Nektar++
Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | List of all members
Nektar::SpatialDomains::Geometry3D Class Reference

3D geometry information More...

#include <Geometry3D.h>

Inheritance diagram for Nektar::SpatialDomains::Geometry3D:
[legend]

Public Member Functions

 Geometry3D ()
 
 Geometry3D (const int coordim)
 
virtual ~Geometry3D ()
 
- Public Member Functions inherited from Nektar::SpatialDomains::Geometry
 Geometry ()
 Default constructor. More...
 
 Geometry (int coordim)
 Constructor when supplied a coordinate dimension. More...
 
virtual ~Geometry ()
 Default destructor. More...
 
int GetCoordim () const
 Return the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded). More...
 
void SetCoordim (int coordim)
 Sets the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded). More...
 
GeomFactorsSharedPtr GetGeomFactors ()
 Get the geometric factors for this object, generating them if required. More...
 
GeomFactorsSharedPtr GetRefGeomFactors (const Array< OneD, const LibUtilities::BasisSharedPtr > &tbasis)
 
GeomFactorsSharedPtr GetMetricInfo ()
 Get the geometric factors for this object. More...
 
LibUtilities::ShapeType GetShapeType (void)
 Get the geometric shape type of this object. More...
 
int GetGlobalID (void) const
 Get the ID of this object. More...
 
void SetGlobalID (int globalid)
 Set the ID of this object. More...
 
int GetVid (int i) const
 Get the ID of vertex i of this object. More...
 
int GetEid (int i) const
 Get the ID of edge i of this object. More...
 
int GetFid (int i) const
 Get the ID of face i of this object. More...
 
int GetTid (int i) const
 Get the ID of trace i of this object. More...
 
PointGeomSharedPtr GetVertex (int i) const
 Returns vertex i of this object. More...
 
Geometry1DSharedPtr GetEdge (int i) const
 Returns edge i of this object. More...
 
Geometry2DSharedPtr GetFace (int i) const
 Returns face i of this object. More...
 
StdRegions::Orientation GetEorient (const int i) const
 Returns the orientation of edge i with respect to the ordering of edges in the standard element. More...
 
StdRegions::Orientation GetForient (const int i) const
 Returns the orientation of face i with respect to the ordering of faces in the standard element. More...
 
int GetNumVerts () const
 Get the number of vertices of this object. More...
 
int GetNumEdges () const
 Get the number of edges of this object. More...
 
int GetNumFaces () const
 Get the number of faces of this object. More...
 
int GetShapeDim () const
 Get the object's shape dimension. More...
 
StdRegions::StdExpansionSharedPtr GetXmap () const
 Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element. More...
 
const Array< OneD, const NekDouble > & GetCoeffs (const int i) const
 Return the coefficients of the transformation Geometry::m_xmap in coordinate direction i. More...
 
void FillGeom ()
 Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements. More...
 
std::array< NekDouble, 6 > GetBoundingBox ()
 Generates the bounding box for the element. More...
 
bool ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
 Determine whether an element contains a particular Cartesian coordinate \((x,y,z)\). More...
 
bool ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol)
 Determine whether an element contains a particular Cartesian coordinate \((x,y,z)\). More...
 
bool ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &dist)
 Determine whether an element contains a particular Cartesian coordinate \(\vec{x} = (x,y,z)\). More...
 
NekDouble GetLocCoords (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
 Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object. More...
 
NekDouble GetCoord (const int i, const Array< OneD, const NekDouble > &Lcoord)
 Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i. More...
 
bool MinMaxCheck (const Array< OneD, const NekDouble > &gloCoord)
 Check if given global coord is within the BoundingBox of the element. More...
 
bool ClampLocCoords (Array< OneD, NekDouble > &locCoord, NekDouble tol)
 Clamp local coords to be within standard regions [-1, 1]^dim. More...
 
int GetVertexEdgeMap (int i, int j) const
 Returns the standard element edge IDs that are connected to a given vertex. More...
 
int GetVertexFaceMap (int i, int j) const
 Returns the standard element face IDs that are connected to a given vertex. More...
 
int GetEdgeFaceMap (int i, int j) const
 Returns the standard element edge IDs that are connected to a given face. More...
 
int GetDir (const int i, const int j=0) const
 Returns the element coordinate direction corresponding to a given face coordinate direction. More...
 
void Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)
 Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors. More...
 
