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 /**
229  * @brief Calculate Frobenius norm \f$ \| A \|_f ^2 \f$ of a matrix \f$ A \f$.
230  *
231  * @param inarray Input matrix \f$ A \f$
232  */
233 template<int DIM>
234 inline NekDouble FrobeniusNorm(NekDouble inarray[][DIM])
235 {
236  return 0.0;
237 }
238 
239 template<>
240 inline NekDouble FrobeniusNorm<2>(NekDouble inarray[][2])
241 {
242  return inarray[0][0] * inarray[0][0]
243  + inarray[0][1] * inarray[0][1]
244  + inarray[1][0] * inarray[1][0]
245  + inarray[1][1] * inarray[1][1] ;
246 }
247 
248 template<>
249 inline NekDouble FrobeniusNorm<3>(NekDouble inarray[][3])
250 {
251  return inarray[0][0] * inarray[0][0]
252  + inarray[0][1] * inarray[0][1]
253  + inarray[0][2] * inarray[0][2]
254  + inarray[1][0] * inarray[1][0]
255  + inarray[1][1] * inarray[1][1]
256  + inarray[1][2] * inarray[1][2]
257  + inarray[2][0] * inarray[2][0]
258  + inarray[2][1] * inarray[2][1]
259  + inarray[2][2] * inarray[2][2] ;
260 }
261 
262 /**
263  * @brief Evaluate functional for elements connected to a node.
264  *
265  * @param minJacNew Stores current minimum Jacobian for the element group
266  * @param gradient If true, calculate gradient.
267  */
268 template <int DIM>
269 NekDouble NodeOpti::GetFunctional(NekDouble &minJacNew, bool gradient)
270 {
271  map<LibUtilities::ShapeType, vector<ElUtilSharedPtr> >::iterator typeIt;
272 
273  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
274  {
275  const int ptsStd = m_derivUtils[typeIt->first]->ptsStd;
276  const int pts = m_derivUtils[typeIt->first]->pts;
277  const int nElmt = typeIt->second.size();
278 
279  NekDouble* X = &m_tmpStore[0];
280 
281  // Store x/y components of each element sequentially in memory
282  for (int i = 0, cnt = 0; i < nElmt; ++i)
283  {
284  for (int j = 0; j < ptsStd; ++j)
285  {
286  for (int d = 0; d < DIM; ++d)
287  {
288  X[cnt + d * ptsStd + j] =
289  *(typeIt->second[i]->nodes[j][d]);
290  }
291  }
292  cnt += DIM * ptsStd;
293  }
294 
295  // Storage for derivatives, ordered by:
296  // - standard coordinate direction
297  // - number of elements
298  // - cartesian coordinate direction
299  // - quadrature points
300 
301  // Calculate x- and y-gradients
302  for (int d = 0; d < DIM; ++d)
303  {
304  Blas::Dgemm('N', 'N', pts, DIM * nElmt, ptsStd, 1.0,
305  m_derivUtils[typeIt->first]->VdmD[d].GetRawPtr(),
306  pts, X, ptsStd, 0.0,
307  &m_derivs[typeIt->first][d][0][0][0], pts);
308  }
309  }
310 
311  minJacNew = std::numeric_limits<double>::max();
312  NekDouble integral = 0.0;
313  NekDouble ep =
314  m_minJac < 0.0 ? sqrt(1e-8 + 0.04 * m_minJac * m_minJac) : 1e-4;
315  NekDouble jacIdeal[DIM][DIM], jacDet;
316  m_grad = Array<OneD, NekDouble>(DIM == 2 ? 5 : 9, 0.0);
317 
318  switch (m_opti)
319  {
320  case eLinEl:
321  {
322  const NekDouble nu = 0.4;
323  const NekDouble mu = 1.0 / 2.0 / (1.0 + nu);
324  const NekDouble K = 1.0 / 3.0 / (1.0 - 2.0 * nu);
325 
326  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
327  {
328  const int nElmt = typeIt->second.size();
329  const int pts = m_derivUtils[typeIt->first]->pts;
330 
331  NekVector<NekDouble> &quadW =
332  m_derivUtils[typeIt->first]->quadW;
333 
334  for (int i = 0; i < nElmt; ++i)
335  {
336  for (int k = 0; k < pts; ++k)
337  {
338  NekDouble phiM[DIM][DIM];
339  for (int l = 0; l < DIM; ++l)
340  {
341  for (int n = 0; n < DIM; ++n)
342  {
343  phiM[n][l] =
344  m_derivs[typeIt->first][l][i][n][k];
345  }
346  }
347 
348  // begin CalcIdealJac
349  for (int m = 0; m < DIM; ++m)
350  {
351  for (int n = 0; n < DIM; ++n)
352  {
353  jacIdeal[n][m] = 0.0;
354  for (int l = 0; l < DIM; ++l)
355  {
356  jacIdeal[n][m] += phiM[n][l] *
357  typeIt->second[i]->maps[k][m * 3 + l];
358  }
359  }
360  }
361  jacDet = Determinant(jacIdeal);
362  // end CalcIdealJac
363 
364  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
365  minJacNew = min(minJacNew, jacDet);
366 
367  NekDouble Emat[DIM][DIM];
368  EMatrix<DIM>(jacIdeal, Emat);
369 
370  NekDouble trEtE = FrobProd<DIM>(Emat, Emat);
371  NekDouble sigma =
372  0.5 *
373  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
374 
375  if (sigma < numeric_limits<double>::min() && !gradient)
376  {
377  return numeric_limits<double>::max();
378  }
379  ASSERTL0(sigma > numeric_limits<double>::min(),
380  std::string("dividing by zero ") +
381  boost::lexical_cast<string>(sigma) + " " +
382  boost::lexical_cast<string>(jacDet) + " " +
383  boost::lexical_cast<string>(ep));
384 
385  NekDouble lsigma = log(sigma);
386  integral += quadW[k] *
387  absIdealMapDet *
388  (K * 0.5 * lsigma * lsigma + mu * trEtE);
389 
390  if (gradient)
391  {
392  NekDouble jacDerivPhi[DIM];
393  NekDouble jacDetDeriv[DIM];
394 
395  NekDouble derivDet = Determinant<DIM>(phiM);
396  NekDouble jacInvTrans[DIM][DIM];
397  InvTrans<DIM>(phiM, jacInvTrans);
398 
399  NekDouble basisDeriv[DIM];
400  for (int m = 0; m < DIM; ++m)
401  {
402  basisDeriv[m] = *(
403  m_derivUtils[typeIt->first]->VdmD[m])(
404  k, typeIt->second[i]->NodeId(m_node->m_id));
405  }
406  // jacDeriv is actually a tensor,
407  // but can be stored as a vector, as 18 out of 27 entries are zero
408  // and the other 9 entries are three triplets
409  // this is due to the delta function in jacDeriv
410  NekDouble jacDeriv[DIM];
411  for (int l = 0; l < DIM; ++l)
412  {
413  jacDeriv[l] = basisDeriv[l];
414  }
415 
416  // jacDerivPhi is actually a tensor,
417  // but can be stored as a vector due to the simple form of jacDeriv
418  for (int n = 0; n < DIM; ++n)
419  {
420  jacDerivPhi[n] = 0.0;
421  for (int l = 0; l < DIM; ++l)
422  {
423  jacDerivPhi[n] += jacDeriv[l] *
424  typeIt->second[i]->maps[k][l + 3 * n];
425  }
426  }
427 
428  for (int m = 0; m < DIM; ++m)
429  {
430  jacDetDeriv[m] = 0.0;
431  for (int n = 0; n < DIM; ++n)
432  {
433  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
434  }
435  jacDetDeriv[m] *= derivDet / absIdealMapDet;
436  }
437  // end of common part to all four versionsNekDouble
438 
439 
440  NekDouble M2[DIM][DIM][DIM];
441  // use the delta function in jacDeriv and do some tensor calculus
442  // to come up with this simplified expression for:
443  // LEM2<DIM>(jacIdeal, jacDerivPhi, M2);
444  for(int d = 0; d < DIM; d++)
445  {
446  for (int m = 0; m < DIM; ++m)
447  {
448  for (int n = 0; n < DIM; ++n)
449  {
450  M2[d][m][n] = 0.5*(jacDerivPhi[m] * jacIdeal[d][n]
451  + jacIdeal[d][m] * jacDerivPhi[n]);
452  }
453  }
454  }
455 
456  for (int m = 0; m < DIM; ++m)
457  {
458  NekDouble frobProdA = FrobProd<DIM>(M2[m], Emat);
459 
460  m_grad[m] +=
461  quadW[k] * absIdealMapDet *
462  (2.0 * mu * frobProdA +
463  K * lsigma * jacDetDeriv[m] /
464  (2.0 * sigma - jacDet));
465  }
466 
467  int ct = 0;
468  for (int m = 0; m < DIM; ++m)
469  {
470  for (int l = m; l < DIM; ++l, ct++)
471  {
472  NekDouble frobProdBC = FrobProd<DIM>(M2[m], M2[l]);
473  NekDouble M3[DIM][DIM];
474  // use the delta function in jacDeriv and do some tensor calculus
475  // to come up with this simplified expression for:
476  // LEM3<DIM>(jacDerivPhi, M3);
477  if (m == l)
478  {
479  for (int p = 0; p < DIM; ++p)
480  {
481  for (int q = 0; q < DIM; ++q)
482  {
483  M3[p][q] = jacDerivPhi[p] * jacDerivPhi[q];
484  }
485  }
486  frobProdBC += FrobProd<DIM>(M3,Emat);
487  }
488 
489  m_grad[ct + DIM] +=
490  quadW[k] * absIdealMapDet *
491  (2.0 * mu * frobProdBC +
492  jacDetDeriv[m] * jacDetDeriv[l] * K /
493  (2.0 * sigma - jacDet) /
494  (2.0 * sigma - jacDet) *
495  (1.0 - jacDet * lsigma /
496  (2.0 * sigma - jacDet)));
497  }
498  }
499  }
500  }
501  }
502  }
503  break;
504  }
505 
506  case eHypEl:
507  {
508  const NekDouble nu = 0.4;
509  const NekDouble mu = 1.0 / 2.0 / (1.0 + nu);
510  const NekDouble K = 1.0 / 3.0 / (1.0 - 2.0 * nu);
511 
512  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
513  {
514  const int nElmt = typeIt->second.size();
515  const int pts = m_derivUtils[typeIt->first]->pts;
516 
517  NekVector<NekDouble> &quadW =
518  m_derivUtils[typeIt->first]->quadW;
519 
520  for (int i = 0; i < nElmt; ++i)
521  {
522  for (int k = 0; k < pts; ++k)
523  {
524  NekDouble phiM[DIM][DIM];
525  for (int l = 0; l < DIM; ++l)
526  {
527  for (int n = 0; n < DIM; ++n)
528  {
529  phiM[n][l] =
530  m_derivs[typeIt->first][l][i][n][k];
531  }
532  }
533  // begin CalcIdealJac
534  for (int m = 0; m < DIM; ++m)
535  {
536  for (int n = 0; n < DIM; ++n)
537  {
538  jacIdeal[n][m] = 0.0;
539  for (int l = 0; l < DIM; ++l)
540  {
541  jacIdeal[n][m] += phiM[n][l] *
542  typeIt->second[i]->maps[k][m * 3 + l];
543  }
544  }
545  }
546  jacDet = Determinant(jacIdeal);
547  // end CalcIdealJac
548 
549  minJacNew = min(minJacNew, jacDet);
550 
551  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
552 
553  NekDouble I1 = FrobeniusNorm(jacIdeal);
554 
555  NekDouble sigma =
556  0.5 *
557  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
558 
559  if (sigma < numeric_limits<double>::min() && !gradient)
560  {
561  return numeric_limits<double>::max();
562  }
563 
564  ASSERTL0(sigma > numeric_limits<double>::min(),
565  std::string("dividing by zero ") +
566  boost::lexical_cast<string>(sigma) + " " +
567  boost::lexical_cast<string>(jacDet) + " " +
568  boost::lexical_cast<string>(ep));
569 
570  NekDouble lsigma = log(sigma);
571  integral += quadW[k] * absIdealMapDet *
572  (0.5 * mu * (I1 - 3.0 - 2.0 * lsigma) +
573  0.5 * K * lsigma * lsigma);
574 
575  // Derivative of basis function in each direction
576  if (gradient)
577  {
578  NekDouble jacDerivPhi[DIM];
579  NekDouble jacDetDeriv[DIM];
580 
581  NekDouble derivDet = Determinant<DIM>(phiM);
582  NekDouble jacInvTrans[DIM][DIM];
583  InvTrans<DIM>(phiM, jacInvTrans);
584 
585  NekDouble basisDeriv[DIM];
586  for (int m = 0; m < DIM; ++m)
587  {
588  basisDeriv[m] =
589  *(m_derivUtils[typeIt->first]->VdmD[m].GetRawPtr() +
590  typeIt->second[i]->NodeId(m_node->m_id) * pts + k);
591  }
592 
593  // jacDeriv is actually a tensor,
594  // but can be stored as a vector, as 18 out of 27 entries are zero
595  // and the other 9 entries are three triplets
596  // this is due to the delta function in jacDeriv
597  NekDouble jacDeriv[DIM];
598  for (int l = 0; l < DIM; ++l)
599  {
600  jacDeriv[l] = basisDeriv[l];
601  }
602 
603  // jacDerivPhi is actually a tensor,
604  // but can be stored as a vector due to the simple form of jacDeriv
605  for (int n = 0; n < DIM; ++n)
606  {
607  jacDerivPhi[n] = 0.0;
608  for (int l = 0; l < DIM; ++l)
609  {
610  jacDerivPhi[n] += jacDeriv[l] *
611  typeIt->second[i]->maps[k][l + 3 * n];
612  }
613  }
614 
615  for (int m = 0; m < DIM; ++m)
616  {
617  jacDetDeriv[m] = 0.0;
618  for (int n = 0; n < DIM; ++n)
619  {
620  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
621  }
622  jacDetDeriv[m] *= derivDet / absIdealMapDet;
623  }
624  // end of common part to all four versionsNekDouble
625 
626  for (int m = 0; m < DIM; ++m)
627  {
628  // because of the zero entries of the tensor jacDerivPhi,
629  // the Frobenius-product becomes a scalar product
630  NekDouble frobProd =
631  ScalarProd<DIM>(jacIdeal[m],jacDerivPhi);
632 
633  m_grad[m] +=
634  quadW[k] * absIdealMapDet *
635  (mu * frobProd +
636  (jacDetDeriv[m] / (2.0 * sigma - jacDet) *
637  (K * lsigma - mu)));
638  }
639 
640  int ct = 0;
641  for (int m = 0; m < DIM; ++m)
642  {
643  for (int l = m; l < DIM; ++l, ct++)
644  {
645  NekDouble frobProdHes = 0.0;
646  // because of the zero entries of the tensor jacDerivPhi,
647  // the matrix frobProdHes has only diagonal entries
648  if (m == l)
649  {
650  // because of the zero entries of the tensor jacDerivPhi,
651  // the Frobenius-product becomes a scalar product
652  frobProdHes = ScalarProd<DIM>(jacDerivPhi,jacDerivPhi);
653  }
654 
655  m_grad[ct + DIM] +=
656  quadW[k] * absIdealMapDet *
657  (mu * frobProdHes +
658  jacDetDeriv[m] * jacDetDeriv[l] /
659  (2.0 * sigma - jacDet) /
660  (2.0 * sigma - jacDet) *
661  (K - jacDet * (K * lsigma - mu) /
662  (2.0 * sigma - jacDet)));
663  }
664  }
665  }
666  }
667  }
668  }
669  break;
670  }
671 
672  case eRoca:
673  {
674  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
675  {
676  const int nElmt = typeIt->second.size();
677  const int pts = m_derivUtils[typeIt->first]->pts;
678 
679  NekVector<NekDouble> &quadW =
680  m_derivUtils[typeIt->first]->quadW;
681 
682  for (int i = 0; i < nElmt; ++i)
683  {
684  for (int k = 0; k < pts; ++k)
685  {
686  NekDouble phiM[DIM][DIM];
687  for (int l = 0; l < DIM; ++l)
688  {
689  for (int n = 0; n < DIM; ++n)
690  {
691  phiM[n][l] =
692  m_derivs[typeIt->first][l][i][n][k];
693  }
694  }
695  // begin CalcIdealJac
696  for (int m = 0; m < DIM; ++m)
697  {
698  for (int n = 0; n < DIM; ++n)
699  {
700  jacIdeal[n][m] = 0.0;
701  for (int l = 0; l < DIM; ++l)
702  {
703  jacIdeal[n][m] += phiM[n][l] *
704  typeIt->second[i]->maps[k][m * 3 + l];
705  }
706  }
707  }
708  jacDet = Determinant(jacIdeal);
709  // end CalcIdealJac
710 
711  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
712  minJacNew = min(minJacNew, jacDet);
713  NekDouble frob = FrobeniusNorm(jacIdeal);
714  NekDouble sigma =
715  0.5 *
716  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
717 
718  if (sigma < numeric_limits<double>::min() && !gradient)
719  {
720  return numeric_limits<double>::max();
721  }
722 
723  ASSERTL0(sigma > numeric_limits<double>::min(),
724  std::string("dividing by zero ") +
725  boost::lexical_cast<string>(sigma) + " " +
726  boost::lexical_cast<string>(jacDet) + " " +
727  boost::lexical_cast<string>(ep));
728 
729  NekDouble W = frob / DIM / pow(fabs(sigma), 2.0 / DIM);
730  integral +=
731  quadW[k] * absIdealMapDet * W;
732 
733  // Derivative of basis function in each direction
734  if (gradient)
735  {
736  NekDouble jacDerivPhi[DIM];
737  NekDouble jacDetDeriv[DIM];
738 
739  NekDouble derivDet = Determinant<DIM>(phiM);
740  NekDouble jacInvTrans[DIM][DIM];
741  InvTrans<DIM>(phiM, jacInvTrans);
742 
743  NekDouble basisDeriv[DIM];
744  for (int m = 0; m < DIM; ++m)
745  {
746  basisDeriv[m] = *(
747  m_derivUtils[typeIt->first]->VdmD[m])(
748  k, typeIt->second[i]->NodeId(m_node->m_id));
749  }
750  // jacDeriv is actually a tensor,
751  // but can be stored as a vector, as 18 out of 27 entries are zero
752  // and the other 9 entries are three triplets
753  // this is due to the delta function in jacDeriv
754  NekDouble jacDeriv[DIM];
755  for (int l = 0; l < DIM; ++l)
756  {
757  jacDeriv[l] = basisDeriv[l];
758  }
759 
760  // jacDerivPhi is actually a tensor,
761  // but can be stored as a vector due to the simple form of jacDeriv
762  for (int n = 0; n < DIM; ++n)
763  {
764  jacDerivPhi[n] = 0.0;
765  for (int l = 0; l < DIM; ++l)
766  {
767  jacDerivPhi[n] += jacDeriv[l] *
768  typeIt->second[i]->maps[k][l + 3 * n];
769  }
770  }
771 
772  for (int m = 0; m < DIM; ++m)
773  {
774  jacDetDeriv[m] = 0.0;
775  for (int n = 0; n < DIM; ++n)
776  {
777  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
778  }
779  jacDetDeriv[m] *= derivDet / absIdealMapDet;
780  }
781  // end of common part to all four versionsNekDouble
782 
783  NekDouble frobProd[DIM];
784  NekDouble inc[DIM];
785  for (int m = 0; m < DIM; ++m)
786  {
787  // because of the zero entries of the tensor jacDerivPhi,
788  // the Frobenius-product becomes a scalar product
789  frobProd[m] = ScalarProd<DIM>(jacIdeal[m],jacDerivPhi);
790 
791  inc[m] =
792  quadW[k] * absIdealMapDet *
793  (2.0 * W * (frobProd[m] / frob -
794  jacDetDeriv[m] / DIM /
795  (2.0 * sigma - jacDet)));
796  m_grad[m] += inc[m];
797  }
798 
799 
800 
801  int ct = 0;
802  for (int m = 0; m < DIM; ++m)
803  {
804  for (int l = m; l < DIM; ++l, ct++)
805  {
806  NekDouble frobProdHes = 0.0;
807  // because of the zero entries of the tensor jacDerivPhi,
808  // the matrix frobProdHes has only diagonal entries
809  if (m == l)
810  {
811  // because of the zero entries of the tensor jacDerivPhi,
812  // the Frobenius-product becomes a scalar product
813  frobProdHes = ScalarProd<DIM>(jacDerivPhi,jacDerivPhi);
814  }
815 
816  m_grad[ct + DIM] +=
817  quadW[k] * absIdealMapDet *
818  (inc[m] * inc[l] / W +
819  2.0 * W *
820  (frobProdHes / frob -
821  2.0 * frobProd[m] * frobProd[l] /
822  frob / frob +
823  jacDetDeriv[m] * jacDetDeriv[l] *
824  jacDet /
825  (2.0 * sigma - jacDet) /
826  (2.0 * sigma - jacDet) /
827  (2.0 * sigma - jacDet) /
828  DIM));
829  }
830  }
831  }
832  }
833  }
834  }
835  break;
836  }
837 
838  case eWins:
839  {
840  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
841  {
842  const int nElmt = typeIt->second.size();
843  const int pts = m_derivUtils[typeIt->first]->pts;
844 
845  NekVector<NekDouble> &quadW =
846  m_derivUtils[typeIt->first]->quadW;
847 
848  for (int i = 0; i < nElmt; ++i)
849  {
850  for (int k = 0; k < pts; ++k)
851  {
852  NekDouble phiM[DIM][DIM];
853  for (int l = 0; l < DIM; ++l)
854  {
855  for (int n = 0; n < DIM; ++n)
856  {
857  phiM[n][l] =
858  m_derivs[typeIt->first][l][i][n][k];
859  }
860  }
861  // begin CalcIdealJac
862  for (int m = 0; m < DIM; ++m)
863  {
864  for (int n = 0; n < DIM; ++n)
865  {
866  jacIdeal[n][m] = 0.0;
867  for (int l = 0; l < DIM; ++l)
868  {
869  jacIdeal[n][m] += phiM[n][l] *
870  typeIt->second[i]->maps[k][m * 3 + l];
871  }
872  }
873  }
874  jacDet = Determinant(jacIdeal);
875  // end CalcIdealJac
876 
877  NekDouble absIdealMapDet = fabs(typeIt->second[i]->maps[k][9]);
878  minJacNew = min(minJacNew, jacDet);
879  NekDouble frob = FrobeniusNorm(jacIdeal);
880  NekDouble sigma =
881  0.5 *
882  (jacDet + sqrt(jacDet * jacDet + 4.0 * ep * ep));
883 
884  if (sigma < numeric_limits<double>::min() && !gradient)
885  {
886  return numeric_limits<double>::max();
887  }
888 
889  ASSERTL0(sigma > numeric_limits<double>::min(),
890  std::string("dividing by zero ") +
891  boost::lexical_cast<string>(sigma) + " " +
892  boost::lexical_cast<string>(jacDet) + " " +
893  boost::lexical_cast<string>(ep));
894 
895  NekDouble W = frob / sigma;
896  integral +=
897  quadW[k] * absIdealMapDet * W;
898 
899  // Derivative of basis function in each direction
900  if (gradient)
901  {
902  NekDouble jacDerivPhi[DIM];
903  NekDouble jacDetDeriv[DIM];
904 
905  NekDouble derivDet = Determinant<DIM>(phiM);
906  NekDouble jacInvTrans[DIM][DIM];
907  InvTrans<DIM>(phiM, jacInvTrans);
908 
909  NekDouble basisDeriv[DIM];
910  for (int m = 0; m < DIM; ++m)
911  {
912  basisDeriv[m] = *(
913  m_derivUtils[typeIt->first]->VdmD[m])(
914  k, typeIt->second[i]->NodeId(m_node->m_id));
915  }
916  // jacDeriv is actually a tensor,
917  // but can be stored as a vector, as 18 out of 27 entries are zero
918  // and the other 9 entries are three triplets
919  // this is due to the delta function in jacDeriv
920  NekDouble jacDeriv[DIM];
921  for (int l = 0; l < DIM; ++l)
922  {
923  jacDeriv[l] = basisDeriv[l];
924  }
925 
926  // jacDerivPhi is actually a tensor,
927  // but can be stored as a vector due to the simple form of jacDeriv
928  for (int n = 0; n < DIM; ++n)
929  {
930  jacDerivPhi[n] = 0.0;
931  for (int l = 0; l < DIM; ++l)
932  {
933  jacDerivPhi[n] += jacDeriv[l] *
934  typeIt->second[i]->maps[k][l + 3 * n];
935  }
936  }
937 
938  for (int m = 0; m < DIM; ++m)
939  {
940  jacDetDeriv[m] = 0.0;
941  for (int n = 0; n < DIM; ++n)
942  {
943  jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
944  }
945  jacDetDeriv[m] *= derivDet / absIdealMapDet;
946  }
947  // end of common part to all four versionsNekDouble
948 
949  NekDouble frobProd[DIM];
950  NekDouble inc[DIM];
951  for (int m = 0; m < DIM; ++m)
952  {
953  // because of the zero entries of the tensor jacDerivPhi,
954  // the Frobenius-product becomes a scalar product
955  frobProd[m] = ScalarProd<DIM>(jacIdeal[m],jacDerivPhi);
956 
957  inc[m] =
958  quadW[k] *
959  absIdealMapDet *
960  (W *
961  (2.0 * frobProd[m] / frob -
962  jacDetDeriv[m] / (2.0 * sigma - jacDet)));
963  m_grad[m] += inc[m];
964  }
965 
966  int ct = 0;
967  for (int m = 0; m < DIM; ++m)
968  {
969  for (int l = m; l < DIM; ++l, ct++)
970  {
971  NekDouble frobProdHes = 0.0;
972  // because of the zero entries of the tensor jacDerivPhi,
973  // the matrix frobProdHes has only diagonal entries
974  if (m == l)
975  {
976  // because of the zero entries of the tensor jacDerivPhi,
977  // the Frobenius-product becomes a scalar product
978  frobProdHes = ScalarProd<DIM>(jacDerivPhi,jacDerivPhi);
979  }
980 
981  m_grad[ct + DIM] +=
982  quadW[k] * absIdealMapDet *
983  (inc[m] * inc[l] / W +
984  2.0 * W *
985  (frobProdHes / frob -
986  2.0 * frobProd[m] * frobProd[l] /
987  frob / frob +
988  0.5 * jacDetDeriv[m] *
989  jacDetDeriv[l] * jacDet /
990  (2.0 * sigma - jacDet) /
991  (2.0 * sigma - jacDet) /
992  (2.0 * sigma - jacDet)));
993  }
994  }
995  }
996  }
997  }
998  }
999  break;
1000  }
1001  }
1002 
1003  // ASSERTL0(std::isfinite(integral),"inf in integral");
1004 
1005  return integral;
1006  // return sqrt(m_grad[0]*m_grad[0] + m_grad[1]*m_grad[1]);
1007 }
1008 }
1009 }
1010 
1011 #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:249
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
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:240
Nektar::Utilities::NodeOpti::GetFunctional
NekDouble GetFunctional(NekDouble &minJacNew, bool gradient=true)
Evaluate functional for elements connected to a node.
Definition: Evaluator.hxx:269
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:234
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 Sat Nov 25 2017 16:18:11 for Nektar++ by   doxygen 1.8.8