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