void Setup ()
 

Static Public Attributes

static const int kDim = 3
 

Protected Member Functions

virtual NekDouble v_GetLocCoords (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
 Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object. More...
 
void NewtonIterationForLocCoord (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, const Array< OneD, const NekDouble > &ptsz, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)
 
virtual void v_FillGeom ()
 Put all quadrature information into face/edge structure and backward transform. More...
 
virtual NekDouble v_GetCoord (const int i, const Array< OneD, const NekDouble > &Lcoord)
 Given local collapsed coordinate Lcoord return the value of physical coordinate in direction i. More...
 
virtual int v_GetShapeDim () const
 Get the object's shape dimension. More...
 
virtual int v_GetNumVerts () const
 Get the number of vertices of this object. More...
 
virtual int v_GetNumEdges () const
 Get the number of edges of this object. More...
 
virtual int v_GetNumFaces () const
 Get the number of faces of this object. More...
 
virtual PointGeomSharedPtr v_GetVertex (int i) const
 
virtual Geometry1DSharedPtr v_GetEdge (int i) const
 Returns edge i of this object. More...
 
virtual Geometry2DSharedPtr v_GetFace (int i) const
 Returns face i of this object. More...
 
virtual StdRegions::Orientation v_GetEorient (const int i) const
 Returns the orientation of edge i with respect to the ordering of edges in the standard element. More...
 
virtual StdRegions::Orientation v_GetForient (const int i) const
 Returns the orientation of face i with respect to the ordering of faces in the standard element. More...
 
- Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
void GenGeomFactors ()
 Handles generation of geometry factors. More...
 
virtual StdRegions::StdExpansionSharedPtr v_GetXmap () const
 Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element. More...
 
virtual bool v_ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &dist)
 
virtual int v_GetVertexEdgeMap (int i, int j) const
 Returns the standard element edge IDs that are connected to a given vertex. More...
 
virtual int v_GetVertexFaceMap (int i, int j) const
 Returns the standard element face IDs that are connected to a given vertex. More...
 
virtual int v_GetEdgeFaceMap (int i, int j) const
 Returns the standard element edge IDs that are connected to a given face. More...
 
virtual int v_GetDir (const int faceidx, const int facedir) const
 Returns the element coordinate direction corresponding to a given face coordinate direction. More...
 
virtual void v_Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)
 Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors. More...
 
virtual void v_Setup ()
 
virtual void v_GenGeomFactors ()=0
 
void SetUpCoeffs (const int nCoeffs)
 Initialise the Geometry::m_coeffs array. More...
 

Protected Attributes

PointGeomVector m_verts
 
SegGeomVector m_edges
 
Geometry2DVector m_faces
 
std::vector< StdRegions::Orientationm_eorient
 
std::vector< StdRegions::Orientationm_forient
 
int m_eid
 
bool m_ownverts
 
- Protected Attributes inherited from Nektar::SpatialDomains::Geometry
int m_coordim
 Coordinate dimension of this geometry object. More...
 
GeomFactorsSharedPtr m_geomFactors
 Geometric factors. More...
 
GeomState m_geomFactorsState
 State of the geometric factors. More...
 
StdRegions::StdExpansionSharedPtr m_xmap
 \(\chi\) mapping containing isoparametric transformation. More...
 
GeomState m_state
 Enumeration to dictate whether coefficients are filled. More...
 
bool m_setupState
 Wether or not the setup routines have been run. More...
 
GeomType m_geomType
 Type of geometry. More...
 
LibUtilities::ShapeType m_shapeType
 Type of shape. More...
 
int m_globalID
 Global ID. More...
 
Array< OneD, Array< OneD, NekDouble > > m_coeffs
 Array containing expansion coefficients of m_xmap. More...
 
Array< OneD, NekDoublem_boundingBox
 Array containing bounding box. More...
 

