geompack_scaffold_mesh.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 #include "geompack_scaffold_mesh.h"
31 
32 namespace oomph
33 {
34 
35 //=====================================================================
36 /// \short Constructor: Pass the filename of the mesh files
37 //=====================================================================
39  const std::string& mesh_file_name,
40  const std::string& curve_file_name)
41 {
42 
43  //Read the output mesh file to find informations about the nodes
44  //and elements of the mesh
45 
46  // Process mesh file
47  //---------------------
48  std::ifstream mesh_file(mesh_file_name.c_str(),std::ios_base::in);
49 
50  // Number of nodes
51  unsigned n_node;
52  mesh_file>>n_node;
53 
54  // Coordinates of nodes and extra information index
55  Vector<double> x_node(n_node);
56  Vector<double> y_node(n_node);
57  Vector<int> vertinfo(n_node);
58  for(unsigned i=0;i<n_node;i++)
59  {
60  mesh_file>>x_node[i];
61  mesh_file>>y_node[i];
62  mesh_file>>vertinfo[i];
63  }
64 
65  // Extra information (nodes lying on a curve)
66  unsigned n_vx;
67  mesh_file>>n_vx;
68 
69  // Dummy storage for the node code
70  int dummy_node_code;
71 
72  // Storage for the curve indice
73  Vector<int> icurv(n_vx); // it's negative if not used
74 
75  // Dummy storage for the location of the point on the curve
76  double dummy_ucurv;
77 
78  for(unsigned i=0;i<n_vx;i++)
79  {
80  mesh_file>>dummy_node_code;
81  mesh_file>>icurv[i];
82  mesh_file>>dummy_ucurv;
83  }
84 
85  // Number of local nodes per element
86  unsigned n_local_node;
87  mesh_file>>n_local_node;
88 
89  // Number of elements
90  unsigned n_element;
91  mesh_file>>n_element;
92 
93  // Storage for global node numbers for all elements
94  Vector<unsigned> global_node(n_local_node*n_element);
95 
96  // Storage for edge information
97  // (needed for a possible construction of midside node
98  // in the following build from scaffold function)
99  Vector<int> edgeinfo(n_local_node*n_element);
100 
101  // Initialize counter
102  unsigned k=0;
103 
104  // Read global node numbers for all elements
105  for(unsigned i=0;i<n_element;i++)
106  {
107  for(unsigned j=0;j<n_local_node;j++)
108  {
109  mesh_file>>global_node[k];
110  k++;
111  }
112  }
113 
114  // Initialize counter
115  unsigned l=0;
116 
117  // Read the edge information
118  for(unsigned i=0;i<n_element;i++)
119  {
120  for(unsigned j=0;j<n_local_node;j++)
121  {
122  mesh_file>>edgeinfo[l];
123  l++;
124  }
125  }
126 
127  mesh_file.close();
128 
129  // Create a vector of boolean so as not to create the same node twice
130  std::vector<bool> done (n_node);
131  for(unsigned i=0;i<n_node;i++)
132  {
133  done[i]=false;
134  }
135 
136  // Resize the Node vector
137  Node_pt.resize(n_node,0);
138 
139  // Resize the Element vector
140  Element_pt.resize(n_element);
141 
142 
143  // Process curve file to extract information about boundaries
144  // ----------------------------------------------------------
145 
146  // Important: the input file must NOT have NURBS curve
147  std::ifstream curve_file(curve_file_name.c_str(),std::ios_base::in);
148 
149  // Number of curves
150  unsigned n_curv;
151  curve_file>>n_curv;
152 
153  // Storage of several information for each curve
154  Vector< Vector<int> > curv;
155 
156  // Resize to n_curv rows
157  curv.resize(n_curv);
158 
159  // Curve type
160  unsigned type;
161 
162  // Loop over the curves to extract information
163  for(unsigned i=0;i<n_curv;i++)
164  {
165  curve_file>>type;
166  if(type==1)
167  {
168  curv[i].resize(4);
169  curv[i][0]=type;
170  for(unsigned j=1;j<4;j++)
171  {
172  curve_file>>curv[i][j];
173  }
174  }
175  else if(type==2)
176  {
177  curv[i].resize(5);
178  curv[i][0]=type;
179  for(unsigned j=1;j<5;j++)
180  {
181  curve_file>>curv[i][j];
182  }
183  }
184  else
185  {
186  std::ostringstream error_stream;
187  error_stream << "Current we can only process curves of\n"
188  << "type 1 (straight lines) and 2 (circular arcs\n"
189  << "You've specified: type " << type << std::endl;
190 
191  throw OomphLibError(
192  error_stream.str(),
193  OOMPH_CURRENT_FUNCTION,
194  OOMPH_EXCEPTION_LOCATION);
195  }
196  }
197 
198  curve_file.close();
199 
200  // Searching the number of boundaries
201  int d=0;
202  for(unsigned i=0;i<n_curv;i++)
203  {
204  if(d<curv[i][1])
205  {
206  d=curv[i][1]; // the boundary code is the 2nd value of each curve
207  }
208  }
209  oomph_info<< "The number of boundaries is "<<d<<std::endl;
210 
211  // Set number of boundaries
212  if(d>0){ set_nboundary(d); }
213 
214  // Search the boundary number of node located on a boundary
215  // If after this, boundary_of_node[j][0] is -1 then node j
216  // is not located on any boundary.
217  // If boundary_of_node[j][0] is positive, the node is located
218  // on the boundary indicated by that number.
219  // If boundary_of_node[j][1] is also positive, the node is also
220  // located on that boundary. Note: We're ignoring the (remote)
221  // possibility that node is located on 3 or more boundaries
222  // as one of them would have to be an internal boundary which
223  // would be odd...
224  Vector<Vector<int> > boundary_of_node;
225  boundary_of_node.resize(n_node);
226  unsigned n;
227  for(unsigned i=0;i<n_node;i++)
228  {
229  n=0;
230  boundary_of_node[i].resize(2);
231  boundary_of_node[i][0]=-1;
232  boundary_of_node[i][1]=-1;
233  if(vertinfo[i]==2) // it's an endpoint
234  {
235  for(unsigned j=0;j<n_curv;j++)
236  {
237  for(unsigned m=2;m<curv[j].size();m++)
238  {
239  if(curv[j][m]==static_cast<int>(i+1)) // node number begins at 1
240  { //in the mesh file !!!
241  boundary_of_node[i][n]=curv[j][1];
242  n++;
243  }
244  }
245  }
246  }
247  if(vertinfo[i]>20)
248  {
249  int a=0;
250  a=(vertinfo[i])/20;
251  int b;
252  b=icurv[a-1]; // 1st value of vector at [0] !!
253  boundary_of_node[i][0]=curv[b-1][1]; // 1st value of vector at [0] !!
254  }
255  }
256 
257 
258  // Create the elements
259  //--------------------
260 
261  unsigned count=0;
262  unsigned c;
263  for(unsigned e=0;e<n_element;e++)
264  {
265  // Build simplex four node quad in the scaffold mesh
267 
268  // Construction of the two first nodes of the element
269  for(unsigned j=0;j<2;j++)
270  {
271  c=global_node[count];
272  if(done[c-1]==false) // c-1 because node number begins
273  // at 1 in the mesh file
274  {
275  //If the node is located on a boundary construct a boundary node
276  if((d>0) && ((boundary_of_node[c-1][0] > 0) ||
277  (boundary_of_node[c-1][1] > 0)))
278  {
279  //Construct a boundary node
281  //Add to the look=up schemes
282  if(boundary_of_node[c-1][0] > 0)
283  {add_boundary_node(boundary_of_node[c-1][0]-1,Node_pt[c-1]);}
284  if(boundary_of_node[c-1][1] > 0)
285  {add_boundary_node(boundary_of_node[c-1][1]-1,Node_pt[c-1]);}
286  }
287  //Otherwise construct a normal node
288  else
289  {
291  }
292  done[c-1]=true;
293  Node_pt[c-1]->x(0)=x_node[c-1];
294  Node_pt[c-1]->x(1)=y_node[c-1];
295  }
296  else
297  {
298  finite_element_pt(e)->node_pt(j) = Node_pt[c-1];
299  }
300  count++;
301  }
302 
303  // Construction of the third node not in anticlockwise order
304  c=global_node[count+1];
305  if(done[c-1]==false) //c-1 because node number begins at 1 in the mesh file
306  {
307  //If the node is on a boundary, construct a boundary node
308  if((d>0) && ((boundary_of_node[c-1][0]>0) ||
309  (boundary_of_node[c-1][1]>0)))
310  {
311  //Construct the node
313  //Add to the look-up schemes
314  if(boundary_of_node[c-1][0]>0)
315  {add_boundary_node(boundary_of_node[c-1][0]-1,Node_pt[c-1]);}
316  if(boundary_of_node[c-1][1]>0)
317  {add_boundary_node(boundary_of_node[c-1][1]-1,Node_pt[c-1]);}
318  }
319  //otherwise construct a normal node
320  else
321  {
323  }
324  done[c-1]=true;
325  Node_pt[c-1]->x(0)=x_node[c-1];
326  Node_pt[c-1]->x(1)=y_node[c-1];
327  }
328  else
329  {
330  finite_element_pt(e)->node_pt(2) = Node_pt[c-1];
331  }
332 
333  count++;
334 
335  // Construction of the fourth node
336  c=global_node[count-1];
337  if(done[c-1]==false) //c-1 because node number begins at 1 in the mesh file
338  {
339  //If the node is on a boundary, constuct a boundary node
340  if((d>0) && ((boundary_of_node[c-1][0]>0) ||
341  (boundary_of_node[c-1][1]>0)))
342  {
343  //Construct the boundary node
345  //Add to the look-up schemes
346  if(boundary_of_node[c-1][0]>0)
347  {add_boundary_node(boundary_of_node[c-1][0]-1,Node_pt[c-1]);}
348  if(boundary_of_node[c-1][1]>0)
349  {add_boundary_node(boundary_of_node[c-1][1]-1,Node_pt[c-1]);}
350  }
351  //Otherwise construct a normal node
352  else
353  {
355  }
356  done[c-1]=true;
357  Node_pt[c-1]->x(0)=x_node[c-1];
358  Node_pt[c-1]->x(1)=y_node[c-1];
359  }
360  else
361  {
362  finite_element_pt(e)->node_pt(3) = Node_pt[c-1];
363  }
364 
365  count++;
366  }
367 
368 }
369 
370 }
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:194
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
GeompackQuadScaffoldMesh()
Empty constructor.
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition: elements.h:2408
OomphInfo oomph_info
e
Definition: cfortran.h:575
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object. ...
Definition: elements.h:2380
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:197
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.