Nektar++
Deform.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Deform.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: Deformation of mesh from fields.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 
42 #include <StdRegions/StdSegExp.h>
43 #include <StdRegions/StdQuadExp.h>
45 #include <MultiRegions/ExpList.h>
46 
47 namespace Nektar {
48 namespace SolverUtils {
49 
50  /**
51  * @brief Update geometry according to displacement that is in current
52  * fields.
53  *
54  * @param graph The MeshGraph of the current geometry.
55  * @param fields The fields containing the displacement.
56  */
60  {
61  // Clear existing curvature.
62  SpatialDomains::CurveMap &curvedEdges = graph->GetCurvedEdges();
63  SpatialDomains::CurveMap &curvedFaces = graph->GetCurvedFaces();
64  curvedEdges.clear();
65  curvedFaces.clear();
66 
67  int i, j, k, l, dim;
68 
69  // Sets to hold IDs of updated vertices to avoid duplicating effort.
70  set<int> updatedVerts, updatedEdges, updatedFaces;
71 
72  dim = graph->GetSpaceDimension();
75 
76  for (i = 0; i < fields[0]->GetExpSize(); ++i)
77  {
78  LocalRegions::ExpansionSharedPtr exp = fields[0]->GetExp(i);
79  int offset = fields[0]->GetPhys_Offset(i);
80  int nquad = exp->GetTotPoints();
81 
82  // Extract displacement for this element, allocate storage for
83  // elemental coordinates.
84  for (j = 0; j < dim; ++j)
85  {
86  phys[j] = Array<OneD, NekDouble>(
87  nquad, fields[j]->UpdatePhys() + offset);
88  coord[j] = Array<OneD, NekDouble>(nquad);
89  }
90 
91  // In 2D loop over edges.
92  if (dim == 2)
93  {
94  exp->GetCoords(coord[0], coord[1]);
95 
97  boost::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
98  exp->GetGeom());
99 
100  for (j = 0; j < exp->GetNedges(); ++j)
101  {
103 
104  // This edge has already been processed.
105  if (updatedEdges.find(edge->GetGlobalID()) !=
106  updatedEdges.end())
107  {
108  continue;
109  }
110 
111  // Extract edge displacement.
112  int nEdgePts = exp->GetEdgeNumPoints(j);
113  Array<OneD, Array<OneD, NekDouble> > edgePhys (dim);
114  Array<OneD, Array<OneD, NekDouble> > edgeCoord(dim);
115 
116  const LibUtilities::BasisKey B(
117  LibUtilities::eModified_A, nEdgePts,
121  StdRegions::StdSegExp>::AllocateSharedPtr(B);
122 
123  for (k = 0; k < dim; ++k)
124  {
125  edgePhys [k] = Array<OneD, NekDouble>(nEdgePts);
126  edgeCoord[k] = Array<OneD, NekDouble>(nEdgePts);
127  exp->GetEdgePhysVals(j, seg, phys [k], edgePhys [k]);
128  exp->GetEdgePhysVals(j, seg, coord[k], edgeCoord[k]);
129  }
130 
131  // Update verts
132  for (k = 0; k < 2; ++k)
133  {
134  int id = edge->GetVid(k);
135  if (updatedVerts.find(id) != updatedVerts.end())
136  {
137  continue;
138  }
139 
141  edge->GetVertex(k);
142 
143  pt->UpdatePosition(
144  (*pt)(0) + edgePhys[0][k*(nEdgePts-1)],
145  (*pt)(1) + edgePhys[1][k*(nEdgePts-1)],
146  (*pt)(2));
147 
148  updatedVerts.insert(id);
149  }
150 
151  // Update curve
153  SpatialDomains::Curve>::AllocateSharedPtr(
154  edge->GetGlobalID(),
156 
157  for (k = 0; k < nEdgePts; ++k)
158  {
162  dim, edge->GetGlobalID(),
163  edgeCoord[0][k] + edgePhys[0][k],
164  edgeCoord[1][k] + edgePhys[1][k], 0.0);
165 
166  curve->m_points.push_back(vert);
167  }
168 
169  curvedEdges[edge->GetGlobalID()] = curve;
170  updatedEdges.insert(edge->GetGlobalID());
171  }
172  }
173  else if (dim == 3)
174  {
175  exp->GetCoords(coord[0], coord[1], coord[2]);
176 
178  boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
179  exp->GetGeom());
180 
181  for (j = 0; j < exp->GetNfaces(); ++j)
182  {
184 
185  // This edge has already been processed.
186  if (updatedFaces.find(face->GetGlobalID()) !=
187  updatedFaces.end())
188  {
189  continue;
190  }
191 
192  // Extract face displacement.
193  LibUtilities::BasisKey B0 = exp->DetFaceBasisKey(j,0);
194  LibUtilities::BasisKey B1 = exp->DetFaceBasisKey(j,1);
195  int nq0 = B0.GetNumPoints();
196  int nq1 = B1.GetNumPoints();
197 
200  B1.GetPointsType()
202  "Deformation requires GLL points in both "
203  "directions on a face.");
204 
206 
208  StdRegions::Orientation orient = exp->GetForient(j);
209 
210  if (face->GetShapeType() == LibUtilities::eTriangle)
211  {
213  AllocateSharedPtr(B0, B1);
214  }
215  else
216  {
218  AllocateSharedPtr(B0, B1);
219  }
220 
221  for (k = 0; k < dim; ++k)
222  {
223  Array<OneD, NekDouble> tmp(nq0*nq1);
224  newPos[k] = Array<OneD, NekDouble>(nq0*nq1);
225  exp->GetFacePhysVals(
226  j, faceexp, phys [k], tmp, orient);
227  exp->GetFacePhysVals(
228  j, faceexp, coord[k], newPos[k], orient);
229  Vmath::Vadd(
230  nq0*nq1, tmp, 1, newPos[k], 1, newPos[k], 1);
231  }
232 
233  // Now interpolate face onto a more reasonable set of
234  // points.
235  int nq = max(nq0, nq1);
236  LibUtilities::PointsKey edgePts(
240 
242 
243  for (k = 0; k < dim; ++k)
244  {
245  intPos[k] = Array<OneD, NekDouble>(nq*nq);
247  faceexp->GetPointsKeys()[0],
248  faceexp->GetPointsKeys()[1],
249  newPos[k], edgePts, edgePts, intPos[k]);
250  }
251 
252  int edgeOff[2][4][2] = {
253  {
254  {0, 1},
255  {nq-1, nq},
256  {nq*(nq-1), -nq},
257  {-1,-1}
258  },
259  {
260  {0, 1},
261  {nq-1, nq},
262  {nq*nq-1, -1},
263  {nq*(nq-1), -nq}
264  }
265  };
266 
267  for (k = 0; k < face->GetNumVerts(); ++k)
268  {
269  // Update verts
270  int id = face->GetVid(k);
271  const int o =
272  face->GetShapeType() - LibUtilities::eTriangle;
273 
274  if (updatedVerts.find(id) == updatedVerts.end())
275  {
277  face->GetVertex(k);
278  pt->UpdatePosition(
279  intPos[0][edgeOff[o][k][0]],
280  intPos[1][edgeOff[o][k][0]],
281  intPos[2][edgeOff[o][k][0]]);
282  updatedVerts.insert(id);
283  }
284 
285  // Update edges
286  id = face->GetEid(k);
287  if (updatedEdges.find(id) == updatedEdges.end())
288  {
290  = face->GetEdge(k);
294  edge->GetGlobalID(),
296 
297  const int offset = edgeOff[o][k][0];
298  const int pos = edgeOff[o][k][1];
299 
300  if (face->GetEorient(k) == StdRegions::eBackwards)
301  {
302  for (l = nq-1; l >= 0; --l)
303  {
304  int m = offset + pos*l;
308  dim, edge->GetGlobalID(),
309  intPos[0][m], intPos[1][m],
310  intPos[2][m]);
311  curve->m_points.push_back(vert);
312  }
313  }
314  else
315  {
316  for (l = 0; l < nq; ++l)
317  {
318  int m = offset + pos*l;
322  dim, edge->GetGlobalID(),
323  intPos[0][m], intPos[1][m],
324  intPos[2][m]);
325  curve->m_points.push_back(vert);
326  }
327  }
328 
329  curvedEdges[edge->GetGlobalID()] = curve;
330  updatedEdges.insert(edge->GetGlobalID());
331  }
332  }
333 
334  // Update face-interior curvature
336  face->GetShapeType() == LibUtilities::eTriangle ?
339 
341  SpatialDomains::Curve>::AllocateSharedPtr(
342  face->GetGlobalID(),
343  pType);
344 
345  if (face->GetShapeType() == LibUtilities::eTriangle)
346  {
347  // This code is probably pretty crappy. Have to go from
348  // GLL-GLL points -> GLL-Gauss-Radau -> nodal triangle
349  // points.
350  const LibUtilities::BasisKey B0(
354  const LibUtilities::BasisKey B1(
358  StdRegions::StdNodalTriExp nodalTri(B0, B1, pType);
359  StdRegions::StdTriExp tri (B0, B1);
360 
361  for (k = 0; k < dim; ++k)
362  {
363  Array<OneD, NekDouble> nodal(nq*nq);
364 
366  faceexp->GetBasis(0)->GetBasisKey(),
367  faceexp->GetBasis(1)->GetBasisKey(),
368  newPos[k], B0, B1, nodal);
369 
370  Array<OneD, NekDouble> tmp1(nq*(nq+1)/2);
371  Array<OneD, NekDouble> tmp2(nq*(nq+1)/2);
372 
373  tri.FwdTrans(nodal, tmp1);
374  nodalTri.ModalToNodal(tmp1, tmp2);
375  newPos[k] = tmp2;
376  }
377 
378  for (l = 0; l < nq*(nq+1)/2; ++l)
379  {
383  dim, face->GetGlobalID(),
384  newPos[0][l], newPos[1][l], newPos[2][l]);
385  curve->m_points.push_back(vert);
386  }
387  }
388  else
389  {
390  for (l = 0; l < nq*nq; ++l)
391  {
395  dim, face->GetGlobalID(),
396  intPos[0][l], intPos[1][l], intPos[2][l]);
397  curve->m_points.push_back(vert);
398  }
399  }
400 
401  curvedFaces[face->GetGlobalID()] = curve;
402  updatedFaces.insert(face->GetGlobalID());
403  }
404  }
405  }
406 
407  // Reset geometry information
408  for (i = 0; i < fields.num_elements(); ++i)
409  {
410  fields[i]->Reset();
411  }
412  }
413 }
414 }
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Principle Modified Functions .
Definition: BasisType.h:49
Geometry2DSharedPtr GetFace(int i)
Return face i in this element.
Definition: Geometry3D.cpp:79
PointsType GetPointsType() const
Return type of quadrature.
Definition: Basis.h:151
boost::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:62
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Principle Orthogonal Functions .
Definition: BasisType.h:47
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:116
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
void UpdateGeometry(SpatialDomains::MeshGraphSharedPtr graph, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Update geometry according to displacement that is in current fields.
Definition: Deform.cpp:57
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:128
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
3D geometry information
Definition: Geometry3D.h:70
2D geometry information
Definition: Geometry2D.h:65
const Geometry1DSharedPtr GetEdge(int i) const
Definition: Geometry2D.cpp:80
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
boost::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:63
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
Describes the specification for a Basis.
Definition: Basis.h:50
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:68
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...