Additional Inherited Members

- Static Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
static GeomFactorsSharedPtr ValidateRegGeomFactor (GeomFactorsSharedPtr geomFactor)
 Check to see if a geometric factor has already been created that contains the same regular information. More...
 
- Static Protected Attributes inherited from Nektar::SpatialDomains::Geometry
static GeomFactorsVector m_regGeomFactorsManager
 

Detailed Description

3D geometry information

Definition at line 67 of file Geometry3D.h.

Constructor & Destructor Documentation

◆ Geometry3D() [1/2]

Nektar::SpatialDomains::Geometry3D::Geometry3D ( )

Definition at line 51 of file Geometry3D.cpp.

52 {
53 }

◆ Geometry3D() [2/2]

Nektar::SpatialDomains::Geometry3D::Geometry3D ( const int  coordim)

Definition at line 55 of file Geometry3D.cpp.

55  : Geometry(coordim)
56 {
57  ASSERTL0(m_coordim > 2,
58  "Coordinate dimension should be at least 3 for a 3D geometry.");
59 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Geometry()
Default constructor.
Definition: Geometry.cpp:54
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:185

References ASSERTL0, and Nektar::SpatialDomains::Geometry::m_coordim.

◆ ~Geometry3D()

Nektar::SpatialDomains::Geometry3D::~Geometry3D ( )
virtual

Definition at line 61 of file Geometry3D.cpp.

62 {
63 }

Member Function Documentation

◆ NewtonIterationForLocCoord()

void Nektar::SpatialDomains::Geometry3D::NewtonIterationForLocCoord ( const Array< OneD, const NekDouble > &  coords,
const Array< OneD, const NekDouble > &  ptsx,
const Array< OneD, const NekDouble > &  ptsy,
const Array< OneD, const NekDouble > &  ptsz,
Array< OneD, NekDouble > &  Lcoords,
NekDouble dist 
)
protected

Definition at line 72 of file Geometry3D.cpp.

79 {
80  // maximum iterations for convergence
81  const int MaxIterations = 51;
82  // |x-xp|^2 < EPSILON error tolerance
83  const NekDouble Tol = 1.e-8;
84  // |r,s| > LcoordDIV stop the search
85  const NekDouble LcoordDiv = 15.0;
86 
87  Array<OneD, const NekDouble> Jac =
88  m_geomFactors->GetJac(m_xmap->GetPointsKeys());
89 
90  NekDouble ScaledTol = Vmath::Vsum(Jac.size(), Jac, 1) /
91  ((NekDouble)Jac.size());
92  ScaledTol *= Tol;
93 
94  NekDouble xmap, ymap, zmap, F1, F2, F3;
95 
96  NekDouble derx_1, derx_2, derx_3, dery_1, dery_2, dery_3, derz_1, derz_2,
97  derz_3, jac;
98 
99  // save intiial guess for later reference if required.
100  NekDouble init0 = Lcoords[0], init1 = Lcoords[1], init2 = Lcoords[2];
101 
102  Array<OneD, NekDouble> DxD1(ptsx.size());
103  Array<OneD, NekDouble> DxD2(ptsx.size());
104  Array<OneD, NekDouble> DxD3(ptsx.size());
105  Array<OneD, NekDouble> DyD1(ptsx.size());
106  Array<OneD, NekDouble> DyD2(ptsx.size());
107  Array<OneD, NekDouble> DyD3(ptsx.size());
108  Array<OneD, NekDouble> DzD1(ptsx.size());
109  Array<OneD, NekDouble> DzD2(ptsx.size());
110  Array<OneD, NekDouble> DzD3(ptsx.size());
111 
112  // Ideally this will be stored in m_geomfactors
113  m_xmap->PhysDeriv(ptsx, DxD1, DxD2, DxD3);
114  m_xmap->PhysDeriv(ptsy, DyD1, DyD2, DyD3);
115  m_xmap->PhysDeriv(ptsz, DzD1, DzD2, DzD3);
116 
117  int cnt = 0;
118  Array<OneD, DNekMatSharedPtr> I(3);
119  Array<OneD, NekDouble> eta(3);
120 
121  F1 = F2 = F3 = 2000; // Starting value of Function
122  NekDouble resid;
123  while (cnt++ < MaxIterations)
124  {
125  // evaluate lagrange interpolant at Lcoords
126  m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
127  I[0] = m_xmap->GetBasis(0)->GetI(eta);
128  I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
129  I[2] = m_xmap->GetBasis(2)->GetI(eta + 2);
130 
131  // calculate the global point `corresponding to Lcoords
132  xmap = m_xmap->PhysEvaluate(I, ptsx);
133  ymap = m_xmap->PhysEvaluate(I, ptsy);
134  zmap = m_xmap->PhysEvaluate(I, ptsz);
135 
136  F1 = coords[0] - xmap;
137  F2 = coords[1] - ymap;
138  F3 = coords[2] - zmap;
139 
140  if (F1 * F1 + F2 * F2 + F3 * F3 < ScaledTol)
141  {
142  resid = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
143  break;
144  }
145 
146  // Interpolate derivative metric at Lcoords
147  derx_1 = m_xmap->PhysEvaluate(I, DxD1);
148  derx_2 = m_xmap->PhysEvaluate(I, DxD2);
149  derx_3 = m_xmap->PhysEvaluate(I, DxD3);
150  dery_1 = m_xmap->PhysEvaluate(I, DyD1);
151  dery_2 = m_xmap->PhysEvaluate(I, DyD2);
152  dery_3 = m_xmap->PhysEvaluate(I, DyD3);
153  derz_1 = m_xmap->PhysEvaluate(I, DzD1);
154  derz_2 = m_xmap->PhysEvaluate(I, DzD2);
155  derz_3 = m_xmap->PhysEvaluate(I, DzD3);
156 
157  jac = derx_1 * (dery_2 * derz_3 - dery_3 * derz_2) -
158  derx_2 * (dery_1 * derz_3 - dery_3 * derz_1) +
159  derx_3 * (dery_1 * derz_2 - dery_2 * derz_1);
160 
161  // use analytical inverse of derivitives which are also similar to
162  // those of metric factors.
163  Lcoords[0] =
164  Lcoords[0] +
165  ((dery_2 * derz_3 - dery_3 * derz_2) * (coords[0] - xmap) -
166  (derx_2 * derz_3 - derx_3 * derz_2) * (coords[1] - ymap) +
167  (derx_2 * dery_3 - derx_3 * dery_2) * (coords[2] - zmap)) /
168  jac;
169 
170  Lcoords[1] =
171  Lcoords[1] -
172  ((dery_1 * derz_3 - dery_3 * derz_1) * (coords[0] - xmap) -
173  (derx_1 * derz_3 - derx_3 * derz_1) * (coords[1] - ymap) +
174  (derx_1 * dery_3 - derx_3 * dery_1) * (coords[2] - zmap)) /
175  jac;
176 
177  Lcoords[2] =
178  Lcoords[2] +
179  ((dery_1 * derz_2 - dery_2 * derz_1) * (coords[0] - xmap) -
180  (derx_1 * derz_2 - derx_2 * derz_1) * (coords[1] - ymap) +
181  (derx_1 * dery_2 - derx_2 * dery_1) * (coords[2] - zmap)) /
182  jac;
183 
184  if( !(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1]) &&
185  std::isfinite(Lcoords[2])) )
186  {
187  dist = 1e16;
188  std::ostringstream ss;
189  ss << "nan or inf found in NewtonIterationForLocCoord in element "
190  << GetGlobalID();
191  WARNINGL1(false, ss.str());
192  return;
193  }
194  if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
195  fabs(Lcoords[0]) > LcoordDiv)
196  {
197  break; // lcoords have diverged so stop iteration
198  }
199  }
200 
201  m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
202  if(ClampLocCoords(eta, 0.))
203  {
204  I[0] = m_xmap->GetBasis(0)->GetI(eta);
205  I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
206  I[2] = m_xmap->GetBasis(2)->GetI(eta + 2);
207  // calculate the global point corresponding to Lcoords
208  xmap = m_xmap->PhysEvaluate(I, ptsx);
209  ymap = m_xmap->PhysEvaluate(I, ptsy);
210  zmap = m_xmap->PhysEvaluate(I, ptsz);
211  F1 = coords[0] - xmap;
212  F2 = coords[1] - ymap;
213  F3 = coords[2] - zmap;
214  dist = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
215  }
216  else
217  {
218  dist = 0.;
219  }
220 
221  if (cnt >= MaxIterations)
222  {
223  Array<OneD, NekDouble> collCoords(3);
224  m_xmap->LocCoordToLocCollapsed(Lcoords, collCoords);
225 
226  // if coordinate is inside element dump error!
227  if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
228  (collCoords[1] >= -1.0 && collCoords[1] <= 1.0) &&
229  (collCoords[2] >= -1.0 && collCoords[2] <= 1.0))
230  {
231  std::ostringstream ss;
232 
233  ss << "Reached MaxIterations (" << MaxIterations
234  << ") in Newton iteration ";
235  ss << "Init value (" << setprecision(4) << init0 << "," << init1
236  << "," << init2 << ") ";
237  ss << "Fin value (" << Lcoords[0] << "," << Lcoords[1] << ","
238  << Lcoords[2] << ") ";
239  ss << "Resid = " << resid << " Tolerance = " << sqrt(ScaledTol);
240 
241  WARNINGL1(cnt < MaxIterations, ss.str());
242  }
243  }
244 }
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:251
bool ClampLocCoords(Array< OneD, NekDouble > &locCoord, NekDouble tol)
Clamp local coords to be within standard regions [-1, 1]^dim.
Definition: Geometry.cpp:486
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:319
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:191
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:187
double NekDouble
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:846
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

