Nektar++
Public Member Functions | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Private Types | List of all members
Nektar::Utilities::NodeOpti Class Referenceabstract

#include <NodeOpti.h>

Inheritance diagram for Nektar::Utilities::NodeOpti:
[legend]

Public Member Functions

 NodeOpti (NodeSharedPtr n, std::vector< ElUtilSharedPtr > e, ResidualSharedPtr r, std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > d, optiType o, int dim)
 
virtual ~NodeOpti ()
 
void CalcMinJac ()
 
virtual void Optimise ()=0
 
NodeOptiJob * GetJob ()
 
template<int DIM>
NekDouble GetFunctional (NekDouble &minJacNew, bool gradient=true)
 Evaluate functional for elements connected to a node. More...
 
template<int DIM>
void MinEigen (NekDouble &val)
 Calculates minimum eigenvalue of Hessian matrix. More...
 
template<>
void MinEigen (NekDouble &val)
 
template<>
void MinEigen (NekDouble &val)
 

Protected Member Functions

template<int DIM>
int IsIndefinite ()
 Returns 1 if Hessian matrix is indefinite and 0 otherwise. More...
 
template<>
int IsIndefinite ()
 
template<>
int IsIndefinite ()
 

Static Protected Member Functions

static NekDouble c1 ()
 
static NekDouble gradTol ()
 
static NekDouble alphaTol ()
 

Protected Attributes

NodeSharedPtr m_node
 
std::mutex mtx
 
std::map< LibUtilities::ShapeType, std::vector< ElUtilSharedPtr > > m_data
 
std::vector< NekDouble > m_grad
 
std::vector< NekDouble > m_tmpStore
 
std::unordered_map< LibUtilities::ShapeType, DerivArray, EnumHash > m_derivs
 
NekDouble m_minJac
 
ResidualSharedPtr m_res
 
std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > m_derivUtils
 
optiType m_opti
 

Private Types

typedef boost::multi_array< NekDouble, 4 > DerivArray
 

Detailed Description

Definition at line 53 of file NodeOpti.h.

Member Typedef Documentation

◆ DerivArray

typedef boost::multi_array<NekDouble, 4> Nektar::Utilities::NodeOpti::DerivArray
private

Definition at line 57 of file NodeOpti.h.

Constructor & Destructor Documentation

◆ NodeOpti()

Nektar::Utilities::NodeOpti::NodeOpti ( NodeSharedPtr  n,
std::vector< ElUtilSharedPtr >  e,
ResidualSharedPtr  r,
std::map< LibUtilities::ShapeType, DerivUtilSharedPtr >  d,
optiType  o,
int  dim 
)
inline

Definition at line 60 of file NodeOpti.h.

64  : m_node(n), m_res(r), m_derivUtils(d), m_opti(o)
65  {
66  // filter element types within d vector
67  for (int i = 0; i < e.size(); i++)
68  {
69  m_data[e[i]->GetEl()->GetShapeType()].push_back(e[i]);
70  }
71 
72  // Set up storage for GetFunctional to avoid reallocation on each call.
73  size_t storageCount = 0;
74 
75  // Count total storage needed.
76  for (auto &typeIt : m_data)
77  {
78  const int pts = m_derivUtils[typeIt.first]->pts;
79  const int nElmt = typeIt.second.size();
80 
81  storageCount = std::max(storageCount,
82  dim * m_derivUtils[typeIt.first]->ptsStd *
83  typeIt.second.size());
84 
85  m_derivs.insert(std::make_pair(
86  typeIt.first,
87  DerivArray(boost::extents[dim][nElmt][dim][pts])));
88  }
89 
90  m_tmpStore.resize(storageCount);
91  }
Nektar::Utilities::NodeOpti::m_derivUtils
std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > m_derivUtils
Definition: NodeOpti.h:118
Nektar::Utilities::NodeOpti::m_tmpStore
std::vector< NekDouble > m_tmpStore
Definition: NodeOpti.h:110
Nektar::Utilities::NodeOpti::m_data
std::map< LibUtilities::ShapeType, std::vector< ElUtilSharedPtr > > m_data
Definition: NodeOpti.h:108
Nektar::Utilities::NodeOpti::m_res
ResidualSharedPtr m_res
Definition: NodeOpti.h:117
Nektar::Utilities::NodeOpti::DerivArray
boost::multi_array< NekDouble, 4 > DerivArray
Definition: NodeOpti.h:57
Nektar::Utilities::NodeOpti::m_derivs
std::unordered_map< LibUtilities::ShapeType, DerivArray, EnumHash > m_derivs
Definition: NodeOpti.h:111

