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 
37 #include <SpatialDomains/SegGeom.h>
38 
39 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
41 #include <StdRegions/StdSegExp.h>
42 
43 namespace Nektar
44 {
45 namespace SpatialDomains
46 {
48 {
50 }
51 
52 SegGeom::SegGeom(int id, const int coordim, const PointGeomSharedPtr vertex[],
53  const CurveSharedPtr curve)
54  : Geometry1D(coordim)
55 {
57  m_globalID = id;
59  m_curve = curve;
60 
61  m_verts[0] = vertex[0];
62  m_verts[1] = vertex[1];
63 }
64 
66 {
67  // From Geometry class
69 
70  // info from EdgeComponent class
72  m_xmap = in.m_xmap;
73  SetUpCoeffs(m_xmap->GetNcoeffs());
74 
75  // info from SegGeom class
76  m_coordim = in.m_coordim;
77  m_verts[0] = in.m_verts[0];
78  m_verts[1] = in.m_verts[1];
79 
80  m_state = in.m_state;
81 }
82 
84 {
85  if (m_curve)
86  {
87  int npts = m_curve->m_points.size();
88  LibUtilities::PointsKey pkey(npts + 1,
92  }
93  else
94  {
95  const LibUtilities::BasisKey B(
99  }
100 }
101 
102 /**
103  * \brief Generate a one dimensional space segment geometry where the vert[0]
104  * has the same x value and vert[1] is set to vert[0] plus the length of the
105  * original segment
106  **/
108 {
110 
111  // info about numbering
112  returnval->m_globalID = m_globalID;
113 
114  // geometric information.
115  returnval->m_coordim = 1;
116  NekDouble x0 = (*m_verts[0])[0];
118  1, m_verts[0]->GetGlobalID(), x0, 0.0, 0.0);
119  vert0->SetGlobalID(vert0->GetGlobalID());
120  returnval->m_verts[0] = vert0;
121 
122  // Get information to calculate length.
124  m_xmap->GetBase();
126  v.push_back(base[0]->GetPointsKey());
128 
129  const Array<OneD, const NekDouble> jac = m_geomFactors->GetJac(v);
130 
131  NekDouble len = 0.0;
132  if (jac.size() == 1)
133  {
134  len = jac[0] * 2.0;
135  }
136  else
137  {
138  Array<OneD, const NekDouble> w0 = base[0]->GetW();
139  len = 0.0;
140 
141  for (int i = 0; i < jac.size(); ++i)
142  {
143  len += jac[i] * w0[i];
144  }
145  }
146  // Set up second vertex.
148  1, m_verts[1]->GetGlobalID(), x0 + len, 0.0, 0.0);
149  vert1->SetGlobalID(vert1->GetGlobalID());
150 
151  returnval->m_verts[1] = vert1;
152 
153  // at present just use previous m_xmap[0];
154  returnval->m_xmap = m_xmap;
155  returnval->SetUpCoeffs(m_xmap->GetNcoeffs());
156  returnval->m_state = eNotFilled;
157 
158  return returnval;
159 }
160 
162 {
163 }
164 
166 {
167  return LibUtilities::eSegment;
168 }
169 
171  const Array<OneD, const NekDouble> &Lcoord)
172 {
173  ASSERTL1(m_state == ePtsFilled, "Geometry is not in physical space");
174 
175  Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
176  m_xmap->BwdTrans(m_coeffs[i], tmp);
177 
178  return m_xmap->PhysEvaluate(Lcoord, tmp);
179 }
180 
181 /**
182  * @brief Get the orientation of @p edge1.
183  *
184  * If @p edge1 is connected to @p edge2 in the same direction as the points
185  * comprising @p edge1 then it is forward, otherwise it is backward.
186  *
187  * For example, assume @p edge1 is comprised of points 1 and 2, and @p edge2 is
188  * comprised of points 2 and 3, then @p edge1 is forward.
189  *
190  * If @p edge1 is comprised of points 2 and 1 and @p edge2 is comprised of
191  * points 3 and 2, then @p edge1 is backward.
192  *
193  * Since both edges are passed, it does not need any information from the
194  * EdgeComponent instance.
195  */
197  const SegGeom &edge2)
198 {
200 
201  if ((*edge1.GetVertex(0) == *edge2.GetVertex(0)) ||
202  (*edge1.GetVertex(0) == *edge2.GetVertex(1)))
203  {
204  // Backward direction. Vertex 0 is connected to edge 2.
205  returnval = StdRegions::eBackwards;
206  }
207  else if ((*edge1.GetVertex(1) != *edge2.GetVertex(0)) &&
208  (*edge1.GetVertex(1) != *edge2.GetVertex(1)))
209  {
210  // Not forward either, then we have a problem.
211  std::ostringstream errstrm;
212  errstrm << "Connected edges do not share a vertex. Edges ";
213  errstrm << edge1.GetGlobalID() << ", " << edge2.GetGlobalID();
214  ASSERTL0(false, errstrm.str());
215  }
216 
217  return returnval;
218 }
219 
221 {
222  if (!m_setupState)
223  {
225  }
226 
228  {
231 
232  if (m_xmap->GetBasisNumModes(0) != 2)
233  {
234  gType = eDeformed;
235  }
236 
238  gType, m_coordim, m_xmap, m_coeffs);
240  }
241 }
242 
244 {
245  if (m_state != ePtsFilled)
246  {
247  int i;
248 
249  if (m_coordim > 0 && m_curve)
250  {
251  int npts = m_curve->m_points.size();
252  LibUtilities::PointsKey pkey(npts + 1,
254  Array<OneD, NekDouble> tmp(npts);
255 
256  if (m_verts[0]->dist(*(m_curve->m_points[0])) >
258  {
259  std::string err =
260  "Vertex 0 is separated from first point by more than ";
261  std::stringstream strstrm;
262  strstrm << NekConstants::kVertexTheSameDouble << " in edge "
263  << m_globalID;
264  err += strstrm.str();
265  NEKERROR(ErrorUtil::ewarning, err.c_str());
266  }
267 
268  if (m_verts[1]->dist(*(m_curve->m_points[npts - 1])) >
270  {
271  std::string err =
272  "Vertex 1 is separated from last point by more than ";
273  std::stringstream strstrm;
274  strstrm << NekConstants::kVertexTheSameDouble << " in edge "
275  << m_globalID;
276  err += strstrm.str();
277  NEKERROR(ErrorUtil::ewarning, err.c_str());
278  }
279 
280  LibUtilities::PointsKey fkey(npts, m_curve->m_ptype);
281  DNekMatSharedPtr I0 =
282  LibUtilities::PointsManager()[fkey]->GetI(pkey);
283  NekVector<NekDouble> out(npts + 1);
284 
285  for (int i = 0; i < m_coordim; ++i)
286  {
287  // Load up coordinate values into tmp
288  for (int j = 0; j < npts; ++j)
289  {
290  tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
291  }
292 
293  // Interpolate to GLL points
294  NekVector<NekDouble> in(npts, tmp, eWrapper);
295  out = (*I0) * in;
296 
297  m_xmap->FwdTrans(out.GetPtr(), m_coeffs[i]);
298  }
299  }
300 
301  for (i = 0; i < m_coordim; ++i)
302  {
303  m_coeffs[i][0] = (*m_verts[0])[i];
304  m_coeffs[i][1] = (*m_verts[1])[i];
305  }
306 
308  }
309 }
310 
311 void SegGeom::v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
312 {
313  Geometry::v_Reset(curvedEdges, curvedFaces);
314  CurveMap::iterator it = curvedEdges.find(m_globalID);
315 
316  if (it != curvedEdges.end())
317  {
318  m_curve = it->second;
319  }
320 
321  SetUpXmap();
322  SetUpCoeffs(m_xmap->GetNcoeffs());
323 }
324 
326 {
327  if (!m_setupState)
328  {
329  SetUpXmap();
330  SetUpCoeffs(m_xmap->GetNcoeffs());
331  m_setupState = true;
332  }
333 }
334 
336 {
337  PointGeomSharedPtr returnval;
338 
339  if (i >= 0 && i < kNverts)
340  {
341  returnval = m_verts[i];
342  }
343 
344  return returnval;
345 }
346 
348 {
349  return kNverts;
350 }
351 
352 } // namespace SpatialDomains
353 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
Describes the specification for a Basis.
Definition: Basis.h:50
Defines a specification for a set of points.
Definition: Points.h:59
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:55
bool m_setupState
Wether or not the setup routines have been run.
Definition: Geometry.h:190
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
Definition: Geometry.h:345
GeomState m_state
Enumeration to dictate whether coefficients are filled.
Definition: Geometry.h:188
void SetUpCoeffs(const int nCoeffs)
Initialise the Geometry::m_coeffs array.
Definition: Geometry.h:671
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:367
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:314
LibUtilities::ShapeType m_shapeType
Type of shape.
Definition: Geometry.h:194
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:198
GeomState m_geomFactorsState
State of the geometric factors.
Definition: Geometry.h:184
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:186
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:182
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:180
virtual int v_GetNumVerts() const
Get the number of vertices of this object.
Definition: SegGeom.cpp:347
virtual PointGeomSharedPtr v_GetVertex(const int i) const
Definition: SegGeom.cpp:335
virtual LibUtilities::ShapeType v_GetShapeType() const
Definition: SegGeom.cpp:165
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: SegGeom.cpp:311
virtual void v_GenGeomFactors()
Definition: SegGeom.cpp:220
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:107
CurveSharedPtr m_curve
Boolean indicating whether object owns the data.
Definition: SegGeom.h:95
SpatialDomains::PointGeomSharedPtr m_verts[kNverts]
Definition: SegGeom.h:79
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:196
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.
Definition: SegGeom.cpp:170
virtual void v_FillGeom()
Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements.
Definition: SegGeom.cpp:243
static const int kNverts
Definition: SegGeom.h:76
PointsManagerT & PointsManager(void)
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
static const NekDouble kVertexTheSameDouble
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:60
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:61
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:62
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:59
@ eNotFilled
Geometric information has not been generated.
@ ePtsFilled
Geometric information has been generated.
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble