Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Friends | List of all members
Nektar::SpatialDomains::GeomFactors Class Reference

Calculation and storage of geometric factors associated with the mapping from StdRegions reference elements to a given LocalRegions physical element in the mesh. More...

#include <GeomFactors.h>

Collaboration diagram for Nektar::SpatialDomains::GeomFactors:
Collaboration graph
[legend]

Public Member Functions

 GeomFactors (const GeomType gtype, const int coordim, const StdRegions::StdExpansionSharedPtr &xmap, const Array< OneD, Array< OneD, NekDouble > > &coords)
 Constructor for GeomFactors class. More...
 
 GeomFactors (const GeomFactors &S)
 Copy constructor. More...
 
 ~GeomFactors ()
 Destructor. More...
 
DerivStorage GetDeriv (const LibUtilities::PointsKeyVector &keyTgt)
 Return the derivative of the mapping with respect to the reference coordinates, $\frac{\partial \chi_i}{\partial \xi_j}$. More...
 
const Array< OneD, const
NekDouble
GetJac (const LibUtilities::PointsKeyVector &keyTgt)
 Return the Jacobian of the mapping and cache the result. More...
 
const Array< TwoD, const
NekDouble
GetGmat (const LibUtilities::PointsKeyVector &keyTgt)
 Return the Laplacian coefficients $g_{ij}$. More...
 
const Array< TwoD, const
NekDouble
GetDerivFactors (const LibUtilities::PointsKeyVector &keyTgt)
 Return the derivative of the reference coordinates with respect to the mapping, $\frac{\partial \xi_i}{\partial \chi_j}$. More...
 
GeomType GetGtype ()
 Returns whether the geometry is regular or deformed. More...
 
bool IsValid () const
 Determine if element is valid and not self-intersecting. More...
 
int GetCoordim () const
 Return the number of dimensions of the coordinate system. More...
 
size_t GetHash ()
 Computes a hash of this GeomFactors element. More...
 

Protected Member Functions

StdRegions::StdExpansionSharedPtrGetXmap ()
 Return the Xmap;. More...
 

Protected Attributes

GeomType m_type
 Type of geometry (e.g. eRegular, eDeformed, eMovingRegular). More...
 
int m_expDim
 Dimension of expansion. More...
 
int m_coordDim
 Dimension of coordinate system. More...
 
bool m_valid
 Validity of element (Jacobian positive) More...
 
StdRegions::StdExpansionSharedPtr m_xmap
 Stores information about the expansion. More...
 
Array< OneD, Array< OneD,
NekDouble > > 
m_coords
 Stores coordinates of the geometry. More...
 
std::map
< LibUtilities::PointsKeyVector,
Array< OneD, NekDouble > > 
m_jacCache
 Jacobian vector cache. More...
 
std::map
< LibUtilities::PointsKeyVector,
Array< TwoD, NekDouble > > 
m_derivFactorCache
 DerivFactors vector cache. More...
 

Private Member Functions

void CheckIfValid ()
 Tests if the element is valid and not self-intersecting. More...
 
DerivStorage ComputeDeriv (const LibUtilities::PointsKeyVector &keyTgt) const
 
Array< OneD, NekDoubleComputeJac (const LibUtilities::PointsKeyVector &keyTgt) const
 Return the Jacobian of the mapping and cache the result. More...
 
Array< TwoD, NekDoubleComputeGmat (const LibUtilities::PointsKeyVector &keyTgt) const
 Computes the Laplacian coefficients $g_{ij}$. More...
 
Array< TwoD, NekDoubleComputeDerivFactors (const LibUtilities::PointsKeyVector &keyTgt) const
 Return the derivative of the reference coordinates with respect to the mapping, $\frac{\partial \xi_i}{\partial \chi_j}$. More...
 
void Interp (const LibUtilities::PointsKeyVector &src_points, const Array< OneD, const NekDouble > &src, const LibUtilities::PointsKeyVector &tgt_points, Array< OneD, NekDouble > &tgt) const
 Perform interpolation of data between two point distributions. More...
 
void Adjoint (const Array< TwoD, const NekDouble > &src, Array< TwoD, NekDouble > &tgt) const
 Compute the transpose of the cofactors matrix. More...
 

Friends

bool operator== (const GeomFactors &lhs, const GeomFactors &rhs)
 Tests if two GeomFactors classes are equal. More...
 

Detailed Description

Calculation and storage of geometric factors associated with the mapping from StdRegions reference elements to a given LocalRegions physical element in the mesh.

This class stores the various geometric factors associated with a specific element, necessary for fundamental integration and differentiation operations as well as edge and surface normals.

Initially, these algorithms are provided with a mapping from the reference region element to the physical element. Practically, this is represented using a corresponding reference region element for each coordinate component. Note that for straight-sided elements, these elements will be of linear order. Curved elements are represented using higher-order coordinate mappings. This geometric order is in contrast to the order of the spectral/hp expansion order on the element.

For application of the chain rule during differentiation we require the partial derivatives

\[\frac{\partial \xi_i}{\partial \chi_j}\]

