Nektar++
Loading...
Searching...
No Matches
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.
 
 GeomFactors (const GeomFactors &S)
 Copy constructor.
 
 ~GeomFactors ()=default
 Destructor.
 
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}\).
 
const Array< OneD, const NekDoubleGetJac (const LibUtilities::PointsKeyVector &keyTgt)
 Return the Jacobian of the mapping and cache the result.
 
const Array< TwoD, const NekDoubleGetGmat (const LibUtilities::PointsKeyVector &keyTgt)
 Return the Laplacian coefficients \(g_{ij}\).
 
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}\).
 
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.
 
GeomType GetGtype ()
 Returns whether the geometry is regular or deformed.
 
bool IsValid () const
 Determine if element is valid and not self-intersecting.
 
int GetCoordim () const
 Return the number of dimensions of the coordinate system.
 
size_t GetHash ()
 Computes a hash of this GeomFactors element.
 

Protected Member Functions

StdRegions::StdExpansionSharedPtrGetXmap ()
 Return the Xmap;.
 

Protected Attributes

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

Private Member Functions

void CheckIfValid ()
 Tests if the element is valid and not self-intersecting.
 
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.
 
Array< TwoD, NekDoubleComputeGmat (const LibUtilities::PointsKeyVector &keyTgt) const
 Computes the Laplacian coefficients \(g_{ij}\).
 
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}\).
 
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.
 
void Adjoint (const Array< TwoD, const NekDouble > &src, Array< TwoD, NekDouble > &tgt) const
 Compute the transpose of the cofactors matrix.
 
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.
 

Friends

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

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

98 : m_type(gtype), m_expDim(xmap->GetShapeDimension()), m_coordDim(coordim),
99 m_valid(true), m_xmap(xmap), m_coords(coords)
100{
101 CheckIfValid();
102}
int m_coordDim
Dimension of coordinate system.
void CheckIfValid()
Tests if the element is valid and not self-intersecting.
int m_expDim
Dimension of expansion.
bool m_valid
Validity of element (Jacobian positive)
StdRegions::StdExpansionSharedPtr m_xmap
Stores information about the expansion.
GeomType m_type
Type of geometry (e.g. eRegular, eDeformed, eMovingRegular).
Array< OneD, Array< OneD, NekDouble > > m_coords
Stores coordinates of the geometry.

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

109 : m_type(S.m_type), m_expDim(S.m_expDim), m_coordDim(S.m_coordDim),
110 m_valid(S.m_valid), m_xmap(S.m_xmap), m_coords(S.m_coords)
111{
112}

◆ ~GeomFactors()

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

Destructor.

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

698{
699 ASSERTL1(src.size() == tgt.size(),
700 "Source matrix is of different size to destination"
701 "matrix for computing adjoint.");
702
703 int n = src[0].size();
704 switch (m_expDim)
705 {
706 case 1:
707 Vmath::Fill(n, 1.0, &tgt[0][0], 1);
708 break;
709 case 2:
710 Vmath::Vcopy(n, &src[3][0], 1, &tgt[0][0], 1);
711 Vmath::Smul(n, -1.0, &src[1][0], 1, &tgt[1][0], 1);
712 Vmath::Smul(n, -1.0, &src[2][0], 1, &tgt[2][0], 1);
713 Vmath::Vcopy(n, &src[0][0], 1, &tgt[3][0], 1);
714 break;
715 case 3:
716 {
717 int a, b, c, d, e, i, j;
718
719 // Compute g^{ij} by computing Cofactors(g_ij)^T
720 for (i = 0; i < m_expDim; ++i)
721 {
722 for (j = 0; j < m_expDim; ++j)
723 {
724 a = ((i + 1) % m_expDim) * m_expDim + ((j + 1) % m_expDim);
725 b = ((i + 1) % m_expDim) * m_expDim + ((j + 2) % m_expDim);
726 c = ((i + 2) % m_expDim) * m_expDim + ((j + 1) % m_expDim);
727 d = ((i + 2) % m_expDim) * m_expDim + ((j + 2) % m_expDim);
728 e = j * m_expDim + i;
729 Vmath::Vvtvvtm(n, &src[a][0], 1, &src[d][0], 1, &src[b][0],
730 1, &src[c][0], 1, &tgt[e][0], 1);
731 }
732 }
733 break;
734 }
735 }
736}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
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, 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 598 of file GeomFactors.cpp.

