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-2022 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
42namespace 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!
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!
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 //##################################################################
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 //##################################################################
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 {
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
941 << "\n\n\n==================================================\n"
942 << "==================================================\n"
943 << "==================================================\n"
944 << "==================================================\n"
945 << "\n\n\n";
946
947 //##################################################################
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 //##################################################################
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 //##################################################################
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 //##################################################################
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
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition: elements.cc:4290
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
SamplePointContainer * sample_point_container_pt() const
Pointer to the sample point container.
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!...
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
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
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
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
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
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 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
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.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
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...