annular_mesh.template.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26
27#ifndef OOMPH_ANNULAR_MESH_HEADER
28#define OOMPH_ANNULAR_MESH_HEADER
29
30// Headers
33
34
35// Include the headers file for domain
36#include "annular_domain.h"
37
38namespace oomph
39{
40 //===================================================================
41 /// 2D annular mesh with a unit circle in the middle and a layer
42 /// of thickness h surrounding it.
43 //==================================================================
44 template<class ELEMENT>
45 class TwoDAnnularMesh : public virtual RectangularQuadMesh<ELEMENT>
46 {
47 public:
48 /// Constructor
49 TwoDAnnularMesh(const bool& periodic,
50 const double& azimuthal_fraction,
51 const unsigned& ntheta,
52 const unsigned& nr,
53 const double& a,
54 const double& h,
55 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
56 : RectangularQuadMesh<ELEMENT>(
57 ntheta, nr, 1.0, 1.0, periodic, time_stepper_pt)
58 {
59 // Mesh can only be built with 2D Qelements.
60 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
61
62 // Wrap mesh into annular shape
63 double phi = 0.0;
64 wrap_into_annular_shape(a, h, azimuthal_fraction, phi);
65 }
66
67 /// Constructor; rotate mesh by angle phi.
68 TwoDAnnularMesh(const bool& periodic,
69 const double& azimuthal_fraction,
70 const unsigned& ntheta,
71 const unsigned& nr,
72 const double& a,
73 const double& h,
74 const double& phi,
75 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
76 : RectangularQuadMesh<ELEMENT>(
77 ntheta, nr, 1.0, 1.0, periodic, time_stepper_pt)
78 {
79 // Mesh can only be built with 2D Qelements.
80 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
81
82 // Wrap mesh into annular shape
83 wrap_into_annular_shape(a, h, azimuthal_fraction, phi);
84 }
85
86
87 private:
88 /// Wrap mesh into annular shape
89 void wrap_into_annular_shape(const double& a,
90 const double& h,
91 const double& azimuthal_fraction,
92 const double& phi);
93 };
94
95
96 /// ///////////////////////////////////////////////////////////////////
97 /// ///////////////////////////////////////////////////////////////////
98 /// ///////////////////////////////////////////////////////////////////
99
100
101 //===================================================================
102 /// Refineable 2D annular mesh with a unit circle in the middle and a layer
103 /// of thickness h surrounding it.
104 //==================================================================
105 template<class ELEMENT>
106 class RefineableTwoDAnnularMesh : public virtual TwoDAnnularMesh<ELEMENT>,
107 public virtual RefineableQuadMesh<ELEMENT>
108 {
109 public:
110 /// Constructor
112 const bool& periodic,
113 const double& azimuthal_fraction,
114 const unsigned& ntheta,
115 const unsigned& nr,
116 const double& a,
117 const double& h,
118 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
119 : RectangularQuadMesh<ELEMENT>(
120 ntheta, nr, 1.0, 1.0, periodic, time_stepper_pt),
121 TwoDAnnularMesh<ELEMENT>(
122 periodic, azimuthal_fraction, ntheta, nr, a, h, time_stepper_pt)
123 {
124 // Mesh can only be built with 2D Qelements.
125 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
126
127 // Build macro element-based domain
128 double phi = 0.0;
129 Domain_pt = new AnnularDomain(azimuthal_fraction, ntheta, nr, a, h, phi);
130
131 // Loop over all elements and set macro element pointer
132 unsigned nel = this->nelement();
133 for (unsigned ielem = 0; ielem < nel; ielem++)
134 {
135 dynamic_cast<RefineableQElement<2>*>(this->element_pt(ielem))
136 ->set_macro_elem_pt(this->Domain_pt->macro_element_pt(ielem));
137 }
138
139 // Update nodal positions based on macro-element representation
140 this->node_update();
141
142 // Nodal positions etc. were created in constructor for
143 // RectangularMesh<...>. Only need to setup quadtree forest
144 this->setup_quadtree_forest();
145
146 // Setup the periodic neighbour information of the TreeRoots
147 // Cast to specific elements
148 if (periodic)
149 {
150 Vector<TreeRoot*> left_root_pt(nr);
151 Vector<TreeRoot*> right_root_pt(nr);
152 for (unsigned i = 0; i < nr; i++)
153 {
154 left_root_pt[i] = dynamic_cast<ELEMENT*>(this->element_pt(i * ntheta))
155 ->tree_pt()
156 ->root_pt();
157
158 right_root_pt[i] =
159 dynamic_cast<ELEMENT*>(this->element_pt((i + 1) * ntheta - 1))
160 ->tree_pt()
161 ->root_pt();
162 }
163
164 // Set the neighbour and periodicity
165 using namespace QuadTreeNames;
166 for (unsigned i = 0; i < nr; i++)
167 {
168 left_root_pt[i]->neighbour_pt(W) = right_root_pt[i];
169 left_root_pt[i]->set_neighbour_periodic(W);
170
171 right_root_pt[i]->neighbour_pt(E) = left_root_pt[i];
172 right_root_pt[i]->set_neighbour_periodic(E);
173 }
174 }
175 }
176
177
178 /// Constructor; rotate mesh by angle phi
180 const bool& periodic,
181 const double& azimuthal_fraction,
182 const unsigned& ntheta,
183 const unsigned& nr,
184 const double& a,
185 const double& h,
186 const double& phi,
187 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
188 : RectangularQuadMesh<ELEMENT>(
189 ntheta, nr, 1.0, 1.0, periodic, time_stepper_pt),
190 TwoDAnnularMesh<ELEMENT>(
191 periodic, azimuthal_fraction, ntheta, nr, a, h, phi, time_stepper_pt)
192 {
193 // Mesh can only be built with 2D Qelements.
194 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
195
196 // Build macro element-based domain
197 Domain_pt = new AnnularDomain(azimuthal_fraction, ntheta, nr, a, h, phi);
198
199 // Loop over all elements and set macro element pointer
200 unsigned nel = this->nelement();
201 for (unsigned ielem = 0; ielem < nel; ielem++)
202 {
203 dynamic_cast<RefineableQElement<2>*>(this->element_pt(ielem))
204 ->set_macro_elem_pt(this->Domain_pt->macro_element_pt(ielem));
205 }
206
207 // Update nodal positions based on macro-element representation
208 this->node_update();
209
210 // Nodal positions etc. were created in constructor for
211 // RectangularMesh<...>. Only need to setup quadtree forest
212 this->setup_quadtree_forest();
213
214 // Setup the periodic neighbour information of the TreeRoots
215 // Cast to specific elements
216 if (periodic)
217 {
218 Vector<TreeRoot*> left_root_pt(nr);
219 Vector<TreeRoot*> right_root_pt(nr);
220 for (unsigned i = 0; i < nr; i++)
221 {
222 left_root_pt[i] = dynamic_cast<ELEMENT*>(this->element_pt(i * ntheta))
223 ->tree_pt()
224 ->root_pt();
225
226 right_root_pt[i] =
227 dynamic_cast<ELEMENT*>(this->element_pt((i + 1) * ntheta - 1))
228 ->tree_pt()
229 ->root_pt();
230 }
231
232 // Set the neighbour and periodicity
233 using namespace QuadTreeNames;
234 for (unsigned i = 0; i < nr; i++)
235 {
236 left_root_pt[i]->neighbour_pt(W) = right_root_pt[i];
237 left_root_pt[i]->set_neighbour_periodic(W);
238
239 right_root_pt[i]->neighbour_pt(E) = left_root_pt[i];
240 right_root_pt[i]->set_neighbour_periodic(E);
241 }
242 }
243 }
244
245 private:
246 /// Pointer to domain
248 };
249
250} // namespace oomph
251
252#endif
cstr elem_len * i
Definition: cfortran.h:603
Annular domain.
MacroElement * macro_element_pt(const unsigned &i)
Access to i-th macro element.
Definition: domain.h:116
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
virtual void node_update(const bool &update_all_solid_nodes=false)
Update nodal positions in response to changes in the domain shape. Uses the FiniteElement::get_x(....
Definition: mesh.cc:287
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
RectangularQuadMesh is a two-dimensional mesh of Quad elements with Nx elements in the "x" (horizonal...
Refineable version of QElement<2,NNODE_1D>.
Intermediate mesh class that implements the mesh adaptation functions specified in the TreeBasedRefin...
void setup_quadtree_forest()
Set up QuadTreeForest. Wipes any existing tree structure below the minimum refinement level and regar...
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
RefineableTwoDAnnularMesh(const bool &periodic, const double &azimuthal_fraction, const unsigned &ntheta, const unsigned &nr, const double &a, const double &h, const double &phi, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor; rotate mesh by angle phi.
AnnularDomain * Domain_pt
Pointer to domain.
RefineableTwoDAnnularMesh(const bool &periodic, const double &azimuthal_fraction, const unsigned &ntheta, const unsigned &nr, const double &a, const double &h, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
2D annular mesh with a unit circle in the middle and a layer of thickness h surrounding it.
TwoDAnnularMesh(const bool &periodic, const double &azimuthal_fraction, const unsigned &ntheta, const unsigned &nr, const double &a, const double &h, const double &phi, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor; rotate mesh by angle phi.
TwoDAnnularMesh(const bool &periodic, const double &azimuthal_fraction, const unsigned &ntheta, const unsigned &nr, const double &a, const double &h, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor.
void wrap_into_annular_shape(const double &a, const double &h, const double &azimuthal_fraction, const double &phi)
Wrap mesh into annular shape.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...