11 #ifndef _CUTCELLMOMENTSIMPLEM_H_
12 #define _CUTCELLMOMENTSIMPLEM_H_
14 #if defined(CH_Darwin) && defined(__GNUC__) && ( __GNUC__ == 3 )
17 #define _GLIBCPP_USE_C99 1
32 #include "NamespaceHeader.H"
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)
60 for (
int hilo = 0; hilo < 2; ++hilo)
62 for (
int idir = 0; idir <
dim; ++idir)
76 if (reducedInfo.m_allVerticesOn)
82 if (reducedInfo.m_badNormal)
98 typename BdCutCellMoments::const_iterator it = m_bdCutCellMoments.find(a_bdId);
100 if (it == m_bdCutCellMoments.end())
114 const int & a_degreePmax,
115 const bool & a_useConstraints,
128 for (
int iDegree = a_degreePmax; iDegree > 0; iDegree--)
135 for (
typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
138 fullCellMap[it->first] = fullCellQuadrature(it->first,m_IFData.m_cellCenterCoord);
142 for (
int iDegree = a_degreePmax; iDegree > 0; iDegree--)
151 for (
typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
153 Real bdMoment = getBdMoment(it->first,a_IFData,a_refinedCenterDelta,fullCellMap);
154 a_coarseBdCutCell.
m_moments[it->first] += bdMoment;
160 for (
int iDegree = a_degreePmax; iDegree > 0; iDegree--)
169 for (
typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
171 Real bdMoment = getBdMoment(it->first,a_IFData,a_refinedCenterDelta);
172 a_coarseBdCutCell.
m_moments[it->first] += bdMoment;
177 for (
int iDegree = a_degreePmax; iDegree >= 0; iDegree--)
185 for (
typename PthMomentLoc::const_iterator it = momMapP.begin(); it != momMapP.end(); ++it)
187 Real bdEBMoment = getBdEBMoment(it->first,a_IFData,a_refinedCenterDelta);
188 a_coarseBdCutCell.
m_EBmoments[it->first] += bdEBMoment;
194 for (
int idir = 0; idir <
dim; idir++)
204 for (
int jdir = 0; jdir <
dim; jdir++)
208 localRefinedCellCenter[jdir] = a_refinedCenterDelta[jdir];
209 localHilo[jdir] = a_hilo[jdir];
211 else if (jdir > idir)
213 localRefinedCellCenter[jdir-1] = a_refinedCenterDelta[jdir];
214 localHilo[jdir-1] = a_hilo[jdir];
220 bdId[1] = a_hilo[idir];
224 (m_bdCutCellMoments[bdId]).addBdMoments(a_coarseBdCutCell.
m_bdCutCellMoments[bdId],m_IFData,a_degreePmax+1,a_useConstraints,localRefinedCellCenter,localHilo);
228 for (
int side = 0; side < 2; side++)
231 (m_bdCutCellMoments[bdId]).addBdMoments(a_coarseBdCutCell.
m_bdCutCellMoments[bdId],m_IFData,a_degreePmax,a_useConstraints,localRefinedCellCenter,localHilo);
242 for (
int idir = 0; idir <
dim; ++idir)
244 Real loPt = a_coord.
convertDir(-0.5*m_IFData.m_cellCenterCoord.m_dx[idir],
245 m_IFData.m_cellCenterCoord,
248 Real hiPt = a_coord.
convertDir(0.5*m_IFData.m_cellCenterCoord.m_dx[idir],
249 m_IFData.m_cellCenterCoord,
252 Real idirInt = pow(hiPt,a_mono[idir] + 1) - pow(loPt,a_mono[idir] + 1);
254 idirInt /= (a_mono[idir] + 1);
275 for (
int i = 0; i <
dim; i++)
277 degree += a_monomial[i];
283 for (
int r = 0; r <= degree; r++)
290 for (
int idir = 0; idir <
dim; ++idir)
292 derivative[idir] = 0;
297 for (
int j = 1; j < dim-1; ++j)
300 for (
int k = j+1; k <
dim; ++k)
305 if (derivative[j] > t)
307 derivative[j+1] += 1;
316 if (derivative[dim-1] > r)
323 for (
int j = 1; j <
dim; ++j)
325 derivative[0] -= derivative[j];
332 for (
int idir = 0; idir <
dim; idir++)
334 for (
int j = 0; j < derivative[idir]; j++)
336 coeff *= a_monomial[idir] - j;
342 if (a_monomial[idir]-derivative[idir] < 0)
346 MayDay::Abort(
"Coeff should be zero when the derivative has higher coefficient than the monomial");
351 coeff *= pow(a_refinedCenterDelta[idir],a_monomial[idir]-derivative[idir]);
356 moment += coeff * a_refinedMomentMap[derivative] /
factorial;
371 delta -= m_IFData.m_parentCoord .m_origin;
374 for (
typename PthMoment::const_iterator it = copyMoments.begin();
375 it != copyMoments.end(); ++it)
377 IvDim mono = it->first;
378 m_moments[mono] = changeMomentCoordinates(copyMoments, mono, delta);
382 for (
typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
384 IvDim mono = it->first;
385 m_EBmoments[mono] = changeMomentCoordinates(copyEBMoments, mono, delta);
393 delta -= m_IFData.m_cellCenterCoord.m_origin;
396 for (
typename PthMoment::const_iterator it = copyMoments.begin();
397 it != copyMoments.end(); ++it)
399 IvDim mono = it->first;
400 m_moments[mono] = changeMomentCoordinates(copyMoments, mono, delta);
404 for (
typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
406 IvDim mono = it->first;
407 m_EBmoments[mono] = changeMomentCoordinates(copyEBMoments, mono, delta);
415 initializeMap(m_EBmoments,a_refinedCutCell.
m_EBmoments);
416 initializeMap(m_moments,a_refinedCutCell.
m_moments);
421 for (
int idir = 0; idir <
dim; idir++)
425 for (
int hilo = 0; hilo < 2; hilo++)
430 m_bdCutCellMoments[bdId].initialize(refinedBdCutCell);
438 for (
typename PthMomentLesserDimension::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
440 a_map1[it->first] = 0.0;
447 for (
typename PthMoment::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
449 a_map1[it->first] = 0.0;
470 moment = changeMomentCoordinates(a_fullCellMap,a_mono,a_refinedCenterDelta);
474 moment = changeMomentCoordinates(m_moments,a_mono,a_refinedCenterDelta);
494 EBmoment = changeMomentCoordinates(m_EBmoments,a_mono,a_refinedCenterDelta);
502 const EBorVol & a_EBorVol)
const
509 typename PthMoment::const_iterator it = m_moments.find(a_mono);
511 if (it != m_moments.end())
523 typename PthMoment::const_iterator it = m_EBmoments.find(a_mono);
525 if (it != m_EBmoments.end())
547 IvDim zeroIv = IvDim::Zero;
548 volume = getMoment(zeroIv,a_EBorVol);
556 Real volume = getVol(a_EBorVol);
561 for (
int idir = 0; idir <
dim; ++idir)
563 IvDim mono = BASISV_TM<int,dim>(idir);
564 centroid[idir] = getMoment(mono,a_EBorVol);
568 MayDay::Warning(
"CutCellMoments::getCentroid: Volume fraction is negative");
570 else if (volume == 0.0)
575 centroid[idir] /= volume;
582 const int & a_normJ)
const
585 if (isCovered() || isRegular())
591 retval = m_residual[a_iDegree][a_normJ];
598 const int & a_iDegree,
601 if (isCovered() || isRegular())
603 MayDay::Abort(
"Invalid assignment to m_residual in covered or full cell");
607 m_residual[a_iDegree][a_normJ] = a_value;
613 if (isCovered() || isRegular())
615 pout()<<
"Dim = "<<
dim << endl;
616 MayDay::Abort(
"Invalid attemtp to slice of m_residual in a covered or full CCM");
620 return m_residual[a_iDegree];
626 return m_IFData.m_allVerticesOut;
631 return (m_IFData.m_allVerticesIn && (!m_bdCCOn));
636 string padding =
" ";
642 a_out << padding <<
"Dim = " << dim <<
"\n";
643 a_out << padding <<
"\n";
645 for (
typename PthMoment::const_iterator it = m_moments.begin(); it != m_moments.end(); ++it)
647 std::ios::fmtflags origFlags = a_out.flags();
648 int origWidth = a_out.width();
649 int origPrecision = a_out.precision();
651 a_out << padding <<
"Integral "
656 << setiosflags(ios::showpoint)
657 << setiosflags(ios::scientific)
661 a_out.flags(origFlags);
662 a_out.width(origWidth);
663 a_out.precision(origPrecision);
666 if (m_moments.size() > 0)
668 a_out << padding <<
"\n";
671 for (
typename PthMoment::const_iterator it = m_EBmoments.begin(); it != m_EBmoments.end(); ++it)
673 std::ios::fmtflags origFlags = a_out.flags();
674 int origWidth = a_out.width();
675 int origPrecision = a_out.precision();
677 a_out << padding <<
"EBIntegral "
682 << setiosflags(ios::showpoint)
683 << setiosflags(ios::scientific)
687 a_out.flags(origFlags);
688 a_out.width(origWidth);
689 a_out.precision(origPrecision);
692 if (m_EBmoments.size() > 0)
694 a_out << padding <<
"\n";
697 a_out << padding <<
"Residuals:" <<
"\n";
698 for (
int i = 0; i < m_residual.size(); i++)
700 a_out << padding <<
" " << i <<
":";
702 std::ios::fmtflags origFlags = a_out.flags();
703 int origWidth = a_out.width();
704 int origPrecision = a_out.precision();
706 for (
int j = 0; j < m_residual[i].size(); j++)
712 << setiosflags(ios::showpoint)
713 << setiosflags(ios::scientific)
717 a_out.flags(origFlags);
718 a_out.width(origWidth);
719 a_out.precision(origPrecision);
721 a_out << padding <<
"\n";
723 a_out << padding <<
"\n";
725 a_out << padding <<
"IFData:" <<
"\n";
726 a_out << padding <<
"\n";
729 a_out << padding <<
"Number of bdCutCellMoments = "
730 << m_bdCutCellMoments.size()
732 a_out << padding <<
"\n";
734 for (
typename BdCutCellMoments::const_iterator it = m_bdCutCellMoments.begin(); it != m_bdCutCellMoments.end(); ++it)
736 a_out << padding <<
"BdId = " << it->first <<
"\n";
750 if (
this != &a_cutCellMoments)
755 m_IFData = a_cutCellMoments.
m_IFData;
756 m_bdCCOn = a_cutCellMoments.
m_bdCCOn;
766 a_cutCellMoments.
print(a_out);
770 #include "NamespaceFooter.H"
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: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
void changeMomentCoordinatesToCellCenter()
Definition: CutCellMomentsImplem.H:367
RvDim getCentroid(const EBorVol &a_EBorVOL) const
Definition: CutCellMomentsImplem.H:553
#define LARGEINTVAL
Definition: Notation.H:76
EBorVol
Definition: Notation.H:101
map< IndexTM< int, dim-1 >, Real > PthMomentLesserDimension
Definition: CutCellMoments.H:40
bool m_badNormal
Definition: CutCellMoments.H:153
void dump() const
Definition: CutCellMomentsImplem.H:741
Definition: ComputeCutCellMoments.H:35
~CutCellMoments()
Definition: CutCellMomentsImplem.H:92
bool isCovered() const
Definition: CutCellMomentsImplem.H:624
Vector< Real > sliceResidual(const int &a_iDegree) const
Definition: CutCellMomentsImplem.H:611
void print(ostream &out) const
Definition: CutCellMomentsImplem.H:634
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
IFData< dim > m_IFData
Definition: CutCellMoments.H:141
Vector< Vector< Real > > m_residual
Definition: CutCellMoments.H:147
#define BDID_DIR
Definition: Notation.H:88
bool m_bdCCOn
Definition: CutCellMoments.H:144
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:96
Real fullCellQuadrature(const IndexTM< int, dim > &a_mono, const CoordinateSystem< dim > &a_coord)
Definition: CutCellMomentsImplem.H:237
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
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:389
PthMoment m_moments
Definition: CutCellMoments.H:132
Real getMoment(const IvDim &a_mono, const EBorVol &a_EBorVOL) const
Definition: CutCellMomentsImplem.H:501
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
Real getVol(const EBorVol &a_EBorVol) const
Definition: CutCellMomentsImplem.H:542
const PthMomentLoc & getMonomialLocMapDegreePLess1() const
Definition: LSProblem.H:99
int m_maxOrder
Definition: IFData.H:67
PthMoment m_EBmoments
Definition: CutCellMoments.H:135
void operator=(const CutCellMoments< dim > &a_cutCellMoments)
Definition: CutCellMomentsImplem.H:747
bool isRegular() const
Definition: CutCellMomentsImplem.H:629
Real getBdEBMoment(const IvDim &a_mono, const IFData< dim+1 > &a_IFData, const IndexTM< Real, dim > &a_refinedCenterDelta)
Definition: CutCellMomentsImplem.H:480
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
Real changeMomentCoordinates(PthMoment &a_refinedMomentMap, const IndexTM< int, dim > &a_monomial, const IndexTM< Real, dim > &a_refinedCenterDelta)
Definition: CutCellMomentsImplem.H:261
Real getResidual(const int &a_iDegree, const int &a_normJ) const
Definition: CutCellMomentsImplem.H:581
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).