Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::NodeOpti:
Collaboration graph
[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
 
boost::mutex mtx
 
std::map
< LibUtilities::ShapeType,
std::vector< ElUtilSharedPtr > > 		m_data
 
Array< OneD, NekDouble > m_grad
 
std::vector< NekDouble > m_tmpStore
 
boost::unordered_map
< LibUtilities::ShapeType,
DerivArray > 		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 52 of file NodeOpti.h.

Member Typedef Documentation

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

Definition at line 56 of file NodeOpti.h.

Constructor & Destructor Documentation

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 59 of file NodeOpti.h.

References Nektar::iterator, m_data, m_derivs, m_derivUtils, and m_tmpStore.

63  : m_node(n), m_res(r), m_derivUtils(d), m_opti(o)
64  {
65  // filter element types within d vector
66  for (int i = 0; i < e.size(); i++)
67  {
68  m_data[e[i]->GetEl()->GetShapeType()].push_back(e[i]);
69  }
70 
71  // Set up storage for GetFunctional to avoid reallocation on each call.
72  size_t storageCount = 0;
73  std::map<LibUtilities::ShapeType, std::vector<ElUtilSharedPtr> >
74  ::iterator typeIt;
75  // Count total storage needed.
76  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
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_data
std::map< LibUtilities::ShapeType, std::vector< ElUtilSharedPtr > > m_data
Definition: NodeOpti.h:108
Nektar::Utilities::NodeOpti::m_derivs
boost::unordered_map< LibUtilities::ShapeType, DerivArray > m_derivs
Definition: NodeOpti.h:111
Nektar::Utilities::NodeOpti::m_tmpStore
std::vector< NekDouble > m_tmpStore
Definition: NodeOpti.h:110
Nektar::Utilities::NodeOpti::DerivArray
boost::multi_array< NekDouble, 4 > DerivArray
Definition: NodeOpti.h:56
Nektar::Utilities::NodeOpti::m_res
ResidualSharedPtr m_res
Definition: NodeOpti.h:117
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Nektar::Utilities::NodeOpti::m_derivUtils
std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > m_derivUtils
Definition: NodeOpti.h:118
virtual Nektar::Utilities::NodeOpti::~NodeOpti ( )
inlinevirtual

Definition at line 93 of file NodeOpti.h.

93 {};

Member Function Documentation

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

Definition at line 129 of file NodeOpti.h.

130  {
131  return 1e-8;
132  }
static NekDouble Nektar::Utilities::NodeOpti::c1 ( )
inlinestaticprotected

Definition at line 121 of file NodeOpti.h.

122  {
123  return 1e-3;
124  }
void Nektar::Utilities::NodeOpti::CalcMinJac ( )

Definition at line 58 of file NodeOpti.cpp.

References Nektar::iterator.

59 {
60  m_minJac = numeric_limits<double>::max();
61  map<LibUtilities::ShapeType, vector<ElUtilSharedPtr> >::iterator typeIt;
62  for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
63  {
64  for (int i = 0; i < typeIt->second.size(); i++)
65  {
66  m_minJac = min(m_minJac, typeIt->second[i]->GetMinJac());
67  }
68  }
69 }
Nektar::Utilities::NodeOpti::m_data
std::map< LibUtilities::ShapeType, std::vector< ElUtilSharedPtr > > m_data
Definition: NodeOpti.h:108
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
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 269 of file Evaluator.hxx.

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

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 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::Utilities::NodeOpti::m_grad
Array< OneD, NekDouble > m_grad
Definition: NodeOpti.h:109
Nektar::Utilities::NodeOpti::m_data
std::map< LibUtilities::ShapeType, std::vector< ElUtilSharedPtr > > m_data
Definition: NodeOpti.h:108
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
Nektar::Utilities::NodeOpti::m_derivs
boost::unordered_map< LibUtilities::ShapeType, DerivArray > m_derivs
Definition: NodeOpti.h:111
Nektar::Utilities::Determinant
NekDouble Determinant(NekDouble jac[][DIM])
Calculate determinant of input matrix.
Definition: Evaluator.hxx:55
Nektar::Utilities::NodeOpti::m_tmpStore
std::vector< NekDouble > m_tmpStore
Definition: NodeOpti.h:110
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Nektar::Utilities::FrobeniusNorm
NekDouble FrobeniusNorm(NekDouble inarray[][DIM])
Calculate Frobenius norm of a matrix .
Definition: Evaluator.hxx:234
Nektar::Utilities::NodeOpti::m_derivUtils
std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > m_derivUtils
Definition: NodeOpti.h:118
NodeOptiJob * Nektar::Utilities::NodeOpti::GetJob ( )

Definition at line 276 of file NodeOpti.cpp.

277 {
278  return new NodeOptiJob(this);
279 }
static NekDouble Nektar::Utilities::NodeOpti::gradTol ( )
inlinestaticprotected

Definition at line 125 of file NodeOpti.h.

126  {
127  return 1e-8;
128  }
template<>
int Nektar::Utilities::NodeOpti::IsIndefinite ( )
protected

Definition at line 56 of file Hessian.hxx.

References ASSERTL0, and Nektar::eDIAGONAL.

57 {
58  Array<OneD, NekDouble> eigR(2);
59  Array<OneD, NekDouble> eigI(2);
60  NekMatrix<NekDouble> H(2, 2);
61  H(0, 0) = m_grad[2];
62  H(1, 0) = m_grad[3];
63  H(0, 1) = H(1, 0);
64  H(1, 1) = m_grad[4];
65 
66  // cout << H << endl << endl;
67 
68  int nVel = 2;
69  char jobvl = 'N', jobvr = 'N';
70  int worklen = 8 * nVel, info;
71  NekDouble dum;
72 
73  DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
74  Array<OneD, NekDouble> vl(nVel * nVel);
75  Array<OneD, NekDouble> work(worklen);
76  Array<OneD, NekDouble> wi(nVel);
77 
78  Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
79  &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
80 
81  ASSERTL0(!info, "dgeev failed");
82 
83  if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0)
84  {
85  if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0)
86  {
87  return 2;
88  }
89  else
90  {
91  return 1;
92  }
93  }
94 
95  return 0;
96 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::Utilities::NodeOpti::m_grad
Array< OneD, NekDouble > m_grad
Definition: NodeOpti.h:109
Nektar::DNekMat
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
template<>
int Nektar::Utilities::NodeOpti::IsIndefinite ( )
protected

Definition at line 98 of file Hessian.hxx.

References ASSERTL0, and Nektar::eDIAGONAL.

99 {
100  Array<OneD, NekDouble> eigR(3);
101  Array<OneD, NekDouble> eigI(3);
102  NekMatrix<NekDouble> H(3, 3);
103  H(0, 0) = m_grad[3];
104  H(1, 0) = m_grad[4];
105  H(0, 1) = H(1, 0);
106  H(2, 0) = m_grad[5];
107  H(0, 2) = H(2, 0);
108  H(1, 1) = m_grad[6];
109  H(2, 1) = m_grad[7];
110  H(1, 2) = H(2, 1);
111  H(2, 2) = m_grad[8];
112 
113  int nVel = 3;
114  char jobvl = 'N', jobvr = 'N';
115  int worklen = 8 * nVel, info;
116  NekDouble dum;
117 
118  DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
119  Array<OneD, NekDouble> vl(nVel * nVel);
120  Array<OneD, NekDouble> work(worklen);
121  Array<OneD, NekDouble> wi(nVel);
122 
123  Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
124  &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
125 
126  ASSERTL0(!info, "dgeev failed");
127 
128  if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0 || eval(2, 2) < 0.0)
129  {
130  if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0 && eval(2, 2))
131  {
132  return 2;
133  }
134  else
135  {
136  return 1;
137  }
138  }
139 
140  return 0;
141 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::Utilities::NodeOpti::m_grad
Array< OneD, NekDouble > m_grad
Definition: NodeOpti.h:109
Nektar::DNekMat
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
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 50 of file Hessian.hxx.

