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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Deformation of mesh from fields.
32//
33///////////////////////////////////////////////////////////////////////////////
34
42
44{
45
46/**
47 * @brief Update geometry according to displacement that is in current
48 * fields.
49 *
50 * @param graph The MeshGraph of the current geometry.
51 * @param fields The fields containing the displacement.
52 */
55 Array<OneD, Array<OneD, NekDouble>> &PhysVals, bool modal)
56{
57 // Clear existing curvature.
58 SpatialDomains::CurveMap &curvedEdges = graph->GetCurvedEdges();
59 SpatialDomains::CurveMap &curvedFaces = graph->GetCurvedFaces();
60 curvedEdges.clear();
61 curvedFaces.clear();
62
63 int i, j, k, l, dim;
64
65 // Sets to hold IDs of updated vertices to avoid duplicating effort.
66 std::set<int> updatedVerts, updatedEdges, updatedFaces;
67
68 dim = graph->GetSpaceDimension();
71
72 for (i = 0; i < fields[0]->GetExpSize(); ++i)
73 {
74 LocalRegions::ExpansionSharedPtr exp = fields[0]->GetExp(i);
75 int offset = fields[0]->GetPhys_Offset(i);
76 int nquad = exp->GetTotPoints();
77
78 // Extract displacement for this element, allocate storage for
79 // elemental coordinates.
80 for (j = 0; j < dim; ++j)
81 {
82 phys[j] = Array<OneD, NekDouble>(nquad, PhysVals[j] + offset);
83 coord[j] = Array<OneD, NekDouble>(nquad);
84 }
85
86 // In 2D loop over edges.
87 if (dim == 2)
88 {
89 exp->GetCoords(coord[0], coord[1]);
90
92 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
93 exp->GetGeom());
94
95 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
96 {
97 SpatialDomains::Geometry1DSharedPtr edge = geom->GetEdge(j);
98
99 // This edge has already been processed.
100 if (updatedEdges.find(edge->GetGlobalID()) !=
101 updatedEdges.end())
102 {
103 continue;
104 }
105
106 // Extract edge displacement.
107 int nEdgePts = exp->GetTraceNumPoints(j);
110
117
118 for (k = 0; k < dim; ++k)
119 {
120 edgePhys[k] = Array<OneD, NekDouble>(nEdgePts);
121 edgeCoord[k] = Array<OneD, NekDouble>(nEdgePts);
122 exp->GetTracePhysVals(j, seg, phys[k], edgePhys[k]);
123 exp->GetTracePhysVals(j, seg, coord[k], edgeCoord[k]);
124 }
125
126 // Update verts
127 for (k = 0; k < 2; ++k)
128 {
129 int id = edge->GetVid(k);
130 if (updatedVerts.find(id) != updatedVerts.end())
131 {
132 continue;
133 }
134
135 SpatialDomains::PointGeomSharedPtr pt = edge->GetVertex(k);
136
137 pt->UpdatePosition(
138 (*pt)(0) + edgePhys[0][k * (nEdgePts - 1)],
139 (*pt)(1) + edgePhys[1][k * (nEdgePts - 1)], (*pt)(2));
140
141 updatedVerts.insert(id);
142 }
143
144 // Update curve
147 edge->GetGlobalID(),
149
150 for (k = 0; k < nEdgePts; ++k)
151 {
154 AllocateSharedPtr(dim, edge->GetGlobalID(),
155 edgeCoord[0][k] + edgePhys[0][k],
156 edgeCoord[1][k] + edgePhys[1][k],
157 0.0);
158
159 curve->m_points.push_back(vert);
160 }
161
162 curvedEdges[edge->GetGlobalID()] = curve;
163 updatedEdges.insert(edge->GetGlobalID());
164 }
165 }
166 else if (dim == 3)
167 {
168 exp->GetCoords(coord[0], coord[1], coord[2]);
169
171 std::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
172 exp->GetGeom());
173
174 for (j = 0; j < exp->GetNtraces(); ++j)
175 {
176 SpatialDomains::Geometry2DSharedPtr face = geom->GetFace(j);
177
179 exp->as<LocalRegions::Expansion3D>();
180
181 // This edge has already been processed.
182 if (updatedFaces.find(face->GetGlobalID()) !=
183 updatedFaces.end())
184 {
185 continue;
186 }
187
188 // Extract face displacement.
189 LibUtilities::BasisKey B0 = exp->GetTraceBasisKey(j, 0);
190 LibUtilities::BasisKey B1 = exp->GetTraceBasisKey(j, 1);
191 int nq0 = B0.GetNumPoints();
192 int nq1 = B1.GetNumPoints();
193
194 ASSERTL1(B0.GetPointsType() ==
196 B1.GetPointsType() ==
198 "Deformation requires GLL points in both "
199 "directions on a face.");
200
202
204 StdRegions::Orientation orient = exp->GetTraceOrient(j);
205
206 if (face->GetShapeType() == LibUtilities::eTriangle)
207 {
208 faceexp =
210 B0, B1);
211 }
212 else
213 {
214 faceexp = MemoryManager<
215 StdRegions::StdQuadExp>::AllocateSharedPtr(B0, B1);
216 }
217
218 for (k = 0; k < dim; ++k)
219 {
220 Array<OneD, NekDouble> tmp(nq0 * nq1);
221 newPos[k] = Array<OneD, NekDouble>(nq0 * nq1);
222 exp3d->GetTracePhysVals(j, faceexp, phys[k], tmp, orient);
223 exp3d->GetTracePhysVals(j, faceexp, coord[k], newPos[k],
224 orient);
225 Vmath::Vadd(nq0 * nq1, tmp, 1, newPos[k], 1, newPos[k], 1);
226 }
227
228 // Now interpolate face onto a more reasonable set of
229 // points.
230 int nq = std::max(nq0, nq1);
231 if (!modal)
232 {
233 nq--;
234 }
235
239
241
242 for (k = 0; k < dim; ++k)
243 {
244 intPos[k] = Array<OneD, NekDouble>(nq * nq);
245 LibUtilities::Interp2D(faceexp->GetPointsKeys()[0],
246 faceexp->GetPointsKeys()[1],
247 newPos[k], edgePts, edgePts,
248 intPos[k]);
249 }
250
251 int edgeOff[2][4][2] = {
252 {{0, 1}, {nq - 1, nq}, {nq * (nq - 1), -nq}, {-1, -1}},
253 {{0, 1},
254 {nq - 1, nq},
255 {nq * nq - 1, -1},
256 {nq * (nq - 1), -nq}}};
257
258 for (k = 0; k < face->GetNumVerts(); ++k)
259 {
260 // Update verts
261 int id = face->GetVid(k);
262 const int o =
263 face->GetShapeType() - LibUtilities::eTriangle;
264
265 if (updatedVerts.find(id) == updatedVerts.end())
266 {
268 face->GetVertex(k);
269 pt->UpdatePosition(intPos[0][edgeOff[o][k][0]],
270 intPos[1][edgeOff[o][k][0]],
271 intPos[2][edgeOff[o][k][0]]);
272 updatedVerts.insert(id);
273 }
274
275 // Update edges
276 id = face->GetEid(k);
277 if (updatedEdges.find(id) == updatedEdges.end())
278 {
280 face->GetEdge(k);
284 edge->GetGlobalID(),
286
287 const int offset = edgeOff[o][k][0];
288 const int pos = edgeOff[o][k][1];
289
290 if (face->GetEorient(k) == StdRegions::eBackwards)
291 {
292 for (l = nq - 1; l >= 0; --l)
293 {
294 int m = offset + pos * l;
298 dim, edge->GetGlobalID(),
299 intPos[0][m], intPos[1][m],
300 intPos[2][m]);
301 curve->m_points.push_back(vert);
302 }
303 }
304 else
305 {
306 for (l = 0; l < nq; ++l)
307 {
308 int m = offset + pos * l;
312 dim, edge->GetGlobalID(),
313 intPos[0][m], intPos[1][m],
314 intPos[2][m]);
315 curve->m_points.push_back(vert);
316 }
317 }
318
319 curvedEdges[edge->GetGlobalID()] = curve;
320 updatedEdges.insert(edge->GetGlobalID());
321 }
322 }
323
324 // Update face-interior curvature
326 face->GetShapeType() == LibUtilities::eTriangle
329
332 face->GetGlobalID(), pType);
333
334 if (face->GetShapeType() == LibUtilities::eTriangle)
335 {
336 // This code is probably pretty crappy. Have to go from
337 // GLL-GLL points -> GLL-Gauss-Radau -> nodal triangle
338 // points.
339 const LibUtilities::BasisKey B0(
343 const LibUtilities::BasisKey B1(
346 nq, LibUtilities::eGaussRadauMAlpha1Beta0));
347 StdRegions::StdNodalTriExp nodalTri(B0, B1, pType);
348 StdRegions::StdTriExp tri(B0, B1);
349
350 for (k = 0; k < dim; ++k)
351 {
352 Array<OneD, NekDouble> nodal(nq * nq);
353
355 faceexp->GetBasis(0)->GetBasisKey(),
356 faceexp->GetBasis(1)->GetBasisKey(), newPos[k], B0,
357 B1, nodal);
358
359 Array<OneD, NekDouble> tmp1(nq * (nq + 1) / 2);
360 Array<OneD, NekDouble> tmp2(nq * (nq + 1) / 2);
361
362 tri.FwdTrans(nodal, tmp1);
363 nodalTri.ModalToNodal(tmp1, tmp2);
364 newPos[k] = tmp2;
365 }
366
367 for (l = 0; l < nq * (nq + 1) / 2; ++l)
368 {
371 AllocateSharedPtr(dim, face->GetGlobalID(),
372 newPos[0][l], newPos[1][l],
373 newPos[2][l]);
374 curve->m_points.push_back(vert);
375 }
376 }
377 else
378 {
379 for (l = 0; l < nq * nq; ++l)
380 {
383 AllocateSharedPtr(dim, face->GetGlobalID(),
384 intPos[0][l], intPos[1][l],
385 intPos[2][l]);
386 curve->m_points.push_back(vert);
387 }
388 }
389
390 curvedFaces[face->GetGlobalID()] = curve;
391 updatedFaces.insert(face->GetGlobalID());
392 }
393 }
394 }
395
396 // Reset geometry information
397 for (i = 0; i < fields.size(); ++i)
398 {
399 fields[i]->Reset();
400 }
401}
402} // namespace Nektar::GlobalMapping
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
Describes the specification for a Basis.
Definition: Basis.h:45
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:120
PointsType GetPointsType() const
Return type of quadrature.
Definition: Basis.h:143
Defines a specification for a set of points.
Definition: Points.h:50
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void UpdateGeometry(SpatialDomains::MeshGraphSharedPtr graph, Array< OneD, MultiRegions::ExpListSharedPtr > &fields, Array< OneD, Array< OneD, NekDouble > > &PhysVals, bool modal)
Update geometry according to displacement that is in current fields.
Definition: Deform.cpp:53
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:101
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:81
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:47
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:58
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:59
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
std::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:50
std::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
std::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
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.hpp:180