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