Nektar++
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>

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 NekDoubleGetJac (const LibUtilities::PointsKeyVector &keyTgt)
 Return the Jacobian of the mapping and cache the result. More...
 
const Array< TwoD, const NekDoubleGetGmat (const LibUtilities::PointsKeyVector &keyTgt)
 Return the Laplacian coefficients \(g_{ij}\). More...
 
const Array< TwoD, const NekDoubleGetDerivFactors (const LibUtilities::PointsKeyVector &keyTgt)
 Return the derivative of the reference coordinates with respect to the mapping, \(\frac{\partial \xi_i}{\partial \chi_j}\). More...
 
void GetMovingFrames (const LibUtilities::PointsKeyVector &keyTgt, const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble >> &outarray)
 Returns moving frames. 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...
 
enum GeomMMF m_MMFDir
 Principle tangent direction for MMF. 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 ComputeMovingFrames (const LibUtilities::PointsKeyVector &keyTgt, const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble >> &movingframes)
 
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...
 
void ComputePrincipleDirection (const LibUtilities::PointsKeyVector &keyTgt, const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble >> &output)
 
void VectorNormalise (Array< OneD, Array< OneD, NekDouble >> &array)
 
void VectorCrossProd (const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &v3)
 Computes the vector cross-product in 3D of v1 and v2, storing the result in v3. 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 73 of file GeomFactors.h.

Constructor & Destructor Documentation

◆ GeomFactors() [1/2]

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 88 of file GeomFactors.cpp.

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

References CheckIfValid().

◆ GeomFactors() [2/2]

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 101 of file GeomFactors.cpp.

102  : m_type(S.m_type), m_expDim(S.m_expDim), m_coordDim(S.m_coordDim),
103  m_valid(S.m_valid), m_xmap(S.m_xmap), m_coords(S.m_coords)
104 {
105 }

◆ ~GeomFactors()

Nektar::SpatialDomains::GeomFactors::~GeomFactors ( )

Destructor.

Definition at line 110 of file GeomFactors.cpp.

111 {
112 }

Member Function Documentation

◆ Adjoint()

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 677 of file GeomFactors.cpp.

679 {
680  ASSERTL1(src.size() == tgt.size(),
681  "Source matrix is of different size to destination"
682  "matrix for computing adjoint.");
683 
684  int n = src[0].size();
685  switch (m_expDim)
686  {
687  case 1:
688  Vmath::Fill(n, 1.0, &tgt[0][0], 1);
689  break;
690  case 2:
691  Vmath::Vcopy(n, &src[3][0], 1, &tgt[0][0], 1);
692  Vmath::Smul(n, -1.0, &src[1][0], 1, &tgt[1][0], 1);
693  Vmath::Smul(n, -1.0, &src[2][0], 1, &tgt[2][0], 1);
694  Vmath::Vcopy(n, &src[0][0], 1, &tgt[3][0], 1);
695  break;
696  case 3:
697  {
698  int a, b, c, d, e, i, j;
699 
700  // Compute g^{ij} by computing Cofactors(g_ij)^T
701  for (i = 0; i < m_expDim; ++i)
702  {
703  for (j = 0; j < m_expDim; ++j)
704  {
705  a = ((i + 1) % m_expDim) * m_expDim + ((j + 1) % m_expDim);
706  b = ((i + 1) % m_expDim) * m_expDim + ((j + 2) % m_expDim);
707  c = ((i + 2) % m_expDim) * m_expDim + ((j + 1) % m_expDim);
708  d = ((i + 2) % m_expDim) * m_expDim + ((j + 2) % m_expDim);
709  e = j * m_expDim + i;
710  Vmath::Vvtvvtm(n, &src[a][0], 1, &src[d][0], 1, &src[b][0],
711  1, &src[c][0], 1, &tgt[e][0], 1);
712  }
713  }
714  break;
715  }
716  }
717 }
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
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:721
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255

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

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

◆ CheckIfValid()

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 579 of file GeomFactors.cpp.