evaluated at the physical points of the expansion basis. We also construct the inverse metric tensor $g^{ij}$ which, in the case of a domain embedded in a higher-dimensional space, supports the conversion of covariant quantities to contravariant quantities. When the expansion dimension is equal to the coordinate dimension the Jacobian of the mapping $\chi_j$ is a square matrix and consequently the required terms are the entries of the inverse of the Jacobian. However, in general this is not the case, so we therefore implement the construction of these terms following the derivation in Cantwell, et. al. [CaYaKiPeSh13.] Given the coordinate maps $\chi_i$, this comprises of five steps

  1. Compute the terms of the Jacobian $\frac{\partial \chi_i}{\partial \xi_j}$.
  2. Compute the metric tensor $g_{ij}=\mathbf{t}_{(i)}\cdot\mathbf{t}_{(j)}$.
  3. Compute the square of the Jacobian determinant $g=|\mathbf{g}|$.
  4. Compute the inverse metric tensor $g^{ij}$.
  5. Compute the terms $\frac{\partial \xi_i}{\partial \chi_j}$.

Definition at line 80 of file GeomFactors.h.

Constructor & Destructor Documentation

Nektar::SpatialDomains::GeomFactors::GeomFactors ( const GeomType  gtype,
const int  coordim,
const StdRegions::StdExpansionSharedPtr xmap,
const Array< OneD, Array< OneD, NekDouble > > &  coords 
)

Constructor for GeomFactors class.

Parameters
gtypeSpecified whether the geometry is regular or deformed.
coordimSpecifies the dimension of the coordinate system.
CoordsCoordinate maps of the element.

Definition at line 89 of file GeomFactors.cpp.

References CheckIfValid().

93  :
94  m_type(gtype),
95  m_expDim(xmap->GetShapeDimension()),
96  m_coordDim(coordim),
97  m_valid(true),
98  m_xmap(xmap),
99  m_coords(coords)
100  {
101  CheckIfValid();
102  }
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
void CheckIfValid()
Tests if the element is valid and not self-intersecting.
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
bool m_valid
Validity of element (Jacobian positive)
Definition: GeomFactors.h:139
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
Array< OneD, Array< OneD, NekDouble > > m_coords
Stores coordinates of the geometry.
Definition: GeomFactors.h:143
StdRegions::StdExpansionSharedPtr m_xmap
Stores information about the expansion.
Definition: GeomFactors.h:141
Nektar::SpatialDomains::GeomFactors::GeomFactors ( const GeomFactors S)

Copy constructor.

Parameters
SAn instance of a GeomFactors class from which to construct a new instance.

Definition at line 109 of file GeomFactors.cpp.

109  :
110  m_type(S.m_type),
111  m_expDim(S.m_expDim),
112  m_coordDim(S.m_coordDim),
113  m_valid(S.m_valid),
114  m_xmap(S.m_xmap),
115  m_coords(S.m_coords)
116  {
117  }
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
bool m_valid
Validity of element (Jacobian positive)
Definition: GeomFactors.h:139
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
Array< OneD, Array< OneD, NekDouble > > m_coords
Stores coordinates of the geometry.
Definition: GeomFactors.h:143
StdRegions::StdExpansionSharedPtr m_xmap
Stores information about the expansion.
Definition: GeomFactors.h:141
Nektar::SpatialDomains::GeomFactors::~GeomFactors ( )

Destructor.

Definition at line 123 of file GeomFactors.cpp.

124  {
125  }

Member Function Documentation

void Nektar::SpatialDomains::GeomFactors::Adjoint ( const Array< TwoD, const NekDouble > &  src,
Array< TwoD, NekDouble > &  tgt 
) const
private

Compute the transpose of the cofactors matrix.

Input and output arrays are of dimension (m_expDim*m_expDim) x num_points. The first index of the input and output arrays are ordered row-by-row.

Parameters
srcInput data array.
tgtStorage for adjoint matrix data.

Definition at line 589 of file GeomFactors.cpp.

References ASSERTL1, Vmath::Fill(), m_expDim, Vmath::Smul(), Vmath::Vcopy(), and Vmath::Vvtvvtm().

Referenced by ComputeDerivFactors(), ComputeGmat(), and ComputeJac().

