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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_ONE_D_MESH_TEMPLATE_CC
31 #define OOMPH_ONE_D_MESH_TEMPLATE_CC
32 
33 
34 
35 // oomph-lib headers
36 #include "one_d_mesh.template.h"
37 
38 
39 
40 namespace oomph
41 {
42 
43 
44 /// The generic mesh construction routine --- this contains all the hard
45 /// work and is called by all constructors
46  template <class ELEMENT>
48  {
49  // Set the length of the domain
50  Length = Xmax - Xmin;
51 
52  // Set the number of boundaries -- there are 2 boundaries in a 1D mesh
53  set_nboundary(2);
54 
55  // Allocate storage for the pointers to the elements
56  Element_pt.resize(N);
57 
58  // Allocate memory for the first element
59  Element_pt[0] = new ELEMENT;
60 
61  // Read out the number of nodes in the element (the member function
62  // nnode_1d() is implemented in QElement)
63  const unsigned n_node
64  = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
65 
66  // We can now allocate storage for the pointers to the nodes in the mesh
67  Node_pt.resize(1+(n_node-1)*N);
68 
69  // Initialise the node counter
70  unsigned node_count = 0;
71 
72  // Initialise the minimum x coordinate in the mesh
73  const double xinit = Xmin;
74 
75  // Calculate the length of the element
76  const double el_length = Length/double(N);
77 
78  // Allocate storage for the local coordinate in the element
79  Vector<double> s_fraction;
80 
81  // If the number of elements is 1, the first element is also the
82  // last element
83  if(N==1)
84  {
85  // Set the first node
86  // ------------------
87 
88  // Allocate memory for the node, using the element's own construct_node
89  // function -- only the element knows what type of nodes it needs!
90  Node_pt[node_count]
91  = finite_element_pt(0)->construct_boundary_node(0,time_stepper_pt);
92 
93  // Set the position of the node
94  node_pt(node_count)->x(0) = xinit;
95 
96  // Add the node to the boundary 0
97  add_boundary_node(0,Node_pt[node_count]);
98 
99  // Increment the counter for the nodes
100  node_count++;
101 
102  // Now build central nodes (ignore last one which needs special
103  // ------------------------------------------------------------
104  // treatment because it's on the boundary)
105  // ---------------------------------------
106  for(unsigned jnod=1;jnod<(n_node-1);jnod++)
107  {
108  // Allocate memory for nodes, as before
109  Node_pt[node_count]
110  = finite_element_pt(0)->construct_node(jnod,time_stepper_pt);
111 
112  // Get the local coordinate of the node
113  finite_element_pt(0)->local_fraction_of_node(jnod,s_fraction);
114 
115  //Set the position of the node (linear mapping)
116  node_pt(node_count)->x(0) = xinit + el_length*s_fraction[0];
117 
118  //Increment the node counter
119  node_count++;
120  }
121 
122  // New final node
123  // --------------
124 
125  // Allocate memory for the node, as before
126  Node_pt[node_count] = finite_element_pt(0)
127  ->construct_boundary_node(n_node-1,time_stepper_pt);
128 
129  // Set the position of the node
130  node_pt(node_count)->x(0) = xinit + Length;
131 
132  // Add the node to the boundary 1
133  add_boundary_node(1,Node_pt[node_count]);
134 
135  // Increment the node counter
136  node_count++;
137  }
138 
139  // Otherwise, i.e. if there is more than one element, build all elements
140  else
141  {
142  // -------------
143  // FIRST ELEMENT
144  // -------------
145 
146  // Set the first node
147  // ------------------
148 
149  // Allocate memory for the node, using the element's own construct_node
150  // function -- only the element knows what type of nodes it needs!
151  Node_pt[node_count]
152  = finite_element_pt(0)->construct_boundary_node(0,time_stepper_pt);
153 
154  // Set the position of the node
155  node_pt(node_count)->x(0) = xinit;
156 
157  // Add the node to the boundary 0
158  add_boundary_node(0,Node_pt[node_count]);
159 
160  //Increment the counter for the nodes
161  node_count++;
162 
163  // Now build the other nodes in the first element
164  // ----------------------------------------------
165 
166  // Loop over the other nodes in the first element
167  for(unsigned jnod=1;jnod<n_node;jnod++)
168  {
169  // Allocate memory for the nodes
170  Node_pt[node_count]
171  = finite_element_pt(0)->construct_node(jnod,time_stepper_pt);
172 
173  // Get the local coordinate of the node
174  finite_element_pt(0)->local_fraction_of_node(jnod,s_fraction);
175 
176  // Set the position of the node (linear mapping)
177  node_pt(node_count)->x(0) = xinit + el_length*s_fraction[0];
178 
179  // Increment the node counter
180  node_count++;
181  }
182 
183  // ----------------
184  // CENTRAL ELEMENTS
185  // ----------------
186 
187  // Loop over central elements in mesh
188  for(unsigned e=1;e<(N-1);e++)
189  {
190  // Allocate memory for the new element
191  Element_pt[e] = new ELEMENT;
192 
193  // The first node is the same as the last node in the neighbouring
194  // element on the left
195  finite_element_pt(e)->node_pt(0)
196  = finite_element_pt(e-1)->node_pt((n_node-1));
197 
198  // Loop over the remaining nodes in the element
199  for(unsigned jnod=1;jnod<n_node;jnod++)
200  {
201  // Allocate memory for the nodes, as before
202  Node_pt[node_count]
203  = finite_element_pt(e)->construct_node(jnod,time_stepper_pt);
204 
205  // Get the local coordinate of the nodes
206  finite_element_pt(e)->local_fraction_of_node(jnod,s_fraction);
207 
208  // Set the position of the node (linear mapping)
209  node_pt(node_count)->x(0) = xinit + el_length*(e + s_fraction[0]);
210 
211  // Increment the node counter
212  node_count++;
213  }
214  } // End of loop over central elements
215 
216 
217  // FINAL ELEMENT
218  //--------------
219 
220  // Allocate memory for element
221  Element_pt[N-1] = new ELEMENT;
222 
223  // The first node is the same as the last node in the neighbouring
224  // element on the left
225  finite_element_pt(N-1)->node_pt(0)
226  = finite_element_pt(N-2)->node_pt(n_node-1);
227 
228  // New central nodes (ignore last one which needs special treatment
229  // because it's on the boundary)
230  for(unsigned jnod=1;jnod<(n_node-1);jnod++)
231  {
232  // Allocate memory for nodes, as before
233  Node_pt[node_count]
234  = finite_element_pt(N-1)->construct_node(jnod,time_stepper_pt);
235 
236  // Get the local coordinate of the node
237  finite_element_pt(N-1)->local_fraction_of_node(jnod,s_fraction);
238 
239  // Set the position of the node
240  node_pt(node_count)->x(0) = xinit + el_length*(N-1 + s_fraction[0]);
241 
242  // Increment the node counter
243  node_count++;
244  }
245 
246  // New final node
247  // --------------
248 
249  // Allocate memory for the node, as before
250  Node_pt[node_count] = finite_element_pt(N-1)
251  ->construct_boundary_node(n_node-1,time_stepper_pt);
252 
253  // Set the position of the node
254  node_pt(node_count)->x(0) = xinit + Length;
255 
256  // Add the node to the boundary 1
257  add_boundary_node(1,Node_pt[node_count]);
258 
259  // Increment the node counter
260  node_count++;
261  }
262  }
263 
264 } // End of namespace
265 
266 #endif
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh constuction routine, called by all constructors.
e
Definition: cfortran.h:575
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219