580 {
581  // Jacobian test only makes sense when expdim = coorddim
582  // If one-dimensional then element is valid.
583  if (m_coordDim != m_expDim || m_expDim == 1)
584  {
585  m_valid = true;
586  return;
587  }
588 
590  int nqtot = 1;
591  for (int i = 0; i < m_expDim; ++i)
592  {
593  p[i] = m_xmap->GetBasis(i)->GetPointsKey();
594  nqtot *= p[i].GetNumPoints();
595  }
596  int pts = (m_type == eRegular || m_type == eMovingRegular) ? 1 : nqtot;
597  Array<OneD, NekDouble> jac(pts, 0.0);
598 
599  DerivStorage deriv = GetDeriv(p);
600 
601  switch (m_expDim)
602  {
603  case 2:
604  {
605  Vmath::Vvtvvtm(pts, &deriv[0][0][0], 1, &deriv[1][1][0], 1,
606  &deriv[1][0][0], 1, &deriv[0][1][0], 1, &jac[0], 1);
607  break;
608  }
609  case 3:
610  {
611  Array<OneD, NekDouble> tmp(pts, 0.0);
612 
613  Vmath::Vvtvvtm(pts, &deriv[1][1][0], 1, &deriv[2][2][0], 1,
614  &deriv[2][1][0], 1, &deriv[1][2][0], 1, &tmp[0], 1);
615  Vmath::Vvtvp(pts, &deriv[0][0][0], 1, &tmp[0], 1, &jac[0], 1,
616  &jac[0], 1);
617 
618  Vmath::Vvtvvtm(pts, &deriv[2][1][0], 1, &deriv[0][2][0], 1,
619  &deriv[0][1][0], 1, &deriv[2][2][0], 1, &tmp[0], 1);
620  Vmath::Vvtvp(pts, &deriv[1][0][0], 1, &tmp[0], 1, &jac[0], 1,
621  &jac[0], 1);
622 
623  Vmath::Vvtvvtm(pts, &deriv[0][1][0], 1, &deriv[1][2][0], 1,
624  &deriv[1][1][0], 1, &deriv[0][2][0], 1, &tmp[0], 1);
625  Vmath::Vvtvp(pts, &deriv[2][0][0], 1, &tmp[0], 1, &jac[0], 1,
626  &jac[0], 1);
627 
628  break;
629  }
630  }
631 
632  if (Vmath::Vmin(pts, &jac[0], 1) < 0)
633  {
634  m_valid = false;
635  }
636 }
DerivStorage GetDeriv(const LibUtilities::PointsKeyVector &keyTgt)
Return the derivative of the mapping with respect to the reference coordinates, .
Definition: GeomFactors.h:218
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eMovingRegular
Currently unused.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > DerivStorage
Storage type for derivative of mapping.
Definition: GeomFactors.h:68
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:1050
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:574

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().

◆ ComputeDeriv()

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 155 of file GeomFactors.cpp.

157 {
158  ASSERTL1(keyTgt.size() == m_expDim,
159  "Dimension of target point distribution does not match "
160  "expansion dimension.");
161 
162  int i = 0, j = 0;
163  int nqtot_map = 1;
164  int nqtot_tbasis = 1;
168 
169  // Allocate storage and compute number of points
170  for (i = 0; i < m_expDim; ++i)
171  {
172  map_points[i] = m_xmap->GetBasis(i)->GetPointsKey();
173  nqtot_map *= map_points[i].GetNumPoints();
174  nqtot_tbasis *= keyTgt[i].GetNumPoints();
175  deriv[i] = Array<OneD, Array<OneD, NekDouble>>(m_coordDim);
176  d_map[i] = Array<OneD, Array<OneD, NekDouble>>(m_coordDim);
177  }
178 
179  // Calculate local derivatives
180  for (i = 0; i < m_coordDim; ++i)
181  {
182  Array<OneD, NekDouble> tmp(nqtot_map);
183  // Transform from coefficient space to physical space
184  m_xmap->BwdTrans(m_coords[i], tmp);
185 
186  // Allocate storage and take the derivative (calculated at the
187  // points as specified in 'Coords')
188  for (j = 0; j < m_expDim; ++j)
189  {
190  d_map[j][i] = Array<OneD, NekDouble>(nqtot_map);
191  deriv[j][i] = Array<OneD, NekDouble>(nqtot_tbasis);
192  m_xmap->StdPhysDeriv(j, tmp, d_map[j][i]);
193  }
194  }
195 
196  for (i = 0; i < m_coordDim; ++i)
197  {
198  // Interpolate the derivatives:
199  // - from the points as defined in the mapping ('Coords')
200  // - to the points at which we want to know the metrics
201  // ('tbasis')
202  bool same = true;
203  for (j = 0; j < m_expDim; ++j)
204  {
205  same = same && (map_points[j] == keyTgt[j]);
206  }
207  if (same)
208  {
209  for (j = 0; j < m_expDim; ++j)
210  {
211  deriv[j][i] = d_map[j][i];
212  }
213  }
214  else
215  {
216  for (j = 0; j < m_expDim; ++j)
217  {
218  Interp(map_points, d_map[j][i], keyTgt, deriv[j][i]);
219  }
220  }
221  }
222 
223  return deriv;
224 }
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.

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

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

◆ ComputeDerivFactors()

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. A 1D example: /f$ Jac =(\partial x/ \partial \xi) /f$ ; /f$ factor = 1/Jac = (\partial \xi/ \partial x) /f$

Definition at line 378 of file GeomFactors.cpp.

