one_d_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-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_ONE_D_MESH_TEMPLATE_CC
27#define OOMPH_ONE_D_MESH_TEMPLATE_CC
28
29
30// oomph-lib headers
31#include "one_d_mesh.template.h"
32
33
34namespace oomph
35{
36 /// The generic mesh construction routine --- this contains all the hard
37 /// work and is called by all constructors
38 template<class ELEMENT>
40 {
41 // Set the length of the domain
42 Length = Xmax - Xmin;
43
44 // Set the number of boundaries -- there are 2 boundaries in a 1D mesh
46
47 // Allocate storage for the pointers to the elements
48 Element_pt.resize(N);
49
50 // Allocate memory for the first element
51 Element_pt[0] = new ELEMENT;
52
53 // Read out the number of nodes in the element (the member function
54 // nnode_1d() is implemented in QElement)
55 const unsigned n_node =
56 dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
57
58 // We can now allocate storage for the pointers to the nodes in the mesh
59 Node_pt.resize(1 + (n_node - 1) * N);
60
61 // Initialise the node counter
62 unsigned node_count = 0;
63
64 // Initialise the minimum x coordinate in the mesh
65 const double xinit = Xmin;
66
67 // Calculate the length of the element
68 const double el_length = Length / double(N);
69
70 // Allocate storage for the local coordinate in the element
72
73 // If the number of elements is 1, the first element is also the
74 // last element
75 if (N == 1)
76 {
77 // Set the first node
78 // ------------------
79
80 // Allocate memory for the node, using the element's own construct_node
81 // function -- only the element knows what type of nodes it needs!
83 finite_element_pt(0)->construct_boundary_node(0, time_stepper_pt);
84
85 // Set the position of the node
86 node_pt(node_count)->x(0) = xinit;
87
88 // Add the node to the boundary 0
90
91 // Increment the counter for the nodes
92 node_count++;
93
94 // Now build central nodes (ignore last one which needs special
95 // ------------------------------------------------------------
96 // treatment because it's on the boundary)
97 // ---------------------------------------
98 for (unsigned jnod = 1; jnod < (n_node - 1); jnod++)
99 {
100 // Allocate memory for nodes, as before
102 finite_element_pt(0)->construct_node(jnod, time_stepper_pt);
103
104 // Get the local coordinate of the node
105 finite_element_pt(0)->local_fraction_of_node(jnod, s_fraction);
106
107 // Set the position of the node (linear mapping)
109
110 // Increment the node counter
111 node_count++;
112 }
113
114 // New final node
115 // --------------
116
117 // Allocate memory for the node, as before
118 Node_pt[node_count] = finite_element_pt(0)->construct_boundary_node(
120
121 // Set the position of the node
122 node_pt(node_count)->x(0) = xinit + Length;
123
124 // Add the node to the boundary 1
126
127 // Increment the node counter
128 node_count++;
129 }
130
131 // Otherwise, i.e. if there is more than one element, build all elements
132 else
133 {
134 // -------------
135 // FIRST ELEMENT
136 // -------------
137
138 // Set the first node
139 // ------------------
140
141 // Allocate memory for the node, using the element's own construct_node
142 // function -- only the element knows what type of nodes it needs!
144 finite_element_pt(0)->construct_boundary_node(0, time_stepper_pt);
145
146 // Set the position of the node
147 node_pt(node_count)->x(0) = xinit;
148
149 // Add the node to the boundary 0
151
152 // Increment the counter for the nodes
153 node_count++;
154
155 // Now build the other nodes in the first element
156 // ----------------------------------------------
157
158 // Loop over the other nodes in the first element
159 for (unsigned jnod = 1; jnod < n_node; jnod++)
160 {
161 // Allocate memory for the nodes
163 finite_element_pt(0)->construct_node(jnod, time_stepper_pt);
164
165 // Get the local coordinate of the node
166 finite_element_pt(0)->local_fraction_of_node(jnod, s_fraction);
167
168 // Set the position of the node (linear mapping)
170
171 // Increment the node counter
172 node_count++;
173 }
174
175 // ----------------
176 // CENTRAL ELEMENTS
177 // ----------------
178
179 // Loop over central elements in mesh
180 for (unsigned e = 1; e < (N - 1); e++)
181 {
182 // Allocate memory for the new element
183 Element_pt[e] = new ELEMENT;
184
185 // The first node is the same as the last node in the neighbouring
186 // element on the left
187 finite_element_pt(e)->node_pt(0) =
188 finite_element_pt(e - 1)->node_pt((n_node - 1));
189
190 // Loop over the remaining nodes in the element
191 for (unsigned jnod = 1; jnod < n_node; jnod++)
192 {
193 // Allocate memory for the nodes, as before
195 finite_element_pt(e)->construct_node(jnod, time_stepper_pt);
196
197 // Get the local coordinate of the nodes
198 finite_element_pt(e)->local_fraction_of_node(jnod, s_fraction);
199
200 // Set the position of the node (linear mapping)
201 node_pt(node_count)->x(0) = xinit + el_length * (e + s_fraction[0]);
202
203 // Increment the node counter
204 node_count++;
205 }
206 } // End of loop over central elements
207
208
209 // FINAL ELEMENT
210 //--------------
211
212 // Allocate memory for element
213 Element_pt[N - 1] = new ELEMENT;
214
215 // The first node is the same as the last node in the neighbouring
216 // element on the left
217 finite_element_pt(N - 1)->node_pt(0) =
218 finite_element_pt(N - 2)->node_pt(n_node - 1);
219
220 // New central nodes (ignore last one which needs special treatment
221 // because it's on the boundary)
222 for (unsigned jnod = 1; jnod < (n_node - 1); jnod++)
223 {
224 // Allocate memory for nodes, as before
226 finite_element_pt(N - 1)->construct_node(jnod, time_stepper_pt);
227
228 // Get the local coordinate of the node
229 finite_element_pt(N - 1)->local_fraction_of_node(jnod, s_fraction);
230
231 // Set the position of the node
232 node_pt(node_count)->x(0) = xinit + el_length * (N - 1 + s_fraction[0]);
233
234 // Increment the node counter
235 node_count++;
236 }
237
238 // New final node
239 // --------------
240
241 // Allocate memory for the node, as before
242 Node_pt[node_count] = finite_element_pt(N - 1)->construct_boundary_node(
244
245 // Set the position of the node
246 node_pt(node_count)->x(0) = xinit + Length;
247
248 // Add the node to the boundary 1
250
251 // Increment the node counter
252 node_count++;
253 }
254 }
255
256} // namespace oomph
257
258#endif
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh constuction routine, called by all constructors.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
////////////////////////////////////////////////////////////////////// //////////////////////////////...