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