Chombo + EB + MF  3.2
DivNormalRefinementImplem.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 _DIVNORMALREFINEMENTIMPLEM_H_
12 #define _DIVNORMALREFINEMENTIMPLEM_H_
13 
14 #include "MayDay.H"
15 
16 #include "NamespaceHeader.H"
17 
19 {
20 }
21 
22 template <int dim> DivNormalRefinement<dim>::DivNormalRefinement(const Real & a_maxChangeThreshold,
23  const int& a_maxNumberOfRefinements)
24 {
25  setMaxChangeThreshold(a_maxChangeThreshold);
26  setMaxNumberOfRefinements(a_maxNumberOfRefinements);
27 }
28 
30 {
31 }
32 
33 template <int dim> bool DivNormalRefinement<dim>::doRefine(IndexTM<int,dim> & a_refineInDir,
34  const CutCellMoments<dim> & a_ccm,
35  const int & a_numberOfRefinements)
36 {
37  bool retval = false;
38  a_refineInDir = IndexTM<int,dim>::Zero;
39 
40  if (a_numberOfRefinements < m_maxNumberOfRefinements)
41  {
42  Real changeInNormal = approximateDivNormal(a_ccm);
43 
44  //check whether normal equals the zero vector or whether constraints were active
45  if (changeInNormal > m_maxChangeThreshold)
46  {
47  retval = true;
48  a_refineInDir = IndexTM<int,dim>::Unit;
49  }
50  }
51 
52  return retval;
53 }
54 
56 {
57  typedef IndexTM<int,dim> IvDim;
58  typedef IndexTM<Real,dim> RvDim;
59  typedef IndexTM<int,dim> EdgeIndex;
60  typedef map<EdgeIndex,Real,LexLT<EdgeIndex> > EdgeIntersections;
61 
62  Real changeInNormal = 0.0;
63  const IFData<dim>& edgeData = a_ccm.m_IFData;
64 
65  //normal at the center of the Taylor Expansion, which always exists.
66  const RvDim& normal = edgeData.m_normalDerivatives.find(IndexTM<int,dim>::Zero)->second;
67 
68  //iterate through intersection points
69  for (typename EdgeIntersections::const_iterator it = edgeData.m_intersections.begin();it != edgeData.m_intersections.end(); ++it)
70  {
71  const EdgeIndex& edgeIndex = it->first;
72  const Real& intercept = it->second;
73 
74  int varyDir = edgeIndex[0];
75 
76  RvDim intersectPt = RvDim::Zero;
77  intersectPt[varyDir] = intercept;
78 
79  for (int i = 1; i < dim; i++)
80  {
81  int curDir;
82  int loHi;
83 
84  if (i <= varyDir)
85  {
86  curDir = i-1;
87  }
88  else
89  {
90  curDir = i;
91  }
92 
93  loHi = edgeIndex[i];
94  intersectPt[curDir] = loHi - 0.5;
95  intersectPt[curDir] *= edgeData.m_globalCoord.m_dx[curDir];
96  }
97 
98  //normal at intersection point is the zeroth derivative of the normal
99  NormalDerivative<dim> normalDerivative;
100  RvDim intersectPtNormal;
101  for (int idir = 0; idir < dim; idir++)
102  {
103  intersectPtNormal[idir] = normalDerivative.evaluate(IndexTM<int,dim>::Zero,
104  idir,
105  edgeData.m_globalCoord.convert(intersectPt,edgeData.m_cellCenterCoord),
106  edgeData.m_function);
107  }
108 
109  // find the angle between the normal at intersection point the normal at the average of interesection points
110  Real dotProduct = 0.0;
111  for (int idir = 0; idir < dim; idir++)
112  {
113  dotProduct += intersectPtNormal[idir]*normal[idir];
114  }
115 
116  if (Abs(dotProduct - 1.0) > changeInNormal)
117  {
118  changeInNormal = Abs(dotProduct -1.0);
119  }
120  }
121  return changeInNormal;
122 }
123 
124 template <int dim> void DivNormalRefinement<dim>::setMaxChangeThreshold(const Real & a_maxChangeThreshold)
125 {
126  if (a_maxChangeThreshold < 0)
127  {
128  MayDay::Abort("DivNormalRefinement<dim>::setMaxChangeThreshold - maxChangeThreshold must be >= 0");
129  }
130 
131  m_maxChangeThreshold = a_maxChangeThreshold;
132 }
133 
135 {
136  return m_maxChangeThreshold;
137 }
138 
139 template <int dim> void DivNormalRefinement<dim>::setMaxNumberOfRefinements(const int & a_maxNumberOfRefinements)
140 {
141  if (a_maxNumberOfRefinements < 0)
142  {
143  MayDay::Abort("DivNormalRefinement<dim>::setMaxNumberOfRefinements - maxNumberOfRefinents must be >= 0");
144  }
145 
146  m_maxNumberOfRefinements = a_maxNumberOfRefinements;
147 }
148 
150 {
151  return m_maxNumberOfRefinements;
152 }
153 
154 #include "NamespaceFooter.H"
155 
156 #endif
NormalDerivatives m_normalDerivatives
Definition: IFData.H:68
virtual bool doRefine(IndexTM< int, dim > &a_refineInDir, const CutCellMoments< dim > &a_ccm, const int &a_numberOfRefinements)
Definition: DivNormalRefinementImplem.H:33
virtual ~DivNormalRefinement()
Destructor.
Definition: DivNormalRefinementImplem.H:29
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
This computes the derivatives of the normal of a sliced implicit function.
Definition: NormalDerivative.H:25
IFSlicer< dim > * m_function
Definition: IFData.H:60
Definition: IndexTM.H:36
IFData< dim > m_IFData
Definition: CutCellMoments.H:141
Definition: IFData.H:42
CoordinateSystem< dim > m_globalCoord
Definition: IFData.H:62
virtual Real approximateDivNormal(const CutCellMoments< dim > &a_ccm)
Definition: DivNormalRefinementImplem.H:55
double Real
Definition: REAL.H:33
T Abs(const T &a_a)
Definition: Misc.H:53
virtual void setMaxNumberOfRefinements(const int &a_maxNumberOfRefinements)
Definition: DivNormalRefinementImplem.H:139
virtual Real getMaxChangeThreshold()
Get threshold.
Definition: DivNormalRefinementImplem.H:134
CoordinateSystem< dim > m_cellCenterCoord
Definition: IFData.H:63
DivNormalRefinement()
Null constructor.
Definition: DivNormalRefinementImplem.H:18
virtual void setMaxChangeThreshold(const Real &a_maxChangeThreshold)
Definition: DivNormalRefinementImplem.H:124
EdgeIntersections m_intersections
Definition: IFData.H:59
virtual int getMaxNumberOfRefinements()
Get threshold.
Definition: DivNormalRefinementImplem.H:149
int dim
Definition: EBInterface.H:146
Definition: CutCellMoments.H:32
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).