Nektar++
QuadGeom.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: QuadGeom.cpp
4//
5// For more information, please see: http://www.nektar.info/
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description:
32//
33//
34////////////////////////////////////////////////////////////////////////////////
35
38
43
45{
46
48{
50}
51
52QuadGeom::QuadGeom(const int id, const SegGeomSharedPtr edges[],
53 const CurveSharedPtr curve)
54 : Geometry2D(edges[0]->GetVertex(0)->GetCoordim(), curve)
55{
56 int j;
57
59 m_globalID = id;
60
61 /// Copy the edge shared pointers.
62 m_edges.insert(m_edges.begin(), edges, edges + QuadGeom::kNedges);
63 m_eorient.resize(kNedges);
64
65 for (j = 0; j < kNedges; ++j)
66 {
67 m_eorient[j] =
68 SegGeom::GetEdgeOrientation(*edges[j], *edges[(j + 1) % kNedges]);
69 m_verts.push_back(
70 edges[j]->GetVertex(m_eorient[j] == StdRegions::eForwards ? 0 : 1));
71 }
72
73 for (j = 2; j < kNedges; ++j)
74 {
78 }
79
80 m_coordim = edges[0]->GetVertex(0)->GetCoordim();
81 ASSERTL0(m_coordim > 1, "Cannot call function with dim == 1");
82}
83
85{
86 // From Geometry
89
90 // From QuadGeom
91 m_verts = in.m_verts;
92 m_edges = in.m_edges;
93 for (int i = 0; i < kNedges; i++)
94 {
95 m_eorient[i] = in.m_eorient[i];
96 }
97}
98
100{
101 int order0 = std::max(m_edges[0]->GetXmap()->GetBasis(0)->GetNumModes(),
102 m_edges[2]->GetXmap()->GetBasis(0)->GetNumModes());
103 int order1 = std::max(m_edges[1]->GetXmap()->GetBasis(0)->GetNumModes(),
104 m_edges[3]->GetXmap()->GetBasis(0)->GetNumModes());
105
106 const LibUtilities::BasisKey B0(
108 LibUtilities::PointsKey(order0 + 1,
110 const LibUtilities::BasisKey B1(
112 LibUtilities::PointsKey(order1 + 1,
114
116}
117
119 const Array<OneD, const NekDouble> &Lcoord)
120{
121 ASSERTL1(m_state == ePtsFilled, "Geometry is not in physical space");
122
123 Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
124 m_xmap->BwdTrans(m_coeffs[i], tmp);
125
126 return m_xmap->PhysEvaluate(Lcoord, tmp);
127}
128
130 const QuadGeom &face2,
131 bool doRot, int dir,
132 NekDouble angle,
133 NekDouble tol)
134{
135 return GetFaceOrientation(face1.m_verts, face2.m_verts, doRot, dir, angle,
136 tol);
137}
138
139/**
140 * Calculate the orientation of face2 to face1 (note this is
141 * not face1 to face2!).
142 */
144 const PointGeomVector &face1, const PointGeomVector &face2, bool doRot,
145 int dir, NekDouble angle, NekDouble tol)
146{
147 int i, j, vmap[4] = {-1, -1, -1, -1};
148
149 if (doRot)
150 {
151 PointGeom rotPt;
152
153 for (i = 0; i < 4; ++i)
154 {
155 rotPt.Rotate((*face1[i]), dir, angle);
156 for (j = 0; j < 4; ++j)
157 {
158 if (rotPt.dist(*face2[j]) < tol)
159 {
160 vmap[j] = i;
161 break;
162 }
163 }
164 }
165 }
166 else
167 {
168
169 NekDouble x, y, z, x1, y1, z1, cx = 0.0, cy = 0.0, cz = 0.0;
170
171 // For periodic faces, we calculate the vector between the centre
172 // points of the two faces. (For connected faces this will be
173 // zero). We can then use this to determine alignment later in the
174 // algorithm.
175 for (i = 0; i < 4; ++i)
176 {
177 cx += (*face2[i])(0) - (*face1[i])(0);
178 cy += (*face2[i])(1) - (*face1[i])(1);
179 cz += (*face2[i])(2) - (*face1[i])(2);
180 }
181 cx /= 4;
182 cy /= 4;
183 cz /= 4;
184
185 // Now construct a mapping which takes us from the vertices of one
186 // face to the other. That is, vertex j of face2 corresponds to
187 // vertex vmap[j] of face1.
188 for (i = 0; i < 4; ++i)
189 {
190 x = (*face1[i])(0);
191 y = (*face1[i])(1);
192 z = (*face1[i])(2);
193 for (j = 0; j < 4; ++j)
194 {
195 x1 = (*face2[j])(0) - cx;
196 y1 = (*face2[j])(1) - cy;
197 z1 = (*face2[j])(2) - cz;
198 if (sqrt((x1 - x) * (x1 - x) + (y1 - y) * (y1 - y) +
199 (z1 - z) * (z1 - z)) < 1e-8)
200 {
201 vmap[j] = i;
202 break;
203 }
204 }
205 }
206 }
207
208 // Use the mapping to determine the eight alignment options between
209 // faces.
210 if (vmap[1] == (vmap[0] + 1) % 4)
211 {
212 switch (vmap[0])
213 {
214 case 0:
216 break;
217 case 1:
219 break;
220 case 2:
222 break;
223 case 3:
225 break;
226 }
227 }
228 else
229 {
230 switch (vmap[0])
231 {
232 case 0:
234 break;
235 case 1:
237 break;
238 case 2:
240 break;
241 case 3:
243 break;
244 }
245 }
246
247 ASSERTL0(false, "unable to determine face orientation");
249}
250
251/**
252 * Set up GeoFac for this geometry using Coord quadrature distribution
253 */
255{
256 if (!m_setupState)
257 {
259 }
260
262 {
263 GeomType Gtype = eRegular;
264
266
267 // We will first check whether we have a regular or deformed
268 // geometry. We will define regular as those cases where the
269 // Jacobian and the metric terms of the derivative are constants
270 // (i.e. not coordinate dependent)
271
272 // Check to see if expansions are linear
273 // If not linear => deformed geometry
274 m_straightEdge = 1;
275 if ((m_xmap->GetBasisNumModes(0) != 2) ||
276 (m_xmap->GetBasisNumModes(1) != 2))
277 {
278 Gtype = eDeformed;
279 m_straightEdge = 0;
280 }
281
282 // For linear expansions, the mapping from standard to local
283 // element is given by the relation:
284 // x_i = 0.25 * [ ( x_i^A + x_i^B + x_i^C + x_i^D) +
285 // (-x_i^A + x_i^B + x_i^C - x_i^D)*xi_1 +
286 // (-x_i^A - x_i^B + x_i^C + x_i^D)*xi_2 +
287 // ( x_i^A - x_i^B + x_i^C - x_i^D)*xi_1*xi_2 ]
288 //
289 // The jacobian of the transformation and the metric terms
290 // dxi_i/dx_j, involve only terms of the form dx_i/dxi_j (both
291 // for coordim == 2 or 3). Inspecting the formula above, it can
292 // be appreciated that the derivatives dx_i/dxi_j will be
293 // constant, if the coefficient of the non-linear term is zero.
294 //
295 // That is why for regular geometry, we require
296 //
297 // x_i^A - x_i^B + x_i^C - x_i^D = 0
298 //
299 // or equivalently
300 //
301 // x_i^A - x_i^B = x_i^D - x_i^C
302 //
303 // This corresponds to quadrilaterals which are paralellograms.
305 m_manifold[0] = 0;
306 m_manifold[1] = 1;
307 if (m_coordim == 3)
308 {
309 PointGeom e01, e21, norm;
310 e01.Sub(*m_verts[0], *m_verts[1]);
311 e21.Sub(*m_verts[3], *m_verts[1]);
312 norm.Mult(e01, e21);
313 int tmpi = 0;
314 double tmp = std::fabs(norm[0]);
315 if (tmp < fabs(norm[1]))
316 {
317 tmp = fabs(norm[1]);
318 tmpi = 1;
319 }
320 if (tmp < fabs(norm[2]))
321 {
322 tmpi = 2;
323 }
324 m_manifold[0] = (tmpi + 1) % 3;
325 m_manifold[1] = (tmpi + 2) % 3;
326 m_manifold[2] = (tmpi + 3) % 3;
327 }
328
329 if (Gtype == eRegular)
330 {
332 for (int i = 0; i < m_verts.size(); ++i)
333 {
334 verts[i] = Array<OneD, NekDouble>(3);
335 m_verts[i]->GetCoords(verts[i]);
336 }
337 // a00 + a01 xi1 + a02 xi2 + a03 xi1 xi2
338 // a10 + a11 xi1 + a12 xi2 + a03 xi1 xi2
340 for (int i = 0; i < 2; i++)
341 {
342 unsigned int d = m_manifold[i];
344 // Karniadakis, Sherwin 2005, Appendix D
345 NekDouble A = verts[0][d];
346 NekDouble B = verts[1][d];
347 NekDouble D = verts[2][d];
348 NekDouble C = verts[3][d];
349 m_isoParameter[i][0] = 0.25 * (A + B + C + D); // 1
350 m_isoParameter[i][1] = 0.25 * (-A + B - C + D); // xi1
351 m_isoParameter[i][2] = 0.25 * (-A - B + C + D); // xi2
352 m_isoParameter[i][3] = 0.25 * (A - B - C + D); // xi1*xi2
353 NekDouble tmp =
354 fabs(m_isoParameter[i][1]) + fabs(m_isoParameter[i][2]);
355 if (fabs(m_isoParameter[i][3]) >
357 {
358 Gtype = eDeformed;
359 }
360 }
361 }
362
363 if (Gtype == eRegular)
364 {
366 }
367 else if (m_straightEdge)
368 {
370 }
371
373 Gtype, m_coordim, m_xmap, m_coeffs);
375 }
376}
377
378/**
379 * Note verts and edges are listed according to anticlockwise
380 * convention but points in _coeffs have to be in array format from
381 * left to right.
382 */
384{
385 // check to see if geometry structure is already filled
386 if (m_state != ePtsFilled)
387 {
388 int i, j, k;
389 int nEdgeCoeffs;
390
391 if (m_curve)
392 {
393 int npts = m_curve->m_points.size();
394 int nEdgePts = (int)sqrt(static_cast<NekDouble>(npts));
395 Array<OneD, NekDouble> tmp(npts);
396 Array<OneD, NekDouble> tmp2(m_xmap->GetTotPoints());
397 LibUtilities::PointsKey curveKey(nEdgePts, m_curve->m_ptype);
398
399 // Sanity checks:
400 // - Curved faces should have square number of points;
401 // - Each edge should have sqrt(npts) points.
402 ASSERTL0(nEdgePts * nEdgePts == npts,
403 "NUMPOINTS should be a square number in"
404 " quadrilteral " +
405 std::to_string(m_globalID));
406
407 for (i = 0; i < kNedges; ++i)
408 {
409 ASSERTL0(m_edges[i]->GetXmap()->GetNcoeffs() == nEdgePts,
410 "Number of edge points does not correspond to "
411 "number of face points in quadrilateral " +
412 std::to_string(m_globalID));
413 }
414
415 for (i = 0; i < m_coordim; ++i)
416 {
417 for (j = 0; j < npts; ++j)
418 {
419 tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
420 }
421
422 // Interpolate m_curve points to GLL points
423 LibUtilities::Interp2D(curveKey, curveKey, tmp,
424 m_xmap->GetBasis(0)->GetPointsKey(),
425 m_xmap->GetBasis(1)->GetPointsKey(),
426 tmp2);
427
428 // Forwards transform to get coefficient space.
429 m_xmap->FwdTrans(tmp2, m_coeffs[i]);
430 }
431 }
432
433 // Now fill in edges.
435 Array<OneD, int> signArray;
436
437 for (i = 0; i < kNedges; i++)
438 {
439 m_edges[i]->FillGeom();
440 m_xmap->GetTraceToElementMap(i, mapArray, signArray, m_eorient[i]);
441
442 nEdgeCoeffs = m_edges[i]->GetXmap()->GetNcoeffs();
443
444 for (j = 0; j < m_coordim; j++)
445 {
446 for (k = 0; k < nEdgeCoeffs; k++)
447 {
448 m_coeffs[j][mapArray[k]] =
449 signArray[k] * (m_edges[i]->GetCoeffs(j))[k];
450 }
451 }
452 }
453
455 }
456}
457
459{
460 int i0, i1, j1, j2;
461 if (fabs(m_isoParameter[0][3]) >= fabs(m_isoParameter[1][3]))
462 {
463 i0 = 0;
464 i1 = 1;
465 }
466 else
467 {
468 i1 = 0;
469 i0 = 1;
470 m_straightEdge |= 2;
471 }
472 NekDouble gamma = m_isoParameter[i1][3] / m_isoParameter[i0][3];
473 std::vector<NekDouble> c(3);
474 for (int i = 0; i < 3; ++i)
475 {
476 c[i] = m_isoParameter[i1][i] - gamma * m_isoParameter[i0][i];
477 }
478 if (fabs(c[1]) >= fabs(c[2]))
479 {
480 j1 = 1;
481 j2 = 2;
482 }
483 else
484 {
485 j1 = 2;
486 j2 = 1;
487 m_straightEdge |= 4;
488 }
489 NekDouble beta = c[j2] / c[j1];
490 if (i0 == 1)
491 {
493 }
494 if (j1 == 2)
495 {
496 NekDouble temp = m_isoParameter[0][j1];
497 m_isoParameter[0][j1] = m_isoParameter[0][j2];
498 m_isoParameter[0][j2] = temp;
499 }
500 m_isoParameter[0][2] -= m_isoParameter[0][1] * beta;
501 m_isoParameter[1][0] = c[0];
502 m_isoParameter[1][1] = 1. / c[j1];
503 m_isoParameter[1][2] = beta;
504 m_isoParameter[1][3] = gamma;
505}
506
507int QuadGeom::v_GetDir(const int i, [[maybe_unused]] const int j) const
508{
509 return i % 2;
510}
511
512void QuadGeom::v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
513{
514 Geometry::v_Reset(curvedEdges, curvedFaces);
515 CurveMap::iterator it = curvedFaces.find(m_globalID);
516
517 if (it != curvedFaces.end())
518 {
519 m_curve = it->second;
520 }
521
522 for (int i = 0; i < 4; ++i)
523 {
524 m_edges[i]->Reset(curvedEdges, curvedFaces);
525 }
526
527 SetUpXmap();
528 SetUpCoeffs(m_xmap->GetNcoeffs());
529}
530
532{
533 if (!m_setupState)
534 {
535 for (int i = 0; i < 4; ++i)
536 {
537 m_edges[i]->Setup();
538 }
539 SetUpXmap();
540 SetUpCoeffs(m_xmap->GetNcoeffs());
541 m_setupState = true;
542 }
543}
544
545} // namespace Nektar::SpatialDomains
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
Describes the specification for a Basis.
Definition: Basis.h:45
Defines a specification for a set of points.
Definition: Points.h:50
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
2D geometry information
Definition: Geometry2D.h:65
Array< OneD, int > m_manifold
Definition: Geometry2D.h:83
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry2D.h:81
void v_CalculateInverseIsoParam() override
Definition: Geometry2D.cpp:629
bool m_setupState
Wether or not the setup routines have been run.
Definition: Geometry.h:198
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
Definition: Geometry.h:357
GeomState m_state
Enumeration to dictate whether coefficients are filled.
Definition: Geometry.h:196
void SetUpCoeffs(const int nCoeffs)
Initialise the Geometry::m_coeffs array.
Definition: Geometry.h:700
Array< OneD, Array< OneD, NekDouble > > m_isoParameter
Definition: Geometry.h:209
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated Geo...
Definition: Geometry.cpp:364
LibUtilities::ShapeType m_shapeType
Type of shape.
Definition: Geometry.h:202
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:206
GeomState m_geomFactorsState
State of the geometric factors.
Definition: Geometry.h:192
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:194
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
Definition: Geometry.h:435
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:190
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:188
void Sub(PointGeom &a, PointGeom &b)
Definition: PointGeom.cpp:119
void Mult(PointGeom &a, PointGeom &b)
_this = a x b
Definition: PointGeom.cpp:128
void Rotate(PointGeom &a, int dir, NekDouble angle)
_this = rotation of a by angle 'angle' around axis dir
Definition: PointGeom.cpp:137
NekDouble dist(PointGeom &a)
return distance between this and input a
Definition: PointGeom.cpp:175
void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces) override
Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated Geo...
Definition: QuadGeom.cpp:512
static StdRegions::Orientation GetFaceOrientation(const QuadGeom &face1, const QuadGeom &face2, bool doRot=false, int dir=0, NekDouble angle=0.0, NekDouble tol=1e-8)
Get the orientation of face1.
Definition: QuadGeom.cpp:129
NekDouble v_GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord) override
Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i.
Definition: QuadGeom.cpp:118
void v_GenGeomFactors() override
Definition: QuadGeom.cpp:254
static const int kNedges
Definition: QuadGeom.h:75
int v_GetDir(const int faceidx, const int facedir) const override
Returns the element coordinate direction corresponding to a given face coordinate direction.
Definition: QuadGeom.cpp:507
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:190
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
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
static const NekDouble kNekZeroTol
std::vector< PointGeomSharedPtr > PointGeomVector
Definition: Geometry2D.h:60
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:58
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:59
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:58
@ ePtsFilled
Geometric information has been generated.
std::vector< double > z(NPUPPER)
std::vector< double > d(NPUPPER *NPUPPER)
double NekDouble
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285