592  {
593  ASSERTL1(src.num_elements() == tgt.num_elements(),
594  "Source matrix is of different size to destination"
595  "matrix for computing adjoint.");
596 
597  int n = src[0].num_elements();
598  switch (m_expDim)
599  {
600  case 1:
601  Vmath::Fill (n, 1.0, &tgt[0][0], 1);
602  break;
603  case 2:
604  Vmath::Vcopy(n, &src[3][0], 1, &tgt[0][0], 1);
605  Vmath::Smul (n, -1.0, &src[1][0], 1, &tgt[1][0], 1);
606  Vmath::Smul (n, -1.0, &src[2][0], 1, &tgt[2][0], 1);
607  Vmath::Vcopy(n, &src[0][0], 1, &tgt[3][0], 1);
608  break;
609  case 3:
610  {
611  int a, b, c, d, e, i, j;
612 
613  // Compute g^{ij} by computing Cofactors(g_ij)^T
614  for (i = 0; i < m_expDim; ++i)
615  {
616  for (j = 0; j < m_expDim; ++j)
617  {
618  a = ((i+1)%m_expDim) * m_expDim + ((j+1)%m_expDim);
619  b = ((i+1)%m_expDim) * m_expDim + ((j+2)%m_expDim);
620  c = ((i+2)%m_expDim) * m_expDim + ((j+1)%m_expDim);
621  d = ((i+2)%m_expDim) * m_expDim + ((j+2)%m_expDim);
622  e = j*m_expDim + i;
623  Vmath::Vvtvvtm(n, &src[a][0], 1, &src[d][0], 1,
624  &src[b][0], 1, &src[c][0], 1,
625  &tgt[e][0], 1);
626  }
627  }
628  break;
629  }
630  }
631  }
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
void Vvtvvtm(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtm (vector times vector minus vector times vector):
Definition: Vmath.cpp:564
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
void Nektar::SpatialDomains::GeomFactors::CheckIfValid ( )
private

Tests if the element is valid and not self-intersecting.

Constructs the Jacobian as per Spencer's book p158 and tests if negative.

Definition at line 481 of file GeomFactors.cpp.

References Nektar::SpatialDomains::eMovingRegular, Nektar::SpatialDomains::eRegular, GetDeriv(), m_coordDim, m_expDim, m_type, m_valid, m_xmap, CellMLToNektar.cellml_metadata::p, Vmath::Vmin(), Vmath::Vvtvp(), and Vmath::Vvtvvtm().

Referenced by GeomFactors().

482  {
483  // Jacobian test only makes sense when expdim = coorddim
484  // If one-dimensional then element is valid.
485  if (m_coordDim != m_expDim || m_expDim == 1)
486  {
487  m_valid = true;
488  return;
489  }
490 
492  int nqtot = 1;
493  for (int i = 0; i < m_expDim; ++i)
494  {
495  p[i] = m_xmap->GetBasis(i)->GetPointsKey();
496  nqtot *= p[i].GetNumPoints();
497  }
498  int pts = (m_type == eRegular || m_type == eMovingRegular)
499  ? 1 : nqtot;
500  Array<OneD, NekDouble> jac(pts, 0.0);
501 
502  DerivStorage deriv = GetDeriv(p);
503 
504  switch (m_expDim)
505  {
506  case 2:
507  {
508  Vmath::Vvtvvtm(pts, &deriv[0][0][0], 1, &deriv[1][1][0], 1,
509  &deriv[1][0][0], 1, &deriv[0][1][0], 1,
510  &jac[0], 1);
511  break;
512  }
513  case 3:
514  {
515  Array<OneD, NekDouble> tmp(pts, 0.0);
516 
517  Vmath::Vvtvvtm(pts, &deriv[1][1][0], 1, &deriv[2][2][0], 1,
518  &deriv[2][1][0], 1, &deriv[1][2][0], 1,
519  &tmp[0], 1);
520  Vmath::Vvtvp (pts, &deriv[0][0][0], 1, &tmp[0], 1,
521  &jac[0], 1, &jac[0], 1);
522 
523  Vmath::Vvtvvtm(pts, &deriv[2][1][0], 1, &deriv[0][2][0], 1,
524  &deriv[0][1][0], 1, &deriv[2][2][0], 1,
525  &tmp[0], 1);
526  Vmath::Vvtvp (pts, &deriv[1][0][0], 1, &tmp[0], 1,
527  &jac[0], 1, &jac[0], 1);
528 
529  Vmath::Vvtvvtm(pts, &deriv[0][1][0], 1, &deriv[1][2][0], 1,
530  &deriv[1][1][0], 1, &deriv[0][2][0], 1,
531  &tmp[0], 1);
532  Vmath::Vvtvp (pts, &deriv[2][0][0], 1, &tmp[0], 1,
533  &jac[0], 1, &jac[0], 1);
534 
535  break;
536  }
537  }
538 
539  if (Vmath::Vmin(pts, &jac[0], 1) < 0)
540  {
541  m_valid = false;
542  }
543  }
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:871
DerivStorage GetDeriv(const LibUtilities::PointsKeyVector &keyTgt)
Return the derivative of the mapping with respect to the reference coordinates, . ...
Definition: GeomFactors.h:207
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > DerivStorage
Storage type for derivative of mapping.
Definition: GeomFactors.h:75
bool m_valid
Validity of element (Jacobian positive)
Definition: GeomFactors.h:139
void Vvtvvtm(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtm (vector times vector minus vector times vector):
Definition: Vmath.cpp:564
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
Geometry is straight-sided with constant geometric factors.
StdRegions::StdExpansionSharedPtr m_xmap
Stores information about the expansion.
Definition: GeomFactors.h:141
DerivStorage Nektar::SpatialDomains::GeomFactors::ComputeDeriv ( const LibUtilities::PointsKeyVector keyTgt) const
private

Derivatives are computed at the geometry point distributions and interpolated to the target point distributions.

Parameters
tpointsTarget point distributions.
Returns
Derivative of coordinate map evaluated at target point distributions.

Definition at line 170 of file GeomFactors.cpp.

References ASSERTL1, Interp(), m_coordDim, m_coords, m_expDim, and m_xmap.

Referenced by ComputeDerivFactors(), ComputeGmat(), ComputeJac(), and GetDeriv().

172  {
173  ASSERTL1(keyTgt.size() == m_expDim,
174  "Dimension of target point distribution does not match "
175  "expansion dimension.");
176 
177  int i = 0, j = 0;
178  int nqtot_map = 1;
179  int nqtot_tbasis = 1;
183 
184  // Allocate storage and compute number of points
185  for (i = 0; i < m_expDim; ++i)
186  {
187  map_points[i] = m_xmap->GetBasis(i)->GetPointsKey();
188  nqtot_map *= map_points[i].GetNumPoints();
189  nqtot_tbasis *= keyTgt[i].GetNumPoints();
190  deriv[i] = Array<OneD, Array<OneD,NekDouble> >(m_coordDim);
191  d_map[i] = Array<OneD, Array<OneD,NekDouble> >(m_coordDim);
192  }
193 
194  // Calculate local derivatives
195  for(i = 0; i < m_coordDim; ++i)
196  {
197  Array<OneD, NekDouble> tmp(nqtot_map);
198  // Transform from coefficient space to physical space
199  m_xmap->BwdTrans(m_coords[i], tmp);
200 
201  // Allocate storage and take the derivative (calculated at the
202  // points as specified in 'Coords')
203  for (j = 0; j < m_expDim; ++j)
204  {
205  d_map[j][i] = Array<OneD,NekDouble>(nqtot_map);
206  deriv[j][i] = Array<OneD,NekDouble>(nqtot_tbasis);
207  m_xmap->StdPhysDeriv(j, tmp, d_map[j][i]);
208  }
209  }
210 
211  for (i = 0; i < m_coordDim; ++i)
212  {
213  // Interpolate the derivatives:
214  // - from the points as defined in the mapping ('Coords')
215  // - to the points we at which we want to know the metrics
216  // ('tbasis')
217  bool same = true;
218  for (j = 0; j < m_expDim; ++j)
219  {
220  same = same && (map_points[j] == keyTgt[j]);
221  }
222  if( same )
223  {
224  for (j = 0; j < m_expDim; ++j)
225  {
226  deriv[j][i] = d_map[j][i];
227  }
228  }
229  else
230  {
231  for (j = 0; j < m_expDim; ++j)
232  {
233  Interp(map_points, d_map[j][i], keyTgt, deriv[j][i]);
234  }
235  }
236  }
237 
238  return deriv;
239  }
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > DerivStorage
Storage type for derivative of mapping.
Definition: GeomFactors.h:75
void Interp(const LibUtilities::PointsKeyVector &src_points, const Array< OneD, const NekDouble > &src, const LibUtilities::PointsKeyVector &tgt_points, Array< OneD, NekDouble > &tgt) const
Perform interpolation of data between two point distributions.
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
Array< OneD, Array< OneD, NekDouble > > m_coords
Stores coordinates of the geometry.
Definition: GeomFactors.h:143
StdRegions::StdExpansionSharedPtr m_xmap
Stores information about the expansion.
Definition: GeomFactors.h:141
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Array< TwoD, NekDouble > Nektar::SpatialDomains::GeomFactors::ComputeDerivFactors ( const LibUtilities::PointsKeyVector keyTgt) const
private

Return the derivative of the reference coordinates with respect to the mapping, $\frac{\partial \xi_i}{\partial \chi_j}$.

Parameters
keyTgtTarget point distributions.
Returns
Derivative factors evaluated at the target point distributions.

Definition at line 398 of file GeomFactors.cpp.

References Adjoint(), ASSERTL1, ComputeDeriv(), Nektar::SpatialDomains::eDeformed, m_coordDim, m_expDim, m_type, Vmath::Vdiv(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by GetDerivFactors().

400  {
401  ASSERTL1(keyTgt.size() == m_expDim,
402  "Dimension of target point distribution does not match "
403  "expansion dimension.");
404 
405  int i = 0, j = 0, k = 0, l = 0;
406  int ptsTgt = 1;
407 
408  if (m_type == eDeformed)
409  {
410  // Allocate storage and compute number of points
411  for (i = 0; i < m_expDim; ++i)
412  {
413  ptsTgt *= keyTgt[i].GetNumPoints();
414  }
415  }
416 
417  // Get derivative at geometry points
418  DerivStorage deriv = ComputeDeriv(keyTgt);
419 
420  Array<TwoD, NekDouble> tmp (m_expDim*m_expDim, ptsTgt, 0.0);
421  Array<TwoD, NekDouble> gmat(m_expDim*m_expDim, ptsTgt, 0.0);
422  Array<OneD, NekDouble> jac (ptsTgt, 0.0);
423  Array<TwoD, NekDouble> factors(m_expDim*m_coordDim, ptsTgt, 0.0);
424 
425  // Compute g_{ij} as t_i \cdot t_j and store in tmp
426  for (i = 0, l = 0; i < m_expDim; ++i)
427  {
428  for (j = 0; j < m_expDim; ++j, ++l)
429  {
430  for (k = 0; k < m_coordDim; ++k)
431  {
432  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
433  &deriv[j][k][0], 1,
434  &tmp[l][0], 1,
435  &tmp[l][0], 1);
436  }
437  }
438  }
439 
440  Adjoint(tmp, gmat);
441 
442  // Compute g = det(g_{ij}) (= Jacobian squared) and store
443  // temporarily in m_jac.
444  for (i = 0; i < m_expDim; ++i)
445  {
446  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i*m_expDim][0], 1,
447  &jac[0], 1, &jac[0], 1);
448  }
449 
450  for (i = 0; i < m_expDim*m_expDim; ++i)
451  {
452  Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
453  }
454 
455  // Compute the Jacobian = sqrt(g)
456  Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
457 
458  // Compute the derivative factors
459  for (k = 0, l = 0; k < m_coordDim; ++k)
460  {
461  for (j = 0; j < m_expDim; ++j, ++l)
462  {
463  for (i = 0; i < m_expDim; ++i)
464  {
465  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
466  &gmat[m_expDim*i+j][0], 1,
467  &factors[l][0], 1,
468  &factors[l][0], 1);
469  }
470  }
471  }
472 
473  return factors;
474  }
void Adjoint(const Array< TwoD, const NekDouble > &src, Array< TwoD, NekDouble > &tgt) const
Compute the transpose of the cofactors matrix.
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:241
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > DerivStorage
Storage type for derivative of mapping.
Definition: GeomFactors.h:75
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Geometry is curved or has non-constant factors.
DerivStorage ComputeDeriv(const LibUtilities::PointsKeyVector &keyTgt) const
Array< TwoD, NekDouble > Nektar::SpatialDomains::GeomFactors::ComputeGmat ( const LibUtilities::PointsKeyVector keyTgt) const
private

