Chombo + EB  3.0
IFDataImplem.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 _IFDATAIMPLEM_H_
12 #define _IFDATAIMPLEM_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 <iostream>
21 #include <iomanip>
22 
23 #include "MultiIndex.H"
24 #include "NormalDerivative.H"
25 #include "ParmParse.H"
26 
27 #include "NamespaceHeader.H"
28 
29 // empty constructor
30 template <int dim> IFData<dim>::IFData()
31 {
32  m_function = NULL;
33 }
34 
35 
36 template <int dim> IFData<dim>::IFData(const IFData<dim>& a_IFData)
37  :m_cornerSigns (a_IFData.m_cornerSigns),
38  m_intersections (a_IFData.m_intersections),
39  m_globalCoord (a_IFData.m_globalCoord),
40  m_cellCenterCoord (a_IFData.m_cellCenterCoord),
41  m_parentCoord (a_IFData.m_parentCoord),
42  m_localCoord (a_IFData.m_localCoord),
43  m_maxOrder (a_IFData.m_maxOrder),
44  m_normalDerivatives (a_IFData.m_normalDerivatives),
45  m_badNormal (a_IFData.m_badNormal),
46  m_allVerticesIn (a_IFData.m_allVerticesIn),
47  m_allVerticesOut (a_IFData.m_allVerticesOut),
48  m_allVerticesOn (a_IFData.m_allVerticesOn)
49 {
50  if (a_IFData.m_function != NULL)
51  {
52  // IFData always owns the data pointed to by m_function
53  m_function = new IFSlicer<dim>(*(a_IFData.m_function));
54  }
55  else
56  {
57  m_function = NULL;
58  }
59 }
60 
61 // Constructor from the implicit function at the highest dimension
62 template <int dim> IFData<dim>::IFData(const BaseIF & a_function,
63  const RvDim & a_dx,
64  const RvDim & a_cellCenter,
65  const int & a_maxOrder)
66  :m_globalCoord (a_cellCenter,a_dx),
67  m_cellCenterCoord(IndexTM<Real,dim>::Zero,a_dx),
68  m_parentCoord (IndexTM<Real,dim>::Zero,a_dx),
69  m_maxOrder (a_maxOrder)
70 {
71  // Copy the implicit function; m_function is a pointer
72  // IFData always owns the data pointed to by m_function
73  m_function = new IFSlicer<dim>(a_function);
74 
75  // makecornerSigns sets the bools m_allVerticesIn and m_allVerticesOut,
76  // and m_allVerticesOn
78 
79  // find an intersection point, if it exists, of the interface with each edge
81 
82  // will ultimately assign m_localCoord.m_origin to the the average of
83  // intersection points. If full or empty cell, m_localCoord.m_origin = m_parentCoord.m_origin
85 
86  // set the normal derivatives
88 }
89 
90 // Constructor using IFSlicer, used for refinement
91 template <int dim> IFData<dim>::IFData(IFSlicer<dim> * a_function,
92  const RvDim & a_dx,
93  const RvDim & a_cellCenter,
94  const int & a_maxOrder)
95  :m_globalCoord (a_cellCenter,a_dx),
96  m_cellCenterCoord(IndexTM<Real,dim>::Zero,a_dx),
97  m_parentCoord (IndexTM<Real,dim>::Zero,a_dx),
98  m_maxOrder (a_maxOrder)
99 {
100  // IFData always owns the data pointed to by m_function
101  m_function = new IFSlicer<dim>(*a_function);
102 
103  // makeCornerSigns sets the bools m_allVerticesIn and m_allVerticesOut,
104  // and m_allVerticesOn
105  makeCornerSigns();
106 
107  // Find an intersection point, if it exists, of the interface with each edge
109 
110  //
112 
113  // Set the normal derivatives
115 }
116 
117 template <int dim> IFData<dim>::IFData(const IFData<dim+1> & a_hIFData,
119  const int & a_maxOrder,
120 #endif
121  const int & a_idir,
122  const int & a_hilo)
123 {
124  typedef IndexTM<int,dim+1>HDEdgeIndex;
125  typedef map<HDEdgeIndex,Real,LexLT<HDEdgeIndex> > HDEdgeIntersections;
126  typedef IndexTM<int,dim+1>HDVertex;
127  typedef map<HDVertex,int,LexLT<HDVertex> > HDCornerSigns;
128 
129  // pm = -1 for lo and 1 for hi
130  int pm = a_hilo;
131  pm *= 2;
132  pm -= 1;
133 
134  // fixedComp and fixed value are used by the slicer to pull points into one higher dimension
135  int fixedComp = a_idir;
136  Real fixedValue = a_hIFData.m_globalCoord.convertDir(pm*0.5*a_hIFData.m_cellCenterCoord.m_dx[fixedComp],
137  a_hIFData.m_cellCenterCoord,
138  fixedComp);
139 
140  // IFData always owns the data pointed to by m_function
141  m_function = new IFSlicer<dim>(a_hIFData.m_function,fixedComp,fixedValue);
142 
143  // reduce the global coordinate system
144  m_globalCoord = CoordinateSystem<dim>(a_hIFData.m_globalCoord,fixedComp);
145 
146  // Construct cell centered coord sys
147  m_cellCenterCoord = CoordinateSystem<dim>(a_hIFData.m_cellCenterCoord,fixedComp);
148 
149  // reduce the parent coordinate system
150  m_parentCoord = CoordinateSystem<dim>(a_hIFData.m_localCoord,fixedComp);
151 
152  m_allVerticesIn = true;
153  m_allVerticesOut = true;
154  m_allVerticesOn = true;
155 
156  for (typename HDCornerSigns::const_iterator it = a_hIFData.m_cornerSigns.begin();
157  it != a_hIFData.m_cornerSigns.end(); ++it)
158  {
159  const HDVertex& hVertex = it->first;
160  if (hVertex[fixedComp] == a_hilo)
161  {
162  Vertex vertex;
163  for (int j = 0; j < dim; ++j)
164  {
165  if (j < fixedComp)
166  {
167  vertex[j] = hVertex[j];
168  }
169  else
170  {
171  vertex[j] = hVertex[j+1];
172  }
173  }
174 
175  m_cornerSigns[vertex] = it->second;
176  if (m_cornerSigns[vertex] == IN)
177  {
178  m_allVerticesOut = false;
179  m_allVerticesOn = false;
180  }
181  else if (m_cornerSigns[vertex] == OUT)
182  {
183  m_allVerticesIn = false;
184  m_allVerticesOn = false;
185  }
186  }
187  }
188 
189  for (typename HDEdgeIntersections::const_iterator it = a_hIFData.m_intersections.begin();
190  it != a_hIFData.m_intersections.end(); ++it)
191  {
192  const HDEdgeIndex& hEdgeIndex = it->first;
193  EdgeIndex edgeIndex;
194  int hEdgeDir = hEdgeIndex[0];
195 
196  // only interested in edgeDir != fixedComp(ie. a_idir)
197  // Among these, need to find edges in the higher dim that are a_hilo in the fixedComp component.
198  // Example: fixedComp = 1; a_hilo = 0. Observe that an edge in the x (N.B 0 < 1) direction
199  // might look like (edgeDir= 0:yhilo,zhilo ). In other words, y information is in index 1.
200  // On the other hand, an edge in the z direction(N.B 2 > 1)
201  // might look like (edgeDir= 0:xhilo,yhilo ). In other words, y information is in index 2.
202  if (hEdgeDir < fixedComp)
203  {
204  if (hEdgeIndex[fixedComp] == a_hilo)
205  {
206  edgeIndex[0] = hEdgeDir;
207 
208  for (int j = 1; j < dim; ++j)
209  {
210  if (j < fixedComp)
211  {
212  edgeIndex[j] = hEdgeIndex[j];
213  }
214  else
215  {
216  edgeIndex[j] = hEdgeIndex[j+1];
217  }
218  }
219  m_intersections[edgeIndex] = it->second;
220  }
221  }
222  else if (hEdgeDir > fixedComp)
223  {
224  if (hEdgeIndex[fixedComp+1] == a_hilo)
225  {
226  edgeIndex[0] = hEdgeDir - 1;
227 
228  for (int j = 1; j < dim; ++j)
229  {
230  if (j < fixedComp+1)
231  {
232  edgeIndex[j] = hEdgeIndex[j];
233  }
234  else
235  {
236  edgeIndex[j] = hEdgeIndex[j+1];
237  }
238  }
239  m_intersections[edgeIndex] = it->second;
240  }
241  }
242  }
243 
244  // assigns m_localCoord.m_origin to the the average of intersection points
245  // or to m_parentCoord.m_origin if there are no interesection points
246  defineLocalCoords();
247 
248 #if RECURSIVE_GEOMETRY_GENERATION == 0
249  m_maxOrder = a_hIFData.m_maxOrder;
250 #else
251  m_maxOrder = a_maxOrder;
252 #endif
253 
254  if (!m_allVerticesOut && !m_allVerticesIn)
255  {
256  setNormalDerivatives();
257  }
258  else
259  {
260  m_badNormal = false;
261  }
262 }
263 
264 // Destructor
265 template <int dim> IFData<dim>::~IFData()
266 {
267  if (m_function != NULL)
268  {
269  delete m_function;
270  }
271 }
272 
273 template <int dim> void IFData<dim>::setNormalDerivatives()
274 {
275  NormalDerivative<dim> normalDerivative;
276 
277  for (int order = 0; order <= m_maxOrder; order++)
278  {
279  Vector<IvDim> derivatives;
280 
281  generateMultiIndices(derivatives,order);
282 
283  for (int i = 0; i < derivatives.size(); ++i)
284  {
285  const IvDim & derivative = derivatives[i];
286 
287  for (int idir = 0; idir < dim; idir++)
288  {
289  m_normalDerivatives[derivative][idir] = normalDerivative.evaluate(derivative,
290  idir,
291  m_globalCoord.convert(RvDim::Zero,m_localCoord),
292  m_function);
293  }
294  }
295  }
296 
297  // if dim < GLOBALDIM, then m_badNormal == true for dim => m_badNormal == true for GLOBALDIM
298  if (normalDerivative.getMagnitudeOfGradient()== 0.0)
299  {
300  m_badNormal = true;
301  }
302  else
303  {
304  m_badNormal = false;
305  }
306 }
307 
308 template <int dim> void IFData<dim>::makeCornerSigns()
309 {
310  m_allVerticesIn = true;
311  m_allVerticesOut = true;
312  m_allVerticesOn = true;
313 
314  Vertex vertex;
315 
316  // there are 2^dim vertices
317  int numVertices = 1 << dim;
318 
319  for (int i = 0; i < numVertices; ++i)
320  {
321  int ii = i;
322 
323  //label the ith vertex by i represented in base 2
324  for (int j = 0; j < dim; ++j)
325  {
326  //convert to base 2 by successivly "anding" with 1 and shifting
327  vertex[j] = (ii & 1);
328  ii = ii >> 1;
329  }
330 
331  // represent the vertex as an RvDim in cell centered coordinates
332  RvDim corner;
333  for (int idir = 0; idir < dim; ++idir)
334  {
335  corner[idir] = vertex[idir] - 0.5;
336  corner[idir] *= m_cellCenterCoord.m_dx[idir];
337  }
338 
339  // compute coordinates of a_vertex in global coordinates
340  RvDim cornerCoord = m_globalCoord.convert(corner,m_cellCenterCoord);
341 
342  // evaluate (zeroth derivative of) the function at the corner
343  Real val = m_function->value(IndexTM<int,dim>::Zero,cornerCoord);
344 
345  // true = negative = in the computational domain = in the fluid
346  if (val < -MACHINEPRECISION)
347  {
348  m_cornerSigns[vertex] = IN;
349  m_allVerticesOut = false;
350  m_allVerticesOn = false;
351  }
352  else if (val > MACHINEPRECISION)
353  {
354  m_cornerSigns[vertex] = OUT;
355  m_allVerticesIn = false;
356  m_allVerticesOn = false;
357  }
358  else
359  {
360  m_cornerSigns[vertex] = ON;
361  }
362  }
363 }
364 
365 template <int dim> void IFData<dim>::findIntersectionPts()
366 {
367  // double iterate through corner signs map,writing to a edge map if corner1 is connected to corner2
368  for (typename CornerSigns::const_iterator it0 = m_cornerSigns.begin();
369  it0 != m_cornerSigns.end(); ++it0)
370  {
371  for (typename CornerSigns::const_iterator it1 = m_cornerSigns.begin();
372  it1 != m_cornerSigns.end(); ++it1)
373  {
374  int edgeDir = LARGEINTVAL;
375  if (isConnected(edgeDir,it0->first,it1->first))
376  {
377  // make edge key:m_intersections[thisEdge] = pt;
378  makeEdgeKey(edgeDir,it0->first,it1->first);
379  }
380  }
381  }
382 }
383 template <int dim> void IFData<dim>::defineLocalCoords()
384 {
385  // default for full cells
386  if (m_intersections.size() == 0)
387  {
388  m_localCoord.m_origin = m_parentCoord.m_origin;
389  }
390  else
391  {
392  m_localCoord.m_origin = RvDim::Zero;
393 
394  for (typename EdgeIntersections::const_iterator it = m_intersections.begin();it != m_intersections.end(); ++it)
395  {
396  const EdgeIndex& edgeIndex = it->first;
397  const Real& intercept = it->second;
398 
399  int varyDir = edgeIndex[0];
400 
401  RvDim intersectPt = RvDim::Zero;
402  intersectPt[varyDir] = intercept;
403 
404  for (int i = 1; i < dim; i++)
405  {
406  int curDir;
407  int loHi;
408 
409  if (i <= varyDir)
410  {
411  curDir = i-1;
412  }
413  else
414  {
415  curDir = i;
416  }
417 
418  loHi = edgeIndex[i];
419  intersectPt[curDir] = loHi - 0.5;
420  intersectPt[curDir] *= m_globalCoord.m_dx[curDir];
421  }
422 
423  m_localCoord.m_origin -= intersectPt;
424  }
425 
426  m_localCoord.m_origin /= m_intersections.size();
427  }
428 
429  m_localCoord.m_dx = m_globalCoord.m_dx;
430 }
431 
432 template <int dim> bool IFData<dim>::isConnected(int & a_edgeDir,
433  const Vertex & a_vertex1,
434  const Vertex & a_vertex2)
435 {
436  bool connected = true;
437  int numDif = 0;
438 
439  // connected = true if and only if (a_vertex1 and a_vertex2 differ in
440  // exactly one coordinate.)
441  for (int idir = 0; idir < dim; ++idir)
442  {
443  if (a_vertex1[idir] != a_vertex2[idir])
444  {
445  // vertices differ in idir direction
446  a_edgeDir = idir;
447  numDif += 1;
448  }
449  }
450 
451  if (numDif == 1)
452  {
453  connected = true;
454  }
455  else
456  {
457  connected = false;
458  a_edgeDir = LARGEINTVAL;
459  }
460 
461  return connected;
462 }
463 
464 template <int dim> void IFData<dim>::makeEdgeKey(const int & a_edgeDir,
465  const Vertex & a_vertex0,
466  const Vertex & a_vertex1)
467 {
468  EdgeIndex thisEdge;
469 
470  thisEdge[0] = a_edgeDir;
471 
472  for (int idir = 0; idir < dim; ++idir)
473  {
474  if (idir < a_edgeDir)
475  {
476  thisEdge[idir + 1] = a_vertex0[idir];
477  }
478  else if (idir > a_edgeDir)
479  {
480  thisEdge[idir] = a_vertex0[idir];
481  }
482  }
483 
484  int lo = LARGEINTVAL;
485  int hi = LARGEINTVAL;
486 
487  if (m_cornerSigns.find(a_vertex0)!= m_cornerSigns.end())
488  {
489  lo = m_cornerSigns.find(a_vertex0)->second;
490  }
491  else
492  {
493  MayDay::Abort("Vertex not well defined in makeEdgeKey");
494  }
495 
496  if (m_cornerSigns.find(a_vertex1)!= m_cornerSigns.end())
497  {
498  hi = m_cornerSigns.find(a_vertex1)->second;
499  }
500  else
501  {
502  MayDay::Abort("Vertex not well defined in makeEdgeKey");
503  }
504 
505  if ((lo == OUT && hi == IN) || (lo == IN && hi == OUT))
506  {
507  // calculate a value between -1.0 and 1.0
508  Real pt = rootFinder(thisEdge);
509 
510  // check whether intersection is at the end point
511  bool hiOn = false;
512  bool loOn = false;
513 
514  checkIntersection(hiOn,loOn,pt);
515 
516  if (hiOn || loOn)
517  {
518  if (loOn)
519  {
520  m_cornerSigns[a_vertex0] = ON;
521  lo = ON;
522  }
523 
524  if (hiOn)
525  {
526  m_cornerSigns[a_vertex1] = ON;
527  hi = ON;
528  }
529 
530  remakeCornerSigns();
531  }
532  else
533  {
534  // Scale the intersection between -dx/2 and dx/2
535  m_intersections[thisEdge] = 0.5*m_cellCenterCoord.m_dx[thisEdge[0]]*pt;
536  }
537  }
538 
539  // If the edge in full and the EB intersects a vertex, record an
540  // intersection point
541  if (lo == IN && hi == ON)
542  {
543  // Intersection at -dx/2
544  m_intersections[thisEdge] = -0.5*m_cellCenterCoord.m_dx[thisEdge[0]];
545  }
546  else if (lo == ON && hi == IN)
547  {
548  // Intersection at dx/2
549  m_intersections[thisEdge] = 0.5*m_cellCenterCoord.m_dx[thisEdge[0]];
550  }
551 }
552 
553 template <int dim> Real IFData<dim>::rootFinder(const EdgeIndex& a_thisEdge)
554 {
555  Real pt = LARGEREALVAL;
556 
557  // the first component of an EdgeIndex gives the direction of the edge
558  // the other components select an edge by using the same (lo-hi) logic as Vertex
559  int edgeDir = a_thisEdge[0];
560 
561  Vertex loEnd;
562  loEnd[edgeDir] = 0;
563 
564  Vertex hiEnd;
565  hiEnd[edgeDir] = 1;
566 
567  for (int idir = 0; idir < dim; ++idir)
568  {
569  if (idir < edgeDir)
570  {
571  loEnd[idir] = a_thisEdge[idir + 1];
572  hiEnd[idir] = a_thisEdge[idir + 1];
573  }
574  else if (idir > edgeDir)
575  {
576  loEnd[idir] = a_thisEdge[idir];
577  hiEnd[idir] = a_thisEdge[idir];
578  }
579  }
580 
581  // computes coordinates of a_vertex relative to the cell-center = (0,0...0)
582  // represent the vertex as an RvDim in cell centered coordinates
583  RvDim loCorner;
584  RvDim hiCorner;
585  for (int idir = 0; idir < dim; ++idir)
586  {
587  loCorner[idir] = (loEnd[idir] - 0.5) * m_cellCenterCoord.m_dx[idir];
588  hiCorner[idir] = (hiEnd[idir] - 0.5) * m_cellCenterCoord.m_dx[idir];
589  }
590 
591  // compute coordinates of a_vertex in global coordinates
592  RvDim loPt = m_globalCoord.convert(loCorner,m_cellCenterCoord);
593  RvDim hiPt = m_globalCoord.convert(hiCorner,m_cellCenterCoord);
594 
595  // range of possible roots
596  // Real smallestRoot = -0.5;
597  // Real biggestRoot = 0.5;
598 
599  // returns a number between -1.0 and 1.0
600  pt = BrentRootFinder(loPt,hiPt,edgeDir); // ,smallestRoot,biggestRoot);
601 
602  // return this result
603  return pt;
604 }
605 
606 template <int dim> Real IFData<dim>::BrentRootFinder(const RvDim & a_loPt,
607  const RvDim & a_hiPt,
608  const int & a_edgeDir) const
609 // const Real & a_smallestRoot,
610 // const Real & a_biggestRoot) const
611 {
612  // Max allowed iterations and floating point precision
613  const unsigned int MAXITER = 100;
614  const Real EPS = 3.0e-15;
615 
616  unsigned int i;
617  Real aPt;
618  Real bPt;
619  Real c, fa, fb, fc;
620  Real d, e;
621  Real tol1, xm;
622  Real p, q, r, s;
623 
624  aPt = -1.0; // a_smallestRoot;
625  bPt = 1.0; // a_biggestRoot;
626 
627  RvDim physCoordAPt = a_loPt;
628  RvDim physCoordBPt = a_hiPt;
629 
630  fa = -m_function->value(IndexTM<int,dim>::Zero,physCoordAPt);
631  fb = -m_function->value(IndexTM<int,dim>::Zero,physCoordBPt);
632 
633  // Init these to be safe
634  c = d = e = 0.0;
635 
636  if (fb*fa > 0)
637  {
638  pout() << "fa " << fa << " fb " << fb <<endl;
639  MayDay::Abort("IFData::BrentRootFinder. Root must be bracketed, but instead the supplied end points have the same sign.");
640  }
641 
642  fc = fb;
643 
644  for (i = 0; i < MAXITER; i++)
645  {
646  if (fb*fc > 0)
647  {
648  // Rename a, b, c and adjust bounding interval d
649  c = aPt;
650  fc = fa;
651  d = bPt - aPt;
652  e = d;
653  }
654 
655  if (Abs(fc) < Abs(fb))
656  {
657  aPt = bPt;
658  bPt = c;
659  c = aPt;
660  fa = fb;
661  fb = fc;
662  fc = fa;
663  }
664 
665  // Convergence check
666  tol1 = 2.0 * EPS * Abs(bPt) + 0.5 * TOLERANCE;
667  xm = 0.5 * (c - bPt);
668 
669  if (Abs(xm) <= tol1 || fb == 0.0)
670  {
671  break;
672  }
673 
674  if (Abs(e) >= tol1 && Abs(fa) > Abs(fb))
675  {
676  // Attempt inverse quadratic interpolation
677  s = fb / fa;
678  if (aPt == c)
679  {
680  p = 2.0 * xm * s;
681  q = 1.0 - s;
682  }
683  else
684  {
685  q = fa / fc;
686  r = fb / fc;
687  p = s * (2.0 * xm * q * (q-r) - (bPt-aPt) * (r-1.0));
688  q = (q-1.0) * (r-1.0) * (s-1.0);
689  }
690 
691  // Check whether in bounds
692  if (p > 0) q = -q;
693 
694  p = Abs(p);
695 
696  if (2.0 * p < Min(((Real)3.0)*xm*q-Abs(tol1*q), Abs(e*q)))
697  {
698  // Accept interpolation
699  e = d;
700  d = p / q;
701  }
702  else
703  {
704  // Interpolation failed, use bisection
705  d = xm;
706  e = d;
707  }
708  }
709  else
710  {
711  // Bounds decreasing too slowly, use bisection
712  d = xm;
713  e = d;
714  }
715 
716  // Move last best guess to a
717  aPt = bPt;
718  fa = fb;
719 
720  // Evaluate new trial root
721  if (Abs(d) > tol1)
722  {
723  bPt = bPt + d;
724  }
725  else
726  {
727  if (xm < 0) bPt = bPt - tol1;
728  else bPt = bPt + tol1;
729  }
730 
731  physCoordBPt[a_edgeDir] = ((1.0 - bPt)/2.0) * a_loPt[a_edgeDir] + ((1.0 + bPt)/2.0) * a_hiPt[a_edgeDir];
732  fb = -m_function->value(IndexTM<int,dim>::Zero,physCoordBPt);
733  }
734 
735  if (i >= MAXITER)
736  {
737  cerr << "BrentRootFinder: exceeding maximum iterations: "
738  << MAXITER
739  << "\n";
740  }
741 
742  return bPt;
743 }
744 
745 template <int dim> void IFData<dim>::checkIntersection(bool & a_hiOn,
746  bool & a_loOn,
747  const Real & a_pt) const
748 {
749  a_hiOn = false;
750  a_loOn = false;
751 
752  if (a_pt >= 1.0 - MACHINEPRECISION)
753  {
754  a_hiOn = true;
755  }
756  else if (a_pt <= -1.0 + MACHINEPRECISION)
757  {
758  a_loOn = true;
759  }
760 }
761 
762 template <int dim> void IFData<dim>::remakeCornerSigns()
763 {
764  m_allVerticesIn = true;
765  m_allVerticesOut = true;
766  m_allVerticesOn = true;
767 
768  for (typename CornerSigns::const_iterator it = m_cornerSigns.begin();
769  it != m_cornerSigns.end(); ++it)
770  {
771  if (it->second == IN)
772  {
773  m_allVerticesOut = false;
774  }
775 
776  if (it->second == OUT)
777  {
778  m_allVerticesIn = false;
779  }
780 
781  if (it->second != ON)
782  {
783  m_allVerticesOn = false;
784  }
785  }
786 }
787 
788 template <int dim> void IFData<dim>::print(ostream& a_out) const
789 {
790  string padding = " ";
791  for (int i = 0; i < GLOBALDIM - dim; i++)
792  {
793  padding += " ";
794  }
795 
796  if (m_function == NULL)
797  {
798  a_out << padding << "undefined IFData\n";
799  }
800  else
801  {
802  for (typename CornerSigns::const_iterator it = m_cornerSigns.begin();
803  it != m_cornerSigns.end(); ++it)
804  {
805  a_out << padding << "Vertex "
806  << "("
807  << it->first
808  << ") = "
809  << it->second
810  << "\n";
811  }
812  a_out << padding << "\n";
813 
814  for (typename EdgeIntersections::const_iterator it = m_intersections.begin();
815  it != m_intersections.end(); ++it)
816  {
817  //save output format
818  std::ios::fmtflags origFlags = a_out.flags();
819  int origWidth = a_out.width();
820  int origPrecision = a_out.precision();
821 
822  a_out << padding << "EdgeIndex "
823  << it->first
824  << " = "
825  << setprecision(16)
826  << setiosflags(ios::showpoint)
827  << setiosflags(ios::scientific)
828  << it->second
829  << "\n";
830 
831  //restore output format
832  a_out.flags(origFlags);
833  a_out.width(origWidth);
834  a_out.precision(origPrecision);
835  }
836  a_out << padding << "\n";
837 
838  a_out << padding << "m_allVerticesIn = " << m_allVerticesIn << "\n";
839  a_out << padding << "m_allVerticesOut = " << m_allVerticesOut << "\n";
840  a_out << padding << "m_allVerticesOn = " << m_allVerticesOn << "\n";
841  a_out << padding << "m_badNormal = " << m_badNormal << "\n";
842  a_out << padding << "\n";
843 
844  if (!m_allVerticesOut)
845  {
846  a_out << padding << "m_globalCoord = " << m_globalCoord ;
847  a_out << padding << "m_cellCenterCoord = " << m_cellCenterCoord;
848  a_out << padding << "m_parentCoord = " << m_parentCoord ;
849  a_out << padding << "m_localCoord = " << m_localCoord ;
850  a_out << padding << "\n";
851 
852  a_out << padding << "m_maxOrder = " << m_maxOrder << "\n";
853 
854  std::ios::fmtflags origFlags = a_out.flags();
855  int origWidth = a_out.width();
856  int origPrecision = a_out.precision();
857 
858  a_out << padding << "m_normalDerivatives = " << "\n";
859 
860  for (typename NormalDerivatives::const_iterator it = m_normalDerivatives.begin();
861  it != m_normalDerivatives.end();
862  ++it)
863  {
864  for (int idir = 0; idir < dim; idir++)
865  {
866  a_out << padding << " "
867  << it->first
868  << ", "
869  << idir
870  << ": "
871  << setprecision(16)
872  << setiosflags(ios::showpoint)
873  << setiosflags(ios::scientific)
874  << (it->second)[idir]
875  << "\n";
876  }
877  a_out << padding << "\n";
878  }
879 
880  a_out.flags(origFlags);
881  a_out.width(origWidth);
882  a_out.precision(origPrecision);
883  }
884  else
885  {
886  a_out << padding << "All vertices out" << "\n";
887  }
888  a_out << padding << "\n";
889  }
890 }
891 
892 template <int dim> ostream& operator<<(ostream & a_out,
893  const IFData<dim> & a_IFData)
894 {
895  a_IFData.print(a_out);
896  return a_out;
897 }
898 
899 // equals operator
900 template <int dim> void IFData<dim>::operator=(const IFData & a_IFData)
901 {
902  // Only do something if the objects are distinct
903  if (this != &a_IFData)
904  {
905  // Delete the existing m_function (if there is one)
906  if (m_function != NULL)
907  {
908  delete m_function;
909  }
910 
911  // IFData always owns the data pointed to by m_function
912  if (a_IFData.m_function != NULL)
913  {
914  m_function = new IFSlicer<dim>(*(a_IFData.m_function));
915  }
916  else
917  {
918  m_function = NULL;
919  }
920 
921  m_cornerSigns = a_IFData.m_cornerSigns;
922  m_intersections = a_IFData.m_intersections;
923 
924  m_globalCoord = a_IFData.m_globalCoord;
925  m_cellCenterCoord = a_IFData.m_cellCenterCoord;
926  m_parentCoord = a_IFData.m_parentCoord;
927  m_localCoord = a_IFData.m_localCoord;
928 
929  m_maxOrder = a_IFData.m_maxOrder;
930  m_normalDerivatives = a_IFData.m_normalDerivatives;
931  m_badNormal = a_IFData.m_badNormal;
932 
933  m_allVerticesIn = a_IFData.m_allVerticesIn;
934  m_allVerticesOut = a_IFData.m_allVerticesOut;
935  m_allVerticesOn = a_IFData.m_allVerticesOn;
936  }
937 }
938 
939 #include "NamespaceFooter.H"
940 
941 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
NormalDerivatives m_normalDerivatives
Definition: IFData.H:61
Definition: CoordinateSystem.H:34
void makeEdgeKey(const int &a_edgeDir, const Vertex &a_vertex1, const Vertex &a_vertex2)
Definition: IFDataImplem.H:464
~IFData()
Definition: IFDataImplem.H:265
#define LARGEREALVAL
Definition: Notation.H:77
#define LARGEINTVAL
Definition: Notation.H:76
one dimensional dynamic array
Definition: Vector.H:52
virtual Real evaluate(const IvDim &a_multiIndex, const int &a_direction, const RvDim &a_point, const IFSlicer< dim > *a_ifSlicer)
Evaluate derivatives of the normal of an IFSlicer class.
Definition: NormalDerivativeImplem.H:27
bool isConnected(int &a_edgeDir, const Vertex &a_vertex1, const Vertex &a_vertex2)
Definition: IFDataImplem.H:432
void print(ostream &out) const
Definition: IFDataImplem.H:788
CoordinateSystem< dim > m_parentCoord
Definition: IFData.H:57
This computes the derivatives of the normal of a sliced implicit function.
Definition: NormalDerivative.H:25
IFSlicer< dim > * m_function
Definition: IFData.H:53
#define TOLERANCE
Definition: Notation.H:78
IFData()
Definition: IFDataImplem.H:30
Real BrentRootFinder(const RvDim &a_loPt, const RvDim &a_hiPt, const int &a_edgeDir) const
Definition: IFDataImplem.H:606
CoordinateSystem< dim > m_localCoord
Definition: IFData.H:58
void operator=(const IFData &a_ifData)
Definition: IFDataImplem.H:900
Definition: BaseIF.H:30
void setNormalDerivatives()
Definition: IFDataImplem.H:273
Definition: IFData.H:35
CoordinateSystem< dim > m_globalCoord
Definition: IFData.H:55
bool m_allVerticesOn
Definition: IFData.H:66
double Real
Definition: REAL.H:33
void makeCornerSigns()
Definition: IFDataImplem.H:308
#define ON
Definition: Notation.H:84
T Abs(const T &a_a)
Definition: Misc.H:53
bool m_allVerticesIn
Definition: IFData.H:64
CoordinateSystem< dim > m_cellCenterCoord
Definition: IFData.H:56
Definition: IFSlicer.H:27
#define MACHINEPRECISION
Definition: Notation.H:79
Real getMagnitudeOfGradient()
Definition: NormalDerivativeImplem.H:64
bool m_allVerticesOut
Definition: IFData.H:65
T Min(const T &a_a, const T &a_b)
Definition: Misc.H:26
CornerSigns m_cornerSigns
Definition: IFData.H:51
size_t size() const
Definition: Vector.H:177
void remakeCornerSigns()
Definition: IFDataImplem.H:762
ostream & operator<<(ostream &a_out, const IFData< dim > &a_IFData)
Definition: IFDataImplem.H:892
void defineLocalCoords()
Definition: IFDataImplem.H:383
int m_maxOrder
Definition: IFData.H:60
EdgeIntersections m_intersections
Definition: IFData.H:52
#define IN
Definition: Notation.H:85
void findIntersectionPts()
Definition: IFDataImplem.H:365
#define OUT
Definition: Notation.H:83
void checkIntersection(bool &a_hiOn, bool &a_loOn, const Real &a_pt) const
Definition: IFDataImplem.H:745
Real rootFinder(const EdgeIndex &a_thisEdge)
Definition: IFDataImplem.H:553
void generateMultiIndices(Vector< IndexTM< int, dim > > &a_indices, const int &a_magnitude)
Definition: MultiIndexImplem.H:18
int dim
Definition: EBInterface.H:146
#define GLOBALDIM
Definition: Notation.H:35
bool m_badNormal
Definition: IFData.H:62
#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).