380 {
381  ASSERTL1(keyTgt.size() == m_expDim,
382  "Dimension of target point distribution does not match "
383  "expansion dimension.");
384 
385  int i = 0, j = 0, k = 0, l = 0;
386  int ptsTgt = 1;
387 
388  if (m_type == eDeformed)
389  {
390  // Allocate storage and compute number of points
391  for (i = 0; i < m_expDim; ++i)
392  {
393  ptsTgt *= keyTgt[i].GetNumPoints();
394  }
395  }
396 
397  // Get derivative at geometry points
398  DerivStorage deriv = ComputeDeriv(keyTgt);
399 
400  Array<TwoD, NekDouble> tmp(m_expDim * m_expDim, ptsTgt, 0.0);
401  Array<TwoD, NekDouble> gmat(m_expDim * m_expDim, ptsTgt, 0.0);
402  Array<OneD, NekDouble> jac(ptsTgt, 0.0);
403  Array<TwoD, NekDouble> factors(m_expDim * m_coordDim, ptsTgt, 0.0);
404 
405  // Compute g_{ij} as t_i \cdot t_j and store in tmp
406  for (i = 0, l = 0; i < m_expDim; ++i)
407  {
408  for (j = 0; j < m_expDim; ++j, ++l)
409  {
410  for (k = 0; k < m_coordDim; ++k)
411  {
412  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1, &deriv[j][k][0], 1,
413  &tmp[l][0], 1, &tmp[l][0], 1);
414  }
415  }
416  }
417 
418  Adjoint(tmp, gmat);
419 
420  // Compute g = det(g_{ij}) (= Jacobian squared) and store
421  // temporarily in m_jac.
422  for (i = 0; i < m_expDim; ++i)
423  {
424  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i * m_expDim][0], 1, &jac[0],
425  1, &jac[0], 1);
426  }
427 
428  for (i = 0; i < m_expDim * m_expDim; ++i)
429  {
430  Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
431  }
432 
433  // Compute the Jacobian = sqrt(g)
434  Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
435 
436  // Compute the derivative factors
437  for (k = 0, l = 0; k < m_coordDim; ++k)
438  {
439  for (j = 0; j < m_expDim; ++j, ++l)
440  {
441  for (i = 0; i < m_expDim; ++i)
442  {
443  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
444  &gmat[m_expDim * i + j][0], 1, &factors[l][0], 1,
445  &factors[l][0], 1);
446  }
447  }
448  }
449 
450  return factors;
451 }
DerivStorage ComputeDeriv(const LibUtilities::PointsKeyVector &keyTgt) const
void Adjoint(const Array< TwoD, const NekDouble > &src, Array< TwoD, NekDouble > &tgt) const
Compute the transpose of the cofactors matrix.
@ eDeformed
Geometry is curved or has non-constant factors.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
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:284

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

Referenced by GetDerivFactors().

◆ ComputeGmat()

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 314 of file GeomFactors.cpp.

316 {
317  ASSERTL1(keyTgt.size() == m_expDim,
318  "Dimension of target point distribution does not match "
319  "expansion dimension.");
320 
321  int i = 0, j = 0, k = 0, l = 0;
322  int ptsTgt = 1;
323 
324  if (m_type == eDeformed)
325  {
326  // Allocate storage and compute number of points
327  for (i = 0; i < m_expDim; ++i)
328  {
329  ptsTgt *= keyTgt[i].GetNumPoints();
330  }
331  }
332 
333  // Get derivative at geometry points
334  DerivStorage deriv = ComputeDeriv(keyTgt);
335 
336  Array<TwoD, NekDouble> tmp(m_expDim * m_expDim, ptsTgt, 0.0);
337  Array<TwoD, NekDouble> gmat(m_expDim * m_expDim, ptsTgt, 0.0);
338  Array<OneD, NekDouble> jac(ptsTgt, 0.0);
339 
340  // Compute g_{ij} as t_i \cdot t_j and store in tmp
341  for (i = 0, l = 0; i < m_expDim; ++i)
342  {
343  for (j = 0; j < m_expDim; ++j, ++l)
344  {
345  for (k = 0; k < m_coordDim; ++k)
346  {
347  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1, &deriv[j][k][0], 1,
348  &tmp[l][0], 1, &tmp[l][0], 1);
349  }
350  }
351  }
352 
353  Adjoint(tmp, gmat);
354 
355  // Compute g = det(g_{ij}) (= Jacobian squared) and store
356  // temporarily in m_jac.
357  for (i = 0; i < m_expDim; ++i)
358  {
359  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i * m_expDim][0], 1, &jac[0],
360  1, &jac[0], 1);
361  }
362 
363  for (i = 0; i < m_expDim * m_expDim; ++i)
364  {
365  Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
366  }
367 
368  return gmat;
369 }

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

Referenced by GetGmat().

◆ ComputeJac()

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 237 of file GeomFactors.cpp.