Computes the Laplacian coefficients $g_{ij}$.

This routine returns a two-dimensional array of values specifying the inverse metric terms associated with the coordinate mapping of the corresponding reference region to the physical element. These terms correspond to the $g^{ij}$ terms in [CaYaKiPeSh13] and, in the case of an embedded manifold, map covariant quantities to contravariant quantities. The leading index of the array is the index of the term in the tensor numbered as

\[\left(\begin{array}{ccc} 0 & 1 & 2 \\ 1 & 3 & 4 \\ 2 & 4 & 5 \end{array}\right)\]

. The second dimension is either of size 1 in the case of elements having GeomType eRegular, or of size equal to the number of quadrature points for eDeformed elements.

See also
Wikipedia "Covariance and Contravariance of Vectors"
Returns
Two-dimensional array containing the inverse metric tensor of the coordinate mapping.

Definition at line 333 of file GeomFactors.cpp.

References Adjoint(), ASSERTL1, ComputeDeriv(), Nektar::SpatialDomains::eDeformed, m_coordDim, m_expDim, m_type, Vmath::Vdiv(), and Vmath::Vvtvp().

Referenced by GetGmat().

335  {
336  ASSERTL1(keyTgt.size() == m_expDim,
337  "Dimension of target point distribution does not match "
338  "expansion dimension.");
339 
340  int i = 0, j = 0, k = 0, l = 0;
341  int ptsTgt = 1;
342 
343  if (m_type == eDeformed)
344  {
345  // Allocate storage and compute number of points
346  for (i = 0; i < m_expDim; ++i)
347  {
348  ptsTgt *= keyTgt[i].GetNumPoints();
349  }
350  }
351 
352  // Get derivative at geometry points
353  DerivStorage deriv = ComputeDeriv(keyTgt);
354 
355  Array<TwoD, NekDouble> tmp (m_expDim*m_expDim, ptsTgt, 0.0);
356  Array<TwoD, NekDouble> gmat(m_expDim*m_expDim, ptsTgt, 0.0);
357  Array<OneD, NekDouble> jac (ptsTgt, 0.0);
358 
359  // Compute g_{ij} as t_i \cdot t_j and store in tmp
360  for (i = 0, l = 0; i < m_expDim; ++i)
361  {
362  for (j = 0; j < m_expDim; ++j, ++l)
363  {
364  for (k = 0; k < m_coordDim; ++k)
365  {
366  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
367  &deriv[j][k][0], 1,
368  &tmp[l][0], 1,
369  &tmp[l][0], 1);
370  }
371  }
372  }
373 
374  Adjoint(tmp, gmat);
375 
376  // Compute g = det(g_{ij}) (= Jacobian squared) and store
377  // temporarily in m_jac.
378  for (i = 0; i < m_expDim; ++i)
379  {
380  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i*m_expDim][0], 1,
381  &jac[0], 1, &jac[0], 1);
382  }
383 
384  for (i = 0; i < m_expDim*m_expDim; ++i)
385  {
386  Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
387  }
388 
389  return gmat;
390  }
void Adjoint(const Array< TwoD, const NekDouble > &src, Array< TwoD, NekDouble > &tgt) const
Compute the transpose of the cofactors matrix.
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:241
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > DerivStorage
Storage type for derivative of mapping.
Definition: GeomFactors.h:75
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Geometry is curved or has non-constant factors.
DerivStorage ComputeDeriv(const LibUtilities::PointsKeyVector &keyTgt) const
Array< OneD, NekDouble > Nektar::SpatialDomains::GeomFactors::ComputeJac ( const LibUtilities::PointsKeyVector keyTgt) const
private

