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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description:
33 //
34 //
35 ////////////////////////////////////////////////////////////////////////////////
36 
37 #include <SpatialDomains/SegGeom.h>
39 
41 #include <StdRegions/StdSegExp.h>
42 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
43 
44 namespace Nektar
45 {
46  namespace SpatialDomains
47  {
49  {
51  }
52 
53  SegGeom::SegGeom(int id, const int coordim):
54  Geometry1D(coordim)
55  {
56 
58  m_eid = id;
59  m_globalID = id;
60  SetUpXmap();
61  SetUpCoeffs(m_xmap->GetNcoeffs());
62  }
63 
65  int id,
66  const int coordim,
67  const PointGeomSharedPtr vertex[]):
68  Geometry1D(coordim)
69  {
71  m_eid = id;
72  m_globalID = id;
74 
75  if (coordim > 0)
76  {
77  SetUpXmap();
78  SetUpCoeffs(m_xmap->GetNcoeffs());
79  }
80 
81  m_verts[0] = vertex[0];
82  m_verts[1] = vertex[1];
83  }
84 
86  int id,
87  const int coordim,
88  const PointGeomSharedPtr vertex[],
89  const CurveSharedPtr& curve):
90  Geometry1D(coordim)
91  {
93  m_eid = id;
94  m_globalID = id;
96  m_curve = curve;
97 
98  SetUpXmap();
99  SetUpCoeffs(m_xmap->GetNcoeffs());
100 
101  m_verts[0] = vertex[0];
102  m_verts[1] = vertex[1];
103  }
104 
105 
107  const int id,
108  const PointGeomSharedPtr& vert1,
109  const PointGeomSharedPtr& vert2):
110  Geometry1D(vert1->GetCoordim())
111  {
113 
114  m_verts[0] = vert1;
115  m_verts[1] = vert2;
116 
118 
119  m_eid = id;
120  m_globalID = id;
121 
122  SetUpXmap();
123  SetUpCoeffs(m_xmap->GetNcoeffs());
124  }
125 
127  {
128  // From Geometry class
130 
131  // info from EdgeComponent class
132  m_eid = in.m_eid;
133  m_globalID = in.m_globalID;
134 
135  std::list<CompToElmt>::const_iterator def;
136  for(def = in.m_elmtMap.begin(); def != in.m_elmtMap.end(); def++)
137  {
138  m_elmtMap.push_back(*def);
139  }
140  m_xmap = in.m_xmap;
141  SetUpCoeffs(m_xmap->GetNcoeffs());
142 
143  // info from SegGeom class
144  m_coordim = in.m_coordim;
145  m_verts[0] = in.m_verts[0];
146  m_verts[1] = in.m_verts[1];
147 
148  m_state = in.m_state;
149  }
150 
152  {
153  if (m_curve)
154  {
155  int npts = m_curve->m_points.size();
159  }
160  else
161  {
162  const LibUtilities::BasisKey B(
166  }
167  }
168 
169 
170  /** \brief Generate a one dimensional space segment geometry where
171  the vert[0] has the same x value and vert[1] is set to
172  vert[0] plus the length of the original segment
173 
174  **/
176  {
178 
179  // info about numbering
180  returnval->m_eid = m_eid;
181  returnval->m_globalID = m_globalID;
182  returnval->m_elmtMap = m_elmtMap;
183 
184 
185  // geometric information.
186  returnval->m_coordim = 1;
187  NekDouble x0 = (*m_verts[0])[0];
189  vert0->SetGlobalID(vert0->GetVid());
190  returnval->m_verts[0] = vert0;
191 
192  // Get information to calculate length.
195  v.push_back(base[0]->GetPointsKey());
197 
198  const Array<OneD, const NekDouble> jac = m_geomFactors->GetJac(v);
199 
200  NekDouble len = 0.0;
201  if(jac.num_elements() == 1)
202  {
203  len = jac[0]*2.0;
204  }
205  else
206  {
207  Array<OneD, const NekDouble> w0 = base[0]->GetW();
208  len = 0.0;
209 
210  for(int i = 0; i < jac.num_elements(); ++i)
211  {
212  len += jac[i]*w0[i];
213  }
214  }
215  // Set up second vertex.
217  vert1->SetGlobalID(vert1->GetVid());
218 
219  returnval->m_verts[1] = vert1;
220 
221  // at present just use previous m_xmap[0];
222  returnval->m_xmap = m_xmap;
223  returnval->SetUpCoeffs(m_xmap->GetNcoeffs());
224  returnval->m_state = eNotFilled;
225 
226  return returnval;
227  }
228 
229 
231  {
232  }
233 
234  /** given local collapsed coordinate Lcoord return the value of
235  physical coordinate in direction i **/
237  const int i,
238  const Array<OneD, const NekDouble>& Lcoord)
239  {
241  "Geometry is not in physical space");
242 
243  Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
244  m_xmap->BwdTrans(m_coeffs[i], tmp);
245 
246  return m_xmap->PhysEvaluate(Lcoord, tmp);
247  }
248 
249  void SegGeom::v_AddElmtConnected(int gvo_id, int locid)
250  {
251  CompToElmt ee(gvo_id,locid);
252  m_elmtMap.push_back(ee);
253  }
254 
256  {
257  return int(m_elmtMap.size());
258  }
259 
260  bool SegGeom::v_IsElmtConnected(int gvo_id, int locid) const
261  {
262  std::list<CompToElmt>::const_iterator def;
263  CompToElmt ee(gvo_id,locid);
264 
265  def = find(m_elmtMap.begin(),m_elmtMap.end(),ee);
266 
267  // Found the element connectivity object in the list
268  if(def != m_elmtMap.end())
269  {
270  return(true);
271  }
272 
273  return(false);
274  }
275 
276 
277 
278  /// \brief Get the orientation of edge1.
279  ///
280  /// If edge1 is connected to edge2 in the same direction as
281  /// the points comprising edge1 then it is forward, otherwise
282  /// it is backward.
283  ///
284  /// For example, assume edge1 is comprised of points 1 and 2,
285  /// and edge2 is comprised of points 2 and 3, then edge1 is
286  /// forward.
287  ///
288  /// If edge1 is comprised of points 2 and 1 and edge2 is
289  /// comprised of points 3 and 2, then edge1 is backward.
290  ///
291  /// Since both edges are passed, it does
292  /// not need any information from the EdgeComponent instance.
294  const SegGeom& edge1,
295  const SegGeom& edge2)
296  {
298 
299  /// Backward direction. Vertex 0 is connected to edge 2.
300  if ((*edge1.GetVertex(0) == *edge2.GetVertex(0)) ||
301  (*edge1.GetVertex(0) == *edge2.GetVertex(1)))
302  {
303  returnval = StdRegions::eBackwards;
304  }
305  // Not forward either, then we have a problem.
306  else if ((*edge1.GetVertex(1) != *edge2.GetVertex(0)) &&
307  (*edge1.GetVertex(1) != *edge2.GetVertex(1)))
308  {
309  std::ostringstream errstrm;
310  errstrm << "Connected edges do not share a vertex. Edges ";
311  errstrm << edge1.GetEid() << ", " << edge2.GetEid();
312  ASSERTL0(false, errstrm.str());
313  }
314 
315  return returnval;
316  }
317 
318  /// Set up GeoFac for this geometry using Coord quadrature distribution
320  {
322  {
324 
326 
327  if(m_xmap->GetBasisNumModes(0)!=2)
328  {
329  gType = eDeformed;
330  }
331 
333  gType, m_coordim, m_xmap, m_coeffs);
335  }
336  }
337 
338 
339  /// \brief put all quadrature information into edge structure and
340  /// backward transform
342  {
343  if(m_state != ePtsFilled)
344  {
345  int i;
346 
347  if (m_coordim > 0 && m_curve)
348  {
349  int npts = m_curve->m_points.size();
351  Array<OneD,NekDouble> tmp(npts);
352 
353  if(m_verts[0]->dist(*(m_curve->m_points[0])) > NekConstants::kVertexTheSameDouble)
354  {
355  std::string err = "Vertex 0 is separated from first point by more than ";
356  std::stringstream strstrm;
358  << " in edge " << m_globalID;
359  err += strstrm.str();
360  NEKERROR(ErrorUtil::ewarning, err.c_str());
361  }
362 
363  if(m_verts[1]->dist(*(m_curve->m_points[npts-1])) > NekConstants::kVertexTheSameDouble)
364  {
365  std::string err = "Vertex 1 is separated from last point by more than ";
366  std::stringstream strstrm;
368  << " in edge " << m_globalID;
369  err += strstrm.str();
370  NEKERROR(ErrorUtil::ewarning, err.c_str());
371  }
372 
373  LibUtilities::PointsKey fkey(npts,m_curve->m_ptype);
374  DNekMatSharedPtr I0 =
375  LibUtilities::PointsManager()[fkey]->GetI(pkey);
376  NekVector<NekDouble> out(npts+1);
377 
378  for(int i = 0; i < m_coordim; ++i)
379  {
380  // Load up coordinate values into tmp
381  for(int j = 0; j < npts; ++j)
382  {
383  tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
384  }
385 
386  // Interpolate to GLL points
387  NekVector<NekDouble> in (npts, tmp, eWrapper);
388  out = (*I0)*in;
389 
390  m_xmap->FwdTrans(out.GetPtr(), m_coeffs[i]);
391  }
392  }
393 
394  for (i = 0; i < m_coordim; ++i)
395  {
396  m_coeffs[i][0] = (*m_verts[0])[i];
397  m_coeffs[i][1] = (*m_verts[1])[i];
398  }
399 
401  }
402  }
403 
404  void SegGeom::v_Reset(CurveMap &curvedEdges,
405  CurveMap &curvedFaces)
406  {
407  Geometry::v_Reset(curvedEdges, curvedFaces);
408  CurveMap::iterator it = curvedEdges.find(m_globalID);
409 
410  if (it != curvedEdges.end())
411  {
412  m_curve = it->second;
413  }
414 
415  SetUpXmap();
416  SetUpCoeffs(m_xmap->GetNcoeffs());
417  }
418 
420  const Array<OneD, const NekDouble>& coords,
421  Array<OneD,NekDouble>& Lcoords)
422  {
423  int i;
424  NekDouble resid = 0.0;
426 
427  // calculate local coordinate for coord
428  if(GetMetricInfo()->GetGtype() == eRegular)
429  {
430  NekDouble len = 0.0;
431  NekDouble xi = 0.0;
432 
433  const int npts = m_xmap->GetTotPoints();
434  Array<OneD, NekDouble> pts(npts);
435 
436  for(i = 0; i < m_coordim; ++i)
437  {
438  m_xmap->BwdTrans(m_coeffs[i], pts);
439  len += (pts[npts-1]-pts[0])*(pts[npts-1]-pts[0]);
440  xi += (coords[i]-pts[0])*(coords[i]-pts[0]);
441  }
442 
443  len = sqrt(len);
444  xi = sqrt(xi);
445 
446  Lcoords[0] = 2*xi/len-1.0;
447  }
448  else
449  {
451  "inverse mapping must be set up to use this call");
452  }
453  return resid;
454  }
455 
456  /**
457  * @brief Determines if a point specified in global coordinates is
458  * located within this tetrahedral geometry.
459  */
461  const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
462  {
463  Array<OneD,NekDouble> locCoord(GetCoordim(),0.0);
464  return v_ContainsPoint(gloCoord,locCoord,tol);
465 
466  }
467 
469  const Array<OneD, const NekDouble>& gloCoord,
470  Array<OneD, NekDouble> &stdCoord,
471  NekDouble tol)
472  {
473  NekDouble resid;
474  return v_ContainsPoint(gloCoord,stdCoord,tol,resid);
475  }
476 
478  const Array<OneD, const NekDouble>& gloCoord,
479  Array<OneD, NekDouble> &stdCoord,
480  NekDouble tol,
481  NekDouble &resid)
482  {
483  resid = GetLocCoords(gloCoord, stdCoord);
484  if (stdCoord[0] >= -(1+tol) && stdCoord[0] <= 1+tol)
485  {
486  return true;
487  }
488  return false;
489  }
490 
491  int SegGeom::v_GetVid(int i) const
492  {
493  ASSERTL2((i >=0) && (i <= 1),"Verted id must be between 0 and 1");
494  return m_verts[i]->GetVid();
495  }
496 
498  {
499  PointGeomSharedPtr returnval;
500 
501  if (i >= 0 && i < kNverts)
502  {
503  returnval = m_verts[i];
504  }
505 
506  return returnval;
507  }
508 
509  int SegGeom::v_GetEid() const
510  {
511  return m_eid;
512  }
513 
515  {
516  return m_xmap->GetBasis(i);
517  }
518 
520  {
521  return m_xmap;
522  }
523 
525  {
526  m_ownData = true;
527  }
528 
530  {
531  //cout << "StdRegions::PointOrientation GetPorient"<<endl;
532  ASSERTL2((i >=0) && (i <= 1),"Point id must be between 0 and 1");
533 
534  if (i % 2 == 0)
535  {
536  return StdRegions::eBwd;
537  }
538  else
539  {
540  return StdRegions::eFwd;
541  }
542  }
543 
544 
546  {
547  return LibUtilities::eSegment;
548  }
549 
551  {
552  return kNverts;
553  }
554 
556  {
557  return kNedges;
558  }
559 
560 
561  }; //end of namespace
562 }; //end of namespace
563 
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
virtual PointGeomSharedPtr v_GetVertex(const int i) const
Definition: SegGeom.cpp:497
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
GeomFactorsSharedPtr m_geomFactors
Definition: Geometry.h:170
Structure holding graphvertexobject id and local element facet id.
virtual int v_GetNumVerts() const
Definition: SegGeom.cpp:550
bool m_ownData
Boolean indicating whether object owns the data.
Definition: SegGeom.h:173
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Definition: SegGeom.cpp:419
virtual StdRegions::StdExpansionSharedPtr v_GetXmap() const
Definition: SegGeom.cpp:519
Principle Modified Functions .
Definition: BasisType.h:49
virtual int v_GetNumEdges() const
Definition: SegGeom.cpp:555
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: Geometry.cpp:307
virtual bool v_IsElmtConnected(int gvoId, int locId) const
Definition: SegGeom.cpp:260
static const NekDouble kVertexTheSameDouble
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
virtual StdRegions::Orientation v_GetPorient(const int i) const
Definition: SegGeom.cpp:529
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:293
boost::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:62
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:175
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: SegGeom.cpp:545
virtual void v_GenGeomFactors()
Set up GeoFac for this geometry using Coord quadrature distribution.
Definition: SegGeom.cpp:319
Geometric information has not been generated.
static const int kNedges
Definition: SegGeom.h:99
virtual int v_NumElmtConnected() const
Definition: SegGeom.cpp:255
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
virtual const LibUtilities::BasisSharedPtr v_GetBasis(const int i)
Definition: SegGeom.cpp:514
static std::string npts
Definition: InputFld.cpp:43
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Definition: Points.h:58
PointGeomSharedPtr GetVertex(const int i) const
Definition: Geometry1D.cpp:56
double NekDouble
SpatialDomains::PointGeomSharedPtr m_verts[kNverts]
Definition: SegGeom.h:104
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
1D geometry information
Definition: Geometry1D.h:56
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
Geometry is straight-sided with constant geometric factors.
int GetVid(int i) const
Definition: Geometry.h:319
std::list< CompToElmt > m_elmtMap
Definition: SegGeom.h:103
NekDouble GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Definition: Geometry.h:450
virtual NekDouble v_GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord)
Definition: SegGeom.cpp:236
LibUtilities::ShapeType m_shapeType
Definition: Geometry.h:177
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:312
virtual void v_FillGeom()
put all quadrature information into edge structure and backward transform
Definition: SegGeom.cpp:341
virtual bool v_ContainsPoint(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
Determines if a point specified in global coordinates is located within this tetrahedral geometry...
Definition: SegGeom.cpp:460
boost::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:63
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: SegGeom.cpp:404
Geometric information has been generated.
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
Definition: Geometry.h:484
GeomFactorsSharedPtr GetMetricInfo()
Definition: Geometry.h:299
virtual void v_AddElmtConnected(int gvo_id, int locid)
Definition: SegGeom.cpp:249
GeomType
Indicates the type of element geometry.
virtual int v_GetEid() const
Definition: SegGeom.cpp:509
boost::shared_ptr< Basis > BasisSharedPtr
GeomState m_state
enum identifier to determine if quad points are filled
Definition: Geometry.h:175
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Geometry is curved or has non-constant factors.
int m_coordim
coordinate dimension
Definition: Geometry.h:169
Describes the specification for a Basis.
Definition: Basis.h:50
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
static const int kNverts
Definition: SegGeom.h:98
virtual int v_GetVid(int i) const
Definition: SegGeom.cpp:491