simple_cubic_tet_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_SIMPLE_CUBIC_TET_MESH_TEMPLATE_CC
27 #define OOMPH_SIMPLE_CUBIC_TET_MESH_TEMPLATE_CC
28 
29 #include <algorithm>
30 
31 // Simple 3D tetrahedral mesh class
32 #include "simple_cubic_tet_mesh.template.h"
33 #include "../generic/map_matrix.h"
34 #include <algorithm>
35 
36 
37 namespace oomph
38 {
39  //====================================================================
40  /// Simple tetrahedral mesh - with 24 tet elements constructed within a
41  /// "brick" form for each element block.
42  //====================================================================
43  template<class ELEMENT>
44  void SimpleCubicTetMesh<ELEMENT>::build_from_scaffold(
45  TimeStepper* time_stepper_pt)
46  {
47  // Mesh can only be built with 3D Telements.
48  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
49 
50  // Create space for elements
51  unsigned nelem = Tmp_mesh_pt->nelement();
52  Element_pt.resize(nelem);
53 
54  // Create space for nodes
55  unsigned nnode_scaffold = Tmp_mesh_pt->nnode();
56  Node_pt.resize(nnode_scaffold);
57 
58  // Set number of boundaries
59  unsigned nbound = Tmp_mesh_pt->nboundary();
60  set_nboundary(nbound);
61 
62  // Loop over elements in scaffold mesh, visit their nodes
63  for (unsigned e = 0; e < nelem; e++)
64  {
65  Element_pt[e] = new ELEMENT;
66  }
67 
68  // In the first instance build all nodes from within all the elements
69  unsigned nnod_el = Tmp_mesh_pt->finite_element_pt(0)->nnode();
70  // Loop over elements in scaffold mesh, visit their nodes
71  for (unsigned e = 0; e < nelem; e++)
72  {
73  // Loop over all nodes in element
74  for (unsigned j = 0; j < nnod_el; j++)
75  {
76  // Create new node, using the NEW element's construct_node
77  // member function
78  finite_element_pt(e)->construct_node(j, time_stepper_pt);
79  }
80  }
81 
82 
83  // Setup map to check the (pseudo-)global node number
84  // Nodes whose number is zero haven't been copied across
85  // into the mesh yet.
86  std::map<Node*, unsigned> global_number;
87  unsigned global_count = 0;
88  // Loop over elements in scaffold mesh, visit their nodes
89  for (unsigned e = 0; e < nelem; e++)
90  {
91  // Loop over all nodes in element
92  for (unsigned j = 0; j < nnod_el; j++)
93  {
94  // Pointer to node in the scaffold mesh
95  Node* scaffold_node_pt = Tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
96 
97  // Get the (pseudo-)global node number in scaffold mesh
98  // (It's zero [=default] if not visited this one yet)
99  unsigned j_global = global_number[scaffold_node_pt];
100 
101  // Haven't done this one yet
102  if (j_global == 0)
103  {
104  // Give it a number (not necessarily the global node
105  // number in the scaffold mesh -- we just need something
106  // to keep track...)
107  global_count++;
108  global_number[scaffold_node_pt] = global_count;
109 
110  // Copy new node, created using the NEW element's construct_node
111  // function into global storage, using the same global
112  // node number that we've just associated with the
113  // corresponding node in the scaffold mesh
114  Node_pt[global_count - 1] = finite_element_pt(e)->node_pt(j);
115 
116  // Assign coordinates
117  Node_pt[global_count - 1]->x(0) = scaffold_node_pt->x(0);
118  Node_pt[global_count - 1]->x(1) = scaffold_node_pt->x(1);
119  Node_pt[global_count - 1]->x(2) = scaffold_node_pt->x(2);
120 
121  // Get pointer to set of mesh boundaries that this
122  // scaffold node occupies; NULL if the node is not on any boundary
123  std::set<unsigned>* boundaries_pt;
124  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
125 
126  // Loop over the mesh boundaries that the node in the scaffold mesh
127  // occupies and assign new node to the same ones.
128  if (boundaries_pt != 0)
129  {
130  this->convert_to_boundary_node(Node_pt[global_count - 1]);
131  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
132  it != boundaries_pt->end();
133  ++it)
134  {
135  add_boundary_node(*it, Node_pt[global_count - 1]);
136  }
137  }
138  }
139  // This one has already been done: Kill it
140  else
141  {
142  // Kill it
143  delete finite_element_pt(e)->node_pt(j);
144 
145  // Overwrite the element's pointer to local node by
146  // pointer to the already existing node -- identified by
147  // the number kept in the map
148  finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
149  }
150  }
151  }
152 
153 
154  // At this point we've created all the elements and
155  // created their vertex nodes. Now we need to create
156  // the additional (midside and internal) nodes!
157 
158 
159  // We'll first create all local nodes for all elements
160  // and then delete the superfluous ones that have
161  // a matching node in an adjacent element.
162 
163  // Get number of nodes along element edge and dimension of element (3)
164  // from first element
165  unsigned nnode_1d = finite_element_pt(0)->nnode_1d();
166 
167  // At the moment we're only able to deal with nnode_1d=2 or 3.
168  if ((nnode_1d != 2) && (nnode_1d != 3))
169  {
170  std::ostringstream error_message;
171  error_message << "Mesh generation currently only works\n";
172  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
173  error_message << "for nnode_1d=" << nnode_1d << std::endl;
174 
175  throw OomphLibError(
176  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
177  }
178 
179  // Spatial dimension of element = number of local coordinates
180  unsigned dim = finite_element_pt(0)->dim();
181 
182  // Storage for the local coordinate of the new node
183  Vector<double> s(dim);
184 
185  // Get number of nodes in the element from first element
186  unsigned nnode = finite_element_pt(0)->nnode();
187 
188  // Loop over all elements
189  for (unsigned e = 0; e < nelem; e++)
190  {
191  // Loop over the new nodes in the element and create them.
192  for (unsigned j = 4; j < nnode; j++)
193  {
194  // Create new node
195  Node* new_node_pt =
196  finite_element_pt(e)->construct_node(j, time_stepper_pt);
197 
198  // What are the node's local coordinates?
199  finite_element_pt(e)->local_coordinate_of_node(j, s);
200 
201  // Find the coordinates of the new node from the existing
202  // and fully-functional element in the scaffold mesh
203  for (unsigned i = 0; i < dim; i++)
204  {
205  new_node_pt->x(i) =
206  Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s, i);
207  }
208  } // end of loop over new nodes
209  } // end of loop over elements
210 
211 
212  // Bracket this away so the edge map goes out of scope
213  // when we're done
214  {
215  // Storage for pointer to mid-edge node
216  MapMatrix<Node*, Node*> central_edge_node_pt;
217  Node* edge_node1_pt = 0;
218  Node* edge_node2_pt = 0;
219 
220  // Loop over elements
221  for (unsigned e = 0; e < nelem; e++)
222  {
223  // Loop over new local nodes
224  for (unsigned j = 4; j < nnode; j++)
225  {
226  // Pointer to the element's local node
227  Node* node_pt = finite_element_pt(e)->node_pt(j);
228 
229  // By default, we assume the node is not new
230  bool is_new = false;
231 
232  // This will have to be changed for higher-order elements
233  //=======================================================
234 
235  // Switch on local node number (all located on edges)
236  switch (j)
237  {
238  // Node 4 is located between nodes 0 and 1
239  case 4:
240 
241  edge_node1_pt = finite_element_pt(e)->node_pt(0);
242  edge_node2_pt = finite_element_pt(e)->node_pt(1);
243  if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
244  {
245  is_new = true;
246  central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
247  central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
248  }
249  break;
250 
251 
252  // Node 5 is located between nodes 0 and 2
253  case 5:
254 
255  edge_node1_pt = finite_element_pt(e)->node_pt(0);
256  edge_node2_pt = finite_element_pt(e)->node_pt(2);
257  if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
258  {
259  is_new = true;
260  central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
261  central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
262  }
263  break;
264 
265 
266  // Node 6 is located between nodes 0 and 3
267  case 6:
268 
269  edge_node1_pt = finite_element_pt(e)->node_pt(0);
270  edge_node2_pt = finite_element_pt(e)->node_pt(3);
271  if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
272  {
273  is_new = true;
274  central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
275  central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
276  }
277  break;
278 
279 
280  // Node 7 is located between nodes 1 and 2
281  case 7:
282 
283  edge_node1_pt = finite_element_pt(e)->node_pt(1);
284  edge_node2_pt = finite_element_pt(e)->node_pt(2);
285  if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
286  {
287  is_new = true;
288  central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
289  central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
290  }
291  break;
292 
293 
294  // Node 8 is located between nodes 2 and 3
295  case 8:
296 
297  edge_node1_pt = finite_element_pt(e)->node_pt(2);
298  edge_node2_pt = finite_element_pt(e)->node_pt(3);
299  if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
300  {
301  is_new = true;
302  central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
303  central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
304  }
305  break;
306 
307 
308  // Node 9 is located between nodes 1 and 3
309  case 9:
310 
311  edge_node1_pt = finite_element_pt(e)->node_pt(1);
312  edge_node2_pt = finite_element_pt(e)->node_pt(3);
313  if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
314  {
315  is_new = true;
316  central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
317  central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
318  }
319  break;
320 
321  default:
322  throw OomphLibError("More than ten nodes in Tet Element",
323  OOMPH_CURRENT_FUNCTION,
324  OOMPH_EXCEPTION_LOCATION);
325  }
326 
327  if (is_new)
328  {
329  // New node: Add it to mesh
330  Node_pt.push_back(node_pt);
331  }
332  else
333  {
334  // Delete local node in element...
335  delete finite_element_pt(e)->node_pt(j);
336 
337  // ... and reset pointer to the existing node
338  finite_element_pt(e)->node_pt(j) =
339  central_edge_node_pt(edge_node1_pt, edge_node2_pt);
340  }
341  }
342  }
343  }
344 
345 
346  // Boundary conditions
347 
348  // Matrix map to check if a node has already been added to
349  // the boundary number b (NOTE: Enumerated by pointer to ORIGINAL
350  // node before transfer to boundary node)
351  MapMatrixMixed<Node*, unsigned, bool> bound_node_done;
352 
353  // Create lookup scheme for pointers to original local nodes
354  // in elements
355  Vector<Vector<Node*>> orig_node_pt(nelem);
356 
357  // Loop over elements
358  for (unsigned e = 0; e < nelem; e++)
359  {
360  orig_node_pt[e].resize(nnode, 0);
361 
362  // Loop over new local nodes
363  for (unsigned j = 4; j < nnode; j++)
364  {
365  // Pointer to the element's local node
366  orig_node_pt[e][j] = finite_element_pt(e)->node_pt(j);
367  }
368  }
369 
370 
371  // Loop over elements
372  for (unsigned e = 0; e < nelem; e++)
373  {
374  // Loop over new local nodes
375  for (unsigned j = 4; j < nnode; j++)
376  {
377  // Loop over the boundaries
378  for (unsigned bo = 0; bo < nbound; bo++)
379  {
380  // Pointer to the element's local node
381  Node* loc_node_pt = finite_element_pt(e)->node_pt(j);
382 
383  // Pointer to original node
384  Node* orig_loc_node_pt = orig_node_pt[e][j];
385 
386  // value of the map for the node and boundary specified
387  bool bound_test = bound_node_done(orig_loc_node_pt, bo);
388 
389  if (bound_test == false)
390  {
391  bound_node_done(orig_loc_node_pt, bo) = true;
392 
393  // This will have to be changed for higher-order elements
394  //=======================================================
395 
396  // Switch on local node number (all located on edges)
397  switch (j)
398  {
399  // Node 4 is located between nodes 0 and 1
400  case 4:
401 
402  if (finite_element_pt(e)->node_pt(0)->is_on_boundary(bo) &&
403  finite_element_pt(e)->node_pt(1)->is_on_boundary(bo))
404  {
405  this->convert_to_boundary_node(loc_node_pt);
406  add_boundary_node(bo, loc_node_pt);
407  }
408  break;
409 
410 
411  // Node 5 is located between nodes 0 and 2
412  case 5:
413 
414  if (finite_element_pt(e)->node_pt(0)->is_on_boundary(bo) &&
415  finite_element_pt(e)->node_pt(2)->is_on_boundary(bo))
416  {
417  this->convert_to_boundary_node(loc_node_pt);
418  add_boundary_node(bo, loc_node_pt);
419  }
420  break;
421 
422 
423  // Node 6 is located between nodes 0 and 3
424  case 6:
425 
426  if (finite_element_pt(e)->node_pt(0)->is_on_boundary(bo) &&
427  finite_element_pt(e)->node_pt(3)->is_on_boundary(bo))
428  {
429  this->convert_to_boundary_node(loc_node_pt);
430  add_boundary_node(bo, loc_node_pt);
431  }
432  break;
433 
434 
435  // Node 7 is located between nodes 1 and 2
436  case 7:
437 
438  if (finite_element_pt(e)->node_pt(1)->is_on_boundary(bo) &&
439  finite_element_pt(e)->node_pt(2)->is_on_boundary(bo))
440  {
441  this->convert_to_boundary_node(loc_node_pt);
442  add_boundary_node(bo, loc_node_pt);
443  }
444  break;
445 
446 
447  // Node 8 is located between nodes 2 and 3
448  case 8:
449 
450  if (finite_element_pt(e)->node_pt(2)->is_on_boundary(bo) &&
451  finite_element_pt(e)->node_pt(3)->is_on_boundary(bo))
452  {
453  this->convert_to_boundary_node(loc_node_pt);
454  add_boundary_node(bo, loc_node_pt);
455  }
456  break;
457 
458 
459  // Node 9 is located between nodes 1 and 3
460  case 9:
461 
462  if (finite_element_pt(e)->node_pt(1)->is_on_boundary(bo) &&
463  finite_element_pt(e)->node_pt(3)->is_on_boundary(bo))
464  {
465  this->convert_to_boundary_node(loc_node_pt);
466  add_boundary_node(bo, loc_node_pt);
467  }
468  break;
469 
470 
471  default:
472  throw OomphLibError("More than ten nodes in Tet Element",
473  OOMPH_CURRENT_FUNCTION,
474  OOMPH_EXCEPTION_LOCATION);
475  }
476  }
477 
478  } // end for bo
479  } // end for j
480  } // end for e
481  }
482 
483 } // namespace oomph
484 
485 #endif
e
e
Definition: cfortran.h:571
s
static char t char * s
Definition: cfortran.h:568
i
cstr elem_len * i
Definition: cfortran.h:603
oomph::MapMatrixMixed
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
oomph::MapMatrix
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: map_matrix.h:508
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
oomph::Node::get_boundaries_pt
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1365
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition: oomph_definitions.h:222
oomph::SimpleCubicTetMesh::build_from_scaffold
void build_from_scaffold(TimeStepper *time_stepper_pt)
Build mesh from scaffold mesh.
Definition: simple_cubic_tet_mesh.template.cc:44
oomph::TimeStepper
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30
simple_cubic_tet_mesh.template.h