References m_data, m_derivs, m_derivUtils, and m_tmpStore.

◆ ~NodeOpti()

virtual Nektar::Utilities::NodeOpti::~NodeOpti ( )
inlinevirtual

Definition at line 93 of file NodeOpti.h.

93 {};

Member Function Documentation

◆ alphaTol()

static NekDouble Nektar::Utilities::NodeOpti::alphaTol ( )
inlinestaticprotected

Definition at line 129 of file NodeOpti.h.

130  {
131  return 1e-8;
132  }

◆ c1()

static NekDouble Nektar::Utilities::NodeOpti::c1 ( )
inlinestaticprotected

Definition at line 121 of file NodeOpti.h.

122  {
123  return 1e-3;
124  }

◆ CalcMinJac()

void Nektar::Utilities::NodeOpti::CalcMinJac ( )

Definition at line 57 of file NodeOpti.cpp.

58 {
59  m_minJac = numeric_limits<double>::max();
60  for (auto &typeIt : m_data)
61  {
62  for (int i = 0; i < typeIt.second.size(); i++)
63  {
64  m_minJac = min(m_minJac, typeIt.second[i]->GetMinJac());
65  }
66  }
67 }

◆ GetFunctional()

template<int DIM>
NekDouble Nektar::Utilities::NodeOpti::GetFunctional ( NekDouble &  minJacNew,
bool  gradient = true 
)

Evaluate functional for elements connected to a node.

Parameters
minJacNewStores current minimum Jacobian for the element group
gradientIf true, calculate gradient.

Definition at line 274 of file Evaluator.hxx.

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 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::Utilities::NodeOpti::m_grad
std::vector< NekDouble > m_grad
Definition: NodeOpti.h:109
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::FrobeniusNorm
NekDouble FrobeniusNorm(NekDouble inarray[][DIM])
Calculate Frobenius norm of a matrix .
Definition: Evaluator.hxx:238
Nektar::Utilities::Determinant
NekDouble Determinant(NekDouble jac[][DIM])
Calculate determinant of input matrix.
Definition: Evaluator.hxx:54
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

References ASSERTL0, Nektar::Utilities::Determinant(), Blas::Dgemm(), Nektar::Utilities::eHypEl, Nektar::Utilities::eLinEl, Nektar::Utilities::eRoca, Nektar::Utilities::eWins, Nektar::Utilities::FrobeniusNorm(), and CellMLToNektar.cellml_metadata::p.

◆ GetJob()

NodeOptiJob * Nektar::Utilities::NodeOpti::GetJob ( )

Definition at line 274 of file NodeOpti.cpp.

275 {
276  return new NodeOptiJob(this);
277 }

◆ gradTol()

static NekDouble Nektar::Utilities::NodeOpti::gradTol ( )
inlinestaticprotected

Definition at line 125 of file NodeOpti.h.

126  {
127  return 1e-8;
128  }

◆ IsIndefinite() [1/3]

template<>
int Nektar::Utilities::NodeOpti::IsIndefinite ( )
protected

Definition at line 55 of file Hessian.hxx.

56 {
57  vector<NekDouble> eigR(2);
58  vector<NekDouble> eigI(2);
59  NekMatrix<NekDouble> H(2, 2);
60  H(0, 0) = m_grad[2];
61  H(1, 0) = m_grad[3];
62  H(0, 1) = H(1, 0);
63  H(1, 1) = m_grad[4];
64 
65  // cout << H << endl << endl;
66 
67  int nVel = 2;
68  char jobvl = 'N', jobvr = 'N';
69  int worklen = 8 * nVel, info;
70  NekDouble dum;
71 
72  DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
73  vector<NekDouble> vl(nVel * nVel);
74  vector<NekDouble> work(worklen);
75  vector<NekDouble> wi(nVel);
76 
77  Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
78  &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
79 
80  ASSERTL0(!info, "dgeev failed");
81 
82  if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0)
83  {
84  if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0)
85  {
86  return 2;
87  }
88  else
89  {
90  return 1;
91  }
92  }
93 
94  return 0;
95 }
Lapack::Dgeev
static void Dgeev(const char &uplo, const char &lrev, const int &n, const double *a, const int &lda, double *wr, double *wi, double *rev, const int &ldr, double *lev, const int &ldv, double *work, const int &lwork, int &info)
Solve general real matrix eigenproblem.
Definition: Lapack.hpp:211
Nektar::DNekMat
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:51

References ASSERTL0, Lapack::Dgeev(), Nektar::eDIAGONAL, and m_grad.

◆ IsIndefinite() [2/3]

