geompack_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_GEOMPACK_MESH_TEMPLATE_CC
27 #define OOMPH_GEOMPACK_MESH_TEMPLATE_CC
28 
29 #include "geompack_mesh.template.h"
30 
31 
32 namespace oomph
33 {
34  //========================================================================
35  /// Quadrilateral mesh generator; Uses input from Geompack++.
36  /// See: http://members.shaw.ca/bjoe/
37  /// Currently only for four-noded quads -- extension to higher-order
38  /// quads should be trivial (see the corresponding classes for
39  /// triangular meshes).
40  //========================================================================
41  template<class ELEMENT>
43  TimeStepper* time_stepper_pt)
44  {
45  // Mesh can only be built with four-noded 2D Qelements.
46  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2, 2);
47 
48  // Create space for elements
49  unsigned nelem = Tmp_mesh_pt->nelement();
50  Element_pt.resize(nelem);
51 
52  // Create space for nodes
53  unsigned nnod = Tmp_mesh_pt->nnode();
54  Node_pt.resize(nnod);
55 
56  // Set number of boundaries
57  unsigned nbound = Tmp_mesh_pt->nboundary();
58  set_nboundary(nbound);
59 
60  // Loop over elements in scaffold mesh, visit their nodes
61  for (unsigned e = 0; e < nelem; e++)
62  {
63  Element_pt[e] = new ELEMENT;
64  }
65 
66 
67  // At the moment we can only do four-noded quads
68  if (finite_element_pt(0)->nnode() != 4)
69  {
70  std::string error_message =
71  "GeompackQuadMesh can currently only deal with\n";
72  error_message += "four noded quads! \n";
73  error_message += "Generalisation to higher-order elements should be \n";
74  error_message += "straightforward but hasn't been done yet.\n";
75  error_message += "Do you want to volunteer? \n";
76 
77  throw OomphLibError(
78  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
79  }
80 
81 
82  // In the first instance build all nodes from within all the elements
83  unsigned nnod_el = Tmp_mesh_pt->finite_element_pt(0)->nnode();
84 
85  // Loop over elements in scaffold mesh, visit their nodes
86  for (unsigned e = 0; e < nelem; e++)
87  {
88  // Loop over all nodes in element
89  for (unsigned j = 0; j < nnod_el; j++)
90  {
91  // Create new node, using the NEW element's construct_node
92  // member function
93  finite_element_pt(e)->construct_node(j, time_stepper_pt);
94  }
95  }
96 
97 
98  std::map<Node*, unsigned> global_number;
99  unsigned global_count = 0;
100  // Loop over elements in scaffold mesh, visit their nodes
101  for (unsigned e = 0; e < nelem; e++)
102  {
103  // Loop over all nodes in element
104  for (unsigned j = 0; j < nnod_el; j++)
105  {
106  // Pointer to node in the scaffold mesh
107  Node* scaffold_node_pt = Tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
108 
109  // Get the (pseudo-)global node number in scaffold mesh
110  // (It's zero [=default] if not visited this one yet)
111  unsigned j_global = global_number[scaffold_node_pt];
112 
113  // Haven't done this one yet
114  if (j_global == 0)
115  {
116  // Give it a number (not necessarily the global node
117  // number in the scaffold mesh -- we just need something
118  // to keep track...)
119  global_count++;
120  global_number[scaffold_node_pt] = global_count;
121 
122  // Copy new node, created using the NEW element's construct_node
123  // function into global storage, using the same global
124  // node number that we've just associated with the
125  // corresponding node in the scaffold mesh
126  Node_pt[global_count - 1] = finite_element_pt(e)->node_pt(j);
127 
128  // Assign coordinates
129  Node_pt[global_count - 1]->x(0) = scaffold_node_pt->x(0);
130  Node_pt[global_count - 1]->x(1) = scaffold_node_pt->x(1);
131 
132 
133  // Get pointer to set of mesh boundaries that this
134  // scaffold node occupies; NULL if the node is not on any boundary
135  std::set<unsigned>* boundaries_pt;
136  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
137 
138  // Loop over the mesh boundaries that the node in the scaffold mesh
139  // occupies and assign new node to the same ones.
140  if (boundaries_pt != 0)
141  {
142  // Convert it to a boundary node
143  this->convert_to_boundary_node(Node_pt[global_count - 1]);
144  // Add the node to the boundaries
145  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
146  it != boundaries_pt->end();
147  ++it)
148  {
149  add_boundary_node(*it, Node_pt[global_count - 1]);
150  }
151  }
152  }
153  // This one has already been done: Kill it
154  else
155  {
156  // Kill it
157  delete finite_element_pt(e)->node_pt(j);
158 
159  // Overwrite the element's pointer to local node by
160  // pointer to the already existing node -- identified by
161  // the number kept in the map
162  finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
163  }
164  }
165  }
166  }
167 
168 } // namespace oomph
169 #endif
e
Definition: cfortran.h:571
void build_from_scaffold(TimeStepper *time_stepper_pt)
Build mesh from scaffold.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
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
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...