rectangle_with_hole_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 #ifndef OOMPH_RECTANGLE_WITH_HOLE_MESH_HEADER
27 #define OOMPH_RECTANGLE_WITH_HOLE_MESH_HEADER
28 
29 // Config header generated by autoconfig
30 #ifdef HAVE_CONFIG_H
31 #include <oomph-lib-config.h>
32 #endif
33 
34 // OOMPH-LIB headers
35 #include "../generic/mesh.h"
36 #include "../generic/quadtree.h"
37 #include "../generic/quad_mesh.h"
38 #include "../generic/refineable_quad_mesh.h"
40 
41 namespace oomph
42 {
43  //=============================================================
44  /// Domain-based mesh for rectangular mesh with circular hole
45  //=============================================================
46  template<class ELEMENT>
47  class RectangleWithHoleMesh : public virtual Mesh
48  {
49  public:
50  /// Constructor: Pass pointer to geometric object that
51  /// represents the cylinder, the length and height of the domain.
52  /// The GeomObject must be parametrised such that
53  /// \f$\zeta \in [0,2\pi]\f$ sweeps around the circumference
54  /// in anticlockwise direction. Timestepper defaults to Steady
55  /// default timestepper.
57  GeomObject* cylinder_pt,
58  const double& length,
59  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
60  {
61  // Create the domain
62  Domain_pt = new RectangleWithHoleDomain(cylinder_pt, length);
63 
64  // Initialise the node counter
65  unsigned long node_count = 0;
66 
67  // Vectors used to get data from domains
68  Vector<double> s(2), r(2);
69 
70  // Setup temporary storage for the Node
71  Vector<Node*> Tmp_node_pt;
72 
73  // Now blindly loop over the macro elements and associate and finite
74  // element with each
75  unsigned nmacro_element = Domain_pt->nmacro_element();
76  for (unsigned e = 0; e < nmacro_element; e++)
77  {
78  // Create the FiniteElement and add to the Element_pt Vector
79  Element_pt.push_back(new ELEMENT);
80 
81  // Read out the number of linear points in the element
82  unsigned np = dynamic_cast<ELEMENT*>(finite_element_pt(e))->nnode_1d();
83 
84  // Loop over nodes in the column
85  for (unsigned l1 = 0; l1 < np; l1++)
86  {
87  // Loop over the nodes in the row
88  for (unsigned l2 = 0; l2 < np; l2++)
89  {
90  // Allocate the memory for the node
91  Tmp_node_pt.push_back(finite_element_pt(e)->construct_node(
92  l1 * np + l2, time_stepper_pt));
93 
94  // Read out the position of the node from the macro element
95  s[0] = -1.0 + 2.0 * (double)l2 / (double)(np - 1);
96  s[1] = -1.0 + 2.0 * (double)l1 / (double)(np - 1);
97  Domain_pt->macro_element_pt(e)->macro_map(s, r);
98 
99  // Set the position of the node
100  Tmp_node_pt[node_count]->x(0) = r[0];
101  Tmp_node_pt[node_count]->x(1) = r[1];
102 
103  // Increment the node number
104  node_count++;
105  }
106  }
107  } // End of loop over macro elements
108 
109 
110  // Now the elements have been created, but there will be nodes in
111  // common, need to loop over the common edges and sort it, by reassigning
112  // pointers and the deleting excess nodes
113 
114  // Read out the number of linear points in the element
115  unsigned np = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
116 
117  // Edge between Elements 0 and 1
118  for (unsigned n = 0; n < np; n++)
119  {
120  // Set the nodes in element 1 to be the same as in element 0
121  finite_element_pt(1)->node_pt(n * np) =
122  finite_element_pt(0)->node_pt((np - 1) * np + np - 1 - n);
123 
124  // Remove the nodes in element 1 from the temporary node list
125  delete Tmp_node_pt[np * np + n * np];
126  Tmp_node_pt[np * np + n * np] = 0;
127  }
128 
129  // Edge between Elements 0 and 3
130  for (unsigned n = 0; n < np; n++)
131  {
132  // Set the nodes in element 3 to be the same as in element 0
133  finite_element_pt(3)->node_pt(n * np) =
134  finite_element_pt(0)->node_pt(n);
135 
136  // Remove the nodes in element 3 from the temporary node list
137  delete Tmp_node_pt[3 * np * np + n * np];
138  Tmp_node_pt[3 * np * np + n * np] = 0;
139  }
140 
141  // Edge between Element 1 and 2
142  for (unsigned n = 0; n < np; n++)
143  {
144  // Set the nodes in element 2 to be the same as in element 1
145  finite_element_pt(2)->node_pt(np * (np - 1) + n) =
146  finite_element_pt(1)->node_pt(np * n + np - 1);
147 
148  // Remove the nodes in element 2 from the temporary node list
149  delete Tmp_node_pt[2 * np * np + np * (np - 1) + n];
150  Tmp_node_pt[2 * np * np + np * (np - 1) + n] = 0;
151  }
152 
153 
154  // Edge between Element 3 and 2
155  for (unsigned n = 0; n < np; n++)
156  {
157  // Set the nodes in element 2 to be the same as in element 3
158  finite_element_pt(2)->node_pt(n) =
159  finite_element_pt(3)->node_pt(np * (np - n - 1) + np - 1);
160 
161  // Remove the nodes in element 2 from the temporary node list
162  delete Tmp_node_pt[2 * np * np + n];
163  Tmp_node_pt[2 * np * np + n] = 0;
164  }
165 
166 
167  // Now set the actual true nodes
168  for (unsigned long n = 0; n < node_count; n++)
169  {
170  if (Tmp_node_pt[n] != 0)
171  {
172  Node_pt.push_back(Tmp_node_pt[n]);
173  }
174  }
175 
176  // Finally set the nodes on the boundaries
177  set_nboundary(5);
178 
179  for (unsigned n = 0; n < np; n++)
180  {
181  // Left hand side
182  Node* nod_pt = finite_element_pt(0)->node_pt(n * np);
183  convert_to_boundary_node(nod_pt);
184  add_boundary_node(3, nod_pt);
185 
186  // Right hand side
187  nod_pt = finite_element_pt(2)->node_pt(n * np + np - 1);
188  convert_to_boundary_node(nod_pt);
189  add_boundary_node(1, nod_pt);
190 
191  // First part of lower boundary
192  nod_pt = finite_element_pt(3)->node_pt(n);
193  convert_to_boundary_node(nod_pt);
194  add_boundary_node(0, nod_pt);
195 
196  // First part of upper boundary
197  nod_pt = finite_element_pt(1)->node_pt(np * (np - 1) + n);
198  convert_to_boundary_node(nod_pt);
199  add_boundary_node(2, nod_pt);
200 
201  // First part of hole boundary
202  nod_pt = finite_element_pt(3)->node_pt(np * (np - 1) + n);
203  convert_to_boundary_node(nod_pt);
204  add_boundary_node(4, nod_pt);
205  }
206 
207  for (unsigned n = 1; n < np; n++)
208  {
209  // Next part of hole
210  Node* nod_pt = finite_element_pt(2)->node_pt(n * np);
211  convert_to_boundary_node(nod_pt);
212  add_boundary_node(4, nod_pt);
213  }
214 
215  for (unsigned n = 1; n < np; n++)
216  {
217  // Next part of hole
218  Node* nod_pt = finite_element_pt(1)->node_pt(np - n - 1);
219  convert_to_boundary_node(nod_pt);
220  add_boundary_node(4, nod_pt);
221  }
222 
223  for (unsigned n = 1; n < np - 1; n++)
224  {
225  // Final part of hole
226  Node* nod_pt =
227  finite_element_pt(0)->node_pt(np * (np - n - 1) + np - 1);
228  convert_to_boundary_node(nod_pt);
229  add_boundary_node(4, nod_pt);
230  }
231  }
232 
233  /// Access function to the domain
235  {
236  return Domain_pt;
237  }
238 
239  protected:
240  /// Pointer to the domain
242  };
243 
244 
245  /// ///////////////////////////////////////////////////////////////////////
246  /// ///////////////////////////////////////////////////////////////////////
247  /// ///////////////////////////////////////////////////////////////////////
248 
249 
250  //===================================================================
251  /// Refineable version of RectangleWithHoleMesh. For some reason
252  /// this needs on uniform refinement to work...
253  //===================================================================
254  template<class ELEMENT>
256  public RefineableQuadMesh<ELEMENT>
257  {
258  public:
259  /// Constructor. Pass pointer to geometric object that
260  /// represents the cylinder, the length and height of the domain.
261  /// The GeomObject must be parametrised such that
262  /// \f$\zeta \in [0,2\pi]\f$ sweeps around the circumference
263  /// in anticlockwise direction. Timestepper defaults to Steady
264  /// default timestepper.
266  GeomObject* cylinder_pt,
267  const double& length,
268  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
269  : RectangleWithHoleMesh<ELEMENT>(cylinder_pt, length, time_stepper_pt)
270  {
271  // Nodal positions etc. were created in constructor for
272  // Cylinder...<...>. Need to setup adaptive information.
273 
274  // Loop over all elements and set macro element pointer
275  for (unsigned e = 0; e < 4; e++)
276  {
277  dynamic_cast<ELEMENT*>(this->element_pt(e))
278  ->set_macro_elem_pt(this->Domain_pt->macro_element_pt(e));
279  }
280 
281  // Setup boundary element lookup schemes
282  this->setup_boundary_element_info();
283 
284  // Nodal positions etc. were created in constructor for
285  // RectangularMesh<...>. Only need to setup quadtree forest
286  this->setup_quadtree_forest();
287  }
288 
289  /// Destructor: Empty
291  };
292 
293 
294 } // namespace oomph
295 
296 
297 #endif
Rectangular domain with circular whole.
Domain-based mesh for rectangular mesh with circular hole.
RectangleWithHoleDomain * Domain_pt
Pointer to the domain.
RectangleWithHoleDomain * domain_pt()
Access function to the domain.
RectangleWithHoleMesh(GeomObject *cylinder_pt, const double &length, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that represents the cylinder, the length and height of ...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
RefineableRectangleWithHoleMesh(GeomObject *cylinder_pt, const double &length, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. Pass pointer to geometric object that represents the cylinder, the length and height of ...
////////////////////////////////////////////////////////////////////// //////////////////////////////...