Chombo + EB  3.2
ComputeCutCellMomentsImplem.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 _COMPUTECUTCELLMOMENTSIMPLEM_H_
12 #define _COMPUTECUTCELLMOMENTSIMPLEM_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 "NoRefinement.H"
31 #include "LSProblem.H"
32 
33 #include "NamespaceHeader.H"
34 
35 // Null constructor
37 {
38 }
39 
40 // Copy constructor
41 template <int dim> ComputeCutCellMoments<dim>::ComputeCutCellMoments(const ComputeCutCellMoments<dim>& a_computeCutCellMoments)
42  :m_cutCellMoments(a_computeCutCellMoments.m_cutCellMoments)
43 {
44 }
45 
46 // Use this for initializing
48  :m_cutCellMoments(a_info)
49 {
50  m_cutCellMoments.m_bdCCOn = false;
51 
52  for (int hilo = 0; hilo < 2; ++hilo)
53  {
54  for (int idir = 0; idir < dim; ++idir)
55  {
56  // Identifier for which boundary cutCellMoment
57  Iv2 bdId;
58  bdId[BDID_DIR] = idir;
59  bdId[BDID_HILO] = hilo;
60 
61 #if RECURSIVE_GEOMETRY_GENERATION == 0
62  IFData<dim-1> reducedInfo(a_info,idir,hilo);
63 #else
64  IFData<dim-1> reducedInfo(a_info,a_info.m_maxOrder+1,idir,hilo);
65 #endif
66  CutCellMoments<dim-1>bdCutCellMoments(reducedInfo);
67 
68  m_cutCellMoments.m_bdCutCellMoments[bdId] = bdCutCellMoments;
69 
70 
71  // Notice whether at least one lower dimensional cutCell is on the
72  // interface
73  if (reducedInfo.m_allVerticesOn)
74  {
75  m_cutCellMoments.m_bdCCOn = true;
76  }
77  }
78  }
79 }
80 
81 // Destructor
83 {
84 }
85 
86 #if RECURSIVE_GEOMETRY_GENERATION == 0
87 template <int dim> void ComputeCutCellMoments<dim>::computeMoments(const int & a_order,
88  const int & a_degreeP,
89  const bool & a_useConstraints,
90  RefinementCriterion<dim> & a_refinementCriterion,
91  const int & a_numberOfRefinements)
92 {
93  CH_assert(m_cutCellMoments.m_IFData.m_maxOrder >= a_order);
94 
95  Vector<Real> RNorm(3);
96  for (int i = 0; i < 3; i++)
97  {
98  RNorm[i] = LARGEREALVAL;
99  }
100 
101  for (int i = 0; i < a_degreeP + 1; i++)
102  {
103  m_cutCellMoments.m_residual.push_back(RNorm);
104  }
105 
106  if (m_cutCellMoments.m_IFData.m_allVerticesOut)
107  {
108  for (int iDegree = a_degreeP; iDegree >= 0; --iDegree)
109  {
110  LSProblem<dim> lsp(iDegree,a_useConstraints);
111 
112  // Fill moments
113  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
114  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
115  {
116  m_cutCellMoments.m_EBmoments[it->first] = 0.0;
117  }
118 
119  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
120  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
121  {
122  m_cutCellMoments.m_moments[it->first] = 0.0;
123  }
124  }
125  }
126  else if (m_cutCellMoments.m_IFData.m_allVerticesIn && !m_cutCellMoments.m_bdCCOn)
127  {
128  for (int iDegree = a_degreeP; iDegree >= 0; --iDegree)
129  {
130  LSProblem<dim> lsp(iDegree,a_useConstraints);
131 
132  // Fill moments of degree P and P-1
133  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
134  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
135  {
136  m_cutCellMoments.m_EBmoments[it->first] = 0.0;
137  }
138 
139  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
140  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
141  {
142  m_cutCellMoments.m_moments[it->first] = m_cutCellMoments.fullCellQuadrature(it->first,m_cutCellMoments.m_IFData.m_parentCoord);
143  }
144  }
145  }
146  else
147  {
148  for (typename BdCutCellMoments::iterator it = m_cutCellMoments.m_bdCutCellMoments.begin(); it != m_cutCellMoments.m_bdCutCellMoments.end(); ++it)
149  {
150  NoRefinement<dim-1> noRefinement;
151 
152  ComputeCutCellMoments<dim-1> subProblem(it->second.m_IFData);
153  subProblem.computeMoments(a_order,a_degreeP+1,a_useConstraints,noRefinement);
154  it->second = subProblem.m_cutCellMoments;
155  //pout()<< "Subproblem constraints: " << it->second.m_numActiveBounds << endl;
156  }
157 
158  for (int iDegree = a_degreeP; iDegree >= 0; --iDegree)
159  {
160  // Make a LSProb
161  IvDim zeroDerivative = IndexTM<int,dim>::Zero;
162  LSProblem<dim> lsp(a_order,iDegree,a_useConstraints,m_cutCellMoments.m_IFData.m_normalDerivatives[zeroDerivative]);
163 
164  Vector<Real> rhs = computeRhs(lsp,a_order);
165 
166  if (a_useConstraints)
167  {
168  lsp.computeBounds(m_cutCellMoments.m_IFData.m_globalCoord.m_dx,m_cutCellMoments);
169  }
170 
171 
172  // Solve the problem and return residual
173  int lsCode = lsp.invertNormalEq(rhs,m_cutCellMoments.m_residual[iDegree]);
174  if (lsCode != 0)
175  {
176  pout() << "Geometry generation least squares problem failed with residual: "
177  << m_cutCellMoments.m_residual[iDegree]<< endl;
178  lsp.print(pout());
179 
180  pout () << endl << "Problem occurred generating geometry for these cut cell moments:\n"
181  << m_cutCellMoments << endl;
182 
183  //MayDay::Error("Geometry generation error.");
184  bool status = false;
185  a_refinementCriterion.setConstrantSuccessStatus(status);
186  }
187 
188  // Record number of constraint violations
189  m_cutCellMoments.m_numActiveBounds += lsp.numActiveBounds();
190 
191  //if (m_cutCellMoments.m_numActiveBounds > 0)
192  // {
193  // pout() <<" Greater than zero. Dim= " << dim << " num="<<m_cutCellMoments.m_numActiveBounds << endl;
194  // }
195 
196  // Fill moments
197  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
198  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin();
199  it != monDegreePLess1.end(); ++it)
200  {
201  m_cutCellMoments.m_moments[it->first] = lsp.getUnknown(it->second + lsp.getNumberDegP());
202  }
203 
204  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
205  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
206  {
207  m_cutCellMoments.m_EBmoments[it->first] = lsp.getUnknown(it->second);
208  }
209  }
210 
211  // Move moments from local coord to parent coord
212  IndexTM<Real,dim> delta = m_cutCellMoments.m_IFData.m_parentCoord.m_origin;
213  delta -= m_cutCellMoments.m_IFData.m_localCoord .m_origin;
214 
215  PthMoment copyMoments = m_cutCellMoments.m_moments;
216  for (typename PthMoment::const_iterator it = copyMoments.begin();
217  it != copyMoments.end(); ++it)
218  {
219  IvDim mono = it->first;
220  m_cutCellMoments.m_moments[mono] = m_cutCellMoments.changeMomentCoordinates(copyMoments, mono, delta);
221  }
222 
223  PthMoment copyEBMoments = m_cutCellMoments.m_EBmoments;
224  for (typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
225  {
226  IvDim mono = it->first;
227  m_cutCellMoments.m_EBmoments[mono] = m_cutCellMoments.changeMomentCoordinates(copyEBMoments, mono, delta);
228  }
229 
230  // Move bdCutCell moments from parent coord to cell center coord. From here on bdCutCell moments will
231  // not be used in any least squares problem
232 
233  for (typename BdCutCellMoments::iterator it = m_cutCellMoments.m_bdCutCellMoments.begin(); it != m_cutCellMoments.m_bdCutCellMoments.end(); ++it)
234  {
235  it->second.changeMomentCoordinatesToCellCenter();
236  }
237  }
238 
239  IndexTM<int,dim> refineInDir;
240 
241  if (!m_cutCellMoments.isCovered() && !m_cutCellMoments.isRegular() && a_refinementCriterion.baseDoRefine(refineInDir,m_cutCellMoments,a_numberOfRefinements))
242  {
243  // Refine creates a vector of 2^dim cutcellmoments corresponding to
244  // the refined cells inside the current cell
245  Vector<CutCellMoments<dim> > refinedCCMoms = refine(a_order,a_degreeP,a_useConstraints,a_refinementCriterion,refineInDir,a_numberOfRefinements);
246 
247  // addMomentMaps add up the refined moments in the current
248  // CutCellMoments maps
249  addMomentMaps(refinedCCMoms,a_degreeP,a_useConstraints);
250 
251  // computeResiduals(a_order,a_degreeP,a_useConstraints);
252  computeResiduals(refinedCCMoms,a_degreeP);
253 
254 
255  // Move the origin of moments from cell center to parent center.
256  IndexTM<Real,dim> delta = m_cutCellMoments.m_IFData.m_parentCoord .m_origin;
257  delta -= m_cutCellMoments.m_IFData.m_cellCenterCoord.m_origin;
258 
259  PthMoment copyMoments = m_cutCellMoments.m_moments;
260  for (typename PthMoment::const_iterator it = copyMoments.begin();
261  it != copyMoments.end(); ++it)
262  {
263  IvDim mono = it->first;
264  m_cutCellMoments.m_moments[mono] = m_cutCellMoments.changeMomentCoordinates(copyMoments, mono, delta);
265  }
266 
267  PthMoment copyEBMoments = m_cutCellMoments.m_EBmoments;
268  for (typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
269  {
270  IvDim mono = it->first;
271  m_cutCellMoments.m_EBmoments[mono] = m_cutCellMoments.changeMomentCoordinates(copyEBMoments, mono, delta);
272  }
273  }
274 }
275 #else
276 // Solve
277 template <int dim> void ComputeCutCellMoments<dim>::computeMoments(const int & a_orderPmax,
278  const int & a_degreePmax,
279  const bool & a_useConstraints,
280  RefinementCriterion<dim> & a_refinementCriterion,
281  const int & a_numberOfRefinements)
282 {
283  int zerothOrderOfAccuracy = 0;
284  int highestDegree = a_orderPmax + a_degreePmax;
285 
286  Vector<Real> RNorm(3);
287  for (int i = 0; i < 3; i++)
288  {
289  RNorm[i] = LARGEREALVAL;
290  }
291 
292  for (int i = 0; i <= highestDegree; i++)
293  {
294  m_cutCellMoments.m_residual.push_back(RNorm);
295  }
296 
297  m_boundaryMomentsComputed = false;
298 
299  computeMomentsRecursively(zerothOrderOfAccuracy,
300  highestDegree,
301  a_useConstraints,
302  a_refinementCriterion,
303  a_numberOfRefinements);
304 
305  IndexTM<int,dim> refineInDir;
306 
307  if (!m_cutCellMoments.isCovered() && !m_cutCellMoments.isRegular() && a_refinementCriterion.doRefine(refineInDir,m_cutCellMoments,a_numberOfRefinements))
308  {
309  // Refine creates a vector of 2^dim cutcellmoments corresponding to
310  // the refined cells inside the current cell
311  Vector<CutCellMoments<dim> > refinedCCMoms = refine(a_orderPmax,a_degreePmax,a_useConstraints,a_refinementCriterion,refineInDir,a_numberOfRefinements);
312 
313  // addMomentMaps add up the refined moments in the current
314  // CutCellMoments maps
315  addMomentMaps(refinedCCMoms,a_degreePmax,a_useConstraints);
316 
317  // computeResiduals(a_orderPmax,a_degreePmax,a_useConstraints);
318  computeResiduals(refinedCCMoms,a_degreePmax);
319 
320  // Move the origin of moments from cell center to parent center.
321  IndexTM<Real,dim> delta = m_cutCellMoments.m_IFData.m_parentCoord.m_origin;
322  delta -= m_cutCellMoments.m_IFData.m_cellCenterCoord .m_origin;
323 
324  PthMoment copyMoments = m_cutCellMoments.m_moments;
325  for (typename PthMoment::const_iterator it = copyMoments.begin();
326  it != copyMoments.end(); ++it)
327  {
328  IvDim mono = it->first;
329  m_cutCellMoments.m_moments[mono] = m_cutCellMoments.changeMomentCoordinates(copyMoments, mono, delta);
330  }
331 
332  PthMoment copyEBMoments = m_cutCellMoments.m_EBmoments;
333  for (typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
334  {
335  IvDim mono = it->first;
336  m_cutCellMoments.m_EBmoments[mono] = m_cutCellMoments.changeMomentCoordinates(copyEBMoments, mono, delta);
337  }
338 
339  // Move bdCutCell moments from parent coord to cell center coord. From here on bdCutCell moments will
340  // not be used in any least squares problem
341 
342  for (typename BdCutCellMoments::iterator it = m_cutCellMoments.m_bdCutCellMoments.begin(); it != m_cutCellMoments.m_bdCutCellMoments.end(); ++it)
343  {
344  it->second.changeMomentCoordinatesToCellCenter();
345  }
346  }
347 
348  // Move moments from local coord to parent coord
349  IndexTM<Real,dim> delta = m_cutCellMoments.m_IFData.m_parentCoord.m_origin;
350  delta -= m_cutCellMoments.m_IFData.m_localCoord .m_origin;
351 
352  PthMoment copyMoments = m_cutCellMoments.m_moments;
353  for (typename PthMoment::const_iterator it = copyMoments.begin();
354  it != copyMoments.end(); ++it)
355  {
356  IvDim mono = it->first;
357  m_cutCellMoments.m_moments[mono] = m_cutCellMoments.changeMomentCoordinates(copyMoments, mono, delta);
358  }
359 
360  PthMoment copyEBMoments = m_cutCellMoments.m_EBmoments;
361  for (typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
362  {
363  IvDim mono = it->first;
364  m_cutCellMoments.m_EBmoments[mono] = m_cutCellMoments.changeMomentCoordinates(copyEBMoments, mono, delta);
365  }
366 }
367 
368 template <int dim> void ComputeCutCellMoments<dim>::computeMomentsRecursively(const int & a_orderPmax,
369  const int & a_degreePmax,
370  const bool & a_useConstraints,
371  RefinementCriterion<dim> & a_refinementCriterion,
372  const int & a_numberOfRefinements)
373 {
374  CH_assert(m_cutCellMoments.m_IFData.m_maxOrder >= a_orderPmax);
375 
376  if (m_cutCellMoments.m_IFData.m_allVerticesOut)
377  {
378  for (int iOrder = 0; iOrder <= a_orderPmax; ++iOrder)
379  {
380  LSProblem<dim> lsp(iOrder + a_degreePmax, a_useConstraints);
381 
382  // Fill moments
383  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
384  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
385  {
386  m_cutCellMoments.m_EBmoments[it->first] = 0.0;
387  }
388 
389  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
390  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
391  {
392  m_cutCellMoments.m_moments[it->first] = 0.0;
393  }
394  }
395  }
396  else if (m_cutCellMoments.m_IFData.m_allVerticesIn && !m_cutCellMoments.m_bdCCOn)
397  {
398  for (int iOrder = 0; iOrder <= a_orderPmax; ++iOrder)
399  {
400  LSProblem<dim> lsp(iOrder + a_degreePmax, a_useConstraints);
401 
402  // Fill moments of degree P and P-1
403  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
404  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
405  {
406  m_cutCellMoments.m_EBmoments[it->first] = 0.0;
407  }
408 
409  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
410  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
411  {
412  m_cutCellMoments.m_moments[it->first] = m_cutCellMoments.fullCellQuadrature(it->first,m_cutCellMoments.m_IFData.m_parentCoord);
413  }
414  }
415  }
416  else
417  {
418  // Only compute the boundary moments if they haven't already been
419  // computed (earlier in the recursion)
420  if (!m_boundaryMomentsComputed)
421  {
422  for (typename BdCutCellMoments::iterator it = m_cutCellMoments.m_bdCutCellMoments.begin(); it != m_cutCellMoments.m_bdCutCellMoments.end(); ++it)
423  {
424  NoRefinement<dim-1> noRefinement;
425 
426  ComputeCutCellMoments<dim-1> subProblem(it->second.m_IFData);
427  subProblem.computeMoments(a_orderPmax,a_degreePmax+1,a_useConstraints,noRefinement);
428  it->second = subProblem.m_cutCellMoments;
429  }
430 
431  m_boundaryMomentsComputed = true;
432  }
433 
434  // Make a LSProb
435  IvDim zeroDerivative = IndexTM<int,dim>::Zero;
436  LSProblem<dim> lsp(a_orderPmax,a_degreePmax,a_useConstraints,m_cutCellMoments.m_IFData.m_normalDerivatives[zeroDerivative]);
437 
438  Vector<Real> rhs = computeRhs(lsp,a_orderPmax);
439 
440  if (a_useConstraints)
441  {
442  lsp.computeBounds(m_cutCellMoments.m_IFData.m_globalCoord.m_dx,m_cutCellMoments);
443  }
444 
445  // Solve the problem and return residual
446  int lsCode=lsp.invertNormalEq(rhs,m_cutCellMoments.m_residual[a_degreePmax]);
447  if (lsCode != 0)
448  {
449  pout() << "Geometry generation least squares problem failed with residual:"
450  << m_cutCellMoments.m_residual[a_degreePmax]<< endl;
451  lsp.print(pout());
452  pout () << "Problem occurred generating geometry for these cut cell moments: " << m_cutCellMoments << endl;
453  MayDay::Error("Geometry generation error.[2]");
454  }
455 
456 
457  // Fill moments
458  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
459  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin();
460  it != monDegreePLess1.end(); ++it)
461  {
462  m_cutCellMoments.m_moments[it->first] = lsp.getUnknown(it->second + lsp.getNumberDegP());
463  }
464 
465  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
466  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
467  {
468  m_cutCellMoments.m_EBmoments[it->first] = lsp.getUnknown(it->second);
469  }
470 
471  }
472 
473  if (a_degreePmax > 0)
474  {
475  computeMomentsRecursively(a_orderPmax + 1,
476  a_degreePmax - 1,
477  a_useConstraints,
478  a_refinementCriterion,
479  a_numberOfRefinements);
480 
481  }
482 }
483 #endif
484 
487  const int & a_order)
488 #else
489  const int & a_orderPmax)
490 #endif
491 {
492  // Resize rhs
493  int numEq = dim*a_lsp.getNumberDegP();
494  Vector<Real> rhs(numEq);
495 
496  // For each moment iterate thru bd CutCellMoments incrementing same comp
497  // of rhs
498  const LocPthMoment& locMap = a_lsp.getLocMonomialMapDegreeP();
499  for (typename LocPthMoment::const_iterator it = locMap.begin(); it != locMap.end(); ++it)
500  {
501  int jth = it->first;
502  IvDim mono = it->second;
503 
504  int hiSide = 1;
505  int loSide = 0;
506  Iv2 bdId;
507 
508  for (int idir = 0; idir < dim; ++idir)
509  {
510  // Which lower dimensional monomial corresponds (mono,j)
511  IndexTM<int,dim-1> mono1Less;
512  for (int jdir = 0; jdir < dim; ++jdir)
513  {
514  if (jdir < idir)
515  {
516  mono1Less[jdir] = mono[jdir];
517  }
518  else if (jdir > idir)
519  {
520  mono1Less[jdir-1] = mono[jdir];
521  }
522  }
523 
524  bdId[0] = idir;
525  bdId[1] = hiSide;
526 
527  Real hiMom = m_cutCellMoments.m_bdCutCellMoments[bdId].m_moments[mono1Less];
528 
529  bdId[1] = loSide;
530 
531  Real loMom = m_cutCellMoments.m_bdCutCellMoments[bdId].m_moments[mono1Less];
532  int exponent = it->second[idir];
533 
534  Real loSideValue;
535  Real hiSideValue;
536 
537  loSideValue = m_cutCellMoments.m_IFData.m_localCoord.convertDir(-0.5*m_cutCellMoments.m_IFData.m_cellCenterCoord.m_dx[idir],
538  m_cutCellMoments.m_IFData.m_cellCenterCoord,
539  idir);
540 
541  hiSideValue = m_cutCellMoments.m_IFData.m_localCoord.convertDir( 0.5*m_cutCellMoments.m_IFData.m_cellCenterCoord.m_dx[idir],
542  m_cutCellMoments.m_IFData.m_cellCenterCoord,
543  idir);
544 
545  Real loFactor = pow(loSideValue,exponent);
546  Real hiFactor = pow(hiSideValue,exponent);
547 
548  rhs[(dim*jth) + idir] = hiMom*hiFactor - loMom*loFactor;
549 
550  // Add the Taylor series terms
551 #if RECURSIVE_GEOMETRY_GENERATION == 0
552  for (int order = 1; order <= a_order; order++)
553 #else
554  for (int order = 1; order <= a_orderPmax; order++)
555 #endif
556  {
557  Vector<IvDim> taylorMonomials;
558 
559  generateMultiIndices(taylorMonomials,order);
560 
561  for (int i = 0; i < taylorMonomials.size(); i++)
562  {
563  const IvDim & taylorMonomial = taylorMonomials[i];
564 
565  IvDim totalMonomial = mono + taylorMonomial;
566 
567  if (m_cutCellMoments.m_EBmoments.find(totalMonomial) !=
568  m_cutCellMoments.m_EBmoments.end())
569  {
570  Real normalDerivative = m_cutCellMoments.m_IFData.m_normalDerivatives[taylorMonomial][idir];
571  Real fact = factorial(taylorMonomial);
572 
573  Real moment = m_cutCellMoments.m_EBmoments[totalMonomial];
574 
575  rhs[(dim*jth) + idir] += normalDerivative * moment / fact;
576  }
577 #if RECURSIVE_GEOMETRY_GENERATION != 0
578  else
579  {
580  MayDay::Error("Unable to find needed monomial for Taylor series");
581  }
582 #endif
583  }
584  }
585  }
586  }
587 
588  return rhs;
589 }
590 
591 #if RECURSIVE_GEOMETRY_GENERATION == 0
592 template <int dim> Vector< CutCellMoments<dim> > ComputeCutCellMoments<dim>::refine(const int & a_order,
593  const int & a_degreeP,
594 #else
595 template <int dim> Vector< CutCellMoments<dim> > ComputeCutCellMoments<dim>::refine(const int & a_orderPmax,
596  const int & a_degreePmax,
597 #endif
598  const bool & a_useConstraints,
599  RefinementCriterion<dim> & a_refinementCriterion,
600  const IndexTM<int,dim> & a_refineInDir,
601  const int & a_numberOfRefinements)
602 {
603  Vector<CutCellMoments<dim> > refinedCutCellMoments;
604 
605  // This function refines the current cell by a factor of two in all
606  // directions where a_refineInDir is non-zero
607  IndexTM<Real,dim> refinedDx = m_cutCellMoments.m_IFData.m_globalCoord.m_dx;
608 
609  // Current cell iterator and maximum iterator
610  IndexTM<int,dim> curIter;
611  IndexTM<int,dim> maxIter;
612 
613  for (int idir = 0; idir < dim; idir++)
614  {
615  if (a_refineInDir[idir] != 0)
616  {
617  // This direction will be divided in two
618  refinedDx[idir] /= 2.0;
619  maxIter[idir] = 1;
620  }
621  else
622  {
623  // This direction is undivided
624  maxIter[idir] = 0;
625  }
626 
627  // Initial the start to zero
628  curIter[idir] = 0;
629  }
630 
631  // Iterate through all the subcells
632  while (1)
633  {
634  // Computes the cell center for the iCell refined cell
635  IndexTM<Real,dim> refinedCellCenter = m_cutCellMoments.m_IFData.m_globalCoord.m_origin;
636 
637  for (int idir = 0; idir < dim; idir++)
638  {
639  if (a_refineInDir[idir] != 0)
640  {
641  int sign = 2*curIter[idir] - 1;
642  refinedCellCenter[idir] += sign * 0.5*refinedDx[idir];
643  }
644  }
645 
646  // Creates the IF Data for the refined cell, it uses the global
647  // constructor of IFData and recomputes all the intersection points
648  // and normal/grad normal from the implicit function
649 #if RECURSIVE_GEOMETRY_GENERATION == 0
650  IFData<dim> refinedIFData(m_cutCellMoments.m_IFData.m_function,refinedDx,refinedCellCenter,a_order);
651 #else
652  IFData<dim> refinedIFData(m_cutCellMoments.m_IFData.m_function,refinedDx,refinedCellCenter,a_orderPmax);
653 #endif
654 #if 0
655  pout() << "refinedIFData" << "\n";
656  pout() << refinedIFData<<endl;
657 #endif
658  // Creates a cutcellmoments for the refined cell
659  ComputeCutCellMoments<dim> refinedComputeCCMom(refinedIFData);
660 
661  // Track current refinement level
662  int numberOfRefinements = a_numberOfRefinements + 1;
663 
664  // Compute the moments for the refined cell
665 #if RECURSIVE_GEOMETRY_GENERATION == 0
666  refinedComputeCCMom.computeMoments(a_order,a_degreeP,a_useConstraints,a_refinementCriterion,numberOfRefinements);
667 #else
668  refinedComputeCCMom.computeMoments(a_orderPmax,a_degreePmax,a_useConstraints,a_refinementCriterion,numberOfRefinements);
669 #endif
670  refinedCutCellMoments.push_back(refinedComputeCCMom.m_cutCellMoments);
671 
672  // Move to the next subcell (if there is one)
673  int idir;
674  for (idir = 0; idir < dim; idir++)
675  {
676  if (curIter[idir] < maxIter[idir])
677  {
678  // If the current index is in range, increment it and continue
679  curIter[idir]++;
680  break;
681  }
682  else
683  {
684  // If the current index at the maximum, reset it and move to the next
685  // index
686  curIter[idir] = 0;
687  }
688  }
689 
690  // All the subcells have been done
691  if (idir == dim)
692  {
693  break;
694  }
695  }
696 
697  return refinedCutCellMoments;
698 }
699 
700 template <int dim> void ComputeCutCellMoments<dim>::addMomentMaps(const Vector<CutCellMoments<dim> > & a_refinedCutCellVector,
702  const int & a_degreeP,
703 #else
704  const int & a_degreePmax,
705 #endif
706  const bool & a_useConstraints)
707 {
708  int numberOfRefinedCC = a_refinedCutCellVector.size();
709  CH_assert(numberOfRefinedCC>0);
710 
711  IndexTM<int,2> bdId;
712 
713  // Pick one refined CutCellMoments and use its maps to initialize the maps
714  // of the current CutCellMoment to zero the chosen cutcellmoment need to
715  // have initialized boundary moments, i.e., be irregular
716  CutCellMoments<dim> cutCell;
717  int k = 0;
718 
719  while (k < numberOfRefinedCC && (a_refinedCutCellVector[k].m_IFData.m_allVerticesOut ||
720  (a_refinedCutCellVector[k].m_IFData.m_allVerticesIn && !m_cutCellMoments.m_bdCCOn)))
721  {
722  k++;
723  }
724 
725  if (k >= numberOfRefinedCC)
726  {
727  MayDay::Abort("Exceeded the number of refined ComputeCutCellMoments in vector");
728  }
729 
730  else
731  {
732  cutCell = a_refinedCutCellVector[k];
733  }
734 
735  m_cutCellMoments.initialize(cutCell);
736 
737  // Add up the moments on the refined cells to obtain the moment in the
738  // coarse cell. Move moments to the m_parentCoord, then add.
739  for (int i = 0; i < numberOfRefinedCC; i++)
740  {
741  // The difference between the refined cell center and the current cell
742  // center is needed to change the coordinate system of the refined moments
743  // to the coordinate system where they are needed.
744  CutCellMoments<dim> refinedCutCell = a_refinedCutCellVector[i];
746  IndexTM<Real,dim> refinedCenterDelta = IndexTM<Real,dim>::Zero;
747 
748  int ii = i;
749  for (int j = 0; j < dim; ++j)
750  {
751  hilo[j] = (ii & 1);
752  ii = ii >> 1;
753  }
754 
755  IndexTM<int,dim> sign = 2*hilo-1;
756 
757  for (int idir = 0; idir < dim; idir++)
758  {
759  refinedCenterDelta[idir] = 0.25*sign[idir]*m_cutCellMoments.m_IFData.m_globalCoord.m_dx[idir];
760  }
761 
762  // Add up the moments in the volume
763  addMoments(m_cutCellMoments.m_moments,refinedCutCell.m_moments,refinedCenterDelta);
764  addMoments(m_cutCellMoments.m_EBmoments,refinedCutCell.m_EBmoments,refinedCenterDelta);
765 
766  // Add up the moments on the boundary
767  for (int idir = 0; idir < dim; idir++)
768  {
769  bdId[0] = idir;
770 
771  // The difference in cell centers in dim-1 is needed to change the
772  // coordinate system of the refined moments to the coordinate
773  // system where they are needed.
774  IndexTM<Real,dim-1> localRefinedCenterDelta;
775  IndexTM<int,dim-1> localHilo;
776 
777  for (int jdir = 0; jdir < dim; jdir++)
778  {
779  if (jdir < idir)
780  {
781  localRefinedCenterDelta[jdir] = refinedCenterDelta[jdir];
782  localHilo[jdir] = hilo[jdir];
783  }
784  else if (jdir > idir)
785  {
786  localRefinedCenterDelta[jdir-1] = refinedCenterDelta[jdir];
787  localHilo[jdir-1] = hilo[jdir];
788  }
789  }
790 
791  if (hilo[idir] != LARGEINTVAL)
792  {
793  bdId[1] = hilo[idir];
794 
795  // Drop one dimension to add up the boundary moments dropping
796  // one dimension allows an easier implementation as 1D is a
797  // particular case
798  (refinedCutCell.m_bdCutCellMoments[bdId]).addBdMoments(
799  m_cutCellMoments.m_bdCutCellMoments[bdId],
800  refinedCutCell.m_IFData,
802  a_degreeP,
803 #else
804  a_degreePmax,
805 #endif
806  a_useConstraints,
807  localRefinedCenterDelta,
808  localHilo);
809  }
810  else
811  {
812  for (int side = 0; side < 2; side++)
813  {
814  bdId[1] = side;
815  (refinedCutCell.m_bdCutCellMoments[bdId]).addBdMoments(
816  m_cutCellMoments.m_bdCutCellMoments[bdId],
817  refinedCutCell.m_IFData,
819  a_degreeP+1,
820 #else
821  a_degreePmax+1,
822 #endif
823  a_useConstraints,
824  localRefinedCenterDelta,
825  localHilo);
826  }
827  }
828  }
829  }
830 }
831 
832 template <int dim> void ComputeCutCellMoments<dim>::addMoments(PthMoment & a_momentMap,
833  PthMoment & a_refinedMomentMap,
834  const IndexTM<Real,dim> & a_refinedCenterDelta)
835 {
836  // Iterate through the REFINEDmap and add the refined moment in the right
837  // coordinate system to the coarse map
838  for (typename PthMoment::const_iterator it = a_refinedMomentMap.begin(); it != a_refinedMomentMap.end(); ++it)
839  {
840  IndexTM<int,dim> monomial = it->first;
841  Real moment = m_cutCellMoments.changeMomentCoordinates(a_refinedMomentMap,monomial,a_refinedCenterDelta);
842 
843  a_momentMap[it->first] += moment;
844  }
845 }
846 
847 #if RECURSIVE_GEOMETRY_GENERATION == 0
848 template <int dim> void ComputeCutCellMoments<dim>::computeResiduals(const int & a_order,
849  const int & a_degreeP,
850 #else
851 template <int dim> void ComputeCutCellMoments<dim>::computeResiduals(const int & a_orderPmax,
852  const int & a_degreePmax,
853 #endif
854  const bool & a_useConstraints)
855 {
856 #if RECURSIVE_GEOMETRY_GENERATION == 0
857  for (int iDegree = a_degreeP; iDegree >= 0; iDegree--)
858 #else
859  for (int iDegree = a_degreePmax; iDegree >= 0; iDegree--)
860 #endif
861  {
862  // Initialize residuals to 0
863  int nNorm = 3;
864  m_cutCellMoments.m_residual[iDegree].resize(nNorm);
865 
866  for (int i = 0; i < nNorm; i++)
867  {
868  m_cutCellMoments.setResidual(0.0,iDegree,i);
869  }
870 
871  // Create a lsp problem
872  IvDim zeroDerivative = IndexTM<int,dim>::Zero;
873 #if RECURSIVE_GEOMETRY_GENERATION == 0
874  LSProblem<dim> lsp(a_order,iDegree,a_useConstraints,m_cutCellMoments.m_IFData.m_normalDerivatives[zeroDerivative]);
875 #else
876  LSProblem<dim> lsp(a_orderPmax,iDegree,a_useConstraints,m_cutCellMoments.m_IFData.m_normalDerivatives[zeroDerivative]);
877 #endif
878 
879 
880  // Compute the right hand side
881 #if RECURSIVE_GEOMETRY_GENERATION == 0
882  Vector<Real> rhs = computeRhs(lsp,a_order);
883 #else
884  Vector<Real> rhs = computeRhs(lsp,a_orderPmax);
885 #endif
886 
887  // Get the unknowns from the moment maps
888  Vector<Real> unknowns(lsp.m_numP+lsp.m_numPLess1);
889 
890  for (typename PthMomentLoc::const_iterator it = lsp.m_monoLocP.begin(); it != lsp.m_monoLocP.end(); ++it)
891  {
892  unknowns[it->second] = m_cutCellMoments.m_EBmoments[it->first];
893  }
894 
895  for (typename PthMomentLoc::const_iterator it = lsp.m_monoLocPLess1.begin(); it != lsp.m_monoLocPLess1.end(); ++it)
896  {
897  unknowns[it->second + lsp.m_numP] = m_cutCellMoments.m_moments[it->first];
898  }
899 
900  // Compute the residuals
901  Real maxRi = 0.0;
902 
903  for (int i = 0; i < lsp.m_numP*dim; i++)
904  {
905  Real AtimeX = 0.0;
906 
907  for (int j = 0; j < lsp.m_numP + lsp.m_numPLess1; j++)
908  {
909  AtimeX += lsp.m_matrix[i][j] * unknowns[j];
910  }
911 
912  Real ri = AtimeX - rhs[i];
913 
914  if (Abs(ri) > maxRi)
915  {
916  m_cutCellMoments.setResidual(Abs(ri),iDegree,0);
917  maxRi = Abs(ri);
918  }
919 
920  m_cutCellMoments.setResidual(m_cutCellMoments.getResidual(iDegree,1) + Abs(ri),iDegree,1);
921  m_cutCellMoments.setResidual(m_cutCellMoments.getResidual(iDegree,2) + ri * ri,iDegree,2);
922  }
923 
924  m_cutCellMoments.setResidual(sqrt(m_cutCellMoments.getResidual(iDegree,2)),iDegree,2);
925  }
926 }
927 
928 template <int dim> void ComputeCutCellMoments<dim>::computeResiduals(const Vector< CutCellMoments<dim> > & a_refinedCCMoms,
930  const int & a_degreeP)
931 #else
932  const int & a_degreePmax)
933 #endif
934 {
935 #if RECURSIVE_GEOMETRY_GENERATION == 0
936  for (int iDegree = 0; iDegree <= a_degreeP; iDegree++)
937 #else
938  for (int iDegree = 0; iDegree <= a_degreePmax; iDegree++)
939 #endif
940  {
941  Real maxRes = 0.0;
942  Real tempRes2 = 0.0;
943  Real tempRes1 = 0.0;
944 
945  for (int i = 0; i < a_refinedCCMoms.size(); i++)
946  {
947  Real res0 = a_refinedCCMoms[i].getResidual(iDegree,0);
948  Real res1 = a_refinedCCMoms[i].getResidual(iDegree,1);
949  Real res2 = a_refinedCCMoms[i].getResidual(iDegree,2);
950 
951  if (res0 > maxRes && res0 != LARGEREALVAL)
952  {
953  maxRes = res0;
954  m_cutCellMoments.setResidual(res0,iDegree,0);
955  }
956 
957  if (res1 != LARGEREALVAL)
958  {
959  tempRes1 += res1;
960  }
961 
962  if (res2 != LARGEREALVAL)
963  {
964  tempRes2 += res2*res2;
965  }
966  }
967 
968  m_cutCellMoments.setResidual(tempRes1,iDegree,1);
969 
970  if (tempRes2 != LARGEREALVAL)
971  {
972  m_cutCellMoments.setResidual(sqrt(tempRes2),iDegree,2);
973  }
974  }
975 }
976 
977 template <int dim> void ComputeCutCellMoments<dim>::print(ostream& a_out) const
978 {
979  m_cutCellMoments.print(a_out);
980 }
981 
982 template <int dim> void ComputeCutCellMoments<dim>::dump() const
983 {
984  print(pout());
985 }
986 
987 // Operators
988 template <int dim> void ComputeCutCellMoments<dim>::operator=(const ComputeCutCellMoments<dim> & a_computeCutCellMoments)
989 {
990  // Only copy if the objects are distinct
991  if (this != &a_computeCutCellMoments)
992  {
993  m_cutCellMoments = a_computeCutCellMoments.m_cutCellMoments;
994  }
995 }
996 
997 template <int dim> ostream& operator<<(ostream & a_out,
998  const ComputeCutCellMoments<dim> & a_computeCutCellMoments)
999 {
1000  a_computeCutCellMoments.print(a_out);
1001  return a_out;
1002 }
1003 
1004 template <int dim> Real ComputeCutCellMoments<dim>::factorial(const IvDim & a_multiIndex) const
1005 {
1006  Real fact = 1.0;
1007 
1008  for (int i = 0; i < dim; i++)
1009  {
1010  for (int j = 2; j <= a_multiIndex[i]; j++)
1011  {
1012  fact *= j;
1013  }
1014  }
1015 
1016  return fact;
1017 }
1018 
1019 #include "NamespaceFooter.H"
1020 
1021 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
int numActiveBounds() const
Definition: LSProblem.H:119
void computeBounds(const IndexTM< Real, dim > &a_dx, const CutCellMoments< dim > &a_ccm)
Definition: LSProblemImplem.H:94
#define BDID_HILO
Definition: Notation.H:89
BdCutCellMoments m_bdCutCellMoments
Definition: CutCellMoments.H:138
ostream & operator<<(ostream &a_out, const ComputeCutCellMoments< dim > &a_computeCutCellMoments)
Definition: ComputeCutCellMomentsImplem.H:997
void addMoments(PthMoment &a_momentMap, PthMoment &a_refinedMomentMap, const IndexTM< Real, dim > &a_refinedCenterDelta)
Definition: ComputeCutCellMomentsImplem.H:832
int m_numP
Definition: LSProblem.H:186
#define CH_assert(cond)
Definition: CHArray.H:37
Box refine(const Box &b, const IntVect &refinement_ratio)
Definition: NoRefinement.H:22
map< int, IvDim > LocPthMoment
Definition: ComputeCutCellMoments.H:48
#define LARGEREALVAL
Definition: Notation.H:77
void operator=(const ComputeCutCellMoments< dim > &a_computeCutCellMoments)
Definition: ComputeCutCellMomentsImplem.H:988
Definition: ComputeCutCellMoments.H:37
#define LARGEINTVAL
Definition: Notation.H:76
Vector< CutCellMoments< dim > > refine(const int &a_order, const int &a_degreeP, const bool &a_useConstraints, RefinementCriterion< dim > &a_refinementCriterion, const IndexTM< int, dim > &a_refineInDir, const int &a_numberOfRefinements)
Definition: ComputeCutCellMomentsImplem.H:592
virtual bool doRefine(IndexTM< int, dim > &a_refineInDir, const CutCellMoments< dim > &a_ccm, const int &a_numberOfRefinements)
Definition: RefinementCriterion.H:101
map< IvDim, Real, LexLT< IvDim > > PthMoment
Definition: ComputeCutCellMoments.H:43
Vector< Real > computeRhs(LSProblem< dim > &a_lsp, const int &a_order)
Definition: ComputeCutCellMomentsImplem.H:485
~ComputeCutCellMoments()
Definition: ComputeCutCellMomentsImplem.H:82
Real getUnknown(int loc)
Definition: LSProblem.H:104
void setConstrantSuccessStatus(const bool &a_status)
Definition: RefinementCriterion.H:110
CutCellMoments< dim > m_cutCellMoments
Definition: ComputeCutCellMoments.H:149
Definition: ComputeCutCellMoments.H:35
Definition: RefinementCriterion.H:27
void dump() const
Definition: ComputeCutCellMomentsImplem.H:982
PthMomentLoc m_monoLocPLess1
Definition: LSProblem.H:189
Definition: IndexTM.H:36
void computeResiduals(const int &a_order, const int &a_degreeP, const bool &a_useConstraints)
Definition: ComputeCutCellMomentsImplem.H:848
void print(ostream &a_out) const
Definition: LSProblemImplem.H:560
void push_back(const T &in)
Definition: Vector.H:295
IFData< dim > m_IFData
Definition: CutCellMoments.H:141
Definition: IFData.H:42
int sign(const Side::LoHiSide &a_side)
#define BDID_DIR
Definition: Notation.H:88
double Real
Definition: REAL.H:33
void addMomentMaps(const Vector< CutCellMoments< dim > > &a_refinedCutCellVector, const int &a_degreeP, const bool &a_useConstraints)
Definition: ComputeCutCellMomentsImplem.H:700
const PthMomentLoc & getMonomialLocMapDegreeP() const
Definition: LSProblem.H:94
int invertNormalEq(const Vector< Real > &a_rhs, Vector< Real > &a_residual)
Definition: LSProblemImplem.H:317
T Abs(const T &a_a)
Definition: Misc.H:53
PthMomentLoc m_monoLocP
Definition: LSProblem.H:182
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
void print(ostream &out) const
Definition: ComputeCutCellMomentsImplem.H:977
void initialize(CutCellMoments< dim > &a_refinedCutCell)
Definition: CutCellMomentsImplem.H:411
void computeMoments(const int &a_order, const int &a_degreeP, const bool &a_useConstraints, RefinementCriterion< dim > &a_refinementCriterion, const int &a_numberOfRefinements=0)
Definition: ComputeCutCellMomentsImplem.H:87
int m_numPLess1
Definition: LSProblem.H:193
PthMoment m_moments
Definition: CutCellMoments.H:132
Real ** m_matrix
Definition: LSProblem.H:197
Real factorial(const IvDim &a_multiIndex) const
Definition: ComputeCutCellMomentsImplem.H:1004
int getNumberDegP()
Definition: LSProblem.H:109
Real factorial(const int n)
Calculates factorial for an integer.
Definition: Factorial.H:25
size_t size() const
Definition: Vector.H:192
const PthMomentLoc & getMonomialLocMapDegreePLess1() const
Definition: LSProblem.H:99
int m_maxOrder
Definition: IFData.H:67
PthMoment m_EBmoments
Definition: CutCellMoments.H:135
ComputeCutCellMoments()
Definition: ComputeCutCellMomentsImplem.H:36
void generateMultiIndices(Vector< IndexTM< int, dim > > &a_indices, const int &a_magnitude)
Definition: MultiIndexImplem.H:18
int dim
Definition: EBInterface.H:146
Definition: CutCellMoments.H:32
virtual bool baseDoRefine(IndexTM< int, dim > &a_refineInDir, const CutCellMoments< dim > &a_ccm, const int &a_numberOfRefinements)
Should a cell be subdivided and in which directions.
Definition: RefinementCriterion.H:63
const LocPthMoment & getLocMonomialMapDegreeP() const
Definition: LSProblem.H:89
#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).
map< IvDim, int, LexLT< IvDim > > PthMomentLoc
Definition: ComputeCutCellMoments.H:49