References Nektar::SpatialDomains::Geometry::ClampLocCoords(), Nektar::SpatialDomains::Geometry::GetGlobalID(), Nektar::SpatialDomains::Geometry::m_geomFactors, Nektar::SpatialDomains::Geometry::m_xmap, tinysimd::sqrt(), Vmath::Vsum(), and WARNINGL1.

Referenced by v_GetLocCoords().

◆ v_FillGeom()

void Nektar::SpatialDomains::Geometry3D::v_FillGeom ( )
protectedvirtual

Put all quadrature information into face/edge structure and backward transform.

Note verts, edges, and faces are listed according to anticlockwise convention but points in _coeffs have to be in array format from left to right.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 365 of file Geometry3D.cpp.

366 {
367  if (m_state == ePtsFilled)
368  {
369  return;
370  }
371 
372  int i, j, k;
373 
374  for (i = 0; i < m_forient.size(); i++)
375  {
376  m_faces[i]->FillGeom();
377 
378  int nFaceCoeffs = m_faces[i]->GetXmap()->GetNcoeffs();
379 
380  Array<OneD, unsigned int> mapArray(nFaceCoeffs);
381  Array<OneD, int> signArray(nFaceCoeffs);
382 
383  if (m_forient[i] < 9)
384  {
385  m_xmap->GetTraceToElementMap(
386  i,
387  mapArray,
388  signArray,
389  m_forient[i],
390  m_faces[i]->GetXmap()->GetTraceNcoeffs(0),
391  m_faces[i]->GetXmap()->GetTraceNcoeffs(1));
392  }
393  else
394  {
395  m_xmap->GetTraceToElementMap(
396  i,
397  mapArray,
398  signArray,
399  m_forient[i],
400  m_faces[i]->GetXmap()->GetTraceNcoeffs(1),
401  m_faces[i]->GetXmap()->GetTraceNcoeffs(0));
402  }
403 
404  for (j = 0; j < m_coordim; j++)
405  {
406  const Array<OneD, const NekDouble> &coeffs =
407  m_faces[i]->GetCoeffs(j);
408 
409  for (k = 0; k < nFaceCoeffs; k++)
410  {
411  NekDouble v = signArray[k] * coeffs[k];
412  m_coeffs[j][mapArray[k]] = v;
413  }
414  }
415  }
416 
418 }
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:84
GeomState m_state
Enumeration to dictate whether coefficients are filled.
Definition: Geometry.h:193
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:203
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
Definition: Geometry.h:426
@ ePtsFilled
Geometric information has been generated.