599{
600 // Jacobian test only makes sense when expdim = coorddim
601 // If one-dimensional then element is valid.
602 if (m_coordDim != m_expDim || m_expDim == 1)
603 {
604 m_valid = true;
605 return;
606 }
607
609 int nqtot = 1;
610 for (int i = 0; i < m_expDim; ++i)
611 {
612 p[i] = m_xmap->GetBasis(i)->GetPointsKey();
613 nqtot *= p[i].GetNumPoints();
614 }
615 int pts = (m_type == eRegular || m_type == eMovingRegular) ? 1 : nqtot;
616 Array<OneD, NekDouble> jac(pts, 0.0);
617
618 DerivStorage deriv = GetDeriv(p);
619
620 switch (m_expDim)
621 {
622 case 2:
623 {
624 Vmath::Vvtvvtm(pts, &deriv[0][0][0], 1, &deriv[1][1][0], 1,
625 &deriv[1][0][0], 1, &deriv[0][1][0], 1, &jac[0], 1);
626 break;
627 }
628 case 3:
629 {
630 Array<OneD, NekDouble> tmp(pts, 0.0);
631
632 Vmath::Vvtvvtm(pts, &deriv[1][1][0], 1, &deriv[2][2][0], 1,
633 &deriv[2][1][0], 1, &deriv[1][2][0], 1, &tmp[0], 1);
634 Vmath::Vvtvp(pts, &deriv[0][0][0], 1, &tmp[0], 1, &jac[0], 1,
635 &jac[0], 1);
636
637 Vmath::Vvtvvtm(pts, &deriv[2][1][0], 1, &deriv[0][2][0], 1,
638 &deriv[0][1][0], 1, &deriv[2][2][0], 1, &tmp[0], 1);
639 Vmath::Vvtvp(pts, &deriv[1][0][0], 1, &tmp[0], 1, &jac[0], 1,
640 &jac[0], 1);
641
642 Vmath::Vvtvvtm(pts, &deriv[0][1][0], 1, &deriv[1][2][0], 1,
643 &deriv[1][1][0], 1, &deriv[0][2][0], 1, &tmp[0], 1);
644 Vmath::Vvtvp(pts, &deriv[2][0][0], 1, &tmp[0], 1, &jac[0], 1,
645 &jac[0], 1);
646
647 break;
648 }
649 }
650
651 if (Vmath::Vmin(pts, &jac[0], 1) < 0)
652 {
653 m_valid = false;
654 }
655}
DerivStorage GetDeriv(const LibUtilities::PointsKeyVector &keyTgt)
Return the derivative of the mapping with respect to the reference coordinates, .
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:64
std::vector< double > p(NPUPPER)
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, 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 }
193
194 switch (m_expDim)
195 {
196 case 1:
197 m_xmap->StdPhysDeriv(tmp, d_map[0][i]);
198 break;
199 case 2:
200 m_xmap->StdPhysDeriv(tmp, d_map[0][i], d_map[1][i]);
201 break;
202 case 3:
203 m_xmap->StdPhysDeriv(tmp, d_map[0][i], d_map[1][i],
204 d_map[2][i]);
205 break;
206 }
207 }
208
209 for (i = 0; i < m_coordDim; ++i)
210 {
211 // Interpolate the derivatives:
212 // - from the points as defined in the mapping ('Coords')
213 // - to the points at which we want to know the metrics
214 // ('tbasis')
215 bool same = true;
216 for (j = 0; j < m_expDim; ++j)
217 {
218 same = same && (map_points[j] == keyTgt[j]);
219 }
220 if (same)
221 {
222 for (j = 0; j < m_expDim; ++j)
223 {
224 deriv[j][i] = d_map[j][i];
225 }
226 }
227 else
228 {
229 for (j = 0; j < m_expDim; ++j)
230 {
231 Interp(map_points, d_map[j][i], keyTgt, deriv[j][i]);
232 }
233 }
234 }
235
236 return deriv;
237}
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 397 of file GeomFactors.cpp.

