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