Nektar++
SegGeom.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: SegGeom.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
39#include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
42
44{
46{
48}
49
50SegGeom::SegGeom(int id, const int coordim, const PointGeomSharedPtr vertex[],
51 const CurveSharedPtr curve)
52 : Geometry1D(coordim)
53{
55 m_globalID = id;
57 m_curve = curve;
58
59 m_verts[0] = vertex[0];
60 m_verts[1] = vertex[1];
61}
62
64{
65 // From Geometry class
67
68 // info from EdgeComponent class
70 m_xmap = in.m_xmap;
71 SetUpCoeffs(m_xmap->GetNcoeffs());
72
73 // info from SegGeom class
75 m_verts[0] = in.m_verts[0];
76 m_verts[1] = in.m_verts[1];
77
78 m_state = in.m_state;
79}
80
82{
83 if (m_curve)
84 {
85 int npts = m_curve->m_points.size();
86 LibUtilities::PointsKey pkey(npts + 1,
90 }
91 else
92 {
97 }
98}
99
100/**
101 * \brief Generate a one dimensional space segment geometry where the vert[0]
102 * has the same x value and vert[1] is set to vert[0] plus the length of the
103 * original segment
104 **/
106{
108
109 // info about numbering
110 returnval->m_globalID = m_globalID;
111
112 // geometric information.
113 returnval->m_coordim = 1;
114 NekDouble x0 = (*m_verts[0])[0];
116 1, m_verts[0]->GetGlobalID(), x0, 0.0, 0.0);
117 vert0->SetGlobalID(vert0->GetGlobalID());
118 returnval->m_verts[0] = vert0;
119
120 // Get information to calculate length.
122 m_xmap->GetBase();
124 v.push_back(base[0]->GetPointsKey());
126
127 const Array<OneD, const NekDouble> jac = m_geomFactors->GetJac(v);
128
129 NekDouble len = 0.0;
130 if (jac.size() == 1)
131 {
132 len = jac[0] * 2.0;
133 }
134 else
135 {
136 Array<OneD, const NekDouble> w0 = base[0]->GetW();
137 len = 0.0;
138
139 for (int i = 0; i < jac.size(); ++i)
140 {
141 len += jac[i] * w0[i];
142 }
143 }
144 // Set up second vertex.
146 1, m_verts[1]->GetGlobalID(), x0 + len, 0.0, 0.0);
147 vert1->SetGlobalID(vert1->GetGlobalID());
148
149 returnval->m_verts[1] = vert1;
150
151 // at present just use previous m_xmap[0];
152 returnval->m_xmap = m_xmap;
153 returnval->SetUpCoeffs(m_xmap->GetNcoeffs());
154 returnval->m_state = eNotFilled;
155
156 return returnval;
157}
158
160{
162}
163
165 const Array<OneD, const NekDouble> &Lcoord)
166{
167 ASSERTL1(m_state == ePtsFilled, "Geometry is not in physical space");
168
169 Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
170 m_xmap->BwdTrans(m_coeffs[i], tmp);
171
172 return m_xmap->PhysEvaluate(Lcoord, tmp);
173}
174
175/**
176 * @brief Get the orientation of @p edge1.
177 *
178 * If @p edge1 is connected to @p edge2 in the same direction as the points
179 * comprising @p edge1 then it is forward, otherwise it is backward.
180 *
181 * For example, assume @p edge1 is comprised of points 1 and 2, and @p edge2 is
182 * comprised of points 2 and 3, then @p edge1 is forward.
183 *
184 * If @p edge1 is comprised of points 2 and 1 and @p edge2 is comprised of
185 * points 3 and 2, then @p edge1 is backward.
186 *
187 * Since both edges are passed, it does not need any information from the
188 * EdgeComponent instance.
189 */
191 const SegGeom &edge2)
192{
194
195 if ((*edge1.GetVertex(0) == *edge2.GetVertex(0)) ||
196 (*edge1.GetVertex(0) == *edge2.GetVertex(1)))
197 {
198 // Backward direction. Vertex 0 is connected to edge 2.
199 returnval = StdRegions::eBackwards;
200 }
201 else if ((*edge1.GetVertex(1) != *edge2.GetVertex(0)) &&
202 (*edge1.GetVertex(1) != *edge2.GetVertex(1)))
203 {
204 // Not forward either, then we have a problem.
205 std::ostringstream errstrm;
206 errstrm << "Connected edges do not share a vertex. Edges ";
207 errstrm << edge1.GetGlobalID() << ", " << edge2.GetGlobalID();
208 ASSERTL0(false, errstrm.str());
209 }
210
211 return returnval;
212}
213
215{
216 if (!m_setupState)
217 {
219 }
220
222 {
225
226 if (m_xmap->GetBasisNumModes(0) != 2)
227 {
228 gType = eDeformed;
229 }
230
232 gType, m_coordim, m_xmap, m_coeffs);
234 }
235}
236
238{
239 if (m_state != ePtsFilled)
240 {
241 int i;
242
243 if (m_coordim > 0 && m_curve)
244 {
245 int npts = m_curve->m_points.size();
246 LibUtilities::PointsKey pkey(npts + 1,
248 Array<OneD, NekDouble> tmp(npts);
249
250 if (m_verts[0]->dist(*(m_curve->m_points[0])) >
252 {
253 std::string err =
254 "Vertex 0 is separated from first point by more than ";
255 std::stringstream strstrm;
256 strstrm << NekConstants::kVertexTheSameDouble << " in edge "
257 << m_globalID;
258 err += strstrm.str();
259 NEKERROR(ErrorUtil::ewarning, err.c_str());
260 }
261
262 if (m_verts[1]->dist(*(m_curve->m_points[npts - 1])) >
264 {
265 std::string err =
266 "Vertex 1 is separated from last point by more than ";
267 std::stringstream strstrm;
268 strstrm << NekConstants::kVertexTheSameDouble << " in edge "
269 << m_globalID;
270 err += strstrm.str();
271 NEKERROR(ErrorUtil::ewarning, err.c_str());
272 }
273
274 LibUtilities::PointsKey fkey(npts, m_curve->m_ptype);
276 LibUtilities::PointsManager()[fkey]->GetI(pkey);
277 NekVector<NekDouble> out(npts + 1);
278
279 for (int i = 0; i < m_coordim; ++i)
280 {
281 // Load up coordinate values into tmp
282 for (int j = 0; j < npts; ++j)
283 {
284 tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
285 }
286
287 // Interpolate to GLL points
288 NekVector<NekDouble> in(npts, tmp, eWrapper);
289 out = (*I0) * in;
290
291 m_xmap->FwdTrans(out.GetPtr(), m_coeffs[i]);
292 }
293 }
294
295 for (i = 0; i < m_coordim; ++i)
296 {
297 m_coeffs[i][0] = (*m_verts[0])[i];
298 m_coeffs[i][1] = (*m_verts[1])[i];
299 }
300
302 }
303}
304
305void SegGeom::v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
306{
307 Geometry::v_Reset(curvedEdges, curvedFaces);
308 CurveMap::iterator it = curvedEdges.find(m_globalID);
309
310 if (it != curvedEdges.end())
311 {
312 m_curve = it->second;
313 }
314
315 SetUpXmap();
316 SetUpCoeffs(m_xmap->GetNcoeffs());
317}
318
320{
321 if (!m_setupState)
322 {
323 SetUpXmap();
324 SetUpCoeffs(m_xmap->GetNcoeffs());
325 m_setupState = true;
326 }
327}
328
330{
331 PointGeomSharedPtr returnval;
332
333 if (i >= 0 && i < kNverts)
334 {
335 returnval = m_verts[i];
336 }
337
338 return returnval;
339}
340
342{
343 return kNverts;
344}
345
348{
349 if (m_geomFactors->GetGtype() == eRegular)
350 {
351 xiOut = Array<OneD, NekDouble>(1, 0.0);
352
353 GetLocCoords(xs, xiOut);
354 ClampLocCoords(xiOut);
355
357 NekDouble tmp = 0;
358 for (int i = 0; i < m_coordim; ++i)
359 {
360 gloCoord[i] = GetCoord(i, xiOut);
361 tmp += (xs[i] - gloCoord[i]) * (xs[i] - gloCoord[i]);
362 }
363
364 return sqrt(tmp);
365 }
366 // If deformed edge then the inverse mapping is non-linear so need to
367 // numerically solve for the local coordinate
368 else if (m_geomFactors->GetGtype() == eDeformed)
369 {
370 Array<OneD, NekDouble> xi(1, 0.0);
371
372 // Armijo constants:
373 // https://en.wikipedia.org/wiki/Backtracking_line_search
374 const NekDouble c1 = 1e-4, c2 = 0.9;
375
376 int dim = GetCoordim();
377 int nq = m_xmap->GetTotPoints();
378
379 Array<OneD, Array<OneD, NekDouble>> x(dim), xder(dim), xder2(dim);
380 // Get x,y,z phys values from coefficients
381 for (int i = 0; i < dim; ++i)
382 {
383 x[i] = Array<OneD, NekDouble>(nq);
384 xder[i] = Array<OneD, NekDouble>(nq);
385 xder2[i] = Array<OneD, NekDouble>(nq);
386
387 m_xmap->BwdTrans(m_coeffs[i], x[i]);
388 }
389
390 NekDouble fx_prev = std::numeric_limits<NekDouble>::max();
391
392 // Minimisation loop (Quasi-newton method)
393 for (int i = 0; i < NekConstants::kNewtonIterations; ++i)
394 {
395 // Compute the objective function, f(x_k) and its derivatives
399 NekDouble fx = 0, fxp = 0, fxp2 = 0, xcDiff = 0;
400 for (int j = 0; j < dim; ++j)
401 {
402 xc[j] = m_xmap->PhysEvaluate(xi, x[j], xc_der[j], xc_der2[j]);
403
404 xcDiff = xc[j] - xs[j];
405 // Objective function is the distance to the search point
406 fx += xcDiff * xcDiff;
407 fxp += xc_der[j][0] * xcDiff;
408 fxp2 += xc_der2[j][0] * xcDiff + xc_der[j][0] * xc_der[j][0];
409 }
410
411 fxp *= 2;
412 fxp2 *= 2;
413
414 // Check for convergence
415 if (std::abs(fx - fx_prev) < 1e-12)
416 {
417 fx_prev = fx;
418 break;
419 }
420 else
421 {
422 fx_prev = fx;
423 }
424
425 NekDouble gamma = 1.0;
426 bool conv = false;
427
428 // Search direction: Newton's method
429 NekDouble pk = -fxp / fxp2;
430
431 // Perform backtracking line search
432 while (gamma > 1e-10)
433 {
434 Array<OneD, NekDouble> xi_pk(1);
435 xi_pk[0] = xi[0] + pk * gamma;
436
437 if (xi_pk[0] < -1.0 || xi_pk[0] > 1.0)
438 {
439 gamma /= 2.0;
440 continue;
441 }
442
443 Array<OneD, NekDouble> xc_pk(dim);
445 NekDouble fx_pk = 0, fxp_pk = 0, xc_pkDiff = 0;
446 for (int j = 0; j < dim; ++j)
447 {
448 xc_pk[j] = m_xmap->PhysEvaluate(xi_pk, x[j], xc_der_pk[j]);
449
450 xc_pkDiff = xc_pk[j] - xs[j];
451 fx_pk += xc_pkDiff * xc_pkDiff;
452 fxp_pk += xc_der_pk[j][0] * xc_pkDiff;
453 }
454
455 fxp_pk *= 2;
456
457 // Check Wolfe conditions using Armijo constants
458 // https://en.wikipedia.org/wiki/Wolfe_conditions
459 if ((fx_pk - (fx + c1 * gamma * pk * fxp)) <
460 std::numeric_limits<NekDouble>::epsilon() &&
461 (-pk * fxp_pk + c2 * pk * fxp) <
462 std::numeric_limits<NekDouble>::epsilon())
463 {
464 conv = true;
465 break;
466 }
467
468 gamma /= 2.0;
469 }
470
471 if (!conv)
472 {
473 break;
474 }
475
476 xi[0] += gamma * pk;
477 }
478
479 xiOut = xi;
480 return sqrt(fx_prev);
481 }
482 else
483 {
484 ASSERTL0(false, "Geometry type unknown")
485 }
486
487 return -1.0;
488}
489
490} // namespace Nektar::SpatialDomains
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#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.
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:217
1D geometry information
Definition: Geometry1D.h:54
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.
Definition: Geometry.h:554
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
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 ge...
Definition: Geometry.h:544
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
bool ClampLocCoords(Array< OneD, NekDouble > &locCoord, NekDouble tol=std::numeric_limits< NekDouble >::epsilon())
Clamp local coords to be within standard regions [-1, 1]^dim.
Definition: Geometry.cpp:530
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:326
int GetCoordim() const
Return the coordinate dimension of this object (i.e. the dimension of the space in which this object ...
Definition: Geometry.h:283
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
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:190
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:188
PointGeomSharedPtr v_GetVertex(const int i) const override
Definition: SegGeom.cpp:329
int v_GetNumVerts() const override
Get the number of vertices of this object.
Definition: SegGeom.cpp:341
virtual LibUtilities::ShapeType v_GetShapeType() const
Definition: SegGeom.cpp:159
void v_FillGeom() override
Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements.
Definition: SegGeom.cpp:237
SegGeomSharedPtr GenerateOneSpaceDimGeom(void)
Generate a one dimensional space segment geometry where the vert[0] has the same x value and vert[1] ...
Definition: SegGeom.cpp:105
CurveSharedPtr m_curve
Boolean indicating whether object owns the data.
Definition: SegGeom.h:94
void v_GenGeomFactors() override
Definition: SegGeom.cpp:214
SpatialDomains::PointGeomSharedPtr m_verts[kNverts]
Definition: SegGeom.h:77
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:190
NekDouble v_FindDistance(const Array< OneD, const NekDouble > &xs, Array< OneD, NekDouble > &xi) override
Definition: SegGeom.cpp:346
static const int kNverts
Definition: SegGeom.h:74
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: SegGeom.cpp:164
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: SegGeom.cpp:305
PointsManagerT & PointsManager(void)
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
static const NekDouble kVertexTheSameDouble
static const unsigned int kNewtonIterations
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
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
@ eNotFilled
Geometric information has not been generated.
@ ePtsFilled
Geometric information has been generated.
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285