Update geometry according to displacement that is in current fields.
Adds a summary item to the summary info list.
62{
63
66 curvedEdges.clear();
67 curvedFaces.clear();
68
69 int i, j, k, l, dim;
70
71
72 set<int> updatedVerts, updatedEdges, updatedFaces;
73
74 dim = graph->GetSpaceDimension();
77
78 for (i = 0; i < fields[0]->GetExpSize(); ++i)
79 {
81 int offset = fields[0]->GetPhys_Offset(i);
82 int nquad = exp->GetTotPoints();
83
84
85
86 for (j = 0; j < dim; ++j)
87 {
90 }
91
92
93 if (dim == 2)
94 {
95 exp->GetCoords(coord[0], coord[1]);
96
98 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
99 exp->GetGeom());
100
101 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
102 {
104
105
106 if (updatedEdges.find(edge->GetGlobalID()) !=
107 updatedEdges.end())
108 {
109 continue;
110 }
111
112
113 int nEdgePts = exp->GetTraceNumPoints(j);
114 Array<OneD, Array<OneD, NekDouble>> edgePhys(dim);
115 Array<OneD, Array<OneD, NekDouble>> edgeCoord(dim);
116
117 const LibUtilities::BasisKey B(
119 LibUtilities::PointsKey(
122 MemoryManager<StdRegions::StdSegExp>::AllocateSharedPtr(B);
123
124 for (k = 0; k < dim; ++k)
125 {
126 edgePhys[k] = Array<OneD, NekDouble>(nEdgePts);
127 edgeCoord[k] = Array<OneD, NekDouble>(nEdgePts);
128 exp->GetTracePhysVals(j, seg, phys[k], edgePhys[k]);
129 exp->GetTracePhysVals(j, seg, coord[k], edgeCoord[k]);
130 }
131
132
133 for (k = 0; k < 2; ++k)
134 {
135 int id = edge->GetVid(k);
136 if (updatedVerts.find(id) != updatedVerts.end())
137 {
138 continue;
139 }
140
142
143 pt->UpdatePosition(
144 (*pt)(0) + edgePhys[0][k * (nEdgePts - 1)],
145 (*pt)(1) + edgePhys[1][k * (nEdgePts - 1)], (*pt)(2));
146
147 updatedVerts.insert(id);
148 }
149
150
152 MemoryManager<SpatialDomains::Curve>::AllocateSharedPtr(
153 edge->GetGlobalID(),
155
156 for (k = 0; k < nEdgePts; ++k)
157 {
159 MemoryManager<SpatialDomains::PointGeom>::
160 AllocateSharedPtr(dim, edge->GetGlobalID(),
161 edgeCoord[0][k] + edgePhys[0][k],
162 edgeCoord[1][k] + edgePhys[1][k],
163 0.0);
164
165 curve->m_points.push_back(vert);
166 }
167
168 curvedEdges[edge->GetGlobalID()] = curve;
169 updatedEdges.insert(edge->GetGlobalID());
170 }
171 }
172 else if (dim == 3)
173 {
174 exp->GetCoords(coord[0], coord[1], coord[2]);
175
177 std::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
178 exp->GetGeom());
179
180 for (j = 0; j < exp->GetNtraces(); ++j)
181 {
183
185 exp->as<LocalRegions::Expansion3D>();
186
187
188 if (updatedFaces.find(face->GetGlobalID()) !=
189 updatedFaces.end())
190 {
191 continue;
192 }
193
194
195 LibUtilities::BasisKey B0 = exp->GetTraceBasisKey(j, 0);
196 LibUtilities::BasisKey B1 = exp->GetTraceBasisKey(j, 1);
197 int nq0 = B0.GetNumPoints();
198 int nq1 = B1.GetNumPoints();
199
202 B1.GetPointsType() ==
204 "Deformation requires GLL points in both "
205 "directions on a face.");
206
207 Array<OneD, Array<OneD, NekDouble>> newPos(dim);
208
211
213 {
214 faceexp =
215 MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(
216 B0, B1);
217 }
218 else
219 {
220 faceexp = MemoryManager<
221 StdRegions::StdQuadExp>::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 exp3d->GetTracePhysVals(j, faceexp, phys[k], tmp, orient);
229 exp3d->GetTracePhysVals(j, faceexp, coord[k], newPos[k],
230 orient);
231 Vmath::Vadd(nq0 * nq1, tmp, 1, newPos[k], 1, newPos[k], 1);
232 }
233
234
235
236 int nq = max(nq0, nq1);
237 if (!modal)
238 {
239 nq--;
240 }
241
242 LibUtilities::PointsKey edgePts(
245
246 Array<OneD, Array<OneD, NekDouble>> intPos(dim);
247
248 for (k = 0; k < dim; ++k)
249 {
250 intPos[k] = Array<OneD, NekDouble>(nq * nq);
252 faceexp->GetPointsKeys()[1],
253 newPos[k], edgePts, edgePts,
254 intPos[k]);
255 }
256
257 int edgeOff[2][4][2] = {
258 {{0, 1}, {nq - 1, nq}, {nq * (nq - 1), -nq}, {-1, -1}},
259 {{0, 1},
260 {nq - 1, nq},
261 {nq * nq - 1, -1},
262 {nq * (nq - 1), -nq}}};
263
264 for (k = 0; k < face->GetNumVerts(); ++k)
265 {
266
267 int id = face->GetVid(k);
268 const int o =
270
271 if (updatedVerts.find(id) == updatedVerts.end())
272 {
274 face->GetVertex(k);
275 pt->UpdatePosition(intPos[0][edgeOff[o][k][0]],
276 intPos[1][edgeOff[o][k][0]],
277 intPos[2][edgeOff[o][k][0]]);
278 updatedVerts.insert(id);
279 }
280
281
282 id = face->GetEid(k);
283 if (updatedEdges.find(id) == updatedEdges.end())
284 {
286 face->GetEdge(k);
288 MemoryManager<SpatialDomains::Curve>::
289 AllocateSharedPtr(
290 edge->GetGlobalID(),
292
293 const int offset = edgeOff[o][k][0];
294 const int pos = edgeOff[o][k][1];
295
297 {
298 for (l = nq - 1; l >= 0; --l)
299 {
300 int m = offset + pos * l;
302 MemoryManager<SpatialDomains::PointGeom>::
303 AllocateSharedPtr(
304 dim, edge->GetGlobalID(),
305 intPos[0][m], intPos[1][m],
306 intPos[2][m]);
307 curve->m_points.push_back(vert);
308 }
309 }
310 else
311 {
312 for (l = 0; l < nq; ++l)
313 {
314 int m = offset + pos * l;
316 MemoryManager<SpatialDomains::PointGeom>::
317 AllocateSharedPtr(
318 dim, edge->GetGlobalID(),
319 intPos[0][m], intPos[1][m],
320 intPos[2][m]);
321 curve->m_points.push_back(vert);
322 }
323 }
324
325 curvedEdges[edge->GetGlobalID()] = curve;
326 updatedEdges.insert(edge->GetGlobalID());
327 }
328 }
329
330
335
337 MemoryManager<SpatialDomains::Curve>::AllocateSharedPtr(
338 face->GetGlobalID(), pType);
339
341 {
342
343
344
345 const LibUtilities::BasisKey B0(
347 LibUtilities::PointsKey(
349 const LibUtilities::BasisKey B1(
351 LibUtilities::PointsKey(
352 nq, LibUtilities::eGaussRadauMAlpha1Beta0));
353 StdRegions::StdNodalTriExp nodalTri(B0, B1, pType);
354 StdRegions::StdTriExp tri(B0, B1);
355
356 for (k = 0; k < dim; ++k)
357 {
358 Array<OneD, NekDouble> nodal(nq * nq);
359
361 faceexp->GetBasis(0)->GetBasisKey(),
362 faceexp->GetBasis(1)->GetBasisKey(), newPos[k], B0,
363 B1, nodal);
364
365 Array<OneD, NekDouble> tmp1(nq * (nq + 1) / 2);
366 Array<OneD, NekDouble> tmp2(nq * (nq + 1) / 2);
367
368 tri.FwdTrans(nodal, tmp1);
369 nodalTri.ModalToNodal(tmp1, tmp2);
370 newPos[k] = tmp2;
371 }
372
373 for (l = 0; l < nq * (nq + 1) / 2; ++l)
374 {
376 MemoryManager<SpatialDomains::PointGeom>::
377 AllocateSharedPtr(dim, face->GetGlobalID(),
378 newPos[0][l], newPos[1][l],
379 newPos[2][l]);
380 curve->m_points.push_back(vert);
381 }
382 }
383 else
384 {
385 for (l = 0; l < nq * nq; ++l)
386 {
388 MemoryManager<SpatialDomains::PointGeom>::
389 AllocateSharedPtr(dim, face->GetGlobalID(),
390 intPos[0][l], intPos[1][l],
391 intPos[2][l]);
392 curve->m_points.push_back(vert);
393 }
394 }
395
396 curvedFaces[face->GetGlobalID()] = curve;
397 updatedFaces.insert(face->GetGlobalID());
398 }
399 }
400 }
401
402
403 for (i = 0; i < fields.size(); ++i)
404 {
405 fields[i]->Reset();
406 }
407}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
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,...
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
@ eOrtho_A
Principle Orthogonal Functions .
@ eOrtho_B
Principle Orthogonal Functions .
@ eModified_A
Principle Modified Functions .
std::shared_ptr< Expansion > ExpansionSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
std::shared_ptr< Curve > CurveSharedPtr
std::unordered_map< int, CurveSharedPtr > CurveMap
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
std::shared_ptr< Geometry3D > Geometry3DSharedPtr
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.