Return the Jacobian of the mapping and cache the result.

This routine returns an array of values specifying the Jacobian of the mapping at quadrature points in the element. The array is either of size 1 in the case of elements having GeomType eRegular, or of size equal to the number of quadrature points for eDeformed elements.

Returns
Array containing the Jacobian of the coordinate mapping at the quadrature points of the element.
See also
GeomType

Definition at line 253 of file GeomFactors.cpp.

References Adjoint(), ASSERTL1, ComputeDeriv(), Nektar::SpatialDomains::eDeformed, m_coordDim, m_expDim, m_type, Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by GetJac(), and Nektar::SpatialDomains::operator==().

255  {
256  ASSERTL1(keyTgt.size() == m_expDim,
257  "Dimension of target point distribution does not match "
258  "expansion dimension.");
259 
260  int i = 0, j = 0, k = 0, l = 0;
261  int ptsTgt = 1;
262 
263  if (m_type == eDeformed)
264  {
265  // Allocate storage and compute number of points
266  for (i = 0; i < m_expDim; ++i)
267  {
268  ptsTgt *= keyTgt[i].GetNumPoints();
269  }
270  }
271 
272  // Get derivative at geometry points
273  DerivStorage deriv = ComputeDeriv(keyTgt);
274 
275  Array<TwoD, NekDouble> tmp (m_expDim*m_expDim, ptsTgt, 0.0);
276  Array<TwoD, NekDouble> gmat(m_expDim*m_expDim, ptsTgt, 0.0);
277  Array<OneD, NekDouble> jac (ptsTgt, 0.0);
278 
279  // Compute g_{ij} as t_i \cdot t_j and store in tmp
280  for (i = 0, l = 0; i < m_expDim; ++i)
281  {
282  for (j = 0; j < m_expDim; ++j, ++l)
283  {
284  for (k = 0; k < m_coordDim; ++k)
285  {
286  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
287  &deriv[j][k][0], 1,
288  &tmp[l][0], 1,
289  &tmp[l][0], 1);
290  }
291  }
292  }
293 
294  Adjoint(tmp, gmat);
295 
296  // Compute g = det(g_{ij}) (= Jacobian squared) and store
297  // temporarily in m_jac.
298  for (i = 0; i < m_expDim; ++i)
299  {
300  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i*m_expDim][0], 1,
301  &jac[0], 1, &jac[0], 1);
302  }
303 
304  // Compute the Jacobian = sqrt(g)
305  Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
306 
307  return jac;
308  }
void Adjoint(const Array< TwoD, const NekDouble > &src, Array< TwoD, NekDouble > &tgt) const
Compute the transpose of the cofactors matrix.
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > DerivStorage
Storage type for derivative of mapping.
Definition: GeomFactors.h:75
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Geometry is curved or has non-constant factors.
DerivStorage ComputeDeriv(const LibUtilities::PointsKeyVector &keyTgt) const
int Nektar::SpatialDomains::GeomFactors::GetCoordim ( ) const
inline

