Memosa-FVM  0.2
phononbase/COMETModel.h
Go to the documentation of this file.
1 // This file os part of FVM
2 // Copyright (c) 2012 FVM Authors
3 // See LICENSE file for terms.
4 
5 #ifndef _COMETMODEL_H_
6 #define _COMETMODEL_H_
7 
8 #include "Model.h"
9 #include "Array.h"
10 #include "Mesh.h"
11 #include "Kspace.h"
12 #include "kvol.h"
13 #include "pmode.h"
14 #include "COMETDiscretizer.h"
15 #include "NumType.h"
16 #include "COMETBC.h"
17 #include "PhononMacro.h"
18 #include "COMETBoundary.h"
19 #include "COMETInterface.h"
20 
21 template<class T>
22 class COMETModel : public Model
23 {
24 
25  public:
26  typedef typename NumTypeTraits<T>::T_Scalar T_Scalar;
27  typedef Vector<T_Scalar,3> VectorT3;
28  typedef Array<VectorT3> VectorT3Array;
29  typedef shared_ptr<VectorT3Array> VT3Ptr;
30  typedef Kspace<T> Tkspace;
31  typedef Tkspace* TkspPtr;
32  typedef vector<Tkspace*> TkspList;
33  typedef kvol<T> Tkvol;
34  typedef pmode<T> Tmode;
35  typedef Array<bool> BArray;
36  typedef Array<int> IntArray;
37  typedef shared_ptr<IntArray> IntArrPtr;
38  typedef Array<T> TArray;
39  typedef shared_ptr<TArray> TArrptr;
40  typedef map<int,COMETBC<T>*> COMETBCMap;
41  typedef vector<COMETIC<T>*> IClist;
42  typedef map<int,int> MeshKspaceMap;
43  typedef COMETModel<T> TCOMET;
44  typedef shared_ptr<TCOMET> TCOMETPtr;
45  typedef COMETModelOptions<T> TCModOpts;
46  typedef typename Tmode::Mode_ptr Mode_ptr;
47  typedef typename Tmode::Reflection Reflection;
48  typedef typename Tmode::Reflptr Reflptr;
49  typedef typename Tmode::Refl_pair Refl_pair;
50  typedef typename Tmode::Refl_Map Refl_Map;
51  typedef shared_ptr<MeshList> MshLstPtr;
52  typedef shared_ptr<Mesh> MeshPtr;
53  typedef shared_ptr<GeomFields> GeoFldsPtr;
54  typedef shared_ptr<Tkspace> KspPtr;
55  typedef shared_ptr<PhononMacro> PMacroPtr;
56  typedef shared_ptr<StorageSite> SSPtr;
57  typedef shared_ptr<CRConnectivity> CRPtr;
58  typedef Array<int> BCfaceArray;
59  typedef shared_ptr<BCfaceArray> BfacePtr;
60  typedef vector<BfacePtr> BCfaceList;
61  typedef Array<int> BCcellArray;
62  typedef shared_ptr<BCcellArray> BCellPtr;
63  typedef vector<BCellPtr> BCcellList;
64  typedef vector<IntArrPtr> MeshICmap;
65  typedef KSConnectivity<T> TKConnectivity;
66  typedef TKConnectivity* TKCptr;
67  typedef vector<TKCptr> TKClist;
68  typedef map<int, TKClist> FgTKClistMap;
69  typedef vector<shared_ptr<Field> > FieldVector;
70  typedef pair<const StorageSite*, const StorageSite*> SSPair;
71  typedef map<SSPair,shared_ptr<Array<int> > > ScatGathMaps;
72  typedef map<const StorageSite*, StorageSite*> SiteMap; //fine to coarse
73 
74  COMETModel(const MeshList& meshes, const int level, GeomFields& geomFields,
75  TkspList& kspaces, PhononMacro& macro):
76  Model(meshes),
77  _level(level),
78  _geomFields(geomFields),
79  _kspaces(kspaces),
80  _macro(macro),
81  _finestLevel(NULL),
82  _coarserLevel(NULL),
83  _finerLevel(NULL),
84  _residual(1.0),
85  _MeshToIC(),
86  _rank(-1)
87  {
88  const int numMeshes = _meshes.size();
89 
90 #ifdef FVM_PARALLEL
91  _rank=MPI::COMM_WORLD.Get_rank();
92 #endif
93 
94  for (int n=0; n<numMeshes; n++)
95  {
96  Mesh& mesh = *_meshes[n];
97  const StorageSite& faces=mesh.getFaces();
98  const StorageSite& cells=mesh.getCells();
99  const int faceCount=faces.getCount();
100  const int cellCount=cells.getCount();
101  _MeshKspaceMap[n]=-1; //mesh map not set
102 
103  BfacePtr BFptr(new BCfaceArray(faceCount));
104  BFptr->zero();
105  _BFaces.push_back(BFptr);
106 
107  BCellPtr BCptr(new BCcellArray(cellCount));
108  _BCells.push_back(BCptr);
109  BCptr->zero();
110 
111  if(_level==0)
112  {
113  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
114  {
115  const FaceGroup& fg = *fgPtr;
116  if (_bcMap.find(fg.id) == _bcMap.end() && fg.id>0)
117  {
118  COMETBC<T> *bc(new COMETBC<T>());
119  _bcMap[fg.id] = bc;
120  }
121  }
122  }
123  }
124  }
125 
126  TCModOpts& getOptions() {return _options;}
127  COMETBCMap& getBCs() {return _bcMap;}
128  MeshKspaceMap& getMKMap() {return _MeshKspaceMap;}
129 
130  void init()
131  {
132  const int numMeshes=_meshes.size();
133  const T Tinit=_options["initialTemperature"];
134 
135  for (int n=0;n<numMeshes;n++)
136  {
137  Mesh& mesh=*_meshes[n];
138  if(_MeshKspaceMap[n]==-1)
139  throw CException("Have not set the Kspace for this Mesh!!");
140  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
141  const int kcount=kspace.gettotmodes();
142  BCfaceArray& BCfArray=*(_BFaces[n]);
143  BCcellArray& BCArray=*(_BCells[n]);
144  const int numK=kspace.getlength();
145  const StorageSite& cells=mesh.getCells();
146  const int numcells=cells.getCount();
147  VectorT3 lamTemp;
148 
149  // set micro parameters
150  //cout<<"Allocating Arrays..."<<endl;
151  shared_ptr<TArray> eArray(new TArray(numcells*kcount));
152  shared_ptr<TArray> e0Array(new TArray(numcells*kcount));
153  shared_ptr<TArray> ResidArray(new TArray(numcells*kcount));
154  shared_ptr<TArray> tauArray(new TArray(numcells*kcount));
155  kspace.seteArray(eArray);
156  kspace.sete0Array(e0Array);
157  kspace.setResArray(ResidArray);
158  kspace.setTauArray(tauArray);
159 
160  if(_options.Source)
161  {
162  shared_ptr<TArray> SrcArray(new TArray(numcells*kcount));
163  SrcArray->zero();
164  kspace.setSourceArray(SrcArray);
165  }
166 
167  shared_ptr<TArray> TLcell(new TArray(numcells));
168  shared_ptr<TArray> deltaTcell(new TArray(numcells));
169  shared_ptr<VectorT3Array> lamArray(new VectorT3Array(numcells));
170  shared_ptr<VectorT3Array> q(new VectorT3Array(numcells));
171  lamTemp.zero();
172  q->zero();
173  *lamArray=lamTemp;
174  *deltaTcell=0.;
175  *TLcell=Tinit;
176  _macro.temperature.addArray(cells,TLcell);
177  _macro.deltaT.addArray(cells,deltaTcell);
178  _macro.lam.addArray(cells,lamArray);
179  _macro.heatFlux.addArray(cells,q);
180 
181  shared_ptr<IntArray> f2c(new IntArray(numcells));
182  *f2c=-1;
183  _geomFields.fineToCoarse.addArray(cells, f2c);
184 
185  int modeCount=kspace.getkvol(0).getmodenum();
186  FieldVector* FieldVecPtr=new FieldVector();
187 
188  for(int mode=0;mode<modeCount;mode++)
189  {
190  shared_ptr<Field> modeField(new Field("mode"));
191  shared_ptr<TArray> modeTemp(new TArray(numcells));
192  *modeTemp=Tinit;
193  modeField->addArray(cells,modeTemp);
194  FieldVecPtr->push_back(modeField);
195  }
196 
197  _macro.BranchTemperatures[n]=FieldVecPtr;
198 
199  //cout<<"Arrays Allocated...Initializing Values..."<<endl;
200 
201  for(int c=0;c<numcells;c++)
202  {
203  int cellIndex=kspace.getGlobalIndex(c,0);
204  for (int k=0;k<numK;k++)
205  {
206  Tkvol& kv=kspace.getkvol(k);
207  const int numM=kv.getmodenum();
208  //const T dk3=kv.getdk3();
209 
210  for (int m=0;m<numM;m++)
211  {
212  Tmode& mode=kv.getmode(m);
213  T tau=mode.gettau();
214  const T einit=mode.calce0(Tinit);
215  /*
216  if(m==0 && k==0)
217  (*eArray)[cellIndex]=1.1*einit;
218  else*/
219 
220  (*eArray)[cellIndex]=einit;
221  (*e0Array)[cellIndex]=einit;
222  (*ResidArray)[cellIndex]=0.;
223  (*tauArray)[cellIndex]=tau;
224  cellIndex++;
225  }
226  }
227  kspace.updateTau(c,Tinit);
228  }
229 
230 
231  shared_ptr<TArray> TlResidCell(new TArray(numcells));
232  *TlResidCell=0.;
233  _macro.TlResidual.addArray(cells,TlResidCell);
234 
235  //cout<<"Values Initialized...Setting Facegroups..."<<endl;
236 
237  //setting facegroups
238 
239  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
240  {
241  FaceGroup& fg = *fgPtr;
242  const StorageSite& faces = fg.site;
243  //const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
244 
245  const int faceCount=faces.getCount();
246  const int offSet=faces.getOffset();
247 
248  for(int i=offSet;i<offSet+faceCount;i++)
249  BCfArray[i]=-1; //implicit boundary
250  }
251 
252 
253  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
254  {
255  FaceGroup& fg = *fgPtr;
256  if(fg.id>0)
257  {
258  if(_bcMap[fg.id]->bcType == "reflecting")
259  {
260  const StorageSite& faces = fg.site;
261  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
262  const int faceCount=faces.getCount();
263  const int offSet=faces.getOffset();
264  const bool Imp=(*(_bcMap[fg.id]))["FullyImplicit"];
265  const T ref=(*(_bcMap[fg.id]))["specifiedReflection"];
266  const Field& AreaMagField=_geomFields.areaMag;
267  const TArray& AreaMag=
268  dynamic_cast<const TArray&>(AreaMagField[faces]);
269  const Field& AreaDirField=_geomFields.area;
270  const VectorT3Array& AreaDir=
271  dynamic_cast<const VectorT3Array&>(AreaDirField[faces]);
272 
273  const VectorT3 norm=AreaDir[0]/AreaMag[0];
274 
275  if(ref==1.)
276  {
277  for (int k=0;k<numK;k++)
278  {
279  Tkvol& kv=kspace.getkvol(k);
280  const int numM=kv.getmodenum();
281  //const T dk3=kv.getdk3();
282  for (int m=0;m<numM;m++)
283  {
284 
285  Tmode& mode=kv.getmode(m);
286  const VectorT3 vg=mode.getv();
287  const T vmag=sqrt(pow(vg[0],2)+pow(vg[1],2)+pow(vg[2],2));
288  VectorT3 si=vg/vmag;
289  VectorT3 so;
290  const T sidotn=si[0]*norm[0]+si[1]*norm[1]+si[2]*norm[2];
291  Refl_Map& rmap=mode.getreflmap();
292 
293  if (sidotn > T_Scalar(0.0))
294  {
295  so=si-2.*(si[0]*norm[0]+si[1]*norm[1]+si[2]*norm[2])*norm;
296  T soMag=sqrt(pow(so[0],2)+pow(so[1],2)+pow(so[2],2));
297  so/=soMag;
298  so*=vmag;
299  Refl_pair refls;
300  Refl_pair reflsFrom;
301  kspace.findSpecs(norm,m,k,refls);
302  rmap[fg.id]=refls;
303  const int k1=refls.first.second;
304  Tmode& mode2=kspace.getkvol(k1).getmode(m);
305  Refl_Map& rmap2=mode2.getreflmap();
306  reflsFrom.first.second=-1;
307  reflsFrom.second.second=k;
308  rmap2[fg.id]=reflsFrom;
309  }
310  }
311  }
312  }
313 
314  if(Imp)
315  {
316  for(int i=offSet;i<offSet+faceCount;i++)
317  BCfArray[i]=2; //implicit boundary
318  }
319  else
320  {
321  for(int i=offSet;i<offSet+faceCount;i++)
322  BCfArray[i]=3; //explicit boundary
323  }
324 
325  if(Imp)
326  {
327  for(int i=0;i<faceCount;i++)
328  {
329  int cell1=BfaceCells(i,0);
330  if(BCArray[cell1]==0)
331  BCArray[cell1]=1; //implicit boundary only
332  else if(BCArray[cell1]==2)
333  BCArray[cell1]=3; //mix implicit/explicit boundary
334  }
335  }
336  else
337  {
338  for(int i=0;i<faceCount;i++)
339  {
340  int cell1=BfaceCells(i,0);
341  if(BCArray[cell1]==0)
342  BCArray[cell1]=2; //explicit boundary only
343  else if (BCArray[cell1]==1)
344  BCArray[cell1]=3; //mix implicit/explicit boundary
345  }
346  }
347  }
348  else if(_bcMap[fg.id]->bcType == "temperature")
349  {
350  const StorageSite& faces = fg.site;
351  const int faceCount=faces.getCount();
352  const int offSet=faces.getOffset();
353 
354  for(int i=offSet;i<offSet+faceCount;i++)
355  BCfArray[i]=1;
356  }
357  else if(_bcMap[fg.id]->bcType == "Interface")
358  {
359  StorageSite& faces0 = fg.site;
360  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces0);
361  const int faceCount=faces0.getCount();
362  const int offSet=faces0.getOffset();
363 
364  for(int i=offSet;i<offSet+faceCount;i++)
365  {
366  BCfArray[i]=4; //Interface boundary
367  int cell0=BfaceCells(i-offSet,0);
368  int cell1=BfaceCells(i-offSet,1);
369  BCArray[cell0]=2; //always treated explicitly
370  BCArray[cell1]=4; //interface ghost cells need to be labeled
371  }
372 
373  bool doneAlready=false;
374  foreach(const COMETIC<T>* icPtr, _IClist)
375  {
376  if(icPtr->FgID1==fg.id && icPtr->MeshID1==mesh.getID())
377  {
378  doneAlready=true;
379  break;
380  }
381  }
382 
383  bool foundMatch=false;
384  if(!doneAlready)
385  {
386  StorageSite* faces1Ptr=NULL;
387  int otherFgID,otherMid;
388  for(int otherMeshID=n+1;otherMeshID<numMeshes;otherMeshID++)
389  {
390  const Mesh& otherMesh=*_meshes[otherMeshID];
391  foreach(const FaceGroupPtr otherfgPtr, otherMesh.getBoundaryFaceGroups())
392  {
393  foundMatch=mesh.COMETfindCommonFaces(faces0, otherfgPtr->site, _geomFields);
394  if(foundMatch)
395  {
396  otherFgID=otherfgPtr->id;
397  otherMid=otherMeshID;
398  faces1Ptr=&(otherfgPtr->site);
399  break;
400  }
401  }
402  if(foundMatch)
403  break;
404  }
405 
406  if(foundMatch)
407  {
408  StorageSite& faces1=*faces1Ptr;
409  const Mesh& mesh1=*_meshes[otherMid];
410  COMETIC<T>* icPtr(new COMETIC<T>(n,fg.id,mesh1.getID(),
411  otherFgID, faces1.getCount()));
412 
413  icPtr->InterfaceModel=_bcMap[fg.id]->InterfaceModel;
414 
415  setLocalScatterMaps(mesh, faces0, mesh1, faces1);
416 
417  _IClist.push_back(icPtr);
418  }
419  else if(!doneAlready && !foundMatch)
420  {
421  cout<<"Face Group: "<<fg.id<<" MeshID: "<<mesh.getID()<<endl;
422  throw CException("Could not find a matching face group!");
423  }
424  }
425  }// end if interface
426  }
427  }//end facegroup loop
428 
429  //cout<<"Facegroups Set...Mesh "<<n<<" Complete."<<endl;
430 
431  }//end meshes loop
432 
433  //cout<<"Mesh Loop Complete...Creating Interfaces (if any)..."<<endl;
434 
435  //Make map from mesh to IC
436  IntArray ICcount(numMeshes);
437  ICcount.zero();
438 
439  foreach(const COMETIC<T>* icPtr, _IClist)
440  {
441  ICcount[icPtr->MeshID0]++;
442  ICcount[icPtr->MeshID1]++;
443  }
444 
445  if(!_IClist.empty())
446  {
447  for(int i=0;i<numMeshes;i++)
448  {
449  IntArrPtr MeshToIC(new IntArray(ICcount[i]));
450  _MeshToIC.push_back(MeshToIC);
451  }
452  }
453 
454  ICcount.zero();
455  const int listSize=_IClist.size();
456  for(int ic=0;ic<listSize;ic++)
457  {
458  const COMETIC<T>* icPtr=_IClist[ic];
459  (*(_MeshToIC[icPtr->MeshID0]))[ICcount[icPtr->MeshID0]]=ic;
460  (*(_MeshToIC[icPtr->MeshID1]))[ICcount[icPtr->MeshID1]]=ic;
461  ICcount[icPtr->MeshID0]++;
462  ICcount[icPtr->MeshID1]++;
463  }
464 
465  COMETInterface<T> ComInt(_meshes,_kspaces,_MeshKspaceMap,_macro,_geomFields);
466 
467  for(int ic=0;ic<listSize;ic++)
468  {
469  COMETIC<T>* icPtr=_IClist[ic];
470  const int mid0=icPtr->MeshID0;
471  const int mid1=icPtr->MeshID1;
472  if(icPtr->InterfaceModel=="DMM")
473  ComInt.makeDMMcoeffs(*icPtr);
474  else if(icPtr->InterfaceModel=="NoInterface")
475  ComInt.makeNoInterfaceCoeffs(*icPtr);
476  ComInt.updateOtherGhost(*icPtr,mid0,false);
477  ComInt.updateOtherGhost(*icPtr,mid1,false);
478  }
479 
480  //cout<<"Interfaces Complete..."<<endl;
481 
482  initializeTemperatureBoundaries();
483  _residual=updateResid(false);
484 
485  if(_options.DomainStats=="Loud")
486  {
487  cout<<"Creating Coarse Levels on rank "<<_rank<<"..."<<endl;
488  cout<<"Level: 0, rank: "<<_rank<<endl<<endl;
489  calcDomainStats();
490  cout<<endl;
491  }
492 
493  MakeCoarseModel(this);
494  setFinestLevel(this);
495  if(_options.DomainStats=="Loud")
496  cout<<"Coarse Levels Completed on rank "<<_rank<<"."<<endl;
497  }
498 
499  void initFromOld()
500  {
501  const int numMeshes=_meshes.size();
502 
503  for (int n=0;n<numMeshes;n++)
504  {
505  Mesh& mesh=*_meshes[n];
506 
507  BCfaceArray& BCfArray=*(_BFaces[n]);
508  BCcellArray& BCArray=*(_BCells[n]);
509  //const StorageSite& cells=mesh.getCells();
510  //const int numcells=cells.getCount();
511 
512  //setting facegroups
513  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
514  {
515  FaceGroup& fg = *fgPtr;
516  if(fg.id>0)
517  {
518  if(_bcMap[fg.id]->bcType == "reflecting")
519  {
520  const StorageSite& faces = fg.site;
521  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
522  const int faceCount=faces.getCount();
523  const int offSet=faces.getOffset();
524  const bool Imp=(*(_bcMap[fg.id]))["FullyImplicit"];
525 
526  if(Imp)
527  {
528  for(int i=offSet;i<offSet+faceCount;i++)
529  BCfArray[i]=2; //implicit boundary
530  }
531  else
532  {
533  for(int i=offSet;i<offSet+faceCount;i++)
534  BCfArray[i]=3; //explicit boundary
535  }
536 
537  if(Imp)
538  {
539  for(int i=0;i<faceCount;i++)
540  {
541  int cell1=BfaceCells(i,0);
542  if(BCArray[cell1]==0)
543  BCArray[cell1]=1; //implicit boundary only
544  else if(BCArray[cell1]==2)
545  BCArray[cell1]=3; //mix implicit/explicit boundary
546  }
547  }
548  else
549  {
550  for(int i=0;i<faceCount;i++)
551  {
552  int cell1=BfaceCells(i,0);
553  if(BCArray[cell1]==0)
554  BCArray[cell1]=2; //explicit boundary only
555  else if (BCArray[cell1]==1)
556  BCArray[cell1]=3; //mix implicit/explicit boundary
557  }
558  }
559  }
560  else if(_bcMap[fg.id]->bcType == "temperature")
561  {
562  const StorageSite& faces = fg.site;
563  const int faceCount=faces.getCount();
564  const int offSet=faces.getOffset();
565 
566  for(int i=offSet;i<offSet+faceCount;i++)
567  BCfArray[i]=1;
568  }
569  else if(_bcMap[fg.id]->bcType == "Interface")
570  {
571  StorageSite& faces0 = fg.site;
572  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces0);
573  const int faceCount=faces0.getCount();
574  const int offSet=faces0.getOffset();
575 
576  for(int i=offSet;i<offSet+faceCount;i++)
577  {
578  BCfArray[i]=4; //Interface boundary
579  int cell0=BfaceCells(i-offSet,0);
580  int cell1=BfaceCells(i-offSet,1);
581  BCArray[cell0]=2; //always treated explicitly
582  BCArray[cell1]=4; //interface ghost cells need to be labeled
583  }
584 
585  bool doneAlready=false;
586  foreach(const COMETIC<T>* icPtr, _IClist)
587  {
588  if(icPtr->FgID1==fg.id && icPtr->MeshID1==mesh.getID())
589  {
590  doneAlready=true;
591  break;
592  }
593  }
594 
595  bool foundMatch=false;
596  if(!doneAlready)
597  {
598  StorageSite* faces1Ptr=NULL;
599  int otherFgID,otherMid;
600  for(int otherMeshID=n+1;otherMeshID<numMeshes;otherMeshID++)
601  {
602  const Mesh& otherMesh=*_meshes[otherMeshID];
603  foreach(const FaceGroupPtr otherfgPtr, otherMesh.getBoundaryFaceGroups())
604  {
605  foundMatch=mesh.COMETfindCommonFaces(faces0, otherfgPtr->site, _geomFields);
606  if(foundMatch)
607  {
608  otherFgID=otherfgPtr->id;
609  otherMid=otherMeshID;
610  faces1Ptr=&(otherfgPtr->site);
611  break;
612  }
613  }
614  if(foundMatch)
615  break;
616  }
617 
618  if(foundMatch)
619  {
620  StorageSite& faces1=*faces1Ptr;
621  const Mesh& mesh1=*_meshes[otherMid];
622  COMETIC<T>* icPtr(new COMETIC<T>(n,fg.id,mesh1.getID(),
623  otherFgID, faces1.getCount()));
624 
625  if(_bcMap[fg.id]->InterfaceModel!="DMM")
626  icPtr->InterfaceModel=_bcMap[fg.id]->InterfaceModel;
627 
628  _IClist.push_back(icPtr);
629  }
630  else if(!doneAlready && !foundMatch)
631  {
632  cout<<"Face Group: "<<fg.id<<" MeshID: "<<mesh.getID()<<endl;
633  throw CException("Could not find a matching face group!");
634  }
635  }
636  }// end if interface
637  }
638  }//end facegroup loop
639 
640  }//end meshes loop
641 
642  //Make map from mesh to IC
643  IntArray ICcount(numMeshes);
644  ICcount.zero();
645 
646  foreach(const COMETIC<T>* icPtr, _IClist)
647  {
648  ICcount[icPtr->MeshID0]++;
649  ICcount[icPtr->MeshID1]++;
650  }
651 
652  for(int i=0;i<numMeshes;i++)
653  {
654  IntArrPtr MeshToIC(new IntArray(ICcount[i]));
655  _MeshToIC.push_back(MeshToIC);
656  }
657 
658  ICcount.zero();
659  const int listSize=_IClist.size();
660  for(int ic=0;ic<listSize;ic++)
661  {
662  const COMETIC<T>* icPtr=_IClist[ic];
663  (*(_MeshToIC[icPtr->MeshID0]))[ICcount[icPtr->MeshID0]]=ic;
664  (*(_MeshToIC[icPtr->MeshID1]))[ICcount[icPtr->MeshID1]]=ic;
665  ICcount[icPtr->MeshID0]++;
666  ICcount[icPtr->MeshID1]++;
667  }
668 
669  COMETInterface<T> ComInt(_meshes,_kspaces,_MeshKspaceMap,_macro,_geomFields);
670 
671  for(int ic=0;ic<listSize;ic++)
672  {
673  COMETIC<T>* icPtr=_IClist[ic];
674  if(icPtr->InterfaceModel=="DMM")
675  ComInt.makeDMMcoeffs(*icPtr);
676  }
677 
678  initializeTemperatureBoundaries();
679  _residual=updateResid(false);
680  //cout<<"Creating Coarse Levels..."<<endl;
681  MakeCoarseModel(this);
682  //cout<<"Coarse Levels Completed."<<endl;
683  }
684 
685  void initCoarse()
686  {
687  const int numMeshes=_meshes.size();
688  const T Tinit=_options["initialTemperature"];
689 
690  for (int n=0;n<numMeshes;n++)
691  {
692  Mesh& mesh=*_meshes[n];
693  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
694  const int numK=kspace.getlength();
695  const StorageSite& cells=mesh.getCells();
696  const int numcells=cells.getCount();
697  const int kcount=kspace.gettotmodes();
698 
699  shared_ptr<TArray> eArray(new TArray(numcells*kcount));
700  shared_ptr<TArray> e0Array(new TArray(numcells*kcount));
701  shared_ptr<TArray> ResidArray(new TArray(numcells*kcount));
702  shared_ptr<TArray> InjArray(new TArray(numcells*kcount));
703  shared_ptr<TArray> FASArray(new TArray(numcells*kcount));
704  shared_ptr<TArray> tauArray(new TArray(numcells*kcount));
705  kspace.seteArray(eArray);
706  kspace.sete0Array(e0Array);
707  kspace.setResArray(ResidArray);
708  kspace.setInjArray(InjArray);
709  kspace.setFASArray(FASArray);
710  kspace.setTauArray(tauArray);
711 
712  if(_options.Source)
713  {
714  shared_ptr<TArray> SrcArray(new TArray(numcells*kcount));
715  SrcArray->zero();
716  kspace.setSourceArray(SrcArray);
717  }
718 
719  shared_ptr<TArray> TLcell(new TArray(numcells));
720  shared_ptr<TArray> deltaTcell(new TArray(numcells));
721  *deltaTcell=0.;
722  *TLcell=Tinit;
723  _macro.temperature.addArray(cells,TLcell);
724  _macro.deltaT.addArray(cells,deltaTcell);
725 
726  shared_ptr<IntArray> f2c(new IntArray(numcells));
727  *f2c=-1;
728  _geomFields.fineToCoarse.addArray(cells, f2c);
729 
730  for(int c=0;c<numcells;c++)
731  {
732  int cellIndex=kspace.getGlobalIndex(c,0);
733  for (int k=0;k<numK;k++)
734  {
735  Tkvol& kv=kspace.getkvol(k);
736  const int numM=kv.getmodenum();
737 
738  for (int m=0;m<numM;m++)
739  {
740  Tmode& mode=kv.getmode(m);
741  const T einit=mode.calce0(Tinit);
742  T tau=mode.gettau();
743  (*eArray)[cellIndex]=einit;
744  (*e0Array)[cellIndex]=einit;
745  (*ResidArray)[cellIndex]=0.;
746  (*InjArray)[cellIndex]=0.;
747  (*FASArray)[cellIndex]=0.;
748  (*tauArray)[cellIndex]=tau;
749  cellIndex++;
750 
751  }
752  }
753  kspace.updateTau(c,Tinit);
754  }
755 
756  BCfaceArray& BCfArray=*(_BFaces[n]);
757  BCcellArray& BCArray=*(_BCells[n]);
758 
759  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
760  {
761  FaceGroup& fg = *fgPtr;
762  const StorageSite& faces = fg.site;
763  //const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
764 
765  const int faceCount=faces.getCount();
766  const int offSet=faces.getOffset();
767 
768  for(int i=offSet;i<offSet+faceCount;i++)
769  BCfArray[i]=-1;
770  }
771 
772  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
773  {
774  FaceGroup& fg = *fgPtr;
775  if(fg.id>0)
776  {
777  if(_bcMap[fg.id]->bcType == "reflecting")
778  {
779  const StorageSite& faces = fg.site;
780  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
781  const int faceCount=faces.getCount();
782  const int offSet=faces.getOffset();
783  const bool Imp=(*(_bcMap[fg.id]))["FullyImplicit"];
784 
785  if(Imp)
786  {
787  for(int i=offSet;i<offSet+faceCount;i++)
788  BCfArray[i]=2; //implicit boundary
789  }
790  else
791  {
792  for(int i=offSet;i<offSet+faceCount;i++)
793  BCfArray[i]=3; //explicit boundary
794  }
795 
796  if(Imp)
797  {
798  for(int i=0;i<faceCount;i++)
799  {
800  int cell1=BfaceCells(i,0);
801  if(BCArray[cell1]==0)
802  BCArray[cell1]=1; //implicit boundary only
803  else if(BCArray[cell1]==2)
804  BCArray[cell1]=3; //mix implicit/explicit boundary
805  }
806  }
807  else
808  {
809  for(int i=0;i<faceCount;i++)
810  {
811  int cell1=BfaceCells(i,0);
812  if(BCArray[cell1]==0)
813  BCArray[cell1]=2; //explicit boundary only
814  else if (BCArray[cell1]==1)
815  BCArray[cell1]=3; //mix implicit/explicit boundary
816  }
817  }
818  }
819  else if(_bcMap[fg.id]->bcType == "temperature")
820  {
821  const StorageSite& faces = fg.site;
822  const int faceCount=faces.getCount();
823  const int offSet=faces.getOffset();
824 
825  for(int i=offSet;i<offSet+faceCount;i++)
826  BCfArray[i]=1;
827  }
828  else if(_bcMap[fg.id]->bcType == "Interface")
829  {
830  StorageSite& faces0 = fg.site;
831  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces0);
832  const int faceCount=faces0.getCount();
833  const int offSet=faces0.getOffset();
834 
835  for(int i=offSet;i<offSet+faceCount;i++)
836  {
837  BCfArray[i]=4; //Interface boundary
838  int cell0=BfaceCells(i-offSet,0);
839  BCArray[cell0]=2; //always treated explicitly
840  }
841 
842 
843 
844  /*
845  bool doneAlready=false;
846  foreach(const COMETIC<T>* icPtr, _IClist)
847  {
848  if(icPtr->FgID1==fg.id && icPtr->MeshID1==mesh.getID())
849  {
850  doneAlready=true;
851  break;
852  }
853  }
854 
855  if(!doneAlready)
856  {
857  bool foundMatch=false;
858  StorageSite* faces1Ptr=NULL;
859  int otherFgID,otherMid;
860  for(int otherMeshID=n+1;otherMeshID<numMeshes;otherMeshID++)
861  {
862  const Mesh& otherMesh=*_meshes[otherMeshID];
863  foreach(const FaceGroupPtr otherfgPtr, otherMesh.getAllFaceGroups())
864  {
865  foundMatch=mesh.COMETfindCommonFaces(faces0, otherfgPtr->site, _geomFields);
866  if(foundMatch)
867  {
868  otherFgID=otherfgPtr->id;
869  otherMid=otherMeshID;
870  faces1Ptr=&(otherfgPtr->site);
871  break;
872  }
873  }
874  if(foundMatch)
875  break;
876  }
877 
878  if(foundMatch)
879  {
880  StorageSite& faces1=*faces1Ptr;
881  const Mesh& mesh1=*_meshes[otherMid];
882  COMETIC<T>* icPtr(new COMETIC<T>(n,fg.id,mesh1.getID(),
883  otherFgID, faces1.getCount()));
884 
885  if(_bcMap[fg.id]->InterfaceModel!="DMM")
886  icPtr->InterfaceModel=_bcMap[fg.id]->InterfaceModel;
887 
888  setLocalScatterMaps(mesh, faces0, mesh1, faces1);
889 
890  _IClist.push_back(icPtr);
891  }
892  else
893  {
894  cout<<"Face Group: "<<fg.id<<" MeshID: "<<mesh.getID()<<endl;
895  throw CException("Could not find a matching face group!");
896  }
897  }
898 */
899  }
900  }
901  }
902 
903  TArrptr TlResidCell(new TArray(numcells));
904  *TlResidCell=0.;
905  _macro.TlResidual.addArray(cells,TlResidCell);
906  TArrptr TlInj(new TArray(numcells));
907  *TlInj=0.;
908  _macro.TlInjected.addArray(cells,TlInj);
909  TArrptr TlFAS(new TArray(numcells));
910  *TlFAS=0.;
911  _macro.TlFASCorrection.addArray(cells,TlFAS);
912  }
913 
914  /*
915  const int listSize=_IClist.size();
916  for(int ic=0;ic<listSize;ic++)
917  {
918  //finer level interface condition
919  COMETIC<T>* icPtr=_IClist[ic];
920  const int meshID0=icPtr->MeshID0;
921  const int meshID1=icPtr->MeshID1;
922  const int fgid0=icPtr->FgID0;
923  const int fgid1=icPtr->FgID1;
924 
925  //this level's interface condition
926  Mesh& mesh0this=*_meshes[meshID0];
927  const FaceGroup& fg0=mesh0this.getFaceGroup(fgid0);
928  const StorageSite& faces0=fg0.site;
929  const int faceCount=faces0.getCount();
930  COMETIC<T>* icPtr2(new COMETIC<T>(meshID0,fgid0,meshID1,
931  fgid1, faceCount));
932  _IClist[ic]=icPtr2;
933 
934  }
935 
936 
937  //Make map from mesh to IC
938  IntArray ICcount(numMeshes);
939  ICcount.zero();
940 
941  foreach(const COMETIC<T>* icPtr, _IClist)
942  {
943  ICcount[icPtr->MeshID0]++;
944  ICcount[icPtr->MeshID1]++;
945  }
946 
947  for(int i=0;i<numMeshes;i++)
948  {
949  IntArrPtr MeshToIC(new IntArray(ICcount[i]));
950  _MeshToIC.push_back(MeshToIC);
951  }
952 
953  ICcount.zero();
954  */
955  const int listSize=_IClist.size();
956  for(int i=0;i<listSize;i++)
957  {
958  COMETIC<T> ic=*_IClist[i];
959  const int m0=ic.MeshID0;
960  const int m1=ic.MeshID1;
961  const int fgid0=ic.FgID0;
962  const int fgid1=ic.FgID1;
963 
964  Mesh& mesh0=*_meshes[m0];
965  Mesh& mesh1=*_meshes[m1];
966 
967  FaceGroup& fg0=const_cast<FaceGroup&>(mesh0.getFaceGroup(fgid0));
968  FaceGroup& fg1=const_cast<FaceGroup&>(mesh1.getFaceGroup(fgid1));
969 
970  setLocalScatterMaps(mesh0,fg0.site,mesh1,fg1.site);
971  }
972 
973  /*
974  COMETInterface<T> ComInt(_meshes,_kspaces,_MeshKspaceMap,_macro,_geomFields);
975 
976  for(int ic=0;ic<listSize;ic++)
977  {
978  COMETIC<T>* icPtr=_IClist[ic];
979  if(icPtr->InterfaceModel=="DMM")
980  ComInt.makeDMMcoeffs(*icPtr);
981 
982  }
983  */
984 
985  applyTemperatureBoundaries();
986  }
987 
988  void setLocalScatterMaps(const Mesh& mesh0, StorageSite& faces0, const Mesh& mesh1,StorageSite& faces1)
989  {//makes the IntArray for the indices of the cells1 to which the cells0 scatter
990 
991  const IntArray& common01=*(faces0.getCommonMap()[&faces1]);
992  IntArrPtr scatter01(new IntArray(faces0.getCount()));
993  IntArrPtr scatter10(new IntArray(faces0.getCount()));
994  const CRConnectivity& faceCells0=mesh0.getFaceCells(faces0);
995  const CRConnectivity& faceCells1=mesh1.getFaceCells(faces1);
996 
997  for(int i=0;i<faces0.getCount();i++)
998  {
999  const int f01=common01[i];
1000  const int c01=faceCells1(f01,1);
1001  const int c10=faceCells0(i,1);
1002  (*scatter01)[i]=c01;
1003  (*scatter10)[f01]=c10;
1004  }
1005 
1006  faces0.getScatterMap()[&faces1]=scatter01;
1007  faces1.getScatterMap()[&faces0]=scatter10;
1008 
1009  }
1010 
1011  void applyTemperatureBoundaries()
1012  {
1013  const int numMeshes=_meshes.size();
1014 
1015  for(int n=0;n<numMeshes;n++)
1016  {
1017  const Mesh& mesh=*_meshes[n];
1018  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
1019  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1020  {
1021  const FaceGroup& fg = *fgPtr;
1022  if(fg.id>0)
1023  {
1024  const StorageSite& faces = fg.site;
1025  const COMETBC<T>& bc = *_bcMap[fg.id];
1026 
1027  COMETBoundary<T> cbc(faces, mesh,_geomFields,kspace,_options,fg.id);
1028 
1029  if(bc.bcType=="temperature")
1030  {
1031  FloatValEvaluator<T>
1032  bTemperature(bc.getVal("specifiedTemperature"),faces);
1033 
1034  if(_level>0 || _options.Convection=="FirstOrder")
1035  cbc.applyTemperatureWallCoarse(bTemperature);
1036  else
1037  cbc.applyTemperatureWallCoarse(bTemperature);
1038 
1039  }
1040  }
1041  }
1042  }
1043  }
1044 
1045  void initializeTemperatureBoundaries()
1046  {
1047  const int numMeshes=_meshes.size();
1048 
1049  for(int n=0;n<numMeshes;n++)
1050  {
1051  const Mesh& mesh=*_meshes[n];
1052  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
1053  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1054  {
1055  const FaceGroup& fg = *fgPtr;
1056  if(fg.id>0)
1057  {
1058  const StorageSite& faces = fg.site;
1059  const COMETBC<T>& bc = *_bcMap[fg.id];
1060 
1061  COMETBoundary<T> cbc(faces, mesh,_geomFields,kspace,_options,fg.id);
1062 
1063  if(bc.bcType=="temperature")
1064  {
1065  FloatValEvaluator<T>
1066  bTemperature(bc.getVal("specifiedTemperature"),faces);
1067 
1068 
1069  cbc.applyTemperatureWallCoarse(bTemperature);
1070 
1071  if(_options.Convection=="SecondOrder")
1072  cbc.applyTemperatureWallFine(bTemperature);
1073 
1074  }
1075  }
1076  }
1077  }
1078  }
1079 
1080  void MakeCoarseModel(TCOMET* finerModel)
1081  {
1082 
1083  if(_options.AgglomerationMethod=="FaceArea")
1084  {
1085  int maxLevs=finerModel->getOptions().maxLevels;
1086  int thisLevel=(finerModel->getLevel())+1;
1087 
1088  if(thisLevel<maxLevs) //assumes # of levels will always work for the mesh
1089  {
1090  MeshList* newMeshesPtr=new MeshList;
1091  TkspList* newKspacesPtr=new TkspList;
1092  PhononMacro* newMacroPtr=new PhononMacro("coarse");
1093 
1094  MakeNewKspaces(finerModel->getKspaces(),*newKspacesPtr);
1095 
1096  IntArray CoarseCounts(_meshes.size());
1097  IntArray CoarseGhost(_meshes.size());
1098  map<const StorageSite*, IntArray*> PreFacePairMap;
1099  SiteMap siteMap;
1100 
1101  MakeInteriorCoarseMesh(finerModel->getMeshList(), finerModel->getGeomFields(),
1102  *newMeshesPtr, CoarseCounts, PreFacePairMap, CoarseGhost, siteMap);
1103 
1104 #ifdef FVM_PARALLEL
1105  _geomFields.fineToCoarse.syncLocal();
1106  ScatGathMaps coarseScatterMaps;
1107  ScatGathMaps coarseGatherMaps;
1108 
1109  syncGhostCoarsening(finerModel->getMeshList(), finerModel->getGeomFields(),
1110  *newMeshesPtr, coarseScatterMaps, coarseGatherMaps, CoarseCounts,
1111  CoarseGhost);
1112 
1113  makeCoarseScatGath(finerModel->getMeshList(), siteMap, coarseScatterMaps, coarseGatherMaps);
1114 #endif
1115 
1116  int newCount=FinishCoarseMesh(finerModel->getMeshList(), finerModel->getGeomFields(),
1117  *newMeshesPtr, CoarseCounts, PreFacePairMap, CoarseGhost);
1118 
1119  TCOMET* newModelPtr=new COMETModel(*newMeshesPtr,thisLevel,
1120  finerModel->getGeomFields(),
1121  *newKspacesPtr,*newMacroPtr);
1122 
1123  newModelPtr->setFinerLevel(finerModel);
1124  finerModel->setCoarserLevel(newModelPtr);
1125  newModelPtr->getOptions()=finerModel->getOptions();
1126  newModelPtr->getBCs()=finerModel->getBCs();
1127  newModelPtr->getMKMap()=finerModel->getMKMap();
1128  newModelPtr->getIClist()=finerModel->getIClist();
1129  newModelPtr->getMeshICmap()=finerModel->getMeshICmap();
1130  COMETInterface<T> ComInt(_meshes,_kspaces,_MeshKspaceMap,_macro,_geomFields);
1131  ComInt.makeCoarseCoeffs(finerModel->getIClist(),newModelPtr->getIClist(),*newMeshesPtr);
1132  newModelPtr->initCoarse();
1133 
1134  if(_options.DomainStats=="Loud")
1135  {
1136  cout<<"Level: "<<newModelPtr->getLevel()<<endl<<endl;
1137  newModelPtr->calcDomainStats();
1138  cout<<endl;
1139  }
1140 
1141  if(newCount>3)
1142  newModelPtr->MakeCoarseModel(newModelPtr);
1143  else
1144  _options.maxLevels=newModelPtr->getLevel();
1145  }
1146  }
1147  else if(_options.AgglomerationMethod=="AMG")
1148  throw CException("Have not implemented AMG agglomeration method.");
1149  else
1150  throw CException("Unknown agglomeration method.");
1151  }
1152 
1153  void setFinestLevel(TCOMET* finest)
1154  {
1155  _finestLevel=finest;
1156  if(_coarserLevel!=NULL)
1157  _coarserLevel->setFinestLevel(finest);
1158  }
1159 
1160  void MakeNewKspaces(TkspList& inList, TkspList& outList)
1161  {
1162  const int len=inList.size();
1163  for(int i=0;i<len;i++)
1164  {
1165  TkspPtr newKspacePtr=TkspPtr(new Tkspace());
1166  newKspacePtr->CopyKspace(*inList[i]);
1167  inList[i]->setCoarseKspace(newKspacePtr);
1168  outList.push_back(newKspacePtr);
1169  }
1170 
1171  for(int i=0;i<len;i++)
1172  inList[i]->giveTransmissions();
1173 
1174  }
1175 
1176  void MakeInteriorCoarseMesh(const MeshList& inMeshes, GeomFields& inGeomFields,
1177  MeshList& outMeshes, IntArray& CoarseCounts,
1178  map<const StorageSite*, IntArray*>& PreFacePairMap, IntArray& CoarseGhost,
1179  SiteMap& siteMap)
1180  {
1181 
1182  const int numMeshes=inMeshes.size();
1183  CoarseCounts=0;
1184 
1185  //coarsen interfaces first
1186 
1187  foreach(COMETIC<T>* icPtr, _IClist)
1188  {
1189  COMETIC<T>& ic=*icPtr;
1190  coarsenInterfaceCells(ic, CoarseCounts, inGeomFields, inMeshes);
1191  }
1192 
1193  //coarsen interiors
1194  for(int m=0;m<numMeshes;m++)
1195  {
1196  const Mesh& mesh=*inMeshes[m];
1197  const int dim=mesh.getDimension();
1198  Mesh* newMeshPtr=new Mesh(dim);
1199  outMeshes.push_back(newMeshPtr);
1200  newMeshPtr->setID(m);
1201  siteMap[&(mesh.getCells())]=&(newMeshPtr->getCells());
1202  int coarseCount=coarsenInterior(m, mesh, inGeomFields, CoarseCounts[m]);
1203  }
1204 
1205  foreach(COMETIC<T>* icPtr, _IClist)
1206  {
1207  COMETIC<T>& ic=*icPtr;
1208  const int Mid0=ic.MeshID0;
1209  const int Mid1=ic.MeshID1;
1210  const int Fid0=ic.FgID0;
1211  const int Fid1=ic.FgID1;
1212  const Mesh& mesh0=*inMeshes[Mid0];
1213  const Mesh& mesh1=*inMeshes[Mid1];
1214  const StorageSite& faces0=mesh0.getFaceGroup(Fid0).site;
1215  const StorageSite& faces1=mesh1.getFaceGroup(Fid1).site;
1216  const StorageSite::CommonMap& ComMap = faces0.getCommonMap();
1217  const IntArray& common01=*(ComMap.find(&faces1)->second);
1218  const StorageSite& cells0=mesh0.getCells();
1219  const StorageSite& cells1=mesh1.getCells();
1220  //const int f1Offset=faces1.getOffset();
1221  const StorageSite& allFaces0=mesh0.getFaces();
1222  const CRConnectivity& faceCells0=mesh0.getFaceCells(faces0);
1223  //const CRConnectivity& cellCells0=mesh0.getCellCells();
1224  //const CRConnectivity& cellCells1=mesh1.getCellCells();
1225  const CRPtr cellFaces0Ptr=faceCells0.getTranspose();
1226  const CRConnectivity& cellFaces0=*cellFaces0Ptr;
1227  const CRConnectivity& faceCells1=mesh1.getFaceCells(faces1);
1228  const CRPtr cellFaces1Ptr=faceCells1.getTranspose();
1229  //const CRConnectivity& cellFaces1=*cellFaces1Ptr;
1230  const int faceCount=faces0.getCount();
1231  IntArray& FineToCoarse0=dynamic_cast<IntArray&>(inGeomFields.fineToCoarse[cells0]);
1232  IntArray& FineToCoarse1=dynamic_cast<IntArray&>(inGeomFields.fineToCoarse[cells1]);
1233  const int cCount0=CoarseCounts[Mid0];
1234  //const int cCount1=CoarseCounts[Mid1];
1235 
1236  const CRPtr cellsIntGhost0Ptr=cellFaces0.multiply(faceCells0,true);
1237  const CRConnectivity& cellsIntGhost0=*cellsIntGhost0Ptr;
1238 
1239  Field& FaceAreaField=inGeomFields.area;
1240  const VectorT3Array& FArea0=
1241  dynamic_cast<const VectorT3Array&>(FaceAreaField[allFaces0]); //global ordering
1242 
1243  IntArray* FaceFine2CoarsePtr0=new IntArray(faceCount);
1244  *FaceFine2CoarsePtr0=-1;
1245  IntArray& FaceFine2Coarse0=*FaceFine2CoarsePtr0;
1246  ic.FineToCoarse0=FaceFine2CoarsePtr0;
1247 
1248  IntArray* FaceFine2CoarsePtr1=new IntArray(faceCount);
1249  *FaceFine2CoarsePtr1=-1;
1250  IntArray& FaceFine2Coarse1=*FaceFine2CoarsePtr1;
1251  ic.FineToCoarse1=FaceFine2CoarsePtr1;
1252 
1253  StorageSite preCoarse0(cCount0);
1254 
1255  CRPtr preCoarseToFineCells0=CRPtr(new CRConnectivity(preCoarse0,cells0));
1256 
1257  preCoarseToFineCells0->initCount();
1258 
1259  for(int c=0;c<cells0.getSelfCount();c++)
1260  preCoarseToFineCells0->addCount(FineToCoarse0[c],1);
1261 
1262  preCoarseToFineCells0->finishCount();
1263 
1264  for(int c=0;c<cells0.getSelfCount();c++)
1265  preCoarseToFineCells0->add(FineToCoarse0[c],c);
1266 
1267  preCoarseToFineCells0->finishAdd();
1268 
1269  //initial pairing, no regard to zero summed area
1270  int coarseFace(0);
1271  for(int f=0;f<faceCount;f++)
1272  {
1273 
1274  if(FaceFine2Coarse0[f]==-1)
1275  {
1276  const int c00=faceCells0(f,0); //interior cells
1277  const int f01=common01[f];
1278  const int c01=faceCells1(f01,0);
1279  const int coarseC00=FineToCoarse0[c00]; //coarse indices
1280  const int coarseC01=FineToCoarse1[c01];
1281 
1282  const int cC00fineCnt=preCoarseToFineCells0->getCount(coarseC00);
1283  for(int cC00fine=0;cC00fine<cC00fineCnt;cC00fine++)
1284  {
1285  const int fineCell=(*preCoarseToFineCells0)(coarseC00,cC00fine);
1286  if(fineCell!=c00)
1287  {
1288  const int fineGhostCount=cellsIntGhost0.getCount(fineCell);
1289  for(int fineGhostNum=0;fineGhostNum<fineGhostCount;fineGhostNum++)
1290  {
1291  const int fineGhost=cellsIntGhost0(fineCell,fineGhostNum);
1292  const int f10=cellFaces0(fineGhost,0);
1293  const int f11=common01[f10];
1294  const int c11=faceCells1(f11,0);
1295 
1296  if(FineToCoarse1[c11]==coarseC01)
1297  {
1298  if(FaceFine2Coarse0[f10]==-1 && FaceFine2Coarse0[f]==-1)
1299  {
1300  FaceFine2Coarse0[f10]=coarseFace;
1301  FaceFine2Coarse0[f]=coarseFace;
1302  FaceFine2Coarse1[f11]=coarseFace;
1303  FaceFine2Coarse1[f01]=coarseFace;
1304  coarseFace++;
1305  }
1306  else if(FaceFine2Coarse0[f10]==-1 && FaceFine2Coarse0[f]!=-1)
1307  {
1308  FaceFine2Coarse0[f10]=FaceFine2Coarse0[f];
1309  FaceFine2Coarse1[f11]=FaceFine2Coarse0[f01];
1310  }
1311  else if(FaceFine2Coarse0[f10]!=-1 && FaceFine2Coarse0[f]==-1)
1312  {
1313  FaceFine2Coarse0[f]=FaceFine2Coarse0[f10];
1314  FaceFine2Coarse1[f01]=FaceFine2Coarse1[f11];
1315  }
1316  }
1317 
1318  }
1319  }
1320  }
1321 
1322  if(FaceFine2Coarse0[f]==-1)
1323  {
1324  FaceFine2Coarse0[f]=coarseFace;
1325  FaceFine2Coarse1[f01]=coarseFace;
1326  coarseFace++;
1327  }
1328  }
1329 
1330  //FOR NO AGGLOMERATION OF INTERFACE FACES
1331  //FaceFine2Coarse0[f]=coarseFace;
1332  //FaceFine2Coarse1[f]=coarseFace;
1333  //coarseFace++;
1334 
1335 
1336  }
1337 
1338 
1339 
1340  //must check for zero summed area
1341  VectorT3Array summedArea(coarseFace);
1342  summedArea.zero();
1343  for(int f=0;f<faceCount;f++)
1344  {
1345  const int offset=faces0.getOffset();
1346  summedArea[FaceFine2Coarse0[f]]+=FArea0[f+offset];
1347  }
1348 
1349  for(int f=0;f<faceCount;f++)
1350  {
1351  const int f01=common01[f];
1352  const int offset=faces0.getOffset();
1353  VectorT3& sumVec=summedArea[FaceFine2Coarse0[f]];
1354  const VectorT3& partVec=FArea0[f+offset];
1355  T sumMag=sqrt(sumVec[0]*sumVec[0]+sumVec[1]*sumVec[1]+
1356  sumVec[2]*sumVec[2]);
1357  T partMag=sqrt(partVec[0]*partVec[0]+partVec[1]*partVec[1]+
1358  partVec[2]*partVec[2]);
1359 
1361  if(sumMag/partMag<1.0)
1362  {
1363  sumVec-=partVec;
1364  FaceFine2Coarse0[f]=coarseFace;
1365  FaceFine2Coarse1[f01]=coarseFace;
1366  coarseFace++;
1367  }
1368 
1369  }
1370 
1371  PreFacePairMap[&faces0]=FaceFine2CoarsePtr0;
1372  PreFacePairMap[&faces1]=FaceFine2CoarsePtr1;
1373 
1374  }
1375 
1376  for(int n=0;n<numMeshes;n++)
1377  {
1378  const Mesh& mesh=*_meshes[n];
1379 
1380  if(!_IClist.empty())
1381  {
1382  CoarseCounts[n]=correctSingleNeighbor(n, mesh, inGeomFields,
1383  CoarseCounts[n], PreFacePairMap);
1384  }
1385 
1386  IntArray& FineToCoarse=dynamic_cast<IntArray&>(
1387  inGeomFields.fineToCoarse[mesh.getCells()]);
1388 
1389  CoarseGhost[n]=CoarseCounts[n];
1390  int outGhost(0);
1391 
1392  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1393  {
1394  const FaceGroup& fg=*fgPtr;
1395  if(fg.id>0)
1396  {
1397  const CRConnectivity& faceCells=mesh.getFaceCells(fg.site);
1398  const int faceCount=fg.site.getCount();
1399 
1400  if(_bcMap[fg.id]->bcType == "Interface")
1401  {
1402  IntArray& FaceFine2Coarse=*PreFacePairMap[&fg.site];
1403 
1404  int coarseFaces(0);
1405  for(int i=0;i<faceCount;i++)
1406  {
1407  const int cghost=faceCells(i,1);
1408  FineToCoarse[cghost]=FaceFine2Coarse[i]+CoarseGhost[n];
1409  if(FaceFine2Coarse[i]>coarseFaces)
1410  coarseFaces=FaceFine2Coarse[i];
1411  }
1412  CoarseGhost[n]+=coarseFaces+1;
1413  outGhost+=coarseFaces+1;
1414  }
1415  else if(_bcMap[fg.id]->bcType == "temperature" ||
1416  _bcMap[fg.id]->bcType == "reflecting" )
1417  {
1418  for(int i=0;i<faceCount;i++)
1419  {
1420  const int cghost=faceCells(i,1);
1421  FineToCoarse[cghost]=CoarseGhost[n];
1422  CoarseGhost[n]++;
1423  outGhost++;
1424  }
1425  }
1426  }
1427  }
1428 
1429  }
1430 
1431  }
1432 
1433  void syncGhostCoarsening(const MeshList& inMeshes, GeomFields& inGeomFields,
1434  MeshList& outMeshes, ScatGathMaps& coarseScatMaps,
1435  ScatGathMaps& coarseGathMaps, IntArray& coarseSizes,
1436  IntArray& CoarseGhost)
1437  {
1438 
1439  const int numMeshes=inMeshes.size();
1440  for(int n=0;n<numMeshes;n++)
1441  {
1442  const Mesh& mesh=*inMeshes[n];
1443 
1444  const StorageSite& inCells=mesh.getCells();
1445  const StorageSite& site=mesh.getCells();
1446  //const int inCellCount=inCells.getSelfCount();
1447  //const int inCellTotal=inCells.getCount();
1448 
1449  Field& FineToCoarseField=inGeomFields.fineToCoarse;
1450  IntArray& coarseIndex=dynamic_cast<IntArray&>(FineToCoarseField[inCells]);
1451  IntArray tempIndex(inCells.getCountLevel1());
1452  for(int c=0;c<inCells.getCountLevel1();c++)
1453  tempIndex[c]=coarseIndex[c];
1454 
1455  //int coarseGhostSize=0;
1456  //int tempGhostSize=0;
1457  //int coarseSize= CoarseGhost[n];
1458 
1459  const StorageSite::GatherMap& gatherMap = site.getGatherMap();
1460 
1461  // collect all the toIndices arrays for each storage site from
1462  // gatherMap
1463 
1464  typedef map<const StorageSite*, vector<const Array<int>* > > IndicesMap;
1465  IndicesMap toIndicesMap;
1466  //IndicesMap tempIndicesMap;
1467 
1468  foreach(const StorageSite::GatherMap::value_type pos, gatherMap)
1469  {
1470  const StorageSite& oSite = *pos.first;
1471  //const Array<int>& tempIndices = *pos.second;
1472  const Array<int>& toIndices = *pos.second;
1473 
1474  //tempIndicesMap[&oSite].push_back(&tempIndices);
1475  toIndicesMap[&oSite].push_back(&toIndices);
1476  }
1477 
1478  /*
1479  foreach(IndicesMap::value_type pos, tempIndicesMap)
1480  {
1481  const StorageSite& oSite = *pos.first;
1482  const vector<const Array<int>* > tempIndicesArrays = pos.second;
1483 
1484  map<int,int> otherToMyMapping;
1485  UnorderedSet gatherSet;
1486 
1487  foreach(const Array<int>* tempIndicesPtr, tempIndicesArrays)
1488  {
1489  const Array<int>& tempIndices = *tempIndicesPtr;
1490  const int nGhostRows = tempIndices.getLength();
1491  for(int ng=0; ng<nGhostRows; ng++)
1492  {
1493  const int fineIndex = tempIndices[ng];
1494  const int coarseOtherIndex = tempIndex[fineIndex];
1495 
1496  if (coarseOtherIndex < 0)
1497  continue;
1498 
1499  if (otherToMyMapping.find(coarseOtherIndex) !=
1500  otherToMyMapping.end())
1501  {
1502 
1503  tempIndex[fineIndex] = otherToMyMapping[coarseOtherIndex];
1504  }
1505  else
1506  {
1507  tempIndex[fineIndex] = tempGhostSize+coarseSize;
1508  otherToMyMapping[coarseOtherIndex] = tempIndex[fineIndex];
1509  gatherSet.insert( tempIndex[fineIndex] );
1510  tempGhostSize++;
1511  }
1512  }
1513  }
1514  }*/
1515 
1516  foreach(IndicesMap::value_type pos, toIndicesMap)
1517  {
1518  const StorageSite& oSite = *pos.first;
1519  const vector<const Array<int>* > toIndicesArrays = pos.second;
1520 
1521 
1522  map<int,int> otherToMyMapping;
1523  UnorderedSet gatherSet;
1524 
1525  foreach(const Array<int>* toIndicesPtr, toIndicesArrays)
1526  {
1527  const Array<int>& toIndices = *toIndicesPtr;
1528  const int nGhostRows = toIndices.getLength();
1529  for(int ng=0; ng<nGhostRows; ng++)
1530  {
1531  const int fineIndex = toIndices[ng];
1532  const int coarseOtherIndex = coarseIndex[fineIndex];
1533 
1534  if (coarseOtherIndex < 0)
1535  continue;
1536 
1537  if (otherToMyMapping.find(coarseOtherIndex) !=
1538  otherToMyMapping.end())
1539  {
1540  coarseIndex[fineIndex] = otherToMyMapping[coarseOtherIndex];
1541  }
1542  else
1543  {
1544  coarseIndex[fineIndex] = CoarseGhost[n];
1545  otherToMyMapping[coarseOtherIndex] = coarseIndex[fineIndex];
1546  gatherSet.insert( coarseIndex[fineIndex] );
1547  CoarseGhost[n]++;
1548  }
1549 
1550  }
1551 
1552  }
1553 
1554  const int coarseMappersSize = otherToMyMapping.size();
1555 
1556  shared_ptr<Array<int> > coarseToIndices(new Array<int>(coarseMappersSize));
1557 
1558  for(int n = 0; n < gatherSet.size(); n++)
1559  (*coarseToIndices)[n] = gatherSet.getData().at(n);
1560 
1561  SSPair sskey(&site,&oSite);
1562  coarseGathMaps [sskey] = coarseToIndices;
1563 
1564  }
1565 
1566  const StorageSite::ScatterMap& scatterMap = site.getScatterMap();
1567 
1568  IndicesMap fromIndicesMap;
1569 
1570  foreach(const StorageSite::GatherMap::value_type pos, scatterMap)
1571  {
1572  const StorageSite& oSite = *pos.first;
1573  const Array<int>& fromIndices = *pos.second;
1574 
1575  fromIndicesMap[&oSite].push_back(&fromIndices);
1576  }
1577 
1578  foreach(IndicesMap::value_type pos, fromIndicesMap)
1579  {
1580  const StorageSite& oSite = *pos.first;
1581  const vector<const Array<int>* > fromIndicesArrays = pos.second;
1582 
1583  UnorderedSet scatterSet;
1584 
1585  foreach(const Array<int>* fromIndicesPtr, fromIndicesArrays)
1586  {
1587  const Array<int>& fromIndices = *fromIndicesPtr;
1588  const int nGhostRows = fromIndices.getLength();
1589  for(int ng=0; ng<nGhostRows; ng++)
1590  {
1591  const int fineIndex = fromIndices[ng];
1592  const int coarseOtherIndex = coarseIndex[fineIndex];
1593  if (coarseOtherIndex >= 0)
1594  scatterSet.insert( coarseOtherIndex );
1595  }
1596 
1597  }
1598 
1599  const int coarseMappersSize = scatterSet.size();
1600 
1601  shared_ptr<Array<int> > coarseFromIndices(new Array<int>(coarseMappersSize));
1602 
1603  for(int n = 0; n < scatterSet.size(); n++ )
1604  (*coarseFromIndices)[n] = scatterSet.getData().at(n);
1605 
1606  SSPair sskey(&site,&oSite);
1607  coarseScatMaps[sskey] = coarseFromIndices;
1608 
1609  }
1610 
1611  }
1612  }
1613 
1614  void makeCoarseScatGath(const MeshList& inMeshes, SiteMap& siteMap,
1615  ScatGathMaps& coarseScatMaps, ScatGathMaps& coarseGathMaps)
1616  {
1617  const int numMeshes =_meshes.size();
1618  for (int n=0; n<numMeshes; n++)
1619  {
1620  const Mesh& mesh = *_meshes[n];
1621  const StorageSite& fineSite = mesh.getCells();
1622  StorageSite& coarseSite = *siteMap[&fineSite];
1623  const StorageSite::ScatterMap& fineScatterMap = fineSite.getScatterMap();
1624  StorageSite::ScatterMap& coarseScatterMap = coarseSite.getScatterMap();
1625 
1626  foreach(const StorageSite::ScatterMap::value_type& pos, fineScatterMap)
1627  {
1628  const StorageSite& fineOSite = *pos.first;
1629 
1630  // the ghost site will not have its corresponding coarse
1631  // site created yet so we create it here
1632  if (siteMap.find(&fineOSite) == siteMap.end())
1633  {
1634  StorageSite* ghostSite=new StorageSite(-1);
1635  ghostSite->setGatherProcID (fineOSite.getGatherProcID());
1636  ghostSite->setScatterProcID(fineOSite.getScatterProcID());
1637  ghostSite->setTag( fineOSite.getTag() );
1638  siteMap[&fineOSite]=ghostSite;
1639  }
1640 
1641  StorageSite* coarseOSite = siteMap[&fineOSite];
1642 
1643  SSPair sskey(&fineSite,&fineOSite);
1644  coarseScatterMap[coarseOSite] = coarseScatMaps[sskey];
1645 
1646  }
1647 
1648  const StorageSite::GatherMap& fineGatherMap = fineSite.getGatherMap();
1649  StorageSite::GatherMap& coarseGatherMap = coarseSite.getGatherMap();
1650 
1651  foreach(const StorageSite::GatherMap::value_type& pos, fineGatherMap)
1652  {
1653  const StorageSite& fineOSite = *pos.first;
1654  StorageSite& coarseOSite = *siteMap[&fineOSite];
1655  SSPair sskey(&fineSite,&fineOSite);
1656 
1657  coarseGatherMap[&coarseOSite] = coarseGathMaps[sskey];
1658  }
1659 
1660  }
1661  }
1662 
1663  int FinishCoarseMesh(const MeshList& inMeshes, GeomFields& inGeomFields,
1664  MeshList& outMeshes, IntArray& CoarseCounts,
1665  map<const StorageSite*, IntArray*>& PreFacePairMap,
1666  IntArray& coarseGhost)
1667  {
1668 
1669  const int numMeshes=inMeshes.size();
1670 
1671  int smallestMesh=-1;
1672 
1673  for(int n=0;n<numMeshes;n++)
1674  {
1675  const Mesh& mesh=*inMeshes[n];
1676  Mesh* newMeshPtr=outMeshes[n];
1677  //int coarseCount=CoarseCounts[n];
1678  const StorageSite& inCells=mesh.getCells();
1679  IntArray& FineToCoarse=dynamic_cast<IntArray&>(inGeomFields.fineToCoarse[inCells]);
1680  StorageSite& outCells=newMeshPtr->getCells();
1681  StorageSite& outFaces=newMeshPtr->getFaces();
1682  const StorageSite& inFaces=mesh.getFaces();
1683  const int inCellTotal=inCells.getCount();
1684  const int inFaceCount=inFaces.getCount();
1685 
1686  const CRConnectivity& inFaceinCells=mesh.getFaceCells(inFaces);
1687  Field& areaMagField=inGeomFields.areaMag;
1688  const TArray& areaMagArray=
1689  dynamic_cast<const TArray&>(areaMagField[inFaces]);
1690 
1691  Field& FaceAreaField=inGeomFields.area;
1692  const VectorT3Array& inFA=
1693  dynamic_cast<const VectorT3Array&>(FaceAreaField[inFaces]);
1694 
1695  //make the coarse cell to fine cell connectivity.
1696  outCells.setCount(CoarseCounts[n],coarseGhost[n]-CoarseCounts[n]);
1697  CRPtr CoarseToFineCells=CRPtr(new CRConnectivity(outCells,inCells));
1698  CoarseToFineCells->initCount();
1699 
1700  for(int c=0;c<inCellTotal;c++)
1701  CoarseToFineCells->addCount(FineToCoarse[c],1);
1702 
1703  CoarseToFineCells->finishCount();
1704 
1705  for(int c=0;c<inCellTotal;c++)
1706  CoarseToFineCells->add(FineToCoarse[c],c);
1707 
1708  CoarseToFineCells->finishAdd();
1709 
1710  //connectivity between itself (cells) and its finer mesh cells.
1711  newMeshPtr->setConnectivity(outCells,inCells,CoarseToFineCells);
1712 
1713  CRPtr FineFacesCoarseCells=CRPtr(new CRConnectivity(inFaces,outCells));
1714  FineFacesCoarseCells->initCount();
1715 
1716  //count surviving faces
1717  for(int f=0;f<inFaceCount;f++)
1718  {
1719  int coarse0=FineToCoarse[inFaceinCells(f,0)];
1720  int coarse1=FineToCoarse[inFaceinCells(f,1)];
1721  if(coarse0!=coarse1)
1722  FineFacesCoarseCells->addCount(f,2);
1723  }
1724 
1725  FineFacesCoarseCells->finishCount();
1726 
1727  //make non-zero's
1728  for(int f=0;f<inFaceCount;f++)
1729  {
1730  int fc0=inFaceinCells(f,0);
1731  int fc1=inFaceinCells(f,1);
1732  int cc0=FineToCoarse[fc0];
1733  int cc1=FineToCoarse[fc1];
1734  if(cc0!=cc1)
1735  {
1736  FineFacesCoarseCells->add(f,cc0);
1737  FineFacesCoarseCells->add(f,cc1);
1738  }
1739  }
1740 
1741  FineFacesCoarseCells->finishAdd();
1742 
1743  CRPtr CoarseCellsFineFaces=FineFacesCoarseCells->getTranspose();
1744  CRPtr CellCellCoarse=CoarseCellsFineFaces->multiply(*FineFacesCoarseCells,true);
1745 
1746  int counter=0; //coarse face counter
1747  BArray counted(outCells.getCount());
1748  counted=false;
1749  for(int c=0;c<outCells.getCount();c++)
1750  {
1751  counted[c]=true;
1752  const int neibs=CellCellCoarse->getCount(c);
1753 
1754  for(int n=0;n<neibs;n++)
1755  {
1756  const int c1=(*CellCellCoarse)(c,n);
1757  if(neibs>1)
1758  {
1759  if(!counted[c1])
1760  {
1761  counter++;
1762  }
1763  }
1764  else //gives two faces to cells with 1 neighbor
1765  {
1766  if(c>=outCells.getSelfCount())
1767  {
1768  if(!counted[c1])
1769  counter++;
1770  }
1771  else //interior cells
1772  {
1773  if(counted[c1])
1774  counter++;
1775  else
1776  counter+=2;
1777  }
1778  }
1779  }
1780  }
1781 
1782  outFaces.setCount(counter);
1783 
1784  CRPtr CoarseCellCoarseFace=CRPtr(new CRConnectivity(outCells,outFaces));
1785  CoarseCellCoarseFace->initCount();
1786 
1787  int tempCount(0);
1788  for(int c=0;c<outCells.getCount();c++)
1789  {
1790  const int neibs=CellCellCoarse->getCount(c);
1791  CoarseCellCoarseFace->addCount(c,neibs);
1792  tempCount+=neibs;
1793  if((neibs==1) && (c<outCells.getSelfCount())) //adding 2 faces to cells with one neighbor
1794  {
1795  CoarseCellCoarseFace->addCount(c,neibs);
1796  tempCount+=neibs;
1797  }
1798  }
1799 
1800  CoarseCellCoarseFace->finishCount();
1801 
1802  //make cell connectivity to interior faces.
1803  counter=0;
1804  counted=false;
1805  for(int c=0;c<outCells.getSelfCount();c++)
1806  {
1807  counted[c]=true;
1808  const int neibs=CellCellCoarse->getCount(c);
1809  for(int n=0;n<neibs;n++)
1810  {
1811  const int c1=(*CellCellCoarse)(c,n);
1812  if(neibs>1)
1813  {
1814  if(!counted[c1] && c1<outCells.getSelfCount())
1815  {
1816  CoarseCellCoarseFace->add(c,counter);
1817  CoarseCellCoarseFace->add(c1,counter);
1818  counter++;
1819  }
1820  }
1821  else
1822  {
1823  if(counted[c1] && c1<outCells.getSelfCount())
1824  {
1825  CoarseCellCoarseFace->add(c,counter);
1826  CoarseCellCoarseFace->add(c1,counter);
1827  counter++;
1828  }
1829  else if(!counted[c1] && c1<outCells.getSelfCount())
1830  {
1831  CoarseCellCoarseFace->add(c,counter);
1832  CoarseCellCoarseFace->add(c1,counter);
1833  counter++;
1834  CoarseCellCoarseFace->add(c,counter);
1835  CoarseCellCoarseFace->add(c1,counter);
1836  counter++;
1837  }
1838  }
1839  }
1840  }
1841 
1842  const int coarseInteriorCount=counter;
1843 
1844  //cell connectivity to boundary faces
1845  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1846  {
1847  const FaceGroup& fg=*fgPtr;
1848  //const int faceCount=fg.site.getCount();
1849  const CRConnectivity& inBfaceCells=mesh.getFaceCells(fg.site);
1850 
1851  if(fg.id>0)
1852  {
1853  int faceCount(0);
1854  if(_bcMap[fg.id]->bcType == "Interface")
1855  {
1856  IntArray& FaceFine2Coarse=*PreFacePairMap[&fg.site];
1857  int coarseFaces(0);
1858  for(int i=0;i<fg.site.getCount();i++)
1859  if(FaceFine2Coarse[i]>coarseFaces)
1860  coarseFaces=FaceFine2Coarse[i];
1861  faceCount=coarseFaces+1;
1862 
1863  BArray countedFaces(faceCount);
1864  countedFaces=false;
1865  for(int i=0;i<fg.site.getCount();i++)
1866  {
1867  const int coarseFace=FaceFine2Coarse[i];
1868  if(!countedFaces[coarseFace])
1869  {
1870  const int globFace=i+fg.site.getOffset();
1871  const int cc1=(*FineFacesCoarseCells)(globFace,0);
1872  const int cc2=(*FineFacesCoarseCells)(globFace,1);
1873  CoarseCellCoarseFace->add(cc1,counter);
1874  CoarseCellCoarseFace->add(cc2,counter);
1875  counter++;
1876  countedFaces[coarseFace]=true;
1877  }
1878  }
1879 
1880  }
1881  else
1882  {
1883  faceCount=fg.site.getCount();
1884  for(int i=0;i<faceCount;i++)
1885  {
1886  const int c1=inBfaceCells(i,1); //fine ghost cell
1887  const int cc1=FineToCoarse[c1]; //coarse ghost cell
1888  const int cc2=(*CellCellCoarse)(cc1,0); //coarse interior cell
1889  CoarseCellCoarseFace->add(cc1,counter);
1890  CoarseCellCoarseFace->add(cc2,counter);
1891  counter++;
1892  }
1893  }
1894  }
1895  }
1896 
1897  // have to count coarse interfaces (parallel) -- !!
1898 
1899  foreach(const StorageSite::GatherMap::value_type pos, outCells.getGatherMap())
1900  {
1901  const StorageSite& oSite = *pos.first;
1902  const Array<int>& toIndices = *pos.second;
1903 
1904  int to_where = oSite.getGatherProcID();
1905  if ( to_where != -1 )
1906  {
1907  const int fgCount=toIndices.getLength();
1908 
1909  for(int i=0;i<fgCount;i++)
1910  {
1911  const int c1=toIndices[i];
1912  const int c1neibs=CellCellCoarse->getCount(c1);
1913  for(int j=0;j<c1neibs;j++)
1914  {
1915  const int c2=(*CellCellCoarse)(toIndices[i],j);
1916  CoarseCellCoarseFace->add(c2,counter);
1917  CoarseCellCoarseFace->add(toIndices[i],counter);
1918  counter++;
1919  }
1920  }
1921  }
1922  }
1923 
1924  CoarseCellCoarseFace->finishAdd();
1925 
1926  CRPtr CoarseFaceCoarseCell=CoarseCellCoarseFace->getTranspose();
1927 
1928  newMeshPtr->setConnectivity(outCells,outFaces,CoarseCellCoarseFace);
1929  newMeshPtr->setConnectivity(outFaces,outCells,CoarseFaceCoarseCell);
1930 
1931  CRPtr CoarseFacesFineFaces=CRPtr(new CRConnectivity(outFaces,inFaces));
1932  CoarseFacesFineFaces->initCount();
1933 
1934  VectorT3Array areaSum(outFaces.getCount());
1935  areaSum.zero();
1936 
1937  for(int f=0;f<inFaceCount;f++)
1938  {
1939  int fc0=inFaceinCells(f,0);
1940  int fc1=inFaceinCells(f,1);
1941  const int cc0=FineToCoarse[fc0];
1942  const int cc1=FineToCoarse[fc1];
1943 
1944  int cc0neibs=CellCellCoarse->getCount(cc0);
1945  int cc1neibs=CellCellCoarse->getCount(cc1);
1946 
1947  if(cc0>=outCells.getSelfCount())
1948  cc0neibs++;
1949  if(cc1>=outCells.getSelfCount())
1950  cc1neibs++;
1951 
1952  if(cc1!=cc0)
1953  {
1954  const int cfaces=CoarseCellCoarseFace->getCount(cc0);
1955 
1956  for(int cf=0;cf<cfaces;cf++)
1957  {
1958  const int face=(*CoarseCellCoarseFace)(cc0,cf);
1959  const int tempc0=(*CoarseFaceCoarseCell)(face,0);
1960  const int tempc1=(*CoarseFaceCoarseCell)(face,1);
1961 
1962  if(((cc0==tempc0)&&(cc1==tempc1))||((cc1==tempc0)&&(cc0==tempc1)))
1963  {
1964  T sign(0);
1965 
1966  if(cc1==tempc0)
1967  sign*=-1.;
1968 
1969  areaSum[face]+=inFA[f]*sign;
1970 
1971  if(cc1neibs>1 && cc0neibs>1)
1972  {
1973  CoarseFacesFineFaces->addCount(face,1);
1974  break;
1975  }
1976  else
1977  {
1978  T sumMag=sqrt(areaSum[face][0]*areaSum[face][0]+
1979  areaSum[face][1]*areaSum[face][1]+
1980  areaSum[face][2]*areaSum[face][2]);
1981 
1982  T partMag=sqrt(inFA[f][0]*inFA[f][0]+
1983  inFA[f][1]*inFA[f][1]+
1984  inFA[f][2]*inFA[f][2]);
1985 
1986  if(sumMag/partMag>.75)
1987  {
1988  CoarseFacesFineFaces->addCount(face,1);
1989  break;
1990  }
1991  else
1992  areaSum[face]-=inFA[f]*sign;
1993 
1994  }
1995  }
1996  }
1997  }
1998  }
1999 
2000  CoarseFacesFineFaces->finishCount();
2001 
2002  areaSum.zero();
2003 
2004  for(int f=0;f<inFaceCount;f++)
2005  {
2006  int fc0=inFaceinCells(f,0);
2007  int fc1=inFaceinCells(f,1);
2008  const int cc0=FineToCoarse[fc0];
2009  const int cc1=FineToCoarse[fc1];
2010 
2011  int cc0neibs=CellCellCoarse->getCount(cc0);
2012  int cc1neibs=CellCellCoarse->getCount(cc1);
2013 
2014  if(cc0>=outCells.getSelfCount())
2015  cc0neibs++;
2016  if(cc1>=outCells.getSelfCount())
2017  cc1neibs++;
2018 
2019  if(cc1!=cc0)
2020  {
2021  const int cfaces=CoarseCellCoarseFace->getCount(cc0);
2022 
2023  for(int cf=0;cf<cfaces;cf++)
2024  {
2025  const int face=(*CoarseCellCoarseFace)(cc0,cf);
2026  const int tempc0=(*CoarseFaceCoarseCell)(face,0);
2027  const int tempc1=(*CoarseFaceCoarseCell)(face,1);
2028 
2029  if(((cc0==tempc0)&&(cc1==tempc1))||((cc1==tempc0)&&(cc0==tempc1)))
2030  {
2031 
2032  T sign(1);
2033  if(cc1==tempc0)
2034  sign*=-1.;
2035 
2036  areaSum[face]+=inFA[f]*sign;
2037 
2038  if(cc1neibs>1 && cc0neibs>1)
2039  {
2040  CoarseFacesFineFaces->add(face,f);
2041  break;
2042  }
2043  else
2044  {
2045  T sumMag=sqrt(areaSum[face][0]*areaSum[face][0]+
2046  areaSum[face][1]*areaSum[face][1]+
2047  areaSum[face][2]*areaSum[face][2]);
2048 
2049  T partMag=sqrt(inFA[f][0]*inFA[f][0]+
2050  inFA[f][1]*inFA[f][1]+
2051  inFA[f][2]*inFA[f][2]);
2052 
2053  if(sumMag/partMag>.75)
2054  {
2055  CoarseFacesFineFaces->add(face,f);
2056  break;
2057  }
2058  else
2059  areaSum[face]-=inFA[f]*sign;
2060  }
2061  }
2062  }
2063  }
2064  }
2065 
2066  CoarseFacesFineFaces->finishAdd();
2067 
2068  //const StorageSite& interiorFaces=
2069  newMeshPtr->createInteriorFaceGroup(coarseInteriorCount);
2070 
2071  int inOffset=coarseInteriorCount;
2072  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2073  {
2074  const FaceGroup& fg=*fgPtr;
2075  const int faceCount=fg.site.getCount();
2076  if(fg.id>0)
2077  {
2078  if(_bcMap[fg.id]->bcType == "Interface")
2079  {
2080  IntArray& FaceFine2Coarse=*PreFacePairMap[&fg.site];
2081  int coarseFaces(0);
2082  for(int i=0;i<faceCount;i++)
2083  if(FaceFine2Coarse[i]>coarseFaces)
2084  coarseFaces=FaceFine2Coarse[i];
2085 
2086  coarseFaces+=1;
2087 
2088  //const StorageSite& newBoundarySite=
2089  newMeshPtr->createBoundaryFaceGroup(coarseFaces,inOffset,fg.id,fg.groupType);
2090  inOffset+=coarseFaces;
2091  }
2092  else
2093  {
2094  const int size=fg.site.getCount();
2095  const StorageSite& newBoundarySite=newMeshPtr->createBoundaryFaceGroup(
2096  size,inOffset,fg.id,fg.groupType);
2097  const VectorT3Array& oldFgCoords=dynamic_cast<const VectorT3Array&>(inGeomFields.coordinate[fg.site]);
2098  VT3Ptr newFgCoordsPtr=VT3Ptr(new VectorT3Array(size));
2099  (*newFgCoordsPtr)=oldFgCoords;
2100  inGeomFields.coordinate.addArray(newBoundarySite,newFgCoordsPtr);
2101  inOffset+=size;
2102  }
2103  }
2104  }
2105 
2106  // have to make coarse interface face groups -- !!
2107 
2108  const StorageSite::GatherMap& coarseGatherMap = outCells.getGatherMap();
2109 
2110  foreach(const StorageSite::GatherMap::value_type pos, coarseGatherMap)
2111  {
2112  const StorageSite& oSite = *pos.first;
2113  const Array<int>& toIndices = *pos.second;
2114 
2115  int to_where = oSite.getGatherProcID();
2116  if ( to_where != -1 )
2117  {
2118  const int cellCount=toIndices.getLength();
2119  int faceCount(0);
2120 
2121  for(int i=0;i<cellCount;i++)
2122  {
2123  const int c1=toIndices[i];
2124  faceCount+=CellCellCoarse->getCount(c1);
2125  }
2126 
2127 
2128  newMeshPtr->createInterfaceGroup(faceCount,inOffset,to_where);
2129  inOffset+=faceCount;
2130  }
2131  }
2132 
2133  //now make the geom fields
2134  const int outCellsCount=outCells.getSelfCount();
2135  const int outCellsTotCount=outCells.getCount();
2136  VT3Ptr outCellCoordsPtr=VT3Ptr(new VectorT3Array(outCellsTotCount));
2137  TArrptr outCellVolumePtr=TArrptr(new TArray(outCellsTotCount));
2138  TArray& outCV=*outCellVolumePtr;
2139  VectorT3Array& outCoords=*outCellCoordsPtr;
2140  outCV=0.;
2141 
2142  Field& VolumeField=inGeomFields.volume;
2143  const TArray& inCV=dynamic_cast<const TArray&>(VolumeField[inCells]);
2144  const VectorT3Array& inCoords=dynamic_cast<const VectorT3Array&>(inGeomFields.coordinate[inCells]);
2145 
2146  for(int c=0;c<outCellsCount;c++)
2147  {
2148  const int fineCount=CoarseToFineCells->getCount(c);
2149  VectorT3 newCoord;
2150  newCoord[0]=0.;
2151  newCoord[1]=0.;
2152  newCoord[2]=0.;
2153  for(int i=0;i<fineCount;i++)
2154  {
2155  int fc=(*CoarseToFineCells)(c,i);
2156  outCV[c]+=inCV[fc];
2157  newCoord+=inCoords[fc]*inCV[fc];
2158  }
2159  outCoords[c]=newCoord/outCV[c];
2160  }
2161 
2162  for(int c=outCellsCount;c<outCellsTotCount;c++)
2163  {
2164  const int fineCount=CoarseToFineCells->getCount(c);
2165  VectorT3 newCoord;
2166  newCoord[0]=0.;
2167  newCoord[1]=0.;
2168  newCoord[2]=0.;
2169  for(int i=0;i<fineCount;i++)
2170  {
2171  int fc=(*CoarseToFineCells)(c,i);
2172  newCoord+=inCoords[fc];
2173  }
2174  outCoords[c]=newCoord;
2175  }
2176 
2177  VolumeField.addArray(outCells,outCellVolumePtr);
2178  inGeomFields.coordinate.addArray(outCells,outCellCoordsPtr);
2179 
2180  const int outFacesCount=outFaces.getCount();
2181  VT3Ptr outFaceAreaPtr=VT3Ptr(new VectorT3Array(outFacesCount));
2182  VectorT3Array& outFA=*outFaceAreaPtr;
2183  TArrptr outFaceAreaMagPtr=TArrptr(new TArray(outFacesCount));
2184  TArray& outFAMag=*outFaceAreaMagPtr;
2185 
2186  VectorT3 myZero;
2187  myZero[0]=0.;
2188  myZero[1]=0.;
2189  myZero[2]=0.;
2190 
2191  outFA=myZero;
2192  outFAMag=0.;
2193  for(int f=0;f<outFacesCount;f++)
2194  {
2195  const int fineCount=CoarseFacesFineFaces->getCount(f);
2196  const int cCell0=(*CoarseFaceCoarseCell)(f,0);
2197  for(int i=0;i<fineCount;i++)
2198  {
2199  const int fFace=(*CoarseFacesFineFaces)(f,i);
2200  const int fCell0=inFaceinCells(fFace,0);
2201  const int CCell0=FineToCoarse[fCell0];
2202 
2203  //must make sure the area vector is pointing
2204  //from c0 to c1
2205  if(CCell0==cCell0)
2206  outFA[f]+=inFA[fFace];
2207  else
2208  outFA[f]-=inFA[fFace];
2209  outFAMag[f]+=areaMagArray[fFace];
2210  }
2211  }
2212 
2213  FaceAreaField.addArray(outFaces,outFaceAreaPtr);
2214  areaMagField.addArray(outFaces,outFaceAreaMagPtr);
2215  //cout<<"Level: "<<_level+1<<" Mesh: "<<n<<" Cells: "<<outCells.getSelfCount()<<endl;
2216 
2217  if(smallestMesh<0)
2218  smallestMesh=outCells.getSelfCount();
2219  else
2220  {
2221  if(outCells.getSelfCount()<smallestMesh)
2222  smallestMesh=outCells.getSelfCount();
2223  }
2224 
2225  /*
2226  //This is for checking purposes only
2227  cout<<"Coarse Faces to Fine Faces"<<endl;
2228  for(int f=0;f<outFaces.getCount();f++)
2229  {
2230  const int neibs=CoarseFacesFineFaces->getCount(f);
2231  for(int n=0;n<neibs;n++)
2232  cout<<f<<" "<<(*CoarseFacesFineFaces)(f,n)<<endl;
2233  cout<<endl;
2234  }
2235  cout<<"Coarse Cells to Coarse Faces"<<endl;
2236  for(int c=0;c<outCells.getCount();c++)
2237  {
2238  const int neibs=CoarseCellCoarseFace->getCount(c);
2239  for(int n=0;n<neibs;n++)
2240  cout<<c<<" "<<(*CoarseCellCoarseFace)(c,n)<<endl;
2241  cout<<endl;
2242  }
2243  */
2244  }
2245 
2246  return smallestMesh;
2247  }
2248 
2249  int coarsenInterior(const int m, const Mesh& mesh,
2250  GeomFields& inGeomFields, int& coarseCount)
2251  {
2252  const StorageSite& inCells=mesh.getCells();
2253  const StorageSite& inFaces=mesh.getFaces();
2254  const int inCellCount=inCells.getSelfCount();
2255  const CRConnectivity& inCellinFaces=mesh.getCellFaces();
2256  const CRConnectivity& inFaceinCells=mesh.getFaceCells(inFaces);
2257  Field& areaMagField=inGeomFields.areaMag;
2258  const TArray& areaMagArray=dynamic_cast<const TArray&>(areaMagField[inFaces]);
2259  const BCfaceArray& inBCfArray=*(_BFaces[m]);
2260 
2261  IntArray& FineToCoarse=dynamic_cast<IntArray&>(inGeomFields.fineToCoarse[inCells]);
2262 
2263  //first sweep to make initial pairing
2264  int pairWith;
2265  for(int c=0;c<inCellCount;c++)
2266  {
2267  if(FineToCoarse[c]<0) //dont bother if im already paired
2268  {
2269  //loop through all neighbors to find pairing
2270  const int neibCount=inCellinFaces.getCount(c);
2271  pairWith=-1;
2272  T maxArea=0.;
2273  int c2;
2274  for(int neib=0;neib<neibCount;neib++)
2275  {
2276  const int f=inCellinFaces(c,neib);
2277 
2278  if(inBCfArray[f]==0) //not a boundary face
2279  {
2280  if(c==inFaceinCells(f,1))
2281  c2=inFaceinCells(f,0);
2282  else
2283  c2=inFaceinCells(f,1);
2284 
2285  if(FineToCoarse[c2]==-1) //not already paired
2286  {
2287  if(areaMagArray[f]>maxArea)
2288  {
2289  pairWith=c2;
2290  maxArea=areaMagArray[f];
2291  }
2292  }
2293 
2294  }
2295  }
2296 
2297  if(pairWith!=-1)
2298  {
2299  FineToCoarse[c]=coarseCount;
2300  FineToCoarse[pairWith]=coarseCount;
2301  coarseCount++;
2302  }
2303  }
2304  }
2305 
2306  //second sweep to group stragglers, or group with self
2307  for(int c=0;c<inCellCount;c++)
2308  {
2309  if(FineToCoarse[c]==-1)
2310  {
2311  const int neibCount=inCellinFaces.getCount(c);
2312  T maxArea=0.;
2313  int c2,c2perm;
2314  pairWith=-2;
2315 
2316  for(int neib=0;neib<neibCount;neib++)
2317  {
2318  const int f=inCellinFaces(c,neib);
2319 
2320  if(inBCfArray[f]==0) //not a boundary face
2321  {
2322  if(c==inFaceinCells(f,1))
2323  c2=inFaceinCells(f,0);
2324  else
2325  c2=inFaceinCells(f,1);
2326 
2327  if(areaMagArray[f]>maxArea)
2328  {
2329  pairWith=FineToCoarse[c2]; //coarse level cell
2330  c2perm=c2; //fine level cell
2331  maxArea=areaMagArray[f];
2332  }
2333  }
2334  }
2335 
2336  if(pairWith==-2) //could not find suitable pair
2337  {
2338  FineToCoarse[c]=coarseCount;
2339  coarseCount++;
2340  }
2341  else //found cell to pair with
2342  {
2343  if(FineToCoarse[c2perm]==-1)
2344  {
2345  FineToCoarse[c]=coarseCount;
2346  FineToCoarse[c2perm]=coarseCount;
2347  coarseCount++;
2348  }
2349  else
2350  FineToCoarse[c]=FineToCoarse[c2perm];
2351  }
2352  }
2353  }
2354  return coarseCount;
2355  }
2356 
2357  int correctSingleNeighbor(const int m, const Mesh& mesh,
2358  GeomFields& inGeomFields, int coarseCount,
2359  map<const StorageSite*, IntArray*> PreFacePairMap)
2360  {
2361 
2362  const StorageSite& inCells=mesh.getCells();
2363  const int inCellTotal=inCells.getCount();
2364  const StorageSite& inFaces=mesh.getFaces();
2365  const int inFaceCount=inFaces.getCount();
2366  const CRConnectivity& inFaceinCells=mesh.getFaceCells(inFaces);
2367  const IntArray& FineToCoarseConst=dynamic_cast<const IntArray&>
2368  (inGeomFields.fineToCoarse[inCells]);
2369 
2370  IntArray FineToCoarse(FineToCoarseConst.getLength());
2371  FineToCoarse=FineToCoarseConst;
2372 
2373  int coarseGhost=coarseCount;
2374  int outGhost(0);
2375  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2376  {
2377  const FaceGroup& fg=*fgPtr;
2378 
2379  if(fg.id>0)
2380  {
2381  const CRConnectivity& faceCells=mesh.getFaceCells(fg.site);
2382  const int faceCount=fg.site.getCount();
2383 
2384  if(_bcMap[fg.id]->bcType == "Interface")
2385  {
2386  IntArray& FaceFine2Coarse=*PreFacePairMap[&fg.site];
2387 
2388  int coarseFaces(0);
2389  for(int i=0;i<faceCount;i++)
2390  {
2391  const int cghost=faceCells(i,1);
2392  FineToCoarse[cghost]=FaceFine2Coarse[i]+coarseGhost;
2393  if(FaceFine2Coarse[i]>coarseFaces)
2394  coarseFaces=FaceFine2Coarse[i];
2395  }
2396  coarseGhost+=coarseFaces+1;
2397  outGhost+=coarseFaces+1;
2398  }
2399  else if(_bcMap[fg.id]->bcType == "temperature" ||
2400  _bcMap[fg.id]->bcType == "reflecting" )
2401  {
2402  for(int i=0;i<faceCount;i++)
2403  {
2404  const int cghost=faceCells(i,1);
2405  FineToCoarse[cghost]=coarseGhost;
2406  coarseGhost++;
2407  outGhost++;
2408  }
2409  }
2410 
2411  }
2412  }
2413 
2414  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2415  {
2416  const FaceGroup& fg=*fgPtr;
2417  const CRConnectivity& faceCells=mesh.getFaceCells(fg.site);
2418  const int faceCount=fg.site.getCount();
2419 
2420  for(int i=0;i<faceCount;i++)
2421  {
2422  const int cghost=faceCells(i,1);
2423  FineToCoarse[cghost]=coarseGhost;
2424  coarseGhost++;
2425  outGhost++;
2426  }
2427 
2428  }
2429 
2430  StorageSite preOutCells(coarseCount,outGhost);
2431  CRPtr CoarseToFineCells=CRPtr(new CRConnectivity(preOutCells,inCells));
2432  CoarseToFineCells->initCount();
2433 
2434  for(int c=0;c<inCellTotal;c++)
2435  CoarseToFineCells->addCount(FineToCoarse[c],1);
2436 
2437  CoarseToFineCells->finishCount();
2438 
2439  for(int c=0;c<inCellTotal;c++)
2440  CoarseToFineCells->add(FineToCoarse[c],c);
2441 
2442  CoarseToFineCells->finishAdd();
2443 
2444  CRPtr FineFacesCoarseCells=CRPtr(new CRConnectivity(inFaces,preOutCells));
2445  FineFacesCoarseCells->initCount();
2446 
2447  //count surviving faces
2448  int survivingFaces=0;
2449  int coarse0, coarse1;
2450  for(int f=0;f<inFaceCount;f++)
2451  {
2452  coarse0=FineToCoarse[inFaceinCells(f,0)];
2453  coarse1=FineToCoarse[inFaceinCells(f,1)];
2454  if(coarse0!=coarse1)
2455  {
2456  survivingFaces++;
2457  FineFacesCoarseCells->addCount(f,2);
2458  }
2459  }
2460 
2461  FineFacesCoarseCells->finishCount();
2462 
2463  //make non-zero's
2464  for(int f=0;f<inFaceCount;f++)
2465  {
2466  int fc0=inFaceinCells(f,0);
2467  int fc1=inFaceinCells(f,1);
2468  int cc0=FineToCoarse[fc0];
2469  int cc1=FineToCoarse[fc1];
2470  if(cc0!=cc1)
2471  {
2472  FineFacesCoarseCells->add(f,cc0);
2473  FineFacesCoarseCells->add(f,cc1);
2474  }
2475  }
2476 
2477  FineFacesCoarseCells->finishAdd();
2478 
2479  CRPtr CoarseCellsFineFaces=FineFacesCoarseCells->getTranspose();
2480  CRPtr CellCellCoarse=CoarseCellsFineFaces->multiply(*FineFacesCoarseCells,true);
2481 
2482 
2483  for(int c=0;c<preOutCells.getSelfCount();c++)
2484  {
2485  if(CellCellCoarse->getCount(c)<2)
2486  {
2487  /*
2488  const int bigCoarse=(*CellCellCoarse)(c,0);
2489  const int fineCount=CoarseToFineCells->getCount(c);
2490 
2491  for(int j=0;j<fineCount;j++)
2492  {
2493  const int fc=(*CoarseToFineCells)(c,j);
2494  FineToCoarse[fc]=bigCoarse;
2495  }
2496  coarseCount-=1;
2497  */
2498 
2499  const int removal=(*CoarseToFineCells)(c,0);
2500  FineToCoarse[removal]=coarseCount;
2501  coarseCount++;
2502 
2503  }
2504  }
2505 
2506  return coarseCount;
2507  }
2508 
2509  void coarsenInterfaceCells(COMETIC<T>& ic, IntArray& coarseCounts,
2510  GeomFields& inGeomFields, const MeshList& inMeshes)
2511  {
2512  const int Mid0=ic.MeshID0;
2513  const int Mid1=ic.MeshID1;
2514  const int Fid0=ic.FgID0;
2515  const int Fid1=ic.FgID1;
2516  const Mesh& mesh0=*inMeshes[Mid0];
2517  const Mesh& mesh1=*inMeshes[Mid1];
2518  const StorageSite& faces0=mesh0.getFaceGroup(Fid0).site;
2519  const StorageSite& faces1=mesh1.getFaceGroup(Fid1).site;
2520  const StorageSite& cells0=mesh0.getCells();
2521  const StorageSite& cells1=mesh1.getCells();
2522  const StorageSite& globFaces0=mesh0.getFaces();
2523  const StorageSite::CommonMap& ComMap = faces0.getCommonMap();
2524  const IntArray& common01=*(ComMap.find(&faces1)->second);
2525  const CRConnectivity& faceCells0=mesh0.getFaceCells(faces0);
2526  const CRConnectivity& faceCells1=mesh1.getFaceCells(faces1);
2527  const CRConnectivity& cellCells0=mesh0.getCellCells();
2528  const CRConnectivity& cellCells1=mesh1.getCellCells();
2529  const CRConnectivity& globFaceCells0=mesh0.getAllFaceCells();
2530  const CRConnectivity& globCellFaces0=mesh0.getCellFaces();
2531  const int inFaceCount=faces0.getCount();
2532  Field& areaMagField=inGeomFields.areaMag;
2533  const TArray& areaMagArray0=dynamic_cast<const TArray&>(areaMagField[globFaces0]);
2534  const BCfaceArray& inBCfArray0=*(_BFaces[Mid0]);
2535 
2536  const CRPtr cellFaces0Ptr=faceCells0.getTranspose();
2537  const CRConnectivity& cellFaces0=*cellFaces0Ptr;
2538  const CRPtr cellsIntGhost0Ptr=cellFaces0.multiply(faceCells0,true);
2539  const CRConnectivity& cellsIntGhost0=*cellsIntGhost0Ptr;
2540 
2541  IntArray CoarseFineCount(cells0.getCount());
2542  CoarseFineCount.zero();
2543 
2544  IntArray& f2c0=dynamic_cast<IntArray&>(inGeomFields.fineToCoarse[cells0]);
2545  IntArray& f2c1=dynamic_cast<IntArray&>(inGeomFields.fineToCoarse[cells1]);
2546  int& count0=coarseCounts[Mid0];
2547  int& count1=coarseCounts[Mid1];
2548 
2549  //const VectorT3Array& cpos0=dynamic_cast<const VectorT3Array&>(inGeomFields.coordinate[cells0]);
2550  //const VectorT3Array& cpos1=dynamic_cast<const VectorT3Array&>(inGeomFields.coordinate[cells1]);
2551 
2552  //agglomerate mesh0 first
2553  int pairWith;
2554 
2555  if(_level>=0)
2556  {
2557  for(int f=0;f<inFaceCount;f++)
2558  {
2559  const int c=faceCells0(f,0);
2560  const int cCoarse=CoarseFineCount[c];
2561  if(f2c0[c]<0 || cCoarse<4) //dont bother if im already paired
2562  {
2563  //loop through all neighbors to find pairing
2564  const int neibCount=cellCells0.getCount(c);
2565  pairWith=-2;
2566  T maxArea=0.;
2567  int c2;
2568  for(int neib=0;neib<neibCount;neib++)
2569  {
2570  const int fglob=globCellFaces0(c,neib);
2571 
2572  if(inBCfArray0[fglob]==0) //not a boundary face
2573  {
2574  if(c==globFaceCells0(fglob,1))
2575  c2=globFaceCells0(fglob,0);
2576  else
2577  c2=globFaceCells0(fglob,1);
2578 
2579  int c2Coarse=f2c0[c2];
2580 
2581  if(_level==0)
2582  {
2583  if(c2Coarse==-1)// || CoarseFineCount[c2Coarse]<5) //not already paired
2584  {
2585  if(areaMagArray0[fglob]>maxArea)
2586  {
2587  pairWith=c2;
2588  maxArea=areaMagArray0[fglob];
2589  if(cellsIntGhost0.getCount(c2)>0) //pick first interface cell
2590  break;
2591  }
2592  }
2593  }
2594  else
2595  {
2596  if(c2Coarse==-1)
2597  {
2598  if(areaMagArray0[fglob]>maxArea)
2599  {
2600  pairWith=c2;
2601  maxArea=areaMagArray0[fglob];
2602  if(cellsIntGhost0.getCount(c2)>0) //pick first interface cell
2603  break;
2604  }
2605  }
2606  }
2607 
2608  }
2609  }
2610 
2611  if(pairWith>=0)
2612  {
2613  int c2Coarse=f2c0[pairWith];
2614  if(c2Coarse==-1 && f2c0[c]==-1)
2615  {
2616  f2c0[pairWith]=count0;
2617  f2c0[c]=count0;
2618  CoarseFineCount[count0]++;
2619  count0++;
2620  }
2621  else if(f2c0[c]==-1 && c2Coarse>=0)
2622  f2c0[c]=c2Coarse;/*
2623  else if(f2c0[c]!=-1 && c2Coarse==-1)
2624  f2c0[pairWith]=f2c0[c];
2625  */
2626  }
2627  }
2628  }
2629  }
2630  else
2631  {
2632  for(int f=0;f<inFaceCount;f++)
2633  {
2634  const int c=faceCells0(f,0);
2635  if(f2c0[c]<0) //dont bother if im already paired
2636  {
2637  //loop through all neighbors to find pairing
2638  const int neibCount=cellCells0.getCount(c);
2639  pairWith=-1;
2640  T maxArea=0.;
2641  int c2;
2642  for(int neib=0;neib<neibCount;neib++)
2643  {
2644  const int fglob=globCellFaces0(c,neib);
2645 
2646  if(inBCfArray0[fglob]==0) //not a boundary face
2647  {
2648  if(c==globFaceCells0(fglob,1))
2649  c2=globFaceCells0(fglob,0);
2650  else
2651  c2=globFaceCells0(fglob,1);
2652 
2653  if(f2c0[c2]==-1) //not already paired
2654  {
2655  if(areaMagArray0[fglob]>maxArea)
2656  {
2657  pairWith=c2;
2658  maxArea=areaMagArray0[fglob];
2659  }
2660  }
2661 
2662  }
2663  }
2664 
2665  if(pairWith!=-1)
2666  {
2667  f2c0[pairWith]=count0;
2668  f2c0[c]=count0;
2669  count0++;
2670  }
2671  }
2672  }
2673  }
2674 
2675  StorageSite preCoarse0(count0);
2676  CRPtr preCoarseToFineCells0=CRPtr(new CRConnectivity(preCoarse0,cells0));
2677  preCoarseToFineCells0->initCount();
2678 
2679  for(int c=0;c<cells0.getSelfCount();c++)
2680  if(f2c0[c]!=-1)
2681  preCoarseToFineCells0->addCount(f2c0[c],1);
2682 
2683 
2684  preCoarseToFineCells0->finishCount();
2685 
2686  for(int c=0;c<cells0.getSelfCount();c++)
2687  if(f2c0[c]!=-1)
2688  preCoarseToFineCells0->add(f2c0[c],c);
2689 
2690  preCoarseToFineCells0->finishAdd();
2691 
2692  //now use mesh0 pairing to pair mesh1
2693  for(int f=0;f<inFaceCount;f++)
2694  {
2695  const int f1=common01[f];
2696  //const int c01ghost=faceCells1(f1,1);
2697  const int c01=faceCells1(f1,0);
2698  const int c00=faceCells0(f,0);
2699  if(f2c0[c00]>-1) //dont bother if im already paired, mesh0 must be paired
2700  {
2701  //const int neibCount0=cellCells0.getCount(c00);
2702  const int coarseC00=f2c0[c00];
2703  const int cC00fineCnt=preCoarseToFineCells0->getCount(coarseC00);
2704 
2705  for(int cC00fine=0;cC00fine<cC00fineCnt;cC00fine++) //who is mesh0 paired with
2706  {
2707  //const int c10=cellCells0(c00,neib0);
2708  const int fineCell=(*preCoarseToFineCells0)(coarseC00,cC00fine);
2709  if(fineCell!=c00)//(f2c0[c00]==f2c0[c10])
2710  {
2711  const int c10gNeibs=cellsIntGhost0.getCount(fineCell); //cell is on an interface
2712  for(int c10ghost=0;c10ghost<c10gNeibs;c10ghost++)
2713  {
2714  const int c10g=cellsIntGhost0(fineCell,c10ghost);
2715  const int f10=cellFaces0(c10g,0);
2716  const int f11=common01[f10];
2717  //const int c11ghost=faceCells1(f11,1);
2718  const int c11=faceCells1(f11,0);
2719  const int neibCount1=cellCells1.getCount(c01);
2720  bool isNeib(false);
2721  for(int c01neibs=0;c01neibs<neibCount1;c01neibs++)
2722  {
2723  if(cellCells1(c01,c01neibs)==c11)
2724  {
2725  isNeib=true;
2726  break;
2727  }
2728  }
2729 
2730  if(isNeib) //direct neighbor
2731  {
2732  if(f2c1[c11]>-1 && f2c1[c01]<0) //other guy paired, im not
2733  f2c1[c01]=f2c1[c11];
2734  else if(f2c1[c01]>-1 && f2c1[c11]<0) //im paired, other guy's not
2735  f2c1[c11]=f2c1[c01];
2736  else if(f2c1[c11]<0 && f2c1[c01]<0) //both not paired
2737  {
2738  f2c1[c11]=count1;
2739  f2c1[c01]=count1;
2740  count1++;
2741  }
2742  }
2743  else if(_level<-1) //check for indirect neighbors
2744  {
2745  const int c11neibCnt=cellCells1.getCount(c11);
2746  int middleMan=-1;
2747  for(int c11nb=0;c11nb<c11neibCnt;c11nb++)
2748  {
2749  const int c11neib=cellCells1(c11,c11nb);
2750  for(int c01neibs=0;c01neibs<neibCount1;c01neibs++)
2751  {
2752  if(cellCells1(c01,c01neibs)==c11neib)
2753  {
2754  middleMan=c11neib;
2755  break;
2756  }
2757  }
2758  }
2759  if(middleMan!=-1)
2760  {
2761  if(f2c1[middleMan]==-1)
2762  {
2763  if(f2c1[c11]<0 && f2c1[c01]<0) //both not paired
2764  {
2765  f2c1[c11]=count1;
2766  f2c1[c01]=count1;
2767  f2c1[middleMan]=count1;
2768  count1++;
2769  }
2770  else if(f2c1[c11]>-1 && f2c1[c01]<0) //other guy paired, im not
2771  {
2772  f2c1[c01]=f2c1[c11];
2773  f2c1[middleMan]=f2c1[c11];
2774  }
2775  else if(f2c1[c01]>-1 && f2c1[c11]<0) //im paired, other guy's not
2776  {
2777  f2c1[c11]=f2c1[c01];
2778  f2c1[middleMan]=f2c1[c01];
2779  }
2780  }
2781  }
2782  }
2783 
2784  }
2785  }
2786  }
2787 
2788  /*
2789  if(f2c1[c01]==-1) //still hasnt been paired -- means mesh0 cell wasnt paired.
2790  {
2791  f2c1[c01]=count1;
2792  count1++;
2793  }*/
2794 
2795  }
2796  }
2797 
2798  }
2799 
2800  void doSweeps(const int sweeps)
2801  {
2802  for(int sweepNo=0;sweepNo<sweeps;sweepNo++)
2803  {
2804  smooth(1);
2805  smooth(-1);
2806  }
2807  //applyTemperatureBoundaries();
2808  }
2809 
2810  void smooth(int dir)
2811  {
2812 
2813  const int numMeshes=_meshes.size();
2814  int start;
2815  if(dir==1)
2816  start=0;
2817  else
2818  start=numMeshes-1;
2819 
2820  COMETInterface<T> ComInt(_meshes,_kspaces,_MeshKspaceMap,_macro,_geomFields);
2821 
2822  for(int msh=start;((msh<numMeshes)&&(msh>-1));msh+=dir)
2823  {
2824  const Mesh& mesh=*_meshes[msh];
2825  Tkspace& kspace=*_kspaces[_MeshKspaceMap[msh]];
2826  FgTKClistMap FgToKsc;
2827 
2828  if(!_MeshToIC.empty())
2829  {
2830  IntArray& ICs=*_MeshToIC[msh];
2831  const int totIC=ICs.getLength();
2832  for(int i=0;i<totIC;i++)
2833  {
2834  COMETIC<T>& ic=*_IClist[ICs[i]];
2835  int fgid=ic.getSelfFaceID(msh);
2836  FgToKsc[fgid]=ic.getKConnectivity(fgid);
2837  }
2838  }
2839 
2840  const BCcellArray& BCArray=*(_BCells[msh]);
2841  const BCfaceArray& BCfArray=*(_BFaces[msh]);
2842 
2843  COMETDiscretizer<T> CDisc(mesh,_geomFields,_macro,
2844  kspace,_bcMap,BCArray,BCfArray,_options,
2845  FgToKsc);
2846 
2847  CDisc.setfgFinder();
2848 
2849 #ifdef FVM_PARALLEL
2850  swapGhostInfo();
2851 #endif
2852 
2853  if(_options.Scattering=="Full")
2854  {
2855  applyTemperatureBoundaries();
2856  if(_level==0)
2857  {
2858  CDisc.COMETSolveFull(dir,_level);
2859  //CDisc.COMETSolveFull(-dir,_level);
2860  }
2861  else
2862  {
2863  CDisc.COMETSolveCoarse(dir,_level);
2864  CDisc.COMETSolveCoarse(-dir,_level);
2865  }
2866  }
2867  else
2868  {
2869  if(_level>0 || _options.Convection=="FirstOrder")
2870  {
2871  applyTemperatureBoundaries();
2872  CDisc.COMETSolveCoarse(dir,_level);
2873  CDisc.COMETSolveCoarse(-dir,_level);
2874  }
2875  else
2876  {
2877  //applyTemperatureBoundaries();
2878  CDisc.COMETSolveFine(dir,_level);
2879  CDisc.COMETSolveFine(-dir,_level);
2880  }
2881  }
2882 
2883  if(!_MeshToIC.empty())
2884  {
2885  IntArray& ICs=*_MeshToIC[msh];
2886  const int totIC=ICs.getLength();
2887  for(int i=0;i<totIC;i++)
2888  {
2889  COMETIC<T>& ic=*_IClist[ICs[i]];
2890  ComInt.updateOtherGhost(ic,msh,_level!=0);
2891  }
2892  }
2893 
2894  }
2895 
2896  }
2897 
2898  T updateResid(const bool addFAS)
2899  {
2900  const int numMeshes=_meshes.size();
2901 
2902  COMETInterface<T> ComInt(_meshes,_kspaces,_MeshKspaceMap,_macro,_geomFields);
2903 
2904  const int listSize=_IClist.size();
2905  for(int ic=0;ic<listSize;ic++)
2906  {
2907  COMETIC<T>* icPtr=_IClist[ic];
2908  //if(icPtr->InterfaceModel=="DMM" && _level==0)
2909  // ComInt.makeDMMcoeffs(*icPtr);
2910  ComInt.updateResid(*icPtr,addFAS);
2911  }
2912 
2913  T highResid=-1.;
2914  T currentResid;
2915 
2916 #ifdef FVM_PARALLEL
2917  swapGhostInfo();
2918 #endif
2919 
2920  for(int msh=0;msh<numMeshes;msh++)
2921  {
2922  const Mesh& mesh=*_meshes[msh];
2923  Tkspace& kspace=*_kspaces[_MeshKspaceMap[msh]];
2924  const BCcellArray& BCArray=*(_BCells[msh]);
2925  const BCfaceArray& BCfArray=*(_BFaces[msh]);
2926  FgTKClistMap FgToKsc;
2927  COMETDiscretizer<T> CDisc(mesh,_geomFields,_macro,
2928  kspace,_bcMap,BCArray,BCfArray,_options,
2929  FgToKsc);
2930 
2931  CDisc.setfgFinder();
2932  if(_options.Scattering=="Full")
2933  {
2934  if(_level==0)
2935  CDisc.findResidFull();
2936  else
2937  CDisc.findResidCoarse(addFAS);
2938  }
2939  else
2940  {
2941  if(_level>0 || _options.Convection=="FirstOrder")
2942  CDisc.findResidCoarse(addFAS);
2943  else
2944  CDisc.findResidFine();
2945  }
2946 
2947  currentResid=CDisc.getAveResid();
2948 
2949  if(highResid<0)
2950  highResid=currentResid;
2951  else
2952  if(currentResid>highResid)
2953  highResid=currentResid;
2954  }
2955 
2956  return highResid;
2957  }
2958 
2959  void swapGhostInfo()
2960  {
2961  const int numMeshes=_meshes.size();
2962  for(int msh=0;msh<numMeshes;msh++)
2963  {
2964  const Mesh& mesh=*_meshes[msh];
2965  const StorageSite& cells=mesh.getCells();
2966  Tkspace& kspace=*_kspaces[_MeshKspaceMap[msh]];
2967  kspace.syncLocal(cells);
2968  }
2969  }
2970 
2971  void cycle()
2972  {
2973  doSweeps(_options.preSweeps);
2974 
2975  if(_level+1<_options.maxLevels)
2976  {
2977  updateResid(_level!=0);
2978  injectResid();
2979  _coarserLevel->sete0();
2980  _coarserLevel->makeFAS();
2981  _coarserLevel->cycle();
2982  correctSolution();
2983  }
2984 
2985  doSweeps(_options.postSweeps);
2986  }
2987 
2988  void injectResid()
2989  {
2990  const int numMeshes = _meshes.size();
2991 
2992  for (int n=0; n<numMeshes; n++)
2993  {
2994  const Mesh& finerMesh=*_meshes[n];
2995  int Knum=_MeshKspaceMap[n];
2996  Tkspace& finerKspace=*_kspaces[Knum];
2997  const Mesh& coarserMesh=*(_coarserLevel->getMeshList())[n];
2998  Tkspace& coarserKspace=_coarserLevel->getKspace(Knum);
2999  PhononMacro& coarserMacro=_coarserLevel->getMacro();
3000  const StorageSite& finerCells=finerMesh.getCells();
3001  const StorageSite& coarserCells=coarserMesh.getCells();
3002  const CRConnectivity& CoarserToFiner=coarserMesh.getConnectivity(coarserCells,finerCells);
3003  const TArray& coarserVol=dynamic_cast<TArray&>(_geomFields.volume[coarserCells]);
3004  const TArray& finerVol=dynamic_cast<TArray&>(_geomFields.volume[finerCells]);
3005 
3006  TArray& coarserVar=coarserKspace.geteArray();
3007  TArray& coarserInj=coarserKspace.getInjArray();
3008  TArray& coarserFAS=coarserKspace.getFASArray();
3009  TArray& finerVar=finerKspace.geteArray();
3010  TArray& finerRes=finerKspace.getResArray();
3011 
3012  const int cellCount=coarserCells.getSelfCount();
3013  const int cellTotCount=coarserCells.getCount();
3014  coarserVar.zero();
3015  coarserFAS.zero();
3016 
3017  for(int c=0;c<cellCount;c++)
3018  {
3019  const int fineCount=CoarserToFiner.getCount(c);
3020 
3021  //sum up contributions from fine cells
3022  for(int fc=0;fc<fineCount;fc++)
3023  {
3024  const int cell=CoarserToFiner(c,fc);
3025  const int klen=finerKspace.getlength();
3026  int coarserCellIndex=coarserKspace.getGlobalIndex(c,0);
3027  int finerCellIndex=finerKspace.getGlobalIndex(cell,0);
3028 
3029  for(int k=0;k<klen;k++)
3030  {
3031  Tkvol& kvol=finerKspace.getkvol(k);
3032  const int numModes=kvol.getmodenum();
3033  for(int m=0;m<numModes;m++)
3034  {
3035  coarserVar[coarserCellIndex]+=finerVar[finerCellIndex]*finerVol[cell];
3036  coarserFAS[coarserCellIndex]+=finerRes[finerCellIndex];
3037  finerCellIndex++;
3038  coarserCellIndex++;
3039  }
3040  }
3041  }
3042 
3043  //make volume average
3044  int coarserCellIndex=coarserKspace.getGlobalIndex(c,0);
3045  const int klen=finerKspace.getlength();
3046  for(int k=0;k<klen;k++)
3047  {
3048  Tkvol& kvol=finerKspace.getkvol(k);
3049  //Tkvol& ckvol=coarserKspace.getkvol(k);
3050  const int numModes=kvol.getmodenum();
3051  for(int m=0;m<numModes;m++)
3052  {
3053  coarserVar[coarserCellIndex]/=coarserVol[c];
3054  coarserInj[coarserCellIndex]=coarserVar[coarserCellIndex];
3055  coarserCellIndex++;
3056  }
3057  }
3058  }
3059 
3060  for(int c=cellCount;c<cellTotCount;c++)
3061  {
3062  const int cell=CoarserToFiner(c,0);
3063  int coarserCellIndex=coarserKspace.getGlobalIndex(c,0);
3064  int finerCellIndex=finerKspace.getGlobalIndex(cell,0);
3065  const int klen=finerKspace.getlength();
3066  for(int k=0;k<klen;k++)
3067  {
3068  Tkvol& kvol=finerKspace.getkvol(k);
3069  const int numModes=kvol.getmodenum();
3070  for(int m=0;m<numModes;m++)
3071  {
3072  coarserVar[coarserCellIndex]=0.;
3073  coarserFAS[coarserCellIndex]=0.;
3074  coarserVar[coarserCellIndex]+=finerVar[finerCellIndex];
3075  coarserFAS[coarserCellIndex]+=finerRes[finerCellIndex];
3076  coarserInj[coarserCellIndex]=coarserVar[coarserCellIndex];
3077  coarserCellIndex++;
3078  finerCellIndex++;
3079  }
3080  }
3081  }
3082 
3083  TArray& coarserVarM=dynamic_cast<TArray&>(coarserMacro.temperature[coarserCells]);
3084  TArray& coarserInjM=dynamic_cast<TArray&>(coarserMacro.TlInjected[coarserCells]);
3085  TArray& coarserFASM=dynamic_cast<TArray&>(coarserMacro.TlFASCorrection[coarserCells]);
3086  TArray& finerVarM=dynamic_cast<TArray&>(_macro.temperature[finerCells]);
3087  TArray& finerResM=dynamic_cast<TArray&>(_macro.TlResidual[finerCells]);
3088 
3089  for(int c=0;c<cellCount;c++)
3090  {
3091  const int fineCount=CoarserToFiner.getCount(c);
3092  coarserVarM[c]=0.;
3093  coarserFASM[c]=0.;
3094 
3095  for(int fc=0;fc<fineCount;fc++)
3096  {
3097  const int cell=CoarserToFiner(c,fc);
3098  coarserVarM[c]+=finerVarM[cell]*finerVol[cell];
3099  coarserFASM[c]+=finerResM[cell];
3100  }
3101  coarserVarM[c]/=coarserVol[c];
3102  coarserInjM[c]=coarserVarM[c];
3103  }
3104  }
3105  }
3106 
3107  void makeFAS()
3108  {
3109  updateResid(false);
3110 
3111  const int numMeshes = _meshes.size();
3112  for (int n=0; n<numMeshes; n++)
3113  {
3114  const Mesh& mesh=*_meshes[n];
3115  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3116  const StorageSite& cells=mesh.getCells();
3117  kspace.makeFAS();
3118 
3119  TArray& resArray=dynamic_cast<TArray&>(_macro.TlResidual[cells]);
3120  TArray& fasArray=dynamic_cast<TArray&>(_macro.TlFASCorrection[cells]);
3121  fasArray-=resArray;
3122  }
3123  }
3124 
3125  void sete0()
3126  {
3127  const int numMeshes = _meshes.size();
3128  for (int n=0; n<numMeshes; n++)
3129  {
3130  const Mesh& mesh=*_meshes[n];
3131  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3132  const StorageSite& cells=mesh.getCells();
3133  const int cellCount=cells.getSelfCount();
3134  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[cells]);
3135  TArray& e0Array=kspace.gete0Array();
3136  const int klen=kspace.getlength();
3137 
3138  for(int c=0;c<cellCount;c++)
3139  {
3140  const T Tlat=Tl[c];
3141  int cellIndex=kspace.getGlobalIndex(c,0);
3142  for(int k=0;k<klen;k++)
3143  {
3144  Tkvol& kvol=kspace.getkvol(k);
3145  const int numModes=kvol.getmodenum();
3146  for(int m=0;m<numModes;m++)
3147  {
3148  Tmode& mode=kvol.getmode(m);
3149  e0Array[cellIndex]=mode.calce0(Tlat);
3150  kspace.updateTau(c,Tl[c]);
3151  cellIndex++;
3152  }
3153  }
3154  }
3155  }
3156  }
3157 
3158  void correctSolution()
3159  {
3160  const int numMeshes = _meshes.size();
3161 
3162  for (int n=0; n<numMeshes; n++)
3163  {
3164  const Mesh& finerMesh=*_meshes[n];
3165  int Knum=_MeshKspaceMap[n];
3166  Tkspace& finerKspace=*_kspaces[Knum];
3167  const Mesh& coarserMesh=*(_coarserLevel->getMeshList())[n];
3168  Tkspace& coarserKspace=_coarserLevel->getKspace(Knum);
3169  PhononMacro& coarserMacro=_coarserLevel->getMacro();
3170  const StorageSite& finerCells=finerMesh.getCells();
3171  const StorageSite& coarserCells=coarserMesh.getCells();
3172  const CRConnectivity& CoarserToFiner=coarserMesh.getConnectivity(coarserCells,finerCells);
3173  const CRConnectivity& coarseCellFaces=coarserMesh.getCellFaces();
3174  BCfaceArray& BCfArray=*(_BFaces[n]);
3175 
3176  TArray& coarserArray=coarserKspace.geteArray();
3177  TArray& finerArray=finerKspace.geteArray();
3178  TArray& injArray=coarserKspace.getInjArray();
3179 
3180  const int cellCount=coarserCells.getSelfCount();
3181  const int cellTotCount=coarserCells.getCount();
3182 
3183  for(int c=0;c<cellCount;c++)
3184  {
3185  const int fineCount=CoarserToFiner.getCount(c);
3186  for(int fc=0;fc<fineCount;fc++)
3187  {
3188  int coarserCellIndex=coarserKspace.getGlobalIndex(c,0);
3189  const int finerCell=CoarserToFiner(c,fc);
3190  int finerCellIndex=finerKspace.getGlobalIndex(finerCell,0);
3191  const int klen=finerKspace.getlength();
3192  for(int k=0;k<klen;k++)
3193  {
3194  Tkvol& kvol=finerKspace.getkvol(k);
3195  const int numModes=kvol.getmodenum();
3196  for(int m=0;m<numModes;m++)
3197  {
3198  const T correction=coarserArray[coarserCellIndex]-injArray[coarserCellIndex];
3199  finerArray[finerCellIndex]+=correction;
3200  coarserCellIndex++;
3201  finerCellIndex++;
3202  }
3203  }
3204  }
3205  }
3206 
3207  for(int c=cellCount;c<cellTotCount;c++)
3208  {
3209  const int f=coarseCellFaces(c,0);
3210  if(BCfArray[f]==4) //correct interface cells only
3211  {
3212  const int fineCount=CoarserToFiner.getCount(c);
3213  for(int fc=0;fc<fineCount;fc++)
3214  {
3215  const int finerCell=CoarserToFiner(c,fc);
3216  int coarserCellIndex=coarserKspace.getGlobalIndex(c,0);
3217  int finerCellIndex=finerKspace.getGlobalIndex(finerCell,0);
3218  const int klen=finerKspace.getlength();
3219 
3220  for(int k=0;k<klen;k++)
3221  {
3222  Tkvol& kvol=finerKspace.getkvol(k);
3223  const int numModes=kvol.getmodenum();
3224  for(int m=0;m<numModes;m++)
3225  {
3226  const T correction=coarserArray[coarserCellIndex]-injArray[coarserCellIndex];
3227  finerArray[finerCellIndex]+=correction;
3228  coarserCellIndex++;
3229  finerCellIndex++;
3230  }
3231  }
3232 
3233  }
3234  }
3235  }
3236 
3237  TArray& coarserArrayM=dynamic_cast<TArray&>(coarserMacro.temperature[coarserCells]);
3238  TArray& injArrayM=dynamic_cast<TArray&>(coarserMacro.TlInjected[coarserCells]);
3239  TArray& finerArrayM=dynamic_cast<TArray&>(_macro.temperature[finerCells]);
3240 
3241  for(int c=0;c<cellCount;c++)
3242  {
3243  const int fineCount=CoarserToFiner.getCount(c);
3244  const T correction=coarserArrayM[c]-injArrayM[c];
3245 
3246  for(int fc=0;fc<fineCount;fc++)
3247  finerArrayM[CoarserToFiner(c,fc)]+=correction;
3248  }
3249  }
3250  }
3251 
3252  void updateTL()
3253  {
3254  const int numMeshes=_meshes.size();
3255  for(int n=0;n<numMeshes;n++)
3256  {
3257  const Mesh& mesh=*_meshes[n];
3258  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3259  const StorageSite& cells = mesh.getCells();
3260  const int numcells = cells.getCount();
3261  const T tauTot=kspace.calcTauTot();
3262  TArray& TL=dynamic_cast<TArray&>(_macro.temperature[cells]);
3263  TArray& e0Array=dynamic_cast<TArray&>(_macro.e0[cells]);
3264 
3265  for(int c=0;c<numcells;c++)
3266  kspace.NewtonSolve(TL[c],e0Array[c]*tauTot);
3267  }
3268  }
3269 
3270  void advance(const int iters)
3271  {
3272  _residual=updateResid(false);
3273  int niters=0;
3274  const T absTol=_options.absTolerance;
3275  const int show=_options.showResidual;
3276 
3277  while((niters<iters) && (_residual>absTol))
3278  {
3279  cycle();
3280  niters++;
3281  _residual=updateResid(false);
3282  if(niters%show==0)
3283  cout<<"Iteration:"<<niters<<" Residual:"<<_residual<<endl;
3284 
3285  }
3286 
3287  //applyTemperatureBoundaries();
3288  //calcModeTemps();
3289  }
3290 
3291  T HeatFluxIntegral(const Mesh& mesh, const int faceGroupId)
3292  {
3293  T r(0.);
3294  bool found = false;
3295  const int n=mesh.getID();
3296  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3297  TArray& eArray=kspace.geteArray();
3298  const T DK3=kspace.getDK3();
3299 
3300  const T hbar=6.582119e-16; // (eV s)
3301 
3302  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3303  {
3304  const FaceGroup& fg = *fgPtr;
3305  if (fg.id == faceGroupId)
3306  {
3307  const StorageSite& faces = fg.site;
3308  const int nFaces = faces.getCount();
3309  //const StorageSite& cells = mesh.getCells();
3310  const CRConnectivity& faceCells=mesh.getFaceCells(faces);
3311  const Field& areaField=_geomFields.area;
3312  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
3313 
3314  for(int f=0; f<nFaces; f++)
3315  {
3316  const VectorT3 An=faceArea[f];
3317  const int c1=faceCells(f,1);
3318  int cellIndex=kspace.getGlobalIndex(c1,0);
3319  for(int k=0;k<kspace.getlength();k++)
3320  {
3321  Tkvol& kv=kspace.getkvol(k);
3322  int modenum=kv.getmodenum();
3323  for(int m=0;m<modenum;m++)
3324  {
3325  VectorT3 vg=kv.getmode(m).getv();
3326  T dk3=kv.getdk3();
3327  //T energy=hbar*kv.getmode(m).getomega();
3328  const T vgdotAn=An[0]*vg[0]+An[1]*vg[1]+An[2]*vg[2];
3329  r += eArray[cellIndex]*vgdotAn*(dk3/DK3);//*energy;
3330  cellIndex++;
3331  }
3332  }
3333  }
3334  found=true;
3335  }
3336  }
3337 
3338 #ifdef FVM_PARALLEL
3339  found=true;
3340  int one=1;
3341  double tempR=r;
3342  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempR, one, MPI::DOUBLE, MPI::SUM);
3343  r=tempR;
3344 #endif
3345 
3346  if (!found)
3347  throw CException("getHeatFluxIntegral: invalid faceGroupID");
3348  return r*DK3;
3349  }
3350 
3351  T HeatFluxIntegralFace(const Mesh& mesh, const int f)
3352  {
3353  T r(0.);
3354  //bool found = false;
3355  const int n=mesh.getID();
3356  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3357  TArray& eArray=kspace.geteArray();
3358  const T DK3=kspace.getDK3();
3359 
3360  const T hbar=6.582119e-16; // (eV s)
3361  //const FaceGroupPtr fgPtr=mesh.getInteriorFaceGroup();
3362 
3363  const FaceGroup& fg = mesh.getInteriorFaceGroup();
3364 
3365  const StorageSite& faces = fg.site;
3366  //const int nFaces = faces.getCount();
3367  //const StorageSite& cells = mesh.getCells();
3368  const CRConnectivity& faceCells=mesh.getFaceCells(faces);
3369  const Field& areaField=_geomFields.area;
3370  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
3371 
3372 
3373  const VectorT3 An=faceArea[f];
3374  const int c1=faceCells(f,1);
3375  const int c0=faceCells(f,0);
3376  int cellIndex=kspace.getGlobalIndex(c1,0);
3377  int cellIndex0=kspace.getGlobalIndex(c0,0);
3378  for(int k=0;k<kspace.getlength();k++)
3379  {
3380  Tkvol& kv=kspace.getkvol(k);
3381  int modenum=kv.getmodenum();
3382  for(int m=0;m<modenum;m++)
3383  {
3384  VectorT3 vg=kv.getmode(m).getv();
3385  T dk3=kv.getdk3();
3386  //T energy=hbar*kv.getmode(m).getomega();
3387  const T vgdotAn=An[0]*vg[0]+An[1]*vg[1]+An[2]*vg[2];
3388  if(vgdotAn>0)
3389  r += eArray[cellIndex0]*vgdotAn*(dk3/DK3);//*energy;
3390  else
3391  r += eArray[cellIndex]*vgdotAn*(dk3/DK3);
3392  cellIndex++;
3393  cellIndex0++;
3394  }
3395  }
3396 
3397  return r*DK3;
3398  }
3399 
3400  ArrayBase* binwiseHeatFluxIntegral(const Mesh& mesh, const int faceGroupId)
3401  {
3402  bool found = false;
3403  const int n=mesh.getID();
3404  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3405  DensityOfStates<T>& dos=*kspace.getDOSptr();
3406  TArray& freqMids=dos.getFreqMidsT();
3407  const T binNos=freqMids.getLength();
3408  TArray& eArray=kspace.geteArray();
3409  TArray* qptr(new TArray(binNos));
3410  TArray& q(*qptr);
3411  q.zero();
3412 
3413  //const T hbar=6.582119e-16; // (eV s)
3414 
3415  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3416  {
3417  const FaceGroup& fg = *fgPtr;
3418  if (fg.id == faceGroupId)
3419  {
3420  const StorageSite& faces = fg.site;
3421  const int nFaces = faces.getCount();
3422  //const StorageSite& cells = mesh.getCells();
3423  const CRConnectivity& faceCells=mesh.getFaceCells(faces);
3424  const Field& areaField=_geomFields.area;
3425  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
3426 
3427  for(int f=0; f<nFaces; f++)
3428  {
3429  const VectorT3 An=faceArea[f];
3430  const int c1=faceCells(f,1);
3431 
3432  for(int binIndx=0;binIndx<binNos;binIndx++)
3433  {
3434  IntArray& kpts=dos.getKIndices(binIndx);
3435  IntArray& mpts=dos.getMIndices(binIndx);
3436 
3437  for(int k=0;k<kpts.getLength();k++)
3438  {
3439  Tkvol& kv=kspace.getkvol(kpts[k]);
3440  T dk3=kv.getdk3();
3441  Tmode& mode=kv.getmode(mpts[k]);
3442  VectorT3 vg=mode.getv();
3443  T vgdotAn=vg[0]*An[0]+vg[1]*An[1]+vg[2]*An[2];
3444  const int index=mode.getIndex()-1;
3445  int cellIndex=kspace.getGlobalIndex(c1,index);
3446  q[binIndx]+= eArray[cellIndex]*vgdotAn*dk3;//*energy;
3447  cellIndex++;
3448  }
3449  }
3450  }
3451  found=true;
3452  }
3453  }
3454 
3455  if (!found)
3456  throw CException("getHeatFluxIntegral: invalid faceGroupID");
3457  return qptr;
3458  }
3459 
3460  ArrayBase* modewiseHeatFluxIntegral(const Mesh& mesh, const int faceGroupId)
3461  {
3462  //bool found = false;
3463  const int n=mesh.getID();
3464  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3465  TArray& eArray=kspace.geteArray();
3466  //const T DK3=kspace.getDK3();
3467  TArray* qptr(new TArray(kspace.gettotmodes()));
3468  TArray& q(*qptr);
3469  q.zero();
3470 
3471  const T hbar=6.582119e-16; // (eV s)
3472 
3473  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3474  {
3475  const FaceGroup& fg = *fgPtr;
3476  if (fg.id == faceGroupId)
3477  {
3478  const StorageSite& faces = fg.site;
3479  const int nFaces = faces.getCount();
3480  //const StorageSite& cells = mesh.getCells();
3481  const CRConnectivity& faceCells=mesh.getFaceCells(faces);
3482  const Field& areaField=_geomFields.area;
3483  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
3484 
3485  for(int f=0; f<nFaces; f++)
3486  {
3487  const VectorT3 An=faceArea[f];
3488  const int c1=faceCells(f,1);
3489  int cellIndex=kspace.getGlobalIndex(c1,0);
3490  for(int k=0;k<kspace.getlength();k++)
3491  {
3492  Tkvol& kv=kspace.getkvol(k);
3493  int modenum=kv.getmodenum();
3494  for(int m=0;m<modenum;m++)
3495  {
3496  VectorT3 vg=kv.getmode(m).getv();
3497  const int index=kv.getmode(m).getIndex()-1;
3498  T dk3=kv.getdk3();
3499  //T energy=hbar*kv.getmode(m).getomega();
3500  const T vgdotAn=An[0]*vg[0]+An[1]*vg[1]+An[2]*vg[2];
3501  q[index]+= eArray[cellIndex]*vgdotAn*dk3;//*energy;
3502  cellIndex++;
3503  }
3504  }
3505  }
3506 
3507  }
3508  }
3509 
3510  return qptr;
3511  }
3512 
3513  T getWallArea(const Mesh& mesh, const int faceGroupId)
3514  {
3515  T r(0.);
3516  bool found = false;
3517  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3518  {
3519  const FaceGroup& fg = *fgPtr;
3520  if (fg.id == faceGroupId)
3521  {
3522  const StorageSite& faces = fg.site;
3523  const int nFaces = faces.getCount();
3524  const Field& areaMagField=_geomFields.areaMag;
3525  const TArray& faceArea=dynamic_cast<const TArray&>(areaMagField[faces]);
3526 
3527  for(int f=0; f<nFaces; f++)
3528  r += faceArea[f];
3529 
3530  found=true;
3531  break;
3532  }
3533  }
3534 
3535  if (!found)
3536  throw CException("getwallArea: invalid faceGroupID");
3537  return r;
3538  }
3539 
3540  VectorT3 getWallAreaVector(const Mesh& mesh, const int faceGroupId)
3541  {
3542  VectorT3 An;
3543  An=0.;
3544  bool found = false;
3545  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3546  {
3547  const FaceGroup& fg = *fgPtr;
3548  if (fg.id == faceGroupId)
3549  {
3550  const StorageSite& faces = fg.site;
3551  const int nFaces = faces.getCount();
3552  const Field& areaField=_geomFields.area;
3553  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
3554  for(int f=0; f<nFaces; f++)
3555  An+=faceArea[f];
3556  found=true;
3557  break;
3558  }
3559  }
3560 
3561 #ifdef FVM_PARALLEL
3562  found=true;
3563  int one=1;
3564  double tempR=An[0];
3565  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempR, one, MPI::DOUBLE, MPI::SUM);
3566  An[0]=tempR;
3567  tempR=An[1];
3568  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempR, one, MPI::DOUBLE, MPI::SUM);
3569  An[1]=tempR;
3570  tempR=An[2];
3571  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempR, one, MPI::DOUBLE, MPI::SUM);
3572  An[2]=tempR;
3573 #endif
3574 
3575  if (!found)
3576  throw CException("getwallArea: invalid faceGroupID");
3577  return An;
3578  }
3579 
3580  ArrayBase* getValueArray(const Mesh& mesh, const int cell)
3581  {
3582  //only returns the e" values, not the lattice temperature
3583  const int n=mesh.getID();
3584  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3585  const int allModes=kspace.gettotmodes();
3586  TArray* vals=new TArray(allModes);
3587  const StorageSite& cells=mesh.getCells();
3588  const int len=kspace.getlength();
3589  int count=0;
3590  for(int k=0;k<len;k++)
3591  {
3592  Tkvol& kvol=kspace.getkvol(k);
3593  const int modes=kvol.getmodenum();
3594  for(int m=0;m<modes;m++)
3595  {
3596  Tmode& mode=kvol.getmode(m);
3597  Field& eField=mode.getfield();
3598  const TArray& eArray=dynamic_cast<const TArray&>(eField[cells]);
3599  (*vals)[count]=eArray[cell];
3600  count++;
3601  }
3602  }
3603  return vals;
3604  }
3605 
3606  void equilibrate()
3607  {
3608  const int numMeshes = _meshes.size();
3609  for (int n=0; n<numMeshes; n++)
3610  {
3611  const Mesh& mesh=*_meshes[n];
3612  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3613  const StorageSite& cells=mesh.getCells();
3614  const int numcells=cells.getCount();
3615  const TArray& Tl=dynamic_cast<const TArray&>(_macro.temperature[cells]);
3616  const int len=kspace.getlength();
3617 
3618  for(int k=0;k<len;k++)
3619  {
3620  Tkvol& kvol=kspace.getkvol(k);
3621  const int modes=kvol.getmodenum();
3622  for(int m=0;m<modes;m++)
3623  {
3624  Tmode& mode=kvol.getmode(m);
3625  Field& eField=mode.getfield();
3626  TArray& eArray=dynamic_cast<TArray&>(eField[cells]);
3627  for(int c=0;c<numcells;c++)
3628  eArray[c]=mode.calce0(Tl[c]);
3629  }
3630  }
3631  }
3632  }
3633 
3634  ArrayBase& getLatticeTemp(const Mesh& mesh)
3635  {
3636  const StorageSite& cells=mesh.getCells();
3637  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[cells]);
3638  return Tl;
3639  }
3640 
3641  void setStraightLine(const T T1, const T T2)
3642  {
3643  const int numMeshes = _meshes.size();
3644 
3645  T xmax(0.);
3646  T xmin(0.);
3647 
3648  for (int n=0; n<numMeshes; n++)
3649  {
3650  const Mesh& mesh=*_meshes[n];
3651  const StorageSite& cells=mesh.getCells();
3652  const int numcells=cells.getCount();
3653  const VectorT3Array& coords=dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
3654 
3655  for(int c=0;c<numcells;c++)
3656  {
3657  T x=coords[c][0];
3658  if(x>xmax)
3659  xmax=x;
3660  if(x<xmin)
3661  xmin=x;
3662  }
3663 
3664  }
3665 
3666  for (int n=0; n<numMeshes; n++)
3667  {
3668  const Mesh& mesh=*_meshes[n];
3669  const StorageSite& cells=mesh.getCells();
3670  //const int numcells=cells.getCount();
3671  const int selfcells=cells.getSelfCount();
3672  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[cells]);
3673  const VectorT3Array& coords=dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
3674 
3675  for(int c=0;c<selfcells;c++)
3676  {
3677  T factor=(coords[c][0]-xmin)/(xmax-xmin);
3678  Tl[c]=T1+factor*(T2-T1);
3679  }
3680 
3681  }
3682 
3683  }
3684 
3685  void calcModeTemps()
3686  {
3687  const int numMeshes=_meshes.size();
3688  for (int n=0;n<numMeshes;n++)
3689  {
3690  Mesh& mesh=*_meshes[n];
3691  if(_MeshKspaceMap[n]==-1)
3692  throw CException("Have not set the Kspace for this Mesh!!");
3693  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3694  const int numK=kspace.getlength();
3695  const StorageSite& cells=mesh.getCells();
3696  const int numcells=cells.getCount();
3697  FieldVector& FieldVec=*_macro.BranchTemperatures[n];
3698  const int modeCount=FieldVec.size();
3699  TArray eSum(numcells);
3700  const TArray& TL=dynamic_cast<const TArray&>(_macro.temperature[cells]);
3701 
3702  for(int m=0;m<modeCount;m++)
3703  {
3704  eSum.zero();
3705  Field& modeField=*FieldVec[m];
3706  TArray& modeTemp=dynamic_cast<TArray&>(modeField[cells]);
3707  for(int c=0;c<numcells;c++)
3708  {
3709  for(int k=0;k<numK;k++)
3710  {
3711  T dk3=kspace.getkvol(k).getdk3();
3712  Tmode& mode=kspace.getkvol(k).getmode(m);
3713  const int index=mode.getIndex()-1;
3714  eSum[c]+=kspace.gete(c,index)*dk3;
3715  }
3716  }
3717  for(int c=0;c<numcells;c++)
3718  modeTemp[c]=kspace.calcModeTemp(TL[c],eSum[c],m);
3719  }
3720 
3721  }
3722  }
3723 
3724  void calcBandTemps()
3725  {
3726  const int numMeshes=_meshes.size();
3727  for (int n=0;n<numMeshes;n++)
3728  {
3729  Mesh& mesh=*_meshes[n];
3730  if(_MeshKspaceMap[n]==-1)
3731  throw CException("Have not set the Kspace for this Mesh!!");
3732  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3733  DensityOfStates<T>& dos=*kspace.getDOSptr();
3734  const StorageSite& cells=mesh.getCells();
3735  const TArray& TL=dynamic_cast<const TArray&>(_macro.temperature[cells]);
3736  const int numcells=cells.getCount();
3737  TArray eSum(numcells);
3738  const int bandCount=dos.getFreqMidsT().getLength();
3739 
3740  FieldVector* FieldVecPtr=new FieldVector();
3741  _macro.BandTemperatures[n]=FieldVecPtr;
3742  for(int b=0;b<bandCount;b++)
3743  {
3744  eSum.zero();
3745  shared_ptr<Field> bandField(new Field("bandTemp"));
3746  TArrptr Tptr(new TArray(numcells));
3747  bandField->addArray(cells,Tptr);
3748  FieldVecPtr->push_back(bandField);
3749  TArray& bandTemp=*Tptr;
3750  bandTemp.zero();
3751  IntArray& kpts=dos.getKIndices(b);
3752  IntArray& mpts=dos.getMIndices(b);
3753  for(int c=0;c<numcells;c++)
3754  {
3755  for(int i=0;i<kpts.getLength();i++)
3756  {
3757  const int k=kpts[i];
3758  const int m=mpts[i];
3759  T dk3=kspace.getkvol(k).getdk3();
3760  Tmode& mode=kspace.getkvol(k).getmode(m);
3761  const int index=mode.getIndex()-1;
3762  eSum[c]+=kspace.gete(c,index)*dk3;
3763  }
3764  }
3765  for(int c=0;c<numcells;c++)
3766  bandTemp[c]=kspace.calcBandTemp(TL[c],eSum[c],kpts,mpts);
3767  }
3768 
3769  }
3770  }
3771 
3772  void calcBandRelEnergy()
3773  {
3774  const int numMeshes=_meshes.size();
3775  for (int n=0;n<numMeshes;n++)
3776  {
3777  Mesh& mesh=*_meshes[n];
3778  if(_MeshKspaceMap[n]==-1)
3779  throw CException("Have not set the Kspace for this Mesh!!");
3780  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3781  const T DK3=kspace.getDK3();
3782  DensityOfStates<T>& dos=*kspace.getDOSptr();
3783  const StorageSite& cells=mesh.getCells();
3784  const int numcells=cells.getCount();
3785  TArray eSum(numcells);
3786  eSum.zero();
3787  const int bandCount=dos.getFreqMidsT().getLength();
3788  const int numK=kspace.getlength();
3789  FieldVector* FieldVecPtr=new FieldVector();
3790  _macro.BandRelEnergy[n]=FieldVecPtr;
3791 
3792  for(int c=0;c<numcells;c++)
3793  {
3794  int index(0);
3795  for(int k=0;k<numK;k++)
3796  {
3797  Tkvol& kv=kspace.getkvol(k);
3798  T dk3=kv.getdk3();
3799  const int numM=kv.getmodenum();
3800  for(int m=0;m<numM;m++)
3801  {
3802  eSum[c]+=kspace.gete(c,index)*(dk3/DK3);
3803  index++;
3804  }
3805  }
3806  }
3807 
3808  for(int b=0;b<bandCount;b++)
3809  {
3810  shared_ptr<Field> bandField(new Field("bandRelEnergy"));
3811  TArrptr Tptr(new TArray(numcells));
3812  bandField->addArray(cells,Tptr);
3813  FieldVecPtr->push_back(bandField);
3814  TArray& bandEn=*Tptr;
3815  bandEn.zero();
3816  IntArray& kpts=dos.getKIndices(b);
3817  IntArray& mpts=dos.getMIndices(b);
3818  for(int c=0;c<numcells;c++)
3819  {
3820  T cellSum(0);
3821  for(int i=0;i<kpts.getLength();i++)
3822  {
3823  const int k=kpts[i];
3824  const int m=mpts[i];
3825  T dk3=kspace.getkvol(k).getdk3();
3826  Tmode& mode=kspace.getkvol(k).getmode(m);
3827  const int index=mode.getIndex()-1;
3828  cellSum+=kspace.gete(c,index)*(dk3/DK3);
3829  }
3830  bandEn[c]=cellSum/eSum[c];
3831  }
3832  }
3833 
3834  }
3835  }
3836 
3837  void calcModeFlux()
3838  {
3839  const int numMeshes=_meshes.size();
3840  const T eVtoJoule=1.60217646e-19;
3841  for (int n=0;n<numMeshes;n++)
3842  {
3843  Mesh& mesh=*_meshes[n];
3844  if(_MeshKspaceMap[n]==-1)
3845  throw CException("Have not set the Kspace for this Mesh!!");
3846  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3847  const int numK=kspace.getlength();
3848  const StorageSite& cells=mesh.getCells();
3849  const int numcells=cells.getCount();
3850  const int modeCount=kspace.getkvol(0).getmodenum();
3851  const T DK3=kspace.getDK3();
3852  VectorT3Array& Q=dynamic_cast<VectorT3Array&>(_macro.heatFlux[cells]);
3853  Q.zero();
3854 
3855  FieldVector* FieldVecPtr=new FieldVector();
3856  _macro.BranchFlux[n]=FieldVecPtr;
3857  for(int m=0;m<modeCount;m++)
3858  {
3859  shared_ptr<Field> modeField(new Field("mode"));
3860  VT3Ptr qptr(new VectorT3Array(numcells));
3861  VectorT3Array& q(*qptr);
3862  q.zero();
3863  modeField->addArray(cells,qptr);
3864  FieldVecPtr->push_back(modeField);
3865 
3866  for(int c=0;c<numcells;c++)
3867  {
3868  for(int k=0;k<numK;k++)
3869  {
3870  Tmode& mode=kspace.getkvol(k).getmode(m);
3871  T dk3=kspace.getkvol(k).getdk3();
3872  const int index=mode.getIndex()-1;
3873  q[c]+=mode.getv()*(dk3/DK3)*kspace.gete(c,index)*eVtoJoule;
3874  Q[c]+=mode.getv()*(dk3/DK3)*kspace.gete(c,index)*eVtoJoule;
3875  }
3876 
3877  }
3878  for(int i=0;i<q.getLength();i++)
3879  q[i]*=DK3;
3880  }
3881  for(int i=0;i<Q.getLength();i++)
3882  Q[i]*=DK3;
3883  }
3884  }
3885 
3886  void calcBandFlux()
3887  {
3888  const int numMeshes=_meshes.size();
3889  const T eVtoJoule=1.60217646e-19;
3890  for (int n=0;n<numMeshes;n++)
3891  {
3892  Mesh& mesh=*_meshes[n];
3893  if(_MeshKspaceMap[n]==-1)
3894  throw CException("Have not set the Kspace for this Mesh!!");
3895  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
3896  //const int numK=kspace.getlength();
3897  const StorageSite& cells=mesh.getCells();
3898  const int numcells=cells.getCount();
3899  //const int modeCount=kspace.getkvol(0).getmodenum();
3900  const T DK3=kspace.getDK3();
3901  DensityOfStates<T>& dos=*kspace.getDOSptr();
3902  const int bandCount=dos.getFreqMidsT().getLength();
3903 
3904  FieldVector* FieldVecPtr=new FieldVector();
3905  _macro.BandFlux[n]=FieldVecPtr;
3906  for(int b=0;b<bandCount;b++)
3907  {
3908  shared_ptr<Field> bandField(new Field("BandFlux"));
3909  VT3Ptr qptr(new VectorT3Array(numcells));
3910  VectorT3Array& q(*qptr);
3911  q.zero();
3912  bandField->addArray(cells,qptr);
3913  FieldVecPtr->push_back(bandField);
3914  IntArray& kpts=dos.getKIndices(b);
3915  IntArray& mpts=dos.getMIndices(b);
3916  for(int c=0;c<numcells;c++)
3917  {
3918  for(int i=0;i<kpts.getLength();i++)
3919  {
3920  const int k=kpts[i];
3921  const int m=mpts[i];
3922  Tmode& mode=kspace.getkvol(k).getmode(m);
3923  T dk3=kspace.getkvol(k).getdk3();
3924  const int index=mode.getIndex()-1;
3925  q[c]+=mode.getv()*(dk3/DK3)*kspace.gete(c,index)*eVtoJoule;
3926  }
3927 
3928  }
3929  for(int i=0;i<q.getLength();i++)
3930  q[i]*=DK3;
3931  }
3932  }
3933  }
3934 
3935  T calcDomainStats()
3936  {
3937  const int numMeshes=_meshes.size();
3938  TArray meshVols(numMeshes);
3939  T totalVol(0.);
3940  meshVols.zero();
3941  T interfaceArea(0.);
3942  T interfaceAreaX(0.);
3943  T interfaceAreaY(0.);
3944  T interfaceAreaZ(0.);
3945  int interfaceFaces(0);
3946 
3947  for (int n=0;n<numMeshes;n++)
3948  {
3949  Mesh& mesh=*_meshes[n];
3950  const StorageSite& cells=mesh.getCells();
3951  const TArray& vol=dynamic_cast<const TArray&>(_geomFields.volume[cells]);
3952 
3953  for(int c=0;c<cells.getSelfCount();c++)
3954  meshVols[n]+=vol[c];
3955 
3956  totalVol+=meshVols[n];
3957 
3958  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3959  {
3960  const FaceGroup& fg = *fgPtr;
3961  if(fg.id>0)
3962  {
3963  if(_bcMap[fg.id]->bcType == "Interface")
3964  {
3965  const StorageSite& faces=fg.site;
3966  const int numFaces=faces.getCount();
3967  const Field& AreaMagField=_geomFields.areaMag;
3968  const TArray& AreaMag=
3969  dynamic_cast<const TArray&>(AreaMagField[faces]);
3970  const VectorT3Array& AreaDir=
3971  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
3972 
3973  interfaceFaces+=numFaces;
3974 
3975  for(int f=0;f<numFaces;f++)
3976  {
3977  interfaceArea+=AreaMag[f];
3978  interfaceAreaX+=fabs(AreaDir[f][0]);
3979  interfaceAreaY+=fabs(AreaDir[f][1]);
3980  interfaceAreaZ+=fabs(AreaDir[f][2]);
3981  }
3982 
3983  }
3984  }
3985  }
3986  }
3987 
3988  interfaceArea/=2.;
3989  interfaceAreaX/=2.;
3990  interfaceAreaY/=2.;
3991  interfaceAreaZ/=2.;
3992  interfaceFaces/=2;
3993 
3994  for (int n=0;n<numMeshes;n++)
3995  {
3996  Mesh& mesh=*_meshes[n];
3997  const StorageSite& cells=mesh.getCells();
3998  cout<<endl;
3999  cout<<"Mesh: "<<n<<endl;
4000  cout<<"Self Cell Count: "<<cells.getSelfCount()<<endl;
4001  cout<<"Ghost Cell Count: "<<cells.getCount()-cells.getSelfCount()<<endl;
4002  cout<<"Mesh Volume: "<<meshVols[n]<<endl;
4003  }
4004 
4005  cout<<endl;
4006  cout<<"Total Volume: "<<totalVol<<endl;
4007  cout<<"Interface Area: "<<interfaceArea<<endl;
4008  cout<<"Interface Area (X): "<<interfaceAreaX<<endl;
4009  cout<<"Interface Area (Y): "<<interfaceAreaY<<endl;
4010  cout<<"Interface Area (Z): "<<interfaceAreaZ<<endl;
4011  cout<<"Area/Volume: "<<interfaceArea/totalVol<<endl;
4012  cout<<"Volume/Area: "<<totalVol/interfaceArea<<endl;
4013  cout<<"Interface Face Count: "<<interfaceFaces<<endl;
4014  return interfaceArea/totalVol;
4015 
4016  }
4017 
4018  void makeCellColors(const int level)
4019  {
4020 
4021  TCOMET* thismodel=getModelPointer(level);
4022  MeshList& meshes=const_cast<MeshList&>(thismodel->getMeshList());
4023  const int numMeshes = meshes.size();
4024  shared_ptr<Field> colorFieldPtr(new Field("Colors"));
4025  _macro.CellColors.push_back(colorFieldPtr);
4026  Field& colorField=*colorFieldPtr;
4027  for (int n=0; n<numMeshes; n++)
4028  {
4029  const Mesh& mesh=*meshes[n];
4030  const StorageSite& cells=mesh.getCells();
4031  const int cellCount=cells.getSelfCount();
4032  const int cellTotCount=cells.getCount();
4033  const CRConnectivity& cellCells=mesh.getCellCells();
4034 
4035  TArrptr colorPtr(new TArray(cellTotCount));
4036  colorField.addArray(cells, colorPtr);
4037  TArray& colorArray=(*colorPtr);
4038  colorArray=-1;
4039 
4040  for(int c=0;c<cellCount;c++)
4041  {
4042  T color=0.;
4043  const int neibs=cellCells.getCount(c);
4044 
4045  while(colorArray[c]==-1)
4046  {
4047 
4048  bool same(false);
4049 
4050  for(int j=0;j<neibs;j++)
4051  {
4052  if(color==colorArray[cellCells(c,j)])
4053  {
4054  same=true;
4055  break;
4056  }
4057  }
4058 
4059  if(same)
4060  color++;
4061  else
4062  colorArray[c]=color;
4063 
4064  }
4065 
4066  }
4067  }
4068 
4069  //if(_coarserLevel!=NULL)
4070  //_coarserLevel->makeCellColors();
4071 
4072  }
4073 
4074  void makePlotColors(const int level)
4075  {
4076 
4077  if(level==0)
4078  {
4079  if(_options.maxLevels>1)
4080  makeFinestToCoarseConn();
4081  shared_ptr<Field> field0ptr(new Field("plotColor"));
4082  _macro.plottingCellColors.push_back(field0ptr);
4083  Field& field0=*field0ptr;
4084  Field& colorField=_macro.getColorField(0);
4085 
4086  const int numMeshes = _meshes.size();
4087  for (int n=0; n<numMeshes; n++)
4088  {
4089  const Mesh& mesh=*_meshes[n];
4090  const StorageSite& cells=mesh.getCells();
4091  //const int cellCount=cells.getSelfCount();
4092  const int cellTotCount=cells.getCount();
4093  //const CRConnectivity& cellCells=mesh.getCellCells();
4094  TArray& colorArray=dynamic_cast<TArray&>(colorField[cells]);
4095 
4096  TArrptr colorPtr(new TArray(cellTotCount));
4097  (*colorPtr)=colorArray;
4098  field0.addArray(cells, colorPtr);
4099  }
4100  //if(_coarserLevel!=NULL)
4101  // _coarserLevel->makePlotColors();
4102  }
4103  else
4104  {
4105  shared_ptr<Field> fieldptr(new Field("plotColor"));
4106  _macro.plottingCellColors.push_back(fieldptr);
4107  Field& field=*fieldptr;
4108 
4109  Field& colorField=_macro.getColorField(level);
4110 
4111  TCOMET* coarseModel=getModelPointer(level);
4112  const MeshList& coarseMeshes=coarseModel->getMeshList();
4113  const int numMeshes = _meshes.size();
4114  for(int n=0;n<numMeshes;n++)
4115  {
4116  const Mesh& finestMesh=*_meshes[n];
4117  const Mesh& coarseMesh=*coarseMeshes[n];
4118  const StorageSite& coarseCells=coarseMesh.getCells();
4119  const StorageSite& finestCells=finestMesh.getCells();
4120  const int finestCellTotCount=finestCells.getCount();
4121  const int finestCellCount=finestCells.getSelfCount();
4122  const CRConnectivity& finestToCoarse=finestMesh.getConnectivity(finestCells,coarseCells);
4123 
4124  TArrptr plotColorPtr(new TArray(finestCellTotCount));
4125  TArray& plotColorArray=*plotColorPtr;
4126  field.addArray(finestCells, plotColorPtr);
4127 
4128  TArray& colorArray=dynamic_cast<TArray&>(colorField[coarseCells]);
4129 
4130  for(int c=0;c<finestCellCount;c++)
4131  plotColorArray[c]=colorArray[finestToCoarse(c,0)];
4132  }
4133 
4134  //if(_coarserLevel!=NULL)
4135  // _coarserLevel->makePlotColors();
4136  }
4137  }
4138 
4139  void makeFinestToCoarseConn()
4140  {
4141  const int numMeshes = _meshes.size();
4142  for (int n=0; n<numMeshes; n++)
4143  {
4144  Mesh& mesh=*_meshes[n];
4145  const StorageSite& cells=mesh.getCells();
4146 
4147  Mesh& mesh1=*(_coarserLevel->getMeshList()[n]);
4148  const StorageSite& cells1=mesh1.getCells();
4149 
4150  const CRConnectivity& firstCRC=mesh1.getConnectivity(cells1,cells);
4151  CRPtr FineToCoarseOld=firstCRC.getTranspose();
4152 
4153  mesh.setConnectivity(cells,cells1,FineToCoarseOld);
4154 
4155  for(int lvl=2;lvl<_options.maxLevels;lvl++)
4156  {
4157  TCOMET* coarseModel=getModelPointer(lvl);
4158  Mesh& coarseMesh=*(coarseModel->getMeshList()[n]);
4159  const StorageSite& coarseCells=coarseMesh.getCells();
4160  TCOMET* finerModel=coarseModel->getFinerModel();
4161  Mesh& finerMesh=*(finerModel->getMeshList()[n]);
4162  const StorageSite& finerCells=finerMesh.getCells();
4163  const CRConnectivity& coarseToFine=coarseMesh.getConnectivity(coarseCells, finerCells);
4164  CRPtr FineToCoarseNew=coarseToFine.getTranspose();
4165  FineToCoarseNew=FineToCoarseOld->multiply(*FineToCoarseNew,true);
4166  mesh.setConnectivity(cells, coarseCells,FineToCoarseNew);
4167  FineToCoarseOld=FineToCoarseNew;
4168  }
4169 
4170  }
4171  }
4172 
4173  TCOMET* getModelPointer(const int level)
4174  {
4175  if(_level==level)
4176  return this;
4177  else
4178  return _coarserLevel->getModelPointer(level);
4179  }
4180 
4181  bool sameFaceGroup(const Mesh& mesh, const int f0, const int f1)
4182  {
4183  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
4184  {
4185  const FaceGroup& fg=*fgPtr;
4186  const int off=fg.site.getOffset();
4187  const int cnt=fg.site.getCount();
4188  const int fin=off+cnt;
4189 
4190  if(f0<fin && f0>=off)
4191  if(f1<fin && f1>=off)
4192  return true;
4193  }
4194  return false;
4195  }
4196 
4197  T getAverageTemperature()
4198  {
4199  const int numMeshes = _meshes.size();
4200  T r(0);
4201  T volTot(0);
4202  for (int n=0; n<numMeshes; n++)
4203  {
4204  const Mesh& mesh=*_meshes[n];
4205  const StorageSite& cells=mesh.getCells();
4206  const int numcells=cells.getSelfCount();
4207  const TArray& Tl=dynamic_cast<const TArray&>(_macro.temperature[cells]);
4208  const TArray& cv=dynamic_cast<const TArray&>(_geomFields.volume[cells]);
4209  for(int c=0;c<numcells;c++)
4210  {
4211  r+=Tl[c]*cv[c];
4212  volTot+=cv[c];
4213  }
4214  }
4215  return r/volTot;
4216  }
4217 
4218  void makeNonEqTemp()
4219  {
4220  const int numMeshes = _meshes.size();
4221  for (int n=0; n<numMeshes; n++)
4222  {
4223  const Mesh& mesh=*_meshes[n];
4224  const StorageSite& cells=mesh.getCells();
4225  const int numcells=cells.getSelfCount();
4226  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[cells]);
4227  Tkspace& kspace=*_kspaces[_MeshKspaceMap[n]];
4228  //const TArray& e=kspace.geteArray();
4229 
4230  //const int len=kspace.getlength();
4231 
4232  for(int c=0;c<numcells;c++)
4233  Tl[c]=kspace.calcLatTemp(c);
4234 
4235  }
4236  }
4237 
4238  void setBCMap(COMETBCMap* bcMap) {_bcMap=*bcMap;}
4239  void setCoarserLevel(TCOMET* cl) {_coarserLevel=cl;}
4240  void setFinerLevel(TCOMET* fl) {_finerLevel=fl;}
4241  TCOMET* getFinerModel() {return _finerLevel;}
4242  int getLevel() {return _level;}
4243  const MeshList& getMeshList() {return _meshes;}
4244  GeomFields& getGeomFields() {return _geomFields;}
4245  TkspList& getKspaces() {return _kspaces;}
4246  Tkspace& getKspace(const int i) {return *_kspaces[i];}
4247  PhononMacro& getMacro() {return _macro;}
4248  T getResidual()
4249  {
4250 #ifdef FVM_PARALLEL
4251  int one=1;
4252  double tempResid=_residual;
4253  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempResid, one, MPI::DOUBLE, MPI::SUM);
4254  _residual=tempResid;
4255 #endif
4256  return _residual;
4257  }
4258  IClist& getIClist() {return _IClist;}
4259  MeshICmap& getMeshICmap() {return _MeshToIC;}
4260 
4261  private:
4262 
4263  const int _level; //0 being the finest level
4264  GeomFields& _geomFields;
4265  TkspList& _kspaces;
4266  PhononMacro& _macro;
4267  TCOMET* _finestLevel;
4268  TCOMET* _coarserLevel;
4269  TCOMET* _finerLevel;
4270  COMETBCMap _bcMap;
4271  TCModOpts _options;
4272  BCfaceList _BFaces;
4273  BCcellList _BCells;
4274  T _residual;
4275  MeshKspaceMap _MeshKspaceMap;
4276  IClist _IClist;
4277  MeshICmap _MeshToIC;
4278  int _rank;
4279 
4280 };
4281 
4282 #endif
COMETModel::modewiseHeatFluxIntegral
ArrayBase * modewiseHeatFluxIntegral(const Mesh &mesh, const int faceGroupId)
Definition: phononbase/COMETModel.h:3460
COMETModel::TKConnectivity
KSConnectivity< T > TKConnectivity
Definition: phononbase/COMETModel.h:65
COMETModel::getKspace
Tkspace & getKspace(const int i)
Definition: phononbase/COMETModel.h:4246
pmode::Reflptr
shared_ptr< Reflection > Reflptr
Definition: pmode.h:28
Mesh::getBoundaryFaceGroups
const FaceGroupList & getBoundaryFaceGroups() const
Definition: Mesh.h:187
pmode::Refl_pair
pair< Reflection, Reflection > Refl_pair
Definition: pmode.h:29
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
Mesh::getInteriorFaceGroup
const FaceGroup & getInteriorFaceGroup() const
Definition: Mesh.h:181
COMETModel::getMacro
PhononMacro & getMacro()
Definition: phononbase/COMETModel.h:4247
COMETModel::getLatticeTemp
ArrayBase & getLatticeTemp(const Mesh &mesh)
Definition: phononbase/COMETModel.h:3634
Mesh::setID
void setID(const int id)
Definition: Mesh.h:321
COMETModel::BCcellArray
Array< int > BCcellArray
Definition: phononbase/COMETModel.h:61
COMETModel::getOptions
TCModOpts & getOptions()
Definition: phononbase/COMETModel.h:126
Kspace::setFASArray
void setFASArray(TArrPtr FASPtr)
Definition: Kspace.h:1222
Kspace::NewtonSolve
void NewtonSolve(T &guess, const T e_sum)
Definition: Kspace.h:425
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
COMETModel::_geomFields
GeomFields & _geomFields
Definition: esbgkbase/COMETModel.h:5652
COMETModel::updateResid
T updateResid(const bool addFAS)
Definition: esbgkbase/COMETModel.h:5038
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
Kspace::syncLocal
void syncLocal(const StorageSite &site)
Definition: Kspace.h:1291
COMETModel::MeshICmap
vector< IntArrPtr > MeshICmap
Definition: phononbase/COMETModel.h:64
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
COMETModel::binwiseHeatFluxIntegral
ArrayBase * binwiseHeatFluxIntegral(const Mesh &mesh, const int faceGroupId)
Definition: phononbase/COMETModel.h:3400
pmode::getv
Tvec getv()
Definition: pmode.h:59
Mesh::getConnectivity
const CRConnectivity & getConnectivity(const StorageSite &from, const StorageSite &to) const
Definition: Mesh.cpp:416
COMETModel::FinishCoarseMesh
int FinishCoarseMesh(const MeshList &inMeshes, GeomFields &inGeomFields, MeshList &outMeshes, IntArray &CoarseCounts, map< const StorageSite *, IntArray * > &PreFacePairMap, IntArray &coarseGhost)
Definition: phononbase/COMETModel.h:1663
Kspace
Definition: Kspace.h:28
Kspace::findSpecs
void findSpecs(const Tvec n, const int m, const int k, Refl_pair &refls)
Definition: Kspace.h:742
COMETModel::TCOMETPtr
shared_ptr< TCOMET > TCOMETPtr
Definition: phononbase/COMETModel.h:44
Kspace::getFASArray
TArray & getFASArray()
Definition: Kspace.h:1229
COMETModel::BCellPtr
shared_ptr< BCcellArray > BCellPtr
Definition: phononbase/COMETModel.h:62
COMETDiscretizer::COMETSolveFull
void COMETSolveFull(const int dir, const int level)
Definition: COMETDiscretizer.h:403
COMETModel::COMETModel
COMETModel(const MeshList &meshes, const int level, GeomFields &geomFields, MacroFields &macroFields, Quadrature< T > &quad, const int ibm=0, GeomFields *finestGeomFields=NULL, const MeshList *finestMeshes=NULL, MacroFields *finestMacroFields=NULL)
Definition: esbgkbase/COMETModel.h:115
COMETModel::COMETBCMap
std::map< int, COMETBC< T > * > COMETBCMap
Definition: esbgkbase/COMETModel.h:90
Kspace::updateTau
void updateTau(const int c, const T Tl)
Definition: Kspace.h:1491
FaceGroup
Definition: Mesh.h:28
COMETModel::coarsenInterfaceCells
void coarsenInterfaceCells(COMETIC< T > &ic, IntArray &coarseCounts, GeomFields &inGeomFields, const MeshList &inMeshes)
Definition: phononbase/COMETModel.h:2509
COMETModel::Refl_Map
Tmode::Refl_Map Refl_Map
Definition: phononbase/COMETModel.h:50
COMETModel::TArrptr
shared_ptr< TArray > TArrptr
Definition: phononbase/COMETModel.h:39
COMETModel::makeCellColors
void makeCellColors(const int level)
Definition: phononbase/COMETModel.h:4018
StorageSite::setScatterProcID
void setScatterProcID(int proc_id)
Definition: StorageSite.h:51
COMETModel::getBCs
COMETBCMap & getBCs()
Definition: phononbase/COMETModel.h:127
Mesh::COMETfindCommonFaces
bool COMETfindCommonFaces(StorageSite &faces, StorageSite &otherFaces, const GeomFields &geomFields)
Definition: Mesh.cpp:1053
COMETModel::MshLstPtr
shared_ptr< MeshList > MshLstPtr
Definition: phononbase/COMETModel.h:51
COMETModel::T_Scalar
NumTypeTraits< T >::T_Scalar T_Scalar
Definition: esbgkbase/COMETModel.h:54
COMETModel::getGeomFields
GeomFields & getGeomFields()
Definition: esbgkbase/COMETModel.h:5641
pmode::calce0
T calce0(T Tl)
Definition: pmode.h:88
Kspace::setTauArray
void setTauArray(TArrPtr TauPtr)
Definition: Kspace.h:1223
COMETInterface::makeCoarseCoeffs
void makeCoarseCoeffs(const IClist &fineList, IClist &coarseList, MeshList &coarseMeshes)
Definition: COMETInterface.h:467
Field
Definition: Field.h:14
Kspace::getkvol
Tkvol & getkvol(int n) const
Definition: Kspace.h:390
COMETModel::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: phononbase/COMETModel.h:70
kvol::getmode
Tmode & getmode(int n) const
Definition: kvol.h:44
Kspace::seteArray
void seteArray(TArrPtr ePtr)
Definition: Kspace.h:1217
Mesh::createInteriorFaceGroup
const StorageSite & createInteriorFaceGroup(const int size)
Definition: Mesh.cpp:259
Kspace::getlength
int getlength() const
Definition: Kspace.h:391
COMETModel::SiteMap
map< const StorageSite *, StorageSite * > SiteMap
Definition: phononbase/COMETModel.h:72
COMETModel::Mode_ptr
Tmode::Mode_ptr Mode_ptr
Definition: phononbase/COMETModel.h:46
COMETModel::getWallAreaVector
VectorT3 getWallAreaVector(const Mesh &mesh, const int faceGroupId)
Definition: phononbase/COMETModel.h:3540
StorageSite::setTag
void setTag(int tag)
Definition: StorageSite.h:53
PhononMacro::deltaT
Field deltaT
Definition: PhononMacro.h:22
PhononMacro::temperature
Field temperature
Definition: PhononMacro.h:21
pmode::Mode_ptr
shared_ptr< pmode< T > > Mode_ptr
Definition: pmode.h:26
COMETModel::MakeNewKspaces
void MakeNewKspaces(TkspList &inList, TkspList &outList)
Definition: phononbase/COMETModel.h:1160
Mesh::getFaceGroup
const FaceGroup & getFaceGroup(const int fgId) const
Definition: Mesh.cpp:1570
Kspace::geteArray
TArray & geteArray()
Definition: Kspace.h:1224
Mesh
Definition: Mesh.h:49
COMETModel::Reflptr
Tmode::Reflptr Reflptr
Definition: phononbase/COMETModel.h:48
COMETModel::correctSingleNeighbor
int correctSingleNeighbor(const int m, const Mesh &mesh, GeomFields &inGeomFields, int coarseCount, map< const StorageSite *, IntArray * > PreFacePairMap)
Definition: phononbase/COMETModel.h:2357
COMETModel::BCcellList
vector< BCellPtr > BCcellList
Definition: phononbase/COMETModel.h:63
kvol::getmodenum
int getmodenum()
Definition: kvol.h:43
Kspace::getDK3
T getDK3() const
Definition: Kspace.h:408
Kspace::calcLatTemp
T calcLatTemp(const int c)
Definition: Kspace.h:461
COMETModel::MakeInteriorCoarseMesh
void MakeInteriorCoarseMesh(const MeshList &inMeshes, GeomFields &inGeomFields, MeshList &outMeshes, IntArray &CoarseCounts, map< const StorageSite *, IntArray * > &PreFacePairMap, IntArray &CoarseGhost, SiteMap &siteMap)
Definition: phononbase/COMETModel.h:1176
StorageSite::getCommonMap
const CommonMap & getCommonMap() const
Definition: StorageSite.h:60
COMETModel::HeatFluxIntegralFace
T HeatFluxIntegralFace(const Mesh &mesh, const int f)
Definition: phononbase/COMETModel.h:3351
COMETModel::FieldVector
vector< shared_ptr< Field > > FieldVector
Definition: phononbase/COMETModel.h:69
COMETModel::_options
COMETModelOptions< T > _options
Definition: esbgkbase/COMETModel.h:5677
COMETModel::getMKMap
MeshKspaceMap & getMKMap()
Definition: phononbase/COMETModel.h:128
COMETModel::TKCptr
TKConnectivity * TKCptr
Definition: phononbase/COMETModel.h:66
pmode::getreflmap
Refl_Map & getreflmap()
Definition: pmode.h:69
StorageSite::setCount
void setCount(const int selfCount, const int nGhost=0)
Definition: StorageSite.h:42
COMETModel::doSweeps
void doSweeps(const int sweeps)
Definition: phononbase/COMETModel.h:2800
UnorderedSet::insert
void insert(int x)
Definition: MultiFieldMatrix.h:25
COMETModel::setBCMap
void setBCMap(COMETBCMap *bcMap)
Definition: phononbase/COMETModel.h:4238
Mesh::setConnectivity
void setConnectivity(const StorageSite &rowSite, const StorageSite &colSite, shared_ptr< CRConnectivity > conn)
Definition: Mesh.cpp:352
COMETModel::setLocalScatterMaps
void setLocalScatterMaps(const Mesh &mesh0, StorageSite &faces0, const Mesh &mesh1, StorageSite &faces1)
Definition: phononbase/COMETModel.h:988
COMETModel::TKClist
vector< TKCptr > TKClist
Definition: phononbase/COMETModel.h:67
COMETIC::FineToCoarse1
Array< int > * FineToCoarse1
Definition: phononbase/COMETBC.h:155
Kspace::gete
T & gete(const int cell, const int count)
Definition: Kspace.h:1211
COMETModel::TCModOpts
COMETModelOptions< T > TCModOpts
Definition: phononbase/COMETModel.h:45
COMETModel::getMeshICmap
MeshICmap & getMeshICmap()
Definition: phononbase/COMETModel.h:4259
Mesh::createBoundaryFaceGroup
const StorageSite & createBoundaryFaceGroup(const int size, const int offset, const int id, const string &boundaryType)
Definition: Mesh.cpp:278
COMETBoundary::applyTemperatureWallFine
void applyTemperatureWallFine(FloatValEvaluator< T > &bTemp) const
Definition: COMETBoundary.h:61
pmode::Refl_Map
map< int, Refl_pair > Refl_Map
Definition: pmode.h:30
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
COMETModel::SiteMap
map< const StorageSite *, StorageSite * > SiteMap
Definition: esbgkbase/COMETModel.h:102
COMETModel::BfacePtr
shared_ptr< BCfaceArray > BfacePtr
Definition: phononbase/COMETModel.h:59
FloatVarDict::getVal
FloatVal< T > getVal(const string varName) const
Definition: FloatVarDict.h:85
COMETModel::setFinerLevel
void setFinerLevel(TCOMET *fl)
Definition: phononbase/COMETModel.h:4240
COMETModel::getModelPointer
TCOMET * getModelPointer(const int level)
Definition: phononbase/COMETModel.h:4173
COMETModel::TCOMET
COMETModel< T > TCOMET
Definition: phononbase/COMETModel.h:43
COMETInterface::updateResid
void updateResid(const COMETIC< T > &ic, const bool plusFAS)
Definition: COMETInterface.h:885
COMETModel::getValueArray
ArrayBase * getValueArray(const Mesh &mesh, const int cell)
Definition: phononbase/COMETModel.h:3580
COMETModel::GeoFldsPtr
shared_ptr< GeomFields > GeoFldsPtr
Definition: phononbase/COMETModel.h:53
COMETModel::syncGhostCoarsening
void syncGhostCoarsening(const MeshList &inMeshes, GeomFields &inGeomFields, MeshList &outMeshes, ScatGathMaps &coarseScatMaps, ScatGathMaps &coarseGathMaps, IntArray &coarseSizes, IntArray &CoarseGhost)
Definition: phononbase/COMETModel.h:1433
GeomFields::fineToCoarse
Field fineToCoarse
Definition: GeomFields.h:41
sqrt
Tangent sqrt(const Tangent &a)
Definition: Tangent.h:317
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
COMETModel::smooth
void smooth(int dir)
Definition: phononbase/COMETModel.h:2810
COMETModel::syncGhostCoarsening
void syncGhostCoarsening(const MeshList &inMeshes, GeomFields &inGeomFields, MeshList &outMeshes)
Definition: esbgkbase/COMETModel.h:3214
Mesh::getCellFaces
const CRConnectivity & getCellFaces() const
Definition: Mesh.cpp:454
COMETIC::getSelfFaceID
int getSelfFaceID(const int mesh)
Definition: phononbase/COMETBC.h:118
FaceGroup::id
const int id
Definition: Mesh.h:41
COMETInterface::makeDMMcoeffs
void makeDMMcoeffs(COMETIC< T > &ic)
Definition: COMETInterface.h:44
COMETModel::MeshKspaceMap
map< int, int > MeshKspaceMap
Definition: phononbase/COMETModel.h:42
COMETBC::bcType
string bcType
Definition: esbgkbase/COMETBC.h:32
COMETModel::applyTemperatureBoundaries
void applyTemperatureBoundaries()
Definition: phononbase/COMETModel.h:1011
Kspace::setResArray
void setResArray(TArrPtr ResPtr)
Definition: Kspace.h:1221
COMETModel::Reflection
Tmode::Reflection Reflection
Definition: phononbase/COMETModel.h:47
COMETModel::getOptions
COMETModelOptions< T > & getOptions()
Definition: esbgkbase/COMETModel.h:2261
pmode::getIndex
int getIndex()
Definition: pmode.h:73
COMETBoundary::applyTemperatureWallCoarse
void applyTemperatureWallCoarse(FloatValEvaluator< T > &bTemp) const
Definition: COMETBoundary.h:72
COMETModel::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: phononbase/COMETModel.h:28
PhononMacro::TlResidual
Field TlResidual
Definition: PhononMacro.h:24
Kspace::getGlobalIndex
int getGlobalIndex(const int cell, const int count)
Definition: Kspace.h:1216
Mesh::createInterfaceGroup
const StorageSite & createInterfaceGroup(const int size, const int offset, const int id)
Definition: Mesh.cpp:268
StorageSite::getTag
int getTag() const
Definition: StorageSite.h:84
COMETModel::BCfaceList
vector< BfacePtr > BCfaceList
Definition: phononbase/COMETModel.h:60
StorageSite::getScatterProcID
int getScatterProcID() const
Definition: StorageSite.h:82
StorageSite::getGatherProcID
int getGatherProcID() const
Definition: StorageSite.h:83
COMETInterface::updateOtherGhost
void updateOtherGhost(const COMETIC< T > &ic, const int Mid0, const bool plusFAS)
Definition: COMETInterface.h:820
COMETModel::FgTKClistMap
map< int, TKClist > FgTKClistMap
Definition: phononbase/COMETModel.h:68
COMETModel::setStraightLine
void setStraightLine(const T T1, const T T2)
Definition: phononbase/COMETModel.h:3641
COMETModel::SSPtr
shared_ptr< StorageSite > SSPtr
Definition: phononbase/COMETModel.h:56
PhononMacro::TlFASCorrection
Field TlFASCorrection
Definition: PhononMacro.h:26
Kspace::calcTauTot
T calcTauTot()
Definition: Kspace.h:409
COMETModel::IntArrPtr
shared_ptr< IntArray > IntArrPtr
Definition: phononbase/COMETModel.h:37
COMETModel::setFinestLevel
void setFinestLevel(TCOMET *finest)
Definition: phononbase/COMETModel.h:1153
Kspace::getResArray
TArray & getResArray()
Definition: Kspace.h:1228
StorageSite::CommonMap
map< const StorageSite *, shared_ptr< Array< int > > > CommonMap
Definition: StorageSite.h:25
DensityOfStates::getMIndices
IntArray & getMIndices(const int fBin)
Definition: DensityOfStates.h:244
Kspace::calcModeTemp
T calcModeTemp(T guess, const T e_sum, const T m)
Definition: Kspace.h:501
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
COMETModel::_macro
PhononMacro & _macro
Definition: phononbase/COMETModel.h:4266
COMETIC::InterfaceModel
string InterfaceModel
Definition: phononbase/COMETBC.h:147
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Mesh::getCellCells
const CRConnectivity & getCellCells() const
Definition: Mesh.cpp:480
Model
Definition: Model.h:13
UnorderedSet::size
int size() const
Definition: MultiFieldMatrix.h:34
COMETModel::getMacro
MacroFields & getMacro()
Definition: esbgkbase/COMETModel.h:5643
PhononMacro::CellColors
FieldVector CellColors
Definition: PhononMacro.h:34
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
DensityOfStates::getFreqMidsT
TArray & getFreqMidsT()
Definition: DensityOfStates.h:220
COMETModel::_MeshKspaceMap
MeshKspaceMap _MeshKspaceMap
Definition: phononbase/COMETModel.h:4275
GeomFields::volume
Field volume
Definition: GeomFields.h:26
PhononMacro::TlInjected
Field TlInjected
Definition: PhononMacro.h:25
COMETModel::makeCoarseScatGath
void makeCoarseScatGath(const MeshList &inMeshes, SiteMap &siteMap, ScatGathMaps &coarseScatMaps, ScatGathMaps &coarseGathMaps)
Definition: phononbase/COMETModel.h:1614
PhononMacro::BandRelEnergy
FieldVectorMap BandRelEnergy
Definition: PhononMacro.h:32
COMETInterface::makeNoInterfaceCoeffs
void makeNoInterfaceCoeffs(COMETIC< T > &ic)
Definition: COMETInterface.h:399
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
COMETModel::IClist
vector< COMETIC< T > * > IClist
Definition: phononbase/COMETModel.h:41
COMETModel::CRPtr
shared_ptr< CRConnectivity > CRPtr
Definition: phononbase/COMETModel.h:57
COMETModel::setCoarserLevel
void setCoarserLevel(TCOMET *cl)
Definition: esbgkbase/COMETModel.h:5637
StorageSite::getScatterMap
const ScatterMap & getScatterMap() const
Definition: StorageSite.h:58
COMETModel::getMeshList
const MeshList & getMeshList()
Definition: esbgkbase/COMETModel.h:5640
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
Kspace::getInjArray
TArray & getInjArray()
Definition: Kspace.h:1227
Mesh::getFaceCells
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
COMETModel::BArray
Array< bool > BArray
Definition: phononbase/COMETModel.h:35
StorageSite::ScatterMap
map< const StorageSite *, shared_ptr< Array< int > > > ScatterMap
Definition: StorageSite.h:23
UnorderedSet::getData
const vector< int > & getData() const
Definition: MultiFieldMatrix.h:33
COMETDiscretizer::COMETSolveFine
void COMETSolveFine(const int dir, const int level)
Definition: COMETDiscretizer.h:91
Kspace::getDOSptr
DensityOfStates< T > * getDOSptr()
Definition: Kspace.h:1207
COMETModel::VT3Ptr
shared_ptr< VectorT3Array > VT3Ptr
Definition: phononbase/COMETModel.h:29
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
Array
Definition: Array.h:14
StorageSite::setGatherProcID
void setGatherProcID(int proc_id)
Definition: StorageSite.h:52
pmode::getfield
Field & getfield()
Definition: pmode.h:74
StorageSite::getOffset
int getOffset() const
Definition: StorageSite.h:87
CRConnectivity::multiply
shared_ptr< CRConnectivity > multiply(const CRConnectivity &b, const bool implicitDiagonal) const
Definition: CRConnectivity.cpp:168
Kspace::CopyKspace
void CopyKspace(Tkspace &copyFromKspace)
Definition: Kspace.h:781
COMETModel::COMETModel
COMETModel(const MeshList &meshes, const int level, GeomFields &geomFields, TkspList &kspaces, PhononMacro &macro)
Definition: phononbase/COMETModel.h:74
Kspace::gete0Array
TArray & gete0Array()
Definition: Kspace.h:1225
StorageSite::getGatherMap
const GatherMap & getGatherMap() const
Definition: StorageSite.h:59
COMETDiscretizer::COMETSolveCoarse
void COMETSolveCoarse(const int dir, const int level)
Definition: COMETDiscretizer.h:244
PhononMacro::plottingCellColors
FieldVector plottingCellColors
Definition: PhononMacro.h:35
PhononMacro::BandTemperatures
FieldVectorMap BandTemperatures
Definition: PhononMacro.h:31
Kspace::gettotmodes
int gettotmodes()
Definition: Kspace.h:393
CRConnectivity::getTranspose
shared_ptr< CRConnectivity > getTranspose() const
Definition: CRConnectivity.cpp:74
COMETModel::TkspList
vector< Tkspace * > TkspList
Definition: phononbase/COMETModel.h:32
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
COMETModel::VectorT3
Vector< T_Scalar, 3 > VectorT3
Definition: phononbase/COMETModel.h:27
pmode::Reflection
pair< T, int > Reflection
Definition: pmode.h:27
COMETIC::FineToCoarse0
Array< int > * FineToCoarse0
Definition: phononbase/COMETBC.h:154
COMETModel::initializeTemperatureBoundaries
void initializeTemperatureBoundaries()
Definition: phononbase/COMETModel.h:1045
PhononMacro::BranchFlux
FieldVectorMap BranchFlux
Definition: PhononMacro.h:30
COMETModel::IntArray
Array< int > IntArray
Definition: phononbase/COMETModel.h:36
DensityOfStates::getKIndices
IntArray & getKIndices(const int fBin)
Definition: DensityOfStates.h:243
GeomFields::area
Field area
Definition: GeomFields.h:23
PhononMacro::BranchTemperatures
FieldVectorMap BranchTemperatures
Definition: PhononMacro.h:29
StorageSite::GatherMap
map< const StorageSite *, shared_ptr< Array< int > > > GatherMap
Definition: StorageSite.h:24
Mesh::getInterfaceGroups
const FaceGroupList & getInterfaceGroups() const
Definition: Mesh.h:190
COMETModel::T_Scalar
NumTypeTraits< T >::T_Scalar T_Scalar
Definition: phononbase/COMETModel.h:26
COMETModel::Tkspace
Kspace< T > Tkspace
Definition: phononbase/COMETModel.h:30
COMETModel::coarsenInterior
int coarsenInterior(const int m, const Mesh &mesh, GeomFields &inGeomFields, int &coarseCount)
Definition: phononbase/COMETModel.h:2249
Kspace::makeFAS
void makeFAS()
Definition: Kspace.h:1289
Vector::zero
void zero()
Definition: Vector.h:156
COMETModel::COMETBCMap
map< int, COMETBC< T > * > COMETBCMap
Definition: phononbase/COMETModel.h:40
PhononMacro::BandFlux
FieldVectorMap BandFlux
Definition: PhononMacro.h:33
COMETModel::getKspaces
TkspList & getKspaces()
Definition: phononbase/COMETModel.h:4245
COMETModel::smooth
void smooth(const int num, const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:4927
PhononMacro::getColorField
Field & getColorField(int level)
Definition: PhononMacro.cpp:101
Kspace::calcBandTemp
T calcBandTemp(const T guess, const T eSum, const IntArray &kpts, const IntArray &mpts)
Definition: Kspace.h:1103
Kspace::setSourceArray
void setSourceArray(TArrPtr SPtr)
Definition: Kspace.h:1219
pmode
Definition: pmode.h:18
COMETModel::Refl_pair
Tmode::Refl_pair Refl_pair
Definition: phononbase/COMETModel.h:49
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
COMETModel::ScatGathMaps
map< SSPair, shared_ptr< Array< int > > > ScatGathMaps
Definition: phononbase/COMETModel.h:71
GeomFields::areaMag
Field areaMag
Definition: GeomFields.h:25
kvol
Definition: kvol.h:14
pmode::gettau
T gettau()
Definition: pmode.h:61
PhononMacro::heatFlux
Field heatFlux
Definition: PhononMacro.h:27
COMETModel::HeatFluxIntegral
T HeatFluxIntegral(const Mesh &mesh, const int faceGroupId)
Definition: phononbase/COMETModel.h:3291
Field::syncLocal
void syncLocal()
Definition: Field.cpp:334
COMETModel::KspPtr
shared_ptr< Tkspace > KspPtr
Definition: phononbase/COMETModel.h:54
Mesh::getID
int getID() const
Definition: Mesh.h:106
Kspace::sete0Array
void sete0Array(TArrPtr e0Ptr)
Definition: Kspace.h:1218
COMETModel::getWallArea
T getWallArea(const Mesh &mesh, const int faceGroupId)
Definition: phononbase/COMETModel.h:3513
COMETModel::MeshPtr
shared_ptr< Mesh > MeshPtr
Definition: phononbase/COMETModel.h:52
COMETModel::PMacroPtr
shared_ptr< PhononMacro > PMacroPtr
Definition: phononbase/COMETModel.h:55
COMETModel::advance
void advance(const int iters)
Definition: phononbase/COMETModel.h:3270
COMETModel::makePlotColors
void makePlotColors(const int level)
Definition: phononbase/COMETModel.h:4074
COMETModel::MakeCoarseModel
void MakeCoarseModel(TCOMET *finerModel)
Definition: esbgkbase/COMETModel.h:575
COMETModel::BCfaceArray
Array< int > BCfaceArray
Definition: phononbase/COMETModel.h:58
COMETModel::doSweeps
void doSweeps(const int sweeps, const int num)
Definition: esbgkbase/COMETModel.h:4915
MeshList
vector< Mesh * > MeshList
Definition: Mesh.h:439
Kspace::setInjArray
void setInjArray(TArrPtr InjPtr)
Definition: Kspace.h:1220
COMETModel::sameFaceGroup
bool sameFaceGroup(const Mesh &mesh, const int f0, const int f1)
Definition: phononbase/COMETModel.h:4181
kvol::getdk3
T getdk3()
Definition: kvol.h:42
COMETModel::getFinerModel
TCOMET * getFinerModel()
Definition: phononbase/COMETModel.h:4241
COMETIC::getKConnectivity
const TKClist & getKConnectivity(const int fgid) const
Definition: phononbase/COMETBC.h:98
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
Array::getLength
int getLength() const
Definition: Array.h:87
COMETDiscretizer::findResidCoarse
void findResidCoarse(const bool plusFAS)
Definition: COMETDiscretizer.h:1314