399{
400 ASSERTL1(keyTgt.size() == m_expDim,
401 "Dimension of target point distribution does not match "
402 "expansion dimension.");
403
404 int i = 0, j = 0, k = 0, l = 0;
405 int ptsTgt = 1;
406
407 if (m_type == eDeformed)
408 {
409 // Allocate storage and compute number of points
410 for (i = 0; i < m_expDim; ++i)
411 {
412 ptsTgt *= keyTgt[i].GetNumPoints();
413 }
414 }
415
416 // Get derivative at geometry points
417 DerivStorage deriv = ComputeDeriv(keyTgt);
418
419 Array<TwoD, NekDouble> tmp(m_expDim * m_expDim, ptsTgt, 0.0);
420 Array<TwoD, NekDouble> gmat(m_expDim * m_expDim, ptsTgt, 0.0);
421 Array<OneD, NekDouble> jac(ptsTgt, 0.0);
422 Array<TwoD, NekDouble> factors(m_expDim * m_coordDim, ptsTgt, 0.0);
423
424 // Compute g_{ij} as t_i \cdot t_j and store in tmp
425 for (i = 0, l = 0; i < m_expDim; ++i)
426 {
427 for (j = 0; j < m_expDim; ++j, ++l)
428 {
429 for (k = 0; k < m_coordDim; ++k)
430 {
431 Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1, &deriv[j][k][0], 1,
432 &tmp[l][0], 1, &tmp[l][0], 1);
433 }
434 }
435 }
436
437 Adjoint(tmp, gmat);
438
439 // Compute g = det(g_{ij}) (= Jacobian squared) and store
440 // temporarily in m_jac.
441 for (i = 0; i < m_expDim; ++i)
442 {
443 Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i * m_expDim][0], 1, &jac[0],
444 1, &jac[0], 1);
445 }
446
447 for (i = 0; i < m_expDim * m_expDim; ++i)
448 {
449 Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
450 }
451
452 // Compute the Jacobian = sqrt(g)
453 Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
454
455 // Compute the derivative factors
456 for (k = 0, l = 0; k < m_coordDim; ++k)
457 {
458 for (j = 0; j < m_expDim; ++j, ++l)
459 {
460 for (i = 0; i < m_expDim; ++i)
461 {
462 Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
463 &gmat[m_expDim * i + j][0], 1, &factors[l][0], 1,
464 &factors[l][0], 1);
465 }
466 }
467 }
468
469 return factors;
470}
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, 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 333 of file GeomFactors.cpp.

335{
336 ASSERTL1(keyTgt.size() == m_expDim,
337 "Dimension of target point distribution does not match "
338 "expansion dimension.");
339
340 int i = 0, j = 0, k = 0, l = 0;
341 int ptsTgt = 1;
342
343 if (m_type == eDeformed)
344 {
345 // Allocate storage and compute number of points
346 for (i = 0; i < m_expDim; ++i)
347 {
348 ptsTgt *= keyTgt[i].GetNumPoints();
349 }
350 }
351
352 // Get derivative at geometry points
353 DerivStorage deriv = ComputeDeriv(keyTgt);
354
355 Array<TwoD, NekDouble> tmp(m_expDim * m_expDim, ptsTgt, 0.0);
356 Array<TwoD, NekDouble> gmat(m_expDim * m_expDim, ptsTgt, 0.0);
357 Array<OneD, NekDouble> jac(ptsTgt, 0.0);
358
359 // Compute g_{ij} as t_i \cdot t_j and store in tmp
360 for (i = 0, l = 0; i < m_expDim; ++i)
361 {
362 for (j = 0; j < m_expDim; ++j, ++l)
363 {
364 for (k = 0; k < m_coordDim; ++k)
365 {
366 Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1, &deriv[j][k][0], 1,
367 &tmp[l][0], 1, &tmp[l][0], 1);
368 }
369 }
370 }
371
372 Adjoint(tmp, gmat);
373
374 // Compute g = det(g_{ij}) (= Jacobian squared) and store
375 // temporarily in m_jac.
376 for (i = 0; i < m_expDim; ++i)
377 {
378 Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i * m_expDim][0], 1, &jac[0],
379 1, &jac[0], 1);
380 }
381
382 for (i = 0; i < m_expDim * m_expDim; ++i)
383 {
384 Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
385 }
386
387 return gmat;
388}

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

252{
253 ASSERTL1(keyTgt.size() == m_expDim,
254 "Dimension of target point distribution does not match "
255 "expansion dimension.");
256
257 // A point always has a unit jacobian
258 if (m_expDim == 0)
259 {
260 return Array<OneD, NekDouble>(1, 1.0);
261 }
262
263 int i = 0, j = 0, k = 0, l = 0;
264 int ptsTgt = 1;
265
266 if (m_type == eDeformed)
267 {
268 // Allocate storage and compute number of points
269 for (i = 0; i < m_expDim; ++i)
270 {
271 ptsTgt *= keyTgt[i].GetNumPoints();
272 }
273 }
274
275 // Get derivative at geometry points
276 DerivStorage deriv = ComputeDeriv(keyTgt);
277
278 Array<TwoD, NekDouble> tmp(m_expDim * m_expDim, ptsTgt, 0.0);
279 Array<TwoD, NekDouble> gmat(m_expDim * m_expDim, ptsTgt, 0.0);
280 Array<OneD, NekDouble> jac(ptsTgt, 0.0);
281
282 // Compute g_{ij} as t_i \cdot t_j and store in tmp
283 for (i = 0, l = 0; i < m_expDim; ++i)
284 {
285 for (j = 0; j < m_expDim; ++j, ++l)
286 {
287 for (k = 0; k < m_coordDim; ++k)
288 {
289 Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1, &deriv[j][k][0], 1,
290 &tmp[l][0], 1, &tmp[l][0], 1);
291 }
292 }
293 }
294
295 Adjoint(tmp, gmat);
296
297 // Compute g = det(g_{ij}) (= Jacobian squared) and store
298 // temporarily in m_jac.
299 for (i = 0; i < m_expDim; ++i)
300 {
301 Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i * m_expDim][0], 1, &jac[0],
302 1, &jac[0], 1);
303 }
304
305 // Compute the Jacobian = sqrt(g)
306 Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
307
308 return jac;
309}

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

