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{
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
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.cpp:716
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:245
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:43
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191

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 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:236
@ 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:1045
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:569

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.
StdRegions::ConstFactorMap factors
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:529
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:280

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 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 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.cpp:513

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 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 }
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 }
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, 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 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;
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:49
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:164
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:103

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:207
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:593

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