Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 
40 #include <GlobalMapping/Deform.h>
42 #include <StdRegions/StdSegExp.h>
43 #include <StdRegions/StdQuadExp.h>
45 #include <MultiRegions/ExpList.h>
46 
47 using namespace std;
48 
49 namespace Nektar {
50 namespace GlobalMapping {
51 
52  /**
53  * @brief Update geometry according to displacement that is in current
54  * fields.
55  *
56  * @param graph The MeshGraph of the current geometry.
57  * @param fields The fields containing the displacement.
58  */
62  bool modal)
63  {
64  // Clear existing curvature.
65  SpatialDomains::CurveMap &curvedEdges = graph->GetCurvedEdges();
66  SpatialDomains::CurveMap &curvedFaces = graph->GetCurvedFaces();
67  curvedEdges.clear();
68  curvedFaces.clear();
69 
70  int i, j, k, l, dim;
71 
72  // Sets to hold IDs of updated vertices to avoid duplicating effort.
73  set<int> updatedVerts, updatedEdges, updatedFaces;
74 
75  dim = graph->GetSpaceDimension();
78 
79  for (i = 0; i < fields[0]->GetExpSize(); ++i)
80  {
81  LocalRegions::ExpansionSharedPtr exp = fields[0]->GetExp(i);
82  int offset = fields[0]->GetPhys_Offset(i);
83  int nquad = exp->GetTotPoints();
84 
85  // Extract displacement for this element, allocate storage for
86  // elemental coordinates.
87  for (j = 0; j < dim; ++j)
88  {
89  phys[j] = Array<OneD, NekDouble>(
90  nquad, fields[j]->UpdatePhys() + offset);
91  coord[j] = Array<OneD, NekDouble>(nquad);
92  }
93 
94  // In 2D loop over edges.
95  if (dim == 2)
96  {
97  exp->GetCoords(coord[0], coord[1]);
98 
100  boost::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
101  exp->GetGeom());
102 
103  for (j = 0; j < exp->GetNedges(); ++j)
104  {
106 
107  // This edge has already been processed.
108  if (updatedEdges.find(edge->GetGlobalID()) !=
109  updatedEdges.end())
110  {
111  continue;
112  }
113 
114  // Extract edge displacement.
115  int nEdgePts = exp->GetEdgeNumPoints(j);
116  Array<OneD, Array<OneD, NekDouble> > edgePhys (dim);
117  Array<OneD, Array<OneD, NekDouble> > edgeCoord(dim);
118 
119  const LibUtilities::BasisKey B(
120  LibUtilities::eModified_A, nEdgePts,
124  StdRegions::StdSegExp>::AllocateSharedPtr(B);
125 
126  for (k = 0; k < dim; ++k)
127  {
128  edgePhys [k] = Array<OneD, NekDouble>(nEdgePts);
129  edgeCoord[k] = Array<OneD, NekDouble>(nEdgePts);
130  exp->GetEdgePhysVals(j, seg, phys [k], edgePhys [k]);
131  exp->GetEdgePhysVals(j, seg, coord[k], edgeCoord[k]);
132  }
133 
134  // Update verts
135  for (k = 0; k < 2; ++k)
136  {
137  int id = edge->GetVid(k);
138  if (updatedVerts.find(id) != updatedVerts.end())
139  {
140  continue;
141  }
142 
144  edge->GetVertex(k);
145 
146  pt->UpdatePosition(
147  (*pt)(0) + edgePhys[0][k*(nEdgePts-1)],
148  (*pt)(1) + edgePhys[1][k*(nEdgePts-1)],
149  (*pt)(2));
150 
151  updatedVerts.insert(id);
152  }
153 
154  // Update curve
156  SpatialDomains::Curve>::AllocateSharedPtr(
157  edge->GetGlobalID(),
159 
160  for (k = 0; k < nEdgePts; ++k)
161  {
165  dim, edge->GetGlobalID(),
166  edgeCoord[0][k] + edgePhys[0][k],
167  edgeCoord[1][k] + edgePhys[1][k], 0.0);
168 
169  curve->m_points.push_back(vert);
170  }
171 
172  curvedEdges[edge->GetGlobalID()] = curve;
173  updatedEdges.insert(edge->GetGlobalID());
174  }
175  }
176  else if (dim == 3)
177  {
178  exp->GetCoords(coord[0], coord[1], coord[2]);
179 
181  boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
182  exp->GetGeom());
183 
184  for (j = 0; j < exp->GetNfaces(); ++j)
185  {
187 
188  // This edge has already been processed.
189  if (updatedFaces.find(face->GetGlobalID()) !=
190  updatedFaces.end())
191  {
192  continue;
193  }
194 
195  // Extract face displacement.
196  LibUtilities::BasisKey B0 = exp->DetFaceBasisKey(j,0);
197  LibUtilities::BasisKey B1 = exp->DetFaceBasisKey(j,1);
198  int nq0 = B0.GetNumPoints();
199  int nq1 = B1.GetNumPoints();
200 
203  B1.GetPointsType()
205  "Deformation requires GLL points in both "
206  "directions on a face.");
207 
209 
211  StdRegions::Orientation orient = exp->GetForient(j);
212 
213  if (face->GetShapeType() == LibUtilities::eTriangle)
214  {
216  AllocateSharedPtr(B0, B1);
217  }
218  else
219  {
221  AllocateSharedPtr(B0, B1);
222  }
223 
224  for (k = 0; k < dim; ++k)
225  {
226  Array<OneD, NekDouble> tmp(nq0*nq1);
227  newPos[k] = Array<OneD, NekDouble>(nq0*nq1);
228  exp->GetFacePhysVals(
229  j, faceexp, phys [k], tmp, orient);
230  exp->GetFacePhysVals(
231  j, faceexp, coord[k], newPos[k], orient);
232  Vmath::Vadd(
233  nq0*nq1, tmp, 1, newPos[k], 1, newPos[k], 1);
234  }
235 
236  // Now interpolate face onto a more reasonable set of
237  // points.
238  int nq = max(nq0, nq1);
239  if(!modal)
240  nq--;
241 
242  LibUtilities::PointsKey edgePts(
246 
248 
249  for (k = 0; k < dim; ++k)
250  {
251  intPos[k] = Array<OneD, NekDouble>(nq*nq);
253  faceexp->GetPointsKeys()[0],
254  faceexp->GetPointsKeys()[1],
255  newPos[k], edgePts, edgePts, intPos[k]);
256  }
257 
258  int edgeOff[2][4][2] = {
259  {
260  {0, 1},
261  {nq-1, nq},
262  {nq*(nq-1), -nq},
263  {-1,-1}
264  },
265  {
266  {0, 1},
267  {nq-1, nq},
268  {nq*nq-1, -1},
269  {nq*(nq-1), -nq}
270  }
271  };
272 
273  for (k = 0; k < face->GetNumVerts(); ++k)
274  {
275  // Update verts
276  int id = face->GetVid(k);
277  const int o =
278  face->GetShapeType() - LibUtilities::eTriangle;
279 
280  if (updatedVerts.find(id) == updatedVerts.end())
281  {
283  face->GetVertex(k);
284  pt->UpdatePosition(
285  intPos[0][edgeOff[o][k][0]],
286  intPos[1][edgeOff[o][k][0]],
287  intPos[2][edgeOff[o][k][0]]);
288  updatedVerts.insert(id);
289  }
290 
291  // Update edges
292  id = face->GetEid(k);
293  if (updatedEdges.find(id) == updatedEdges.end())
294  {
296  = face->GetEdge(k);
300  edge->GetGlobalID(),
302 
303  const int offset = edgeOff[o][k][0];
304  const int pos = edgeOff[o][k][1];
305 
306  if (face->GetEorient(k) == StdRegions::eBackwards)
307  {
308  for (l = nq-1; l >= 0; --l)
309  {
310  int m = offset + pos*l;
314  dim, edge->GetGlobalID(),
315  intPos[0][m], intPos[1][m],
316  intPos[2][m]);
317  curve->m_points.push_back(vert);
318  }
319  }
320  else
321  {
322  for (l = 0; l < nq; ++l)
323  {
324  int m = offset + pos*l;
328  dim, edge->GetGlobalID(),
329  intPos[0][m], intPos[1][m],
330  intPos[2][m]);
331  curve->m_points.push_back(vert);
332  }
333  }
334 
335  curvedEdges[edge->GetGlobalID()] = curve;
336  updatedEdges.insert(edge->GetGlobalID());
337  }
338  }
339 
340  // Update face-interior curvature
342  face->GetShapeType() == LibUtilities::eTriangle ?
345 
347  SpatialDomains::Curve>::AllocateSharedPtr(
348  face->GetGlobalID(),
349  pType);
350 
351  if (face->GetShapeType() == LibUtilities::eTriangle)
352  {
353  // This code is probably pretty crappy. Have to go from
354  // GLL-GLL points -> GLL-Gauss-Radau -> nodal triangle
355  // points.
356  const LibUtilities::BasisKey B0(
360  const LibUtilities::BasisKey B1(
364  StdRegions::StdNodalTriExp nodalTri(B0, B1, pType);
365  StdRegions::StdTriExp tri (B0, B1);
366 
367  for (k = 0; k < dim; ++k)
368  {
369  Array<OneD, NekDouble> nodal(nq*nq);
370 
372  faceexp->GetBasis(0)->GetBasisKey(),
373  faceexp->GetBasis(1)->GetBasisKey(),
374  newPos[k], B0, B1, nodal);
375 
376  Array<OneD, NekDouble> tmp1(nq*(nq+1)/2);
377  Array<OneD, NekDouble> tmp2(nq*(nq+1)/2);
378 
379  tri.FwdTrans(nodal, tmp1);
380  nodalTri.ModalToNodal(tmp1, tmp2);
381  newPos[k] = tmp2;
382  }
383 
384  for (l = 0; l < nq*(nq+1)/2; ++l)
385  {
389  dim, face->GetGlobalID(),
390  newPos[0][l], newPos[1][l], newPos[2][l]);
391  curve->m_points.push_back(vert);
392  }
393  }
394  else
395  {
396  for (l = 0; l < nq*nq; ++l)
397  {
401  dim, face->GetGlobalID(),
402  intPos[0][l], intPos[1][l], intPos[2][l]);
403  curve->m_points.push_back(vert);
404  }
405  }
406 
407  curvedFaces[face->GetGlobalID()] = curve;
408  updatedFaces.insert(face->GetGlobalID());
409  }
410  }
411  }
412 
413  // Reset geometry information
414  for (i = 0; i < fields.num_elements(); ++i)
415  {
416  fields[i]->Reset();
417  }
418  }
419 }
420 }
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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:82
STL namespace.
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:59
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
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< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void UpdateGeometry(SpatialDomains::MeshGraphSharedPtr graph, Array< OneD, MultiRegions::ExpListSharedPtr > &fields, bool modal)
Update geometry according to displacement that is in current fields.
Definition: Deform.cpp:59
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:83
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:228
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
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:52
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:299
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:70
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...