Chombo + EB  3.0
NormalDerivativeImplem.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 _NORMALDERIVATIVEIMPLEM_H_
12 #define _NORMALDERIVATIVEIMPLEM_H_
13 
14 #include "NamespaceHeader.H"
15 
16 // Null constructor
18 {
19 }
20 
21 // Destructor
23 {
24 }
25 
26 // Evaluate derivatives of the normal of an IFSlicer class
27 template <int dim> Real NormalDerivative<dim>::evaluate(const IvDim & a_multiIndex,
28  const int & a_direction,
29  const RvDim & a_point,
30  const IFSlicer<dim> * a_implicitFunction)
31 {
32  // This is the numerator of the "a_direction" component of the normal, i.e.
33  // the "a_direction" component of the gradient of the function.
34  DerivativeProduct gradientComponent;
35  gradientComponent[BASISV_TM<int,dim>(a_direction)] = 1;
36 
37  // Represent the "a_direction" component of the normal as the "a_direction"
38  // component of gradient of the function, "gradientComponent", divided by
39  // the magnitude of the gradient.
40  int magnitudeOfGradientPower = 1;
41  PartialDerivativeTerm normalComponent(gradientComponent,
42  magnitudeOfGradientPower);
43 
44  // Compute and store the magnitude of the gradient of the function.
45  m_magnitudeOfGradient = 0.0;
46 
47  for (int idir = 0; idir < dim; idir++)
48  {
49  Real firstPartial = a_implicitFunction->value(BASISV_TM<int,dim>(idir),a_point);
50  m_magnitudeOfGradient += firstPartial*firstPartial;
51  }
52 
53  m_magnitudeOfGradient = sqrt(m_magnitudeOfGradient);
54 
55  Real value = expand(a_multiIndex,
56  normalComponent,
57  a_point,
58  a_implicitFunction);
59 
60  return value;
61 }
62 
63 //get m_magnitudeOfGradient
65 {
66  return m_magnitudeOfGradient;
67 }
68 
69 // Expand and evaluate the multi-index partial derivative of a
70 // PartialDerivativeTerm recursively. If the multi-index is zero (i.e., no
71 // further derivatives) then simply evaluate the PartialDerivateTerm at
72 // "a_point" using the implicit function. If the multi-index isn't zerom,
73 // explicitly compute one partial derivative which is a sum of
74 // PartialDerivativeTerm's and call "expand" with each of these terms and a
75 // reduced multi-index (which will eventually be zero). The sum the results
76 // and return that sum.
77 template <int dim> Real NormalDerivative<dim>::expand(const IvDim & a_multiIndex,
78  const PartialDerivativeTerm & a_term,
79  const RvDim & a_point,
80  const IFSlicer<dim> * a_implicitFunction) const
81 {
82  Real value;
83 
84  int firstNonZero;
85  for (firstNonZero = 0; firstNonZero < dim; firstNonZero++)
86  {
87  if (a_multiIndex[firstNonZero] != 0)
88  {
89  break;
90  }
91  }
92 
93  // No more derivatives, evaluate the current term
94  if (firstNonZero == dim)
95  {
96  value = 1.0;
97 
98  // Evalute the needed partial derivatives and take the product of the
99  // specified powers
100  const DerivativeProduct& curDerivativeProduct = a_term.first;
101 
102  for (typename DerivativeProduct::const_iterator it=curDerivativeProduct.begin();
103  it != curDerivativeProduct.end();
104  ++it)
105  {
106  const IvDim& curMultiIndex = it->first;
107  const int& curExponent = it->second;
108 
109  // Evaluate a single partial derivative to its power
110  Real curValue = pow(a_implicitFunction->value(curMultiIndex,a_point),curExponent);
111 
112  value *= curValue;
113  }
114 
115  if (m_magnitudeOfGradient != 0.0)
116  {
117  // Divide by the magnitude of the gradient including the exponent
118  int curExponent = a_term.second;
119  value /= pow(m_magnitudeOfGradient,curExponent);
120  }
121  }
122  else
123  {
124  value = 0.0;
125 
126  // This is the (current) partial derivative we are going to apply to the
127  // current product of partial derivatives
128  IvDim curPartialDerivative = BASISV_TM<int,dim>(firstNonZero);
129 
130  // This is the remaining multi-index that we haven't applied
131  IvDim reducedMultiIndex = a_multiIndex;
132  reducedMultiIndex[firstNonZero]--;
133 
134  // Distribute the current partial derivative over the product of
135  // derivatives using the product rule
136  const DerivativeProduct& curDerivativeProduct = a_term.first;
137  int curExponentOfMagnitudeOfGradient = a_term.second;
138 
139  // Loop through each term in the product
140  for (typename DerivativeProduct::const_iterator it=curDerivativeProduct.begin();
141  it != curDerivativeProduct.end();
142  ++it)
143  {
144  // Get the current derivative multi-index and exponent
145  const IvDim& curMultiIndex = it->first;
146  const int& curExponent = it->second;
147 
148  // Create the next term in the product rule by copying the current
149  // product and take the current partial derivative of the current
150  // product term (including the exponent).
151  DerivativeProduct newDerivativeProduct = curDerivativeProduct;
152 
153  // Handle the exponent of the current product term
154  Real multiplier = curExponent;
155 
156  if (curExponent == 1)
157  {
158  // Erase the current product term if its current exponent is one
159  newDerivativeProduct.erase(curMultiIndex);
160  }
161  else
162  {
163  // Otherwise, decrement the exponent by one
164  newDerivativeProduct[curMultiIndex] -= 1;
165  }
166 
167  // Generate the new product term
168  IvDim newMultiIndex = curMultiIndex;
169  newMultiIndex += curPartialDerivative;
170 
171  // Put it into the product
172  newDerivativeProduct[newMultiIndex] += 1;
173 
174  // Put the new product together with magnitude of the gradient term
175  PartialDerivativeTerm newTerm(newDerivativeProduct,
176  curExponentOfMagnitudeOfGradient);
177 
178  // Evaluate this term in the product rule (recursively)
179  Real curValue = multiplier * expand(reducedMultiIndex,
180  newTerm,
181  a_point,
182  a_implicitFunction);
183 
184  // Add the result into the overall product rule sum
185  value += curValue;
186  }
187 
188  // Now handle the last term in the overall product (and product rule)
189  // which is the inverse of the magnitude of the gradient with an exponent.
190 
191  // The derivative of the magnitude of the gradient results in a sum which
192  // has to be handled term by term
193  for (int idir = 0; idir < dim; idir++)
194  {
195  // Copy the current overall product
196  DerivativeProduct newDerivativeProduct = curDerivativeProduct;
197 
198  // Create the two new terms in the product
199  IvDim firstPartial = BASISV_TM<int,dim>(idir);
200  IvDim secondPartial = firstPartial + curPartialDerivative;
201 
202  // Put them in the new overall product
203  newDerivativeProduct[firstPartial] += 1;
204  newDerivativeProduct[secondPartial] += 1;
205 
206  // Generate the new exponent for the magnitude of the gradient
207  int newExponentOfMagnitudeOfGradient = curExponentOfMagnitudeOfGradient + 2;
208 
209  // The multiplier due to the exponent
210  Real multiplier = -curExponentOfMagnitudeOfGradient;
211 
212  // Put the new product together with magnitude of the gradient term
213  PartialDerivativeTerm newTerm(newDerivativeProduct,
214  newExponentOfMagnitudeOfGradient);
215 
216  // Evaluate this term in the product rule (recursively)
217  Real curValue = multiplier * expand(reducedMultiIndex,
218  newTerm,
219  a_point,
220  a_implicitFunction);
221 
222  // Add the result into the overall product rule sum
223  value += curValue;
224  }
225  }
226 
227  return value;
228 }
229 
230 #include "NamespaceFooter.H"
231 
232 #endif
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
pair< DerivativeProduct, int > PartialDerivativeTerm
Definition: NormalDerivative.H:46
Definition: IndexTM.H:36
NormalDerivative()
Null constructor.
Definition: NormalDerivativeImplem.H:17
double Real
Definition: REAL.H:33
virtual ~NormalDerivative()
Destructor.
Definition: NormalDerivativeImplem.H:22
Definition: IFSlicer.H:27
virtual Real value(const IndexTM< int, dim > &a_partialDerivative, const IndexTM< Real, dim > &a_point) const
Return the partial derivative evaluated at a_point.
Definition: IFSlicerImplem.H:60
Real getMagnitudeOfGradient()
Definition: NormalDerivativeImplem.H:64
map< IvDim, int, LexLT< IvDim > > DerivativeProduct
Definition: NormalDerivative.H:37
Real expand(const IvDim &a_multiIndex, const PartialDerivativeTerm &a_term, const RvDim &a_point, const IFSlicer< dim > *a_ifSlicer) const
Definition: NormalDerivativeImplem.H:77
int dim
Definition: EBInterface.H:146