Return the number of dimensions of the coordinate system.

This is greater than or equal to the expansion dimension.

Returns
The dimension of the coordinate system.

Definition at line 297 of file GeomFactors.h.

References m_coordDim.

298  {
299  return m_coordDim;
300  }
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
DerivStorage Nektar::SpatialDomains::GeomFactors::GetDeriv ( const LibUtilities::PointsKeyVector tpoints)
inline

Return the derivative of the mapping with respect to the reference coordinates, $\frac{\partial \chi_i}{\partial \xi_j}$.

Parameters
tpointsTarget point distributions.
Returns
Derivative evaluated at target point distributions.
See also
GeomFactors::ComputeDeriv

Definition at line 207 of file GeomFactors.h.

References ComputeDeriv().

Referenced by CheckIfValid().

209  {
210  return ComputeDeriv(tpoints);
211  }
DerivStorage ComputeDeriv(const LibUtilities::PointsKeyVector &keyTgt) const
const Array< TwoD, const NekDouble > Nektar::SpatialDomains::GeomFactors::GetDerivFactors ( const LibUtilities::PointsKeyVector keyTgt)
inline

Return the derivative of the reference coordinates with respect to the mapping, $\frac{\partial \xi_i}{\partial \chi_j}$.

Returns cached value if available, otherwise computes derivative factors and stores result in cache.

Parameters
keyTgtTarget point distributions.
Returns
Derivative factors evaluated at target point distributions.
See also
GeomFactors::ComputeDerivFactors

Definition at line 262 of file GeomFactors.h.