477{
478 ASSERTL1(keyTgt.size() == m_expDim,
479 "Dimension of target point distribution does not match "
480 "expansion dimension.");
481
482 int i = 0, k = 0;
483 int ptsTgt = 1;
484 int nq = 1;
485
486 for (i = 0; i < m_expDim; ++i)
487 {
488 nq *= keyTgt[i].GetNumPoints();
489 }
490
491 if (m_type == eDeformed)
492 {
493 // Allocate storage and compute number of points
494 for (i = 0; i < m_expDim; ++i)
495 {
496 ptsTgt *= keyTgt[i].GetNumPoints();
497 }
498 }
499
500 // Get derivative at geometry points
501 DerivStorage deriv = ComputeDeriv(keyTgt);
502
503 // number of moving frames is requited to be 3, even for surfaces
504 int MFdim = 3;
505
506 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtmp(MFdim);
507
508 // Compute g_{ij} as t_i \cdot t_j and store in tmp
509 for (i = 0; i < MFdim; ++i)
510 {
511 MFtmp[i] = Array<OneD, Array<OneD, NekDouble>>(m_coordDim);
512 for (k = 0; k < m_coordDim; ++k)
513 {
514 MFtmp[i][k] = Array<OneD, NekDouble>(nq);
515 }
516 }
517
518 // Compute g_{ij} as t_i \cdot t_j and store in tmp
519 for (i = 0; i < MFdim - 1; ++i)
520 {
521 for (k = 0; k < m_coordDim; ++k)
522 {
523 if (m_type == eDeformed)
524 {
525 Vmath::Vcopy(ptsTgt, &deriv[i][k][0], 1, &MFtmp[i][k][0], 1);
526 }
527 else
528 {
529 Vmath::Fill(nq, deriv[i][k][0], MFtmp[i][k], 1);
530 }
531 }
532 }
533
534 // Direction of MF1 is preserved: MF2 is considered in the same
535 // tangent plane as MF1. MF3 is computed by cross product of MF1
536 // and MF2. MF2 is consequently computed as the cross product of
537 // MF3 and MF1.
538 Array<OneD, Array<OneD, NekDouble>> PrincipleDir(m_coordDim);
539 for (k = 0; k < m_coordDim; k++)
540 {
541 PrincipleDir[k] = Array<OneD, NekDouble>(nq);
542 }
543
544 if (!(MMFdir == eLOCAL))
545 {
546 ComputePrincipleDirection(keyTgt, MMFdir, factors, PrincipleDir);
547 }
548
549 // MF3 = MF1 \times MF2
550 VectorCrossProd(MFtmp[0], MFtmp[1], MFtmp[2]);
551
552 // Normalizing MF3
553 VectorNormalise(MFtmp[2]);
554
555 if (!(MMFdir == eLOCAL))
556 {
557 Array<OneD, NekDouble> temp(nq, 0.0);
558
559 // Reorient MF1 along the PrincipleDir
560 for (i = 0; i < m_coordDim; ++i)
561 {
562 Vmath::Vvtvp(nq, MFtmp[2][i], 1, PrincipleDir[i], 1, temp, 1, temp,
563 1);
564 }
565 Vmath::Neg(nq, temp, 1);
566
567 // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
568 for (i = 0; i < m_coordDim; ++i)
569 {
570 Vmath::Vvtvp(nq, temp, 1, MFtmp[2][i], 1, PrincipleDir[i], 1,
571 MFtmp[0][i], 1);
572 }
573 }
574
575 // Normalizing MF1
576 VectorNormalise(MFtmp[0]);
577
578 // MF2 = MF3 \times MF1
579 VectorCrossProd(MFtmp[2], MFtmp[0], MFtmp[1]);
580
581 // Normalizing MF2
582 VectorNormalise(MFtmp[1]);
583
584 for (i = 0; i < MFdim; ++i)
585 {
586 for (k = 0; k < m_coordDim; ++k)
587 {
588 Vmath::Vcopy(nq, &MFtmp[i][k][0], 1,
589 &movingframes[i * m_coordDim + k][0], 1);
590 }
591 }
592}
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, 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 741 of file GeomFactors.cpp.

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

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

