Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Evaluator.hxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Evaluator.hxx
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Inline header used to evaluate functional.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #ifndef UTILITIES_NEKMESH_NODEOPTI_EVALUATOR
37 #define UTILITIES_NEKMESH_NODEOPTI_EVALUATOR
38 
39 namespace Nektar
40 {
41 namespace Utilities
42 {
43 
44 using namespace std;
45 
46 /**
47  * @brief Calculate determinant of input matrix.
48  *
49  * Specialised versions of this function exist only for 2x2 and 3x3 matrices.
50  *
51  * @param jac Input matrix
52  *
53  * @return Jacobian of @p jac.
54  */
55 template <int DIM> inline NekDouble Determinant(NekDouble jac[][DIM])
56 {
57  return 0.0;
58 }
59 
60 template <> inline NekDouble Determinant<2>(NekDouble jac[][2])
61 {
62  return jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
63 }
64 
65 template <> inline NekDouble Determinant<3>(NekDouble jac[][3])
66 {
67  return jac[0][0] * (jac[1][1] * jac[2][2] - jac[2][1] * jac[1][2]) -
68  jac[0][1] * (jac[1][0] * jac[2][2] - jac[1][2] * jac[2][0]) +
69  jac[0][2] * (jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0]);
70 }
71 
72 /**
73  * @brief Calculate inverse transpose of input matrix.
74  *
75  * Specialised versions of this function exist only for 2x2 and 3x3 matrices.
76  *
77  * @param in Input matrix \f$ A \f$
78  * @param out Output matrix \f$ A^{-\top} \f$
79  */
80 template <int DIM>
81 inline void InvTrans(NekDouble in[][DIM], NekDouble out[][DIM])
82 {
83 }
84 
85 template <> inline void InvTrans<2>(NekDouble in[][2], NekDouble out[][2])
86 {
87  NekDouble invDet = 1.0 / Determinant<2>(in);
88 
89  out[0][0] = in[1][1] * invDet;
90  out[1][0] = -in[0][1] * invDet;
91  out[0][1] = -in[1][0] * invDet;
92  out[1][1] = in[0][0] * invDet;
93 }
94 
95 template <> inline void InvTrans<3>(NekDouble in[][3], NekDouble out[][3])
96 {
97  NekDouble invdet = 1.0 / Determinant<3>(in);
98 
99  out[0][0] = (in[1][1] * in[2][2] - in[2][1] * in[1][2]) * invdet;
100  out[1][0] = -(in[0][1] * in[2][2] - in[0][2] * in[2][1]) * invdet;
101  out[2][0] = (in[0][1] * in[1][2] - in[0][2] * in[1][1]) * invdet;
102  out[0][1] = -(in[1][0] * in[2][2] - in[1][2] * in[2][0]) * invdet;
103  out[1][1] = (in[0][0] * in[2][2] - in[0][2] * in[2][0]) * invdet;
104  out[2][1] = -(in[0][0] * in[1][2] - in[1][0] * in[0][2]) * invdet;
105  out[0][2] = (in[1][0] * in[2][1] - in[2][0] * in[1][1]) * invdet;
106  out[1][2] = -(in[0][0] * in[2][1] - in[2][0] * in[0][1]) * invdet;
107  out[2][2] = (in[0][0] * in[1][1] - in[1][0] * in[0][1]) * invdet;
108 }
109 
110 /**
111  * @brief Calculate Scalar product of input vectors.
112  */
113 /*template<int DIM>
114 inline NekDouble ScalarProd(NekDouble in1[DIM], NekDouble in2[DIM])
115 {
116  return 0.0;
117 }
118 
119 template <> inline NekDouble ScalarProd<2>(NekDouble in1[2], NekDouble in2[2])
120 {
121  return in1[0] * in2[0]
122  + in1[1] * in2[1];
123 }
124 
125 template <> inline NekDouble ScalarProd<3>(NekDouble in1[3], NekDouble in2[3])
126 {
127  return in1[0] * in2[0]
128  + in1[1] * in2[1]
129  + in1[2] * in2[2];
130 }*/
131 template<int DIM>
132 inline NekDouble ScalarProd(NekDouble (&in1)[DIM], NekDouble (&in2)[DIM])
133 {
134  return 0.0;
135 }
136 
137 template<>
138 inline NekDouble ScalarProd<2>(NekDouble (&in1)[2], NekDouble (&in2)[2])
139 {
140  return in1[0] * in2[0]
141  + in1[1] * in2[1];
142 }
143 template<>
144 inline NekDouble ScalarProd<3>(NekDouble (&in1)[3], NekDouble (&in2)[3])
145 {
146  return in1[0] * in2[0]
147  + in1[1] * in2[1]
148  + in1[2] * in2[2];
149 }
150 
151 
152 /**
153  * @brief Calculate \f$ E = F^\top F - I \f$ tensor used in derivation of linear
154  * elasticity gradients.
155  *
156  * Specialised versions of this function exist only for 2x2 and 3x3 matrices.
157  *
158  * @param in Input matrix \f$ F \f$
159  * @param out Output matrix \f$ F^\top F - I \f$
160  */
161 template <int DIM>
162 inline void EMatrix(NekDouble in[][DIM], NekDouble out[][DIM])
163 {
164 }
165 
166 template <> inline void EMatrix<2>(NekDouble in[][2], NekDouble out[][2])
167 {
168  out[0][0] = 0.5 * (in[0][0] * in[0][0] + in[1][0] * in[1][0] - 1.0);
169  out[1][0] = 0.5 * (in[0][0] * in[0][1] + in[1][0] * in[1][1]);
170  out[0][1] = 0.5 * (in[0][0] * in[0][1] + in[1][0] * in[1][1]);
171  out[1][1] = 0.5 * (in[1][1] * in[1][1] + in[0][1] * in[0][1] - 1.0);
172 }
173 
174 template <> inline void EMatrix<3>(NekDouble in[][3], NekDouble out[][3])
175 {
176  out[0][0] = 0.5 * (in[0][0] * in[0][0] + in[1][0] * in[1][0] +
177  in[2][0] * in[2][0] - 1.0);
178  out[1][0] = 0.5 * (in[0][0] * in[1][0] + in[1][0] * in[1][1] +
179  in[2][0] * in[2][1]);
180  out[0][1] = out[1][0];
181  out[2][0] = 0.5 * (in[0][0] * in[0][2] + in[1][0] * in[1][2] +
182  in[2][0] * in[2][2]);
183  out[0][2] = out[2][0];
184  out[1][1] = 0.5 * (in[0][1] * in[0][1] + in[1][1] * in[1][1] +
185  in[2][1] * in[2][1] - 1.0);
186  out[1][2] = 0.5 * (in[0][1] * in[0][2] + in[1][1] * in[1][2] +
187  in[2][1] * in[2][2]);
188  out[2][1] = out[1][2];
189  out[2][2] = 0.5 * (in[0][2] * in[0][2] + in[1][2] * in[1][2] +
190  in[2][2] * in[2][2] - 1.0);
191 }
192 
193 
194 
195 /**
196  * @brief Calculate Frobenius inner product of input matrices.
197  */
198 template<int DIM>
199 inline NekDouble FrobProd(NekDouble in1[][DIM],
200  NekDouble in2[][DIM])
201 {
202  return 0.0;
203 }
204 
205 template<>
206 inline NekDouble FrobProd<2>(NekDouble in1[][2], NekDouble in2[][2])
207 {
208  return in1[0][0] * in2[0][0]
209  + in1[0][1] * in2[0][1]
210  + in1[1][0] * in2[1][0]
211  + in1[1][1] * in2[1][1] ;
212 }
213 
214 template<>
215 inline NekDouble FrobProd<3>(NekDouble in1[][3], NekDouble in2[][3])
216 {
217  return in1[0][0] * in2[0][0]
218  + in1[0][1] * in2[0][1]
219  + in1[0][2] * in2[0][2]
220  + in1[1][0] * in2[1][0]
221  + in1[1][1] * in2[1][1]
222  + in1[1][2] * in2[1][2]
223  + in1[2][0] * in2[2][0]
224  + in1[2][1] * in2[2][1]
225  + in1[2][2] * in2[2][2] ;
226 }
227 
228 // Typedef for derivative storage, we use boost::multi_array so we can pass this
229 // to functions easily
230 typedef boost::multi_array<NekDouble, 4> DerivArray;
231 
232 
233 /**
234  * @brief Calculate Frobenius norm \f$ \| A \|_f ^2 \f$ of a matrix \f$ A \f$.
235  *
236  * @param inarray Input matrix \f$ A \f$
237  */
238 template<int DIM>
239 inline NekDouble FrobeniusNorm(NekDouble inarray[][DIM])
240 {
241  return 0.0;
242 }
243 
244 template<>
245 inline NekDouble FrobeniusNorm<2>(NekDouble inarray[][2])
246 {
247  return inarray[0][0] * inarray[0][0]
248  + inarray[0][1] * inarray[0][1]
249  + inarray[1][0] * inarray[1][0]
250  + inarray[1][1] * inarray[1][1] ;
251 }
252 
253 template<>
254 inline NekDouble FrobeniusNorm<3>(NekDouble inarray[][3])
255 {
256  return inarray[0][0] * inarray[0][0]
257  + inarray[0][1] * inarray[0][1]
258  + inarray[0][2] * inarray[0][2]
259  + inarray[1][0] * inarray[1][0]
260  + inarray[1][1] * inarray[1][1]
261  + inarray[1][2] * inarray[1][2]
262  + inarray[2][0] * inarray[2][0]
263  + inarray[2][1] * inarray[2][1]
264  + inarray[2][2] * inarray[2][2] ;
265 }
266 
267 /**
268  * @brief Evaluate functional for elements connected to a node.
269  *
270  * @param minJacNew Stores current minimum Jacobian for the element group
271  * @param gradient If true, calculate gradient.
272  */
273 template <int DIM>
274 NekDouble NodeOpti::GetFunctional(NekDouble &minJacNew, bool gradient)
275 {
276  map<LibUtilities::ShapeType, vector<ElUtilSharedPtr> >::iterator typeIt;
277  map<LibUtilities::ShapeType, DerivArray> derivs;
278 
279  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
280  {
281  int nElmt = typeIt->second.size();
282  int totpts = m_derivUtils[typeIt->first]->ptsStd * nElmt;
283  NekDouble* X = new NekDouble[DIM*totpts];
284 
285  // Store x/y components of each element sequentially in memory
286  for (int i = 0, cnt = 0; i < nElmt; ++i)
287  {
288  for (int j = 0; j < m_derivUtils[typeIt->first]->ptsStd; ++j)
289  {
290  for (int d = 0; d < DIM; ++d)
291  {
292  X[cnt + d * m_derivUtils[typeIt->first]->ptsStd + j] =
293  *(typeIt->second[i]->nodes[j][d]);
294  }
295  }
296  cnt += DIM * m_derivUtils[typeIt->first]->ptsStd;
297  }
298 
299  // Storage for derivatives, ordered by:
300  // - standard coordinate direction
301  // - number of elements
302  // - cartesian coordinate direction
303  // - quadrature points
304  derivs.insert(std::make_pair(
305  typeIt->first,
306  DerivArray(boost::extents[DIM][nElmt][DIM]
307  [m_derivUtils[typeIt->first]->pts])));
308 
309  // Calculate x- and y-gradients
310  for (int d = 0; d < DIM; ++d)
311  {
312  Blas::Dgemm('N', 'N', m_derivUtils[typeIt->first]->pts, DIM * nElmt,
313  m_derivUtils[typeIt->first]->ptsStd, 1.0,
314  m_derivUtils[typeIt->first]->VdmD[d].GetRawPtr(),
315  m_derivUtils[typeIt->first]->pts, X,
316  m_derivUtils[typeIt->first]->ptsStd, 0.0,
317  &derivs[typeIt->first][d][0][0][0],
318  m_derivUtils[typeIt->first]->pts);
319  }
320  }
321 
322  minJacNew = std::numeric_limits<double>::max();
323  NekDouble integral = 0.0;
324  NekDouble ep =
325  m_minJac < 0.0 ? sqrt(1e-8 + 0.04 * m_minJac * m_minJac) : 1e-4;
326  NekDouble jacIdeal[DIM][DIM], jacDet;
327  m_grad = Array<OneD, NekDouble>(DIM == 2 ? 5 : 9, 0.0);
328 
329  switch (m_opti)
330  {
331  case eLinEl:
332  {
333  const NekDouble nu = 0.45;
334  const NekDouble mu = 1.0 / 2.0 / (1.0 + nu);
335  const NekDouble K = 1.0 / 3.0 / (1.0 - 2.0 * nu);
336 
337  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
338  {
339  NekVector<NekDouble> &quadW =
340  m_derivUtils[typeIt->first]->quadW;
341  for (int i = 0; i < typeIt->second.size(); ++i)
342  {
343  for (int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
344  {
345  NekDouble phiM[DIM][DIM];
346  for (int l = 0; l < DIM; ++l)
347  {
348  for (int n = 0; n < DIM; ++n)
349  {
350  phiM[n][l] =
351  derivs[typeIt->first][l][i][n][k];
352  }
353  }
354  // begin CalcIdealJac
355  for (int m = 0; m < DIM; ++m)
356  {
357  for (int n = 0; n < DIM; ++n)
358  {
359  jacIdeal[n][m] = 0.0;
360  for (int l = 0; l < DIM; ++l)
361  {
362  jacIdeal[n][m] += phiM[n][l] *
363  typeIt->second[i]->maps[k][m * 3 + l];
364  }
365  }
366  }
367  jacDet = Determinant(jacIdeal);
368  // end CalcIdealJac
369 
370  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
371  minJacNew = min(minJacNew, jacDet);
372 
373  NekDouble Emat[DIM][DIM];
374  EMatrix<DIM>(jacIdeal, Emat);
375 
376  NekDouble trEtE = FrobProd<DIM>(Emat, Emat);
377  NekDouble sigma =
378  0.5 *
379  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
380 
381  if (sigma < numeric_limits<double>::min() && !gradient)
382  {
383  return numeric_limits<double>::max();
384  }
385  ASSERTL0(sigma > numeric_limits<double>::min(),
386  std::string("dividing by zero ") +
387  boost::lexical_cast<string>(sigma) + " " +
388  boost::lexical_cast<string>(jacDet) + " " +
389  boost::lexical_cast<string>(ep));
390 
391  NekDouble lsigma = log(sigma);
392  integral += quadW[k] *
393  absIdealMapDet *
394  (K * 0.5 * lsigma * lsigma + mu * trEtE);
395 
396  if (gradient)
397  {
398  NekDouble jacDerivPhi[DIM];
399  NekDouble jacDetDeriv[DIM];
400 
401  NekDouble derivDet = Determinant<DIM>(phiM);
402  NekDouble jacInvTrans[DIM][DIM];
403  InvTrans<DIM>(phiM, jacInvTrans);
404 
405  NekDouble basisDeriv[DIM];
406  for (int m = 0; m < DIM; ++m)
407  {
408  basisDeriv[m] = *(
409  m_derivUtils[typeIt->first]->VdmD[m])(
410  k, typeIt->second[i]->NodeId(m_node->m_id));
411  }
412  // jacDeriv is actually a tensor,
413  // but can be stored as a vector, as 18 out of 27 entries are zero
414  // and the other 9 entries are three triplets
415  // this is due to the delta function in jacDeriv
416  NekDouble jacDeriv[DIM];
417  for (int l = 0; l < DIM; ++l)
418  {
419  jacDeriv[l] = basisDeriv[l];
420  }
421 
422  // jacDerivPhi is actually a tensor,
423  // but can be stored as a vector due to the simple form of jacDeriv
424  for (int n = 0; n < DIM; ++n)
425  {
426  jacDerivPhi[n] = 0.0;
427  for (int l = 0; l < DIM; ++l)
428  {
429  jacDerivPhi[n] += jacDeriv[l] *
430  typeIt->second[i]->maps[k][l + 3 * n];
431  }
432  }
433 
434  for (int m = 0; m < DIM; ++m)
435  {
436  jacDetDeriv[m] = 0.0;
437  for (int n = 0; n < DIM; ++n)
438  {
439  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
440  }
441  jacDetDeriv[m] *= derivDet / absIdealMapDet;
442  }
443  // end of common part to all four versionsNekDouble
444 
445 
446  NekDouble M2[DIM][DIM][DIM];
447  // use the delta function in jacDeriv and do some tensor calculus
448  // to come up with this simplified expression for:
449  // LEM2<DIM>(jacIdeal, jacDerivPhi, M2);
450  for(int d = 0; d < DIM; d++)
451  {
452  for (int m = 0; m < DIM; ++m)
453  {
454  for (int n = 0; n < DIM; ++n)
455  {
456  M2[d][m][n] = 0.5*(jacDerivPhi[m] * jacIdeal[d][n]
457  + jacIdeal[d][m] * jacDerivPhi[n]);
458  }
459  }
460  }
461 
462  for (int m = 0; m < DIM; ++m)
463  {
464  NekDouble frobProdA = FrobProd<DIM>(M2[m], Emat);
465 
466  m_grad[m] +=
467  quadW[k] * absIdealMapDet *
468  (2.0 * mu * frobProdA +
469  K * lsigma * jacDetDeriv[m] /
470  (2.0 * sigma - jacDet));
471  }
472 
473  int ct = 0;
474  for (int m = 0; m < DIM; ++m)
475  {
476  for (int l = m; l < DIM; ++l, ct++)
477  {
478  NekDouble frobProdBC = FrobProd<DIM>(M2[m], M2[l]);
479  NekDouble M3[DIM][DIM];
480  // use the delta function in jacDeriv and do some tensor calculus
481  // to come up with this simplified expression for:
482  // LEM3<DIM>(jacDerivPhi, M3);
483  if (m == l)
484  {
485  for (int p = 0; p < DIM; ++p)
486  {
487  for (int q = 0; q < DIM; ++q)
488  {
489  M3[p][q] = jacDerivPhi[p] * jacDerivPhi[q];
490  }
491  }
492  frobProdBC += FrobProd<DIM>(M3,Emat);
493  }
494 
495  m_grad[ct + DIM] +=
496  quadW[k] * absIdealMapDet *
497  (2.0 * mu * frobProdBC +
498  jacDetDeriv[m] * jacDetDeriv[l] * K /
499  (2.0 * sigma - jacDet) /
500  (2.0 * sigma - jacDet) *
501  (1.0 - jacDet * lsigma /
502  (2.0 * sigma - jacDet)));
503  }
504  }
505  }
506  }
507  }
508  }
509  break;
510  }
511 
512  case eHypEl:
513  {
514  const NekDouble nu = 0.45;
515  const NekDouble mu = 1.0 / 2.0 / (1.0 + nu);
516  const NekDouble K = 1.0 / 3.0 / (1.0 - 2.0 * nu);
517 
518  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
519  {
520  NekVector<NekDouble> &quadW =
521  m_derivUtils[typeIt->first]->quadW;
522  for (int i = 0; i < typeIt->second.size(); ++i)
523  {
524  for (int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
525  {
526  NekDouble phiM[DIM][DIM];
527  for (int l = 0; l < DIM; ++l)
528  {
529  for (int n = 0; n < DIM; ++n)
530  {
531  phiM[n][l] =
532  derivs[typeIt->first][l][i][n][k];
533  }
534  }
535  // begin CalcIdealJac
536  for (int m = 0; m < DIM; ++m)
537  {
538  for (int n = 0; n < DIM; ++n)
539  {
540  jacIdeal[n][m] = 0.0;
541  for (int l = 0; l < DIM; ++l)
542  {
543  jacIdeal[n][m] += phiM[n][l] *
544  typeIt->second[i]->maps[k][m * 3 + l];
545  }
546  }
547  }
548  jacDet = Determinant(jacIdeal);
549  // end CalcIdealJac
550 
551  minJacNew = min(minJacNew, jacDet);
552 
553  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
554 
555  NekDouble I1 = FrobeniusNorm(jacIdeal);
556 
557  NekDouble sigma =
558  0.5 *
559  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
560 
561  if (sigma < numeric_limits<double>::min() && !gradient)
562  {
563  return numeric_limits<double>::max();
564  }
565 
566  ASSERTL0(sigma > numeric_limits<double>::min(),
567  std::string("dividing by zero ") +
568  boost::lexical_cast<string>(sigma) + " " +
569  boost::lexical_cast<string>(jacDet) + " " +
570  boost::lexical_cast<string>(ep));
571 
572  NekDouble lsigma = log(sigma);
573  integral += quadW[k] * absIdealMapDet *
574  (0.5 * mu * (I1 - 3.0 - 2.0 * lsigma) +
575  0.5 * K * lsigma * lsigma);
576 
577  // Derivative of basis function in each direction
578  if (gradient)
579  {
580  NekDouble jacDerivPhi[DIM];
581  NekDouble jacDetDeriv[DIM];
582 
583  NekDouble derivDet = Determinant<DIM>(phiM);
584  NekDouble jacInvTrans[DIM][DIM];
585  InvTrans<DIM>(phiM, jacInvTrans);
586 
587  NekDouble basisDeriv[DIM];
588  for (int m = 0; m < DIM; ++m)
589  {
590  basisDeriv[m] = *(
591  m_derivUtils[typeIt->first]->VdmD[m])(
592  k, typeIt->second[i]->NodeId(m_node->m_id));
593  }
594  // jacDeriv is actually a tensor,
595  // but can be stored as a vector, as 18 out of 27 entries are zero
596  // and the other 9 entries are three triplets
597  // this is due to the delta function in jacDeriv
598  NekDouble jacDeriv[DIM];
599  for (int l = 0; l < DIM; ++l)
600  {
601  jacDeriv[l] = basisDeriv[l];
602  }
603 
604  // jacDerivPhi is actually a tensor,
605  // but can be stored as a vector due to the simple form of jacDeriv
606  for (int n = 0; n < DIM; ++n)
607  {
608  jacDerivPhi[n] = 0.0;
609  for (int l = 0; l < DIM; ++l)
610  {
611  jacDerivPhi[n] += jacDeriv[l] *
612  typeIt->second[i]->maps[k][l + 3 * n];
613  }
614  }
615 
616  for (int m = 0; m < DIM; ++m)
617  {
618  jacDetDeriv[m] = 0.0;
619  for (int n = 0; n < DIM; ++n)
620  {
621  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
622  }
623  jacDetDeriv[m] *= derivDet / absIdealMapDet;
624  }
625  // end of common part to all four versionsNekDouble
626 
627  for (int m = 0; m < DIM; ++m)
628  {
629  // because of the zero entries of the tensor jacDerivPhi,
630  // the Frobenius-product becomes a scalar product
631  NekDouble frobProd =
632  ScalarProd<DIM>(jacIdeal[m],jacDerivPhi);
633 
634  m_grad[m] +=
635  quadW[k] * absIdealMapDet *
636  (mu * frobProd +
637  (jacDetDeriv[m] / (2.0 * sigma - jacDet) *
638  (K * lsigma - mu)));
639  }
640 
641  int ct = 0;
642  for (int m = 0; m < DIM; ++m)
643  {
644  for (int l = m; l < DIM; ++l, ct++)
645  {
646  NekDouble frobProdHes = 0.0;
647  // because of the zero entries of the tensor jacDerivPhi,
648  // the matrix frobProdHes has only diagonal entries
649  if (m == l)
650  {
651  // because of the zero entries of the tensor jacDerivPhi,
652  // the Frobenius-product becomes a scalar product
653  frobProdHes = ScalarProd<DIM>(jacDerivPhi,jacDerivPhi);
654  }
655 
656  m_grad[ct + DIM] +=
657  quadW[k] * absIdealMapDet *
658  (mu * frobProdHes +
659  jacDetDeriv[m] * jacDetDeriv[l] /
660  (2.0 * sigma - jacDet) /
661  (2.0 * sigma - jacDet) *
662  (K - jacDet * (K * lsigma - mu) /
663  (2.0 * sigma - jacDet)));
664  }
665  }
666  }
667  }
668  }
669  }
670  break;
671  }
672 
673  case eRoca:
674  {
675  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
676  {
677  NekVector<NekDouble> &quadW =
678  m_derivUtils[typeIt->first]->quadW;
679  for (int i = 0; i < typeIt->second.size(); ++i)
680  {
681  for (int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
682  {
683  NekDouble phiM[DIM][DIM];
684  for (int l = 0; l < DIM; ++l)
685  {
686  for (int n = 0; n < DIM; ++n)
687  {
688  phiM[n][l] =
689  derivs[typeIt->first][l][i][n][k];
690  }
691  }
692  // begin CalcIdealJac
693  for (int m = 0; m < DIM; ++m)
694  {
695  for (int n = 0; n < DIM; ++n)
696  {
697  jacIdeal[n][m] = 0.0;
698  for (int l = 0; l < DIM; ++l)
699  {
700  jacIdeal[n][m] += phiM[n][l] *
701  typeIt->second[i]->maps[k][m * 3 + l];
702  }
703  }
704  }
705  jacDet = Determinant(jacIdeal);
706  // end CalcIdealJac
707 
708  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
709  minJacNew = min(minJacNew, jacDet);
710  NekDouble frob = FrobeniusNorm(jacIdeal);
711  NekDouble sigma =
712  0.5 *
713  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
714 
715  if (sigma < numeric_limits<double>::min() && !gradient)
716  {
717  return numeric_limits<double>::max();
718  }
719 
720  ASSERTL0(sigma > numeric_limits<double>::min(),
721  std::string("dividing by zero ") +
722  boost::lexical_cast<string>(sigma) + " " +
723  boost::lexical_cast<string>(jacDet) + " " +
724  boost::lexical_cast<string>(ep));
725 
726  NekDouble W = frob / DIM / pow(fabs(sigma), 2.0 / DIM);
727  integral +=
728  quadW[k] * absIdealMapDet * W;
729 
730  // Derivative of basis function in each direction
731  if (gradient)
732  {
733  NekDouble jacDerivPhi[DIM];
734  NekDouble jacDetDeriv[DIM];
735 
736  NekDouble derivDet = Determinant<DIM>(phiM);
737  NekDouble jacInvTrans[DIM][DIM];
738  InvTrans<DIM>(phiM, jacInvTrans);
739 
740  NekDouble basisDeriv[DIM];
741  for (int m = 0; m < DIM; ++m)
742  {
743  basisDeriv[m] = *(
744  m_derivUtils[typeIt->first]->VdmD[m])(
745  k, typeIt->second[i]->NodeId(m_node->m_id));
746  }
747  // jacDeriv is actually a tensor,
748  // but can be stored as a vector, as 18 out of 27 entries are zero
749  // and the other 9 entries are three triplets
750  // this is due to the delta function in jacDeriv
751  NekDouble jacDeriv[DIM];
752  for (int l = 0; l < DIM; ++l)
753  {
754  jacDeriv[l] = basisDeriv[l];
755  }
756 
757  // jacDerivPhi is actually a tensor,
758  // but can be stored as a vector due to the simple form of jacDeriv
759  for (int n = 0; n < DIM; ++n)
760  {
761  jacDerivPhi[n] = 0.0;
762  for (int l = 0; l < DIM; ++l)
763  {
764  jacDerivPhi[n] += jacDeriv[l] *
765  typeIt->second[i]->maps[k][l + 3 * n];
766  }
767  }
768 
769  for (int m = 0; m < DIM; ++m)
770  {
771  jacDetDeriv[m] = 0.0;
772  for (int n = 0; n < DIM; ++n)
773  {
774  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
775  }
776  jacDetDeriv[m] *= derivDet / absIdealMapDet;
777  }
778  // end of common part to all four versionsNekDouble
779 
780  NekDouble frobProd[DIM];
781  NekDouble inc[DIM];
782  for (int m = 0; m < DIM; ++m)
783  {
784  // because of the zero entries of the tensor jacDerivPhi,
785  // the Frobenius-product becomes a scalar product
786  frobProd[m] = ScalarProd<DIM>(jacIdeal[m],jacDerivPhi);
787 
788  inc[m] =
789  quadW[k] * absIdealMapDet *
790  (2.0 * W * (frobProd[m] / frob -
791  jacDetDeriv[m] / DIM /
792  (2.0 * sigma - jacDet)));
793  m_grad[m] += inc[m];
794  }
795 
796 
797 
798  int ct = 0;
799  for (int m = 0; m < DIM; ++m)
800  {
801  for (int l = m; l < DIM; ++l, ct++)
802  {
803  NekDouble frobProdHes = 0.0;
804  // because of the zero entries of the tensor jacDerivPhi,
805  // the matrix frobProdHes has only diagonal entries
806  if (m == l)
807  {
808  // because of the zero entries of the tensor jacDerivPhi,
809  // the Frobenius-product becomes a scalar product
810  frobProdHes = ScalarProd<DIM>(jacDerivPhi,jacDerivPhi);
811  }
812 
813  m_grad[ct + DIM] +=
814  quadW[k] * absIdealMapDet *
815  (inc[m] * inc[l] / W +
816  2.0 * W *
817  (frobProdHes / frob -
818  2.0 * frobProd[m] * frobProd[l] /
819  frob / frob +
820  jacDetDeriv[m] * jacDetDeriv[l] *
821  jacDet /
822  (2.0 * sigma - jacDet) /
823  (2.0 * sigma - jacDet) /
824  (2.0 * sigma - jacDet) /
825  DIM));
826  }
827  }
828  }
829  }
830  }
831  }
832  break;
833  }
834 
835  case eWins:
836  {
837  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
838  {
839  NekVector<NekDouble> &quadW =
840  m_derivUtils[typeIt->first]->quadW;
841  for (int i = 0; i < typeIt->second.size(); ++i)
842  {
843  for (int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
844  {
845  NekDouble phiM[DIM][DIM];
846  for (int l = 0; l < DIM; ++l)
847  {
848  for (int n = 0; n < DIM; ++n)
849  {
850  phiM[n][l] =
851  derivs[typeIt->first][l][i][n][k];
852  }
853  }
854  // begin CalcIdealJac
855  for (int m = 0; m < DIM; ++m)
856  {
857  for (int n = 0; n < DIM; ++n)
858  {
859  jacIdeal[n][m] = 0.0;
860  for (int l = 0; l < DIM; ++l)
861  {
862  jacIdeal[n][m] += phiM[n][l] *
863  typeIt->second[i]->maps[k][m * 3 + l];
864  }
865  }
866  }
867  jacDet = Determinant(jacIdeal);
868  // end CalcIdealJac
869 
870  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
871  minJacNew = min(minJacNew, jacDet);
872  NekDouble frob = FrobeniusNorm(jacIdeal);
873  NekDouble sigma =
874  0.5 *
875  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
876 
877  if (sigma < numeric_limits<double>::min() && !gradient)
878  {
879  return numeric_limits<double>::max();
880  }
881 
882  ASSERTL0(sigma > numeric_limits<double>::min(),
883  std::string("dividing by zero ") +
884  boost::lexical_cast<string>(sigma) + " " +
885  boost::lexical_cast<string>(jacDet) + " " +
886  boost::lexical_cast<string>(ep));
887 
888  NekDouble W = frob / sigma;
889  integral +=
890  quadW[k] * absIdealMapDet * W;
891 
892  // Derivative of basis function in each direction
893  if (gradient)
894  {
895  NekDouble jacDerivPhi[DIM];
896  NekDouble jacDetDeriv[DIM];
897 
898  NekDouble derivDet = Determinant<DIM>(phiM);
899  NekDouble jacInvTrans[DIM][DIM];
900  InvTrans<DIM>(phiM, jacInvTrans);
901 
902  NekDouble basisDeriv[DIM];
903  for (int m = 0; m < DIM; ++m)
904  {
905  basisDeriv[m] = *(
906  m_derivUtils[typeIt->first]->VdmD[m])(
907  k, typeIt->second[i]->NodeId(m_node->m_id));
908  }
909  // jacDeriv is actually a tensor,
910  // but can be stored as a vector, as 18 out of 27 entries are zero
911  // and the other 9 entries are three triplets
912  // this is due to the delta function in jacDeriv
913  NekDouble jacDeriv[DIM];
914  for (int l = 0; l < DIM; ++l)
915  {
916  jacDeriv[l] = basisDeriv[l];
917  }
918 
919  // jacDerivPhi is actually a tensor,
920  // but can be stored as a vector due to the simple form of jacDeriv
921  for (int n = 0; n < DIM; ++n)
922  {
923  jacDerivPhi[n] = 0.0;
924  for (int l = 0; l < DIM; ++l)
925  {
926  jacDerivPhi[n] += jacDeriv[l] *
927  typeIt->second[i]->maps[k][l + 3 * n];
928  }
929  }
930 
931  for (int m = 0; m < DIM; ++m)
932  {
933  jacDetDeriv[m] = 0.0;
934  for (int n = 0; n < DIM; ++n)
935  {
936  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
937  }
938  jacDetDeriv[m] *= derivDet / absIdealMapDet;
939  }
940  // end of common part to all four versionsNekDouble
941 
942  NekDouble frobProd[DIM];
943  NekDouble inc[DIM];
944  for (int m = 0; m < DIM; ++m)
945  {
946  // because of the zero entries of the tensor jacDerivPhi,
947  // the Frobenius-product becomes a scalar product
948  frobProd[m] = ScalarProd<DIM>(jacIdeal[m],jacDerivPhi);
949 
950  inc[m] =
951  quadW[k] *
952  absIdealMapDet *
953  (W *
954  (2.0 * frobProd[m] / frob -
955  jacDetDeriv[m] / (2.0 * sigma - jacDet)));
956  m_grad[m] += inc[m];
957  }
958 
959  int ct = 0;
960  for (int m = 0; m < DIM; ++m)
961  {
962  for (int l = m; l < DIM; ++l, ct++)
963  {
964  NekDouble frobProdHes = 0.0;
965  // because of the zero entries of the tensor jacDerivPhi,
966  // the matrix frobProdHes has only diagonal entries
967  if (m == l)
968  {
969  // because of the zero entries of the tensor jacDerivPhi,
970  // the Frobenius-product becomes a scalar product
971  frobProdHes = ScalarProd<DIM>(jacDerivPhi,jacDerivPhi);
972  }
973 
974  m_grad[ct + DIM] +=
975  quadW[k] * absIdealMapDet *
976  (inc[m] * inc[l] / W +
977  2.0 * W *
978  (frobProdHes / frob -
979  2.0 * frobProd[m] * frobProd[l] /
980  frob / frob +
981  0.5 * jacDetDeriv[m] *
982  jacDetDeriv[l] * jacDet /
983  (2.0 * sigma - jacDet) /
984  (2.0 * sigma - jacDet) /
985  (2.0 * sigma - jacDet)));
986  }
987  }
988  }
989  }
990  }
991  }
992  break;
993  }
994  }
995 
996  // ASSERTL0(std::isfinite(integral),"inf in integral");
997 
998  return integral;
999  // return sqrt(m_grad[0]*m_grad[0] + m_grad[1]*m_grad[1]);
1000 }
1001 }
1002 }
1003 
1004 #endif
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::Utilities::InvTrans< 3 >
void InvTrans< 3 >(NekDouble in[][3], NekDouble out[][3])
Definition: Evaluator.hxx:95
Nektar::Utilities::InvTrans< 2 >
void InvTrans< 2 >(NekDouble in[][2], NekDouble out[][2])
Definition: Evaluator.hxx:85
Nektar::Utilities::FrobeniusNorm< 3 >
NekDouble FrobeniusNorm< 3 >(NekDouble inarray[][3])
Definition: Evaluator.hxx:254
Nektar::Utilities::ScalarProd< 3 >
NekDouble ScalarProd< 3 >(NekDouble(&in1)[3], NekDouble(&in2)[3])
Definition: Evaluator.hxx:144
Nektar::Utilities::EMatrix< 2 >
void EMatrix< 2 >(NekDouble in[][2], NekDouble out[][2])
Definition: Evaluator.hxx:166
std
STL namespace.
Nektar::Utilities::FrobProd< 3 >
NekDouble FrobProd< 3 >(NekDouble in1[][3], NekDouble in2[][3])
Definition: Evaluator.hxx:215
Nektar::Utilities::DerivArray
boost::multi_array< NekDouble, 4 > DerivArray
Definition: Evaluator.hxx:230
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
Nektar::Utilities::Determinant< 2 >
NekDouble Determinant< 2 >(NekDouble jac[][2])
Definition: Evaluator.hxx:60
Nektar::Utilities::Determinant
NekDouble Determinant(NekDouble jac[][DIM])
Calculate determinant of input matrix.
Definition: Evaluator.hxx:55
Nektar::Utilities::InvTrans
void InvTrans(NekDouble in[][DIM], NekDouble out[][DIM])
Calculate inverse transpose of input matrix.
Definition: Evaluator.hxx:81
Nektar::Utilities::FrobProd
NekDouble FrobProd(NekDouble in1[][DIM], NekDouble in2[][DIM])
Calculate Frobenius inner product of input matrices.
Definition: Evaluator.hxx:199
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::Utilities::FrobProd< 2 >
NekDouble FrobProd< 2 >(NekDouble in1[][2], NekDouble in2[][2])
Definition: Evaluator.hxx:206
Nektar::Utilities::FrobeniusNorm< 2 >
NekDouble FrobeniusNorm< 2 >(NekDouble inarray[][2])
Definition: Evaluator.hxx:245
Nektar::Utilities::NodeOpti::GetFunctional
NekDouble GetFunctional(NekDouble &minJacNew, bool gradient=true)
Evaluate functional for elements connected to a node.
Definition: Evaluator.hxx:274
Nektar::Utilities::ScalarProd
NekDouble ScalarProd(NekDouble(&in1)[DIM], NekDouble(&in2)[DIM])
Calculate Scalar product of input vectors.
Definition: Evaluator.hxx:132
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Nektar::Utilities::ScalarProd< 2 >
NekDouble ScalarProd< 2 >(NekDouble(&in1)[2], NekDouble(&in2)[2])
Definition: Evaluator.hxx:138
Nektar::Utilities::EMatrix
void EMatrix(NekDouble in[][DIM], NekDouble out[][DIM])
Calculate tensor used in derivation of linear elasticity gradients.
Definition: Evaluator.hxx:162
Nektar::Utilities::FrobeniusNorm
NekDouble FrobeniusNorm(NekDouble inarray[][DIM])
Calculate Frobenius norm of a matrix .
Definition: Evaluator.hxx:239
Nektar::Utilities::Determinant< 3 >
NekDouble Determinant< 3 >(NekDouble jac[][3])
Definition: Evaluator.hxx:65
Nektar::Utilities::EMatrix< 3 >
void EMatrix< 3 >(NekDouble in[][3], NekDouble out[][3])
Definition: Evaluator.hxx:174
Generated on Wed Mar 8 2017 12:28:54 for Nektar++ by   doxygen 1.8.8