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

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

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

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

◆ ~GeomFactors()

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

Destructor.

Definition at line 108 of file GeomFactors.cpp.

109{
110}

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

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

References ASSERTL1, Nektar::UnitTests::d(), 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 577 of file GeomFactors.cpp.

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

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

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

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

References Adjoint(), ASSERTL1, ComputeDeriv(), Nektar::SpatialDomains::eDeformed, Nektar::VarcoeffHashingTest::factors, 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"] (http://en.wikipedia.org/wiki/Covariance_and_contravariance_of_vectors)
Returns
Two-dimensional array containing the inverse metric tensor of the coordinate mapping.

Definition at line 312 of file GeomFactors.cpp.

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

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

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

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

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

References ASSERTL1, ComputeDeriv(), ComputePrincipleDirection(), Nektar::SpatialDomains::eDeformed, Nektar::SpatialDomains::eLOCAL, Nektar::VarcoeffHashingTest::factors, 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 720 of file GeomFactors.cpp.

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

308{
309 return m_coordDim;
310}

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

218{
219 return ComputeDeriv(tpoints);
220}

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

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

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

254{
255 return ComputeGmat(keyTgt);
256}
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 298 of file GeomFactors.h.

299{
300 return m_type;
301}

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

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

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

233{
234 auto x = m_jacCache.find(keyTgt);
235
236 if (x != m_jacCache.end())
237 {
238 return x->second;
239 }
240
241 m_jacCache[keyTgt] = ComputeJac(keyTgt);
242
243 return m_jacCache[keyTgt];
244}
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:143

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

287{
288 ComputeMovingFrames(keyTgt, MMFdir, CircCentre, outarray);
289}
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 345 of file GeomFactors.h.

346{
347 return m_xmap;
348}

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

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

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

318{
319 return m_valid;
320}

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

929{
930 ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
931 ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
932 ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
933
934 int nq = v1[0].size();
935 Array<OneD, NekDouble> temp(nq);
936
937 Vmath::Vmul(nq, v1[2], 1, v2[1], 1, temp, 1);
938 Vmath::Vvtvm(nq, v1[1], 1, v2[2], 1, temp, 1, v3[0], 1);
939
940 Vmath::Vmul(nq, v1[0], 1, v2[2], 1, temp, 1);
941 Vmath::Vvtvm(nq, v1[2], 1, v2[0], 1, temp, 1, v3[1], 1);
942
943 Vmath::Vmul(nq, v1[1], 1, v2[0], 1, temp, 1);
944 Vmath::Vvtvm(nq, v1[0], 1, v2[1], 1, temp, 1, v3[2], 1);
945}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
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.hpp:72
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.hpp:381

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

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

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

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

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 141 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 146 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 143 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 136 of file GeomFactors.h.

◆ m_type

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

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

Definition at line 127 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 133 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 139 of file GeomFactors.h.

Referenced by CheckIfValid(), ComputeDeriv(), ComputePrincipleDirection(), GetHash(), and GetXmap().