239 {
240  ASSERTL1(keyTgt.size() == m_expDim,
241  "Dimension of target point distribution does not match "
242  "expansion dimension.");
243 
244  int i = 0, j = 0, k = 0, l = 0;
245  int ptsTgt = 1;
246 
247  if (m_type == eDeformed)
248  {
249  // Allocate storage and compute number of points
250  for (i = 0; i < m_expDim; ++i)
251  {
252  ptsTgt *= keyTgt[i].GetNumPoints();
253  }
254  }
255 
256  // Get derivative at geometry points
257  DerivStorage deriv = ComputeDeriv(keyTgt);
258 
259  Array<TwoD, NekDouble> tmp(m_expDim * m_expDim, ptsTgt, 0.0);
260  Array<TwoD, NekDouble> gmat(m_expDim * m_expDim, ptsTgt, 0.0);
261  Array<OneD, NekDouble> jac(ptsTgt, 0.0);
262 
263  // Compute g_{ij} as t_i \cdot t_j and store in tmp
264  for (i = 0, l = 0; i < m_expDim; ++i)
265  {
266  for (j = 0; j < m_expDim; ++j, ++l)
267  {
268  for (k = 0; k < m_coordDim; ++k)
269  {
270  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1, &deriv[j][k][0], 1,
271  &tmp[l][0], 1, &tmp[l][0], 1);
272  }
273  }
274  }
275 
276  Adjoint(tmp, gmat);
277 
278  // Compute g = det(g_{ij}) (= Jacobian squared) and store
279  // temporarily in m_jac.
280  for (i = 0; i < m_expDim; ++i)
281  {
282  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i * m_expDim][0], 1, &jac[0],
283  1, &jac[0], 1);
284  }
285 
286  // Compute the Jacobian = sqrt(g)
287  Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
288 
289  return jac;
290 }

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

Referenced by GetJac().

◆ ComputeMovingFrames()

void Nektar::SpatialDomains::GeomFactors::ComputeMovingFrames ( const LibUtilities::PointsKeyVector keyTgt,
const SpatialDomains::GeomMMF  MMFdir,
const Array< OneD, const NekDouble > &  CircCentre,
Array< OneD, Array< OneD, NekDouble >> &  movingframes 
)
private

Definition at line 453 of file GeomFactors.cpp.

458 {
459  ASSERTL1(keyTgt.size() == m_expDim,
460  "Dimension of target point distribution does not match "
461  "expansion dimension.");
462 
463  int i = 0, k = 0;
464  int ptsTgt = 1;
465  int nq = 1;
466 
467  for (i = 0; i < m_expDim; ++i)
468  {
469  nq *= keyTgt[i].GetNumPoints();
470  }
471 
472  if (m_type == eDeformed)
473  {
474  // Allocate storage and compute number of points
475  for (i = 0; i < m_expDim; ++i)
476  {
477  ptsTgt *= keyTgt[i].GetNumPoints();
478  }
479  }
480 
481  // Get derivative at geometry points
482  DerivStorage deriv = ComputeDeriv(keyTgt);
483 
484  // number of moving frames is requited to be 3, even for surfaces
485  int MFdim = 3;
486 
487  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtmp(MFdim);
488 
489  // Compute g_{ij} as t_i \cdot t_j and store in tmp
490  for (i = 0; i < MFdim; ++i)
491  {
492  MFtmp[i] = Array<OneD, Array<OneD, NekDouble>>(m_coordDim);
493  for (k = 0; k < m_coordDim; ++k)
494  {
495  MFtmp[i][k] = Array<OneD, NekDouble>(nq);
496  }
497  }
498 
499  // Compute g_{ij} as t_i \cdot t_j and store in tmp
500  for (i = 0; i < MFdim - 1; ++i)
501  {
502  for (k = 0; k < m_coordDim; ++k)
503  {
504  if (m_type == eDeformed)
505  {
506  Vmath::Vcopy(ptsTgt, &deriv[i][k][0], 1, &MFtmp[i][k][0], 1);
507  }
508  else
509  {
510  Vmath::Fill(nq, deriv[i][k][0], MFtmp[i][k], 1);
511  }
512  }
513  }
514 
515  // Direction of MF1 is preserved: MF2 is considered in the same
516  // tangent plane as MF1. MF3 is computed by cross product of MF1
517  // and MF2. MF2 is consequently computed as the cross product of
518  // MF3 and MF1.
519  Array<OneD, Array<OneD, NekDouble>> PrincipleDir(m_coordDim);
520  for (int k = 0; k < m_coordDim; k++)
521  {
522  PrincipleDir[k] = Array<OneD, NekDouble>(nq);
523  }
524 
525  if (!(MMFdir == eLOCAL))
526  {
527  ComputePrincipleDirection(keyTgt, MMFdir, factors, PrincipleDir);
528  }
529 
530  // MF3 = MF1 \times MF2
531  VectorCrossProd(MFtmp[0], MFtmp[1], MFtmp[2]);
532 
533  // Normalizing MF3
534  VectorNormalise(MFtmp[2]);
535 
536  if (!(MMFdir == eLOCAL))
537  {
538  Array<OneD, NekDouble> temp(nq, 0.0);
539 
540  // Reorient MF1 along the PrincipleDir
541  for (i = 0; i < m_coordDim; ++i)
542  {
543  Vmath::Vvtvp(nq, MFtmp[2][i], 1, PrincipleDir[i], 1, temp, 1, temp,
544  1);
545  }
546  Vmath::Neg(nq, temp, 1);
547 
548  // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
549  for (i = 0; i < m_coordDim; ++i)
550  {
551  Vmath::Vvtvp(nq, temp, 1, MFtmp[2][i], 1, PrincipleDir[i], 1,
552  MFtmp[0][i], 1);
553  }
554  }
555 
556  // Normalizing MF1
557  VectorNormalise(MFtmp[0]);
558 
559  // MF2 = MF3 \times MF1
560  VectorCrossProd(MFtmp[2], MFtmp[0], MFtmp[1]);
561 
562  // Normalizing MF2
563  VectorNormalise(MFtmp[1]);
564 
565  for (i = 0; i < MFdim; ++i)
566  {
567  for (k = 0; k < m_coordDim; ++k)
568  {
569  Vmath::Vcopy(nq, &MFtmp[i][k][0], 1,
570  &movingframes[i * m_coordDim + k][0], 1);
571  }
572  }
573 }
void VectorCrossProd(const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &v3)
Computes the vector cross-product in 3D of v1 and v2, storing the result in v3.
void VectorNormalise(Array< OneD, Array< OneD, NekDouble >> &array)
void ComputePrincipleDirection(const LibUtilities::PointsKeyVector &keyTgt, const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble >> &output)
@ eLOCAL
No Principal direction.
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518

