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
41namespace 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
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);
185
186 // Right hand side
187 nod_pt = finite_element_pt(2)->node_pt(n * np + np - 1);
190
191 // First part of lower boundary
192 nod_pt = finite_element_pt(3)->node_pt(n);
195
196 // First part of upper boundary
197 nod_pt = finite_element_pt(1)->node_pt(np * (np - 1) + n);
200
201 // First part of hole boundary
202 nod_pt = finite_element_pt(3)->node_pt(np * (np - 1) + n);
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);
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);
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);
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)
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
Quadrilateral mesh generator; Uses input from Geompack++. See: http://members.shaw....
Rectangular domain with circular whole.
Domain-based mesh for rectangular mesh with circular hole.
RectangleWithHoleDomain * Domain_pt
Pointer 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 ...
RectangleWithHoleDomain * domain_pt()
Access function to the domain.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
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 ...
////////////////////////////////////////////////////////////////////// //////////////////////////////...