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