References ASSERTL1, ComputeDeriv(), ComputePrincipleDirection(), Nektar::SpatialDomains::eDeformed, Nektar::SpatialDomains::eLOCAL, Vmath::Fill(), m_coordDim, m_expDim, m_type, Vmath::Neg(), Vmath::Vcopy(), VectorCrossProd(), VectorNormalise(), and Vmath::Vvtvp().

Referenced by GetMovingFrames().

◆ ComputePrincipleDirection()

void Nektar::SpatialDomains::GeomFactors::ComputePrincipleDirection ( const LibUtilities::PointsKeyVector keyTgt,
const SpatialDomains::GeomMMF  MMFdir,
const Array< OneD, const NekDouble > &  CircCentre,
Array< OneD, Array< OneD, NekDouble >> &  output 
)
private

Definition at line 722 of file GeomFactors.cpp.

727 {
728  int nq = output[0].size();
729 
730  output = Array<OneD, Array<OneD, NekDouble>>(m_coordDim);
731  for (int i = 0; i < m_coordDim; ++i)
732  {
733  output[i] = Array<OneD, NekDouble>(nq, 0.0);
734  }
735 
736  // Construction of Connection
737  switch (MMFdir)
738  {
739  // projection to x-axis
740  case eTangentX:
741  {
742  Vmath::Fill(nq, 1.0, output[0], 1);
743  break;
744  }
745  case eTangentY:
746  {
747  Vmath::Fill(nq, 1.0, output[1], 1);
748  break;
749  }
750  case eTangentXY:
751  {
752  Vmath::Fill(nq, sqrt(2.0), output[0], 1);
753  Vmath::Fill(nq, sqrt(2.0), output[1], 1);
754  break;
755  }
756  case eTangentZ:
757  {
758  Vmath::Fill(nq, 1.0, output[2], 1);
759  break;
760  }
761  case eTangentCircular:
762  {
763  // Tangent direction depends on spatial location.
764  Array<OneD, Array<OneD, NekDouble>> x(m_coordDim);
765  for (int k = 0; k < m_coordDim; k++)
766  {
767  x[k] = Array<OneD, NekDouble>(nq);
768  }
769 
770  // m_coords are StdExpansions which store the mapping
771  // between the std element and the local element. Bwd
772  // transforming the std element minimum basis gives a
773  // minimum physical basis for geometry. Need to then
774  // interpolate this up to the quadrature basis.
775  int nqtot_map = 1;
777  for (int i = 0; i < m_expDim; ++i)
778  {
779  map_points[i] = m_xmap->GetBasis(i)->GetPointsKey();
780  nqtot_map *= map_points[i].GetNumPoints();
781  }
782  Array<OneD, NekDouble> tmp(nqtot_map);
783  for (int k = 0; k < m_coordDim; k++)
784  {
785  m_xmap->BwdTrans(m_coords[k], tmp);
786  Interp(map_points, tmp, keyTgt, x[k]);
787  }
788 
789  // circular around the center of the domain
790  NekDouble radius, xc = 0.0, yc = 0.0, xdis, ydis;
791  NekDouble la, lb;
792 
793  ASSERTL1(factors.size() >= 4, "factors is too short.");
794 
795  la = factors[0];
796  lb = factors[1];
797  xc = factors[2];
798  yc = factors[3];
799 
800  for (int i = 0; i < nq; i++)
801  {
802  xdis = x[0][i] - xc;
803  ydis = x[1][i] - yc;
804  radius = sqrt(xdis * xdis / la / la + ydis * ydis / lb / lb);
805  output[0][i] = ydis / radius;
806  output[1][i] = -1.0 * xdis / radius;
807  }
808  break;
809  }
810  case eTangentIrregular:
811  {
812  // Tangent direction depends on spatial location.
813  Array<OneD, Array<OneD, NekDouble>> x(m_coordDim);
814  for (int k = 0; k < m_coordDim; k++)
815  {
816  x[k] = Array<OneD, NekDouble>(nq);
817  }
818 
819  int nqtot_map = 1;
821  for (int i = 0; i < m_expDim; ++i)
822  {
823  map_points[i] = m_xmap->GetBasis(i)->GetPointsKey();
824  nqtot_map *= map_points[i].GetNumPoints();
825  }
826  Array<OneD, NekDouble> tmp(nqtot_map);
827  for (int k = 0; k < m_coordDim; k++)
828  {
829  m_xmap->BwdTrans(m_coords[k], tmp);
830  Interp(map_points, tmp, keyTgt, x[k]);
831  }
832 
833  // circular around the center of the domain
834  NekDouble xtan, ytan, mag;
835  for (int i = 0; i < nq; i++)
836  {
837  xtan = -1.0 * (x[1][i] * x[1][i] * x[1][i] + x[1][i]);
838  ytan = 2.0 * x[0][i];
839  mag = sqrt(xtan * xtan + ytan * ytan);
840  output[0][i] = xtan / mag;
841  output[1][i] = ytan / mag;
842  }
843  break;
844  }
845  case eTangentNonconvex:
846  {
847  // Tangent direction depends on spatial location.
848  Array<OneD, Array<OneD, NekDouble>> x(m_coordDim);
849  for (int k = 0; k < m_coordDim; k++)
850  {
851  x[k] = Array<OneD, NekDouble>(nq);
852  }
853 
854  int nqtot_map = 1;
856  for (int i = 0; i < m_expDim; ++i)
857  {
858  map_points[i] = m_xmap->GetBasis(i)->GetPointsKey();
859  nqtot_map *= map_points[i].GetNumPoints();
860  }
861  Array<OneD, NekDouble> tmp(nqtot_map);
862  for (int k = 0; k < m_coordDim; k++)
863  {
864  m_xmap->BwdTrans(m_coords[k], tmp);
865  Interp(map_points, tmp, keyTgt, x[k]);
866  }
867 
868  // circular around the center of the domain
869  NekDouble xtan, ytan, mag;
870  for (int i = 0; i < nq; i++)
871  {
872  xtan = -2.0 * x[1][i] * x[1][i] * x[1][i] + x[1][i];
873  ytan = sqrt(3.0) * x[0][i];
874  mag = sqrt(xtan * xtan + ytan * ytan);
875  output[0][i] = xtan / mag;
876  output[1][i] = ytan / mag;
877  }
878  break;
879  }
880  default:
881  {
882  break;
883  }
884  }
885 }
@ eTangentIrregular
Circular around the centre of domain.
@ eTangentX
X coordinate direction.
@ eTangentCircular
Circular around the centre of domain.
@ eTangentNonconvex
Circular around the centre of domain.
@ eTangentXY
XY direction.
@ eTangentZ
Z coordinate direction.
@ eTangentY
Y coordinate direction.
double NekDouble
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References ASSERTL1, Nektar::SpatialDomains::eTangentCircular, Nektar::SpatialDomains::eTangentIrregular, Nektar::SpatialDomains::eTangentNonconvex, Nektar::SpatialDomains::eTangentX, Nektar::SpatialDomains::eTangentXY, Nektar::SpatialDomains::eTangentY, Nektar::SpatialDomains::eTangentZ, Vmath::Fill(), Interp(), m_coordDim, m_coords, m_expDim, m_xmap, and tinysimd::sqrt().