References Nektar::SpatialDomains::ePtsFilled, Nektar::SpatialDomains::Geometry::GetXmap(), Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_coordim, m_faces, m_forient, Nektar::SpatialDomains::Geometry::m_state, and Nektar::SpatialDomains::Geometry::m_xmap.

Referenced by Nektar::SpatialDomains::HexGeom::v_GenGeomFactors(), Nektar::SpatialDomains::PrismGeom::v_GenGeomFactors(), Nektar::SpatialDomains::PyrGeom::v_GenGeomFactors(), Nektar::SpatialDomains::TetGeom::v_GenGeomFactors(), and v_GetLocCoords().

◆ v_GetCoord()

NekDouble Nektar::SpatialDomains::Geometry3D::v_GetCoord ( const int  i,
const Array< OneD, const NekDouble > &  Lcoord 
)
protectedvirtual

Given local collapsed coordinate Lcoord return the value of physical coordinate in direction i.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 424 of file Geometry3D.cpp.

426 {
427  ASSERTL1(m_state == ePtsFilled, "Geometry is not in physical space");
428 
429  Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
430  m_xmap->BwdTrans(m_coeffs[i], tmp);
431 
432  return m_xmap->PhysEvaluate(Lcoord, tmp);
433 }
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250

References ASSERTL1, Nektar::SpatialDomains::ePtsFilled, Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_state, and Nektar::SpatialDomains::Geometry::m_xmap.