References ASSERTL0.

51 {
52  ASSERTL0(false, "DIM error");
53  return 0;
54 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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 148 of file Hessian.hxx.

References ASSERTL0.

149 {
150  ASSERTL0(false, "DIM error");
151 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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 }
Nektar::Utilities::NodeOpti::m_grad
Array< OneD, NekDouble > m_grad
Definition: NodeOpti.h:109
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
template<>
void Nektar::Utilities::NodeOpti::MinEigen ( NekDouble &  val)

Definition at line 168 of file Hessian.hxx.

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

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::NodeOpti::m_grad
Array< OneD, NekDouble > m_grad
Definition: NodeOpti.h:109
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::Utilities::Determinant< 3 >
NekDouble Determinant< 3 >(NekDouble jac[][3])
Definition: Evaluator.hxx:65
virtual void Nektar::Utilities::NodeOpti::Optimise ( )
pure virtual

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

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

Member Data Documentation

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

Definition at line 108 of file NodeOpti.h.

Referenced by NodeOpti().

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

Definition at line 111 of file NodeOpti.h.

Referenced by NodeOpti().

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

Definition at line 118 of file NodeOpti.h.

Referenced by NodeOpti().

Array<OneD, NekDouble> Nektar::Utilities::NodeOpti::m_grad
protected

Definition at line 109 of file NodeOpti.h.

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

Definition at line 116 of file NodeOpti.h.

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

Definition at line 106 of file NodeOpti.h.

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

Definition at line 119 of file NodeOpti.h.

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

Definition at line 117 of file NodeOpti.h.

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

Definition at line 110 of file NodeOpti.h.

Referenced by NodeOpti().

boost::mutex Nektar::Utilities::NodeOpti::mtx
protected

Definition at line 107 of file NodeOpti.h.

Generated on Sat Nov 25 2017 16:19:12 for Nektar++ by   doxygen 1.8.8