Referenced by ComputeMovingFrames().

◆ GetCoordim()

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 309 of file GeomFactors.h.

310 {
311  return m_coordDim;
312 }

References m_coordDim.

◆ GetDeriv()

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 218 of file GeomFactors.h.

220 {
221  return ComputeDeriv(tpoints);
222 }

References ComputeDeriv().

Referenced by CheckIfValid().

◆ GetDerivFactors()

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 269 of file GeomFactors.h.

271 {
272  auto x = m_derivFactorCache.find(keyTgt);
273 
274  if (x != m_derivFactorCache.end())
275  {
276  return x->second;
277  }
278 
279  m_derivFactorCache[keyTgt] = ComputeDerivFactors(keyTgt);
280 
281  return m_derivFactorCache[keyTgt];
282 }
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:148

References ComputeDerivFactors(), and m_derivFactorCache.

◆ GetGmat()

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 254 of file GeomFactors.h.

256 {
257  return ComputeGmat(keyTgt);
258 }
Array< TwoD, NekDouble > ComputeGmat(const LibUtilities::PointsKeyVector &keyTgt) const
Computes the Laplacian coefficients .

References ComputeGmat().

◆ GetGtype()

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 300 of file GeomFactors.h.

301 {
302  return m_type;
303 }

References m_type.

