Chombo + EB + MF  3.2
CutCellMomentsImplem.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _CUTCELLMOMENTSIMPLEM_H_
12 #define _CUTCELLMOMENTSIMPLEM_H_
13 
14 #if defined(CH_Darwin) && defined(__GNUC__) && ( __GNUC__ == 3 )
15 // deal with the broken isnan()/isinf() in GCC on MacOS
16 #include <unistd.h>
17 #define _GLIBCPP_USE_C99 1
18 #endif
19 
20 #include <cmath>
21 #include <cstdio>
22 #include <cstdlib>
23 #include <fstream>
24 #include <iostream>
25 #include <iomanip>
26 #include <string>
27 
28 #include "MayDay.H"
29 
30 #include "LSProblem.H"
31 
32 #include "NamespaceHeader.H"
33 
34 // Null constructor
36 {
37 }
38 
39 // Copy constructor
40 template <int dim> CutCellMoments<dim>::CutCellMoments(const CutCellMoments<dim>& a_cutCellMoments)
41  :m_moments (a_cutCellMoments.m_moments),
42  m_EBmoments (a_cutCellMoments.m_EBmoments),
43  m_bdCutCellMoments(a_cutCellMoments.m_bdCutCellMoments),
44  m_IFData (a_cutCellMoments.m_IFData),
45  m_bdCCOn (a_cutCellMoments.m_bdCCOn),
46  m_residual (a_cutCellMoments.m_residual),
47  m_numActiveBounds (a_cutCellMoments.m_numActiveBounds),
48  m_badNormal (a_cutCellMoments.m_badNormal)
49 {
50 }
51 
52 // Use this for initializing
53 template <int dim> CutCellMoments<dim>::CutCellMoments(const IFData<dim>& a_info)
54  :m_IFData(a_info),
56 {
57  m_bdCCOn = false;
58  m_badNormal = false;
59 
60  for (int hilo = 0; hilo < 2; ++hilo)
61  {
62  for (int idir = 0; idir < dim; ++idir)
63  {
64  //identifier for which boundary cutCellMoment
65  Iv2 bdId;
66  bdId[BDID_DIR] = idir;
67  bdId[BDID_HILO] = hilo;
68 
69  IFData<dim-1> reducedInfo(a_info,a_info.m_maxOrder+1,idir,hilo);
70  CutCellMoments<dim-1>bdCutCellMoments(reducedInfo);
71 
72  m_bdCutCellMoments[bdId] = bdCutCellMoments;
73 
74 
75  // Notice whether reducedInfo coincides with the interface
76  if (reducedInfo.m_allVerticesOn)
77  {
78  m_bdCCOn = true;
79  }
80 
81  // Notice whether reducedInfo.m_normal = zero vector
82  if (reducedInfo.m_badNormal)
83  {
84  m_badNormal = true;
85  }
86 
87  }
88  }
89 }
90 
91 // Destructor
93 {
94 }
95 
96 template <int dim> const CutCellMoments<dim - 1> CutCellMoments<dim>::getBdCutCellMoments(const Iv2 & a_bdId) const
97 {
98  typename BdCutCellMoments::const_iterator it = m_bdCutCellMoments.find(a_bdId);
99 
100  if (it == m_bdCutCellMoments.end())
101  {
102  MayDay::Abort("Can't find this bdId in m_bdCutCellMoments");
103  }
104 
105  return it->second;
106 }
107 
108 // This function is called to fill the boundary moments from the higher
109 // dimension: it fills the moments in a_coarseCutCell using either the
110 // moments in the current cutcell map or the full cell/covered cell answers
111 // depending on the value of a_IFData.m_allVerticesIn/Out
112 template <int dim> void CutCellMoments<dim>::addBdMoments(CutCellMoments<dim> & a_coarseBdCutCell,
113  const IFData<dim+1> & a_IFData,
114  const int & a_degreePmax,
115  const bool & a_useConstraints,
116  const IndexTM<Real,dim> & a_refinedCenterDelta,
117  const IndexTM<int,dim> & a_hilo)
118 {
119  typedef map<IvDim,int > PthMomentLoc;
120 
121  if (a_IFData.m_allVerticesIn)
122  {
123  // The boundary map doesn't exist, we need to compute the full cell
124  // quadrature for all monomial in order to use the changecoordinates
125  // function. Generate a full cell quadrature map.
126  PthMoment fullCellMap;
127 
128  for (int iDegree = a_degreePmax; iDegree > 0; iDegree--)
129  {
130  // Creating a lsproblem generates the monomial map that is
131  // needed to fill the moment maps
132  LSProblem<dim> lsp(iDegree,a_useConstraints);
133  const PthMomentLoc& momMap = lsp.getMonomialLocMapDegreePLess1();
134 
135  for (typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
136  {
137  // pout()<<"mono = "<<it->first<<endl;
138  fullCellMap[it->first] = fullCellQuadrature(it->first,m_IFData.m_cellCenterCoord);
139  }
140  }
141 
142  for (int iDegree = a_degreePmax; iDegree > 0; iDegree--)
143  {
144  LSProblem<dim> lsp(iDegree,a_useConstraints);
145 
146  // Change the coordinates of the refined moments using the full
147  // cell quadrature map. Add up the obtained moment to the
148  // coarse map
149  const PthMomentLoc& momMap = lsp.getMonomialLocMapDegreePLess1();
150 
151  for (typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
152  {
153  Real bdMoment = getBdMoment(it->first,a_IFData,a_refinedCenterDelta,fullCellMap);
154  a_coarseBdCutCell.m_moments[it->first] += bdMoment;
155  }
156  }
157  }
158  else
159  {
160  for (int iDegree = a_degreePmax; iDegree > 0; iDegree--)
161  {
162  LSProblem<dim> lsp(iDegree,a_useConstraints);
163 
164  // Get the moments from the refined cell map + change coordinates
165  // function if the refined cell is irregular or they are 0 as all
166  // vertices are out
167  const PthMomentLoc& momMap = lsp.getMonomialLocMapDegreePLess1();
168 
169  for (typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
170  {
171  Real bdMoment = getBdMoment(it->first,a_IFData,a_refinedCenterDelta);
172  a_coarseBdCutCell.m_moments[it->first] += bdMoment;
173  }
174  }
175  }
176 
177  for (int iDegree = a_degreePmax; iDegree >= 0; iDegree--)
178  {
179  LSProblem<dim> lsp(iDegree,a_useConstraints);
180 
181  // EB moments are either found in the EB moment map or equals to 0
182  // (allIn or allOut)
183  const PthMomentLoc& momMapP = lsp.getMonomialLocMapDegreeP();
184 
185  for (typename PthMomentLoc::const_iterator it = momMapP.begin(); it != momMapP.end(); ++it)
186  {
187  Real bdEBMoment = getBdEBMoment(it->first,a_IFData,a_refinedCenterDelta);
188  a_coarseBdCutCell.m_EBmoments[it->first] += bdEBMoment;
189  }
190  }
191 
192  IndexTM<int,2> bdId;
193 
194  for (int idir = 0; idir < dim; idir++)
195  {
196  bdId[0] = idir;
197 
198  // The refined cell center in dim-1 is needed to change the coordinate
199  // system of the refined moments to the coordinate system where they are
200  // needed.
201  IndexTM<Real,dim-1> localRefinedCellCenter;
202  IndexTM<int,dim-1> localHilo;
203 
204  for (int jdir = 0; jdir < dim; jdir++)
205  {
206  if (jdir < idir)
207  {
208  localRefinedCellCenter[jdir] = a_refinedCenterDelta[jdir];
209  localHilo[jdir] = a_hilo[jdir];
210  }
211  else if (jdir > idir)
212  {
213  localRefinedCellCenter[jdir-1] = a_refinedCenterDelta[jdir];
214  localHilo[jdir-1] = a_hilo[jdir];
215  }
216  }
217 
218  if (a_hilo[idir] != LARGEINTVAL)
219  {
220  bdId[1] = a_hilo[idir];
221  // Drop one dimension to add up the boundary moments. Dropping one
222  // dimension allows an easier implementation as 1D is a particular
223  // case
224  (m_bdCutCellMoments[bdId]).addBdMoments(a_coarseBdCutCell.m_bdCutCellMoments[bdId],m_IFData,a_degreePmax+1,a_useConstraints,localRefinedCellCenter,localHilo);
225  }
226  else
227  {
228  for (int side = 0; side < 2; side++)
229  {
230  bdId[1]=side;
231  (m_bdCutCellMoments[bdId]).addBdMoments(a_coarseBdCutCell.m_bdCutCellMoments[bdId],m_IFData,a_degreePmax,a_useConstraints,localRefinedCellCenter,localHilo);
232  }
233  }
234  }
235 }
236 
238  const CoordinateSystem<dim> & a_coord)
239 {
240  Real moment = 1.0;
241 
242  for (int idir = 0; idir < dim; ++idir)
243  {
244  Real loPt = a_coord.convertDir(-0.5*m_IFData.m_cellCenterCoord.m_dx[idir],
245  m_IFData.m_cellCenterCoord,
246  idir);
247 
248  Real hiPt = a_coord.convertDir(0.5*m_IFData.m_cellCenterCoord.m_dx[idir],
249  m_IFData.m_cellCenterCoord,
250  idir);
251 
252  Real idirInt = pow(hiPt,a_mono[idir] + 1) - pow(loPt,a_mono[idir] + 1);
253 
254  idirInt /= (a_mono[idir] + 1);
255  moment *= idirInt;
256  }
257 
258  return moment;
259 }
260 
261 template <int dim> Real CutCellMoments<dim>::changeMomentCoordinates(PthMoment & a_refinedMomentMap,
262  const IndexTM<int,dim> & a_monomial,
263  const IndexTM<Real,dim> & a_refinedCenterDelta)
264 {
265  // This function outputs the moment in the coordinate system where the
266  // origin is the current cell center, given the moment map of the refined
267  // moments computed in the coordinate system where the origin is the
268  // refined cell center it uses a Taylor development around 0 of the function
269  // (x+xc)^p1*(y+yc)^p2*...
270  Real moment = 0.0;
271 
272  // Degree is the degree of the monomial
273  int degree = 0;
274 
275  for (int i = 0; i < dim; i++)
276  {
277  degree += a_monomial[i];
278  }
279 
280  // Loop over the possible values of r in the development of the Taylor series
281  // The Taylor series has a finite expansion as the function's derivatives
282  // are zero for r>degree
283  for (int r = 0; r <= degree; r++)
284  {
285  if (r >= 0)
286  {
287  // Generate all the possible monomials of degree r and add up the moment
288  IndexTM<int,dim> derivative;
289 
290  for (int idir = 0; idir < dim; ++idir)
291  {
292  derivative[idir] = 0;
293  }
294 
295  while (true)
296  {
297  for (int j = 1; j < dim-1; ++j)
298  {
299  int t = r;
300  for (int k = j+1; k < dim; ++k)
301  {
302  t -= derivative[k];
303  }
304 
305  if (derivative[j] > t)
306  {
307  derivative[j+1] += 1;
308  derivative[j ] = 0;
309  }
310  else
311  {
312  break;
313  }
314  }
315 
316  if (derivative[dim-1] > r)
317  {
318  break;
319  }
320 
321  derivative[0] = r;
322 
323  for (int j = 1; j < dim; ++j)
324  {
325  derivative[0] -= derivative[j];
326  }
327 
328  // Add up the relevant term of the refined moment map to the map
329  Real coeff = 1.0;
330  Real factorial = 1.0;
331 
332  for (int idir = 0; idir < dim; idir++)
333  {
334  for (int j = 0; j < derivative[idir]; j++)
335  {
336  coeff *= a_monomial[idir] - j;
337  factorial *= j+1;
338  }
339 
340  // This is to prevent the computation of 0*0^-1 which result
341  // in a nan
342  if (a_monomial[idir]-derivative[idir] < 0)
343  {
344  if (coeff != 0)
345  {
346  MayDay::Abort("Coeff should be zero when the derivative has higher coefficient than the monomial");
347  }
348  }
349  else
350  {
351  coeff *= pow(a_refinedCenterDelta[idir],a_monomial[idir]-derivative[idir]);
352  }
353  }
354 
355  //pout()<<a_refinedMomentMap[derivative]<<endl;
356  moment += coeff * a_refinedMomentMap[derivative] / factorial;
357 
358  // Increments to get to the next derivative
359  derivative[1] += 1;
360  }
361  }
362  }
363 
364  return moment;
365 }
366 
368 {
369  // Move moments from parent coord to cell center coord
370  IndexTM<Real,dim> delta = m_IFData.m_cellCenterCoord.m_origin;
371  delta -= m_IFData.m_parentCoord .m_origin;
372 
373  PthMoment copyMoments = m_moments;
374  for (typename PthMoment::const_iterator it = copyMoments.begin();
375  it != copyMoments.end(); ++it)
376  {
377  IvDim mono = it->first;
378  m_moments[mono] = changeMomentCoordinates(copyMoments, mono, delta);
379  }
380 
381  PthMoment copyEBMoments = m_EBmoments;
382  for (typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
383  {
384  IvDim mono = it->first;
385  m_EBmoments[mono] = changeMomentCoordinates(copyEBMoments, mono, delta);
386  }
387 }
388 
390 {
391  // Move moments from cell center coord to parent coord
392  IndexTM<Real,dim> delta = m_IFData.m_parentCoord .m_origin;
393  delta -= m_IFData.m_cellCenterCoord.m_origin;
394 
395  PthMoment copyMoments = m_moments;
396  for (typename PthMoment::const_iterator it = copyMoments.begin();
397  it != copyMoments.end(); ++it)
398  {
399  IvDim mono = it->first;
400  m_moments[mono] = changeMomentCoordinates(copyMoments, mono, delta);
401  }
402 
403  PthMoment copyEBMoments = m_EBmoments;
404  for (typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
405  {
406  IvDim mono = it->first;
407  m_EBmoments[mono] = changeMomentCoordinates(copyEBMoments, mono, delta);
408  }
409 }
410 
411 template <int dim> void CutCellMoments<dim>::initialize(CutCellMoments<dim> & a_refinedCutCell)
412 {
413  // This initialize the maps of the current cut cell moment to zero by
414  // iterating through the refined maps
415  initializeMap(m_EBmoments,a_refinedCutCell.m_EBmoments);
416  initializeMap(m_moments,a_refinedCutCell.m_moments);
417 
418  // Initialize maps on the boundary
419  IndexTM<int,2> bdId;
420 
421  for (int idir = 0; idir < dim; idir++)
422  {
423  bdId[0] = idir;
424 
425  for (int hilo = 0; hilo < 2; hilo++)
426  {
427  bdId[1] = hilo;
428  CutCellMoments<dim-1> refinedBdCutCell = a_refinedCutCell.getBdCutCellMoments(bdId);
429 
430  m_bdCutCellMoments[bdId].initialize(refinedBdCutCell);
431  }
432  }
433 }
434 
436  PthMomentLesserDimension & a_map2)
437 {
438  for (typename PthMomentLesserDimension::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
439  {
440  a_map1[it->first] = 0.0;
441  }
442 }
443 
444 template <int dim> void CutCellMoments<dim>::initializeMap(PthMoment & a_map1,
445  PthMoment & a_map2)
446 {
447  for (typename PthMoment::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
448  {
449  a_map1[it->first] = 0.0;
450  }
451 }
452 
453 
454 template <int dim> Real CutCellMoments<dim>::getBdMoment(const IvDim & a_mono,
455  const IFData<dim+1> & a_IFData,
456  const IndexTM<Real,dim> & a_refinedCenterDelta,
457  PthMoment a_fullCellMap)
458 {
459  // This computes the moments on the fly in the case the refined cell (in
460  // dim+1) was all in or all out (in that case the boundary moments were
461  // never computed) or looks for them in the moment map
462  Real moment = LARGEREALVAL;
463 
464  if (a_IFData.m_allVerticesOut)
465  {
466  moment = 0.0;
467  }
468  else if (a_IFData.m_allVerticesIn)
469  {
470  moment = changeMomentCoordinates(a_fullCellMap,a_mono,a_refinedCenterDelta);
471  }
472  else
473  {
474  moment = changeMomentCoordinates(m_moments,a_mono,a_refinedCenterDelta);
475  }
476 
477  return moment;
478 }
479 
480 template <int dim> Real CutCellMoments<dim>::getBdEBMoment(const IvDim & a_mono,
481  const IFData<dim+1> & a_IFData,
482  const IndexTM<Real,dim> & a_refinedCenterDelta)
483 {
484  // EB moments are 0 if the refined cell (in dim+1) is all in or all out,
485  // otherwise they are in the current cell EB moment maps
486  Real EBmoment = LARGEREALVAL;
487 
488  if (a_IFData.m_allVerticesOut || a_IFData.m_allVerticesIn)
489  {
490  EBmoment = 0.0;
491  }
492  else
493  {
494  EBmoment = changeMomentCoordinates(m_EBmoments,a_mono,a_refinedCenterDelta);
495  }
496 
497  return EBmoment;
498 }
499 
500 // Methods for reading moments that check the moment exists
501 template <int dim> Real CutCellMoments<dim>::getMoment(const IvDim & a_mono,
502  const EBorVol & a_EBorVol) const
503 {
504  Real moment = LARGEREALVAL;
505 
506  if (a_EBorVol == VolMoment)
507  {
508  // Find the vol in the map
509  typename PthMoment::const_iterator it = m_moments.find(a_mono);
510 
511  if (it != m_moments.end())
512  {
513  moment = it->second;
514  }
515  else
516  {
517  MayDay::Abort("No volume moment in m_moments");
518  }
519  }
520  else if (a_EBorVol == EBMoment)
521  {
522  // Find the vol in the EBmap
523  typename PthMoment::const_iterator it = m_EBmoments.find(a_mono);
524 
525  if (it != m_EBmoments.end())
526  {
527  moment = it->second;
528  }
529  else
530  {
531  MayDay::Abort("No volume moment in m_bdMoments");
532  }
533  }
534  else
535  {
536  MayDay::Abort("Must ask for EBMoment or VolMoment");
537  }
538 
539  return moment;
540 }
541 
542 template <int dim> Real CutCellMoments<dim>::getVol(const EBorVol& a_EBorVol) const
543 {
544  Real volume;
545 
546  // Volume calculated by integrating x^0 = 1
547  IvDim zeroIv = IvDim::Zero;
548  volume = getMoment(zeroIv,a_EBorVol);
549 
550  return volume;
551 }
552 
553 template <int dim> IndexTM<Real,dim> CutCellMoments<dim>::getCentroid(const EBorVol& a_EBorVol) const
554 {
555  // This call checks the string to be either "VolMoment" or "EBMoment"
556  Real volume = getVol(a_EBorVol);
557 
558  RvDim centroid;
559 
560  // Mono= 0,0,...1,0,0
561  for (int idir = 0; idir < dim; ++idir)
562  {
563  IvDim mono = BASISV_TM<int,dim>(idir);
564  centroid[idir] = getMoment(mono,a_EBorVol);
565 
566  if (volume < 0.0)
567  {
568  MayDay::Warning("CutCellMoments::getCentroid: Volume fraction is negative");
569  }
570  else if (volume == 0.0)
571  {
572  //MayDay::Abort("CutCellMoments::getCentroid: Volume fraction is zero");
573  }
574 
575  centroid[idir] /= volume;
576  }
577 
578  return centroid;
579 }
580 
581 template <int dim> Real CutCellMoments<dim>::getResidual(const int & a_iDegree,
582  const int & a_normJ) const
583 {
584  Real retval;
585  if (isCovered() || isRegular())
586  {
587  retval = 0.0;
588  }
589  else
590  {
591  retval = m_residual[a_iDegree][a_normJ];
592  }
593 
594  return retval;
595 }
596 
597 template <int dim> void CutCellMoments<dim>::setResidual(const Real& a_value,
598  const int & a_iDegree,
599  const int & a_normJ)
600 {
601  if (isCovered() || isRegular())
602  {
603  MayDay::Abort("Invalid assignment to m_residual in covered or full cell");
604  }
605  else
606  {
607  m_residual[a_iDegree][a_normJ] = a_value;
608  }
609 }
610 
611 template <int dim> Vector<Real> CutCellMoments<dim>::sliceResidual(const int & a_iDegree) const
612 {
613  if (isCovered() || isRegular())
614  {
615  pout()<<"Dim = "<< dim << endl;
616  MayDay::Abort("Invalid attemtp to slice of m_residual in a covered or full CCM");
617  }
618  else
619  {
620  return m_residual[a_iDegree];
621  }
622 }
623 
624 template <int dim> bool CutCellMoments<dim>::isCovered() const
625 {
626  return m_IFData.m_allVerticesOut;
627 }
628 
629 template <int dim> bool CutCellMoments<dim>::isRegular() const
630 {
631  return (m_IFData.m_allVerticesIn && (!m_bdCCOn));
632 }
633 
634 template <int dim> void CutCellMoments<dim>::print(ostream& a_out) const
635 {
636  string padding = " ";
637  for (int i = 0; i < GLOBALDIM - dim; i++)
638  {
639  padding += " ";
640  }
641 
642  a_out << padding << "Dim = " << dim << "\n";
643  a_out << padding << "\n";
644 
645  for (typename PthMoment::const_iterator it = m_moments.begin(); it != m_moments.end(); ++it)
646  {
647  std::ios::fmtflags origFlags = a_out.flags();
648  int origWidth = a_out.width();
649  int origPrecision = a_out.precision();
650 
651  a_out << padding << "Integral "
652  << it->first
653  << " = "
654  << setw(23)
655  << setprecision(16)
656  << setiosflags(ios::showpoint)
657  << setiosflags(ios::scientific)
658  << it->second
659  << "\n";
660 
661  a_out.flags(origFlags);
662  a_out.width(origWidth);
663  a_out.precision(origPrecision);
664  }
665 
666  if (m_moments.size() > 0)
667  {
668  a_out << padding << "\n";
669  }
670 
671  for (typename PthMoment::const_iterator it = m_EBmoments.begin(); it != m_EBmoments.end(); ++it)
672  {
673  std::ios::fmtflags origFlags = a_out.flags();
674  int origWidth = a_out.width();
675  int origPrecision = a_out.precision();
676 
677  a_out << padding << "EBIntegral "
678  << it->first
679  << " = "
680  << setw(23)
681  << setprecision(16)
682  << setiosflags(ios::showpoint)
683  << setiosflags(ios::scientific)
684  << it->second
685  << "\n";
686 
687  a_out.flags(origFlags);
688  a_out.width(origWidth);
689  a_out.precision(origPrecision);
690  }
691 
692  if (m_EBmoments.size() > 0)
693  {
694  a_out << padding << "\n";
695  }
696 
697  a_out << padding << "Residuals:" << "\n";
698  for (int i = 0; i < m_residual.size(); i++)
699  {
700  a_out << padding << " " << i << ":";
701 
702  std::ios::fmtflags origFlags = a_out.flags();
703  int origWidth = a_out.width();
704  int origPrecision = a_out.precision();
705 
706  for (int j = 0; j < m_residual[i].size(); j++)
707  {
708 
709  a_out << " "
710  << setw(23)
711  << setprecision(16)
712  << setiosflags(ios::showpoint)
713  << setiosflags(ios::scientific)
714  << m_residual[i][j];
715  }
716 
717  a_out.flags(origFlags);
718  a_out.width(origWidth);
719  a_out.precision(origPrecision);
720 
721  a_out << padding << "\n";
722  }
723  a_out << padding << "\n";
724 
725  a_out << padding << "IFData:" << "\n";
726  a_out << padding << "\n";
727  a_out << m_IFData;
728 
729  a_out << padding << "Number of bdCutCellMoments = "
730  << m_bdCutCellMoments.size()
731  << "\n";
732  a_out << padding << "\n";
733 
734  for (typename BdCutCellMoments::const_iterator it = m_bdCutCellMoments.begin(); it != m_bdCutCellMoments.end(); ++it)
735  {
736  a_out << padding << "BdId = " << it->first << "\n";
737  a_out << it->second;
738  }
739 }
740 
741 template <int dim> void CutCellMoments<dim>::dump() const
742 {
743  print(pout());
744 }
745 
746 // Operators
747 template <int dim> void CutCellMoments<dim>::operator=(const CutCellMoments<dim> & a_cutCellMoments)
748 {
749  // Only copy if the objects are distinct
750  if (this != &a_cutCellMoments)
751  {
752  m_moments = a_cutCellMoments.m_moments;
753  m_EBmoments = a_cutCellMoments.m_EBmoments;
754  m_bdCutCellMoments = a_cutCellMoments.m_bdCutCellMoments;
755  m_IFData = a_cutCellMoments.m_IFData;
756  m_bdCCOn = a_cutCellMoments.m_bdCCOn;
757  m_residual = a_cutCellMoments.m_residual;
758  m_numActiveBounds = a_cutCellMoments.m_numActiveBounds;
759  m_badNormal = a_cutCellMoments.m_badNormal;
760  }
761 }
762 
763 template <int dim> ostream& operator<<(ostream & a_out,
764  const CutCellMoments<dim> & a_cutCellMoments)
765 {
766  a_cutCellMoments.print(a_out);
767  return a_out;
768 }
769 
770 #include "NamespaceFooter.H"
771 
772 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
Real convertDir(const Real &a_coord, const CoordinateSystem< dim > &a_system, const int &a_dir) const
Definition: CoordinateSystemImplem.H:91
#define BDID_HILO
Definition: Notation.H:89
bool isCovered() const
Definition: CutCellMomentsImplem.H:624
Real getVol(const EBorVol &a_EBorVol) const
Definition: CutCellMomentsImplem.H:542
BdCutCellMoments m_bdCutCellMoments
Definition: CutCellMoments.H:138
void addBdMoments(CutCellMoments< dim > &a_coarseCutCell, const IFData< dim+1 > &a_IFData, const int &a_degreePmax, const bool &a_useConstraints, const IndexTM< Real, dim > &a_refinedCenterDelta, const IndexTM< int, dim > &a_localHilo)
Definition: CutCellMomentsImplem.H:112
Definition: CoordinateSystem.H:34
#define LARGEREALVAL
Definition: Notation.H:77
Real getResidual(const int &a_iDegree, const int &a_normJ) const
Definition: CutCellMomentsImplem.H:581
void changeMomentCoordinatesToCellCenter()
Definition: CutCellMomentsImplem.H:367
#define LARGEINTVAL
Definition: Notation.H:76
EBorVol
Definition: Notation.H:101
map< IndexTM< int, dim-1 >, Real > PthMomentLesserDimension
Definition: CutCellMoments.H:40
Vector< Real > sliceResidual(const int &a_iDegree) const
Definition: CutCellMomentsImplem.H:611
Real getMoment(const IvDim &a_mono, const EBorVol &a_EBorVOL) const
Definition: CutCellMomentsImplem.H:501
bool m_badNormal
Definition: CutCellMoments.H:153
Definition: ComputeCutCellMoments.H:35
~CutCellMoments()
Definition: CutCellMomentsImplem.H:92
void print(ostream &out) const
Definition: CutCellMomentsImplem.H:634
Definition: IndexTM.H:36
CutCellMoments()
Definition: CutCellMomentsImplem.H:35
void setResidual(const Real &a_value, const int &a_iDegree, const int &a_normJ)
Definition: CutCellMomentsImplem.H:597
Definition: Notation.H:103
const PthMomentLoc & getMonomialLocMapDegreeP() const
Definition: LSProblem.H:94
IFData< dim > m_IFData
Definition: CutCellMoments.H:141
Vector< Vector< Real > > m_residual
Definition: CutCellMoments.H:147
Definition: IFData.H:42
const PthMomentLoc & getMonomialLocMapDegreePLess1() const
Definition: LSProblem.H:99
#define BDID_DIR
Definition: Notation.H:88
bool m_bdCCOn
Definition: CutCellMoments.H:144
double Real
Definition: REAL.H:33
Real fullCellQuadrature(const IndexTM< int, dim > &a_mono, const CoordinateSystem< dim > &a_coord)
Definition: CutCellMomentsImplem.H:237
bool isRegular() const
Definition: CutCellMomentsImplem.H:629
size_t size() const
Definition: Vector.H:192
bool m_allVerticesIn
Definition: IFData.H:71
Definition: Notation.H:104
void initialize(CutCellMoments< dim > &a_refinedCutCell)
Definition: CutCellMomentsImplem.H:411
static void Warning(const char *const a_msg=m_nullString)
Print out message to cerr and continue.
int m_numActiveBounds
Definition: CutCellMoments.H:150
void changeMomentCoordinatesToParentCenter()
Definition: CutCellMomentsImplem.H:389
PthMoment m_moments
Definition: CutCellMoments.H:132
map< IvDim, Real > PthMoment
Definition: CutCellMoments.H:38
bool m_allVerticesOut
Definition: IFData.H:72
Real getBdMoment(const IvDim &a_mono, const IFData< dim+1 > &a_IFData, const IndexTM< Real, dim > &a_refinedCenterDelta, PthMoment a_fullCellMap=PthMoment())
Definition: CutCellMomentsImplem.H:454
Real factorial(const int n)
Calculates factorial for an integer.
Definition: Factorial.H:25
void initializeMap(PthMoment &a_map1, PthMoment &a_map2)
Definition: CutCellMomentsImplem.H:444
RvDim getCentroid(const EBorVol &a_EBorVOL) const
Definition: CutCellMomentsImplem.H:553
int m_maxOrder
Definition: IFData.H:67
PthMoment m_EBmoments
Definition: CutCellMoments.H:135
const CutCellMoments< dim - 1 > getBdCutCellMoments(const Iv2 &a_bdId) const
Definition: CutCellMomentsImplem.H:96
void operator=(const CutCellMoments< dim > &a_cutCellMoments)
Definition: CutCellMomentsImplem.H:747
Real getBdEBMoment(const IvDim &a_mono, const IFData< dim+1 > &a_IFData, const IndexTM< Real, dim > &a_refinedCenterDelta)
Definition: CutCellMomentsImplem.H:480
void dump() const
Definition: CutCellMomentsImplem.H:741
map< IvDim, int > PthMomentLoc
Definition: CutCellMoments.H:44
int dim
Definition: EBInterface.H:146
ostream & operator<<(ostream &a_out, const CutCellMoments< dim > &a_cutCellMoments)
Definition: CutCellMomentsImplem.H:763
#define GLOBALDIM
Definition: Notation.H:35
Definition: CutCellMoments.H:32
static const IndexTM Zero
Definition: IndexTM.H:467
Real changeMomentCoordinates(PthMoment &a_refinedMomentMap, const IndexTM< int, dim > &a_monomial, const IndexTM< Real, dim > &a_refinedCenterDelta)
Definition: CutCellMomentsImplem.H:261
static void Abort(const char *const a_msg=m_nullString)
Print out message to cerr and exit via abort() (if serial) or MPI_Abort() (if parallel).