Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  const LibUtilities::BasisKey B(
59 
61  m_eid = id;
62  m_globalID = id;
64  SetUpCoeffs(m_xmap->GetNcoeffs());
65  }
66 
68  int id,
69  const int coordim,
70  const PointGeomSharedPtr vertex[]):
71  Geometry1D(coordim)
72  {
74  m_eid = id;
75  m_globalID = id;
77 
78  if (coordim > 0)
79  {
83  )
84  );
85 
87  SetUpCoeffs(m_xmap->GetNcoeffs());
88  }
89 
90  m_verts[0] = vertex[0];
91  m_verts[1] = vertex[1];
92  }
93 
95  int id,
96  const int coordim,
97  const PointGeomSharedPtr vertex[],
98  const CurveSharedPtr& curve):
99  Geometry1D(coordim)
100  {
102  m_eid = id;
103  m_globalID = id;
105 
106  if (coordim > 0)
107  {
108  int npts = curve->m_points.size();
111 
112  Array<OneD,NekDouble> tmp(npts);
113 
114  if(vertex[0]->dist(*(curve->m_points[0])) > NekConstants::kVertexTheSameDouble)
115  {
116  cout<<"edge="<<id<<endl;
117  std::string err = "Vertex 0 is separated from first point by more than ";
118  std::stringstream strstrm;
120  err += strstrm.str();
121  NEKERROR(ErrorUtil::ewarning, err.c_str());
122  }
123 
124  if(vertex[1]->dist(*(curve->m_points[npts-1])) > NekConstants::kVertexTheSameDouble)
125  {
126  cout<<"edge="<<id<<endl;
127  std::string err = "Vertex 1 is separated from last point by more than ";
128  std::stringstream strstrm;
130  err += strstrm.str();
131  NEKERROR(ErrorUtil::ewarning, err.c_str());
132  }
133 
135  SetUpCoeffs(m_xmap->GetNcoeffs());
136 
137  for(int i = 0; i < m_coordim; ++i)
138  {
139  // Load up coordinate values into tmp
140  for(int j = 0; j < npts; ++j)
141  {
142  tmp[j] = (curve->m_points[j]->GetPtr())[i];
143  }
144 
145  // Interpolate to GLL points
146  DNekMatSharedPtr I0;
147 
148  LibUtilities::PointsKey fkey(npts,curve->m_ptype);
149  I0 = LibUtilities::PointsManager()[fkey]->GetI(pkey);
150 
151  NekVector<NekDouble> in (npts, tmp, eWrapper);
152  NekVector<NekDouble> out(npts+1);
153  out = (*I0)*in;
154 
155  m_xmap->FwdTrans(out.GetPtr(), m_coeffs[i]);
156  }
157  }
158 
159  m_verts[0] = vertex[0];
160  m_verts[1] = vertex[1];
161  }
162 
163 
165  const int id,
166  const PointGeomSharedPtr& vert1,
167  const PointGeomSharedPtr& vert2):
168  Geometry1D(vert1->GetCoordim())
169  {
171 
172  m_verts[0] = vert1;
173  m_verts[1] = vert2;
174 
176 
180  )
181  );
182  m_eid = id;
183  m_globalID = id;
185  SetUpCoeffs(m_xmap->GetNcoeffs());
186  }
187 
189  {
190  // From Geometry class
192 
193  // info from EdgeComponent class
194  m_eid = in.m_eid;
195  m_globalID = in.m_globalID;
196 
197  std::list<CompToElmt>::const_iterator def;
198  for(def = in.m_elmtMap.begin(); def != in.m_elmtMap.end(); def++)
199  {
200  m_elmtMap.push_back(*def);
201  }
202  m_xmap = in.m_xmap;
203  SetUpCoeffs(m_xmap->GetNcoeffs());
204 
205  // info from SegGeom class
206  m_coordim = in.m_coordim;
207  m_verts[0] = in.m_verts[0];
208  m_verts[1] = in.m_verts[1];
209 
210  m_state = in.m_state;
211  }
212 
213 
214  /** \brief Generate a one dimensional space segment geometry where
215  the vert[0] has the same x value and vert[1] is set to
216  vert[0] plus the length of the original segment
217 
218  **/
220  {
222 
223  // info about numbering
224  returnval->m_eid = m_eid;
225  returnval->m_globalID = m_globalID;
226  returnval->m_elmtMap = m_elmtMap;
227 
228 
229  // geometric information.
230  returnval->m_coordim = 1;
231  NekDouble x0 = (*m_verts[0])[0];
233  vert0->SetGlobalID(vert0->GetVid());
234  returnval->m_verts[0] = vert0;
235 
236  // Get information to calculate length.
237  const Array<OneD, const LibUtilities::BasisSharedPtr> base = m_xmap->GetBase();
239  v.push_back(base[0]->GetPointsKey());
241 
242  const Array<OneD, const NekDouble> jac = m_geomFactors->GetJac(v);
243 
244  NekDouble len;
245  if(jac.num_elements() == 1)
246  {
247  len = jac[0]*2.0;
248  }
249  else
250  {
251  Array<OneD, const NekDouble> w0 = base[0]->GetW();
252  len = 0.0;
253 
254  for(int i = 0; i < jac.num_elements(); ++i)
255  {
256  len += jac[i]*w0[i];
257  }
258  }
259  // Set up second vertex.
261  vert1->SetGlobalID(vert1->GetVid());
262 
263  returnval->m_verts[1] = vert1;
264 
265  // at present just use previous m_xmap[0];
266  returnval->m_xmap = m_xmap;
267  returnval->SetUpCoeffs(m_xmap->GetNcoeffs());
268  returnval->m_state = eNotFilled;
269 
270  return returnval;
271  }
272 
273 
275  {
276  }
277 
278  /** given local collapsed coordinate Lcoord return the value of
279  physical coordinate in direction i **/
281  const int i,
282  const Array<OneD, const NekDouble>& Lcoord)
283  {
285  "Geometry is not in physical space");
286 
287  Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
288  m_xmap->BwdTrans(m_coeffs[i], tmp);
289 
290  return m_xmap->PhysEvaluate(Lcoord, tmp);
291  }
292 
293  void SegGeom::v_AddElmtConnected(int gvo_id, int locid)
294  {
295  CompToElmt ee(gvo_id,locid);
296  m_elmtMap.push_back(ee);
297  }
298 
300  {
301  return int(m_elmtMap.size());
302  }
303 
304  bool SegGeom::v_IsElmtConnected(int gvo_id, int locid) const
305  {
306  std::list<CompToElmt>::const_iterator def;
307  CompToElmt ee(gvo_id,locid);
308 
309  def = find(m_elmtMap.begin(),m_elmtMap.end(),ee);
310 
311  // Found the element connectivity object in the list
312  if(def != m_elmtMap.end())
313  {
314  return(true);
315  }
316 
317  return(false);
318  }
319 
320 
321 
322  /// \brief Get the orientation of edge1.
323  ///
324  /// If edge1 is connected to edge2 in the same direction as
325  /// the points comprising edge1 then it is forward, otherwise
326  /// it is backward.
327  ///
328  /// For example, assume edge1 is comprised of points 1 and 2,
329  /// and edge2 is comprised of points 2 and 3, then edge1 is
330  /// forward.
331  ///
332  /// If edge1 is comprised of points 2 and 1 and edge2 is
333  /// comprised of points 3 and 2, then edge1 is backward.
334  ///
335  /// Since both edges are passed, it does
336  /// not need any information from the EdgeComponent instance.
338  const SegGeom& edge1,
339  const SegGeom& edge2)
340  {
342 
343  /// Backward direction. Vertex 0 is connected to edge 2.
344  if ((*edge1.GetVertex(0) == *edge2.GetVertex(0)) ||
345  (*edge1.GetVertex(0) == *edge2.GetVertex(1)))
346  {
347  returnval = StdRegions::eBackwards;
348  }
349  // Not forward either, then we have a problem.
350  else if ((*edge1.GetVertex(1) != *edge2.GetVertex(0)) &&
351  (*edge1.GetVertex(1) != *edge2.GetVertex(1)))
352  {
353  std::ostringstream errstrm;
354  errstrm << "Connected edges do not share a vertex. Edges ";
355  errstrm << edge1.GetEid() << ", " << edge2.GetEid();
356  ASSERTL0(false, errstrm.str());
357  }
358 
359  return returnval;
360  }
361 
362  /// Set up GeoFac for this geometry using Coord quadrature distribution
364  {
366  {
368 
370 
371  if(m_xmap->GetBasisNumModes(0)!=2)
372  {
373  gType = eDeformed;
374  }
375 
377  gType, m_coordim, m_xmap, m_coeffs);
379  }
380  }
381 
382 
383  /// \brief put all quadrature information into edge structure and
384  /// backward transform
386  {
387  if(m_state != ePtsFilled)
388  {
389  int i;
390 
391  for (i = 0; i < m_coordim; ++i)
392  {
393  m_coeffs[i][0] = (*m_verts[0])[i];
394  m_coeffs[i][1] = (*m_verts[1])[i];
395  }
396 
398  }
399  }
400 
402  const Array<OneD, const NekDouble>& coords,
403  Array<OneD,NekDouble>& Lcoords)
404  {
405  int i;
406  NekDouble resid = 0.0;
408 
409  // calculate local coordinate for coord
410  if(GetMetricInfo()->GetGtype() == eRegular)
411  {
412  NekDouble len = 0.0;
413  NekDouble xi = 0.0;
414 
415  const int npts = m_xmap->GetTotPoints();
416  Array<OneD, NekDouble> pts(npts);
417 
418  for(i = 0; i < m_coordim; ++i)
419  {
420  m_xmap->BwdTrans(m_coeffs[i], pts);
421  len += (pts[npts-1]-pts[0])*(pts[npts-1]-pts[0]);
422  xi += (coords[i]-pts[0])*(coords[i]-pts[0]);
423  }
424 
425  len = sqrt(len);
426  xi = sqrt(xi);
427 
428  Lcoords[0] = 2*xi/len-1.0;
429  }
430  else
431  {
433  "inverse mapping must be set up to use this call");
434  }
435  return resid;
436  }
437 
438  /**
439  * @brief Determines if a point specified in global coordinates is
440  * located within this tetrahedral geometry.
441  */
443  const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
444  {
445  Array<OneD,NekDouble> locCoord(GetCoordim(),0.0);
446  return v_ContainsPoint(gloCoord,locCoord,tol);
447 
448  }
449 
451  const Array<OneD, const NekDouble>& gloCoord,
452  Array<OneD, NekDouble> &stdCoord,
453  NekDouble tol)
454  {
455  NekDouble resid;
456  return v_ContainsPoint(gloCoord,stdCoord,tol,resid);
457  }
458 
460  const Array<OneD, const NekDouble>& gloCoord,
461  Array<OneD, NekDouble> &stdCoord,
462  NekDouble tol,
463  NekDouble &resid)
464  {
465  resid = GetLocCoords(gloCoord, stdCoord);
466  if (stdCoord[0] >= -(1+tol) && stdCoord[0] <= 1+tol)
467  {
468  return true;
469  }
470  return false;
471  }
472 
473  int SegGeom::v_GetVid(int i) const
474  {
475  ASSERTL2((i >=0) && (i <= 1),"Verted id must be between 0 and 1");
476  return m_verts[i]->GetVid();
477  }
478 
480  {
481  PointGeomSharedPtr returnval;
482 
483  if (i >= 0 && i < kNverts)
484  {
485  returnval = m_verts[i];
486  }
487 
488  return returnval;
489  }
490 
491  int SegGeom::v_GetEid() const
492  {
493  return m_eid;
494  }
495 
497  {
498  return m_xmap->GetBasis(i);
499  }
500 
502  {
503  return m_xmap;
504  }
505 
507  {
508  m_ownData = true;
509  }
510 
512  {
513  //cout << "StdRegions::PointOrientation GetPorient"<<endl;
514  ASSERTL2((i >=0) && (i <= 1),"Point id must be between 0 and 1");
515 
516  if (i % 2 == 0)
517  {
518  return StdRegions::eBwd;
519  }
520  else
521  {
522  return StdRegions::eFwd;
523  }
524  }
525 
526 
528  {
529  return LibUtilities::eSegment;
530  }
531 
533  {
534  return kNverts;
535  }
536 
538  {
539  return kNedges;
540  }
541 
542 
543  }; //end of namespace
544 }; //end of namespace
545