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-2023 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 
38 namespace 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
Annular domain.
RectangularQuadMesh is a two-dimensional mesh of Quad elements with Nx elements in the "x" (horizonal...
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
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.
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.
////////////////////////////////////////////////////////////////////// //////////////////////////////...