◆ GetHash()

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 329 of file GeomFactors.h.

330 {
331  LibUtilities::PointsKeyVector ptsKeys = m_xmap->GetPointsKeys();
332  const Array<OneD, const NekDouble> jac = GetJac(ptsKeys);
333 
334  size_t hash = 0;
335  hash_combine(hash, (int)m_type, m_expDim, m_coordDim);
336  if (m_type == eDeformed)
337  {
338  hash_range(hash, jac.begin(), jac.end());
339  }
340  else
341  {
342  hash_combine(hash, jac[0]);
343  }
344  return hash;
345 }
const Array< OneD, const NekDouble > GetJac(const LibUtilities::PointsKeyVector &keyTgt)
Return the Jacobian of the mapping and cache the result.
Definition: GeomFactors.h:233
std::size_t hash_range(Iter first, Iter last)
Definition: HashUtils.hpp:68
void hash_combine(std::size_t &seed)
Definition: HashUtils.hpp:46

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

◆ GetJac()

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 233 of file GeomFactors.h.

235 {
236  auto x = m_jacCache.find(keyTgt);
237 
238  if (x != m_jacCache.end())
239  {
240  return x->second;
241  }
242 
243  m_jacCache[keyTgt] = ComputeJac(keyTgt);
244 
245  return m_jacCache[keyTgt];
246 }
Array< OneD, NekDouble > ComputeJac(const LibUtilities::PointsKeyVector &keyTgt) const
Return the Jacobian of the mapping and cache the result.
std::map< LibUtilities::PointsKeyVector, Array< OneD, NekDouble > > m_jacCache
Jacobian vector cache.
Definition: GeomFactors.h:145

References ComputeJac(), and m_jacCache.

Referenced by GetHash().

◆ GetMovingFrames()

void Nektar::SpatialDomains::GeomFactors::GetMovingFrames ( const LibUtilities::PointsKeyVector keyTgt,
const SpatialDomains::GeomMMF  MMFdir,
const Array< OneD, const NekDouble > &  CircCentre,
Array< OneD, Array< OneD, NekDouble >> &  outarray 
)
inline

Returns moving frames.

Definition at line 284 of file GeomFactors.h.

289 {
290  ComputeMovingFrames(keyTgt, MMFdir, CircCentre, outarray);
291 }
void ComputeMovingFrames(const LibUtilities::PointsKeyVector &keyTgt, const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble >> &movingframes)

References ComputeMovingFrames().

◆ GetXmap()

StdRegions::StdExpansionSharedPtr & Nektar::SpatialDomains::GeomFactors::GetXmap ( void  )
inlineprotected

Return the Xmap;.

Definition at line 347 of file GeomFactors.h.

348 {
349  return m_xmap;
350 }

References m_xmap.

◆ Interp()

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 644 of file GeomFactors.cpp.

648 {
649  ASSERTL1(src_points.size() == tgt_points.size(),
650  "Dimension of target point distribution does not match "
651  "expansion dimension.");
652 
653  switch (m_expDim)
654  {
655  case 1:
656  LibUtilities::Interp1D(src_points[0], src, tgt_points[0], tgt);
657  break;
658  case 2:
659  LibUtilities::Interp2D(src_points[0], src_points[1], src,
660  tgt_points[0], tgt_points[1], tgt);
661  break;
662  case 3:
663  LibUtilities::Interp3D(src_points[0], src_points[1], src_points[2],
664  src, tgt_points[0], tgt_points[1],
665  tgt_points[2], tgt);
666  break;
667  }
668 }
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:52
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:167
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:106

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

Referenced by ComputeDeriv(), and ComputePrincipleDirection().

◆ IsValid()

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 319 of file GeomFactors.h.

320 {
321  return m_valid;
322 }

References m_valid.

◆ VectorCrossProd()

void Nektar::SpatialDomains::GeomFactors::VectorCrossProd ( const Array< OneD, const Array< OneD, NekDouble >> &  v1,
const Array< OneD, const Array< OneD, NekDouble >> &  v2,
Array< OneD, Array< OneD, NekDouble >> &  v3 
)
private