References ComputeDerivFactors(), and m_derivFactorCache.

264  {
266  Array<TwoD, NekDouble> >::const_iterator x;
267 
268  if ((x = m_derivFactorCache.find(keyTgt)) != m_derivFactorCache.end())
269  {
270  return x->second;
271  }
272 
273  m_derivFactorCache[keyTgt] = ComputeDerivFactors(keyTgt);
274 
275  return m_derivFactorCache[keyTgt];
276 
277  }
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
Array< TwoD, NekDouble > ComputeDerivFactors(const LibUtilities::PointsKeyVector &keyTgt) const
Return the derivative of the reference coordinates with respect to the mapping, . ...
std::map< LibUtilities::PointsKeyVector, Array< TwoD, NekDouble > > m_derivFactorCache
DerivFactors vector cache.
Definition: GeomFactors.h:149
const Array< TwoD, const NekDouble > Nektar::SpatialDomains::GeomFactors::GetGmat ( const LibUtilities::PointsKeyVector keyTgt)
inline

Return the Laplacian coefficients $g_{ij}$.

Parameters
keyTgtTarget point distributions.
Returns
Inverse metric tensor evaluated at target point distributions.
See also
GeomFactors::ComputeGmat

Definition at line 246 of file GeomFactors.h.

References ComputeGmat().

248  {
249  return ComputeGmat(keyTgt);
250  }
Array< TwoD, NekDouble > ComputeGmat(const LibUtilities::PointsKeyVector &keyTgt) const
Computes the Laplacian coefficients .
GeomType Nektar::SpatialDomains::GeomFactors::GetGtype ( )
inline

Returns whether the geometry is regular or deformed.

A geometric shape is considered regular if it has constant geometric information, and deformed if this information changes throughout the shape.

Returns
The type of geometry.
See also
GeomType

Definition at line 287 of file GeomFactors.h.

References m_type.

288  {
289  return m_type;
290  }
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
size_t Nektar::SpatialDomains::GeomFactors::GetHash ( )
inline

Computes a hash of this GeomFactors element.

The hash is computed from the geometry type, expansion dimension, coordinate dimension and Jacobian.

Returns
Hash of this GeomFactors object.

Definition at line 317 of file GeomFactors.h.

References Nektar::SpatialDomains::eDeformed, GetJac(), m_coordDim, m_expDim, m_type, and m_xmap.

318  {
319  LibUtilities::PointsKeyVector ptsKeys = m_xmap->GetPointsKeys();
320  const Array<OneD, const NekDouble> jac = GetJac(ptsKeys);
321 
322  size_t hash = 0;
323  boost::hash_combine(hash, (int)m_type);
324  boost::hash_combine(hash, m_expDim);
325  boost::hash_combine(hash, m_coordDim);
326  if (m_type == eDeformed)
327  {
328  boost::hash_range(hash, jac.begin(), jac.end());
329  }
330  else
331  {
332  boost::hash_combine(hash, jac[0]);
333  }
334  return hash;
335  }
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Definition: GeomFactors.h:133
int m_coordDim
Dimension of coordinate system.
Definition: GeomFactors.h:137
const Array< OneD, const NekDouble > GetJac(const LibUtilities::PointsKeyVector &keyTgt)
Return the Jacobian of the mapping and cache the result.
Definition: GeomFactors.h:223
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
StdRegions::StdExpansionSharedPtr m_xmap
Stores information about the expansion.
Definition: GeomFactors.h:141
Geometry is curved or has non-constant factors.
const Array< OneD, const NekDouble > Nektar::SpatialDomains::GeomFactors::GetJac ( const LibUtilities::PointsKeyVector keyTgt)
inline

Return the Jacobian of the mapping and cache the result.

Returns cached value if available, otherwise computes Jacobian and stores result in cache.

Parameters
keyTgtTarget point distributions.
Returns
Jacobian evaluated at target point distributions.
See also
GeomFactors::ComputeJac

Definition at line 223 of file GeomFactors.h.

References ComputeJac(), and m_jacCache.

Referenced by GetHash().

225  {
227  Array<OneD, NekDouble> >::const_iterator x;
228 
229  if ((x = m_jacCache.find(keyTgt)) != m_jacCache.end())
230  {
231  return x->second;
232  }
233 
234  m_jacCache[keyTgt] = ComputeJac(keyTgt);
235 
236  return m_jacCache[keyTgt];
237 
238  }
Array< OneD, NekDouble > ComputeJac(const LibUtilities::PointsKeyVector &keyTgt) const
Return the Jacobian of the mapping and cache the result.
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
std::map< LibUtilities::PointsKeyVector, Array< OneD, NekDouble > > m_jacCache
Jacobian vector cache.
Definition: GeomFactors.h:146
StdRegions::StdExpansionSharedPtr & Nektar::SpatialDomains::GeomFactors::GetXmap ( void  )
inlineprotected

Return the Xmap;.

Definition at line 337 of file GeomFactors.h.

References m_xmap.

