Nektar++
Movement/Zones.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: Zones.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: Zones used in the non-conformal interfaces
33// and ALE implementations
34//
35////////////////////////////////////////////////////////////////////////////////
36
40#include <tinyxml.h>
41
43{
44
45ZoneBase::ZoneBase(MovementType type, int indx, int domainID,
46 CompositeMap domain, int coordDim)
47 : m_type(type), m_id(indx), m_domainID(domainID), m_domain(domain),
48 m_coordDim(coordDim)
49{
50 for (auto &comp : domain)
51 {
52 for (auto &geom : comp.second->m_geomVec)
53 {
54 m_elements.emplace_back(geom);
55
56 // Fill constituent elements (i.e. faces and edges)
57 switch (geom->GetShapeDim())
58 {
59 case 3:
60 for (int i = 0; i < geom->GetNumFaces(); ++i)
61 {
62 m_constituentElements[2].insert(geom->GetFace(i));
63 }
64 /* fall through */
65 case 2:
66 for (int i = 0; i < geom->GetNumEdges(); ++i)
67 {
68 m_constituentElements[1].insert(geom->GetEdge(i));
69 }
70 /* fall through */
71 case 1:
72 m_constituentElements[0].insert(geom);
73 /* fall through */
74 default:
75 break;
76 }
77 }
78 }
79
80 // The seenVerts/Edges/Faces keeps track of geometry so duplicates aren't
81 // emplaced in to the storage vector
82 std::unordered_set<int> seenVerts, seenEdges, seenFaces;
83
84 // Fill verts/edges/faces vector storage that are to be moved each timestep
85 for (auto &comp : m_domain)
86 {
87 for (auto &geom : comp.second->m_geomVec)
88 {
89 for (int i = 0; i < geom->GetNumVerts(); ++i)
90 {
91 PointGeomSharedPtr vert = geom->GetVertex(i);
92
93 if (seenVerts.find(vert->GetGlobalID()) != seenVerts.end())
94 {
95 continue;
96 }
97
98 seenVerts.insert(vert->GetGlobalID());
99 m_verts.emplace_back(vert);
100 }
101
102 for (int i = 0; i < geom->GetNumEdges(); ++i)
103 {
104 SegGeomSharedPtr edge =
105 std::static_pointer_cast<SegGeom>(geom->GetEdge(i));
106
107 if (seenEdges.find(edge->GetGlobalID()) != seenEdges.end())
108 {
109 continue;
110 }
111
112 seenEdges.insert(edge->GetGlobalID());
113
114 CurveSharedPtr curve = edge->GetCurve();
115 if (!curve)
116 {
117 continue;
118 }
119
120 m_curves.emplace_back(curve);
121 }
122
123 for (int i = 0; i < geom->GetNumFaces(); ++i)
124 {
126 std::static_pointer_cast<Geometry2D>(geom->GetFace(i));
127
128 if (seenFaces.find(face->GetGlobalID()) != seenFaces.end())
129 {
130 continue;
131 }
132
133 seenFaces.insert(face->GetGlobalID());
134
135 CurveSharedPtr curve = face->GetCurve();
136 if (!curve)
137 {
138 continue;
139 }
140
141 m_curves.emplace_back(curve);
142 }
143 }
144 }
145
146 // Copy points so we know original positions.
147 for (auto &pt : m_verts)
148 {
149 m_origVerts.emplace_back(*pt);
150 }
151
152 for (auto &curve : m_curves)
153 {
154 for (auto &pt : curve->m_points)
155 {
156 m_origVerts.emplace_back(*pt);
157 }
158 }
159}
160
161ZoneRotate::ZoneRotate(int id, int domainID, const CompositeMap &domain,
162 const int coordDim, const NekPoint<NekDouble> &origin,
163 const DNekVec &axis,
164 const LibUtilities::EquationSharedPtr &angularVelEqn)
165 : ZoneBase(MovementType::eRotate, id, domainID, domain, coordDim),
166 m_origin(origin), m_axis(axis), m_angularVelEqn(angularVelEqn)
167{
168 // Construct rotation matrix
169 m_W(0, 1) = -m_axis[2];
170 m_W(0, 2) = m_axis[1];
171 m_W(1, 0) = m_axis[2];
172 m_W(1, 2) = -m_axis[0];
173 m_W(2, 0) = -m_axis[1];
174 m_W(2, 1) = m_axis[0];
175
176 m_W2 = m_W * m_W;
177}
178
180{
181 // Clear bboxes (these will be regenerated next time GetBoundingBox is
182 // called)
183 for (auto &el : m_elements)
184 {
185 el->ClearBoundingBox();
186
187 int nfaces = el->GetNumFaces();
188 for (int i = 0; i < nfaces; ++i)
189 {
190 el->GetFace(i)->ClearBoundingBox();
191 }
192
193 int nedges = el->GetNumEdges();
194 for (int i = 0; i < nedges; ++i)
195 {
196 el->GetEdge(i)->ClearBoundingBox();
197 }
198 }
199}
200
202{
203 NekDouble rampTime = 1;
204 if (time < rampTime)
205 {
206 return m_angularVelEqn->Evaluate(0, 0, 0, time) * time / rampTime;
207 }
208 else
209 {
210 return m_angularVelEqn->Evaluate(0, 0, 0, time);
211 }
212}
213
214// Calculate new location of points using Rodrigues formula
216{
217 // This works for a linear ramp
218 NekDouble rampTime = 1;
219 NekDouble angle;
220 if (time < rampTime)
221 {
222 angle = GetAngularVel(time) * (time / rampTime) / 2;
223 }
224 else
225 {
226 angle = GetAngularVel(rampTime) * rampTime / 2 +
227 GetAngularVel(time) * (time - rampTime);
228 }
229
230 // TODO: I want to integrate m_angularVelEqn up to current time
231
232 // Identity matrix
233 DNekMat rot(3, 3, 0.0);
234 rot(0, 0) = 1.0;
235 rot(1, 1) = 1.0;
236 rot(2, 2) = 1.0;
237
238 // Rodrigues' rotation formula in matrix form
239 rot = rot + sin(angle) * m_W + (1 - cos(angle)) * m_W2;
240
241 int cnt = 0;
242 for (auto &vert : m_verts)
243 {
245 DNekVec pntVec = {pnt[0], pnt[1], pnt[2]};
246
247 DNekVec newLoc = rot * pntVec;
248
249 vert->UpdatePosition(newLoc(0) + m_origin[0], newLoc(1) + m_origin[1],
250 newLoc(2) + m_origin[2]);
251 cnt++;
252 }
253
254 for (auto &curve : m_curves)
255 {
256 for (auto &vert : curve->m_points)
257 {
259 DNekVec pntVec = {pnt[0], pnt[1], pnt[2]};
260
261 DNekVec newLoc = rot * pntVec;
262
263 vert->UpdatePosition(newLoc(0) + m_origin[0],
264 newLoc(1) + m_origin[1],
265 newLoc(2) + m_origin[2]);
266 cnt++;
267 }
268 }
269
271
272 return true;
273}
274
275std::vector<NekDouble> ZoneTranslate::GetVel(NekDouble &time) const
276{
277 std::vector<NekDouble> vel(m_coordDim);
278 for (int i = 0; i < m_coordDim; ++i)
279 {
280 vel[i] = m_velocityEqns[i]->Evaluate(0, 0, 0, time);
281 }
282
283 return vel;
284}
285
286std::vector<NekDouble> ZoneTranslate::GetDisp(NekDouble &time)
287{
288 m_disp = std::vector<NekDouble>(m_coordDim);
289 for (int i = 0; i < m_coordDim; ++i)
290 {
291 m_disp[i] = m_displacementEqns[i]->Evaluate(0, 0, 0, time);
292 }
293 return m_disp;
294}
295
297{
298 auto disp = GetDisp(time);
299
300 int cnt = 0;
301 for (auto &vert : m_verts)
302 {
303 Array<OneD, NekDouble> newLoc(3, 0.0);
304 auto pnt = m_origVerts[cnt];
305
306 for (int i = 0; i < m_coordDim; ++i)
307 {
308 newLoc[i] = pnt(i) + disp[i];
309 }
310
311 vert->UpdatePosition(newLoc[0], newLoc[1], newLoc[2]);
312 cnt++;
313 }
314
315 for (auto &curve : m_curves)
316 {
317 for (auto &vert : curve->m_points)
318 {
319 Array<OneD, NekDouble> newLoc(3, 0.0);
320 auto pnt = m_origVerts[cnt];
321
322 for (int i = 0; i < m_coordDim; ++i)
323 {
324 newLoc[i] = pnt(i) + disp[i];
325 }
326
327 vert->UpdatePosition(newLoc[0], newLoc[1], newLoc[2]);
328 cnt++;
329 }
330 }
331
333
334 return true;
335}
336
337bool ZoneFixed::v_Move([[maybe_unused]] NekDouble time)
338{
339 return false;
340}
341
342std::vector<NekDouble> ZoneFixed::v_GetDisp() const
343{
344 std::vector<NekDouble> disp(m_coordDim);
345 for (int i = 0; i < m_coordDim; ++i)
346 {
347 disp[i] = 0.0;
348 }
349
350 return disp;
351}
352
353} // namespace Nektar::SpatialDomains
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
const NekPoint< NekDouble > origin
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:58
MovementType
Enum of zone movement type.
Definition: Zones.h:47
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:59
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:136
double NekDouble
Zone base: Contains the shared functions and variables.
Definition: Zones.h:61
ZoneBase(MovementType type, int indx, int domainID, CompositeMap domain, int coordDim)
Constructor.
std::array< std::set< GeometrySharedPtr >, 3 > m_constituentElements
Array of all dimension elements i.e. faces = [2], edges = [1], geom = [0].
Definition: Zones.h:151
std::vector< PointGeomSharedPtr > m_verts
Vector of all points in the zone.
Definition: Zones.h:157
int m_coordDim
Coordinate dimension.
Definition: Zones.h:155
void ClearBoundingBoxes()
Clears all bounding boxes associated with the zones elements.
std::vector< PointGeom > m_origVerts
Vector of all points in the zone at initialisation.
Definition: Zones.h:161
std::vector< CurveSharedPtr > m_curves
Vector of all curves in the zone.
Definition: Zones.h:159
CompositeMap m_domain
Zone domain.
Definition: Zones.h:146
std::vector< GeometrySharedPtr > m_elements
Vector of highest dimension zone elements.
Definition: Zones.h:148
std::vector< NekDouble > v_GetDisp() const override
Returns the displacement of the zone.
bool v_Move(NekDouble time) final
Virtual function for movement of the zone at.
DNekVec m_axis
Axis rotation is performed around.
Definition: Zones.h:220
NekDouble GetAngularVel(NekDouble &time) const
Return the angular velocity of the zone at.
DNekMat m_W
W matrix Rodrigues' rotation formula, cross product of axis.
Definition: Zones.h:224
NekPoint< NekDouble > m_origin
Origin point rotation is performed around.
Definition: Zones.h:218
LibUtilities::EquationSharedPtr m_angularVelEqn
Equation defining angular velocity as a function of time.
Definition: Zones.h:222
bool v_Move(NekDouble time) final
Virtual function for movement of the zone at.
DNekMat m_W2
W^2 matrix Rodrigues' rotation formula, cross product of axis squared.
Definition: Zones.h:226
ZoneRotate(int id, int domainID, const CompositeMap &domain, const int coordDim, const NekPoint< NekDouble > &origin, const DNekVec &axis, const LibUtilities::EquationSharedPtr &angularVelEqn)
std::vector< NekDouble > GetDisp(NekDouble &time)
Returns the displacement of the zone.
std::vector< NekDouble > m_disp
Definition: Zones.h:285
bool v_Move(NekDouble time) final
Virtual function for movement of the zone at.
Array< OneD, LibUtilities::EquationSharedPtr > m_displacementEqns
Definition: Zones.h:284
Array< OneD, LibUtilities::EquationSharedPtr > m_velocityEqns
Definition: Zones.h:283
std::vector< NekDouble > GetVel(NekDouble &time) const
Returns the velocity of the zone.