◆ v_GetEdge()

Geometry1DSharedPtr Nektar::SpatialDomains::Geometry3D::v_GetEdge ( int  i) const
protectedvirtual

Returns edge i of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 464 of file Geometry3D.cpp.

465 {
466  ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
467  "Edge ID must be between 0 and " +
468  boost::lexical_cast<string>(m_edges.size() - 1));
469  return m_edges[i];
470 }
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274

References ASSERTL2, and m_edges.

◆ v_GetEorient()

StdRegions::Orientation Nektar::SpatialDomains::Geometry3D::v_GetEorient ( const int  i) const
inlineprotectedvirtual

Returns the orientation of edge i with respect to the ordering of edges in the standard element.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 478 of file Geometry3D.cpp.

479 {
480  ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
481  "Edge ID must be between 0 and " +
482  boost::lexical_cast<string>(m_edges.size() - 1));
483  return m_eorient[i];
484 }
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:83

References ASSERTL2, m_edges, and m_eorient.

◆ v_GetFace()

Geometry2DSharedPtr Nektar::SpatialDomains::Geometry3D::v_GetFace ( int  i) const
protectedvirtual

Returns face i of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 472 of file Geometry3D.cpp.

473 {
474  ASSERTL2((i >= 0) && (i <= 5), "Edge id must be between 0 and 4");
475  return m_faces[i];
476 }

References ASSERTL2, and m_faces.

◆ v_GetForient()

StdRegions::Orientation Nektar::SpatialDomains::Geometry3D::v_GetForient ( const int  i) const
protectedvirtual

Returns the orientation of face i with respect to the ordering of faces in the standard element.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 486 of file Geometry3D.cpp.

487 {
488  ASSERTL2(i >= 0 && i <= m_faces.size() - 1,
489  "Face ID must be between 0 and " +
490  boost::lexical_cast<string>(m_faces.size() - 1));
491  return m_forient[i];
492 }

References ASSERTL2, m_faces, and m_forient.

◆ v_GetLocCoords()

NekDouble Nektar::SpatialDomains::Geometry3D::v_GetLocCoords ( const Array< OneD, const NekDouble > &  coords,
Array< OneD, NekDouble > &  Lcoords 
)
protectedvirtual

Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object.

For curvilinear and non-affine elements (i.e. where the Jacobian varies as a function of the standard element coordinates), this is a non-linear optimisation problem that requires the use of a Newton iteration. Note therefore that this can be an expensive operation.

