Chombo + EB  3.2
NormalDerivativeNew.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 _NORMALDERIVATIVENEW_H_
12 #define _NORMALDERIVATIVENEW_H_
13 
14 #include "MayDay.H"
15 
16 #include "Notation.H"
17 #include "BoxIterator.H"
18 #include "IFSlicer.H"
19 
20 #include "NamespaceHeader.H"
21 
22 /// This computes the derivatives of the normal of a sliced implicit function
23 /**
24  This computes the derivatives of the normal of a sliced implicit function
25  */
26 template <int dim> class NormalDerivativeNew
27 {
28 public:
31 
32  // These are maps to calculate partial derivatives, based on
33  // multi-index IvDim, and of the normal, with components RvDim
34  typedef map<IvDim,Real,LexLT<IvDim> > ScalarPowerMap;
35  typedef map<IvDim,RvDim,LexLT<IvDim> > NormalDerivativeMap;
36 
37  /// Null constructor
38  /**
39  Null constructor
40  */
42 
43  /// Destructor
44  /**
45  Destructor
46  */
47  virtual ~NormalDerivativeNew();
48 
49  /// Evaluate the normal of a BaseIF subclass
50  /**
51  Evaluate the normal.
52  a_direction specifies which component of the normal use,
53  a_point is the point in space to evaluate the derivative, and
54  a_ifSlicer is a sliced function whose gradient is the normal.
55  */
56  Real normal(const int & a_direction,
57  const RvDim & a_point,
58  const IFSlicer<dim> * a_implicitFunction) const
59  {
60  // Compute and store the magnitude of the gradient of the function.
61  Real lenGradIF = 0.0;
62  Real value = 0.0;
63 
64  for (int idir = 0; idir < dim; idir++)
65  {
66  Real dxd = a_implicitFunction->
67  value(BASISV_TM<int,dim>(idir), a_point);
68  if (idir == a_direction)
69  value = dxd;
70  lenGradIF += dxd*dxd;
71  }
72 
73  return value / sqrt(lenGradIF);
74  }
75 
76  /// Evaluate all derivatives of the normal of an IFSlicer class
77  /**
78  Evaluate all derivatives of the normal,
79  a_maxP specifies the max sum of derivatives to
80  take in any coordinate direction,
81  a_point is the point in space to evaluate the derivatives, and
82  a_ifSlicer is a sliced function whose gradient is the normal.
83  */
84  NormalDerivativeMap calculateAll(const int& a_maxP,
85  const RvDim& a_point,
86  const IFSlicer<dim>* a_implicitFunction)
87  {
88  // pout() << "Normal derivatives<" << dim << "> for x=" << a_point << endl;
89 
90  ScalarPowerMap phiDerivs; // partials of \phi
91  ScalarPowerMap lDerivs; // partials of l = |\grad \phi|
92  NormalDerivativeMap nDerivs; // partials of the normal, n = \grad \phi / l
93 
94  // We only want the powers that would be less than
95  // the power in each direction and less than the total
96  Box b(IntVect::Zero, a_maxP*IntVect::Unit);
97  BoxIterator bit(b);
98  for (bit.begin(); bit.ok(); ++bit)
99  {
100  IvDim pow;
101  for(int idir = 0; idir < dim; idir++)
102  {
103  pow[idir] = bit()[idir];
104  }
105  if (pow.sum() <= a_maxP)
106  nDerivs[pow] = RvDim::Zero;
107  }
108 
109  // Calculate the phi derivatives we'll need
110  calculatePhiDerivs(phiDerivs, a_point, a_implicitFunction, a_maxP+1);
111 
112  // In order of increasing derivatives, in slices
113  typename NormalDerivativeMap::const_iterator it;
114  for (it = nDerivs.begin(); it != nDerivs.end(); ++it)
115  {
116  IvDim p = it->first;
117  addLDeriv(lDerivs, phiDerivs, nDerivs, p);
118  addNDeriv(nDerivs, phiDerivs, lDerivs, p);
119  }
120  /*
121  for (it = nDerivs.begin(); it != nDerivs.end(); ++it)
122  pout() << "New algorithm Normal derivative: p=" << it->first << ", " <<
123  it->second << endl;
124  */
125 
126  return nDerivs;
127  }
128 
129 protected:
130 
132  ScalarPowerMap& a_phiDerivs,
133  const RvDim& a_point,
134  const IFSlicer<dim>* a_implicitFunction,
135  const int& a_maxPartial)
136  {
137  Box b(IntVect::Zero, a_maxPartial*IntVect::Unit);
138  BoxIterator bit(b);
139  // We only want the powers that would be less than
140  // the power in each direction and less than the total
141  for (bit.begin(); bit.ok(); ++bit)
142  {
143  if (bit().sum() <= a_maxPartial)
144  {
145  IvDim deriv;
146  for(int idir = 0; idir < dim; idir++)
147  {
148  deriv[idir]= bit()[idir];
149  }
150 
151  a_phiDerivs[deriv] = a_implicitFunction->value(deriv, a_point);
152  }
153  }
154  }
155 
156 
157  void addLDeriv(
158  ScalarPowerMap& a_lDerivs,
159  ScalarPowerMap& a_phiDerivs,
160  NormalDerivativeMap& a_nDerivs,
161  const IvDim& a_deriv)
162  {
163  // If this is the first call, just calc l = |\grad phi|
164  if (a_deriv == IvDim::Zero)
165  {
166  RvDim dphi;
167  Real len = 0;
168  for (int d=0; d < dim; ++d)
169  {
170  Real val = a_phiDerivs[BASISV_TM<int,dim>(d)];
171  dphi[d] = val;
172  len += val*val;
173  }
174  a_lDerivs[a_deriv] = sqrt(len);
175  return;
176  }
177 
178  // Otherwise, assume we are moving in increasing x, then y...
179  int dir;
180  for (dir=0; dir < dim-1; ++dir)
181  if (a_deriv[dir] > 0)
182  break;
183 
184  Real dl = 0;
185  IvDim q = a_deriv - BASISV_TM<int,dim>(dir);
186  typename NormalDerivativeMap::const_iterator it;
187  for (it = a_nDerivs.begin(); it != a_nDerivs.end(); ++it)
188  {
189  IvDim r = it->first;
190  if (r <= q)
191  {
192  Real coef = nChoosek(q,r);
193  RvDim dn = a_nDerivs[r];
194  for (int d=0; d < dim; ++d)
195  {
196  IvDim dphip = q - r + BASISV_TM<int,dim>(d) +
197  BASISV_TM<int,dim>(dir);
198  dl += coef * dn[d] * a_phiDerivs[dphip];
199  }
200  }
201  }
202  a_lDerivs[a_deriv] = dl;
203  }
204 
205 
206  void addNDeriv(
207  NormalDerivativeMap& a_nDerivs,
208  ScalarPowerMap& a_phiDerivs,
209  ScalarPowerMap& a_lDerivs,
210  const IvDim& a_deriv)
211  {
212  RvDim dn = RvDim::Zero;
213  typename NormalDerivativeMap::const_iterator it;
214  for (it = a_nDerivs.begin(); it != a_nDerivs.end(); ++it)
215  {
216  IvDim q = it->first;
217  if ((q <= a_deriv) && (q != a_deriv))
218  {
219  Real coef = nChoosek(a_deriv,q);
220  dn += (coef * a_lDerivs[a_deriv-q]) * it->second;
221  }
222  }
223 
224  Real l = a_lDerivs[IvDim::Zero];
225  RvDim dphi;
226  for (int d=0; d < dim; ++d)
227  dphi[d] = a_phiDerivs[a_deriv + BASISV_TM<int,dim>(d)];
228  a_nDerivs[a_deriv] = (dphi - dn) / l;
229  }
230 
231 
232  // Calculates the multi-index binomial coefficient,
233  // product of "n choose k" for each index
234  static int nChoosek(const IvDim& a_n, const IvDim& a_k)
235  {
236  int numer = 1;
237  int denom = 1;
238 
239  IvDim diff = a_n - a_k;
240  CH_assert(diff >= IvDim::Zero);
241  for (int d = 0; d < dim; ++d)
242  {
243  for (int i=diff[d]+1; i <= a_n[d]; ++i)
244  numer *= i;
245  for (int i=2; i <= a_k[d]; ++i)
246  denom *= i;
247  }
248 
249  return numer / denom;
250  }
251 
252 private:
254  {
255  MayDay::Abort("NormalDerivativeNew doesn't allow copy construction");
256  }
257 
258  void operator=(const NormalDerivativeNew& a_input)
259  {
260  MayDay::Abort("NormalDerivativeNew doesn't allow assignment");
261  }
262 };
263 
264 #include "NamespaceFooter.H"
265 
267 
268 #endif
bool ok()
Definition: BoxIterator.H:281
#define CH_assert(cond)
Definition: CHArray.H:37
Real normal(const int &a_direction, const RvDim &a_point, const IFSlicer< dim > *a_implicitFunction) const
Evaluate the normal of a BaseIF subclass.
Definition: NormalDerivativeNew.H:56
NormalDerivativeNew(const NormalDerivativeNew &a_input)
Definition: NormalDerivativeNew.H:253
IndexTM< int, dim > IvDim
Definition: NormalDerivativeNew.H:29
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
map< IvDim, Real, LexLT< IvDim > > ScalarPowerMap
Definition: NormalDerivativeNew.H:34
Definition: IndexTM.H:36
void operator=(const NormalDerivativeNew &a_input)
Definition: NormalDerivativeNew.H:258
NormalDerivativeNew()
Null constructor.
Definition: NormalDerivativeNewImplem.H:17
static const IntVect Unit
Definition: IntVect.H:659
void addNDeriv(NormalDerivativeMap &a_nDerivs, ScalarPowerMap &a_phiDerivs, ScalarPowerMap &a_lDerivs, const IvDim &a_deriv)
Definition: NormalDerivativeNew.H:206
This computes the derivatives of the normal of a sliced implicit function.
Definition: NormalDerivativeNew.H:26
double Real
Definition: REAL.H:33
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
virtual ~NormalDerivativeNew()
Destructor.
Definition: NormalDerivativeNewImplem.H:22
static const IntVect Zero
Definition: IntVect.H:654
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
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
void begin()
Definition: BoxIterator.H:150
Definition: IFSlicer.H:27
T sum() const
Definition: IndexTMI.H:260
static int nChoosek(const IvDim &a_n, const IvDim &a_k)
Definition: NormalDerivativeNew.H:234
map< IvDim, RvDim, LexLT< IvDim > > NormalDerivativeMap
Definition: NormalDerivativeNew.H:35
IndexTM< Real, dim > RvDim
Definition: NormalDerivativeNew.H:30
int dim
Definition: EBInterface.H:146
void addLDeriv(ScalarPowerMap &a_lDerivs, ScalarPowerMap &a_phiDerivs, NormalDerivativeMap &a_nDerivs, const IvDim &a_deriv)
Definition: NormalDerivativeNew.H:157
static const IndexTM Zero
Definition: IndexTM.H:467
void calculatePhiDerivs(ScalarPowerMap &a_phiDerivs, const RvDim &a_point, const IFSlicer< dim > *a_implicitFunction, const int &a_maxPartial)
Definition: NormalDerivativeNew.H:131
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).