tet_mesh.h
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// Common base class for all Tet Meshes
27#ifndef OOMPH_TET_MESH_HEADER
28#define OOMPH_TET_MESH_HEADER
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35// oomph-lib includes
36#include "Vector.h"
37#include "nodes.h"
38#include "matrices.h"
39#include "mesh.h"
41
42namespace oomph
43{
44 //=======================================================================
45 /// Vertex for Tet mesh generation. Can lie on multiple boundaries
46 /// (identified via one-based enumeration!) and can have intrinisic
47 /// coordinates in a DiskLikeGeomObjectWithBoundaries.
48 //=======================================================================
50 {
51 public:
52 /// Only friends can set boundary ID -- the facet is my only friend!
53 friend class TetMeshFacet;
54
55 /// Constructor: Pass coordinates (length 3!)
57 {
58#ifdef PARANOID
59 if (X.size() != 3)
60 {
61 std::ostringstream error_stream;
62 error_stream << "TetMeshVertex should only be used in 3D!\n"
63 << "Your Vector of coordinates, contains data for "
64 << x.size() << "-dimensional coordinates." << std::endl;
65 throw OomphLibError(
66 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
67 }
68#endif
69 }
70
71 // Constructor: Take coordinates from a node
72 TetMeshVertex(Node* const& node_pt)
73 {
74 const unsigned n_dim = node_pt->ndim();
75#ifdef PARANOID
76 if (n_dim != 3)
77 {
78 std::ostringstream error_stream;
79 error_stream << "TetMeshVertex should only be used in 3D!\n"
80 << "Your Node contains data for " << n_dim
81 << "-dimensional coordinates." << std::endl;
82 throw OomphLibError(
83 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
84 }
85#endif
86
87 // Copy over the positions from the node
88 X.resize(n_dim);
89 for (unsigned i = 0; i < n_dim; ++i)
90 {
91 X[i] = node_pt->x(i);
92 }
93 }
94
95
96 /// Set intrinisic coordinates in GeomObject
98 {
99#ifdef PARANOID
100 if (zeta.size() != 2)
101 {
102 std::ostringstream error_stream;
103 error_stream
104 << "TetMeshVertex should only be used in 3D!\n"
105 << "Your Vector of intrinisic coordinates, contains data for "
106 << zeta.size() << "-dimensional coordinates but should be 2!"
107 << std::endl;
108 throw OomphLibError(
109 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
110 }
111#endif
112
113 Zeta_in_geom_object = zeta;
114 }
115
116 /// Get intrinisic coordinates in GeomObject (returns zero sized
117 /// vector if no such coordinates have been specified)
119 {
120 return Zeta_in_geom_object;
121 }
122
123 /// i-th coordinate
124 double x(const unsigned& i) const
125 {
126 return X[i];
127 }
128
129 /// First (of possibly multiple) one-based boundary id
130 unsigned one_based_boundary_id() const
131 {
132 if (One_based_boundary_id.size() == 0)
133 {
134 return 0;
135 }
136 return *(One_based_boundary_id.begin());
137 }
138
139 private:
140 /// Set of (one-based!) boundary IDs this vertex lives on
141 void set_one_based_boundary_id(const unsigned& id)
142 {
143 One_based_boundary_id.insert(id);
144 }
145
146 /// Coordinate vector
148
149 /// Set of (one-based!) boundary IDs this vertex lives on
150 std::set<unsigned> One_based_boundary_id;
151
152 /// Intrinisic coordinates in GeomObject (zero sized if there
153 /// isn't one.
155 };
156
157 /// /////////////////////////////////////////////////////////////////////
158 /// /////////////////////////////////////////////////////////////////////
159 /// /////////////////////////////////////////////////////////////////////
160
161
162 //=======================================================================
163 /// Facet for Tet mesh generation. Can lie on boundary
164 /// (identified via one-based enumeration!) and can have
165 /// GeomObject associated with those boundaries.
166 //=======================================================================
168 {
169 public:
170 /// Constructor: Specify number of vertices
171 TetMeshFacet(const unsigned& nvertex)
172 : One_based_boundary_id(0), // Initialision implies not on any boundary
174 0) // Initialisation implies
175 // not embedded in any region
176 {
177 Vertex_pt.resize(nvertex, 0);
178 }
179
180 /// Number of vertices
181 unsigned nvertex() const
182 {
183 return Vertex_pt.size();
184 }
185
186 /// Pointer to j-th vertex (const)
187 TetMeshVertex* vertex_pt(const unsigned& j) const
188 {
189 return Vertex_pt[j];
190 }
191
192 /// Set pointer to j-th vertex and pass one-based boundary id that
193 /// may already have been set for this facet.
194 void set_vertex_pt(const unsigned& j, TetMeshVertex* vertex_pt)
195 {
196 Vertex_pt[j] = vertex_pt;
197 // If not set yet, this is harmless since it simply over-writes
198 // the dummy value in vertex
199 TetMeshVertex* v_pt = Vertex_pt[j];
200 if (v_pt != 0)
201 {
203 }
204 }
205
206 /// Constant access to (one-based!) boundary IDs this facet lives on
207 unsigned one_based_boundary_id() const
208 {
210 }
211
212 /// Set (one-based!) boundary IDs this facet lives on. Passed to any
213 /// existing vertices and to any future ones
214 void set_one_based_boundary_id(const unsigned& one_based_id)
215 {
216 One_based_boundary_id = one_based_id;
217 unsigned nv = Vertex_pt.size();
218 for (unsigned j = 0; j < nv; j++)
219 {
220 TetMeshVertex* v_pt = Vertex_pt[j];
221 if (v_pt != 0)
222 {
223 v_pt->set_one_based_boundary_id(one_based_id);
224 }
225 }
226 }
227
228 /// Set (one-based!) region ID this facet is adjacent to.
229 /// Specification of zero argument is harmless as it indicates that
230 /// that facet is not adjacent to any region.
231 void set_one_based_adjacent_region_id(const unsigned& one_based_id)
232 {
233 One_based_adjacent_region_id.insert(one_based_id);
234 }
235
236 /// Return set of (one-based!) region IDs this facet is adjacent to
237 std::set<unsigned> one_based_adjacent_region_id() const
238 {
240 }
241
242 /// Boolean indicating that facet is embedded in a specified region
244 {
246 }
247
248 /// Facet is to be embedded in specified one-based region
250 const unsigned& one_based_region_id)
251 {
253 }
254
255 /// Which (one-based) region is facet embedded in (if zero, it's not
256 /// embedded in any region)
258 {
260 }
261
262 /// Output
263 void output(std::ostream& outfile) const
264 {
265 unsigned n = Vertex_pt.size();
266 outfile << "ZONE I=" << n + 1 << std::endl;
267 for (unsigned j = 0; j < n; j++)
268 {
269 outfile << Vertex_pt[j]->x(0) << " " << Vertex_pt[j]->x(1) << " "
270 << Vertex_pt[j]->x(2) << " " << One_based_boundary_id
271 << std::endl;
272 }
273 outfile << Vertex_pt[0]->x(0) << " " << Vertex_pt[0]->x(1) << " "
274 << Vertex_pt[0]->x(2) << " " << One_based_boundary_id
275 << std::endl;
276 }
277
278 private:
279 /// Pointer to vertices
281
282 /// (One-based!) boundary IDs this facet lives on
284
285 /// Set of one-based adjacent region ids; no adjacent region if
286 /// it's zero.
288
289
290 /// Facet is to be embedded in specified one-based region.
291 /// Defaults to zero, indicating that its not embedded.
293 };
294
295
296 /// /////////////////////////////////////////////////////////////////////
297 /// /////////////////////////////////////////////////////////////////////
298 /// /////////////////////////////////////////////////////////////////////
299
300
301 //========================================================================
302 /// Base class for tet mesh boundary defined by polygonal
303 /// planar facets
304 //========================================================================
306 {
307 public:
308 /// Constructor:
312 {
313 }
314
315 /// Empty destructor
317
318 /// Number of vertices
319 unsigned nvertex() const
320 {
321 return Vertex_pt.size();
322 }
323
324 /// Number of facets
325 unsigned nfacet() const
326 {
327 return Facet_pt.size();
328 }
329
330 /// One-based boundary id of j-th facet
331 unsigned one_based_facet_boundary_id(const unsigned& j) const
332 {
333 return Facet_pt[j]->one_based_boundary_id();
334 }
335
336 /// First (of possibly multiple) one-based boundary id of j-th vertex
337 unsigned one_based_vertex_boundary_id(const unsigned& j) const
338 {
339 return Vertex_pt[j]->one_based_boundary_id();
340 }
341
342 /// i-th coordinate of j-th vertex
343 double vertex_coordinate(const unsigned& j, const unsigned& i) const
344 {
345 return Vertex_pt[j]->x(i);
346 }
347
348 /// Number of vertices defining the j-th facet
349 unsigned nvertex_on_facet(const unsigned& j) const
350 {
351 return Facet_pt[j]->nvertex();
352 }
353
354 /// Test whether boundary can be split in tetgen
356 {
358 }
359
360 /// Test whether boundaries can be split in tetgen
362 {
364 }
365
366 /// Test whether boundaries can be split in tetgen
368 {
370 }
371
372 /// Pointer to j-th facet
373 TetMeshFacet* facet_pt(const unsigned& j) const
374 {
375 return Facet_pt[j];
376 }
377
378 /// Pointer to j-th vertex
379 TetMeshVertex* vertex_pt(const unsigned& j) const
380 {
381 return Vertex_pt[j];
382 }
383
384 /// Access to GeomObject with boundaries associated with this
385 /// surface (Null if there isn't one!)
387 {
389 }
390
391 /// Output
392 void output(std::ostream& outfile) const
393 {
394 unsigned n = Facet_pt.size();
395 for (unsigned j = 0; j < n; j++)
396 {
397 Facet_pt[j]->output(outfile);
398 }
399 }
400
401
402 /// Output
403 void output(const std::string& filename) const
404 {
405 std::ofstream outfile;
406 outfile.open(filename.c_str());
407 output(outfile);
408 outfile.close();
409 }
410
411 /// Virtual function that specifies the variation of the
412 /// zeta coordinates in the GeomObject along the
413 /// boundary connecting vertices 0 and 1 in
414 /// facet facet_id. Default implementation: Linear interpolation
415 /// between edges. zeta_boundary=0.0: we're on vertex 0;
416 /// zeta_boundary=1.0: we're on vertex 1.
417 virtual void boundary_zeta01(const unsigned& facet_id,
418 const double& zeta_boundary,
419 Vector<double>& zeta)
420 {
421 Vector<Vector<double>> zeta_vertex(2);
422 zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(0)->zeta_in_geom_object();
423 zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(1)->zeta_in_geom_object();
424 zeta[0] = zeta_vertex[0][0] +
425 (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
426 zeta[1] = zeta_vertex[0][1] +
427 (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
428 }
429
430 /// Virtual function that specifies the variation of the
431 /// zeta coordinates in the GeomObject along the
432 /// boundary connecting vertices 1 and 2 in
433 /// facet facet_id. Default implementation: Linear interpolation
434 /// between edges. zeta_boundary=0.0: we're on vertex 1;
435 /// zeta_boundary=1.0: we're on vertex 2.
436 virtual void boundary_zeta12(const unsigned& facet_id,
437 const double& zeta_boundary,
438 Vector<double>& zeta)
439 {
440 Vector<Vector<double>> zeta_vertex(2);
441 zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(1)->zeta_in_geom_object();
442 zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(2)->zeta_in_geom_object();
443 zeta[0] = zeta_vertex[0][0] +
444 (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
445 zeta[1] = zeta_vertex[0][1] +
446 (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
447 }
448
449 /// Virtual function that specifies the variation of the
450 /// zeta coordinates in the GeomObject along the
451 /// boundary connecting vertices 2 and 0 in
452 /// facet facet_id. Default implementation: Linear interpolation
453 /// between edges. zeta_boundary=0.0: we're on vertex 2;
454 /// zeta_boundary=1.0: we're on vertex 0.
455 virtual void boundary_zeta20(const unsigned& facet_id,
456 const double& zeta_boundary,
457 Vector<double>& zeta)
458 {
459 Vector<Vector<double>> zeta_vertex(2);
460 zeta_vertex[0] = Facet_pt[facet_id]->vertex_pt(2)->zeta_in_geom_object();
461 zeta_vertex[1] = Facet_pt[facet_id]->vertex_pt(0)->zeta_in_geom_object();
462 zeta[0] = zeta_vertex[0][0] +
463 (zeta_vertex[1][0] - zeta_vertex[0][0]) * zeta_boundary;
464 zeta[1] = zeta_vertex[0][1] +
465 (zeta_vertex[1][1] - zeta_vertex[0][1]) * zeta_boundary;
466 }
467
468
469 /// Facet connectivity: vertex_index[j] is the index of the
470 /// j-th vertex (in the Vertex_pt vector) in facet f. Bit of an obscure
471 /// functionality that's only needed for setup tetgen_io.
473 {
474 if (Facet_vertex_index_in_tetgen.size() != nfacet())
475 {
477 }
479 }
480
481 protected:
482 /// Vector pointers to vertices
484
485 /// Vector of pointers to facets
487
488 /// Boolean to indicate whether extra vertices can be added
489 /// on the boundary in tetgen
491
492 /// Facet connectivity: Facet_vertex_index[f][j] is the index of the
493 /// j-th vertex (in the Vertex_pt vector) in facet f.
495
496 /// GeomObject with boundaries associated with this surface
498
499
500 private:
501 /// Setup facet connectivity for tetgen
503 {
504 unsigned nv_overall = Vertex_pt.size();
505 unsigned nf = nfacet();
507 for (unsigned f = 0; f < nf; f++)
508 {
509 unsigned nv = Facet_pt[f]->nvertex();
510 Facet_vertex_index_in_tetgen[f].resize(nv);
511 for (unsigned v = 0; v < nv; v++)
512 {
513 TetMeshVertex* my_vertex_pt = Facet_pt[f]->vertex_pt(v);
514 for (unsigned j = 0; j < nv_overall; j++)
515 {
516 if (my_vertex_pt == Vertex_pt[j])
517 {
519 break;
520 }
521 }
522 }
523 }
524 }
525 };
526
527
528 /// /////////////////////////////////////////////////////////////////////
529 /// /////////////////////////////////////////////////////////////////////
530 /// /////////////////////////////////////////////////////////////////////
531
532
533 //========================================================================
534 /// Base class for closed tet mesh boundary bounded by polygonal
535 /// planar facets
536 //========================================================================
538 {
539 public:
540 /// Constructor:
543 {
544 }
545
546 /// Empty destructor
548
549 /// Declare closed surface to represent hole for gmsh
551 {
553 }
554
555 /// Declare closed surface NOT to represent hole for gmsh
557 {
559 }
560
561 /// Does closed surface represent hole for gmsh?
563 {
565 }
566
567 /// i=th coordinate of the j-th internal point for tetgen
568 const double& internal_point_for_tetgen(const unsigned& j,
569 const unsigned& i) const
570 {
571 return (Internal_point_for_tetgen[j].first)[i];
572 }
573
574 /// Specify coordinate of hole for tetgen
575 void set_hole_for_tetgen(const Vector<double>& hole_point)
576 {
577 Internal_point_for_tetgen.push_back(std::make_pair(hole_point, -1));
578 }
579
580 /// Specify a region; pass (zero-based) region ID and coordinate
581 /// of point in region for tetgen
582 void set_region_for_tetgen(const unsigned& region_id,
583 const Vector<double>& region_point)
584 {
586 std::make_pair(region_point, region_id));
587 }
588
589 /// Number of internal points (identifying either regions or holes)
590 /// for tetgen
592 {
593 return Internal_point_for_tetgen.size();
594 }
595
596 /// Return the (zero-based) region ID of j-th internal point for
597 /// tetgen. Negative if it's actually a hole.
598 const int& region_id_for_tetgen(const unsigned& j) const
599 {
600 return Internal_point_for_tetgen[j].second;
601 }
602
603
604 /// Is j-th internal point for tetgen associated with a hole?
606 {
607 return (Internal_point_for_tetgen[j].second < 0);
608 }
609
610 /// Is j-th internal point for tetgen associated with a region?
612 {
613 return (Internal_point_for_tetgen[j].second >= 0);
614 }
615
616
617 private:
618 /// Storage for internal points for tetgen. Stores pair of:
619 /// -- Vector containing coordinates of internal point
620 /// -- region ID (negative if it's a hole)
622
623 /// Does closed surface represent hole for gmsh?
625 };
626
627 /// ////////////////////////////////////////////////////////////////////
628 /// ////////////////////////////////////////////////////////////////////
629 /// ////////////////////////////////////////////////////////////////////
630
631
634 {
635 public:
636 // Constructor, which requires node, connectivity and boundary information
638 Vector<Node*> const& vertex_node_pt,
639 Vector<Vector<unsigned>> const& facet_connectivity,
640 Vector<unsigned> const& facet_boundary_id);
641
642 // Destructor
644 };
645
646
647 /// ////////////////////////////////////////////////////////////////////
648 /// ////////////////////////////////////////////////////////////////////
649 /// ////////////////////////////////////////////////////////////////////
650
651
652 /// ////////////////////////////////////////////////////////////////////
653 /// ////////////////////////////////////////////////////////////////////
654 /// ////////////////////////////////////////////////////////////////////
655
656
657 //================================================================
658 /// Base class for tet meshes (meshes made of 3D tet elements).
659 //================================================================
660 class TetMeshBase : public virtual Mesh
661 {
662 public:
663 /// Constructor
665
666 /// Broken copy constructor
667 TetMeshBase(const TetMeshBase& node) = delete;
668
669 /// Broken assignment operator
670 void operator=(const TetMeshBase&) = delete;
671
672 /// Destructor (empty)
673 virtual ~TetMeshBase() {}
674
675 /// Global static data that specifies the permitted
676 /// error in the setup of the boundary coordinates
678
679 /// Assess mesh quality: Ratio of max. edge length to min. height,
680 /// so if it's very large it's BAAAAAD.
681 void assess_mesh_quality(std::ofstream& some_file);
682
683 /// Setup boundary coordinate on boundary b which is
684 /// assumed to be planar. Boundary coordinates are the
685 /// x-y coordinates in the plane of that boundary, with the
686 /// x-axis along the line from the (lexicographically)
687 /// "lower left" to the "upper right" node. The y axis
688 /// is obtained by taking the cross-product of the positive
689 /// x direction with the outer unit normal computed by
690 /// the face elements (or its negative if switch_normal is set
691 /// to true). Doc faces in output file (if it's open).
692 ///
693 /// Note 1: Setup of boundary coordinates is not done if the boundary in
694 /// question turns out to be nonplanar.
695 ///
696 /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
697 /// also
698 /// store the boundary coordinates of its vertices. They are needed
699 /// to interpolated intrinsic coordinates of an associated
700 /// GeomObject (if any) into the interior.
701 template<class ELEMENT>
702 void setup_boundary_coordinates(const unsigned& b)
703 {
704 // Dummy file
705 std::ofstream some_file;
706
707 // Don't switch the normal
708 bool switch_normal = false;
709 this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
710 }
711
712 /// Setup boundary coordinate on boundary b which is
713 /// assumed to be planar. Boundary coordinates are the
714 /// x-y coordinates in the plane of that boundary, with the
715 /// x-axis along the line from the (lexicographically)
716 /// "lower left" to the "upper right" node. The y axis
717 /// is obtained by taking the cross-product of the positive
718 /// x direction with the outer unit normal computed by
719 /// the face elements (or its negative if switch_normal is set
720 /// to true). Doc faces in output file (if it's open).
721 ///
722 /// Note 1: Setup of boundary coordinates is not done if the boundary in
723 /// question turns out to be nonplanar.
724 ///
725 /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
726 /// also
727 /// store the boundary coordinates of its vertices. They are needed
728 /// to interpolated intrinsic coordinates of an associated
729 /// GeomObject (if any) into the interior.
730 /// Final boolean argument allows switching of the direction of the outer
731 /// unit normal.
732 template<class ELEMENT>
733 void setup_boundary_coordinates(const unsigned& b,
734 const bool& switch_normal)
735 {
736 // Dummy file
737 std::ofstream some_file;
738
739 this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, some_file);
740 }
741
742
743 /// Setup boundary coordinate on boundary b which is
744 /// assumed to be planar. Boundary coordinates are the
745 /// x-y coordinates in the plane of that boundary, with the
746 /// x-axis along the line from the (lexicographically)
747 /// "lower left" to the "upper right" node. The y axis
748 /// is obtained by taking the cross-product of the positive
749 /// x direction with the outer unit normal computed by
750 /// the face elements (or its negative if switch_normal is set
751 /// to true). Doc faces in output file (if it's open).
752 ///
753 /// Note 1: Setup of boundary coordinates is not done if the boundary in
754 /// question turns out to be nonplanar.
755 ///
756 /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
757 /// also
758 /// store the boundary coordinates of its vertices. They are needed
759 /// to interpolated intrinsic coordinates of an associated
760 /// GeomObject (if any) into the interior.
761 /// Boolean argument allows switching of the direction of the outer
762 /// unit normal. Output file for doc.
763 template<class ELEMENT>
764 void setup_boundary_coordinates(const unsigned& b,
765 const bool& switch_normal,
766 std::ofstream& outfile);
767
768 /// Setup boundary coordinate on boundary b which is
769 /// assumed to be planar. Boundary coordinates are the
770 /// x-y coordinates in the plane of that boundary, with the
771 /// x-axis along the line from the (lexicographically)
772 /// "lower left" to the "upper right" node. The y axis
773 /// is obtained by taking the cross-product of the positive
774 /// x direction with the outer unit normal computed by
775 /// the face elements (or its negative if switch_normal is set
776 /// to true). Doc faces in output file (if it's open).
777 ///
778 /// Note 1: Setup of boundary coordinates is not done if the boundary in
779 /// question turns out to be nonplanar.
780 ///
781 /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
782 /// also
783 /// store the boundary coordinates of its vertices. They are needed
784 /// to interpolated intrinsic coordinates of an associated
785 /// GeomObject (if any) into the interior.
786 /// Output file for doc.
787 template<class ELEMENT>
788 void setup_boundary_coordinates(const unsigned& b, std::ofstream& outfile)
789 {
790 // Don't switch the normal
791 bool switch_normal = false;
792 this->setup_boundary_coordinates<ELEMENT>(b, switch_normal, outfile);
793 }
794
795
796 /// Return the number of elements adjacent to boundary b in region r
797 inline unsigned nboundary_element_in_region(const unsigned& b,
798 const unsigned& r) const
799 {
800 // Need to use a constant iterator here to keep the function "const"
801 // Return an iterator to the appropriate entry, if we find it
802 std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
804 if (it != Boundary_region_element_pt[b].end())
805 {
806 return (it->second).size();
807 }
808 // Otherwise there are no elements adjacent to boundary b in the region r
809 else
810 {
811 return 0;
812 }
813 }
814
815 /// Return pointer to the e-th element adjacent to boundary b in region r
817 const unsigned& r,
818 const unsigned& e) const
819 {
820 // Use a constant iterator here to keep function "const" overall
821 std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
823 if (it != Boundary_region_element_pt[b].end())
824 {
825 return (it->second)[e];
826 }
827 else
828 {
829 return 0;
830 }
831 }
832
833 /// Return face index of the e-th element adjacent to boundary b in region r
834 int face_index_at_boundary_in_region(const unsigned& b,
835 const unsigned& r,
836 const unsigned& e) const
837 {
838 // Use a constant iterator here to keep function "const" overall
839 std::map<unsigned, Vector<int>>::const_iterator it =
841 if (it != Face_index_region_at_boundary[b].end())
842 {
843 return (it->second)[e];
844 }
845 else
846 {
847 std::ostringstream error_message;
848 error_message << "Face indices not set up for boundary " << b
849 << " in region " << r << "\n";
850 error_message << "This probably means that the boundary is not "
851 "adjacent to region\n";
852 throw OomphLibError(error_message.str(),
853 OOMPH_CURRENT_FUNCTION,
854 OOMPH_EXCEPTION_LOCATION);
855 }
856 }
857
858
859 /// Return the number of regions specified by attributes
860 unsigned nregion()
861 {
862 return Region_element_pt.size();
863 }
864
865 /// Return the number of elements in region r
866 unsigned nregion_element(const unsigned& r)
867 {
868 unsigned entry = 0;
869 bool found = false;
870 unsigned n = Region_attribute.size();
871 for (unsigned i = 0; i < n; i++)
872 {
873#ifdef PARANOID
874 double diff =
875 fabs(Region_attribute[i] -
876 static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
877 if (diff > 0.0)
878 {
879 std::ostringstream error_message;
880 error_message << "Region attributes should be unsigneds because we \n"
881 << "only use them to set region ids\n";
882 throw OomphLibError(error_message.str(),
883 OOMPH_CURRENT_FUNCTION,
884 OOMPH_EXCEPTION_LOCATION);
885 }
886#endif
887 if (static_cast<unsigned>(Region_attribute[i]) == r)
888 {
889 entry = i;
890 found = true;
891 break;
892 }
893 }
894 if (found)
895 {
896 return Region_element_pt[entry].size();
897 }
898 else
899 {
900 return 0;
901 }
902 }
903
904 /// Return the i-th region attribute (here only used as the
905 /// (assumed to be unsigned) region id
906 double region_attribute(const unsigned& i)
907 {
908 return Region_attribute[i];
909 }
910
911 /// Return the e-th element in the r-th region
912 FiniteElement* region_element_pt(const unsigned& r, const unsigned& e)
913 {
914 unsigned entry = 0;
915 bool found = false;
916 unsigned n = Region_attribute.size();
917 for (unsigned i = 0; i < n; i++)
918 {
919#ifdef PARANOID
920 double diff =
921 fabs(Region_attribute[i] -
922 static_cast<double>(static_cast<unsigned>(Region_attribute[i])));
923 if (diff > 0.0)
924 {
925 std::ostringstream error_message;
926 error_message << "Region attributes should be unsigneds because we \n"
927 << "only use the to set region ids\n";
928 throw OomphLibError(error_message.str(),
929 OOMPH_CURRENT_FUNCTION,
930 OOMPH_EXCEPTION_LOCATION);
931 }
932#endif
933 if (static_cast<unsigned>(Region_attribute[i]) == r)
934 {
935 entry = i;
936 found = true;
937 break;
938 }
939 }
940 if (found)
941 {
942 return Region_element_pt[entry][e];
943 }
944 else
945 {
946 return 0;
947 }
948 }
949
950 /// Snap boundaries specified by the IDs listed in boundary_id to
951 /// a quadratric surface, specified in the file
952 /// quadratic_surface_file_name. This is usually used with vmtk-based
953 /// meshes for which oomph-lib's xda to poly conversion code produces the
954 /// files "quadratic_fsi_boundary.dat" and
955 /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
956 /// boundary (for the fluid and the solid) and the quadratic representation
957 /// of the outer boundary of the solid. When used with these files, the flag
958 /// switch_normal should be set to true when calling the function for the
959 /// outer boundary of the solid. The DocInfo object can be used to label
960 /// optional output files. (Uses directory and label).
961 template<class ELEMENT>
963 const Vector<unsigned>& boundary_id,
964 const std::string& quadratic_surface_file_name,
965 const bool& switch_normal,
966 DocInfo& doc_info);
967
968 /// Snap boundaries specified by the IDs listed in boundary_id to
969 /// a quadratric surface, specified in the file
970 /// quadratic_surface_file_name. This is usually used with vmtk-based
971 /// meshes for which oomph-lib's xda to poly conversion code produces the
972 /// files "quadratic_fsi_boundary.dat" and
973 /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
974 /// boundary (for the fluid and the solid) and the quadratic representation
975 /// of the outer boundary of the solid. When used with these files, the flag
976 /// switch_normal should be set to true when calling the function for the
977 /// outer boundary of the solid.
978 template<class ELEMENT>
980 const Vector<unsigned>& boundary_id,
981 const std::string& quadratic_surface_file_name,
982 const bool& switch_normal)
983 {
984 // Dummy doc info
985 DocInfo doc_info;
986 doc_info.disable_doc();
987 snap_to_quadratic_surface<ELEMENT>(
988 boundary_id, quadratic_surface_file_name, switch_normal, doc_info);
989 }
990
991 /// Move the nodes on boundaries with associated GeomObjects so
992 /// that they exactly coincide with the geometric object. This requires
993 /// that the boundary coordinates are set up consistently.
995
996
997 /// Non-Delaunay split elements that have three faces on a boundary
998 /// into sons. Timestepper species timestepper for new nodes; defaults
999 /// to to steady timestepper.
1000 template<class ELEMENT>
1002 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper);
1003
1004
1005 /// Setup lookup schemes which establish which elements are located
1006 /// next to mesh's boundaries (wrapper to suppress doc).
1008 {
1009 std::ofstream outfile;
1010 this->setup_boundary_element_info(outfile);
1011 }
1012
1013
1014 /// Setup lookup schemes which establish which elements are located
1015 /// next to mesh's boundaries. Doc in outfile (if it's open).
1016 void setup_boundary_element_info(std::ostream& outfile);
1017
1018
1019 protected:
1020 /// Vectors of vectors of elements in each region (note: this just
1021 /// stores them; the region IDs are contained in Region_attribute!)
1023
1024 /// Vector of attributes associated with the elements in each region
1025 /// NOTE: double is enforced on us by tetgen. We use it as an unsigned
1026 /// to indicate the actual (zero-based) region ID
1028
1029 /// Storage for elements adjacent to a boundary in a particular region
1032
1033 /// Storage for the face index adjacent to a boundary in a particular
1034 /// region
1036
1037 /// Faceted surface that defines outer boundaries
1039
1040 /// Vector to faceted surfaces that define internal boundaries
1042
1043 /// Reverse lookup scheme: Pointer to faceted surface (if any!)
1044 /// associated with boundary b
1045 std::map<unsigned, TetMeshFacetedSurface*> Tet_mesh_faceted_surface_pt;
1046
1047 /// Reverse lookup scheme: Pointer to facet (if any!)
1048 /// associated with boundary b
1049 std::map<unsigned, TetMeshFacet*> Tet_mesh_facet_pt;
1050
1051 /// Boundary coordinates of vertices in triangular facets
1052 /// associated with given boundary. Is only set up for triangular facets!
1053 std::map<unsigned, Vector<Vector<double>>>
1055
1056 /// Timestepper used to build nodes
1058 };
1059
1060
1061 /// ///////////////////////////////////////////////////////////////////////////
1062 /// ///////////////////////////////////////////////////////////////////////////
1063 /// ///////////////////////////////////////////////////////////////////////////
1064
1065
1066 //###########################################################################
1067 // Templated member functions
1068 //###########################################################################
1069
1070
1071 //======================================================================
1072 /// Setup boundary coordinate on boundary b which is
1073 /// assumed to be planar. Boundary coordinates are the
1074 /// x-y coordinates in the plane of that boundary, with the
1075 /// x-axis along the line from the (lexicographically)
1076 /// "lower left" to the "upper right" node. The y axis
1077 /// is obtained by taking the cross-product of the positive
1078 /// x direction with the outer unit normal computed by
1079 /// the face elements (or its negative if switch_normal is set
1080 /// to true). Doc faces in output file (if it's open).
1081 ///
1082 /// Note 1: Setup of boundary coordinates is not done if the boundary in
1083 /// question turns out to be nonplanar.
1084 ///
1085 /// Note 2: If a triangular TetMeshFacet is associated with a boundary we
1086 /// also
1087 /// store the boundary coordinates of its vertices. They are needed
1088 /// to interpolated intrinsic coordinates of an associated GeomObject
1089 /// (if any) into the interior.
1090 //======================================================================
1091 template<class ELEMENT>
1093 const bool& switch_normal,
1094 std::ofstream& outfile)
1095 {
1096 Node* lower_left_node_pt = 0;
1097 Node* upper_right_node_pt = 0;
1098
1099 // Unit vector connecting lower left and upper right nodes
1100 Vector<double> b0(3);
1101
1102 // ...and a vector normal to it
1103 Vector<double> b1(3);
1104
1105
1106 // Facet?
1107 TetMeshFacet* f_pt = 0;
1108 std::map<unsigned, TetMeshFacet*>::iterator it = Tet_mesh_facet_pt.find(b);
1109 if (it != Tet_mesh_facet_pt.end())
1110 {
1111 f_pt = (*it).second;
1112 }
1113
1114 // std::cout << "Debug " << b; f_pt->output(std::cout);
1115
1116 // Number of vertices
1117 unsigned nv = 0;
1118 if (f_pt != 0)
1119 {
1120 nv = f_pt->nvertex();
1121 }
1122
1123 // Check for planarity and bail out if nodes are not in the same plane
1124 bool vertices_are_in_same_plane = true;
1125 for (unsigned do_for_real = 0; do_for_real < 2; do_for_real++)
1126 {
1127 // Temporary storage for face elements
1128 Vector<FiniteElement*> face_el_pt;
1129
1130 // Loop over all elements on boundaries
1131 unsigned nel = this->nboundary_element(b);
1132
1133 if (nel > 0)
1134 {
1135 // Loop over the bulk elements adjacent to boundary b
1136 for (unsigned e = 0; e < nel; e++)
1137 {
1138 // Get pointer to the bulk element that is adjacent to boundary b
1139 FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
1140
1141 // Find the index of the face of element e along boundary b
1142 int face_index = this->face_index_at_boundary(b, e);
1143
1144 // Create new face element
1145 face_el_pt.push_back(
1146 new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index));
1147
1148 // Output faces?
1149 if (outfile.is_open())
1150 {
1151 face_el_pt[face_el_pt.size() - 1]->output(outfile);
1152 }
1153 }
1154
1155 // Loop over all nodes to find the lower left and upper right ones
1156 lower_left_node_pt = this->boundary_node_pt(b, 0);
1157 upper_right_node_pt = this->boundary_node_pt(b, 0);
1158
1159 // Loop over all nodes on the boundary
1160 unsigned nnod = this->nboundary_node(b);
1161 for (unsigned j = 0; j < nnod; j++)
1162 {
1163 // Get node
1164 Node* nod_pt = this->boundary_node_pt(b, j);
1165
1166 // Primary criterion for lower left: Does it have a lower
1167 // z-coordinate?
1168 if (nod_pt->x(2) < lower_left_node_pt->x(2))
1169 {
1170 lower_left_node_pt = nod_pt;
1171 }
1172 // ...or is it a draw?
1173 else if (nod_pt->x(2) == lower_left_node_pt->x(2))
1174 {
1175 // If it's a draw: Does it have a lower y-coordinate?
1176 if (nod_pt->x(1) < lower_left_node_pt->x(1))
1177 {
1178 lower_left_node_pt = nod_pt;
1179 }
1180 // ...or is it a draw?
1181 else if (nod_pt->x(1) == lower_left_node_pt->x(1))
1182 {
1183 // If it's a draw: Does it have a lower x-coordinate?
1184 if (nod_pt->x(0) < lower_left_node_pt->x(0))
1185 {
1186 lower_left_node_pt = nod_pt;
1187 }
1188 }
1189 }
1190
1191 // Primary criterion for upper right: Does it have a higher
1192 // z-coordinate?
1193 if (nod_pt->x(2) > upper_right_node_pt->x(2))
1194 {
1195 upper_right_node_pt = nod_pt;
1196 }
1197 // ...or is it a draw?
1198 else if (nod_pt->x(2) == upper_right_node_pt->x(2))
1199 {
1200 // If it's a draw: Does it have a higher y-coordinate?
1201 if (nod_pt->x(1) > upper_right_node_pt->x(1))
1202 {
1203 upper_right_node_pt = nod_pt;
1204 }
1205 // ...or is it a draw?
1206 else if (nod_pt->x(1) == upper_right_node_pt->x(1))
1207 {
1208 // If it's a draw: Does it have a higher x-coordinate?
1209 if (nod_pt->x(0) > upper_right_node_pt->x(0))
1210 {
1211 upper_right_node_pt = nod_pt;
1212 }
1213 }
1214 }
1215 }
1216
1217 // Prepare storage for boundary coordinates
1218 Vector<double> zeta(2);
1219 /*std::cout << upper_right_node_pt->x(0) << " "
1220 << upper_right_node_pt->x(1) << " "
1221 << upper_right_node_pt->x(2) << " L ";
1222 std::cout << lower_left_node_pt->x(0) << " "
1223 << lower_left_node_pt->x(1) << " "
1224 << lower_left_node_pt->x(2) << "\n ";*/
1225
1226
1227 // Unit vector connecting lower left and upper right nodes
1228 b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
1229 b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
1230 b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
1231
1232 // Normalise
1233 double inv_length =
1234 1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
1235 b0[0] *= inv_length;
1236 b0[1] *= inv_length;
1237 b0[2] *= inv_length;
1238
1239 // std::cout << "B0 ";
1240 // std::cout << b0[0] << " " << b0[1] << " " << b0[2] << "\n";
1241
1242 // Get (outer) unit normal to first face element
1243 Vector<double> normal(3);
1244 Vector<double> s(2, 0.0);
1245 if (nv != 3)
1246 {
1247 dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[0])
1248 ->outer_unit_normal(s, normal);
1249 }
1250 else
1251 {
1252 double t1x = f_pt->vertex_pt(1)->x(0) - f_pt->vertex_pt(0)->x(0);
1253
1254 double t1y = f_pt->vertex_pt(1)->x(1) - f_pt->vertex_pt(0)->x(1);
1255
1256 double t1z = f_pt->vertex_pt(1)->x(2) - f_pt->vertex_pt(0)->x(2);
1257
1258 double t2x = f_pt->vertex_pt(2)->x(0) - f_pt->vertex_pt(0)->x(0);
1259
1260 double t2y = f_pt->vertex_pt(2)->x(1) - f_pt->vertex_pt(0)->x(1);
1261
1262 double t2z = f_pt->vertex_pt(2)->x(2) - f_pt->vertex_pt(0)->x(2);
1263
1264 normal[0] = t1y * t2z - t1z * t2y;
1265 normal[1] = t1z * t2x - t1x * t2z;
1266 normal[2] = t1x * t2y - t1y * t2x;
1267 double inv_length =
1268 1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
1269 normal[2] * normal[2]);
1270 normal[0] *= inv_length;
1271 normal[1] *= inv_length;
1272 normal[2] *= inv_length;
1273 }
1274
1275 if (switch_normal)
1276 {
1277 normal[0] = -normal[0];
1278 normal[1] = -normal[1];
1279 normal[2] = -normal[2];
1280 }
1281
1282 // std::cout << "Normal ";
1283 // std::cout << normal[0] << " " << normal[1] << " " << normal[2] <<
1284 // "\n";
1285
1286
1287 // Check that all elements have the same normal
1288 for (unsigned e = 0; e < nel; e++)
1289 {
1290 // Get (outer) unit normal to face element
1291 Vector<double> my_normal(3);
1292 dynamic_cast<DummyFaceElement<ELEMENT>*>(face_el_pt[e])
1293 ->outer_unit_normal(s, my_normal);
1294
1295 // Dot product should be one!
1296 double dot_prod = normal[0] * my_normal[0] +
1297 normal[1] * my_normal[1] + normal[2] * my_normal[2];
1298
1299
1300 double error = abs(dot_prod) - 1.0;
1301 if (abs(error) > Tolerance_for_boundary_finding)
1302 {
1303 vertices_are_in_same_plane = false;
1304 }
1305 }
1306
1307 // Bail out?
1308 if (vertices_are_in_same_plane)
1309 {
1310 // Cross-product to get second in-plane vector, normal to b0
1311 b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
1312 b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
1313 b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
1314
1315 // std::cout << "B1 ";
1316 // std::cout << b1[0] << " " << b1[1] << " " << b1[2] << "\n";
1317
1318
1319 // Assign boundary coordinates: projection onto the axes
1320 for (unsigned j = 0; j < nnod; j++)
1321 {
1322 // Get node
1323 Node* nod_pt = this->boundary_node_pt(b, j);
1324
1325 // Difference vector to lower left corner
1326 Vector<double> delta(3);
1327 delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
1328 delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
1329 delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
1330
1331 /*std::cout << j << ": "
1332 << nod_pt->x(0) << " " << nod_pt->x(1)
1333 << " " << nod_pt->x(2) << " Delta ";
1334 std::cout << delta[0] << " " << delta[1] << " " << delta[2] << "\n";
1335 */
1336
1337 // Project
1338 zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1339 zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1340
1341 // Check:
1342 double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
1343 zeta[1] * b1[0] - nod_pt->x(0),
1344 2) +
1345 pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
1346 zeta[1] * b1[1] - nod_pt->x(1),
1347 2) +
1348 pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
1349 zeta[1] * b1[2] - nod_pt->x(2),
1350 2);
1351
1352 if (sqrt(error) > Tolerance_for_boundary_finding)
1353 {
1354 std::ostringstream error_message;
1355 error_message
1356 << "Error in projection of boundary coordinates = "
1357 << sqrt(error) << " > Tolerance_for_boundary_finding = "
1358 << Tolerance_for_boundary_finding << std::endl
1359 << "nv = " << nv << std::endl
1360 << "Lower left: " << lower_left_node_pt->x(0) << " "
1361 << lower_left_node_pt->x(1) << " " << lower_left_node_pt->x(2)
1362 << " " << std::endl
1363 << "Upper right: " << upper_right_node_pt->x(0) << " "
1364 << upper_right_node_pt->x(1) << " " << upper_right_node_pt->x(2)
1365 << " " << std::endl
1366 << "Nodes: ";
1367 for (unsigned j = 0; j < nnod; j++)
1368 {
1369 // Get node
1370 Node* nod_pt = this->boundary_node_pt(b, j);
1371 error_message << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1372 << nod_pt->x(2) << " " << std::endl;
1373 }
1374 throw OomphLibError(error_message.str(),
1375 OOMPH_CURRENT_FUNCTION,
1376 OOMPH_EXCEPTION_LOCATION);
1377 }
1378
1379 // Set boundary coordinate
1380 if (do_for_real == 1)
1381 {
1382 nod_pt->set_coordinates_on_boundary(b, zeta);
1383 }
1384 }
1385 } // End of if vertices are in the same plane
1386 }
1387
1388
1389 // Indicate that boundary coordinate has been set up
1390 if (do_for_real == 1)
1391 {
1393
1394 // Facet associated with this boundary?
1395 if (f_pt != 0)
1396 {
1397 // Triangular facet: Set coordinates at vertices
1398 if (nv == 3)
1399 {
1401 for (unsigned j = 0; j < 3; j++)
1402 {
1403 // Two surface coordinates
1405
1406 // Difference vector to lower left corner
1407 Vector<double> delta(3);
1408 delta[0] = f_pt->vertex_pt(j)->x(0) - lower_left_node_pt->x(0);
1409 delta[1] = f_pt->vertex_pt(j)->x(1) - lower_left_node_pt->x(1);
1410 delta[2] = f_pt->vertex_pt(j)->x(2) - lower_left_node_pt->x(2);
1411
1412 // Project
1413 Vector<double> zeta(2);
1414 zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1415 zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1416
1417 for (unsigned ii = 0; ii < 2; ii++)
1418 {
1420 zeta[ii];
1421 }
1422 }
1423 }
1424 }
1425 }
1426
1427 // Cleanup
1428 unsigned n = face_el_pt.size();
1429 for (unsigned e = 0; e < n; e++)
1430 {
1431 delete face_el_pt[e];
1432 }
1433
1434 // Bail out?
1435 if (!vertices_are_in_same_plane)
1436 {
1437 /* oomph_info << "Vertices on boundary " << b */
1438 /* << " are not in the same plane; bailing out\n"; */
1439 return;
1440 }
1441 }
1442 }
1443
1444
1445 //======================================================================
1446 /// Snap boundaries specified by the IDs listed in boundary_id to
1447 /// a quadratric surface, specified in the file
1448 /// quadratic_surface_file_name. This is usually used with vmtk-based
1449 /// meshes for which oomph-lib's xda to poly conversion code produces the
1450 /// files "quadratic_fsi_boundary.dat" and
1451 /// "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI
1452 /// boundary (for the fluid and the solid) and the quadratic representation of
1453 /// the outer boundary of the solid. When used with these files, the flag
1454 /// switch_normal should be set to true when calling the function for the
1455 /// outer boundary of the solid. The DocInfo object can be used to label
1456 /// optional output files. (Uses directory and label).
1457 //======================================================================
1458 template<class ELEMENT>
1460 const Vector<unsigned>& boundary_id,
1461 const std::string& quadratic_surface_file_name,
1462 const bool& switch_normal,
1463 DocInfo& doc_info)
1464 {
1465 // Aux storage for processing input
1466 char dummy[101];
1467
1468 // Prepare to doc nodes that couldn't be snapped
1469 std::ofstream the_file_non_snapped;
1470 std::string non_snapped_filename =
1471 "non_snapped_nodes_" + doc_info.label() + ".dat";
1472
1473 // Read the number of nodes and elements (quadratic facets)
1474 std::ifstream infile(quadratic_surface_file_name.c_str(),
1475 std::ios_base::in);
1476 unsigned n_node;
1477 infile >> n_node;
1478
1479 // Ignore rest of line
1480 infile.getline(dummy, 100);
1481
1482 // Number of quadratic facets
1483 unsigned nel;
1484 infile >> nel;
1485
1486 // Ignore rest of line
1487 infile.getline(dummy, 100);
1488
1489 // Check that the number of elements matches
1490 if (nel != boundary_id.size())
1491 {
1492 std::ostringstream error_message;
1493 error_message << "Number of quadratic facets specified in "
1494 << quadratic_surface_file_name << "is: " << nel
1495 << "\nThis doesn't match the number of planar boundaries \n"
1496 << "specified in boundary_id which is "
1497 << boundary_id.size() << std::endl;
1498 throw OomphLibError(
1499 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1500 }
1501
1502 // Temporary storage for face elements
1504
1505 // Loop over quadratic face elements -- one for each facet
1506 for (unsigned e = 0; e < nel; e++)
1507 {
1508 face_el_pt.push_back(new FreeStandingFaceElement<ELEMENT>);
1509 }
1510
1511
1512 // Now build nodes
1513 unsigned n_dim = 3;
1514 unsigned n_position_type = 1;
1515 unsigned initial_n_value = 0;
1516 Vector<Node*> nod_pt(n_node);
1517 unsigned node_nmbr;
1518 std::map<unsigned, unsigned> node_number;
1519 std::ofstream node_file;
1520 for (unsigned j = 0; j < n_node; j++)
1521 {
1522 nod_pt[j] =
1523 new BoundaryNode<Node>(n_dim, n_position_type, initial_n_value);
1524 infile >> nod_pt[j]->x(0);
1525 infile >> nod_pt[j]->x(1);
1526 infile >> nod_pt[j]->x(2);
1527 infile >> node_nmbr;
1528 node_number[node_nmbr] = j;
1529 }
1530
1531
1532 // Now assign nodes to elements -- each element represents
1533 // distinct boundary; assign enumeration as specified by
1534 // boundary_id.
1535 for (unsigned e = 0; e < nel; e++)
1536 {
1537 unsigned nnod_el = face_el_pt[e]->nnode();
1538 unsigned j_global;
1539 for (unsigned j = 0; j < nnod_el; j++)
1540 {
1541 infile >> j_global;
1542 face_el_pt[e]->node_pt(j) = nod_pt[node_number[j_global]];
1543 face_el_pt[e]->node_pt(j)->add_to_boundary(boundary_id[e]);
1544 }
1545 face_el_pt[e]->set_boundary_number_in_bulk_mesh(boundary_id[e]);
1546 face_el_pt[e]->set_nodal_dimension(3);
1547 }
1548
1549
1550 // Setup boundary coordinates for each facet, using
1551 // the same strategy as for the simplex boundaries
1552 // (there's some code duplication here but it doesn't
1553 // seem worth it to rationlise this as the interface would
1554 // be pretty clunky).
1555 for (unsigned e = 0; e < nel; e++)
1556 {
1557 Vector<Vector<double>> vertex_boundary_coord(3);
1558
1559 // Loop over all nodes to find the lower left and upper right ones
1560 Node* lower_left_node_pt = face_el_pt[e]->node_pt(0);
1561 Node* upper_right_node_pt = face_el_pt[e]->node_pt(0);
1562
1563 // Loop over all vertex nodes
1564 for (unsigned j = 0; j < 3; j++)
1565 {
1566 // Get node
1567 Node* nod_pt = face_el_pt[e]->node_pt(j);
1568
1569 // Primary criterion for lower left: Does it have a lower z-coordinate?
1570 if (nod_pt->x(2) < lower_left_node_pt->x(2))
1571 {
1572 lower_left_node_pt = nod_pt;
1573 }
1574 // ...or is it a draw?
1575 else if (nod_pt->x(2) == lower_left_node_pt->x(2))
1576 {
1577 // If it's a draw: Does it have a lower y-coordinate?
1578 if (nod_pt->x(1) < lower_left_node_pt->x(1))
1579 {
1580 lower_left_node_pt = nod_pt;
1581 }
1582 // ...or is it a draw?
1583 else if (nod_pt->x(1) == lower_left_node_pt->x(1))
1584 {
1585 // If it's a draw: Does it have a lower x-coordinate?
1586 if (nod_pt->x(0) < lower_left_node_pt->x(0))
1587 {
1588 lower_left_node_pt = nod_pt;
1589 }
1590 }
1591 }
1592
1593 // Primary criterion for upper right: Does it have a higher
1594 // z-coordinate?
1595 if (nod_pt->x(2) > upper_right_node_pt->x(2))
1596 {
1597 upper_right_node_pt = nod_pt;
1598 }
1599 // ...or is it a draw?
1600 else if (nod_pt->x(2) == upper_right_node_pt->x(2))
1601 {
1602 // If it's a draw: Does it have a higher y-coordinate?
1603 if (nod_pt->x(1) > upper_right_node_pt->x(1))
1604 {
1605 upper_right_node_pt = nod_pt;
1606 }
1607 // ...or is it a draw?
1608 else if (nod_pt->x(1) == upper_right_node_pt->x(1))
1609 {
1610 // If it's a draw: Does it have a higher x-coordinate?
1611 if (nod_pt->x(0) > upper_right_node_pt->x(0))
1612 {
1613 upper_right_node_pt = nod_pt;
1614 }
1615 }
1616 }
1617 }
1618
1619 // Prepare storage for boundary coordinates
1620 Vector<double> zeta(2);
1621
1622 // Unit vector connecting lower left and upper right nodes
1623 Vector<double> b0(3);
1624 b0[0] = upper_right_node_pt->x(0) - lower_left_node_pt->x(0);
1625 b0[1] = upper_right_node_pt->x(1) - lower_left_node_pt->x(1);
1626 b0[2] = upper_right_node_pt->x(2) - lower_left_node_pt->x(2);
1627
1628 // Normalise
1629 double inv_length =
1630 1.0 / sqrt(b0[0] * b0[0] + b0[1] * b0[1] + b0[2] * b0[2]);
1631 b0[0] *= inv_length;
1632 b0[1] *= inv_length;
1633 b0[2] *= inv_length;
1634
1635 // Get (outer) unit normal to face element -- note that
1636 // with the enumeration chosen in oomph-lib's xda to poly
1637 // conversion code the sign below is correct for the outer
1638 // unit normal on the FSI interface.
1639 Vector<double> tang1(3);
1640 Vector<double> tang2(3);
1641 Vector<double> normal(3);
1642 tang1[0] =
1643 face_el_pt[e]->node_pt(1)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
1644 tang1[1] =
1645 face_el_pt[e]->node_pt(1)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
1646 tang1[2] =
1647 face_el_pt[e]->node_pt(1)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
1648 tang2[0] =
1649 face_el_pt[e]->node_pt(2)->x(0) - face_el_pt[e]->node_pt(0)->x(0);
1650 tang2[1] =
1651 face_el_pt[e]->node_pt(2)->x(1) - face_el_pt[e]->node_pt(0)->x(1);
1652 tang2[2] =
1653 face_el_pt[e]->node_pt(2)->x(2) - face_el_pt[e]->node_pt(0)->x(2);
1654 normal[0] = -(tang1[1] * tang2[2] - tang1[2] * tang2[1]);
1655 normal[1] = -(tang1[2] * tang2[0] - tang1[0] * tang2[2]);
1656 normal[2] = -(tang1[0] * tang2[1] - tang1[1] * tang2[0]);
1657
1658 // Normalise
1659 inv_length = 1.0 / sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
1660 normal[2] * normal[2]);
1661 normal[0] *= inv_length;
1662 normal[1] *= inv_length;
1663 normal[2] *= inv_length;
1664
1665 // Change direction -- usually for outer boundary of solid
1666 if (switch_normal)
1667 {
1668 normal[0] = -normal[0];
1669 normal[1] = -normal[1];
1670 normal[2] = -normal[2];
1671 }
1672
1673 // Cross-product to get second in-plane vector, normal to b0
1674 Vector<double> b1(3);
1675 b1[0] = b0[1] * normal[2] - b0[2] * normal[1];
1676 b1[1] = b0[2] * normal[0] - b0[0] * normal[2];
1677 b1[2] = b0[0] * normal[1] - b0[1] * normal[0];
1678
1679 // Assign boundary coordinates for vertex nodes: projection onto the axes
1680 for (unsigned j = 0; j < 3; j++)
1681 {
1682 // Get node
1683 Node* nod_pt = face_el_pt[e]->node_pt(j);
1684
1685 // Difference vector to lower left corner
1686 Vector<double> delta(3);
1687 delta[0] = nod_pt->x(0) - lower_left_node_pt->x(0);
1688 delta[1] = nod_pt->x(1) - lower_left_node_pt->x(1);
1689 delta[2] = nod_pt->x(2) - lower_left_node_pt->x(2);
1690
1691 // Project
1692 zeta[0] = delta[0] * b0[0] + delta[1] * b0[1] + delta[2] * b0[2];
1693 zeta[1] = delta[0] * b1[0] + delta[1] * b1[1] + delta[2] * b1[2];
1694
1695 vertex_boundary_coord[j].resize(2);
1696 vertex_boundary_coord[j][0] = zeta[0];
1697 vertex_boundary_coord[j][1] = zeta[1];
1698
1699#ifdef PARANOID
1700
1701 // Check:
1702 double error = pow(lower_left_node_pt->x(0) + zeta[0] * b0[0] +
1703 zeta[1] * b1[0] - nod_pt->x(0),
1704 2) +
1705 pow(lower_left_node_pt->x(1) + zeta[0] * b0[1] +
1706 zeta[1] * b1[1] - nod_pt->x(1),
1707 2) +
1708 pow(lower_left_node_pt->x(2) + zeta[0] * b0[2] +
1709 zeta[1] * b1[2] - nod_pt->x(2),
1710 2);
1711
1712 if (sqrt(error) > Tolerance_for_boundary_finding)
1713 {
1714 std::ostringstream error_message;
1715 error_message
1716 << "Error in setup of boundary coordinate " << sqrt(error)
1717 << std::endl
1718 << "exceeds tolerance specified by the static member data\n"
1719 << "TetMeshBase::Tolerance_for_boundary_finding = "
1720 << Tolerance_for_boundary_finding << std::endl
1721 << "This usually means that the boundary is not planar.\n\n"
1722 << "You can suppress this error message by recompiling \n"
1723 << "recompiling without PARANOID or by changing the tolerance.\n";
1724 throw OomphLibError(error_message.str(),
1725 OOMPH_CURRENT_FUNCTION,
1726 OOMPH_EXCEPTION_LOCATION);
1727 }
1728#endif
1729
1730 // Set boundary coordinate
1731 nod_pt->set_coordinates_on_boundary(boundary_id[e], zeta);
1732 }
1733
1734 // Assign boundary coordinates of three midside nodes by linear
1735 // interpolation (corresponding to a flat facet)
1736
1737 // Node 3 is between 0 and 1
1738 zeta[0] =
1739 0.5 * (vertex_boundary_coord[0][0] + vertex_boundary_coord[1][0]);
1740 zeta[1] =
1741 0.5 * (vertex_boundary_coord[0][1] + vertex_boundary_coord[1][1]);
1742 face_el_pt[e]->node_pt(3)->set_coordinates_on_boundary(boundary_id[e],
1743 zeta);
1744
1745 // Node 4 is between 1 and 2
1746 zeta[0] =
1747 0.5 * (vertex_boundary_coord[1][0] + vertex_boundary_coord[2][0]);
1748 zeta[1] =
1749 0.5 * (vertex_boundary_coord[1][1] + vertex_boundary_coord[2][1]);
1750 face_el_pt[e]->node_pt(4)->set_coordinates_on_boundary(boundary_id[e],
1751 zeta);
1752
1753 // Node 5 is between 2 and 0
1754 zeta[0] =
1755 0.5 * (vertex_boundary_coord[2][0] + vertex_boundary_coord[0][0]);
1756 zeta[1] =
1757 0.5 * (vertex_boundary_coord[2][1] + vertex_boundary_coord[0][1]);
1758 face_el_pt[e]->node_pt(5)->set_coordinates_on_boundary(boundary_id[e],
1759 zeta);
1760 }
1761
1762
1763 // Loop over elements/facets = boundaries to snap
1764 bool success = true;
1765 for (unsigned b = 0; b < nel; b++)
1766 {
1767 // Doc boundary coordinates on quadratic patches
1768 std::ofstream the_file;
1769 std::ofstream the_file_before;
1770 std::ofstream the_file_after;
1771 if (doc_info.is_doc_enabled())
1772 {
1773 std::ostringstream filename;
1774 filename << doc_info.directory() << "/quadratic_coordinates_"
1775 << doc_info.label() << b << ".dat";
1776 the_file.open(filename.str().c_str());
1777
1778 std::ostringstream filename1;
1779 filename1 << doc_info.directory() << "/quadratic_nodes_before_"
1780 << doc_info.label() << b << ".dat";
1781 the_file_before.open(filename1.str().c_str());
1782
1783 std::ostringstream filename2;
1784 filename2 << doc_info.directory() << "/quadratic_nodes_after_"
1785 << doc_info.label() << b << ".dat";
1786 the_file_after.open(filename2.str().c_str());
1787 }
1788
1789 // Cast the element pointer
1790 FreeStandingFaceElement<ELEMENT>* el_pt = face_el_pt[b];
1791
1792 // Doc boundary coordinate on quadratic facet representation
1793 if (doc_info.is_doc_enabled())
1794 {
1795 Vector<double> s(2);
1796 Vector<double> zeta(2);
1797 Vector<double> x(3);
1798 unsigned n_plot = 3;
1799 the_file << el_pt->tecplot_zone_string(n_plot);
1800
1801 // Loop over plot points
1802 unsigned num_plot_points = el_pt->nplot_points(n_plot);
1803 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1804 {
1805 // Get local coordinates of plot point
1806 el_pt->get_s_plot(iplot, n_plot, s);
1807 el_pt->interpolated_zeta(s, zeta);
1808 el_pt->interpolated_x(s, x);
1809 for (unsigned i = 0; i < 3; i++)
1810 {
1811 the_file << x[i] << " ";
1812 }
1813 for (unsigned i = 0; i < 2; i++)
1814 {
1815 the_file << zeta[i] << " ";
1816 }
1817 the_file << std::endl;
1818 }
1819 el_pt->write_tecplot_zone_footer(the_file, n_plot);
1820
1821 // std::cout << "Facet " << b << " corresponds to quadratic
1822 // boundary "
1823 // << boundary_id[b] << " which contains "
1824 // << this->nboundary_element(boundary_id[b])
1825 // << " element[s] " << std::endl;
1826 }
1827
1828
1829 // Loop over bulk elements that are adjacent to quadratic boundary
1830 Vector<double> boundary_zeta(2);
1831 Vector<double> quadratic_coordinates(3);
1832 GeomObject* geom_obj_pt = 0;
1833 Vector<double> s_geom_obj(2);
1834 unsigned nel = this->nboundary_element(boundary_id[b]);
1835 for (unsigned e = 0; e < nel; e++)
1836 {
1837 // Get pointer to the bulk element that is adjacent to boundary
1838 FiniteElement* bulk_elem_pt =
1839 this->boundary_element_pt(boundary_id[b], e);
1840
1841 // Loop over nodes
1842 unsigned nnod = bulk_elem_pt->nnode();
1843 for (unsigned j = 0; j < nnod; j++)
1844 {
1845 Node* nod_pt = bulk_elem_pt->node_pt(j);
1846 if (nod_pt->is_on_boundary(boundary_id[b]))
1847 {
1848 nod_pt->get_coordinates_on_boundary(boundary_id[b], boundary_zeta);
1849
1850 // Doc it?
1851 if (doc_info.is_doc_enabled())
1852 {
1853 the_file_before << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1854 << nod_pt->x(2) << " " << boundary_zeta[0] << " "
1855 << boundary_zeta[1] << " " << std::endl;
1856 }
1857
1858 // Find local coordinate in quadratic facet
1859 el_pt->locate_zeta(boundary_zeta, geom_obj_pt, s_geom_obj);
1860 if (geom_obj_pt != 0)
1861 {
1862 geom_obj_pt->position(s_geom_obj, quadratic_coordinates);
1863 nod_pt->x(0) = quadratic_coordinates[0];
1864 nod_pt->x(1) = quadratic_coordinates[1];
1865 nod_pt->x(2) = quadratic_coordinates[2];
1866 }
1867 else
1868 {
1869 // Get ready to bail out below...
1870 success = false;
1871
1872 std::ostringstream error_message;
1873 error_message
1874 << "Warning: Couldn't find GeomObject during snapping to\n"
1875 << "quadratic surface for boundary " << boundary_id[b]
1876 << ". I'm leaving the node where it was. Will bail out "
1877 "later.\n";
1878 OomphLibWarning(error_message.str(),
1879 "TetgenMesh::snap_to_quadratic_surface()",
1880 OOMPH_EXCEPTION_LOCATION);
1881 if (!the_file_non_snapped.is_open())
1882 {
1883 the_file_non_snapped.open(non_snapped_filename.c_str());
1884 }
1885 the_file_non_snapped << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1886 << nod_pt->x(2) << " " << boundary_zeta[0]
1887 << " " << boundary_zeta[1] << " "
1888 << std::endl;
1889 }
1890
1891 // Doc it?
1892 if (doc_info.is_doc_enabled())
1893 {
1894 the_file_after << nod_pt->x(0) << " " << nod_pt->x(1) << " "
1895 << nod_pt->x(2) << " " << boundary_zeta[0] << " "
1896 << boundary_zeta[1] << " " << std::endl;
1897 }
1898 }
1899 }
1900 }
1901
1902 // Close doc file
1903 the_file.close();
1904 the_file_before.close();
1905 the_file_after.close();
1906 }
1907
1908 // Bail out?
1909 if (!success)
1910 {
1911 std::ostringstream error_message;
1912 error_message
1913 << "Warning: Couldn't find GeomObject during snapping to\n"
1914 << "quadratic surface. Bailing out.\n"
1915 << "Nodes that couldn't be snapped are contained in \n"
1916 << "file: " << non_snapped_filename << ".\n"
1917 << "This problem may arise because the switch_normal flag was \n"
1918 << "set wrongly.\n";
1919 throw OomphLibError(
1920 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1921 }
1922
1923 // Close
1924 if (!the_file_non_snapped.is_open())
1925 {
1926 the_file_non_snapped.close();
1927 }
1928
1929 // Kill auxiliary FreeStandingFaceElements
1930 for (unsigned e = 0; e < nel; e++)
1931 {
1932 delete face_el_pt[e];
1933 face_el_pt[e] = 0;
1934 }
1935
1936 // Kill boundary nodes
1937 unsigned nn = nod_pt.size();
1938 for (unsigned j = 0; j < nn; j++)
1939 {
1940 delete nod_pt[j];
1941 }
1942 }
1943
1944
1945 //========================================================================
1946 /// Non-delaunay split elements that have three faces on a boundary
1947 /// into sons.
1948 //========================================================================
1949 template<class ELEMENT>
1951 {
1952 // Setup boundary lookup scheme if required
1954 {
1956 }
1957
1958 // Find out how many nodes we have along each element edge
1959 unsigned nnode_1d = finite_element_pt(0)->nnode_1d();
1960 // Find out the total number of nodes
1961 unsigned nnode = this->finite_element_pt(0)->nnode();
1962
1963 // At the moment we're only able to deal with nnode_1d=2 or 3.
1964 if ((nnode_1d != 2) && (nnode_1d != 3))
1965 {
1966 std::ostringstream error_message;
1967 error_message << "Mesh generation from tetgen currently only works\n";
1968 error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
1969 error_message << "for nnode_1d=" << nnode_1d << std::endl;
1970
1971 throw OomphLibError(
1972 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1973 }
1974
1975 // Map to store how many boundaries elements are located on
1976 std::map<FiniteElement*, unsigned> count;
1977
1978 // Loop over all boundaries
1979 unsigned nb = this->nboundary();
1980 for (unsigned b = 0; b < nb; b++)
1981 {
1982 // Loop over elements next to that boundary
1983 unsigned nel = this->nboundary_element(b);
1984 for (unsigned e = 0; e < nel; e++)
1985 {
1986 // Get pointer to element
1987 FiniteElement* el_pt = boundary_element_pt(b, e);
1988
1989 // Bump up counter
1990 count[el_pt]++;
1991 }
1992 }
1993
1994 // To avoid having to check the map for all elements (which will
1995 // inflate it to the size of all elements!), move offending elements
1996 // into set
1997 std::set<FiniteElement*> elements_to_be_split;
1998 for (std::map<FiniteElement*, unsigned>::iterator it = count.begin();
1999 it != count.end();
2000 it++)
2001 {
2002 // Get pointer to element: Key
2003 FiniteElement* el_pt = it->first;
2004
2005 // Does it have to be split?
2006 if (it->second > 2)
2007 {
2008 elements_to_be_split.insert(el_pt);
2009 }
2010 }
2011
2012 // Vector for retained or newly built elements
2013 unsigned nel = this->nelement();
2014 Vector<FiniteElement*> new_or_retained_el_pt;
2015 new_or_retained_el_pt.reserve(nel);
2016
2017 // Map which returns the 4 newly created elements for each old corner
2018 // element
2019 std::map<FiniteElement*, Vector<FiniteElement*>>
2020 old_to_new_corner_element_map;
2021
2022 // Now loop over all elements
2023 for (unsigned e = 0; e < nel; e++)
2024 {
2025 // Get pointer to element
2026 FiniteElement* el_pt = this->finite_element_pt(e);
2027
2028 // Does it have to be split? I.e. is it in the set?
2029 std::set<FiniteElement*>::iterator it = std::find(
2030 elements_to_be_split.begin(), elements_to_be_split.end(), el_pt);
2031
2032 // It's not in the set, so iterator has reached end
2033 if (it == elements_to_be_split.end())
2034 {
2035 // Carry it across
2036 new_or_retained_el_pt.push_back(el_pt);
2037 }
2038 // It's in the set of elements to be split
2039 else
2040 {
2041 // Storage for new nodes -- Note: All new nodes are interior and
2042 // therefore cannot be boundary nodes!
2043 Node* node0_pt = 0;
2044 Node* node1_pt = 0;
2045 Node* node2_pt = 0;
2046 Node* node3_pt = 0;
2047 Node* node4_pt = 0;
2048 Node* node5_pt = 0;
2049 Node* node6_pt = 0;
2050 Node* node7_pt = 0;
2051 Node* node8_pt = 0;
2052 Node* node9_pt = 0;
2053 Node* node10_pt = 0;
2054
2055 // Create first new element
2056 FiniteElement* el1_pt = new ELEMENT;
2057
2058 // Copy existing nodes
2059 el1_pt->node_pt(0) = el_pt->node_pt(0);
2060 el1_pt->node_pt(1) = el_pt->node_pt(1);
2061 el1_pt->node_pt(3) = el_pt->node_pt(3);
2062 if (nnode_1d == 3)
2063 {
2064 el1_pt->node_pt(4) = el_pt->node_pt(4);
2065 el1_pt->node_pt(6) = el_pt->node_pt(6);
2066 el1_pt->node_pt(9) = el_pt->node_pt(9);
2067 }
2068
2069 // Create new nodes
2070 // If we have an enriched element then don't need to construct centroid
2071 // node
2072 if (nnode == 15)
2073 {
2074 node0_pt = el_pt->node_pt(14);
2075 el1_pt->node_pt(2) = node0_pt;
2076 node5_pt = el1_pt->construct_node(11, time_stepper_pt);
2077 node6_pt = el1_pt->construct_node(13, time_stepper_pt);
2078 node9_pt = el1_pt->construct_node(12, time_stepper_pt);
2079
2080 // Copy others over
2081 el1_pt->node_pt(10) = el_pt->node_pt(10);
2082 }
2083 // If not enriched we do
2084 else
2085 {
2086 node0_pt = el1_pt->construct_node(2, time_stepper_pt);
2087 }
2088 if (nnode_1d == 3)
2089 {
2090 node1_pt = el1_pt->construct_boundary_node(5, time_stepper_pt);
2091 node2_pt = el1_pt->construct_boundary_node(7, time_stepper_pt);
2092 node4_pt = el1_pt->construct_boundary_node(8, time_stepper_pt);
2093 }
2094
2095
2096 // Create second new element
2097 FiniteElement* el2_pt = new ELEMENT;
2098
2099 // Copy existing nodes
2100 el2_pt->node_pt(0) = el_pt->node_pt(0);
2101 el2_pt->node_pt(1) = el_pt->node_pt(1);
2102 el2_pt->node_pt(2) = el_pt->node_pt(2);
2103 if (nnode_1d == 3)
2104 {
2105 el2_pt->node_pt(4) = el_pt->node_pt(4);
2106 el2_pt->node_pt(5) = el_pt->node_pt(5);
2107 el2_pt->node_pt(7) = el_pt->node_pt(7);
2108 }
2109
2110 // Create new node
2111 if (nnode_1d == 3)
2112 {
2113 node3_pt = el2_pt->construct_boundary_node(8, time_stepper_pt);
2114 }
2115
2116 // Copy existing new nodes
2117 el2_pt->node_pt(3) = node0_pt;
2118 if (nnode_1d == 3)
2119 {
2120 el2_pt->node_pt(6) = node1_pt;
2121 el2_pt->node_pt(9) = node2_pt;
2122 }
2123
2124 // Copy and constuct other nodes for enriched elements
2125 if (nnode == 15)
2126 {
2127 el2_pt->node_pt(10) = node5_pt;
2128 el2_pt->node_pt(11) = el_pt->node_pt(11);
2129 node8_pt = el2_pt->construct_node(12, time_stepper_pt);
2130 node10_pt = el2_pt->construct_node(13, time_stepper_pt);
2131 }
2132
2133 // Create third new element
2134 FiniteElement* el3_pt = new ELEMENT;
2135
2136 // Copy existing nodes
2137 el3_pt->node_pt(1) = el_pt->node_pt(1);
2138 el3_pt->node_pt(2) = el_pt->node_pt(2);
2139 el3_pt->node_pt(3) = el_pt->node_pt(3);
2140 if (nnode_1d == 3)
2141 {
2142 el3_pt->node_pt(7) = el_pt->node_pt(7);
2143 el3_pt->node_pt(8) = el_pt->node_pt(8);
2144 el3_pt->node_pt(9) = el_pt->node_pt(9);
2145 }
2146
2147 // Copy existing new nodes
2148 el3_pt->node_pt(0) = node0_pt;
2149 if (nnode_1d == 3)
2150 {
2151 el3_pt->node_pt(4) = node2_pt;
2152 el3_pt->node_pt(5) = node3_pt;
2153 el3_pt->node_pt(6) = node4_pt;
2154 }
2155
2156 // Copy and constuct other nodes for enriched elements
2157 if (nnode == 15)
2158 {
2159 el3_pt->node_pt(10) = node6_pt;
2160 el3_pt->node_pt(11) = node10_pt;
2161 node7_pt = el3_pt->construct_node(12, time_stepper_pt);
2162 el3_pt->node_pt(13) = el_pt->node_pt(13);
2163 }
2164
2165
2166 // Create fourth new element
2167 FiniteElement* el4_pt = new ELEMENT;
2168
2169 // Copy existing nodes
2170 el4_pt->node_pt(0) = el_pt->node_pt(0);
2171 el4_pt->node_pt(2) = el_pt->node_pt(2);
2172 el4_pt->node_pt(3) = el_pt->node_pt(3);
2173 if (nnode_1d == 3)
2174 {
2175 el4_pt->node_pt(5) = el_pt->node_pt(5);
2176 el4_pt->node_pt(6) = el_pt->node_pt(6);
2177 el4_pt->node_pt(8) = el_pt->node_pt(8);
2178 }
2179
2180 // Copy existing new nodes
2181 el4_pt->node_pt(1) = node0_pt;
2182 if (nnode_1d == 3)
2183 {
2184 el4_pt->node_pt(4) = node1_pt;
2185 el4_pt->node_pt(7) = node3_pt;
2186 el4_pt->node_pt(9) = node4_pt;
2187 }
2188
2189 // Copy all other nodes
2190 if (nnode == 15)
2191 {
2192 el4_pt->node_pt(10) = node9_pt;
2193 el4_pt->node_pt(11) = node8_pt;
2194 el4_pt->node_pt(12) = el_pt->node_pt(12);
2195 el4_pt->node_pt(13) = node7_pt;
2196 ;
2197 }
2198
2199
2200 // Add new elements and nodes
2201 new_or_retained_el_pt.push_back(el1_pt);
2202 new_or_retained_el_pt.push_back(el2_pt);
2203 new_or_retained_el_pt.push_back(el3_pt);
2204 new_or_retained_el_pt.push_back(el4_pt);
2205
2206 // create a vector of the newly created elements
2207 Vector<FiniteElement*> temp_vec(4);
2208 temp_vec[0] = el1_pt;
2209 temp_vec[1] = el2_pt;
2210 temp_vec[2] = el3_pt;
2211 temp_vec[3] = el4_pt;
2212
2213 // add the vector to the map
2214 old_to_new_corner_element_map.insert(
2215 std::pair<FiniteElement*, Vector<FiniteElement*>>(el_pt, temp_vec));
2216
2217 if (nnode != 15)
2218 {
2219 this->add_node_pt(node0_pt);
2220 }
2221 this->add_node_pt(node1_pt);
2222 this->add_node_pt(node2_pt);
2223 this->add_node_pt(node3_pt);
2224 this->add_node_pt(node4_pt);
2225 if (nnode == 15)
2226 {
2227 this->add_node_pt(node5_pt);
2228 this->add_node_pt(node6_pt);
2229 this->add_node_pt(node7_pt);
2230 this->add_node_pt(node8_pt);
2231 this->add_node_pt(node9_pt);
2232 }
2233
2234 // Set nodal positions
2235 for (unsigned i = 0; i < 3; i++)
2236 {
2237 // Only bother to set centroid if does not already exist
2238 if (nnode != 15)
2239 {
2240 node0_pt->x(i) =
2241 0.25 * (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
2242 el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i));
2243 }
2244
2245 if (nnode_1d == 3)
2246 {
2247 node1_pt->x(i) = 0.5 * (el_pt->node_pt(0)->x(i) + node0_pt->x(i));
2248 node2_pt->x(i) = 0.5 * (el_pt->node_pt(1)->x(i) + node0_pt->x(i));
2249 node3_pt->x(i) = 0.5 * (el_pt->node_pt(2)->x(i) + node0_pt->x(i));
2250 node4_pt->x(i) = 0.5 * (el_pt->node_pt(3)->x(i) + node0_pt->x(i));
2251 }
2252 }
2253
2254
2255 // Construct the four interior nodes if needed
2256 // and add to the list
2257 if (nnode == 15)
2258 {
2259 // Set the positions of the newly created mid-face nodes
2260 // New Node 5 lies in the plane between original nodes 0 1 centroid
2261 for (unsigned i = 0; i < 3; ++i)
2262 {
2263 node5_pt->x(i) =
2264 (el_pt->node_pt(0)->x(i) + el_pt->node_pt(1)->x(i) +
2265 el_pt->node_pt(14)->x(i)) /
2266 3.0;
2267 }
2268
2269 // New Node 6 lies in the plane between original nodes 1 3 centroid
2270 for (unsigned i = 0; i < 3; ++i)
2271 {
2272 node6_pt->x(i) =
2273 (el_pt->node_pt(1)->x(i) + el_pt->node_pt(3)->x(i) +
2274 el_pt->node_pt(14)->x(i)) /
2275 3.0;
2276 }
2277
2278 // New Node 7 lies in the plane between original nodes 2 3 centroid
2279 for (unsigned i = 0; i < 3; ++i)
2280 {
2281 node7_pt->x(i) =
2282 (el_pt->node_pt(2)->x(i) + el_pt->node_pt(3)->x(i) +
2283 el_pt->node_pt(14)->x(i)) /
2284 3.0;
2285 }
2286
2287 // New Node 8 lies in the plane between original nodes 0 2 centroid
2288 for (unsigned i = 0; i < 3; ++i)
2289 {
2290 node8_pt->x(i) =
2291 (el_pt->node_pt(0)->x(i) + el_pt->node_pt(2)->x(i) +
2292 el_pt->node_pt(14)->x(i)) /
2293 3.0;
2294 }
2295
2296 // New Node 9 lies in the plane between original nodes 0 3 centroid
2297 for (unsigned i = 0; i < 3; ++i)
2298 {
2299 node9_pt->x(i) =
2300 (el_pt->node_pt(0)->x(i) + el_pt->node_pt(3)->x(i) +
2301 el_pt->node_pt(14)->x(i)) /
2302 3.0;
2303 }
2304
2305 // New Node 10 lies in the plane between original nodes 1 2 centroid
2306 for (unsigned i = 0; i < 3; ++i)
2307 {
2308 node10_pt->x(i) =
2309 (el_pt->node_pt(1)->x(i) + el_pt->node_pt(2)->x(i) +
2310 el_pt->node_pt(14)->x(i)) /
2311 3.0;
2312 }
2313
2314 // Now create the new centroid nodes
2315
2316 // First element
2317 Node* temp_nod_pt = el1_pt->construct_node(14, time_stepper_pt);
2318 for (unsigned i = 0; i < 3; ++i)
2319 {
2320 double av_pos = 0.0;
2321 for (unsigned j = 0; j < 4; j++)
2322 {
2323 av_pos += el1_pt->node_pt(j)->x(i);
2324 }
2325
2326 temp_nod_pt->x(i) = 0.25 * av_pos;
2327 }
2328
2329 this->add_node_pt(temp_nod_pt);
2330
2331 // Second element
2332 temp_nod_pt = el2_pt->construct_node(14, time_stepper_pt);
2333 for (unsigned i = 0; i < 3; ++i)
2334 {
2335 double av_pos = 0.0;
2336 for (unsigned j = 0; j < 4; j++)
2337 {
2338 av_pos += el2_pt->node_pt(j)->x(i);
2339 }
2340 temp_nod_pt->x(i) = 0.25 * av_pos;
2341 }
2342 this->add_node_pt(temp_nod_pt);
2343
2344 // Third element
2345 temp_nod_pt = el3_pt->construct_node(14, time_stepper_pt);
2346 for (unsigned i = 0; i < 3; ++i)
2347 {
2348 double av_pos = 0.0;
2349 for (unsigned j = 0; j < 4; j++)
2350 {
2351 av_pos += el3_pt->node_pt(j)->x(i);
2352 }
2353 temp_nod_pt->x(i) = 0.25 * av_pos;
2354 }
2355 this->add_node_pt(temp_nod_pt);
2356
2357 // Fourth element
2358 temp_nod_pt = el4_pt->construct_node(14, time_stepper_pt);
2359 for (unsigned i = 0; i < 3; ++i)
2360 {
2361 double av_pos = 0.0;
2362 for (unsigned j = 0; j < 4; j++)
2363 {
2364 av_pos += el4_pt->node_pt(j)->x(i);
2365 }
2366 temp_nod_pt->x(i) = 0.25 * av_pos;
2367 }
2368 this->add_node_pt(temp_nod_pt);
2369 }
2370
2371 // Kill the old element
2372 delete el_pt;
2373 }
2374 }
2375
2376 // Flush element storage
2377 Element_pt.clear();
2378
2379 // Copy across
2380 nel = new_or_retained_el_pt.size();
2381 Element_pt.resize(nel);
2382 for (unsigned e = 0; e < nel; e++)
2383 {
2384 Element_pt[e] = new_or_retained_el_pt[e];
2385 }
2386
2387 // Setup boundary lookup scheme again
2389
2390 // -------------------------------------------------------------------------
2391 // The various boundary/region lookups now need updating to account for the
2392 // newly added/removed elements. This will be done in two stages:
2393 // Step 1: Add the new elements to the vector of elements in the same region
2394 // as the original corner element, and then delete the originals.
2395 // Updating this lookup makes things easier in the following step.
2396 // Step 2: Regenerate the two more specific lookups: One which gives the
2397 // elements on a given boundary in a given region, and the other
2398 // which maps elements on a given boundary in a given region to the
2399 // element's face index on that boundary.
2400 //
2401 // N.B. the lookup Triangular_facet_vertex_boundary_coordinate is setup in
2402 // the call to setup_boundary_element_info() above so doesn't need
2403 // additional work.
2404
2405 // if we have no regions then we have no lookups to update so we're done
2406 // here
2407 if (Region_attribute.size() == 0)
2408 {
2409 return;
2410 }
2411 // if we haven't had to split any corner elements then don't need to fiddle
2412 // with the lookups
2413 if (old_to_new_corner_element_map.size() == 0)
2414 {
2415 oomph_info << "\nNo corner elements need splitting\n\n";
2416 return;
2417 }
2418
2419 // ------------------------------------------
2420 // Step 1: update the region element lookup
2421
2422 // loop over the map of old corner elements which have been split
2423 for (std::map<FiniteElement*, Vector<FiniteElement*>>::iterator map_it =
2424 old_to_new_corner_element_map.begin();
2425 map_it != old_to_new_corner_element_map.end();
2426 map_it++)
2427 {
2428 // extract the old and new elements from the map
2429 FiniteElement* original_el_pt = map_it->first;
2430 Vector<FiniteElement*> new_el_pt = map_it->second;
2431
2432 unsigned original_region_index = 0;
2433
2434#ifdef PARANOID
2435 // flag for paranoia, if for some reason we don't find the original corner
2436 // element in any of the regions
2437 bool found = false;
2438#endif
2439
2440 Vector<FiniteElement*>::iterator region_element_it;
2441
2442 // loop over the regions and look for this original corner element to find
2443 // out which region it used to be in, so we can add the new elements to
2444 // the same region.
2445 for (unsigned region_index = 0; region_index < Region_element_pt.size();
2446 region_index++)
2447 {
2448 // for each region, search the vector of elements in this region for the
2449 // original corner element
2450 region_element_it = std::find(Region_element_pt[region_index].begin(),
2451 Region_element_pt[region_index].end(),
2452 original_el_pt);
2453
2454 // if the iterator hasn't reached the end then we've found it
2455 if (region_element_it != Region_element_pt[region_index].end())
2456 {
2457 // save the region index we're at
2458 original_region_index = region_index;
2459
2460#ifdef PARANOID
2461 // update the paranoid flag
2462 found = true;
2463#endif
2464
2465 // regions can't overlap, so once we've found one we're done
2466 break;
2467 }
2468 }
2469
2470#ifdef PARANOID
2471 if (!found)
2472 {
2473 std::ostringstream error_message;
2474 error_message
2475 << "The corner element being split does not appear to be in any "
2476 << "region, so something has gone wrong with the region lookup "
2477 "scheme\n";
2478
2479 throw OomphLibError(error_message.str(),
2480 OOMPH_CURRENT_FUNCTION,
2481 OOMPH_EXCEPTION_LOCATION);
2482 }
2483#endif
2484
2485 // Now update the basic region lookup. The iterator will still point to
2486 // the original corner element in the lookup, so we can delete this easily
2487 Region_element_pt[original_region_index].erase(region_element_it);
2488 for (unsigned i = 0; i < 4; i++)
2489 {
2490 Region_element_pt[original_region_index].push_back(new_el_pt[i]);
2491 }
2492 }
2493 // ------------------------------------------
2494 // Step 2: Clear and regenerate lookups
2495
2498
2501
2502 for (unsigned b = 0; b < nboundary(); b++)
2503 {
2504 // Loop over elements next to that boundary
2505 unsigned nel = this->nboundary_element(b);
2506 for (unsigned e = 0; e < nel; e++)
2507 {
2508 FiniteElement* el_pt = boundary_element_pt(b, e);
2509
2510 // now search for it in each region
2511 for (unsigned r_index = 0; r_index < Region_attribute.size(); r_index++)
2512 {
2513 unsigned region_id = static_cast<unsigned>(Region_attribute[r_index]);
2514
2516 std::find(Region_element_pt[r_index].begin(),
2517 Region_element_pt[r_index].end(),
2518 el_pt);
2519
2520 // if we find this element in the current region, then update our
2521 // lookups
2522 if (it != Region_element_pt[r_index].end())
2523 {
2524 Boundary_region_element_pt[b][region_id].push_back(el_pt);
2525
2526 unsigned face_index = face_index_at_boundary(b, e);
2527 Face_index_region_at_boundary[b][region_id].push_back(face_index);
2528 }
2529 }
2530 }
2531 }
2532
2533 oomph_info << "\nNumber of outer corner elements split: "
2534 << old_to_new_corner_element_map.size() << "\n\n";
2535
2536 } // end split_elements_in_corners()
2537} // namespace oomph
2538
2539#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh....
Definition: nodes.h:2242
Base class for upgraded disk-like GeomObject (i.e. 2D surface in 3D space) with specification of boun...
Information for documentation of results: Directory and file number to enable output in the form RESL...
std::string & label()
String used (e.g.) for labeling output files.
bool is_doc_enabled() const
Are we documenting?
void disable_doc()
Disable documentation.
std::string directory() const
Output directory.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5014
A general Finite Element class.
Definition: elements.h:1313
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition: elements.h:2509
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition: elements.h:2538
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5119
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
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
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
bool Lookup_for_elements_next_boundary_is_setup
Flag to indicate that the lookup schemes for elements that are adjacent to the boundaries has been se...
Definition: mesh.h:87
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:190
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
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
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
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
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2394
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....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: tet_mesh.h:661
virtual ~TetMeshBase()
Destructor (empty)
Definition: tet_mesh.h:673
void snap_nodes_onto_geometric_objects()
Move the nodes on boundaries with associated GeomObjects so that they exactly coincide with the geome...
Definition: tet_mesh.cc:483
void snap_to_quadratic_surface(const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal, DocInfo &doc_info)
Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface,...
Definition: tet_mesh.h:1459
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
Vector< double > Region_attribute
Vector of attributes associated with the elements in each region NOTE: double is enforced on us by te...
Definition: tet_mesh.h:1027
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
void setup_boundary_coordinates(const unsigned &b, const bool &switch_normal)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Definition: tet_mesh.h:733
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
TimeStepper * Time_stepper_pt
Timestepper used to build nodes.
Definition: tet_mesh.h:1057
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
std::map< unsigned, TetMeshFacetedSurface * > Tet_mesh_faceted_surface_pt
Reverse lookup scheme: Pointer to faceted surface (if any!) associated with boundary b.
Definition: tet_mesh.h:1045
Vector< TetMeshFacetedSurface * > Internal_surface_pt
Vector to faceted surfaces that define internal boundaries.
Definition: tet_mesh.h:1041
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
Storage for elements adjacent to a boundary in a particular region.
Definition: tet_mesh.h:1031
void operator=(const TetMeshBase &)=delete
Broken assignment operator.
std::map< unsigned, TetMeshFacet * > Tet_mesh_facet_pt
Reverse lookup scheme: Pointer to facet (if any!) associated with boundary b.
Definition: tet_mesh.h:1049
void setup_boundary_coordinates(const unsigned &b)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Definition: tet_mesh.h:702
TetMeshFacetedClosedSurface * Outer_boundary_pt
Faceted surface that defines outer boundaries.
Definition: tet_mesh.h:1038
static double Tolerance_for_boundary_finding
Global static data that specifies the permitted error in the setup of the boundary coordinates.
Definition: tet_mesh.h:677
void assess_mesh_quality(std::ofstream &some_file)
Assess mesh quality: Ratio of max. edge length to min. height, so if it's very large it's BAAAAAD.
Definition: tet_mesh.cc:373
TetMeshBase()
Constructor.
Definition: tet_mesh.h:664
void snap_to_quadratic_surface(const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal)
Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface,...
Definition: tet_mesh.h:979
void setup_boundary_element_info()
Setup lookup schemes which establish which elements are located next to mesh's boundaries (wrapper to...
Definition: tet_mesh.h:1007
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Storage for the face index adjacent to a boundary in a particular region.
Definition: tet_mesh.h:1035
std::map< unsigned, Vector< Vector< double > > > Triangular_facet_vertex_boundary_coordinate
Boundary coordinates of vertices in triangular facets associated with given boundary....
Definition: tet_mesh.h:1054
TetMeshBase(const TetMeshBase &node)=delete
Broken copy constructor.
Vector< Vector< FiniteElement * > > Region_element_pt
Vectors of vectors of elements in each region (note: this just stores them; the region IDs are contai...
Definition: tet_mesh.h:1022
void split_elements_in_corners(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Non-Delaunay split elements that have three faces on a boundary into sons. Timestepper species timest...
Definition: tet_mesh.h:1950
void setup_boundary_coordinates(const unsigned &b, std::ofstream &outfile)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Definition: tet_mesh.h:788
unsigned nregion_element(const unsigned &r)
Return the number of elements in region r.
Definition: tet_mesh.h:866
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:168
Vector< TetMeshVertex * > Vertex_pt
Pointer to vertices.
Definition: tet_mesh.h:280
void set_vertex_pt(const unsigned &j, TetMeshVertex *vertex_pt)
Set pointer to j-th vertex and pass one-based boundary id that may already have been set for this fac...
Definition: tet_mesh.h:194
void output(std::ostream &outfile) const
Output.
Definition: tet_mesh.h:263
void set_one_based_region_that_facet_is_embedded_in(const unsigned &one_based_region_id)
Facet is to be embedded in specified one-based region.
Definition: tet_mesh.h:249
std::set< unsigned > One_based_adjacent_region_id
Set of one-based adjacent region ids; no adjacent region if it's zero.
Definition: tet_mesh.h:287
bool facet_is_embedded_in_a_specified_region()
Boolean indicating that facet is embedded in a specified region.
Definition: tet_mesh.h:243
void set_one_based_boundary_id(const unsigned &one_based_id)
Set (one-based!) boundary IDs this facet lives on. Passed to any existing vertices and to any future ...
Definition: tet_mesh.h:214
unsigned one_based_boundary_id() const
Constant access to (one-based!) boundary IDs this facet lives on.
Definition: tet_mesh.h:207
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:181
unsigned One_based_boundary_id
(One-based!) boundary IDs this facet lives on
Definition: tet_mesh.h:283
TetMeshVertex * vertex_pt(const unsigned &j) const
Pointer to j-th vertex (const)
Definition: tet_mesh.h:187
unsigned One_based_region_id_that_facet_is_embedded_in
Facet is to be embedded in specified one-based region. Defaults to zero, indicating that its not embe...
Definition: tet_mesh.h:292
void set_one_based_adjacent_region_id(const unsigned &one_based_id)
Set (one-based!) region ID this facet is adjacent to. Specification of zero argument is harmless as i...
Definition: tet_mesh.h:231
unsigned one_based_region_that_facet_is_embedded_in()
Which (one-based) region is facet embedded in (if zero, it's not embedded in any region)
Definition: tet_mesh.h:257
TetMeshFacet(const unsigned &nvertex)
Constructor: Specify number of vertices.
Definition: tet_mesh.h:171
std::set< unsigned > one_based_adjacent_region_id() const
Return set of (one-based!) region IDs this facet is adjacent to.
Definition: tet_mesh.h:237
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: tet_mesh.h:634
virtual ~TetMeshFacetedClosedSurfaceForRemesh()
Destructor. Delete allocated memory.
Definition: tet_mesh.cc:73
TetMeshFacetedClosedSurfaceForRemesh(Vector< Node * > const &vertex_node_pt, Vector< Vector< unsigned > > const &facet_connectivity, Vector< unsigned > const &facet_boundary_id)
Constructor for a FacetedSurface created from a list of nodes and connectivity information....
Definition: tet_mesh.cc:40
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:538
TetMeshFacetedClosedSurface()
Constructor:
Definition: tet_mesh.h:541
bool Faceted_volume_represents_hole_for_gmsh
Does closed surface represent hole for gmsh?
Definition: tet_mesh.h:624
void enable_faceted_volume_represents_hole_for_gmsh()
Declare closed surface to represent hole for gmsh.
Definition: tet_mesh.h:550
bool faceted_volume_represents_hole_for_gmsh() const
Does closed surface represent hole for gmsh?
Definition: tet_mesh.h:562
Vector< std::pair< Vector< double >, int > > Internal_point_for_tetgen
Storage for internal points for tetgen. Stores pair of: – Vector containing coordinates of internal p...
Definition: tet_mesh.h:621
virtual ~TetMeshFacetedClosedSurface()
Empty destructor.
Definition: tet_mesh.h:547
unsigned ninternal_point_for_tetgen()
Number of internal points (identifying either regions or holes) for tetgen.
Definition: tet_mesh.h:591
bool internal_point_identifies_region_for_tetgen(const unsigned &j)
Is j-th internal point for tetgen associated with a region?
Definition: tet_mesh.h:611
const int & region_id_for_tetgen(const unsigned &j) const
Return the (zero-based) region ID of j-th internal point for tetgen. Negative if it's actually a hole...
Definition: tet_mesh.h:598
bool internal_point_identifies_hole_for_tetgen(const unsigned &j)
Is j-th internal point for tetgen associated with a hole?
Definition: tet_mesh.h:605
void disable_faceted_volume_represents_hole_for_gmsh()
Declare closed surface NOT to represent hole for gmsh.
Definition: tet_mesh.h:556
void set_region_for_tetgen(const unsigned &region_id, const Vector< double > &region_point)
Specify a region; pass (zero-based) region ID and coordinate of point in region for tetgen.
Definition: tet_mesh.h:582
void set_hole_for_tetgen(const Vector< double > &hole_point)
Specify coordinate of hole for tetgen.
Definition: tet_mesh.h:575
const double & internal_point_for_tetgen(const unsigned &j, const unsigned &i) const
i=th coordinate of the j-th internal point for tetgen
Definition: tet_mesh.h:568
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:306
TetMeshFacet * facet_pt(const unsigned &j) const
Pointer to j-th facet.
Definition: tet_mesh.h:373
unsigned nvertex_on_facet(const unsigned &j) const
Number of vertices defining the j-th facet.
Definition: tet_mesh.h:349
bool Boundaries_can_be_split_in_tetgen
Boolean to indicate whether extra vertices can be added on the boundary in tetgen.
Definition: tet_mesh.h:490
TetMeshVertex * vertex_pt(const unsigned &j) const
Pointer to j-th vertex.
Definition: tet_mesh.h:379
void disable_boundaries_can_be_split_in_tetgen()
Test whether boundaries can be split in tetgen.
Definition: tet_mesh.h:367
virtual void boundary_zeta01(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:417
Vector< TetMeshVertex * > Vertex_pt
Vector pointers to vertices.
Definition: tet_mesh.h:483
unsigned nfacet() const
Number of facets.
Definition: tet_mesh.h:325
Vector< Vector< unsigned > > Facet_vertex_index_in_tetgen
Facet connectivity: Facet_vertex_index[f][j] is the index of the j-th vertex (in the Vertex_pt vector...
Definition: tet_mesh.h:494
DiskLikeGeomObjectWithBoundaries * Geom_object_with_boundaries_pt
GeomObject with boundaries associated with this surface.
Definition: tet_mesh.h:497
Vector< unsigned > vertex_index_in_tetgen(const unsigned &f)
Facet connectivity: vertex_index[j] is the index of the j-th vertex (in the Vertex_pt vector) in face...
Definition: tet_mesh.h:472
virtual void boundary_zeta12(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:436
bool boundaries_can_be_split_in_tetgen()
Test whether boundary can be split in tetgen.
Definition: tet_mesh.h:355
Vector< TetMeshFacet * > Facet_pt
Vector of pointers to facets.
Definition: tet_mesh.h:486
unsigned one_based_facet_boundary_id(const unsigned &j) const
One-based boundary id of j-th facet.
Definition: tet_mesh.h:331
void enable_boundaries_can_be_split_in_tetgen()
Test whether boundaries can be split in tetgen.
Definition: tet_mesh.h:361
DiskLikeGeomObjectWithBoundaries * geom_object_with_boundaries_pt()
Access to GeomObject with boundaries associated with this surface (Null if there isn't one!...
Definition: tet_mesh.h:386
virtual void boundary_zeta20(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:455
void setup_facet_connectivity_for_tetgen()
Setup facet connectivity for tetgen.
Definition: tet_mesh.h:502
double vertex_coordinate(const unsigned &j, const unsigned &i) const
i-th coordinate of j-th vertex
Definition: tet_mesh.h:343
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:319
TetMeshFacetedSurface()
Constructor:
Definition: tet_mesh.h:309
void output(const std::string &filename) const
Output.
Definition: tet_mesh.h:403
void output(std::ostream &outfile) const
Output.
Definition: tet_mesh.h:392
virtual ~TetMeshFacetedSurface()
Empty destructor.
Definition: tet_mesh.h:316
unsigned one_based_vertex_boundary_id(const unsigned &j) const
First (of possibly multiple) one-based boundary id of j-th vertex.
Definition: tet_mesh.h:337
Vertex for Tet mesh generation. Can lie on multiple boundaries (identified via one-based enumeration!...
Definition: tet_mesh.h:50
void set_zeta_in_geom_object(const Vector< double > &zeta)
Set intrinisic coordinates in GeomObject.
Definition: tet_mesh.h:97
Vector< double > zeta_in_geom_object() const
Get intrinisic coordinates in GeomObject (returns zero sized vector if no such coordinates have been ...
Definition: tet_mesh.h:118
Vector< double > X
Coordinate vector.
Definition: tet_mesh.h:147
TetMeshVertex(Node *const &node_pt)
Definition: tet_mesh.h:72
std::set< unsigned > One_based_boundary_id
Set of (one-based!) boundary IDs this vertex lives on.
Definition: tet_mesh.h:150
unsigned one_based_boundary_id() const
First (of possibly multiple) one-based boundary id.
Definition: tet_mesh.h:130
Vector< double > Zeta_in_geom_object
Intrinisic coordinates in GeomObject (zero sized if there isn't one.
Definition: tet_mesh.h:154
void set_one_based_boundary_id(const unsigned &id)
Set of (one-based!) boundary IDs this vertex lives on.
Definition: tet_mesh.h:141
double x(const unsigned &i) const
i-th coordinate
Definition: tet_mesh.h:124
TetMeshVertex(const Vector< double > &x)
Constructor: Pass coordinates (length 3!)
Definition: tet_mesh.h:56
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...