template<>
int Nektar::Utilities::NodeOpti::IsIndefinite ( )
protected

Definition at line 97 of file Hessian.hxx.

98 {
99  vector<NekDouble> eigR(3);
100  vector<NekDouble> eigI(3);
101  NekMatrix<NekDouble> H(3, 3);
102  H(0, 0) = m_grad[3];
103  H(1, 0) = m_grad[4];
104  H(0, 1) = H(1, 0);
105  H(2, 0) = m_grad[5];
106  H(0, 2) = H(2, 0);
107  H(1, 1) = m_grad[6];
108  H(2, 1) = m_grad[7];
109  H(1, 2) = H(2, 1);
110  H(2, 2) = m_grad[8];
111 
112  int nVel = 3;
113  char jobvl = 'N', jobvr = 'N';
114  int worklen = 8 * nVel, info;
115  NekDouble dum;
116 
117  DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
118  vector<NekDouble> vl(nVel * nVel);
119  vector<NekDouble> work(worklen);
120  vector<NekDouble> wi(nVel);
121 
122  Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
123  &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
124 
125  ASSERTL0(!info, "dgeev failed");
126 
127  if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0 || eval(2, 2) < 0.0)
128  {
129  if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0 && eval(2, 2))
130  {
131  return 2;
132  }
133  else
134  {
135  return 1;
136  }
137  }
138 
139  return 0;
140 }

References ASSERTL0, Lapack::Dgeev(), Nektar::eDIAGONAL, and m_grad.

◆ IsIndefinite() [3/3]

template<int DIM>
int Nektar::Utilities::NodeOpti::IsIndefinite
protected

Returns 1 if Hessian matrix is indefinite and 0 otherwise.

Specialised versions of this function exist only for 2x2 and 3x3 matrices.

Definition at line 49 of file Hessian.hxx.

50 {
51  ASSERTL0(false, "DIM error");
52  return 0;
53 }

References ASSERTL0.

◆ MinEigen() [1/3]

template<>
void Nektar::Utilities::NodeOpti::MinEigen ( NekDouble &  val)

Definition at line 153 of file Hessian.hxx.

154 {
155  NekDouble H[2][2];
156  H[0][0] = m_grad[2];
157  H[1][0] = m_grad[3];
158  //H[0][1] = H[1][0];
159  H[1][1] = m_grad[4];
160 
161  NekDouble D = (H[0][0] - H[1][1]) * (H[0][0] - H[1][1]) + 4.0 * H[1][0] * H[1][0];
162  NekDouble Dsqrt = sqrt(D);
163 
164  //eval[0] = (H[0][0] + H[1][1] + Dsqrt ) / 2.0;
165  val = (H[0][0] + H[1][1] - Dsqrt ) / 2.0; // the minimum Eigenvalue
166 }

References m_grad.

◆ MinEigen() [2/3]

template<>
void Nektar::Utilities::NodeOpti::MinEigen ( NekDouble &  val)

Definition at line 168 of file Hessian.hxx.

169 {
170  NekDouble H[3][3];
171  H[0][0] = m_grad[3];
172  H[1][0] = m_grad[4];
173  H[0][1] = H[1][0];
174  H[2][0] = m_grad[5];
175  H[0][2] = H[2][0];
176  H[1][1] = m_grad[6];
177  H[2][1] = m_grad[7];
178  H[1][2] = H[2][1];
179  H[2][2] = m_grad[8];
180 
181  //double eval[3]; // the eigenvalues
182 
183  NekDouble p1 = H[0][1] * H[0][1] + H[0][2] * H[0][2] + H[1][2] * H[1][2];
184  if (p1 == 0.0) // H is diagonal
185  {
186  // find the minimum Eigenvalue
187  if(H[0][0] < H[1][1])
188  {
189  if(H[0][0] < H[2][2])
190  {
191  val = H[0][0];
192  }
193  else
194  {
195  val = H[2][2];
196  }
197  }
198  else
199  {
200  if(H[1][1] < H[2][2])
201  {
202  val = H[1][1];
203  }
204  else
205  {
206  val = H[2][2];
207  }
208  }
209  }
210  else
211  {
212  NekDouble q = (H[0][0] + H[1][1] + H[2][2]) / 3.0;
213  NekDouble p2 = (H[0][0] - q)*(H[0][0] - q)
214  + (H[1][1] - q)*(H[1][1] - q)
215  + (H[2][2] - q)*(H[2][2] - q)
216  + 2.0 * p1;
217  NekDouble p = sqrt(p2 / 6.0);
218 
219  NekDouble B[3][3]; // B = (1.0 / p) * (H - q * I) with I being the identity matrix
220  NekDouble pinv = 1.0 / p;
221  B[0][0] = pinv * (H[0][0] - q);
222  B[1][1] = pinv * (H[1][1] - q);
223  B[2][2] = pinv * (H[2][2] - q);
224  B[0][1] = pinv * H[0][1];
225  B[1][0] = B[0][1];
226  B[0][2] = pinv * H[0][2];
227  B[2][0] = B[0][2];
228  B[1][2] = pinv * H[1][2];
229  B[2][1] = B[1][2];
230 
231  NekDouble r = Determinant<3>(B) / 2.0;
232 
233  // In exact arithmetic for h symmetric matrix -1 <= r <= 1
234  // but computation error can leave it slightly outside this range.
235  NekDouble phi;
236  if (r <= -1)
237  {
238  phi = M_PI / 3.0;
239  }
240  else if (r >= 1)
241  {
242  phi = 0.0;
243  }
244  else
245  {
246  phi = acos(r) / 3.0;
247  }
248 
249  // the eigenvalues satisfy eval[2] <= eval[1] <= eval[0]
250  //eval[0] = q + 2.0 * p * cos(phi);
251  val = q + 2.0 * p * cos(phi + (2.0*M_PI/3.0));
252  //eval[1] = 3.0 * q - eval[0] - eval[2]; // since trace(H) = eval[0] + eval[1] + eval[2]
253  }
254 }
Nektar::Utilities::Determinant< 3 >
NekDouble Determinant< 3 >(NekDouble jac[][3])
Definition: Evaluator.hxx:65

