xda_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_XDA_TET_MESH_TEMPLATE_CC
31 #define OOMPH_XDA_TET_MESH_TEMPLATE_CC
32 
33 
34 #include "../generic/Telements.h"
35 #include "xda_tet_mesh.template.h"
36 
37 
38 namespace oomph
39 {
40 
41 
42  //======================================================================
43  /// \short Constructor: Pass name of xda file. Note that all
44  /// boundary elements get their own ID -- this is required for
45  /// FSI problems. The vector containing the oomph-lib
46  /// boundary IDs of all oomph-lib boundaries that collectively form
47  /// a given boundary as specified in the xda input file can be
48  /// obtained from the access function oomph_lib_boundary_ids(...).
49  /// Timestepper defaults to steady pseudo-timestepper.
50  //======================================================================
51  template <class ELEMENT>
53  TimeStepper* time_stepper_pt)
54  {
55  // Mesh can only be built with 3D Telements.
56  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
57 
58  // Open and process xda input file
59  std::ifstream infile(xda_file_name.c_str(),std::ios_base::in);
60  unsigned n_node;
61  unsigned n_element;
62  unsigned n_bound_face;
63 
64  if (!infile.is_open())
65  {
66  std::ostringstream error_stream;
67  error_stream
68  << "Failed to open " << xda_file_name
69  << "\n";
70  throw OomphLibError(error_stream.str(),
71  OOMPH_CURRENT_FUNCTION,
72  OOMPH_EXCEPTION_LOCATION);
73  }
74 
75  // Dummy storage to jump lines
76  char dummy[101];
77 
78  // Ignore file format
79  infile.getline(dummy, 100);
80 
81  // Get number of elements
82  infile>>n_element;
83 
84  // Ignore rest of line
85  infile.getline(dummy, 100);
86 
87  // Get number of nodes
88  infile>>n_node;
89 
90  // Ignore rest of line
91  infile.getline(dummy, 100);
92 
93  // Ignore sum of element weights (whatever that is...)
94  infile.getline(dummy, 100);
95 
96  // Get number of enumerated boundary faces on which boundary conditions
97  // are applied.
98  infile>>n_bound_face;
99 
100  // Keep reading until "Title String"
101  while (dummy[0]!= 'T') {infile.getline(dummy, 100);}
102 
103  // Make space for nodes and elements
104  Node_pt.resize(n_node);
105  Element_pt.resize(n_element);
106 
107  // Read first line with node labels and count them
108  std::string line;
109  std::getline(infile,line);
110  std::istringstream ostr(line);
111  std::istream_iterator<std::string> it(ostr);
112  std::istream_iterator<std::string> end;
113  unsigned nnod_el = 0;
114  Vector<unsigned> first_node;
115  while (it != end)
116  {
117  first_node.push_back(atoi((*it).c_str()));
118  it++;
119  nnod_el++;
120  }
121 
122  // Check
123  if (nnod_el!=10)
124  {
125  std::ostringstream error_stream;
126  error_stream
127  << "XdaTetMesh can currently only be built with quadratic tets "
128  << "with 10 nodes. The specified mesh file has " << nnod_el
129  << "nodes per element.\n";
130  throw OomphLibError(error_stream.str(),
131  OOMPH_CURRENT_FUNCTION,
132  OOMPH_EXCEPTION_LOCATION);
133  }
134 
135  // Storage for the global node numbers listed element-by-element
136  Vector<unsigned> global_node(n_element*nnod_el);
137 
138  // Read in nodes
139  unsigned k=0;
140 
141  // Copy across first nodes
142  for(unsigned j=0;j<nnod_el;j++)
143  {
144  global_node[k]=first_node[k];
145  k++;
146  }
147 
148  // Read the other ones
149  for(unsigned i=1;i<n_element;i++)
150  {
151  for(unsigned j=0;j<nnod_el;j++)
152  {
153  infile >> global_node[k];
154  k++;
155  }
156  infile.getline(dummy, 100);
157  }
158 
159 
160  // Create storage for coordinates
161  Vector<double> x_node(n_node);
162  Vector<double> y_node(n_node);
163  Vector<double> z_node(n_node);
164 
165  // Get nodal coordinates
166  for(unsigned i=0;i<n_node;i++)
167  {
168  infile>>x_node[i];
169  infile>>y_node[i];
170  infile>>z_node[i];
171  }
172 
173 
174  // Read in boundaries for faces
175  unsigned element_nmbr;
176  unsigned bound_id;
177  unsigned side_nmbr;
178  unsigned max_bound=0;
179 
180  // Make space for enumeration of sub-boundaries
181  Boundary_id.resize(n_bound_face+1);
182 
183  // Counter for number of separate boundary faces
184  unsigned count=0;
185  Vector<std::set<unsigned> > boundary_id(n_node);
186  for(unsigned i=0;i<n_bound_face;i++)
187  {
188  // Number of the element
189  infile>> element_nmbr ;
190 
191  // Which side/face on the tet are we dealing with (xda enumeratation)?
192  infile>> side_nmbr ;
193 
194  // What's the boundary ID?
195  infile>> bound_id ;
196 
197  // Turn into zero-based oomph-lib mesh boundary id
198  unsigned oomph_lib_bound_id=bound_id-1;
199  oomph_lib_bound_id=count;
200  Boundary_id[bound_id].push_back(count);
201 
202  // Increment number of separate boundary faces
203  count++;
204 
205  // Get ready for allocation of total number of boundaries
206  if (oomph_lib_bound_id>max_bound) max_bound=oomph_lib_bound_id;
207 
208  // Identify the "side nodes" (i.e. the nodes on the faces of
209  // the bulk tet) according to the
210  // conventions in '.xda' mesh files so that orientation of the
211  // faces is always the same (required for computation of
212  // outer unit normals
213  Vector<unsigned> side_node(6);
214  switch(side_nmbr)
215  {
216  case 0:
217  side_node[0]=global_node[nnod_el*element_nmbr+1];
218  side_node[1]=global_node[nnod_el*element_nmbr];
219  side_node[2]=global_node[nnod_el*element_nmbr+2];
220  side_node[3]=global_node[nnod_el*element_nmbr+4];
221  side_node[4]=global_node[nnod_el*element_nmbr+6];
222  side_node[5]=global_node[nnod_el*element_nmbr+5];
223  break;
224 
225  case 1:
226  side_node[0]=global_node[nnod_el*element_nmbr];
227  side_node[1]=global_node[nnod_el*element_nmbr+1];
228  side_node[2]=global_node[nnod_el*element_nmbr+3];
229  side_node[3]=global_node[nnod_el*element_nmbr+4];
230  side_node[4]=global_node[nnod_el*element_nmbr+8];
231  side_node[5]=global_node[nnod_el*element_nmbr+7];
232  break;
233 
234  case 2:
235  side_node[0]=global_node[nnod_el*element_nmbr+1];
236  side_node[1]=global_node[nnod_el*element_nmbr+2];
237  side_node[2]=global_node[nnod_el*element_nmbr+3];
238  side_node[3]=global_node[nnod_el*element_nmbr+5];
239  side_node[4]=global_node[nnod_el*element_nmbr+9];
240  side_node[5]=global_node[nnod_el*element_nmbr+8];
241  break;
242 
243  case 3:
244  side_node[0]=global_node[nnod_el*element_nmbr+2];
245  side_node[1]=global_node[nnod_el*element_nmbr];
246  side_node[2]=global_node[nnod_el*element_nmbr+3];
247  side_node[3]=global_node[nnod_el*element_nmbr+6];
248  side_node[4]=global_node[nnod_el*element_nmbr+7];
249  side_node[5]=global_node[nnod_el*element_nmbr+9];
250  break;
251 
252  default :
253  throw OomphLibError(
254  "Unexcpected side number in your '.xda' input file\n",
255  OOMPH_CURRENT_FUNCTION,
256  OOMPH_EXCEPTION_LOCATION);
257  }
258 
259  // Associate boundaries with nodes
260  for (unsigned j=0;j<6;j++)
261  {
262  boundary_id[side_node[j]].insert(oomph_lib_bound_id);
263  }
264  }
265 
266  // Set number of boundaries
267  set_nboundary(max_bound+1);
268 
269  // Vector of bools, will tell us if we already visited a node
270  std::vector<bool> done(n_node,false);
271 
272  Vector<unsigned> translate(n_node);
273  translate[0]=0;
274  translate[1]=2;
275  translate[2]=1;
276  translate[3]=3;
277  translate[4]=5;
278  translate[5]=7;
279  translate[6]=4;
280  translate[7]=6;
281  translate[8]=8;
282  translate[9]=9;
283 
284  // Create the elements
285  unsigned node_count=0;
286  for(unsigned e=0;e<n_element;e++)
287  {
288  Element_pt[e]=new ELEMENT;
289 
290  // Loop over all nodes
291  for(unsigned j=0;j<nnod_el;j++)
292  {
293  unsigned j_global=global_node[node_count];
294  if(done[j_global]==false)
295  {
296  if (boundary_id[j_global].size()==0)
297  {
298  Node_pt[j_global]=
299  finite_element_pt(e)->construct_node(translate[j],time_stepper_pt);
300  }
301  else
302  {
303  Node_pt[j_global]=
304  finite_element_pt(e)->construct_boundary_node(translate[j],
305  time_stepper_pt);
306  for (std::set<unsigned>::iterator it=boundary_id[j_global].begin();
307  it!=boundary_id[j_global].end();it++)
308  {
309  add_boundary_node(*it,Node_pt[j_global]);
310  }
311  }
312  done[j_global]=true;
313  Node_pt[j_global]->x(0)=x_node[j_global];
314  Node_pt[j_global]->x(1)=y_node[j_global];
315  Node_pt[j_global]->x(2)=z_node[j_global];
316  }
317  else
318  {
319  finite_element_pt(e)->node_pt(translate[j]) = Node_pt[j_global];
320  }
321  node_count++;
322  }
323  }
324 
325  // Figure out which elements are next to the boundaries.
326  setup_boundary_element_info();
327 
328 
329  // Setup boundary coordinates
330  unsigned nb=nboundary();
331  for (unsigned b=0;b<nb;b++)
332  {
333  bool switch_normal=false;
334  setup_boundary_coordinates(b,switch_normal);
335  }
336 
337  }
338 
339 
340 //======================================================================
341 /// Setup boundary coordinate on boundary b while is
342 /// temporarily flattened to simplex faces. Boundary coordinates are the
343 /// x-y coordinates in the plane of that boundary with the
344 /// x-axis along the line from the (lexicographically)
345 /// "lower left" to the "upper right" node. The y axis
346 /// is obtained by taking the cross-product of the positive
347 /// x direction with the outer unit normal computed by
348 /// the face elements (or its negative if switch_normal is set
349 /// to true). Doc faces in output file.
350 //======================================================================
351  template <class ELEMENT>
353  const bool& switch_normal,
354  std::ofstream& outfile)
355  {
356 
357  // Temporary storage for face elements
358  Vector<FiniteElement*> face_el_pt;
359 
360  // Backup for nodal positions
361  std::map<Node*,Vector<double> > backup_position;
362 
363  // Loop over all elements on boundaries
364  unsigned nel=this->nboundary_element(b);
365  if (nel>0)
366  {
367  // Loop over the bulk elements adjacent to boundary b
368  for(unsigned e=0;e<nel;e++)
369  {
370  // Get pointer to the bulk element that is adjacent to boundary b
371  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b,e);
372 
373  //Find the index of the face of element e along boundary b
374  int face_index = this->face_index_at_boundary(b,e);
375 
376  // Create new face element
378  new DummyFaceElement<ELEMENT>(bulk_elem_pt,face_index);
379  face_el_pt.push_back(el_pt);
380 
381  // Backup nodal position
382  Vector<double> pos(3);
383  for (unsigned j=3;j<6;j++)
384  {
385  if (backup_position[el_pt->node_pt(j)].size()==0)
386  {
387  el_pt->node_pt(j)->position(pos);
388  backup_position[el_pt->node_pt(j)]=pos;
389  }
390  }
391 
392  // Temporarily flatten the element to a simplex
393  for (unsigned i=0;i<3;i++)
394  {
395  // Node 3 is between vertex nodes 0 and 1
396  el_pt->node_pt(3)->x(i)=
397  0.5*(el_pt->node_pt(0)->x(i)+el_pt->node_pt(1)->x(i));
398 
399  // Node 4 is between vertex nodes 1 and 2
400  el_pt->node_pt(4)->x(i)=
401  0.5*(el_pt->node_pt(1)->x(i)+el_pt->node_pt(2)->x(i));
402 
403  // Node 5 is between vertex nodes 2 and 0
404  el_pt->node_pt(5)->x(i)=
405  0.5*(el_pt->node_pt(2)->x(i)+el_pt->node_pt(0)->x(i));
406  }
407 
408 
409  // Output faces?
410  if (outfile.is_open())
411  {
412  face_el_pt[face_el_pt.size()-1]->output(outfile);
413  }
414  }
415 
416 
417  // Loop over all nodes to find the lower left and upper right ones
418  Node* lower_left_node_pt=this->boundary_node_pt(b,0);
419  Node* upper_right_node_pt=this->boundary_node_pt(b,0);
420 
421  // Loop over all nodes on the boundary
422  unsigned nnod=this->nboundary_node(b);
423  for (unsigned j=0;j<nnod;j++)
424  {
425 
426  //Get node
427  Node* nod_pt=this->boundary_node_pt(b,j);
428 
429  // Primary criterion for lower left: Does it have a lower z-coordinate?
430  if (nod_pt->x(2)<lower_left_node_pt->x(2))
431  {
432  lower_left_node_pt=nod_pt;
433  }
434  // ...or is it a draw?
435  else if (nod_pt->x(2)==lower_left_node_pt->x(2))
436  {
437  // If it's a draw: Does it have a lower y-coordinate?
438  if (nod_pt->x(1)<lower_left_node_pt->x(1))
439  {
440  lower_left_node_pt=nod_pt;
441  }
442  // ...or is it a draw?
443  else if (nod_pt->x(1)==lower_left_node_pt->x(1))
444  {
445 
446  // If it's a draw: Does it have a lower x-coordinate?
447  if (nod_pt->x(0)<lower_left_node_pt->x(0))
448  {
449  lower_left_node_pt=nod_pt;
450  }
451  }
452  }
453 
454  // Primary criterion for upper right: Does it have a higher
455  // z-coordinate?
456  if (nod_pt->x(2)>upper_right_node_pt->x(2))
457  {
458  upper_right_node_pt=nod_pt;
459  }
460  // ...or is it a draw?
461  else if (nod_pt->x(2)==upper_right_node_pt->x(2))
462  {
463  // If it's a draw: Does it have a higher y-coordinate?
464  if (nod_pt->x(1)>upper_right_node_pt->x(1))
465  {
466  upper_right_node_pt=nod_pt;
467  }
468  // ...or is it a draw?
469  else if (nod_pt->x(1)==upper_right_node_pt->x(1))
470  {
471 
472  // If it's a draw: Does it have a higher x-coordinate?
473  if (nod_pt->x(0)>upper_right_node_pt->x(0))
474  {
475  upper_right_node_pt=nod_pt;
476  }
477  }
478  }
479  }
480 
481  // Prepare storage for boundary coordinates
482  Vector<double> zeta(2);
483 
484  // Unit vector connecting lower left and upper right nodes
485  Vector<double> b0(3);
486  b0[0]=upper_right_node_pt->x(0)-lower_left_node_pt->x(0);
487  b0[1]=upper_right_node_pt->x(1)-lower_left_node_pt->x(1);
488  b0[2]=upper_right_node_pt->x(2)-lower_left_node_pt->x(2);
489 
490  // Normalise
491  double inv_length=1.0/sqrt(b0[0]*b0[0]+b0[1]*b0[1]+b0[2]*b0[2]);
492  b0[0]*=inv_length;
493  b0[1]*=inv_length;
494  b0[2]*=inv_length;
495 
496  // Get (outer) unit normal to first face element
497  Vector<double> normal(3);
498  Vector<double> s(2,0.0);
499  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])->
500  outer_unit_normal(s,normal);
501 
502  if (switch_normal)
503  {
504  normal[0]=-normal[0];
505  normal[1]=-normal[1];
506  normal[2]=-normal[2];
507  }
508 
509 #ifdef PARANOID
510 
511  // Check that all elements have the same normal
512  for(unsigned e=0;e<nel;e++)
513  {
514  // Get (outer) unit normal to face element
515  Vector<double> my_normal(3);
516  dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])->
517  outer_unit_normal(s,my_normal);
518 
519  // Dot product should be one!
520  double error=
521  normal[0]*my_normal[0]+
522  normal[1]*my_normal[1]+
523  normal[2]*my_normal[2];
524 
525  error-=1.0;
526  if (switch_normal) error+=1.0;
527 
528  if (error>Tolerance_for_boundary_finding)
529  {
530  std::ostringstream error_message;
531  error_message
532  << "Error in alingment of normals (dot product-1) "
533  << error << " for element " << e << std::endl
534  << "exceeds tolerance specified by the static member data\n"
535  << "TetMeshBase::Tolerance_for_boundary_finding = "
536  << Tolerance_for_boundary_finding << std::endl
537  << "This usually means that the boundary is not planar.\n\n"
538  << "You can suppress this error message by recompiling \n"
539  << "recompiling without PARANOID or by changing the tolerance.\n";
540  throw OomphLibError(error_message.str(),
541  OOMPH_CURRENT_FUNCTION,
542  OOMPH_EXCEPTION_LOCATION);
543  }
544  }
545 #endif
546 
547  // Cross-product to get second in-plane vector, normal to b0
548  Vector<double> b1(3);
549  b1[0]=b0[1]*normal[2]-b0[2]*normal[1];
550  b1[1]=b0[2]*normal[0]-b0[0]*normal[2];
551  b1[2]=b0[0]*normal[1]-b0[1]*normal[0];
552 
553  // Assign boundary coordinates: projection onto the axes
554  for (unsigned j=0;j<nnod;j++)
555  {
556  //Get node
557  Node* nod_pt=this->boundary_node_pt(b,j);
558 
559  // Difference vector to lower left corner
560  Vector<double> delta(3);
561  delta[0]=nod_pt->x(0)-lower_left_node_pt->x(0);
562  delta[1]=nod_pt->x(1)-lower_left_node_pt->x(1);
563  delta[2]=nod_pt->x(2)-lower_left_node_pt->x(2);
564 
565  // Project
566  zeta[0]=delta[0]*b0[0]+delta[1]*b0[1]+delta[2]*b0[2];
567  zeta[1]=delta[0]*b1[0]+delta[1]*b1[1]+delta[2]*b1[2];
568 
569 #ifdef PARANOID
570 
571  // Check:
572  double error=
573  pow(lower_left_node_pt->x(0)+zeta[0]*b0[0]+zeta[1]*b1[0]-
574  nod_pt->x(0),2)+
575  pow(lower_left_node_pt->x(1)+zeta[0]*b0[1]+zeta[1]*b1[1]-
576  nod_pt->x(1),2)+
577  pow(lower_left_node_pt->x(2)+zeta[0]*b0[2]+zeta[1]*b1[2]-
578  nod_pt->x(2),2);
579 
580  if (sqrt(error)>Tolerance_for_boundary_finding)
581  {
582  std::ostringstream error_message;
583  error_message
584  << "Error in setup of boundary coordinate "
585  << sqrt(error) << std::endl
586  << "exceeds tolerance specified by the static member data\n"
587  << "TetMeshBase::Tolerance_for_boundary_finding = "
588  << Tolerance_for_boundary_finding << std::endl
589  << "This usually means that the boundary is not planar.\n\n"
590  << "You can suppress this error message by recompiling \n"
591  << "recompiling without PARANOID or by changing the tolerance.\n";
592  throw OomphLibError(error_message.str(),
593  OOMPH_CURRENT_FUNCTION,
594  OOMPH_EXCEPTION_LOCATION);
595  }
596 #endif
597 
598  // Set boundary coordinate
599  nod_pt->set_coordinates_on_boundary(b,zeta);
600  }
601 
602  }
603 
604  // Indicate that boundary coordinate has been set up
605  Boundary_coordinate_exists[b]=true;
606 
607  // Cleanup
608  unsigned n=face_el_pt.size();
609  for (unsigned e=0;e<n;e++)
610  {
611  delete face_el_pt[e];
612  }
613 
614 
615  // Reset nodal position
616  for (std::map<Node*,Vector<double> >::iterator it=backup_position.begin();
617  it!=backup_position.end();it++)
618  {
619  Node* nod_pt=(*it).first;
620  Vector<double> pos((*it).second);
621  for (unsigned i=0;i<3;i++)
622  {
623  nod_pt->x(i)=pos[i];
624  }
625  }
626 
627  }
628 
629 
630 }
631 
632 
633 #endif
void setup_boundary_coordinates(const unsigned &b, const bool &switch_normal)
Setup boundary coordinate on boundary b while is temporarily flattened to simplex faces...
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2301
cstr elem_len * i
Definition: cfortran.h:607
A general Finite Element class.
Definition: elements.h:1271
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
XdaTetMesh(const std::string xda_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass name of xda file. Note that all boundary elements get their own ID – this is requir...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
static char t char * s
Definition: cfortran.h:572
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.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void position(Vector< double > &pos) const
Compute Vector of nodal positions either directly or via hanging node representation.
Definition: nodes.cc:2399