338  {
339  return m_xmap;
340  }
StdRegions::StdExpansionSharedPtr m_xmap
Stores information about the expansion.
Definition: GeomFactors.h:141
void Nektar::SpatialDomains::GeomFactors::Interp ( const LibUtilities::PointsKeyVector src_points,
const Array< OneD, const NekDouble > &  src,
const LibUtilities::PointsKeyVector tgt_points,
Array< OneD, NekDouble > &  tgt 
) const
private

Perform interpolation of data between two point distributions.

Parameters
map_pointsSource data point distribution.
srcSource data to be interpolated.
tpointsTarget data point distribution.
tgtTarget data storage.

Definition at line 552 of file GeomFactors.cpp.

References ASSERTL1, Nektar::LibUtilities::Interp1D(), Nektar::LibUtilities::Interp2D(), Nektar::LibUtilities::Interp3D(), and m_expDim.

Referenced by ComputeDeriv().

557  {
558  ASSERTL1(src_points.size() == tgt_points.size(),
559  "Dimension of target point distribution does not match "
560  "expansion dimension.");
561 
562  switch (m_expDim)
563  {
564  case 1:
565  LibUtilities::Interp1D(src_points[0], src,
566  tgt_points[0], tgt);
567  break;
568  case 2:
569  LibUtilities::Interp2D(src_points[0], src_points[1], src,
570  tgt_points[0], tgt_points[1], tgt);
571  break;
572  case 3:
573  LibUtilities::Interp3D(src_points[0], src_points[1],
574  src_points[2], src,
575  tgt_points[0], tgt_points[1],
576  tgt_points[2], tgt);
577  break;
578  }
579  }
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:116
int m_expDim
Dimension of expansion.
Definition: GeomFactors.h:135
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis...
Definition: Interp.cpp:186
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:54
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
bool Nektar::SpatialDomains::GeomFactors::IsValid ( ) const
inline

Determine if element is valid and not self-intersecting.

The validity test is performed by testing if the Jacobian is negative at any point in the shape.

Returns
True if the element is not self-intersecting.

Definition at line 307 of file GeomFactors.h.

References m_valid.

308  {
309  return m_valid;
310  }
bool m_valid
Validity of element (Jacobian positive)
Definition: GeomFactors.h:139

Friends And Related Function Documentation

bool operator== ( const GeomFactors lhs,
const GeomFactors rhs 
)
friend

Tests if two GeomFactors classes are equal.

Member data equivalence is tested in the following order: shape type, expansion dimension, coordinate dimension and coordinates.

Definition at line 132 of file GeomFactors.cpp.

133  {
134  if(!(lhs.m_type == rhs.m_type))
135  {
136  return false;
137  }
138 
139  if(!(lhs.m_expDim == rhs.m_expDim))
140  {
141  return false;
142  }
143 
144  if(!(lhs.m_coordDim == rhs.m_coordDim))
145  {
146  return false;
147  }
148 
149  const Array<OneD, const NekDouble> jac_lhs =
150  lhs.ComputeJac(lhs.m_xmap->GetPointsKeys());
151  const Array<OneD, const NekDouble> jac_rhs =
152  rhs.ComputeJac(rhs.m_xmap->GetPointsKeys());
153  if(!(jac_lhs == jac_rhs))
154  {
155  return false;
156  }
157 
158  return true;
159  }
StandardMatrixTag & lhs

Member Data Documentation

int Nektar::SpatialDomains::GeomFactors::m_coordDim
protected
Array<OneD, Array<OneD, NekDouble> > Nektar::SpatialDomains::GeomFactors::m_coords
protected

Stores coordinates of the geometry.

Definition at line 143 of file GeomFactors.h.

Referenced by ComputeDeriv().

std::map<LibUtilities::PointsKeyVector, Array<TwoD, NekDouble> > Nektar::SpatialDomains::GeomFactors::m_derivFactorCache
protected

DerivFactors vector cache.

Definition at line 149 of file GeomFactors.h.

Referenced by GetDerivFactors().

int Nektar::SpatialDomains::GeomFactors::m_expDim
protected
std::map<LibUtilities::PointsKeyVector, Array<OneD, NekDouble> > Nektar::SpatialDomains::GeomFactors::m_jacCache
protected

Jacobian vector cache.

Definition at line 146 of file GeomFactors.h.

Referenced by GetJac().

GeomType Nektar::SpatialDomains::GeomFactors::m_type
protected

Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).

Definition at line 133 of file GeomFactors.h.

Referenced by CheckIfValid(), ComputeDerivFactors(), ComputeGmat(), ComputeJac(), GetGtype(), GetHash(), and Nektar::SpatialDomains::operator==().

bool Nektar::SpatialDomains::GeomFactors::m_valid
protected

Validity of element (Jacobian positive)

Definition at line 139 of file GeomFactors.h.

Referenced by CheckIfValid(), and IsValid().

StdRegions::StdExpansionSharedPtr Nektar::SpatialDomains::GeomFactors::m_xmap
protected

Stores information about the expansion.

Definition at line 141 of file GeomFactors.h.

Referenced by CheckIfValid(), ComputeDeriv(), GetHash(), GetXmap(), and Nektar::SpatialDomains::operator==().