single_layer_cubic_spine_mesh.template.cc
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#ifndef OOMPH_SINGLE_LAYER_CUBIC_SPINE_MESH_TEMPLATE_CC
27#define OOMPH_SINGLE_LAYER_CUBIC_SPINE_MESH_TEMPLATE_CC
28
31
32namespace oomph
33{
34 //===========================================================================
35 /// Constructor for spine 3D mesh: Pass number of elements in x-direction,
36 /// number of elements in y-direction, number elements in z-direction,
37 /// length, width and height of layer,
38 /// and pointer to timestepper (defaults to Static timestepper).
39 //===========================================================================
40 template<class ELEMENT>
42 const unsigned& nx,
43 const unsigned& ny,
44 const unsigned& nz,
45 const double& lx,
46 const double& ly,
47 const double& h,
48 TimeStepper* time_stepper_pt)
49 : SimpleCubicMesh<ELEMENT>(nx, ny, nz, lx, ly, h, time_stepper_pt)
50 {
51 // Mesh can only be built with 3D Qelements.
52 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
53
54 // Mesh can only be built with spine elements
55 MeshChecker::assert_geometric_element<SpineFiniteElement, ELEMENT>(3);
56
57 // Now build the single layer mesh on top of the existing mesh
58 build_single_layer_mesh(time_stepper_pt);
59 }
60
61 //===========================================================================
62 /// Helper function that actually builds the single-layer spine mesh
63 /// based on the parameters set in the various constructors
64 //===========================================================================
65 template<class ELEMENT>
67 TimeStepper* time_stepper_pt)
68 {
69 // Read out the number of elements in the x-direction
70 unsigned n_x = this->Nx;
71 // Read out the number of elements in the y-direction
72 unsigned n_y = this->Ny;
73 // Read out the number of elements in the z-direction
74 unsigned n_z = this->Nz;
75
76 // Allocate memory for the spines and fractions along spines
77
78 // Read out number of linear points in the element
79 unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
80
81 // Allocate store for the spines: (different in the case of periodic meshes
82 // !!)
83 Spine_pt.reserve(((n_p - 1) * n_x + 1) * ((n_p - 1) * n_y + 1));
84
85 // Now we loop over all the elements and attach the spines
86
87 // FIRST ELEMENT: Element 0
88 // Loop over the nodes on the base of the element
89 for (unsigned l1 = 0; l1 < n_p; l1++) // y loop over the nodes
90 {
91 for (unsigned l2 = 0; l2 < n_p; l2++) // x loop over the nodes
92 {
93 // Assign the new spine with length h
94 Spine* new_spine_pt = new Spine(1.0);
95 Spine_pt.push_back(new_spine_pt);
96
97 // Get pointer to node
98 SpineNode* nod_pt = element_node_pt(0, l2 + l1 * n_p);
99 // Set the pointer to the spine
100 nod_pt->spine_pt() = new_spine_pt;
101 // Set the fraction
102 nod_pt->fraction() = 0.0;
103 // Pointer to the mesh that implements the update fct
104 nod_pt->spine_mesh_pt() = this;
105
106 // Loop vertically along the spine
107 // Loop over the elements
108 for (unsigned long k = 0; k < n_z; k++)
109 {
110 // Loop over the vertical nodes, apart from the first
111 for (unsigned l3 = 1; l3 < n_p; l3++)
112 {
113 // Get pointer to node
114 SpineNode* nod_pt =
115 element_node_pt(k * n_x * n_y, l3 * n_p * n_p + l2 + l1 * n_p);
116 // Set the pointer to the spine
117 nod_pt->spine_pt() = new_spine_pt;
118 // Set the fraction
119 nod_pt->fraction() =
120 (double(k) + double(l3) / double(n_p - 1)) / double(n_z);
121 // Pointer to the mesh that implements the update fct
122 nod_pt->spine_mesh_pt() = this;
123 }
124 }
125 }
126 }
127
128
129 // LOOP OVER OTHER ELEMENTS IN THE FIRST ROW
130 //-----------------------------------------
131
132 // The procedure is the same but we have to identify the
133 // before defined spines for not defining them two times
134
135 for (unsigned j = 1; j < n_x; j++) // loop over the elements in the first
136 // row
137 {
138 for (unsigned l1 = 0; l1 < n_p; l1++) // y loop over the nodes
139 {
140 // First we copy the last row of nodes into the
141 // first row of the new element (and extend to the third dimension)
142 for (unsigned l2 = 1; l2 < n_p; l2++) // x loop over the nodes
143 {
144 // Node j + i*np
145 // Assign the new spine with unit length
146 Spine* new_spine_pt = new Spine(1.0);
147 Spine_pt.push_back(new_spine_pt);
148
149 // Get pointer to node
150 SpineNode* nod_pt = element_node_pt(j, l2 + l1 * n_p);
151
152 // Set the pointer to the spine
153 nod_pt->spine_pt() = new_spine_pt;
154 // Set the fraction
155 nod_pt->fraction() = 0.0;
156 // Pointer to the mesh that implements the update fct
157 nod_pt->spine_mesh_pt() = this;
158
159 // Loop vertically along the spine
160 // Loop over the elements
161 for (unsigned long k = 0; k < n_z; k++)
162 {
163 // Loop over the vertical nodes, apart from the first
164 for (unsigned l3 = 1; l3 < n_p; l3++)
165 {
166 // Get pointer to node
167 SpineNode* nod_pt = element_node_pt(
168 j + k * n_x * n_y, l3 * n_p * n_p + l2 + l1 * n_p);
169 // Set the pointer to the spine
170 nod_pt->spine_pt() = new_spine_pt;
171 // Set the fraction
172 nod_pt->fraction() =
173 (double(k) + double(l3) / double(n_p - 1)) / double(n_z);
174 // Pointer to the mesh that implements the update fct
175 nod_pt->spine_mesh_pt() = this;
176 }
177 }
178 }
179 }
180 }
181
182 // REST OF THE ELEMENTS
183 // Now we loop over the rest of the elements.
184 // We will separate the first of each row being al the rest equal
185 for (unsigned long i = 1; i < n_y; i++)
186 {
187 // FIRST ELEMENT OF THE ROW
188
189 // First line of nodes is copied from the element of the bottom
190 for (unsigned l1 = 1; l1 < n_p; l1++) // y loop over the nodes
191 {
192 for (unsigned l2 = 0; l2 < n_p; l2++) // x loop over the nodes
193 {
194 // Node j + i*np
195 // Assign the new spine with unit length
196 Spine* new_spine_pt = new Spine(1.0);
197 Spine_pt.push_back(new_spine_pt);
198
199
200 // Get pointer to node
201 // Element i*n_x; node l2 + l1*n_p
202 SpineNode* nod_pt = element_node_pt(i * n_x, l2 + l1 * n_p);
203 // Set the pointer to the spine
204 nod_pt->spine_pt() = new_spine_pt;
205 // Set the fraction
206 nod_pt->fraction() = 0.0;
207 // Pointer to the mesh that implements the update fct
208 nod_pt->spine_mesh_pt() = this;
209
210 // Loop vertically along the spine
211 // Loop over the elements
212 for (unsigned long k = 0; k < n_z; k++)
213 {
214 // Loop over the vertical nodes, apart from the first
215 for (unsigned l3 = 1; l3 < n_p; l3++)
216 {
217 // Get pointer to node
218 SpineNode* nod_pt = element_node_pt(
219 i * n_x + k * n_x * n_y, l3 * n_p * n_p + l2 + l1 * n_p);
220 // Set the pointer to the spine
221 nod_pt->spine_pt() = new_spine_pt;
222 // Set the fraction
223 nod_pt->fraction() =
224 (double(k) + double(l3) / double(n_p - 1)) / double(n_z);
225 // Pointer to the mesh that implements the update fct
226 nod_pt->spine_mesh_pt() = this;
227 }
228 }
229 }
230 }
231
232
233 // REST OF THE ELEMENTS OF THE ROW
234 for (unsigned j = 1; j < n_x; j++)
235 {
236 // First line of nodes is copied from the element of the bottom
237 for (unsigned l1 = 1; l1 < n_p; l1++) // y loop over the nodes
238 {
239 for (unsigned l2 = 1; l2 < n_p; l2++) // x loop over the nodes
240 {
241 // Node j + i*np
242 // Assign the new spine with unit length
243 Spine* new_spine_pt = new Spine(1.0);
244 Spine_pt.push_back(new_spine_pt);
245
246
247 // Get pointer to node
248 // Element j + i*n_x; node l2 + l1*n_p
249 SpineNode* nod_pt = element_node_pt(j + i * n_x, l2 + l1 * n_p);
250 // Set the pointer to the spine
251 nod_pt->spine_pt() = new_spine_pt;
252 // Set the fraction
253 nod_pt->fraction() = 0.0;
254 // Pointer to the mesh that implements the update fct
255 nod_pt->spine_mesh_pt() = this;
256
257 // Loop vertically along the spine
258 // Loop over the elements
259 for (unsigned long k = 0; k < n_z; k++)
260 {
261 // Loop over the vertical nodes, apart from the first
262 for (unsigned l3 = 1; l3 < n_p; l3++)
263 {
264 // Get pointer to node
265 SpineNode* nod_pt = element_node_pt(
266 j + i * n_x + k * n_x * n_y, l3 * n_p * n_p + l2 + l1 * n_p);
267 // Set the pointer to the spine
268 nod_pt->spine_pt() = new_spine_pt;
269 // Set the fraction
270 nod_pt->fraction() =
271 (double(k) + double(l3) / double(n_p - 1)) / double(n_z);
272 // Pointer to the mesh that implements the update fct
273 nod_pt->spine_mesh_pt() = this;
274 }
275 }
276 }
277 }
278 }
279 }
280 }
281
282} // namespace oomph
283#endif
cstr elem_len * i
Definition: cfortran.h:603
Simple cubic 3D Brick mesh class.
virtual void build_single_layer_mesh(TimeStepper *time_stepper_pt)
Helper function to actually build the single-layer spine mesh (called from various constructors)
SingleLayerCubicSpineMesh(const unsigned &nx, const unsigned &ny, const unsigned &nz, const double &lx, const double &ly, const double &h, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction,...
Class for nodes that live on spines. The assumption is that each Node lies at a fixed fraction on a s...
Definition: spines.h:328
double & fraction()
Set reference to fraction along spine.
Definition: spines.h:378
SpineMesh *& spine_mesh_pt()
Access function to Pointer to SpineMesh that this node is a part of and which implements the node upd...
Definition: spines.h:391
Spine *& spine_pt()
Access function to spine.
Definition: spines.h:372
Spines are used for algebraic node update operations in free-surface fluid problems: They form the ba...
Definition: spines.h:64
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
//////////////////////////////////////////////////////////////////// ////////////////////////////////...