References Nektar::Utilities::Determinant< 3 >(), m_grad, and CellMLToNektar.cellml_metadata::p.

◆ MinEigen() [3/3]

template<int DIM>
void Nektar::Utilities::NodeOpti::MinEigen ( NekDouble &  val)

Calculates minimum eigenvalue of Hessian matrix.

Specialised versions of this function exist only for 2x2 and 3x3 matrices.

Definition at line 147 of file Hessian.hxx.

148 {
149  boost::ignore_unused(val);
150  ASSERTL0(false, "DIM error");
151 }

References ASSERTL0.

◆ Optimise()

virtual void Nektar::Utilities::NodeOpti::Optimise ( )
pure virtual

Implemented in Nektar::Utilities::NodeOpti2D3D, Nektar::Utilities::NodeOpti1D3D, Nektar::Utilities::NodeOpti2D2D, and Nektar::Utilities::NodeOpti3D3D.

Referenced by Nektar::Utilities::NodeOptiJob::Run().

Member Data Documentation

◆ m_data

std::map<LibUtilities::ShapeType, std::vector<ElUtilSharedPtr> > Nektar::Utilities::NodeOpti::m_data
protected

Definition at line 108 of file NodeOpti.h.

Referenced by NodeOpti().

◆ m_derivs

std::unordered_map<LibUtilities::ShapeType, DerivArray, EnumHash> Nektar::Utilities::NodeOpti::m_derivs
protected

Definition at line 111 of file NodeOpti.h.

Referenced by NodeOpti().

◆ m_derivUtils

std::map<LibUtilities::ShapeType, DerivUtilSharedPtr> Nektar::Utilities::NodeOpti::m_derivUtils
protected

Definition at line 118 of file NodeOpti.h.

Referenced by NodeOpti().

◆ m_grad

std::vector<NekDouble> Nektar::Utilities::NodeOpti::m_grad
protected

Definition at line 109 of file NodeOpti.h.

Referenced by IsIndefinite(), and MinEigen().

◆ m_minJac

NekDouble Nektar::Utilities::NodeOpti::m_minJac
protected

Definition at line 116 of file NodeOpti.h.

◆ m_node

NodeSharedPtr Nektar::Utilities::NodeOpti::m_node
protected

Definition at line 106 of file NodeOpti.h.

◆ m_opti

optiType Nektar::Utilities::NodeOpti::m_opti
protected

Definition at line 119 of file NodeOpti.h.

◆ m_res

ResidualSharedPtr Nektar::Utilities::NodeOpti::m_res
protected

Definition at line 117 of file NodeOpti.h.

◆ m_tmpStore

std::vector<NekDouble> Nektar::Utilities::NodeOpti::m_tmpStore
protected

Definition at line 110 of file NodeOpti.h.

Referenced by NodeOpti().

◆ mtx

std::mutex Nektar::Utilities::NodeOpti::mtx
protected

Definition at line 107 of file NodeOpti.h.

Generated on Wed Aug 18 2021 14:16:20 for Nektar++ by doxygen 1.9.1