Referenced by ComputeMovingFrames().

◆ GetCoordim()

int Nektar::SpatialDomains::GeomFactors::GetCoordim ( ) const
inline

Return the number of dimensions of the coordinate system.

This is greater than or equal to the expansion dimension.

Returns
The dimension of the coordinate system.

Definition at line 305 of file GeomFactors.h.

306{
307 return m_coordDim;
308}

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

216{
217 return ComputeDeriv(tpoints);
218}

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

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

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

252{
253 return ComputeGmat(keyTgt);
254}
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 296 of file GeomFactors.h.

297{
298 return m_type;
299}

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

326{
327 LibUtilities::PointsKeyVector ptsKeys = m_xmap->GetPointsKeys();
328 const Array<OneD, const NekDouble> jac = GetJac(ptsKeys);
329
330 size_t hash = 0;
332 if (m_type == eDeformed)
333 {
334 hash_range(hash, jac.begin(), jac.end());
335 }
336 else
337 {
338 hash_combine(hash, jac[0]);
339 }
340 return hash;
341}
const Array< OneD, const NekDouble > GetJac(const LibUtilities::PointsKeyVector &keyTgt)
Return the Jacobian of the mapping and cache the result.
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 229 of file GeomFactors.h.

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

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

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

344{
345 return m_xmap;
346}

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

667{
668 ASSERTL1(src_points.size() == tgt_points.size(),
669 "Dimension of target point distribution does not match "
670 "expansion dimension.");
671
672 switch (m_expDim)
673 {
674 case 1:
675 LibUtilities::Interp1D(src_points[0], src, tgt_points[0], tgt);
676 break;
677 case 2:
678 LibUtilities::Interp2D(src_points[0], src_points[1], src,
679 tgt_points[0], tgt_points[1], tgt);
680 break;
681 case 3:
682 LibUtilities::Interp3D(src_points[0], src_points[1], src_points[2],
683 src, tgt_points[0], tgt_points[1],
684 tgt_points[2], tgt);
685 break;
686 }
687}
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 315 of file GeomFactors.h.

316{
317 return m_valid;
318}

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

950{
951 ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
952 ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
953 ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
954
955 int nq = v1[0].size();
956 Array<OneD, NekDouble> temp(nq);
957
958 Vmath::Vmul(nq, v1[2], 1, v2[1], 1, temp, 1);
959 Vmath::Vvtvm(nq, v1[1], 1, v2[2], 1, temp, 1, v3[0], 1);
960
961 Vmath::Vmul(nq, v1[0], 1, v2[2], 1, temp, 1);
962 Vmath::Vvtvm(nq, v1[2], 1, v2[0], 1, temp, 1, v3[1], 1);
963
964 Vmath::Vmul(nq, v1[1], 1, v2[0], 1, temp, 1);
965 Vmath::Vvtvm(nq, v1[0], 1, v2[1], 1, temp, 1, v3[2], 1);
966}
#define ASSERTL0(condition, msg)
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 909 of file GeomFactors.cpp.

910{
911 int ndim = array.size();
912 ASSERTL0(ndim > 0, "Number of components must be > 0.");
913 for (int i = 1; i < ndim; ++i)
914 {
915 ASSERTL0(array[i].size() == array[0].size(),
916 "Array size mismatch in coordinates.");
917 }
918
919 int nq = array[0].size();
920 Array<OneD, NekDouble> norm(nq, 0.0);
921
922 // Compute the norm of each vector.
923 for (int i = 0; i < ndim; ++i)
924 {
925 Vmath::Vvtvp(nq, array[i], 1, array[i], 1, norm, 1, norm, 1);
926 }
927
928 Vmath::Vsqrt(nq, norm, 1, norm, 1);
929
930 // Normalise the vectors by the norm
931 for (int i = 0; i < ndim; ++i)
932 {
933 Vmath::Vdiv(nq, array[i], 1, norm, 1, array[i], 1);
934 }
935}

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

Referenced by ComputeMovingFrames().

Friends And Related Symbol 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 139 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 144 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 141 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 134 of file GeomFactors.h.

◆ m_type

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

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

Definition at line 125 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 131 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 137 of file GeomFactors.h.

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