Note that, clearly, the provided Cartesian coordinate lie outside the element. The function therefore returns the minimum distance from some position in the element to . Lcoords will also be constrained to fit within the range \([-1,1]^d\) where \( d \) is the dimension of the element.

Parameters
coordsInput Cartesian global coordinates
LcoordsCorresponding local coordinates
Returns
Distance between obtained coordinates and provided ones.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 246 of file Geometry3D.cpp.

248 {
249  NekDouble dist = 0.;
250  if (GetMetricInfo()->GetGtype() == eRegular)
251  {
252  int v1, v2, v3;
256  {
257  v1 = 1;
258  v2 = 3;
259  v3 = 4;
260  }
262  {
263  v1 = 1;
264  v2 = 2;
265  v3 = 3;
266  }
267  else
268  {
269  v1 = 1;
270  v2 = 2;
271  v3 = 3;
272  ASSERTL0(false, "unrecognized 3D element type");
273  }
274  // Point inside tetrahedron
275  PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
276 
277  // Edges
278  PointGeom er0, e10, e20, e30;
279  er0.Sub(r, *m_verts[0]);
280  e10.Sub(*m_verts[v1], *m_verts[0]);
281  e20.Sub(*m_verts[v2], *m_verts[0]);
282  e30.Sub(*m_verts[v3], *m_verts[0]);
283 
284  // Cross products (Normal times area)
285  PointGeom cp1020, cp2030, cp3010;
286  cp1020.Mult(e10, e20);
287  cp2030.Mult(e20, e30);
288  cp3010.Mult(e30, e10);
289 
290  // Barycentric coordinates (relative volume)
291  NekDouble iV = 2. / e30.dot(cp1020); // Hex Volume = {(e30)dot(e10)x(e20)}
292  Lcoords[0] = er0.dot(cp2030) * iV - 1.0;
293  Lcoords[1] = er0.dot(cp3010) * iV - 1.0;
294  Lcoords[2] = er0.dot(cp1020) * iV - 1.0;
295 
296  // Set distance
297  Array<OneD, NekDouble> eta(3, 0.);
298  m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
299  if(ClampLocCoords(eta, 0.))
300  {
301  Array<OneD, NekDouble> xi(3, 0.);
302  m_xmap->LocCollapsedToLocCoord(eta, xi);
303  xi[0] = (xi[0] + 1.) * 0.5; //re-scaled to ratio [0, 1]
304  xi[1] = (xi[1] + 1.) * 0.5;
305  xi[2] = (xi[2] + 1.) * 0.5;
306  for (int i = 0; i < m_coordim; ++i)
307  {
308  NekDouble tmp = xi[0]*e10[i] + xi[1]*e20[i]
309  + xi[2]*e30[i] - er0[i];
310  dist += tmp * tmp;
311  }
312  dist = sqrt(dist);
313  }
314  }
315  else
316  {
317  v_FillGeom();
318  // Determine nearest point of coords to values in m_xmap
319  int npts = m_xmap->GetTotPoints();
320  Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
321  Array<OneD, NekDouble> tmp1(npts), tmp2(npts);
322 
323  m_xmap->BwdTrans(m_coeffs[0], ptsx);
324  m_xmap->BwdTrans(m_coeffs[1], ptsy);
325  m_xmap->BwdTrans(m_coeffs[2], ptsz);
326 
327  const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
328  const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
329  const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);
330 
331  // guess the first local coords based on nearest point
332  Vmath::Sadd(npts, -coords[0], ptsx, 1, tmp1, 1);
333  Vmath::Vmul(npts, tmp1, 1, tmp1, 1, tmp1, 1);
334  Vmath::Sadd(npts, -coords[1], ptsy, 1, tmp2, 1);
335  Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
336  Vmath::Sadd(npts, -coords[2], ptsz, 1, tmp2, 1);
337  Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
338 
339  int min_i = Vmath::Imin(npts, tmp1, 1);
340 
341  // Get Local coordinates
342  int qa = za.size(), qb = zb.size();
343 
344  Array<OneD, NekDouble> eta(3, 0.);
345  eta[2] = zc[min_i / (qa * qb)];
346  min_i = min_i % (qa * qb);
347  eta[1] = zb[min_i / qa];
348  eta[0] = za[min_i % qa];
349  m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
350 
351  // Perform newton iteration to find local coordinates
352  NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords, dist);
353  }
354  return dist;
355 }
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:365
void NewtonIterationForLocCoord(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, const Array< OneD, const NekDouble > &ptsz, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)
Definition: Geometry3D.cpp:72
LibUtilities::ShapeType m_shapeType
Type of shape.
Definition: Geometry.h:199
GeomFactorsSharedPtr GetMetricInfo()
Get the geometric factors for this object.
Definition: Geometry.h:303
@ eRegular
Geometry is straight-sided with constant geometric factors.
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:192
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:513
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:966
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:341

