refineable_tetgen_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 
27 #ifndef OOMPH_REFINEABLE_TETGEN_MESH_TEMPLATE_CC
28 #define OOMPH_REFINEABLE_TETGEN_MESH_TEMPLATE_CC
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 // OOMPH-LIB Headers
37 #include "../generic/sample_point_parameters.h"
38 #include "../generic/mesh_as_geometric_object.h"
39 #include "../generic/projection.h"
40 #include "../generic/face_element_as_geometric_object.h"
41 
42 namespace oomph
43 {
44  /// Helper function that updates the input faceted surface
45  /// by using the flattened elements from FaceMesh(es) that are
46  /// constructed for the boundaries associated with the segments of the
47  /// polygon.
48  template<class ELEMENT>
50  TetMeshFacetedSurface*& faceted_surface_pt)
51  {
52  // The easiest thing to do is simply to update the
53  // positions of the key control nodes, leaving the connectivity alone,
54  // but that doesn't allow for any surface remeshing
55 
56  /// List of vertex nodes
57  std::list<Node*> new_vertex_nodes;
58  // List of facets and boundary ids
59  std::list<std::pair<Vector<unsigned>, unsigned>> new_facet_list;
60 
61  // Loop over the number of old facets
62  unsigned n_old_facet = faceted_surface_pt->nfacet();
63  for (unsigned f = 0; f < n_old_facet; f++)
64  {
65  // Get the boundary id of the facet. Need to subtract one,
66  // which is confusing now I think about it.
67  // ALH: Should fix this.
68  unsigned bound = faceted_surface_pt->one_based_facet_boundary_id(f) - 1;
69 
70  // Create a face mesh adjacent to the fluid mesh's bound-th boundary.
71  // The face mesh consists of FaceElements that may also be
72  // interpreted as GeomObjects
73  Mesh* face_mesh_pt = new Mesh;
74  this->template build_face_mesh<ELEMENT, FaceElementAsGeomObject>(
75  bound, face_mesh_pt);
76 
77  // Loop over these new face elements and tell them the boundary number
78  // from the bulk fluid mesh -- this is required to they can
79  // get access to the boundary coordinates!
80  unsigned n_face_element = face_mesh_pt->nelement();
81  for (unsigned e = 0; e < n_face_element; e++)
82  {
83  // Cast the element pointer to the correct thing!
85  dynamic_cast<FaceElementAsGeomObject<ELEMENT>*>(
86  face_mesh_pt->element_pt(e));
87 
88  // Set bulk boundary number
90  }
91 
92  // In order to set up the new faceted representation
93  // Need to know the positions of the corner nodes
94  // and the connectivity
95 
96  // Storage for the connectivity information
97  Vector<unsigned> new_local_facet(3);
98 
99  // Now we have the face mesh loop over the face elements and
100  // print out the end points
101  for (unsigned e = 0; e < n_face_element; ++e)
102  {
103  // Cache pointer to the element
104  FiniteElement const* elem_pt = face_mesh_pt->finite_element_pt(e);
105 
106  // Just use the three primary (corner) nodes to define the facet
107  unsigned n_vertex_node = 3;
108  for (unsigned n = 0; n < n_vertex_node; n++)
109  {
110  // Cache the pointer to the node
111  Node* const nod_pt = elem_pt->node_pt(n);
112  // If the vertex node is in the list return the number
113  unsigned counter = 0;
114  bool found_it = false;
115  for (std::list<Node*>::iterator it = new_vertex_nodes.begin();
116  it != new_vertex_nodes.end();
117  ++it, ++counter)
118  {
119  // If we have found the node then store it
120  if (*it == nod_pt)
121  {
122  new_local_facet[n] = counter;
123  found_it = true;
124  break;
125  }
126  }
127 
128  // If we haven't found it
129  // then add the node to the list and fill in the counter
130  if (!found_it)
131  {
132  new_vertex_nodes.push_back(nod_pt);
133  // We know where it is
134  new_local_facet[n] = counter;
135  }
136  }
137 
138  // Add the new facet connectivity to the list
139  new_facet_list.push_back(std::make_pair(new_local_facet, bound + 1));
140  }
141  } // End of loop over old facets
142 
143 
144  // Probably want the surface mesh in a better data structure so
145  // that we can perform refinement or coarsening on it
146  // i.e. want neighbours, edge flips all that fun stuff
147  // That will go here!
148 
149  // Now we need to put the information into the appropriate data structures
150  // for the Surface
151 
152  // Now sort out the facet nodes
153  unsigned n_facet_vertex = new_vertex_nodes.size();
154  // Vector<Vector<double> > facet_point(n_facet_vertex);
155  Vector<Node*> facet_nodes_pt(n_facet_vertex);
156  unsigned count = 0;
157  for (std::list<Node*>::iterator it = new_vertex_nodes.begin();
158  it != new_vertex_nodes.end();
159  ++it)
160  {
161  facet_nodes_pt[count] = *it;
162  // facet_point[count].resize(3);
163  // for(unsigned i=0;i<3;i++)
164  // {
165  // facet_point[count][i] = (*it)->x(i);
166  // }
167  ++count;
168  }
169 
170  // And also the facets
171  unsigned n_facet = new_facet_list.size();
172  Vector<Vector<unsigned>> new_facet(n_facet);
173  Vector<unsigned> new_facet_boundary_id(n_facet);
174  count = 0;
175  for (std::list<std::pair<Vector<unsigned>, unsigned>>::iterator it =
176  new_facet_list.begin();
177  it != new_facet_list.end();
178  ++it)
179  {
180  new_facet[count] = (*it).first;
181  new_facet_boundary_id[count] = (*it).second;
182  ++count;
183  }
184 
185  for (unsigned f = 0; f < n_facet; f++)
186  {
187  for (unsigned i = 0; i < 3; i++)
188  {
189  oomph_info << new_facet[f][i] << " ";
190  }
191  oomph_info << " : ";
192  oomph_info << new_facet_boundary_id[f] << "\n";
193  }
194 
195  // Now just create the new boundary
196  delete faceted_surface_pt;
197  faceted_surface_pt = new TetMeshFacetedClosedSurfaceForRemesh(
198  facet_nodes_pt, new_facet, new_facet_boundary_id);
199 
200  // Also need to setup the reverse look-up scheme again
201  this->setup_reverse_lookup_schemes_for_faceted_surface(faceted_surface_pt);
202 
203  // Take average to get a new hole position (Won't always work)
204  Vector<double> inner_point(3, 0.0);
205  for (unsigned n = 0; n < n_facet_vertex; n++)
206  {
207  for (unsigned i = 0; i < 3; i++)
208  {
209  inner_point[i] += facet_nodes_pt[n]->x(i);
210  }
211  }
212 
213  for (unsigned i = 0; i < 3; i++)
214  {
215  inner_point[i] /= n_facet_vertex;
216  }
217 
218  // Now set the hole if required
219  dynamic_cast<TetMeshFacetedClosedSurface*>(faceted_surface_pt)
220  ->set_hole_for_tetgen(inner_point);
221  }
222 
223 
224  /// Generate a new faceted representation of the inner hole
225  /// boundaries
226  template<class ELEMENT>
228  {
229  // Loop over the number of internal boundarys
230  unsigned n_hole = this->Internal_surface_pt.size();
231  for (unsigned ihole = 0; ihole < n_hole; ihole++)
232  {
233  // Now can the surface update its own representation goes in here
234 
235  // If not we have to generate it from the new mesh
236  {
237  // Update the faceted surface associated with the ihole-th hole
238  this->update_faceted_surface_using_face_mesh(
239  this->Internal_surface_pt[ihole]);
240  }
241  }
242  }
243 
244  /// Snap the boundary nodes onto any curvilinear boundaries
245  template<class ELEMENT>
247  RefineableTetgenMesh<ELEMENT>*& new_mesh_pt, const unsigned& b)
248  {
249  // Quick return
250  if (!Boundary_coordinate_exists[b])
251  {
252  oomph_info << "Not snapping nodes on boundary " << b
253  << " because no boundary coordinate has been set up.\n";
254  return;
255  }
256 
257  // Firstly we set the boundary coordinates of the new nodes
258  // In case the mapping between the geometric object's intrinsic coordiante
259  // and the arc-length coordinate is nonlinear.
260  // This is only an approximation,
261  // but it will ensure that the nodes that were input to triangle will
262  // retain exactly the same boundary coordinates and
263  // then linear interpolation
264  // is used between those values for any newly created nodes.
265 
266 
267  // Create a face mesh adjacent to the fluid mesh's b-th boundary.
268  // The face mesh consists of FaceElements that may also be
269  // interpreted as GeomObjects
270  Mesh* face_mesh_pt = new Mesh;
271  this->template build_face_mesh<ELEMENT, FaceElementAsGeomObject>(
272  b, face_mesh_pt);
273 
274  // Loop over these new face elements and tell them the boundary number
275  // from the bulk fluid mesh -- this is required to they can
276  // get access to the boundary coordinates!
277  unsigned n_face_element = face_mesh_pt->nelement();
278  for (unsigned e = 0; e < n_face_element; e++)
279  {
280  // Cast the element pointer to the correct thing!
282  dynamic_cast<FaceElementAsGeomObject<ELEMENT>*>(
283  face_mesh_pt->element_pt(e));
284 
285  // Set bulk boundary number
287  }
288 
289 
290  // Compare face meshes
291  {
292  Mesh* new_face_mesh_pt = new Mesh;
293  new_mesh_pt->template build_face_mesh<ELEMENT, FaceElementAsGeomObject>(
294  b, new_face_mesh_pt);
295 
296  std::ostringstream filestring;
297  filestring << "old_mesh_boundary" << b << ".dat";
298  face_mesh_pt->output(filestring.str().c_str());
299  filestring.clear();
300  filestring << "new_mesh_boundary" << b << ".dat";
301  new_face_mesh_pt->output(filestring.str().c_str());
302 
303  Vector<double> b_coo(2);
304  std::cout << "OLD---\n";
305  // Now let's look at the boundary coordinates
306  unsigned n_tmp_node = this->nboundary_node(b);
307  for (unsigned n = 0; n < n_tmp_node; ++n)
308  {
309  Node* const nod_pt = this->boundary_node_pt(b, n);
310  nod_pt->get_coordinates_on_boundary(b, b_coo);
311  std::cout << nod_pt->x(0) << " " << nod_pt->x(1) << " " << nod_pt->x(2)
312  << " " << b_coo[0] << " " << b_coo[1] << "\n";
313  }
314 
315  std::cout << "NEW-----------\n";
316  // Now let's look at the boundary coordinates
317  n_tmp_node = new_mesh_pt->nboundary_node(b);
318  for (unsigned n = 0; n < n_tmp_node; ++n)
319  {
320  Node* const nod_pt = new_mesh_pt->boundary_node_pt(b, n);
321  nod_pt->get_coordinates_on_boundary(b, b_coo);
322  std::cout << nod_pt->x(0) << " " << nod_pt->x(1) << " " << nod_pt->x(2)
323  << " " << b_coo[0] << " " << b_coo[1] << "\n";
324  }
325  }
326 
327 
328  // Now make the mesh as geometric object
329  MeshAsGeomObject* mesh_geom_obj_pt = new MeshAsGeomObject(face_mesh_pt);
330 
331  // Now assign the new nodes positions based on the old meshes
332  // potentially curvilinear boundary (its geom object incarnation)
333  Vector<double> new_x(3);
334  Vector<double> b_coord(2);
335  unsigned n_new_boundary_node = new_mesh_pt->nboundary_node(b);
336  for (unsigned n = 0; n < n_new_boundary_node; n++)
337  {
338  // Get the boundary coordinate of all new nodes
339  Node* const nod_pt = new_mesh_pt->boundary_node_pt(b, n);
340  nod_pt->get_coordinates_on_boundary(b, b_coord);
341  // Let's find boundary coordinates of the new node
342  mesh_geom_obj_pt->position(b_coord, new_x);
343  // Now snap to the boundary
344  for (unsigned i = 0; i < 3; i++)
345  {
346  nod_pt->x(i) = new_x[i];
347  }
348  }
349 
350 
351  // Delete the allocated memory for the geometric object and face mesh
352  delete mesh_geom_obj_pt;
353  face_mesh_pt->flush_element_and_node_storage();
354  delete face_mesh_pt;
355 
356  // Loop over the elements adjacent to the boundary
357  unsigned n_element = new_mesh_pt->nboundary_element(b);
358  if (n_element > 0)
359  {
360  // Make a dummy simplex element
361  TElement<3, 2> dummy_four_node_element;
362  // Make a dummy quadratic element
363  TElement<3, 3> dummy_ten_node_element;
364  Vector<double> s(3);
365  Vector<double> x_new(3);
366  for (unsigned n = 0; n < 4; n++)
367  {
368  dummy_four_node_element.construct_node(n);
369  }
370  for (unsigned n = 0; n < 10; n++)
371  {
372  dummy_ten_node_element.construct_node(n);
373  }
374  for (unsigned e = 0; e < n_element; e++)
375  {
376  // Cache the element pointer
377  ELEMENT* elem_pt =
378  dynamic_cast<ELEMENT*>(new_mesh_pt->boundary_element_pt(b, e));
379 
380  // Find the number of nodes
381  const unsigned n_node = elem_pt->nnode();
382  // Only do something if not simplex element
383  if (n_node > 4)
384  {
385  // Copy the nodes into the dummy
386  for (unsigned n = 0; n < 4; n++)
387  {
388  for (unsigned i = 0; i < 3; i++)
389  {
390  dummy_four_node_element.node_pt(n)->x(i) =
391  elem_pt->node_pt(n)->x(i);
392  }
393  }
394 
395  // Now sort out the mid-side nodes
396  for (unsigned n = 4; n < 10; n++)
397  {
398  // If it's not on a boundary then reset to be interpolated
399  // from the simplex
400  if (!elem_pt->node_pt(n)->is_on_boundary())
401  {
402  elem_pt->local_coordinate_of_node(n, s);
403  dummy_four_node_element.interpolated_x(s, x_new);
404  for (unsigned i = 0; i < 3; i++)
405  {
406  elem_pt->node_pt(n)->x(i) = x_new[i];
407  }
408  }
409  }
410  }
411 
412  // If we have more than 10 nodes interpolate from the quadratic shape
413  if (n_node > 10)
414  {
415  // Copy the nodes into the dummy
416  for (unsigned n = 0; n < 10; n++)
417  {
418  for (unsigned i = 0; i < 3; i++)
419  {
420  dummy_ten_node_element.node_pt(n)->x(i) =
421  elem_pt->node_pt(n)->x(i);
422  }
423  }
424 
425  // Now sort out the mid-face and central nodes
426  for (unsigned n = 10; n < n_node; n++)
427  {
428  // If it's not on a boundary then reset to be interpolated
429  // from the simplex
430  if (!elem_pt->node_pt(n)->is_on_boundary())
431  {
432  elem_pt->local_coordinate_of_node(n, s);
433  dummy_ten_node_element.interpolated_x(s, x_new);
434  for (unsigned i = 0; i < 3; i++)
435  {
436  elem_pt->node_pt(n)->x(i) = x_new[i];
437  }
438  }
439  }
440  }
441  }
442  } // End of fix up of elements
443  }
444 
445 
446  //======================================================================
447  /// Adapt problem based on specified elemental error estimates
448  //======================================================================
449  template<class ELEMENT>
451  {
452  double t_start = 0.0;
453  //###################################
454  t_start = TimingHelpers::timer();
455  //###################################
456 
457  // Get refinement targets
458  Vector<double> target_size(elem_error.size());
459  double max_edge_ratio =
460  this->compute_volume_target(elem_error, target_size);
461  // Get maximum target volume
462  unsigned n = target_size.size();
463  double max_size = 0.0;
464  double min_size = DBL_MAX;
465  for (unsigned e = 0; e < n; e++)
466  {
467  if (target_size[e] > max_size) max_size = target_size[e];
468  if (target_size[e] < min_size) min_size = target_size[e];
469  }
470 
471  oomph_info << "Maximum target size: " << max_size << std::endl;
472  oomph_info << "Minimum target size: " << min_size << std::endl;
473  oomph_info << "Number of elements to be refined " << this->Nrefined
474  << std::endl;
475  oomph_info << "Number of elements to be unrefined " << this->Nunrefined
476  << std::endl;
477  oomph_info << "Max edge ratio " << max_edge_ratio << std::endl;
478 
479  double orig_max_size, orig_min_size;
480  this->max_and_min_element_size(orig_max_size, orig_min_size);
481  oomph_info << "Max/min element size in original mesh: " << orig_max_size
482  << " " << orig_min_size << std::endl;
483 
484  //##################################################################
485  oomph_info
486  << "adapt: Time for getting volume targets : "
487  << TimingHelpers::timer() - t_start << " sec " << std::endl;
488  //##################################################################
489 
490  //===============================================================
491  // END: Compute target volumes
492  //===============================================================
493 
494  // Check if boundaries need to be updated (regardless of
495  // requirements of bulk error estimator) but don't do anything!
496  bool inner_boundary_update_necessary = false; // true;
497 
498  // Should we bother to adapt?
499  if ((Nrefined > 0) || (Nunrefined > this->max_keep_unrefined()) ||
500  (max_edge_ratio > this->max_permitted_edge_ratio()))
501  {
502  if (!((Nrefined > 0) || (Nunrefined > max_keep_unrefined())))
503  {
504  oomph_info << "Mesh regeneration triggered by edge ratio criterion\n";
505  }
506 
507  // Generate remesh of inner surface
508  if (inner_boundary_update_necessary)
509  {
510  this->surface_remesh_for_inner_hole_boundaries();
511 
512 
513  // If there is not a geometric object associated with the boundary
514  // the reset the boundary coordinates so that the lengths are consistent
515  // in the new mesh and the old mesh.
516  const unsigned n_boundary = this->nboundary();
517  for (unsigned b = 0; b < n_boundary; ++b)
518  {
519  // if(this->boundary_geom_object_pt(b)==0)
520  {
521  this->template setup_boundary_coordinates<ELEMENT>(b);
522  }
523  }
524  }
525 
526 
527  //###################################
528  t_start = TimingHelpers::timer();
529  //###################################
530 
531  // Are we dealing with a solid mesh?
532  SolidMesh* solid_mesh_pt = dynamic_cast<SolidMesh*>(this);
533 
534  // Build temporary uniform background mesh
535  //----------------------------------------
536  // with volume set by maximum required volume
537  //---------------------------------------
538  RefineableTetgenMesh<ELEMENT>* tmp_new_mesh_pt = 0;
539  if (solid_mesh_pt != 0)
540  {
541  /*throw OomphLibError("Solid RefineableTetgenMesh not done yet.",
542  OOMPH_CURRENT_FUNCTION,
543  OOMPH_EXCEPTION_LOCATION);*/
544  tmp_new_mesh_pt = new RefineableSolidTetgenMesh<ELEMENT>(
545  this->Outer_boundary_pt,
546  this->Internal_surface_pt,
547  max_size,
548  this->Time_stepper_pt,
549  this->Use_attributes,
550  this->Corner_elements_must_be_split);
551  }
552  else
553  {
554  tmp_new_mesh_pt = new RefineableTetgenMesh<ELEMENT>(
555  this->Outer_boundary_pt,
556  this->Internal_surface_pt,
557  max_size,
558  this->Time_stepper_pt,
559  this->Use_attributes,
560  this->Corner_elements_must_be_split);
561  }
562 
563 
564  //##################################################################
565  oomph_info
566  << "adapt: Time for building temp mesh : "
567  << TimingHelpers::timer() - t_start << " sec " << std::endl;
568  //##################################################################
569 
570  // tmp_new_mesh_pt->output("pre_mesh_nodes_snapped_0.dat");
571 
572  // //Move the nodes on the new boundary onto the
573  // //old curvilinear boundary
574  // //If the boundary is straight this will do precisely nothing
575  // //but will be somewhat inefficient
576  // {
577  // const unsigned n_boundary = this->nboundary();
578  // for(unsigned b=0;b<n_boundary;b++)
579  // {
580  // std::cout << "Snapping to boundary " << b << "\n";
581  // this->snap_nodes_onto_boundary(tmp_new_mesh_pt,b);
582  // }
583  // }
584 
585  // tmp_new_mesh_pt->output("mesh_nodes_snapped_0.dat");
586 
587 
588  // Get the tetgenio object associated with that mesh
589  tetgenio* tmp_new_tetgenio_pt = tmp_new_mesh_pt->tetgenio_pt();
590  RefineableTetgenMesh<ELEMENT>* new_mesh_pt = 0;
591 
592  // If the mesh is a solid mesh then do the mapping based on the
593  // Eulerian coordinates
594  bool use_eulerian_coords = false;
595  if (solid_mesh_pt != 0)
596  {
597  use_eulerian_coords = true;
598  }
599 
600 #ifdef OOMPH_HAS_CGAL
601 
602  // Make cgal-based bin
603  CGALSamplePointContainerParameters cgal_params(this);
604  if (use_eulerian_coords)
605  {
607  }
608  MeshAsGeomObject* mesh_geom_obj_pt = new MeshAsGeomObject(&cgal_params);
609 
610 #else
611 
612  std::ostringstream error_message;
613  error_message << "Non-CGAL-based target size transfer not implemented.\n";
614  throw OomphLibError(
615  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
616 
617  // Here's the relevant construction from the triangle mesh. Update...
618 
619  // // Make nonrefineable bin
620  // NonRefineableBinArrayParameters params(this);
621  // if (use_eulerian_coords)
622  // {
623  // params.enable_use_eulerian_coordinates_during_setup();
624  // }
625  // Vector<unsigned> bin_dim(2);
626  // bin_dim[0]=Nbin_x_for_area_transfer;
627  // bin_dim[1]=Nbin_y_for_area_transfer;
628  // params.dimensions_of_bin_array()=bin_dim;
629  // MeshAsGeomObject* mesh_geom_obj_pt=new MeshAsGeomObject(&params);
630 
631 #endif
632 
633  // Set up a map from pointer to element to its number
634  // in the mesh
635  std::map<GeneralisedElement*, unsigned> element_number;
636  unsigned nelem = this->nelement();
637  for (unsigned e = 0; e < nelem; e++)
638  {
639  element_number[this->element_pt(e)] = e;
640  }
641 
642  // Now start iterating to refine mesh recursively
643  //-----------------------------------------------
644  bool done = false;
645  unsigned iter = 0;
646  while (!done)
647  {
648  // Accept by default but overwrite if things go wrong below
649  done = true;
650 
651  // Loop over elements in new (tmp) mesh and visit all
652  // its integration points. Check where it's located in the bin
653  // structure of the current mesh and pass the target size
654  // to the new element
655  unsigned nelem = tmp_new_mesh_pt->nelement();
656 
657  // Store the target sizes for elements in the temporary
658  // mesh
659  Vector<double> new_transferred_target_size(nelem, 0.0);
660  for (unsigned e = 0; e < nelem; e++)
661  {
662  ELEMENT* el_pt =
663  dynamic_cast<ELEMENT*>(tmp_new_mesh_pt->element_pt(e));
664  unsigned nint = el_pt->integral_pt()->nweight();
665  for (unsigned ipt = 0; ipt < nint; ipt++)
666  {
667  // Get the coordinate of current point
668  Vector<double> s(3);
669  for (unsigned i = 0; i < 3; i++)
670  {
671  s[i] = el_pt->integral_pt()->knot(ipt, i);
672  }
673 
674  Vector<double> x(3);
675  el_pt->interpolated_x(s, x);
676 
677 #if OOMPH_HAS_CGAL
678 
679  // Try the five nearest sample points for Newton search
680  // then just settle on the nearest one
681  GeomObject* geom_obj_pt = 0;
682  unsigned max_sample_points = 5;
683  dynamic_cast<CGALSamplePointContainer*>(
684  mesh_geom_obj_pt->sample_point_container_pt())
685  ->limited_locate_zeta(x, max_sample_points, geom_obj_pt, s);
686 #ifdef PARANOID
687  if (geom_obj_pt == 0)
688  {
689  std::stringstream error_message;
690  error_message << "Limited locate zeta failed for zeta = [ "
691  << x[0] << " " << x[1] << " " << x[2]
692  << " ]. Makes no sense!\n";
693  throw OomphLibError(error_message.str(),
694  OOMPH_CURRENT_FUNCTION,
695  OOMPH_EXCEPTION_LOCATION);
696  }
697  else
698  {
699 #endif
700  FiniteElement* fe_pt = dynamic_cast<FiniteElement*>(geom_obj_pt);
701 #ifdef PARANOID
702  if (fe_pt == 0)
703  {
704  std::stringstream error_message;
705  error_message << "Cast to FE for GeomObject returned by "
706  "limited locate zeta failed for zeta = [ "
707  << x[0] << " " << x[1] << " " << x[2]
708  << " ]. Makes no sense!\n";
709  throw OomphLibError(error_message.str(),
710  OOMPH_CURRENT_FUNCTION,
711  OOMPH_EXCEPTION_LOCATION);
712  }
713  else
714  {
715 #endif
716  // What's the target size of the element that contains this
717  // point
718  double tg_size = target_size[element_number[fe_pt]];
719 
720  // Go for smallest target size over all integration
721  // points in new element
722  // to force "one level" of refinement (the one-level-ness
723  // is enforced below by limiting the actual reduction in
724  // size
725  if (new_transferred_target_size[e] != 0.0)
726  {
727  new_transferred_target_size[e] =
728  std::min(new_transferred_target_size[e], tg_size);
729  }
730  else
731  {
732  new_transferred_target_size[e] = tg_size;
733  }
734 #ifdef PARANOID
735  }
736  }
737 #endif
738 
739 // Non-CGAL
740 #else
741 
742  std::ostringstream error_message;
743  error_message
744  << "Non-CGAL-based target size transfer not implemented.\n";
745  throw OomphLibError(error_message.str(),
746  OOMPH_CURRENT_FUNCTION,
747  OOMPH_EXCEPTION_LOCATION);
748 
749  // Here's the relevant construction from the triangle mesh.
750  // Update...
751 
752  // // Find the bin that contains that point and its contents
753  // int bin_number=0;
754  // bin_array_pt->get_bin(x,bin_number);
755 
756  // // Did we find it?
757  // if (bin_number<0)
758  // {
759  // // Not even within bin boundaries... odd
760  // std::stringstream error_message;
761  // error_message
762  // << "Very odd -- we're looking for a point[ "
763  // << x[0] << " " << x[1] << " " << x[2] << " ] that's not even
764  // \n"
765  // << "located within the bin boundaries.\n";
766  // throw OomphLibError(error_message.str(),
767  // OOMPH_CURRENT_FUNCTION,
768  // OOMPH_EXCEPTION_LOCATION);
769  // } // if (bin_number<0)
770  // else
771  // {
772  // // Go for smallest target size of any element in this bin
773  // // to force "one level" of refinement (the one-level-ness
774  // // is enforced below by limiting the actual reduction in
775  // // size
776  // if (new_transferred_target_size[e]!=0)
777  // {
778  // std::min(new_transferred_target_size[e],
779  // bin_min_target_size[bin_number]);
780  // }
781  // else
782  // {
783  // new_transferred_target_size[e]=bin_min_target_size[bin_number];
784  // }
785 
786  // }
787 
788 #endif
789 
790  } // for (ipt<nint)
791 
792  } // for (e<nelem)
793 
794 
795  // do some output (keep it alive!)
796  const bool output_target_sizes = false;
797  if (output_target_sizes)
798  {
799  unsigned length = new_transferred_target_size.size();
800  for (unsigned u = 0; u < length; u++)
801  {
802  oomph_info << "Element" << u
803  << ",target size: " << new_transferred_target_size[u]
804  << std::endl;
805  }
806  }
807 
808  // Now copy into target size for temporary mesh but limit to
809  // the equivalent of one sub-division per iteration
810 #ifdef OOMPH_HAS_MPI
811  unsigned n_ele_need_refinement_iter = 0;
812 #endif
813 
814 
815  // Don't delete! Keep these around for debugging
816  // ofstream tmp_mesh_file;
817  // tmp_mesh_file.open("tmp_mesh_file.dat");
818  // tmp_new_mesh_pt->output(tmp_mesh_file);
819  // tmp_mesh_file.close();
820 
821  std::ofstream target_sizes_file;
822  char filename[100];
823  sprintf(filename, "target_sizes%i.dat", iter);
824  if (output_target_sizes)
825  {
826  target_sizes_file.open(filename);
827  }
828 
829  const unsigned nel_new = tmp_new_mesh_pt->nelement();
830  Vector<double> new_target_size(nel_new);
831  for (unsigned e = 0; e < nel_new; e++)
832  {
833  // The finite element
834  FiniteElement* f_ele_pt = tmp_new_mesh_pt->finite_element_pt(e);
835 
836  // Transferred target size
837  const double new_size = new_transferred_target_size[e];
838  if (new_size <= 0.0)
839  {
840  std::ostringstream error_stream;
841  error_stream
842  << "This shouldn't happen! Element whose centroid is at "
843  << (f_ele_pt->node_pt(0)->x(0) + f_ele_pt->node_pt(1)->x(0) +
844  f_ele_pt->node_pt(2)->x(0) + f_ele_pt->node_pt(3)->x(0)) /
845  4.0
846  << " "
847  << (f_ele_pt->node_pt(0)->x(1) + f_ele_pt->node_pt(1)->x(1) +
848  f_ele_pt->node_pt(2)->x(1) + f_ele_pt->node_pt(3)->x(1)) /
849  4.0
850  << " "
851  << (f_ele_pt->node_pt(0)->x(2) + f_ele_pt->node_pt(1)->x(2) +
852  f_ele_pt->node_pt(2)->x(2) + f_ele_pt->node_pt(3)->x(2)) /
853  4.0
854  << " "
855  << " has no target size assigned\n";
856  throw OomphLibError(error_stream.str(),
857  OOMPH_CURRENT_FUNCTION,
858  OOMPH_EXCEPTION_LOCATION);
859  }
860  else
861  {
862  // Limit target size to the equivalent of uniform refinement
863  // during this stage of the iteration
864  new_target_size[e] = new_size;
865  if (new_target_size[e] < f_ele_pt->size() / 4.0)
866  {
867  new_target_size[e] = f_ele_pt->size() / 4.0;
868 
869  // ALH: It seems that tetgen "enlarges" the volume constraint
870  // so this criterion can never be met unless dividing by 1.2
871  // as well. MH: Is this the reason why Andrew's version of
872  // adaptation never converges? Took it out.
873  // new_target_size[e] /= 1.2;
874 
875  // We'll need to give it another go later
876  done = false;
877  }
878 
879  // Don't delete! Keep around for debugging
880  if (output_target_sizes)
881  {
882  target_sizes_file << "ZONE N=4, E=1, F=FEPOINT, ET=TETRAHEDRON\n";
883  for (unsigned j = 0; j < 4; j++)
884  {
885  target_sizes_file << f_ele_pt->node_pt(j)->x(0) << " "
886  << f_ele_pt->node_pt(j)->x(1) << " "
887  << f_ele_pt->node_pt(j)->x(2) << " "
888  << new_size << " " << new_target_size[e]
889  << std::endl;
890  }
891  target_sizes_file << "1 2 3 4\n"; // connectivity
892 
893  // Keep around; just doc at centroid
894  /* target_sizes_file */
895  /* << (f_ele_pt->node_pt(0)->x(0)+ */
896  /* f_ele_pt->node_pt(1)->x(0)+ */
897  /* f_ele_pt->node_pt(2)->x(0)+ */
898  /* f_ele_pt->node_pt(3)->x(0))/4.0 << " " */
899  /* << (f_ele_pt->node_pt(0)->x(1)+ */
900  /* f_ele_pt->node_pt(1)->x(1)+ */
901  /* f_ele_pt->node_pt(2)->x(1)+ */
902  /* f_ele_pt->node_pt(3)->x(1))/4.0 << " " */
903  /* << (f_ele_pt->node_pt(0)->x(2)+ */
904  /* f_ele_pt->node_pt(1)->x(2)+ */
905  /* f_ele_pt->node_pt(2)->x(2)+ */
906  /* f_ele_pt->node_pt(3)->x(2))/4.0 << " " */
907  /* << new_size << " " */
908  /* << new_target_size[e] << std::endl; */
909  }
910 
911 #ifdef OOMPH_HAS_MPI
912  // Keep track of the elements that require (un)refinement
913  n_ele_need_refinement_iter++;
914 #endif
915 
916  } // else if (new_size <= 0.0)
917 
918  } // for (e < nel_new)
919 
920  // Don't delete! Keep around for debugging
921  if (output_target_sizes)
922  {
923  target_sizes_file.close();
924  }
925 
926  if (done)
927  {
928  oomph_info
929  << "All size adjustments accommodated by max. permitted size"
930  << " reduction during iter " << iter << "\n";
931  }
932  else
933  {
934  oomph_info << "NOT all size adjustments accommodated by max. "
935  << "permitted size reduction during iter " << iter
936  << "\n";
937  }
938 
939 
940  oomph_info
941  << "\n\n\n==================================================\n"
942  << "==================================================\n"
943  << "==================================================\n"
944  << "==================================================\n"
945  << "\n\n\n";
946 
947  //##################################################################
948  oomph_info
949  << "adapt: Time for new_target_size[.] : "
950  << TimingHelpers::timer() - t_start << " sec " << std::endl;
951  //##################################################################
952 
953 
954  // Now create the new mesh from TriangulateIO structure
955  //-----------------------------------------------------
956  // associated with uniform background mesh and the
957  //------------------------------------------------
958  // associated target element sizes.
959  //---------------------------------
960 
961  //###################################
962  t_start = TimingHelpers::timer();
963  //###################################
964 
965  // Solid mesh?
966  if (solid_mesh_pt != 0)
967  {
968  /*std::ostringstream error_message;
969  error_message
970  << "RefineableSolidTetgenMesh not implemented yet.\n";
971  throw OomphLibError(error_message.str(),
972  OOMPH_CURRENT_FUNCTION,
973  OOMPH_EXCEPTION_LOCATION);*/
974  new_mesh_pt =
975  new RefineableSolidTetgenMesh<ELEMENT>(new_target_size,
976  tmp_new_tetgenio_pt,
977  this->Outer_boundary_pt,
978  this->Internal_surface_pt,
979  this->Time_stepper_pt,
980  this->Use_attributes);
981  }
982  // No solid mesh
983  else
984  {
985  new_mesh_pt =
986  new RefineableTetgenMesh<ELEMENT>(new_target_size,
987  tmp_new_tetgenio_pt,
988  this->Outer_boundary_pt,
989  this->Internal_surface_pt,
990  this->Time_stepper_pt,
991  this->Use_attributes);
992  }
993 
994  //##################################################################
995  oomph_info
996  << "adapt: Time for new_mesh_pt : "
997  << TimingHelpers::timer() - t_start << " sec " << std::endl;
998  //##################################################################
999 
1000 
1001  // Not done: get ready for another iteration
1002  iter++;
1003 
1004 
1005  // Check if the new mesh actually differs from the old one
1006  // If not, we're done.
1007  if (!done)
1008  {
1009  unsigned nnod = tmp_new_mesh_pt->nnode();
1010  if (nnod == new_mesh_pt->nnode())
1011  {
1012  unsigned nel = tmp_new_mesh_pt->nelement();
1013  if (nel == new_mesh_pt->nelement())
1014  {
1015  bool nodes_identical = true;
1016  for (unsigned j = 0; j < nnod; j++)
1017  {
1018  bool coords_identical = true;
1019  for (unsigned i = 0; i < 3; i++)
1020  {
1021  if (new_mesh_pt->node_pt(j)->x(i) !=
1022  tmp_new_mesh_pt->node_pt(j)->x(i))
1023  {
1024  coords_identical = false;
1025  }
1026  }
1027  if (!coords_identical)
1028  {
1029  nodes_identical = false;
1030  break;
1031  }
1032  }
1033  if (nodes_identical)
1034  {
1035  done = true;
1036  }
1037  }
1038  }
1039  }
1040 
1041  // Delete the temporary mesh
1042  delete tmp_new_mesh_pt;
1043 
1044  // Now transfer over the pointers
1045  if (!done)
1046  {
1047  tmp_new_mesh_pt = new_mesh_pt;
1048  tmp_new_tetgenio_pt = new_mesh_pt->tetgenio_pt();
1049  }
1050 
1051  } // end of iteration
1052 
1053  // Move the nodes on the new boundary onto the
1054  // old curvilinear boundary
1055  // If the boundary is straight this will do precisely nothing
1056  // but will be somewhat inefficient
1057  /*{
1058  const unsigned n_boundary = this->nboundary();
1059  for(unsigned b=0;b<n_boundary;b++)
1060  {
1061  this->snap_nodes_onto_boundary(new_mesh_pt,b);
1062  }
1063  }*/
1064 
1065 
1066  // Check that the projection step is not disabled
1067  if (!Projection_is_disabled)
1068  {
1069  //###################################
1070  t_start = TimingHelpers::timer();
1071  //###################################
1072 
1073  // Project current solution onto new mesh
1074  //---------------------------------------
1075  ProjectionProblem<ELEMENT>* project_problem_pt =
1077  project_problem_pt->mesh_pt() = new_mesh_pt;
1078  project_problem_pt->project(this);
1079  delete project_problem_pt;
1080 
1081  //##################################################################
1082  oomph_info
1083  << "adapt: Time for project soln onto new mesh : "
1084  << TimingHelpers::timer() - t_start << " sec " << std::endl;
1085  //##################################################################
1086  }
1087 
1088  //###################################
1089  t_start = TimingHelpers::timer();
1090  //###################################
1091 
1092  // this->output("pre_proj",5);
1093  // new_mesh_pt->output("post_proj.dat",5);
1094 
1095  // Flush the old mesh
1096  unsigned nnod = nnode();
1097  for (unsigned j = nnod; j > 0; j--)
1098  {
1099  delete Node_pt[j - 1];
1100  Node_pt[j - 1] = 0;
1101  }
1102  unsigned nel = nelement();
1103  for (unsigned e = nel; e > 0; e--)
1104  {
1105  delete Element_pt[e - 1];
1106  Element_pt[e - 1] = 0;
1107  }
1108 
1109  // Now copy back to current mesh
1110  //------------------------------
1111  nnod = new_mesh_pt->nnode();
1112  Node_pt.resize(nnod);
1113  nel = new_mesh_pt->nelement();
1114  Element_pt.resize(nel);
1115  for (unsigned j = 0; j < nnod; j++)
1116  {
1117  Node_pt[j] = new_mesh_pt->node_pt(j);
1118  }
1119  for (unsigned e = 0; e < nel; e++)
1120  {
1121  Element_pt[e] = new_mesh_pt->element_pt(e);
1122  }
1123 
1124  // Copy the boundary schemes
1125  unsigned nbound = new_mesh_pt->nboundary();
1126  Boundary_element_pt.resize(nbound);
1127  Face_index_at_boundary.resize(nbound);
1128  Boundary_node_pt.resize(nbound);
1129  for (unsigned b = 0; b < nbound; b++)
1130  {
1131  unsigned nel = new_mesh_pt->nboundary_element(b);
1132  Boundary_element_pt[b].resize(nel);
1133  Face_index_at_boundary[b].resize(nel);
1134  for (unsigned e = 0; e < nel; e++)
1135  {
1136  Boundary_element_pt[b][e] = new_mesh_pt->boundary_element_pt(b, e);
1137  Face_index_at_boundary[b][e] =
1138  new_mesh_pt->face_index_at_boundary(b, e);
1139  }
1140  unsigned nnod = new_mesh_pt->nboundary_node(b);
1141  Boundary_node_pt[b].resize(nnod);
1142  for (unsigned j = 0; j < nnod; j++)
1143  {
1144  Boundary_node_pt[b][j] = new_mesh_pt->boundary_node_pt(b, j);
1145  }
1146  }
1147 
1148  // Also copy over the new boundary and region information
1149  unsigned n_region = new_mesh_pt->nregion();
1150 
1151  // Only bother if we have regions
1152  if (n_region > 1)
1153  {
1154  // Deal with the region information first
1155  this->Region_element_pt.resize(n_region);
1156  this->Region_attribute.resize(n_region);
1157  for (unsigned i = 0; i < n_region; i++)
1158  {
1159  // Copy across region attributes (region ids!)
1160  this->Region_attribute[i] = new_mesh_pt->region_attribute(i);
1161 
1162  // Find the number of elements in the region
1163  unsigned r = this->Region_attribute[i];
1164  unsigned n_region_element = new_mesh_pt->nregion_element(r);
1165  this->Region_element_pt[i].resize(n_region_element);
1166  for (unsigned e = 0; e < n_region_element; e++)
1167  {
1168  this->Region_element_pt[i][e] =
1169  new_mesh_pt->region_element_pt(r, e);
1170  }
1171  }
1172 
1173  // Now the boundary region information
1174  this->Boundary_region_element_pt.resize(nbound);
1175  this->Face_index_region_at_boundary.resize(nbound);
1176 
1177  // Now loop over the boundaries
1178  for (unsigned b = 0; b < nbound; ++b)
1179  {
1180  // Loop over the regions
1181  for (unsigned i = 0; i < n_region; ++i)
1182  {
1183  unsigned r = this->Region_attribute[i];
1184 
1185  unsigned n_boundary_el_in_region =
1186  new_mesh_pt->nboundary_element_in_region(b, r);
1187  if (n_boundary_el_in_region > 0)
1188  {
1189  // Allocate storage in the map
1190  this->Boundary_region_element_pt[b][r].resize(
1191  n_boundary_el_in_region);
1192  this->Face_index_region_at_boundary[b][r].resize(
1193  n_boundary_el_in_region);
1194 
1195  // Copy over the information
1196  for (unsigned e = 0; e < n_boundary_el_in_region; ++e)
1197  {
1198  this->Boundary_region_element_pt[b][r][e] =
1199  new_mesh_pt->boundary_element_in_region_pt(b, r, e);
1200  this->Face_index_region_at_boundary[b][r][e] =
1201  new_mesh_pt->face_index_at_boundary_in_region(b, r, e);
1202  }
1203  }
1204  }
1205  } // End of loop over boundaries
1206 
1207  } // End of case when more than one region
1208 
1209  // Copy TriangulateIO representation
1210  this->set_deep_copy_tetgenio_pt(new_mesh_pt->tetgenio_pt());
1211 
1212  // Flush the mesh
1213  new_mesh_pt->flush_element_and_node_storage();
1214 
1215  // Delete the mesh and the problem
1216  delete new_mesh_pt;
1217 
1218 
1219  //##################################################################
1220  oomph_info
1221  << "adapt: Time for moving nodes etc. to actual mesh : "
1222  << TimingHelpers::timer() - t_start << " sec " << std::endl;
1223  //##################################################################
1224 
1225  // Solid mesh?
1226  if (solid_mesh_pt != 0)
1227  {
1228  // Warning
1229  std::stringstream error_message;
1230  error_message
1231  << "Lagrangian coordinates are currently not projected but are\n"
1232  << "are re-set during adaptation. This is not appropriate for\n"
1233  << "real solid mechanics problems!\n";
1234  OomphLibWarning(error_message.str(),
1235  OOMPH_CURRENT_FUNCTION,
1236  OOMPH_EXCEPTION_LOCATION);
1237 
1238  // Reset Lagrangian coordinates
1239  dynamic_cast<SolidMesh*>(this)->set_lagrangian_nodal_coordinates();
1240  }
1241 
1242  double max_area;
1243  double min_area;
1244  this->max_and_min_element_size(max_area, min_area);
1245  oomph_info << "Max/min element size in adapted mesh: " << max_area << " "
1246  << min_area << std::endl;
1247  }
1248  else
1249  {
1250  oomph_info << "Not enough benefit in adaptation.\n";
1251  Nrefined = 0;
1252  Nunrefined = 0;
1253  }
1254  }
1255 
1256 } // namespace oomph
1257 
1258 #endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
Class that is used to create FaceElement from bulk elements and to provide these FaceElement with a g...
void set_boundary_number_in_bulk_mesh(const unsigned &b)
Set function for the boundary number in bulk mesh.
Definition: elements.h:4482
A general Finite Element class.
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition: elements.cc:4290
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void position(const Vector< double > &zeta, Vector< double > &r) const
Return the position as a function of the intrinsic coordinate zeta. This provides an (expensive!...
SamplePointContainer * sample_point_container_pt() const
Pointer to the sample point container.
A general mesh class.
Definition: mesh.h:67
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
void flush_element_and_node_storage()
Flush storage for elements and nodes by emptying the vectors that store the pointers to them....
Definition: mesh.h:407
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
Definition: mesh.h:896
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
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_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b....
Definition: nodes.cc:2379
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
Projection problem. This is created during the adaptation of unstructured meshes and it is assumed th...
Definition: projection.h:695
void project(Mesh *base_mesh_pt, const bool &dont_project_positions=false)
Project from base into the problem's own mesh.
Definition: projection.h:748
////////////////////////////////////////////////////////////////////////// //////////////////////////...
void surface_remesh_for_inner_hole_boundaries()
Generate a new faceted representation of the inner hole boundaries.
void update_faceted_surface_using_face_mesh(TetMeshFacetedSurface *&faceted_surface_pt)
Helper function that updates the input faceted surface by using the flattened elements from FaceMesh(...
void snap_nodes_onto_boundary(RefineableTetgenMesh< ELEMENT > *&new_mesh_pt, const unsigned &b)
Snap the boundary nodes onto any curvilinear boundaries.
void adapt(const Vector< double > &elem_error)
Adapt mesh, based on elemental error provided.
void enable_use_eulerian_coordinates_during_setup()
Enable use of eulerian coordinates (via interpolated_x) during setup (otherwise use interpolated_zeta...
General SolidMesh class.
Definition: mesh.h:2562
General TElement class.
Definition: Telements.h:1208
int face_index_at_boundary_in_region(const unsigned &b, const unsigned &r, const unsigned &e) const
Return face index of the e-th element adjacent to boundary b in region r.
Definition: tet_mesh.h:834
double region_attribute(const unsigned &i)
Return the i-th region attribute (here only used as the (assumed to be unsigned) region id.
Definition: tet_mesh.h:906
unsigned nregion()
Return the number of regions specified by attributes.
Definition: tet_mesh.h:860
unsigned nboundary_element_in_region(const unsigned &b, const unsigned &r) const
Return the number of elements adjacent to boundary b in region r.
Definition: tet_mesh.h:797
FiniteElement * boundary_element_in_region_pt(const unsigned &b, const unsigned &r, const unsigned &e) const
Return pointer to the e-th element adjacent to boundary b in region r.
Definition: tet_mesh.h:816
FiniteElement * region_element_pt(const unsigned &r, const unsigned &e)
Return the e-th element in the r-th region.
Definition: tet_mesh.h:912
unsigned nregion_element(const unsigned &r)
Return the number of elements in region r.
Definition: tet_mesh.h:866
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: tet_mesh.h:634
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:538
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:306
unsigned nfacet() const
Number of facets.
Definition: tet_mesh.h:325
unsigned one_based_facet_boundary_id(const unsigned &j) const
One-based boundary id of j-th facet.
Definition: tet_mesh.h:331
tetgenio *& tetgenio_pt()
Access to the triangulateio representation of the mesh.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...