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 | 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)
 
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
 
NekDouble m_minJac
 
ResidualSharedPtr m_res
 
std::map
< LibUtilities::ShapeType,
DerivUtilSharedPtr > 		m_derivUtils
 
optiType m_opti
 

Detailed Description

Definition at line 52 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 
)
inline

Definition at line 55 of file NodeOpti.h.

References m_data.

59  : m_node(n), m_res(r), m_derivUtils(d), m_opti(o)
60  {
61  // filter element types within d vector
62  for (int i = 0; i < e.size(); i++)
63  {
64  m_data[e[i]->GetEl()->GetShapeType()].push_back(e[i]);
65  }
66  }
Nektar::Utilities::NodeOpti::m_data
std::map< LibUtilities::ShapeType, std::vector< ElUtilSharedPtr > > m_data
Definition: NodeOpti.h:83
Nektar::Utilities::NodeOpti::m_res
ResidualSharedPtr m_res
Definition: NodeOpti.h:89
Nektar::Utilities::NodeOpti::m_node
NodeSharedPtr m_node
Definition: NodeOpti.h:81
Nektar::Utilities::NodeOpti::m_derivUtils
std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > m_derivUtils
Definition: NodeOpti.h:90
virtual Nektar::Utilities::NodeOpti::~NodeOpti ( )
inlinevirtual

Definition at line 68 of file NodeOpti.h.

68 {};

Member Function Documentation

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

Definition at line 101 of file NodeOpti.h.

102  {
103  return 1e-8;
104  }
static NekDouble Nektar::Utilities::NodeOpti::c1 ( )
inlinestaticprotected

Definition at line 93 of file NodeOpti.h.

94  {
95  return 1e-3;
96  }
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:83
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 274 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.

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

98  {
99  return 1e-8;
100  }
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:84
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 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:84
Nektar::DNekMat
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
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:84
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:84
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 83 of file NodeOpti.h.

Referenced by NodeOpti().

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

Definition at line 90 of file NodeOpti.h.

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

Definition at line 84 of file NodeOpti.h.

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

Definition at line 88 of file NodeOpti.h.

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

Definition at line 81 of file NodeOpti.h.

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

Definition at line 91 of file NodeOpti.h.

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

Definition at line 89 of file NodeOpti.h.

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

Definition at line 82 of file NodeOpti.h.

Generated on Wed Mar 8 2017 12:30:09 for Nektar++ by   doxygen 1.8.8