References ASSERTL0, Nektar::SpatialDomains::Geometry::ClampLocCoords(), Nektar::SpatialDomains::PointGeom::dot(), Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePyramid, Nektar::SpatialDomains::eRegular, Nektar::LibUtilities::eTetrahedron, Nektar::SpatialDomains::Geometry::GetMetricInfo(), Vmath::Imin(), Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_coordim, Nektar::SpatialDomains::Geometry::m_shapeType, m_verts, Nektar::SpatialDomains::Geometry::m_xmap, Nektar::SpatialDomains::PointGeom::Mult(), NewtonIterationForLocCoord(), Vmath::Sadd(), tinysimd::sqrt(), Nektar::SpatialDomains::PointGeom::Sub(), v_FillGeom(), Vmath::Vmul(), and Vmath::Vvtvp().

◆ v_GetNumEdges()

int Nektar::SpatialDomains::Geometry3D::v_GetNumEdges ( ) const
protectedvirtual

Get the number of edges of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 449 of file Geometry3D.cpp.

450 {
451  return m_edges.size();
452 }

References m_edges.

◆ v_GetNumFaces()

int Nektar::SpatialDomains::Geometry3D::v_GetNumFaces ( ) const
protectedvirtual

Get the number of faces of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 454 of file Geometry3D.cpp.

455 {
456  return m_faces.size();
457 }

References m_faces.

◆ v_GetNumVerts()

int Nektar::SpatialDomains::Geometry3D::v_GetNumVerts ( ) const
protectedvirtual

Get the number of vertices of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 444 of file Geometry3D.cpp.

445 {
446  return m_verts.size();
447 }

References m_verts.

◆ v_GetShapeDim()

int Nektar::SpatialDomains::Geometry3D::v_GetShapeDim ( ) const
protectedvirtual

Get the object's shape dimension.

For example, a segment is one dimensional and quadrilateral is two dimensional.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 439 of file Geometry3D.cpp.

440 {
441  return 3;
442 }

◆ v_GetVertex()

PointGeomSharedPtr Nektar::SpatialDomains::Geometry3D::v_GetVertex ( int  i) const
protectedvirtual

Implements Nektar::SpatialDomains::Geometry.

Definition at line 459 of file Geometry3D.cpp.

460 {
461  return m_verts[i];
462 }

References m_verts.

Member Data Documentation

◆ kDim

const int Nektar::SpatialDomains::Geometry3D::kDim = 3
static

Definition at line 77 of file Geometry3D.h.

◆ m_edges

SegGeomVector Nektar::SpatialDomains::Geometry3D::m_edges
protected

◆ m_eid

int Nektar::SpatialDomains::Geometry3D::m_eid
protected

Definition at line 85 of file Geometry3D.h.

◆ m_eorient

std::vector<StdRegions::Orientation> Nektar::SpatialDomains::Geometry3D::m_eorient
protected

◆ m_faces

Geometry2DVector Nektar::SpatialDomains::Geometry3D::m_faces
protected

◆ m_forient

std::vector<StdRegions::Orientation> Nektar::SpatialDomains::Geometry3D::m_forient
protected

◆ m_ownverts

bool Nektar::SpatialDomains::Geometry3D::m_ownverts
protected

Definition at line 86 of file Geometry3D.h.

◆ m_verts

PointGeomVector Nektar::SpatialDomains::Geometry3D::m_verts
protected