Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Face.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Face.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: Mesh face object.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 
39 
41 
42 namespace Nektar
43 {
44 namespace NekMeshUtils
45 {
46 
48  std::vector<NodeSharedPtr> &nodeList)
49 {
50  // Treat 2D point distributions differently to 3D.
54  {
55  int n = m_edgeList[0]->GetNodeCount();
56  int n2 = m_edgeList[1]->GetNodeCount();
57  int n3 = m_edgeList[2]->GetNodeCount();
58 
59  bool same = (n == n2 ? (n2 == n3) : false);
60  ASSERTL0(same, "Edges are not consistent");
61 
62  nodeList.insert(
63  nodeList.end(), m_vertexList.begin(), m_vertexList.end());
64  for (int k = 0; k < m_edgeList.size(); ++k)
65  {
66  nodeList.insert(nodeList.end(),
67  m_edgeList[k]->m_edgeNodes.begin(),
68  m_edgeList[k]->m_edgeNodes.end());
69  if (m_edgeList[k]->m_n1 != m_vertexList[k])
70  {
71  // If edge orientation is reversed relative to node
72  // ordering, we need to reverse order of nodes.
73  std::reverse(nodeList.begin() + 3 + k * (n - 2),
74  nodeList.begin() + 3 + (k + 1) * (n - 2));
75  }
76  }
77  nodeList.insert(
78  nodeList.end(), m_faceNodes.begin(), m_faceNodes.end());
79  }
80  else
81  {
82  // Write out in 2D tensor product order.
83  ASSERTL0(m_vertexList.size() == 4,
84  "Face nodes of tensor product only supported "
85  "for quadrilaterals.");
86 
87  int n = (int)sqrt((NekDouble)GetNodeCount());
88  nodeList.resize(n * n);
89 
90  ASSERTL0(n * n == GetNodeCount(), "Wrong number of modes?");
91 
92  // Write vertices
93  nodeList[0] = m_vertexList[0];
94  nodeList[n - 1] = m_vertexList[1];
95  nodeList[n * n - 1] = m_vertexList[2];
96  nodeList[n * (n - 1)] = m_vertexList[3];
97 
98  // Write edge-interior
99  int skips[4][2] = {
100  {0, 1}, {n - 1, n}, {n * n - 1, -1}, {n * (n - 1), -n}};
101  for (int i = 0; i < 4; ++i)
102  {
103  bool reverseEdge = m_edgeList[i]->m_n1 == m_vertexList[i];
104 
105  if (!reverseEdge)
106  {
107  for (int j = 1; j < n - 1; ++j)
108  {
109  nodeList[skips[i][0] + j * skips[i][1]] =
110  m_edgeList[i]->m_edgeNodes[n - 2 - j];
111  }
112  }
113  else
114  {
115  for (int j = 1; j < n - 1; ++j)
116  {
117  nodeList[skips[i][0] + j * skips[i][1]] =
118  m_edgeList[i]->m_edgeNodes[j - 1];
119  }
120  }
121  }
122 
123  // Write interior
124  for (int i = 1; i < n - 1; ++i)
125  {
126  for (int j = 1; j < n - 1; ++j)
127  {
128  nodeList[i * n + j] =
129  m_faceNodes[(i - 1) * (n - 2) + (j - 1)];
130  }
131  }
132  }
133 }
134 
135 
137 {
138  std::stringstream s;
139  std::string str;
140  std::vector<NodeSharedPtr> nodeList;
141 
142  // assemble listof nodes
143  GetCurvedNodes(nodeList);
144 
145  // put them into a string
146  for (int k = 0; k < nodeList.size(); ++k)
147  {
148  s << std::scientific << std::setprecision(8) << " "
149  << nodeList[k]->m_x << " " << nodeList[k]->m_y << " "
150  << nodeList[k]->m_z << " ";
151  }
152 
153  return s.str();
154 }
155 
156 void Face::MakeOrder(int order,
159  int coordDim,
160  int &id)
161 {
162  if (m_vertexList.size() == 3)
163  {
164  // Triangles of order < 3 have no interior volume points.
165  if (order < 3)
166  {
167  m_faceNodes.clear();
168  return;
169  }
170 
171  int nPoints = order + 1;
172  StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
173 
174  Array<OneD, NekDouble> px, py;
175  LibUtilities::PointsKey pKey(nPoints, pType);
176  ASSERTL1(pKey.GetPointsDim() == 2, "Points distribution must be 2D");
177  LibUtilities::PointsManager()[pKey]->GetPoints(px, py);
178 
179  Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
180 
181  for (int i = 0; i < coordDim; ++i)
182  {
183  phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
184  xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
185  }
186 
187  const int nTriPts = nPoints * (nPoints + 1) / 2;
188  const int nTriIntPts = (nPoints - 3) * (nPoints - 2) / 2;
189  m_faceNodes.resize(nTriIntPts);
190 
191  for (int i = 3 + 3*(nPoints-2), cnt = 0; i < nTriPts; ++i, ++cnt)
192  {
194  xp[0] = px[i];
195  xp[1] = py[i];
196 
197  Array<OneD, NekDouble> x(3, 0.0);
198  for (int j = 0; j < coordDim; ++j)
199  {
200  x[j] = xmap->PhysEvaluate(xp, phys[j]);
201  }
202 
203  m_faceNodes[cnt] = boost::shared_ptr<Node>(
204  new Node(id++, x[0], x[1], x[2]));
205  }
206  m_curveType = pType;
207  }
208  else if (m_vertexList.size() == 4)
209  {
210  // Quads of order < 2 have no interior volume points.
211  if (order < 2)
212  {
213  m_faceNodes.clear();
214  return;
215  }
216 
217  int nPoints = order + 1;
218  StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
219 
221  LibUtilities::PointsKey pKey(nPoints, pType);
222  ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D");
223  LibUtilities::PointsManager()[pKey]->GetPoints(px);
224 
225  Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
226 
227  for (int i = 0; i < coordDim; ++i)
228  {
229  phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
230  xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
231  }
232 
233  int nQuadIntPts = (nPoints - 2) * (nPoints - 2);
234  m_faceNodes.resize(nQuadIntPts);
235 
236  for (int i = 1, cnt = 0; i < nPoints-1; ++i)
237  {
238  for (int j = 1; j < nPoints-1; ++j, ++cnt)
239  {
241  xp[0] = px[j];
242  xp[1] = px[i];
243 
244  Array<OneD, NekDouble> x(3, 0.0);
245  for (int k = 0; k < coordDim; ++k)
246  {
247  x[k] = xmap->PhysEvaluate(xp, phys[k]);
248  }
249 
250  m_faceNodes[cnt] = boost::shared_ptr<Node>(
251  new Node(id++, x[0], x[1], x[2]));
252  }
253  }
254 
255  m_curveType = pType;
256  }
257  else
258  {
259  ASSERTL0(false, "Unknown number of vertices");
260  }
261 
262  if(m_parentCAD)
263  {
264  CADSurfSharedPtr s = boost::dynamic_pointer_cast<CADSurf>(m_parentCAD);
265  for(int i = 0; i < m_faceNodes.size(); i++)
266  {
267  Array<OneD, NekDouble> loc(3);
268  loc[0] = m_faceNodes[i]->m_x;
269  loc[1] = m_faceNodes[i]->m_y;
270  loc[2] = m_faceNodes[i]->m_z;
272  s->ProjectTo(loc,uv);
273  loc = s->P(uv);
274  m_faceNodes[i]->m_x = loc[0];
275  m_faceNodes[i]->m_y = loc[1];
276  m_faceNodes[i]->m_z = loc[2];
277  m_faceNodes[i]->SetCADSurf(s->GetId(),s,uv);
278  }
279  }
280 }
281 
283 {
284  int nEdge = m_edgeList.size();
285 
288  StdRegions::Orientation edgeo[4];
289 
290  for (int i = 0; i < nEdge; ++i)
291  {
292  edges[i] = m_edgeList[i]->GetGeom(coordDim);
293  }
294 
295  for (int i = 0; i < nEdge; ++i)
296  {
297  edgeo[i] = m_edgeList[i]->m_n1 == m_vertexList[i]
300  }
301 
302  if (m_faceNodes.size() > 0)
303  {
304  if (nEdge == 3)
305  {
308  m_id, m_curveType);
309 
310  for (int j = 0; j < m_vertexList.size(); j++)
311  {
312  c->m_points.push_back(m_vertexList[j]->GetGeom(coordDim));
313  }
314  for (int j = 0; j < nEdge; j++)
315  {
316  std::vector<NodeSharedPtr> ed = m_edgeList[j]->m_edgeNodes;
317  if (edgeo[j] == StdRegions::eBackwards)
318  {
319  for (int k = ed.size() - 1; k >= 0; k--)
320  {
321  c->m_points.push_back(ed[k]->GetGeom(coordDim));
322  }
323  }
324  else
325  {
326  for (int k = 0; k < ed.size(); k++)
327  {
328  c->m_points.push_back(ed[k]->GetGeom(coordDim));
329  }
330  }
331  }
332  for (int j = 0; j < m_faceNodes.size(); j++)
333  {
334  c->m_points.push_back(m_faceNodes[j]->GetGeom(coordDim));
335  }
336 
338  m_id, edges, edgeo, c);
339  }
340  else
341  {
344  m_id, m_curveType);
345 
346  ASSERTL0(m_vertexList.size() == 4,
347  "Face nodes of tensor product only supported "
348  "for quadrilaterals.");
349 
350  int n = (int)sqrt((NekDouble)GetNodeCount());
351  std::vector<NodeSharedPtr> tmp(n * n);
352 
353  ASSERTL0(n * n == GetNodeCount(), "Wrong number of modes?");
354 
355  // Write vertices
356  tmp[0] = m_vertexList[0];
357  tmp[n - 1] = m_vertexList[1];
358  tmp[n * n - 1] = m_vertexList[2];
359  tmp[n * (n - 1)] = m_vertexList[3];
360 
361  // Write edge-interior
362  int skips[4][2] = {
363  {0, 1}, {n - 1, n}, {n * n - 1, -1}, {n * (n - 1), -n}};
364  for (int i = 0; i < 4; ++i)
365  {
366  bool reverseEdge = edgeo[i] == StdRegions::eBackwards;
367 
368  if (reverseEdge)
369  {
370  for (int j = 1; j < n - 1; ++j)
371  {
372  tmp[skips[i][0] + j * skips[i][1]] =
373  m_edgeList[i]->m_edgeNodes[n - 2 - j];
374  }
375  }
376  else
377  {
378  for (int j = 1; j < n - 1; ++j)
379  {
380  tmp[skips[i][0] + j * skips[i][1]] =
381  m_edgeList[i]->m_edgeNodes[j - 1];
382  }
383  }
384  }
385 
386  // Write interior
387  for (int i = 1; i < n - 1; ++i)
388  {
389  for (int j = 1; j < n - 1; ++j)
390  {
391  tmp[i * n + j] =
392  m_faceNodes[(i - 1) * (n - 2) + (j - 1)];
393  }
394  }
395 
396  for (int k = 0; k < tmp.size(); ++k)
397  {
398  c->m_points.push_back(tmp[k]->GetGeom(coordDim));
399  }
400 
401  ret =
403  m_id, edges, edgeo, c);
404  }
405  }
406  else
407  {
408  if (nEdge == 3)
409  {
411  m_id, edges, edgeo);
412  }
413  else
414  {
415  ret =
417  m_id, edges, edgeo);
418  }
419  }
420 
421  return ret;
422 }
423 
424 }
425 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
NEKMESHUTILS_EXPORT void GetCurvedNodes(std::vector< NodeSharedPtr > &nodeList)
Assemble a list of nodes on curved face.
Definition: Face.cpp:47
void MakeOrder(int order, SpatialDomains::GeometrySharedPtr geom, LibUtilities::PointsType pType, int coordDim, int &id)
Make this face an order order face.
Definition: Face.cpp:156
std::vector< NodeSharedPtr > m_faceNodes
List of face-interior nodes defining the shape of the face.
Definition: Face.h:137
NEKMESHUTILS_EXPORT unsigned int GetNodeCount() const
Returns the total number of nodes (vertices, edge nodes and face nodes).
Definition: Face.h:100
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
std::vector< EdgeSharedPtr > m_edgeList
List of corresponding edges.
Definition: Face.h:135
unsigned int m_id
ID of the face.
Definition: Face.h:131
boost::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:62
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
base class for a cad surface
Definition: CADSurf.h:55
unsigned int GetPointsDim() const
Definition: Points.h:149
PointsManagerT & PointsManager(void)
NEKMESHUTILS_EXPORT std::string GetXmlCurveString()
Generates a string listing the coordinates of all nodes associated with this face.
Definition: Face.cpp:136
Defines a specification for a set of points.
Definition: Points.h:58
LibUtilities::PointsType m_curveType
Distribution of points in this face.
Definition: Face.h:139
double NekDouble
2D Nodal Fekete Points on a Triangle
Definition: PointsType.h:71
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
NEKMESHUTILS_EXPORT SpatialDomains::Geometry2DSharedPtr GetGeom(int coordDim)
Generate either SpatialDomains::TriGeom or SpatialDomains::QuadGeom for this element.
Definition: Face.cpp:282
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
2D Evenly-spaced points on a Triangle
Definition: PointsType.h:72
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:228
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
CADObjectSharedPtr m_parentCAD
Definition: Face.h:145
std::vector< NodeSharedPtr > m_vertexList
List of vertex nodes.
Definition: Face.h:133
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:70