Computes the vector cross-product in 3D of v1 and v2, storing the result in v3.

Parameters
v1First input vector.
v2Second input vector.
v3Output vector computed to be orthogonal to both v1 and v2.

Definition at line 927 of file GeomFactors.cpp.

931 {
932  ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
933  ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
934  ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
935 
936  int nq = v1[0].size();
937  Array<OneD, NekDouble> temp(nq);
938 
939  Vmath::Vmul(nq, v1[2], 1, v2[1], 1, temp, 1);
940  Vmath::Vvtvm(nq, v1[1], 1, v2[2], 1, temp, 1, v3[0], 1);
941 
942  Vmath::Vmul(nq, v1[0], 1, v2[2], 1, temp, 1);
943  Vmath::Vvtvm(nq, v1[2], 1, v2[0], 1, temp, 1, v3[1], 1);
944 
945  Vmath::Vmul(nq, v1[1], 1, v2[0], 1, temp, 1);
946  Vmath::Vvtvm(nq, v1[0], 1, v2[1], 1, temp, 1, v3[2], 1);
947 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
void Vmul(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:209
void Vvtvm(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)
vvtvm (vector times vector minus vector): z = w*x - y
Definition: Vmath.cpp:598

References ASSERTL0, Vmath::Vmul(), and Vmath::Vvtvm().

Referenced by ComputeMovingFrames().

◆ VectorNormalise()

void Nektar::SpatialDomains::GeomFactors::VectorNormalise ( Array< OneD, Array< OneD, NekDouble >> &  array)
private

Definition at line 890 of file GeomFactors.cpp.

891 {
892  int ndim = array.size();
893  ASSERTL0(ndim > 0, "Number of components must be > 0.");
894  for (int i = 1; i < ndim; ++i)
895  {
896  ASSERTL0(array[i].size() == array[0].size(),
897  "Array size mismatch in coordinates.");
898  }
899 
900  int nq = array[0].size();
901  Array<OneD, NekDouble> norm(nq, 0.0);
902 
903  // Compute the norm of each vector.
904  for (int i = 0; i < ndim; ++i)
905  {
906  Vmath::Vvtvp(nq, array[i], 1, array[i], 1, norm, 1, norm, 1);
907  }
908 
909  Vmath::Vsqrt(nq, norm, 1, norm, 1);
910 
911  // Normalise the vectors by the norm
912  for (int i = 0; i < ndim; ++i)
913  {
914  Vmath::Vdiv(nq, array[i], 1, norm, 1, array[i], 1);
915  }
916 }

References ASSERTL0, Vmath::Vdiv(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by ComputeMovingFrames().

Friends And Related Function Documentation

◆ operator==

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 118 of file GeomFactors.cpp.

119 {
120  if (!(lhs.m_type == rhs.m_type))
121  {
122  return false;
123  }
124 
125  if (!(lhs.m_expDim == rhs.m_expDim))
126  {
127  return false;
128  }
129 
130  if (!(lhs.m_coordDim == rhs.m_coordDim))
131  {
132  return false;
133  }
134 
135  const Array<OneD, const NekDouble> jac_lhs =
136  lhs.ComputeJac(lhs.m_xmap->GetPointsKeys());
137  const Array<OneD, const NekDouble> jac_rhs =
138  rhs.ComputeJac(rhs.m_xmap->GetPointsKeys());
139  if (!(jac_lhs == jac_rhs))
140  {
141  return false;
142  }
143 
144  return true;
145 }

Member Data Documentation

◆ m_coordDim

int Nektar::SpatialDomains::GeomFactors::m_coordDim
protected

◆ m_coords

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(), and ComputePrincipleDirection().

◆ m_derivFactorCache

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

DerivFactors vector cache.

Definition at line 148 of file GeomFactors.h.

Referenced by GetDerivFactors().

◆ m_expDim

int Nektar::SpatialDomains::GeomFactors::m_expDim
protected

◆ m_jacCache

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

Jacobian vector cache.

Definition at line 145 of file GeomFactors.h.

Referenced by GetJac().

◆ m_MMFDir

enum GeomMMF Nektar::SpatialDomains::GeomFactors::m_MMFDir
protected

Principle tangent direction for MMF.

Definition at line 135 of file GeomFactors.h.

◆ m_type

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

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

Definition at line 129 of file GeomFactors.h.

Referenced by CheckIfValid(), ComputeDerivFactors(), ComputeGmat(), ComputeJac(), ComputeMovingFrames(), GetGtype(), and GetHash().

◆ m_valid

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

Validity of element (Jacobian positive)

Definition at line 135 of file GeomFactors.h.

Referenced by CheckIfValid(), and IsValid().

◆ m_xmap

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(), ComputePrincipleDirection(), GetHash(), and GetXmap().