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