gmsh_tet_mesh.template.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
27#ifndef OOMPH_GMSH_TET_MESH_HEADER
28#define OOMPH_GMSH_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#include <algorithm>
36#include <iterator>
37
38// OOMPH-LIB Headers
40#include "../generic/sample_point_parameters.h"
41#include "../generic/mesh_as_geometric_object.h"
42#include "../generic/projection.h"
43
44namespace oomph
45{
46 //=========================================================================
47 /// Class to collate parameters for Gmsh mesh generation
48 //=========================================================================
50 {
51 public:
52 /// Specify outer boundary of domain to be meshed.
53 /// Other parameters get default values and can be set via
54 /// member functions
58 Element_volume(1.0),
59 Target_size_file_name(".gmsh_target_size_for_adaptation.dat"),
60 Geo_and_msh_file_stem(".geo_file"),
63 Min_element_size(0.01),
67 5),
73 {
74 }
75
76 /// Outer boundary
78 {
79 return Outer_boundary_pt;
80 }
81
82 /// Internal boundaries
84 {
86 }
87
88 /// Uniform target element volume
90 {
91 return Element_volume;
92 }
93
94 /// Filename for target volumes (for system-call based transfer to
95 /// gmsh during mesh adaptation). Default:
96 /// .gmsh_target_size_for_adaptation.dat
98 {
100 }
101
102 /// String to be issued via system command to activate gmsh
104 {
106 }
107
108 /// Stem for geo and msh files (input/output to command-line gmsh
109 /// invocation)
111 {
113 }
114
115
116 /// Max. element size during refinement
118 {
119 return Max_element_size;
120 }
121
122 /// Min. element size during refinement
124 {
125 return Min_element_size;
126 }
127
128 /// Max. permitted edge ratio
130 {
132 }
133
134 /// (Isotropic) grid spacing for target size transfer
136 {
138 }
139
140 /// Target size is transferred onto regular grid (for gmsh) by
141 /// locate zeta. We try to find the exact point in the existing
142 /// mesh but if we fail to converge from the nearest specified number
143 /// of sample points we use the nearest of those.
145 {
147 }
148
149 /// Stem for filename used to doc target element sizes on
150 /// gmsh grid. No doc if stem is equal to empty string (or counter
151 /// is negative)
153 {
155 }
156
157 /// Counter for filename used to doc target element sizes on
158 /// gmsh grid. No doc if stem is equal to empty string (or counter
159 /// is negative)
161 {
163 }
164
165 /// Is projection of old solution onto new mesh disabled?
167 {
169 }
170
171 /// Disable projection of old solution onto new mesh
173 {
175 }
176
177 /// Disable projection of old solution onto new mesh
179 {
181 }
182
183 /// Output filename for gmsh on-screen output
185 {
187 }
188
189 /// Counter for marker that indicates where we are
190 /// in gmsh on-screen output
192 {
194 }
195
196 private:
197 /// Pointer to outer boundary
199
200 /// Internal boundaries
202
203 /// Uniform element volume
205
206 /// Filename for target volume (for system-call based transfer
207 /// to gmsh during mesh adaptation)
209
210 /// Stem for geo and msh files (input/output to command-line gmsh
211 /// invocation)
213
214 /// Gmsh command line invocation string
216
217 /// Max. element size during refinement
219
220 /// Min. element size during refinement
222
223 /// Max edge ratio before remesh gets triggered
225
226 /// (Isotropic) grid spacing for target size transfer
228
229 /// Target size is transferred onto regular grid (for gmsh) by
230 /// locate zeta. We try to find the exact point in the existing
231 /// mesh but if we fail to converge from the nearest specified number of
232 /// sample points we use the nearest of those.
233 unsigned
235
236 /// Stem for filename used to doc target element sizes on
237 /// gmsh grid. No doc if stem is equal to empty string (or counter
238 /// is negative)
240
241 /// Counter for filename used to doc target element sizes on
242 /// gmsh grid. No doc if stem is equal to empty string (or counter
243 /// is negative)
245
246 /// Is projection of old solution onto new mesh disabled?
248
249 /// Output filename for gmsh on-screen output
251
252 /// Counter for marker that indicates where we are
253 /// in gmsh on-screen output
255 };
256
257 /// //////////////////////////////////////////////////////////////////
258 /// //////////////////////////////////////////////////////////////////
259 /// //////////////////////////////////////////////////////////////////
260
261
262 //=========================================================================
263 /// Helper class to keep track of edges in tet mesh generation
264 //=========================================================================
266 {
267 public:
268 /// Constructor: Pass two vertices, identified by their indices
269 /// Edge "direction" is from lower vertex to higher vertex id so
270 /// can compare if we're dealing with the same one...
271 TetEdge(const unsigned& vertex1, const unsigned& vertex2)
272 {
273 if (vertex1 > vertex2)
274 {
275 Reversed = true;
276 Vertex_pair = std::make_pair(vertex2, vertex1);
277 }
278 else if (vertex1 < vertex2)
279 {
280 Reversed = false;
281 Vertex_pair = std::make_pair(vertex1, vertex2);
282 }
283 else
284 {
285 throw OomphLibError(
286 "Muppet! You can't build an edge from one vertex to the same vertex!",
287 OOMPH_CURRENT_FUNCTION,
288 OOMPH_EXCEPTION_LOCATION);
289 }
290 }
291
292
293 /// First vertex id
294 unsigned first_vertex_id() const
295 {
296 return Vertex_pair.first;
297 }
298
299 /// Second vertex id
300 unsigned second_vertex_id() const
301 {
302 return Vertex_pair.second;
303 }
304
305
306 /// Edge is reversed in the sense that vertex1 actually has a higher
307 /// id than vertex2 (when specified in the constructor)
308 bool is_reversed() const
309 {
310 return Reversed;
311 }
312
313 /// Comparison operator: Edges are identical if their sorted
314 /// (and therefore possibly reversed) vertex ids agree
315 bool operator==(const TetEdge& tet_edge) const
316 {
317 return ((tet_edge.first_vertex_id() == Vertex_pair.first) &&
318 (tet_edge.second_vertex_id() == Vertex_pair.second));
319 }
320
321 /// Comparison operator. Lexicographic comparison based on
322 /// vertex ids
323 bool operator<(const TetEdge& tet_edge) const
324 {
325 if ((tet_edge.first_vertex_id() == Vertex_pair.first) &&
326 (tet_edge.second_vertex_id() == Vertex_pair.second))
327 {
328 return false;
329 }
330 else
331 {
332 if (tet_edge.first_vertex_id() == Vertex_pair.first)
333 {
334 return (tet_edge.second_vertex_id() < Vertex_pair.second);
335 }
336 else
337 {
338 return (tet_edge.first_vertex_id() < Vertex_pair.first);
339 }
340 }
341 }
342
343
344 private:
345 /// The vertices (sorted by vertex ids)
346 std::pair<unsigned, unsigned> Vertex_pair;
347
348 /// Is it reversed? I.e. is the first input vertex stored after
349 /// the second one?
351 };
352
353
354 /// ///////////////////////////////////////////////////////////////////////
355 /// ///////////////////////////////////////////////////////////////////////
356 /// ///////////////////////////////////////////////////////////////////////
357
358
359 /// Forward declaration
360 template<class ELEMENT>
361 class GmshTetMesh;
362
363
364 //=========================================================================
365 // Gmsh-based Tet scaffold mesh
366 //=========================================================================
367 class GmshTetScaffoldMesh : public virtual TetMeshBase, public virtual Mesh
368 {
369 public:
370 /// We're friends with the actual mesh
371 template<class ELEMENT>
372 friend class GmshTetMesh;
373
374 /// Build mesh, based on specified parameters. If boolean is set
375 /// to true, the target element sizes are read from file (used during
376 /// adaptation; otherwise uniform target size is used).
378 const bool& use_mesh_grading_from_file)
379 : Gmsh_parameters_pt(gmsh_parameters_pt)
380 {
381 double t_start = TimingHelpers::timer();
382
383 // Create .geo file
384 write_geo_file(use_mesh_grading_from_file);
385
386 oomph_info << "Time for writing geo file : "
387 << TimingHelpers::timer() - t_start << " sec " << std::endl;
388 t_start = TimingHelpers::timer();
389
390 // Execute gmsh on command line
391 std::string gmsh_command_line_string = "";
392
393 std::string gmsh_onscreen_output_file_name =
394 gmsh_parameters_pt->gmsh_onscreen_output_file_name();
395 std::ofstream gmsh_on_screen_output_file;
396 std::stringstream marker;
397 if (gmsh_onscreen_output_file_name != "")
398 {
399 marker << "\n\n====================================================\n"
400 << " gmsh invocation: "
401 << gmsh_parameters_pt->gmsh_onscreen_output_counter()
402 << std::endl
403 << "====================================================\n\n\n"
404 << std::endl;
405 gmsh_parameters_pt->gmsh_onscreen_output_counter()++;
406 oomph_info << marker.str();
407 gmsh_on_screen_output_file.open(gmsh_onscreen_output_file_name.c_str(),
408 std::ofstream::app);
409 gmsh_on_screen_output_file << marker.str();
410 gmsh_on_screen_output_file.close();
411 }
412
413 gmsh_command_line_string +=
416 if (gmsh_onscreen_output_file_name != "")
417 {
418 gmsh_command_line_string += " >> " + gmsh_onscreen_output_file_name;
419 }
420
421 // Note return flag isn't particularly well defined but we report it
422 // anyway to aid detection of problems...
423 int return_flag = system(gmsh_command_line_string.c_str());
424 oomph_info << "fyi: return from system command: " << return_flag
425 << std::endl;
426 oomph_info << "Time for gmsh system call : "
427 << TimingHelpers::timer() - t_start << " sec " << std::endl;
428 t_start = TimingHelpers::timer();
429
430
431 // Create the mesh
433
434 oomph_info << "Time for creating mesh from msh file: "
435 << TimingHelpers::timer() - t_start << " sec " << std::endl;
436 }
437
438 private:
439 /// Create mesh from msh file (created internally via disk-based operations)
441 {
442 // Create filename from stem
443 std::string mesh_file_name =
445
446 // Keep around if we ever write a version where name is passed in.
447 /* // Check extension */
448 /* #ifdef PARANOID */
449 /* std::string mesh_file_stem= */
450 /* mesh_file_name.substr(0,mesh_file_name.length()-4); */
451 /* std::string test=mesh_file_stem+".msh"; */
452 /* if (test!=mesh_file_name) */
453 /* { */
454 /* std::string error_msg("msh file has wrong extension: "); */
455 /* error_msg += " " + mesh_file_name; */
456 /* throw OomphLibError(error_msg, */
457 /* OOMPH_CURRENT_FUNCTION, */
458 /* OOMPH_EXCEPTION_LOCATION); */
459 /* } */
460 /* #endif */
461
462 // Open wide...
463 std::ifstream mesh_file(mesh_file_name.c_str(), std::ios_base::in);
464
465// Check that the file actually opened correctly
466#ifdef PARANOID
467 if (!mesh_file.is_open())
468 {
469 std::string error_msg("Failed to open mesh file: ");
470 error_msg += "\"" + mesh_file_name + "\".";
471 throw OomphLibError(
472 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
473 }
474#endif
475
476 // First line: Must be "$MeshFormat"
477 //----------------------------------
478 std::string line;
479 mesh_file >> line;
480 if (line != "$MeshFormat")
481 {
482 std::string error_msg(
483 "First line has to contain the string \"$MeshFormat\"; ");
484 error_msg += " yours contains: " + line;
485 throw OomphLibError(
486 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
487 }
488
489 // Rest of line
490 mesh_file >> line;
491
492
493 // Nodes
494 //------
495
496 // Now keep reading until we find the nodes
497 bool keep_going = true;
498 while (keep_going)
499 {
500 std::string line;
501 mesh_file >> line;
502 if (line == "$Nodes")
503 {
504 keep_going = false;
505 }
506 }
507 unsigned nnod = 0;
508 mesh_file >> nnod;
509 std::map<unsigned, Vector<double>> node_coordinate;
510 for (unsigned j = 0; j < nnod; j++)
511 {
512 unsigned node_number = 0;
513 mesh_file >> node_number;
514
515 // Read rest of line word by word
517 std::getline(mesh_file, s);
518 std::istringstream iss(s);
519 std::string sub;
520 while (iss >> sub)
521 {
522 node_coordinate[node_number].push_back(atof(sub.c_str()));
523 }
524 }
525
526 mesh_file >> line;
527 if (line != "$EndNodes")
528 {
529 std::string error_msg("Line has to contain the string \"$EndNodes\"; ");
530 error_msg += " yours contains: " + line;
531 throw OomphLibError(
532 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
533 }
534
535 // Rewind
536 mesh_file.clear();
537 mesh_file.seekg(0);
538
539 mesh_file >> line;
540 if (line != "$MeshFormat")
541 {
542 std::string error_msg(
543 "First line has to contain the string \"$MeshFormat\"; ");
544 error_msg += " yours contains: " + line;
545 throw OomphLibError(
546 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
547 }
548
549 // Rest of line
550 mesh_file >> line;
551
552 // Storage for boundaries the nodes are on
553 std::map<unsigned, std::set<unsigned>> one_based_boundaries_of_node;
554
555 // Elements
556 //---------
557
558 // node number = boundary_node[one_based_bound_id][...]
559 std::map<unsigned, std::set<unsigned>> boundary_node;
560
561 // element number = region_element[one_based_region_id][...]
562 std::map<unsigned, std::set<unsigned>> region_element;
563
564 // Now keep reading until we find the nodes
565 keep_going = true;
566 while (keep_going)
567 {
568 std::string line;
569 mesh_file >> line;
570 if (line == "$Elements")
571 {
572 keep_going = false;
573 }
574 }
575
576
577 unsigned highest_one_based_boundary_id = 0;
578 unsigned n_tet_el = 0;
579
580 unsigned nel = 0;
581 mesh_file >> nel;
582 std::map<unsigned, Vector<unsigned>> element_node_index;
583 for (unsigned e = 0; e < nel; e++)
584 {
585 // For each element the msh file provides:
586 //
587 // elm-number elm-type number-of-tags < tag > ... node-number-list
588 //
589 // By default, the first tag is the number of the physical entity to
590 // which the element belongs; the second is the number of the elementary
591 // geometrical entity to which the element belongs; the third is the
592 // number of mesh partitions to which the element belongs, followed by
593 // the partition ids (negative partition ids indicate ghost cells). A
594 // zero tag is equivalent to no tag. Gmsh and most codes using the MSH 2
595 // format require at least the first two tags (physical and elementary
596 // tags).
597 unsigned el_number = 0;
598 mesh_file >> el_number;
599
600
601 unsigned el_type = 0;
602 mesh_file >> el_type;
603
604 switch (el_type)
605 {
606 case 1:
607 // oomph_info << "Line element" << std::endl;
608 break;
609
610 case 2:
611 // oomph_info << "Triangle" << std::endl;
612 break;
613
614 case 4:
615 n_tet_el++;
616 // oomph_info << "Tet" << std::endl;
617 break;
618
619 case 15:
620 // oomph_info << "Point" << std::endl;
621 break;
622
623 default:
624 std::string error_msg("Can't handle element type: ");
625 error_msg += oomph::StringConversion::to_string(el_type);
626 throw OomphLibError(
627 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
628 }
629
630 unsigned ntags;
631 mesh_file >> ntags;
632
633 // Gmesh produces at least two flags:
634
635 // Physical tag = what we use as boundary or region IDs.
636 int physical_tag = 0;
637 mesh_file >> physical_tag;
638
639 // Geometric tag; not something we use
640 int geom_tag = 0;
641 mesh_file >> geom_tag;
642
643 Vector<int> other_tags;
644 for (unsigned i = 2; i < ntags; i++)
645 {
646 int tag = 0;
647 mesh_file >> tag;
648 other_tags.push_back(tag);
649 }
650
651
652 // Now read the rest: node numbers
653 // https://stackoverflow.com/questions/16991002/getting-a-line-from-a-file-and-then-reading-word-by-word-c
655 std::getline(mesh_file, s);
656 std::istringstream iss(s);
657 std::string sub;
658 Vector<int> other_ints;
659 while (iss >> sub)
660 {
661 other_ints.push_back(atoi(sub.c_str()));
662 }
663 unsigned n_el_nod = other_ints.size();
664 for (unsigned j = 0; j < n_el_nod; j++)
665 {
666 unsigned node_number = unsigned(other_ints[j]);
667 element_node_index[el_number].push_back(node_number);
668
669 // If the element is a triangle, add node to boundary
670 // lookup scheme (if tag is not zero, i.e. hasn't been specified)
671 if (el_type == 2)
672 {
673 if (physical_tag != 0)
674 {
675 boundary_node[unsigned(physical_tag)].insert(node_number);
676 one_based_boundaries_of_node[node_number].insert(
677 unsigned(physical_tag));
678 if (unsigned(physical_tag) > highest_one_based_boundary_id)
679 {
680 highest_one_based_boundary_id = physical_tag;
681 }
682 }
683 }
684 }
685 // If it's a bulk element (tet) and the physical tag (region id)
686 // is not equal to 0 add it to region list
687 if (el_type == 4)
688 {
689 if (physical_tag != 0)
690 {
691 region_element[unsigned(physical_tag)].insert(n_tet_el);
692 }
693 }
694 }
695
696 mesh_file >> line;
697 if (line != "$EndElements")
698 {
699 std::string error_msg(
700 "Line has to contain the string \"$EndElements\"; ");
701 error_msg += " yours contains: " + line;
702 throw OomphLibError(
703 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
704 }
705
706 // Rewind
707 mesh_file.clear();
708 mesh_file.seekg(0);
709
710 mesh_file >> line;
711 if (line != "$MeshFormat")
712 {
713 std::string error_msg(
714 "First line has to contain the string \"$MeshFormat\"; ");
715 error_msg += " yours contains: " + line;
716 throw OomphLibError(
717 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
718 }
719
720 // Rest of line
721 mesh_file >> line;
722 mesh_file.close();
723
724
725 // Identify elements next to boundaries
726 //-------------------------------------
727
728 // Now loop over tet elements and check if three their nodes are
729 // on given boundary
730 std::map<unsigned, std::set<unsigned>> element_next_to_boundary;
731 for (std::map<unsigned, Vector<unsigned>>::iterator it =
732 element_node_index.begin();
733 it != element_node_index.end();
734 it++)
735 {
736 unsigned el_number = (*it).first;
737 unsigned nnod = ((*it).second).size();
738 if (nnod == 4)
739 {
740 std::map<unsigned, unsigned> boundary_node_count;
741 for (unsigned j = 0; j < nnod; j++)
742 {
743 unsigned node_number = ((*it).second)[j];
744 std::map<unsigned, std::set<unsigned>>::iterator itt =
745 one_based_boundaries_of_node.find(node_number);
746 if (itt != one_based_boundaries_of_node.end())
747 {
748 for (std::set<unsigned>::iterator ittt = (itt->second).begin();
749 ittt != (itt->second).end();
750 ittt++)
751 {
752 unsigned one_based_boundary_id = (*ittt);
753 boundary_node_count[one_based_boundary_id]++;
754 }
755 }
756 }
757 for (std::map<unsigned, unsigned>::iterator itt =
758 boundary_node_count.begin();
759 itt != boundary_node_count.end();
760 itt++)
761 {
762 if ((*itt).second == 3)
763 {
764 element_next_to_boundary[(*itt).first].insert(el_number);
765 }
766 }
767 }
768 }
769
770 this->set_nboundary(highest_one_based_boundary_id);
771
772 // Done reading/processing; now move across
773 // ----------------------------------------
774
775 // Translate enumeration
776 Vector<unsigned> oomph_lib_node_number(4);
777 oomph_lib_node_number[0] = 0;
778 oomph_lib_node_number[1] = 3;
779 oomph_lib_node_number[2] = 2;
780 oomph_lib_node_number[3] = 1;
781
782 // make space to accomodate all this in the actual mesh
783 Element_pt.resize(n_tet_el);
784 Node_pt.resize(nnod);
785
786 // Existing nodes
787 Vector<Node*> existing_node_pt(nnod, 0);
788
789 // Now process the simplex tets only (they have four nodes!)
790 unsigned e = 0;
791 for (std::map<unsigned, Vector<unsigned>>::iterator it =
792 element_node_index.begin();
793 it != element_node_index.end();
794 it++)
795 {
796 unsigned el_nnod = (*it).second.size();
797 if (el_nnod == 4)
798 {
799 // Here comes the new element
800 TElement<3, 2>* el_pt = new TElement<3, 2>;
801
802 // Store it
803 Element_pt[e] = el_pt;
804 e++;
805
806 // Make/get nodes
807 for (unsigned j = 0; j < el_nnod; j++)
808 {
809 // Get global, one-based node number
810 unsigned node_number = (*it).second[j];
811
812 // Does it already exist?
813 if (existing_node_pt[node_number - 1] != 0)
814 {
815 el_pt->node_pt(oomph_lib_node_number[j]) =
816 existing_node_pt[node_number - 1];
817 }
818 // Make new node
819 else
820 {
821 Node* nod_pt = 0;
822
823 // Is it on the boundary?
824 if (one_based_boundaries_of_node[node_number].size() == 0)
825 {
826 // Make normal node
827 nod_pt = el_pt->construct_node(oomph_lib_node_number[j]);
828 Node_pt[node_number - 1] = nod_pt;
829 existing_node_pt[node_number - 1] = nod_pt;
830 }
831 // Make boundary node
832 else
833 {
834 nod_pt =
835 el_pt->construct_boundary_node(oomph_lib_node_number[j]);
836 Node_pt[node_number - 1] = nod_pt;
837 existing_node_pt[node_number - 1] = nod_pt;
838
839 // Add to boundary lookup scheme
840 for (std::set<unsigned>::iterator it =
841 one_based_boundaries_of_node[node_number].begin();
842 it != one_based_boundaries_of_node[node_number].end();
843 it++)
844 {
845 add_boundary_node((*it) - 1, nod_pt);
846 }
847 }
848
849 // Assign coordinates of new node
850 Vector<double> x(node_coordinate[node_number]);
851 for (unsigned i = 0; i < 3; i++)
852 {
853 nod_pt->x(i) = x[i];
854 }
855 }
856 }
857 } // end for tet element
858
859 } // End of loop over all elements
860
861
862 // Setup region info. This is ugly because we're using
863 // a lookup scheme that was originally designed for tetgen...
864 unsigned n_region = region_element.size();
865 this->Region_element_pt.resize(n_region);
866 this->Region_attribute.resize(n_region);
867 unsigned region_count = 0;
868 for (std::map<unsigned, std::set<unsigned>>::iterator it =
869 region_element.begin();
870 it != region_element.end();
871 it++)
872 {
873 this->Region_element_pt[region_count].resize((*it).second.size());
874 unsigned one_based_region_id = (*it).first;
875 this->Region_attribute[region_count] = one_based_region_id - 1;
876 unsigned count = 0;
877 for (std::set<unsigned>::iterator itt = (*it).second.begin();
878 itt != (*it).second.end();
879 itt++)
880 {
881 unsigned one_based_el_number = (*itt);
882 this->Region_element_pt[region_count][count] =
883 dynamic_cast<FiniteElement*>(Element_pt[one_based_el_number - 1]);
884 count++;
885 }
886 region_count++;
887 }
888
889
890 // Keep alive for debugging
891 bool plot_for_debugging = false;
892 if (plot_for_debugging)
893 {
894 std::ofstream outfile;
895 std::string outfile_name;
896 std::string mesh_file_stem = "shite";
897
898 // Output element nodes
899 //----------------------
900 outfile_name = mesh_file_stem + "_element_nodes.dat";
901 outfile.open(outfile_name.c_str());
902 for (std::map<unsigned, Vector<unsigned>>::iterator it =
903 element_node_index.begin();
904 it != element_node_index.end();
905 it++)
906 {
907 std::set<unsigned> shite_nodes;
908 unsigned el_id = (*it).first;
909 std::stringstream tmp_out;
910 for (Vector<unsigned>::iterator itt = (*it).second.begin();
911 itt != (*it).second.end();
912 itt++)
913 {
914 unsigned node_number = (*itt);
915 shite_nodes.insert(node_number);
916 Vector<double> x(node_coordinate[node_number]);
917 unsigned n = x.size();
918 for (unsigned i = 0; i < n; i++)
919 {
920 tmp_out << x[i] << " ";
921 }
922 tmp_out << node_number << std::endl;
923 }
924 std::string prefix = ", N=4, E=1, F=FEPOINT, ET=TETRAHEDRON";
925 std::string postfix = "1 2 3 4";
926 if ((*it).second.size() == 3)
927 {
928 prefix = ", N=3, E=1, F=FEPOINT, ET=TRIANGLE";
929 postfix = "1 2 3";
930 }
931 outfile << "ZONE T=\"one-based element id=" << el_id << "\"" << prefix
932 << std::endl;
933 outfile << tmp_out.str();
934 outfile << postfix << std::endl;
935 }
936 outfile.close();
937
938 // Output boundary nodes
939 //----------------------
940 outfile_name = mesh_file_stem + "_boundary_nodes.dat";
941 outfile.open(outfile_name.c_str());
942 for (std::map<unsigned, std::set<unsigned>>::iterator it =
943 boundary_node.begin();
944 it != boundary_node.end();
945 it++)
946 {
947 unsigned one_based_boundary_id = (*it).first;
948 outfile << "ZONE T=\"one-based boundary id=" << one_based_boundary_id
949 << "\"" << std::endl;
950 for (std::set<unsigned>::iterator itt = (*it).second.begin();
951 itt != (*it).second.end();
952 itt++)
953 {
954 unsigned node_number = (*itt);
955 Vector<double> x(node_coordinate[node_number]);
956 unsigned n = x.size();
957 for (unsigned i = 0; i < n; i++)
958 {
959 outfile << x[i] << " ";
960 }
961 outfile << node_number << std::endl;
962 }
963 }
964 outfile.close();
965
966
967 // Output elements next to boundaries
968 //-----------------------------------
969 for (std::map<unsigned, std::set<unsigned>>::iterator it =
970 element_next_to_boundary.begin();
971 it != element_next_to_boundary.end();
972 it++)
973 {
974 unsigned one_based_boundary_id = (*it).first;
975 outfile_name =
976 mesh_file_stem + "_elements_next_to_boundary_" +
977 oomph::StringConversion::to_string(one_based_boundary_id) + ".dat";
978 outfile.open(outfile_name.c_str());
979 for (std::set<unsigned>::iterator itt = (*it).second.begin();
980 itt != (*it).second.end();
981 itt++)
982 {
983 outfile << "ZONE T=\"one-based boundary " << one_based_boundary_id
984 << "\", N=4, E=1, F=FEPOINT, ET=TETRAHEDRON" << std::endl;
985 unsigned el_number = (*itt);
986 unsigned nnod = element_node_index[el_number].size();
987 for (unsigned j = 0; j < nnod; j++)
988 {
989 unsigned node_number = element_node_index[el_number][j];
990 Vector<double> x(node_coordinate[node_number]);
991 unsigned n = x.size();
992 for (unsigned i = 0; i < n; i++)
993 {
994 outfile << x[i] << " ";
995 }
996 outfile << std::endl;
997 }
998 outfile << "1 2 3 4" << std::endl;
999 }
1000 outfile.close();
1001 }
1002
1003
1004 // Compute/report total volume for sanity check
1005 {
1006 double vol = 0.0;
1007 unsigned nel = this->nelement();
1008 for (unsigned e = 0; e < nel; e++)
1009 {
1010 vol += this->finite_element_pt(e)->size();
1011 }
1012 oomph_info << "Total volume of all elements in scaffold mesh: " << vol
1013 << std::endl;
1014 }
1015 }
1016 }
1017
1018
1019 /// Write geo file for gmsh
1020 void write_geo_file(const bool& use_mesh_grading_from_file)
1021 {
1022 // Transfer data from parameters
1023 TetMeshFacetedClosedSurface* outer_boundary_pt =
1025
1026 Vector<TetMeshFacetedSurface*> internal_surface_pt =
1028
1029 double element_volume = Gmsh_parameters_pt->element_volume();
1030
1031 std::string& target_size_file_name =
1033
1034 // Write gmsh geo file:
1035 std::string filename =
1037 std::ofstream geo_file;
1038 geo_file.open(filename.c_str());
1039
1040 geo_file << "// Uniform element size" << std::endl;
1041 geo_file << "//---------------------" << std::endl;
1042 geo_file << "lc=" << pow(element_volume, 1.0 / 3.0) << ";" << std::endl;
1043 geo_file << std::endl;
1044
1045 // Outer boundary
1046 //===============
1047
1048 // Create vertices
1049 geo_file << "// Outer box" << std::endl;
1050 geo_file << "//==========" << std::endl;
1051 geo_file << std::endl;
1052 geo_file << "// Vertices" << std::endl;
1053 geo_file << "//---------" << std::endl;
1054 std::map<TetMeshVertex*, unsigned> vertex_number;
1055 unsigned nv = outer_boundary_pt->nvertex();
1056 for (unsigned j = 0; j < nv; j++)
1057 {
1058 TetMeshVertex* vertex_pt = outer_boundary_pt->vertex_pt(j);
1059 vertex_number[vertex_pt] = j;
1060 geo_file << "Point(" << j + 1 << ")={";
1061 for (unsigned i = 0; i < 3; i++)
1062 {
1063 geo_file << vertex_pt->x(i) << ",";
1064 }
1065 geo_file << "lc};" << std::endl;
1066 }
1067
1068
1069 // Determine unique edges
1070 std::set<TetEdge> tet_edge_set;
1071 unsigned nfacet = outer_boundary_pt->nfacet();
1072 for (unsigned f = 0; f < nfacet; f++)
1073 {
1074 TetMeshFacet* facet_pt = outer_boundary_pt->facet_pt(f);
1075 unsigned nv = facet_pt->nvertex();
1076 for (unsigned j = 0; j < nv - 1; j++)
1077 {
1078 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1079 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1080 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1081 vertex_number[second_vertex_pt] + 1);
1082 tet_edge_set.insert(my_tet_edge);
1083 }
1084 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1085 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1086 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1087 vertex_number[second_vertex_pt] + 1);
1088 tet_edge_set.insert(my_tet_edge);
1089 }
1090
1091
1092 geo_file << std::endl;
1093 geo_file << "// Edge of outer box" << std::endl;
1094 geo_file << "//------------------" << std::endl;
1095 unsigned count = 0;
1096 std::map<TetEdge, unsigned> tet_edge;
1097 for (std::set<TetEdge>::iterator it = tet_edge_set.begin();
1098 it != tet_edge_set.end();
1099 it++)
1100 {
1101 tet_edge.insert(std::make_pair((*it), count));
1102 geo_file << "Line(" << count + 1 << ")={" << (*it).first_vertex_id()
1103 << "," << (*it).second_vertex_id() << "};" << std::endl;
1104
1105 count++;
1106 }
1107
1108
1109 geo_file << std::endl;
1110 geo_file << "// Faces of outer box" << std::endl;
1111 geo_file << "//-------------------" << std::endl;
1112 for (unsigned f = 0; f < nfacet; f++)
1113 {
1114 geo_file << "Line Loop(" << f + 1 << ")={";
1115
1116 TetMeshFacet* facet_pt = outer_boundary_pt->facet_pt(f);
1117 unsigned nv = facet_pt->nvertex();
1118 for (unsigned j = 0; j < nv - 1; j++)
1119 {
1120 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1121 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1122 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1123 vertex_number[second_vertex_pt] + 1);
1124
1125 std::map<TetEdge, unsigned>::iterator it = tet_edge.find(my_tet_edge);
1126 if (my_tet_edge.is_reversed())
1127 {
1128 geo_file << -int((it->second) + 1) << ",";
1129 }
1130 else
1131 {
1132 geo_file << ((it->second) + 1) << ",";
1133 }
1134 }
1135
1136 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1137 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1138 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1139 vertex_number[second_vertex_pt] + 1);
1140
1141 std::map<TetEdge, unsigned>::iterator it = tet_edge.find(my_tet_edge);
1142 if (my_tet_edge.is_reversed())
1143 {
1144 geo_file << -int((it->second) + 1) << "};" << std::endl;
1145 }
1146 else
1147 {
1148 geo_file << ((it->second) + 1) << "};" << std::endl;
1149 }
1150 geo_file << "Plane Surface(" << f + 1 << ")={" << f + 1 << "};"
1151 << std::endl;
1152 }
1153
1154
1155 geo_file << std::endl;
1156 geo_file << "// Define Plane Surfaces bounding the volume" << std::endl;
1157 geo_file << "//------------------------------------------" << std::endl;
1158 geo_file << "Surface Loop(1) = {";
1159 for (unsigned f = 0; f < nfacet - 1; f++)
1160 {
1161 geo_file << f + 1 << ",";
1162 }
1163 geo_file << nfacet << "};" << std::endl;
1164
1165
1166 geo_file << std::endl;
1167 geo_file << "// Define one-based boundary IDs" << std::endl;
1168 geo_file << "//------------------------------" << std::endl;
1169 for (unsigned f = 0; f < nfacet; f++)
1170 {
1171 unsigned one_based_boundary_id =
1172 outer_boundary_pt->one_based_facet_boundary_id(f);
1173 geo_file << "Physical Surface(" << one_based_boundary_id << ") = {"
1174 << f + 1 << "};" << std::endl;
1175 }
1176
1177
1178 // Offsets before we start adding the various internal volumes/surfaces
1179 Vector<unsigned> volume_id_to_be_subtracted_off;
1180 unsigned nvertex_offset = outer_boundary_pt->nvertex();
1181 unsigned nfacet_offset = outer_boundary_pt->nfacet();
1182 unsigned nvolume_offset = 1;
1183 unsigned nline_offset = tet_edge.size();
1184
1185
1186 // Storage for one-based ids of surfaces to be embedded in
1187 // main volume
1188 Vector<unsigned> surfaces_to_be_embedded_in_main_volume;
1189 std::map<unsigned, Vector<unsigned>>
1190 surfaces_to_be_embedded_in_specified_one_based_region;
1191
1192 // Now deal with internal faceted objects
1193 //=======================================
1194 unsigned n_internal = internal_surface_pt.size();
1195 for (unsigned i_internal = 0; i_internal < n_internal; i_internal++)
1196 {
1197 // Closed surface?
1198 TetMeshFacetedClosedSurface* closed_srf_pt =
1199 dynamic_cast<TetMeshFacetedClosedSurface*>(
1200 internal_surface_pt[i_internal]);
1201 bool inner_surface_is_closed = true;
1202 if (closed_srf_pt == 0)
1203 {
1204 inner_surface_is_closed = false;
1205 }
1206
1207 // What it says
1208 unsigned number_of_volumes_created_for_this_internal_object = 0;
1209
1210 // Create vertices
1211 geo_file << std::endl;
1212 geo_file << std::endl;
1213 geo_file << "// Inner faceted surface " << i_internal << std::endl;
1214 geo_file << "//==========================" << std::endl;
1215 geo_file << std::endl;
1216 geo_file << "// Vertices" << std::endl;
1217 geo_file << "//---------" << std::endl;
1218 std::map<TetMeshVertex*, unsigned> vertex_number;
1219 unsigned nv = internal_surface_pt[i_internal]->nvertex();
1220 for (unsigned j = 0; j < nv; j++)
1221 {
1222 TetMeshVertex* vertex_pt =
1223 internal_surface_pt[i_internal]->vertex_pt(j);
1224 vertex_number[vertex_pt] = nvertex_offset + j;
1225 geo_file << "Point(" << nvertex_offset + j + 1 << ")={";
1226 for (unsigned i = 0; i < 3; i++)
1227 {
1228 geo_file << vertex_pt->x(i) << ",";
1229 }
1230 geo_file << "lc};" << std::endl;
1231 }
1232
1233 // Determine unique edges
1234 std::set<TetEdge> tet_edge_set;
1235 unsigned nfacet = internal_surface_pt[i_internal]->nfacet();
1236 for (unsigned f = 0; f < nfacet; f++)
1237 {
1238 TetMeshFacet* facet_pt = internal_surface_pt[i_internal]->facet_pt(f);
1239 unsigned nv = facet_pt->nvertex();
1240 for (unsigned j = 0; j < nv - 1; j++)
1241 {
1242 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1243 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1244 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1245 vertex_number[second_vertex_pt] + 1);
1246 tet_edge_set.insert(my_tet_edge);
1247 }
1248 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1249 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1250 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1251 vertex_number[second_vertex_pt] + 1);
1252 tet_edge_set.insert(my_tet_edge);
1253 }
1254
1255 geo_file << std::endl;
1256 geo_file << "// Edge of inner faceted surface" << std::endl;
1257 geo_file << "//------------------------------" << std::endl;
1258 unsigned count = 0;
1259 std::map<TetEdge, unsigned> tet_edge;
1260 for (std::set<TetEdge>::iterator it = tet_edge_set.begin();
1261 it != tet_edge_set.end();
1262 it++)
1263 {
1264 tet_edge.insert(std::make_pair((*it), nline_offset + count));
1265 geo_file << "Line(" << nline_offset + count + 1 << ")={"
1266 << (*it).first_vertex_id() << "," << (*it).second_vertex_id()
1267 << "};" << std::endl;
1268
1269 count++;
1270 }
1271
1272
1273 geo_file << std::endl;
1274 geo_file << "// Faces of inner faceted surface" << std::endl;
1275 geo_file << "//-------------------------------" << std::endl;
1276 for (unsigned f = 0; f < nfacet; f++)
1277 {
1278 geo_file << "Line Loop(" << nfacet_offset + f + 1 << ")={";
1279
1280 TetMeshFacet* facet_pt = internal_surface_pt[i_internal]->facet_pt(f);
1281 unsigned nv = facet_pt->nvertex();
1282 for (unsigned j = 0; j < nv - 1; j++)
1283 {
1284 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(j);
1285 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(j + 1);
1286 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1287 vertex_number[second_vertex_pt] + 1);
1288
1289 std::map<TetEdge, unsigned>::iterator it =
1290 tet_edge.find(my_tet_edge);
1291 if (my_tet_edge.is_reversed())
1292 {
1293 geo_file << -int((it->second) + 1) << ",";
1294 }
1295 else
1296 {
1297 geo_file << ((it->second) + 1) << ",";
1298 }
1299 }
1300
1301 TetMeshVertex* first_vertex_pt = facet_pt->vertex_pt(nv - 1);
1302 TetMeshVertex* second_vertex_pt = facet_pt->vertex_pt(0);
1303 TetEdge my_tet_edge(vertex_number[first_vertex_pt] + 1,
1304 vertex_number[second_vertex_pt] + 1);
1305
1306 std::map<TetEdge, unsigned>::iterator it = tet_edge.find(my_tet_edge);
1307 if (my_tet_edge.is_reversed())
1308 {
1309 geo_file << -int((it->second) + 1) << "};" << std::endl;
1310 }
1311 else
1312 {
1313 geo_file << ((it->second) + 1) << "};" << std::endl;
1314 }
1315 geo_file << "Plane Surface(" << nfacet_offset + f + 1 << ")={"
1316 << nfacet_offset + f + 1 << "};" << std::endl;
1317
1318 // Keep track of plane surfaces that we've created so we
1319 // can embed them in volume for non-closed surfaces
1320 bool facet_is_embedded_in_a_volume =
1322 if (facet_is_embedded_in_a_volume)
1323 {
1324 unsigned one_based_region_id =
1326 if (one_based_region_id == 1)
1327 {
1328 surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset +
1329 f + 1);
1330 }
1331 else
1332 {
1333 surfaces_to_be_embedded_in_specified_one_based_region
1334 [one_based_region_id]
1335 .push_back(nfacet_offset + f + 1);
1336 }
1337 }
1338 else
1339 {
1340 // By default put all surfaces into main volume
1341 if (!inner_surface_is_closed)
1342 {
1343 surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset +
1344 f + 1);
1345 }
1346 }
1347 }
1348
1349 // Store all region IDs defined by bounding facets
1350 std::set<unsigned> all_regions_id;
1351
1352 // region_bounding_facet[r][i] returns the facet id (in gmsh
1353 // counting) of i-th facet bounding region r
1354 std::map<unsigned, Vector<unsigned>> region_bounding_facet;
1355
1356 // outer_bounding_facet[i] returns the facet id (in gmsh
1357 // counting) that encloses the actual regions and acts as a
1358 // hole in the main volume (volume 1)
1359 Vector<unsigned> outer_bounding_facet;
1360
1361 // Loop pver all facets to figure out which regions are bounded
1362 // by them
1363 for (unsigned f = 0; f < nfacet; f++)
1364 {
1365 TetMeshFacet* facet_pt = internal_surface_pt[i_internal]->facet_pt(f);
1366 std::set<unsigned> region_id(
1367 facet_pt->one_based_adjacent_region_id());
1368 unsigned nr = region_id.size();
1369 if (nr == 1) outer_bounding_facet.push_back(nfacet_offset + f + 1);
1370 if ((nr == 0) && inner_surface_is_closed)
1371 {
1372 // Add to list of plane surfaces that don't bound regions so we
1373 // can embed them in volume for non-closed surfaces
1374 surfaces_to_be_embedded_in_main_volume.push_back(nfacet_offset + f +
1375 1);
1376 }
1377 for (std::set<unsigned>::iterator it = region_id.begin();
1378 it != region_id.end();
1379 it++)
1380 {
1381 all_regions_id.insert((*it));
1382 region_bounding_facet[(*it)].push_back(nfacet_offset + f + 1);
1383 }
1384 }
1385
1386 // Number of regions bounded by facets
1387 unsigned n_regions_bounded_by_facets = all_regions_id.size();
1388
1389 // No bounded regions
1390 if (n_regions_bounded_by_facets == 0)
1391 {
1392 if (inner_surface_is_closed)
1393 {
1394 std::ostringstream error_message;
1395 error_message
1396 << "Something fishy going on! "
1397 << "Internal faceted surface " << i_internal
1398 << " is closed but does not bound any regions!\n"
1399 << "Specify one-based region ID for all facets using\n\n"
1400 << " TetMeshFacet::set_one_based_adjacent_region_id(...)\n\n";
1401 throw OomphLibError(error_message.str(),
1402 OOMPH_CURRENT_FUNCTION,
1403 OOMPH_EXCEPTION_LOCATION);
1404 }
1405 }
1406 else
1407 {
1408 // Do we need to create a hole in the main volume that contains
1409 // the multiple sub-volumes/regions?
1410 unsigned offset_for_extra_hole = 0;
1411 if (n_regions_bounded_by_facets > 1)
1412 {
1413 geo_file << std::endl;
1414 geo_file << "// Define Plane Surfaces bounding compound volume"
1415 << std::endl
1416 << "//-----------------------------------------------"
1417 << std::endl;
1418 geo_file << "// that will have to be treated as hole in main volume"
1419 << std::endl
1420 << "//-----------------------------------------------"
1421 << std::endl;
1422 offset_for_extra_hole = 1;
1423 geo_file << "Surface Loop(" << nvolume_offset + 1 << ") = {";
1424 unsigned n = outer_bounding_facet.size();
1425 for (unsigned f = 0; f < n - 1; f++)
1426 {
1427 geo_file << outer_bounding_facet[f] << ",";
1428 }
1429 geo_file << outer_bounding_facet[n - 1] << "};" << std::endl;
1430
1431 // Bump
1432 number_of_volumes_created_for_this_internal_object++;
1433 }
1434
1435 // Need to remove the internal volume from the outer volume
1436 volume_id_to_be_subtracted_off.push_back(nvolume_offset + 1);
1437
1438 // Deal with actual regions that fill the hole
1439 unsigned count = 0;
1440 for (std::map<unsigned, Vector<unsigned>>::iterator it =
1441 region_bounding_facet.begin();
1442 it != region_bounding_facet.end();
1443 it++)
1444 {
1445 geo_file << std::endl;
1446 geo_file << "// Define Plane Surfaces bounding the region volume "
1447 << (*it).first << std::endl;
1448 geo_file << "//----------------------------------------------------"
1449 << std::endl;
1450 geo_file << "Surface Loop("
1451 << nvolume_offset + 1 + offset_for_extra_hole + count
1452 << ") = {";
1453 unsigned n = (*it).second.size();
1454 for (unsigned f = 0; f < n - 1; f++)
1455 {
1456 geo_file << ((*it).second)[f] << ",";
1457 }
1458 geo_file << ((*it).second)[n - 1] << "};" << std::endl;
1459
1460 geo_file << std::endl;
1461 geo_file << "// Define volume "
1462 << nvolume_offset + 1 + offset_for_extra_hole + count
1463 << " as the volume bounded by Surface Loop "
1464 << nvolume_offset + 1 + offset_for_extra_hole + count
1465 << std::endl;
1466 geo_file
1467 << "//--------------------------------------------------------"
1468 << std::endl;
1469 geo_file << "Volume("
1470 << nvolume_offset + 1 + offset_for_extra_hole + count
1471 << ")={"
1472 << nvolume_offset + 1 + offset_for_extra_hole + count
1473 << "};" << std::endl;
1474 geo_file << std::endl;
1475
1476 // Bump
1477 number_of_volumes_created_for_this_internal_object++;
1478
1479
1480 // Add as physical (to be meshed) volume if it's not a volume
1481 bool mesh_the_volume = true;
1482 if (closed_srf_pt != 0)
1483 {
1484 if (closed_srf_pt->faceted_volume_represents_hole_for_gmsh())
1485 {
1486 mesh_the_volume = false;
1487 }
1488 }
1489 if (mesh_the_volume)
1490 {
1491 geo_file << "// Define one-based region IDs" << std::endl;
1492 geo_file << "//----------------------------" << std::endl;
1493 geo_file << "Physical Volume(" << (*it).first << ")={"
1494 << nvolume_offset + 1 + offset_for_extra_hole + count
1495 << "};" << std::endl;
1496
1497 unsigned ns_embedded =
1498 surfaces_to_be_embedded_in_specified_one_based_region[(*it)
1499 .first]
1500 .size();
1501 if (ns_embedded > 0)
1502 {
1503 geo_file << "// This region has " << ns_embedded
1504 << " embedded surfaces\n";
1505 geo_file << "Surface{";
1506 for (unsigned i = 0; i < ns_embedded - 1; i++)
1507 {
1508 geo_file
1509 << surfaces_to_be_embedded_in_specified_one_based_region
1510 [(*it).first][i]
1511 << ",";
1512 }
1513 geo_file
1514 << surfaces_to_be_embedded_in_specified_one_based_region
1515 [(*it).first][ns_embedded - 1]
1516 << "}In Volume {"
1517 << nvolume_offset + 1 + offset_for_extra_hole + count << "};"
1518 << std::endl;
1519 }
1520 }
1521
1522 count++;
1523 }
1524 }
1525
1526 geo_file << std::endl;
1527 geo_file << "// Define one-based boundary IDs" << std::endl;
1528 geo_file << "//------------------------------" << std::endl;
1529 for (unsigned f = 0; f < nfacet; f++)
1530 {
1531 unsigned one_based_boundary_id =
1532 internal_surface_pt[i_internal]->one_based_facet_boundary_id(f);
1533 geo_file << "Physical Surface(" << one_based_boundary_id << ") = {"
1534 << nfacet_offset + f + 1 << "};" << std::endl;
1535 }
1536
1537 // Bump
1538 nvertex_offset += internal_surface_pt[i_internal]->nvertex();
1539 nfacet_offset += internal_surface_pt[i_internal]->nfacet();
1540 nvolume_offset += number_of_volumes_created_for_this_internal_object;
1541 nline_offset += tet_edge.size();
1542 }
1543
1544
1545 // Done with the internal regions; write the actual volume
1546 // and specify embedded surfaces
1547 geo_file << std::endl;
1548 geo_file << "// Define volume 1 as the volume bounded by Surface Loop 1"
1549 << std::endl;
1550 geo_file << "//--------------------------------------------------------"
1551 << std::endl;
1552 unsigned n = volume_id_to_be_subtracted_off.size();
1553 if (n > 0)
1554 {
1555 geo_file << "// with volume[s]: ";
1556 for (unsigned i = 0; i < n; i++)
1557 {
1558 geo_file << volume_id_to_be_subtracted_off[i] << " ";
1559 }
1560 geo_file << "removed." << std::endl;
1561 geo_file << "//--------------------------------------------------------"
1562 << std::endl;
1563 }
1564
1565 // Add initial volume
1566 geo_file << "Volume(1)={1";
1567
1568 // Remove volumes that has separate IDs
1569 for (unsigned i = 0; i < n; i++)
1570 {
1571 geo_file << "," << volume_id_to_be_subtracted_off[i];
1572 }
1573 geo_file << "};" << std::endl;
1574 geo_file << std::endl;
1575
1576
1577 geo_file << "// Define one-based region IDs" << std::endl;
1578 geo_file << "//----------------------------" << std::endl;
1579 geo_file << "Physical Volume(1"
1580 << ")={1};" << std::endl;
1581
1582
1583 // Now embed any surfaces that don't bound volumes
1584 unsigned ns = surfaces_to_be_embedded_in_main_volume.size();
1585 if (ns > 0)
1586 {
1587 geo_file << std::endl;
1588 geo_file << "// Embed Plane Surfaces in main volume (volume 1)"
1589 << std::endl;
1590 geo_file << "//-----------------------------------------------"
1591 << std::endl;
1592 geo_file << "Surface{";
1593 for (unsigned s = 0; s < ns - 1; s++)
1594 {
1595 geo_file << surfaces_to_be_embedded_in_main_volume[s] << ",";
1596 }
1597 geo_file << surfaces_to_be_embedded_in_main_volume[ns - 1]
1598 << "} In Volume{1};" << std::endl;
1599 geo_file << std::endl;
1600 }
1601
1602
1603 // Mesh grading
1604 {
1605 if (use_mesh_grading_from_file)
1606 {
1607#ifdef PARANOID
1608
1609 // Open wide...
1610 std::ifstream file(target_size_file_name.c_str(), std::ios_base::in);
1611
1612 // Check that the file actually opened correctly
1613 if (!file.is_open())
1614 {
1615 std::string error_msg("Failed to open target volume file: ");
1616 error_msg += "\"" + target_size_file_name;
1617 throw OomphLibError(
1618 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1619 }
1620#endif
1621
1622 geo_file << "Field[1]=Structured;" << std::endl;
1623 geo_file << "Field[1].FileName=\"" << target_size_file_name << "\";"
1624 << std::endl;
1625 geo_file << "Field[1].TextFormat=1;" << std::endl;
1626 geo_file << "Background Field = 1;" << std::endl;
1627 }
1628 }
1629
1630 // Mesh the bloody thing
1631 geo_file << std::endl;
1632 geo_file << "Mesh 3;" << std::endl;
1633
1634
1635 /* // hierher Christmas attempt at optimisation */
1636 /* geo_file << std::endl; */
1637 /* geo_file << "// Attempt at optimisation:
1638 * http://onelab.info/pipermail/gmsh/2015/010126.html" << std::endl; */
1639 /* geo_file << "//----------------------------" << std::endl; */
1640 /* geo_file << "Mesh.Optimize=1;" << std::endl; */
1641 /* geo_file << "Mesh.OptimizeNetgen=1;" << std::endl; */
1642
1643
1644 geo_file.close();
1645 }
1646
1647 /// Parameters
1649 };
1650
1651
1652 /// ///////////////////////////////////////////////////////////////////////
1653 /// ///////////////////////////////////////////////////////////////////////
1654 /// ///////////////////////////////////////////////////////////////////////
1655
1656
1657 //=========================================================================
1658 // Gmsh-based Tet Mesh
1659 //=========================================================================
1660 template<class ELEMENT>
1661 class GmshTetMesh : public virtual TetMeshBase, public virtual Mesh
1662 {
1663 public:
1664 /// Constructor
1665 GmshTetMesh(GmshParameters* gmsh_parameters_pt,
1666 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1667 : Gmsh_parameters_pt(gmsh_parameters_pt)
1668 {
1669 bool use_mesh_grading_from_file = false;
1670 build_it(time_stepper_pt, use_mesh_grading_from_file);
1671 }
1672
1673 /// Constructor. If boolean is set
1674 /// to true, the target element sizes are read from file (used during
1675 /// adaptation; otherwise uniform target size is used).
1676 GmshTetMesh(GmshParameters* gmsh_parameters_pt,
1677 const bool& use_mesh_grading_from_file,
1678 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1679 : Gmsh_parameters_pt(gmsh_parameters_pt)
1680 {
1681 build_it(time_stepper_pt, use_mesh_grading_from_file);
1682 }
1683
1684
1685 private:
1686 // Build function
1687 void build_it(TimeStepper* time_stepper_pt,
1688 const bool& use_mesh_grading_from_file)
1689 {
1690 // Mesh can only be built with 3D Telements.
1691 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
1692
1693 // Transfer data from parameters
1694 TetMeshFacetedClosedSurface* outer_boundary_pt =
1696
1697 Vector<TetMeshFacetedSurface*> internal_surface_pt =
1699
1700 // Remember timestepper
1701 Time_stepper_pt = time_stepper_pt;
1702
1703 // Store the boundary
1704 Outer_boundary_pt = outer_boundary_pt;
1705
1706 // Setup reverse lookup scheme
1707 {
1708 unsigned n_facet = Outer_boundary_pt->nfacet();
1709 for (unsigned f = 0; f < n_facet; f++)
1710 {
1712 if (b != 0)
1713 {
1716 }
1717 else
1718 {
1719 std::ostringstream error_message;
1720 error_message << "Boundary IDs have to be one-based. Yours is " << b
1721 << "\n";
1722 throw OomphLibError(error_message.str(),
1723 OOMPH_CURRENT_FUNCTION,
1724 OOMPH_EXCEPTION_LOCATION);
1725 }
1726 }
1727 }
1728
1729
1730 // Store the internal boundary
1731 Internal_surface_pt = internal_surface_pt;
1732
1733 // Setup reverse lookup scheme
1734 {
1735 unsigned n = Internal_surface_pt.size();
1736 for (unsigned i = 0; i < n; i++)
1737 {
1738 unsigned n_facet = Internal_surface_pt[i]->nfacet();
1739 for (unsigned f = 0; f < n_facet; f++)
1740 {
1741 unsigned b = Internal_surface_pt[i]->one_based_facet_boundary_id(f);
1742 if (b != 0)
1743 {
1745 Tet_mesh_facet_pt[b - 1] = Internal_surface_pt[i]->facet_pt(f);
1746 }
1747 else
1748 {
1749 std::ostringstream error_message;
1750 error_message << "Boundary IDs have to be one-based. Yours is "
1751 << b << "\n";
1752 throw OomphLibError(error_message.str(),
1753 OOMPH_CURRENT_FUNCTION,
1754 OOMPH_EXCEPTION_LOCATION);
1755 }
1756 }
1757 }
1758 }
1759
1760 // Build scaffold
1761 GmshTetScaffoldMesh* tmp_scaffold_mesh_pt =
1762 new GmshTetScaffoldMesh(Gmsh_parameters_pt, use_mesh_grading_from_file);
1763
1764 // Convert mesh from scaffold to actual mesh
1765 build_from_scaffold(tmp_scaffold_mesh_pt, time_stepper_pt);
1766
1767 // Kill the scaffold
1768 delete tmp_scaffold_mesh_pt;
1769 tmp_scaffold_mesh_pt = 0;
1770
1771 // Setup boundary coordinates
1772 unsigned nb = nboundary();
1773 for (unsigned b = 0; b < nb; b++)
1774 {
1775 bool switch_normal = false;
1776 setup_boundary_coordinates<ELEMENT>(b, switch_normal);
1777 }
1778
1779 // Now snap onto geometric objects associated with triangular facets
1780 // (if any!)
1782 }
1783
1784
1785 protected:
1786 /// Parameters
1788
1789
1790 private:
1791 /// Build unstructured tet gmesh mesh based on output from scaffold
1792 void build_from_scaffold(GmshTetScaffoldMesh* tmp_scaffold_mesh_pt,
1793 TimeStepper* time_stepper_pt)
1794 {
1795 // Mesh can only be built with 3D Telements.
1796 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
1797
1798 // Create space for elements
1799 unsigned nelem = tmp_scaffold_mesh_pt->nelement();
1800 Element_pt.resize(nelem);
1801
1802 // Relation between elements pointers and numbers in old mesh
1803 std::map<FiniteElement*, unsigned> scaffold_mesh_element_number;
1804
1805 // Create space for nodes
1806 unsigned nnode_scaffold = tmp_scaffold_mesh_pt->nnode();
1807 Node_pt.resize(nnode_scaffold);
1808
1809 // Set number of boundaries
1810 unsigned nbound = tmp_scaffold_mesh_pt->nboundary();
1811 set_nboundary(nbound);
1812
1813 // Resize boundary info stuff
1814 Boundary_element_pt.resize(nbound); // hierher shouldn't this be a map?
1815 Face_index_at_boundary.resize(nbound); // hierher shouldn't this be a map?
1816 Boundary_region_element_pt.resize(nbound);
1817 Face_index_region_at_boundary.resize(nbound);
1818
1819 // Build elements
1820 for (unsigned e = 0; e < nelem; e++)
1821 {
1822 Element_pt[e] = new ELEMENT;
1823 }
1824
1825 // Number of nodes per element
1826 unsigned nnod_el = tmp_scaffold_mesh_pt->finite_element_pt(0)->nnode();
1827
1828 // Setup map to check the (pseudo-)global node number
1829 // Nodes whose number is zero haven't been copied across
1830 // into the mesh yet.
1831 std::map<Node*, unsigned> global_number;
1832 unsigned global_count = 0;
1833
1834 // Loop over elements in scaffold mesh, visit their nodes
1835 for (unsigned e = 0; e < nelem; e++)
1836 {
1837 // Setup reverse lookup scheme to decipher lookup schemes from
1838 // scaffold mesh
1839 scaffold_mesh_element_number[tmp_scaffold_mesh_pt->finite_element_pt(
1840 e)] = e;
1841
1842 // Loop over all nodes in element
1843 for (unsigned j = 0; j < nnod_el; j++)
1844 {
1845 // Pointer to node in the scaffold mesh
1846 Node* scaffold_node_pt =
1847 tmp_scaffold_mesh_pt->finite_element_pt(e)->node_pt(j);
1848
1849 // Get the (pseudo-)global node number in scaffold mesh
1850 // (It's zero [=default] if not visited this one yet)
1851 unsigned j_global = global_number[scaffold_node_pt];
1852
1853 // Haven't done this one yet
1854 if (j_global == 0)
1855 {
1856 // Get pointer to set of mesh boundaries that this
1857 // scaffold node occupies; NULL if the node is not on any boundary
1858 std::set<unsigned>* boundaries_pt;
1859 scaffold_node_pt->get_boundaries_pt(boundaries_pt);
1860
1861 // Is it on boundaries?
1862 if (boundaries_pt != 0)
1863 {
1864 // Create new boundary node
1866 j, time_stepper_pt);
1867
1868 // Give it a number (not necessarily the global node
1869 // number in the scaffold mesh -- we just need something
1870 // to keep track...)
1871 global_count++;
1872 global_number[scaffold_node_pt] = global_count;
1873
1874 // Add to boundaries
1875 for (std::set<unsigned>::iterator it = boundaries_pt->begin();
1876 it != boundaries_pt->end();
1877 ++it)
1878 {
1879 add_boundary_node(*it, new_node_pt);
1880 }
1881 }
1882 // Build normal node
1883 else
1884 {
1885 // Create new normal node
1886 finite_element_pt(e)->construct_node(j, time_stepper_pt);
1887
1888 // Give it a number (not necessarily the global node
1889 // number in the scaffold mesh -- we just need something
1890 // to keep track...)
1891 global_count++;
1892 global_number[scaffold_node_pt] = global_count;
1893 }
1894
1895 // Copy new node, created using the NEW element's construct_node
1896 // function into global storage, using the same global
1897 // node number that we've just associated with the
1898 // corresponding node in the scaffold mesh
1899 Node_pt[global_count - 1] = finite_element_pt(e)->node_pt(j);
1900
1901 // Assign coordinates
1902 Node_pt[global_count - 1]->x(0) = scaffold_node_pt->x(0);
1903 Node_pt[global_count - 1]->x(1) = scaffold_node_pt->x(1);
1904 Node_pt[global_count - 1]->x(2) = scaffold_node_pt->x(2);
1905 }
1906 // This one has already been done: Copy across
1907 else
1908 {
1909 finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
1910 }
1911 }
1912
1913
1914 // Now figure out which boundaries the faces are on
1916 for (unsigned f = 0; f < 4; f++)
1917 {
1918 Node* face_node0_pt = 0;
1919 Node* face_node1_pt = 0;
1920 Node* face_node2_pt = 0;
1921
1922 switch (f)
1923 {
1924 case 0:
1925 face_node0_pt = fe_pt->node_pt(1);
1926 face_node1_pt = fe_pt->node_pt(2);
1927 face_node2_pt = fe_pt->node_pt(3);
1928 break;
1929
1930 case 1:
1931 face_node0_pt = fe_pt->node_pt(0);
1932 face_node1_pt = fe_pt->node_pt(2);
1933 face_node2_pt = fe_pt->node_pt(3);
1934 break;
1935
1936 case 2:
1937 face_node0_pt = fe_pt->node_pt(0);
1938 face_node1_pt = fe_pt->node_pt(1);
1939 face_node2_pt = fe_pt->node_pt(3);
1940 break;
1941
1942 case 3:
1943 face_node0_pt = fe_pt->node_pt(0);
1944 face_node1_pt = fe_pt->node_pt(1);
1945 face_node2_pt = fe_pt->node_pt(2);
1946 break;
1947
1948 default:
1949
1950 std::ostringstream error_message;
1951 error_message << "Wrong face number " << f << std::endl;
1952 throw OomphLibError(error_message.str(),
1953 OOMPH_CURRENT_FUNCTION,
1954 OOMPH_EXCEPTION_LOCATION);
1955 }
1956
1957 // If any of the boundary sets are empty we don't have
1958 // three nodes on the boundary...
1959 std::set<unsigned>* bc0_pt;
1960 face_node0_pt->get_boundaries_pt(bc0_pt);
1961 if (bc0_pt != 0)
1962 {
1963 std::set<unsigned>* bc1_pt;
1964 face_node1_pt->get_boundaries_pt(bc1_pt);
1965 if (bc1_pt != 0)
1966 {
1967 std::set<unsigned>* bc2_pt;
1968 face_node2_pt->get_boundaries_pt(bc2_pt);
1969 if (bc2_pt != 0)
1970 {
1971 std::set<unsigned> common_bound_0_and_1;
1972 std::set_intersection(
1973 bc0_pt->begin(),
1974 bc0_pt->end(),
1975 bc1_pt->begin(),
1976 bc1_pt->end(),
1977 std::inserter(common_bound_0_and_1,
1978 common_bound_0_and_1.begin()));
1979 std::set<unsigned> common_bound_0_and_1_and_2;
1980 std::set_intersection(
1981 common_bound_0_and_1.begin(),
1982 common_bound_0_and_1.end(),
1983 bc2_pt->begin(),
1984 bc2_pt->end(),
1985 std::inserter(common_bound_0_and_1_and_2,
1986 common_bound_0_and_1_and_2.begin()));
1987 for (std::set<unsigned>::iterator it =
1988 common_bound_0_and_1_and_2.begin();
1989 it != common_bound_0_and_1_and_2.end();
1990 it++)
1991 {
1992 Boundary_element_pt[(*it)].push_back(fe_pt);
1993 Face_index_at_boundary[(*it)].push_back(f);
1994 }
1995 }
1996 }
1997 }
1998 }
1999 }
2000
2001 // Copy across region information (scaffold mesh is a friend)
2002 unsigned nr = tmp_scaffold_mesh_pt->Region_element_pt.size();
2003 Region_attribute.resize(nr);
2004 Region_element_pt.resize(nr);
2005 for (unsigned i = 0; i < nr; i++)
2006 { //--
2007 Region_attribute[i] = tmp_scaffold_mesh_pt->Region_attribute[i];
2008 unsigned nel = tmp_scaffold_mesh_pt->Region_element_pt[i].size();
2009 Region_element_pt[i].resize(nel);
2010 for (unsigned e = 0; e < nel; e++)
2011 {
2012 FiniteElement* scaff_el_pt =
2013 tmp_scaffold_mesh_pt->Region_element_pt[i][e];
2014 unsigned scaff_el_number = scaffold_mesh_element_number[scaff_el_pt];
2016 dynamic_cast<FiniteElement*>(Element_pt[scaff_el_number]);
2017
2018
2019 // Now figure out which boundaries the faces are on (again!)
2021 for (unsigned f = 0; f < 4; f++)
2022 {
2023 Node* face_node0_pt = 0;
2024 Node* face_node1_pt = 0;
2025 Node* face_node2_pt = 0;
2026
2027 switch (f)
2028 {
2029 case 0:
2030 face_node0_pt = fe_pt->node_pt(1);
2031 face_node1_pt = fe_pt->node_pt(2);
2032 face_node2_pt = fe_pt->node_pt(3);
2033 break;
2034
2035 case 1:
2036 face_node0_pt = fe_pt->node_pt(0);
2037 face_node1_pt = fe_pt->node_pt(2);
2038 face_node2_pt = fe_pt->node_pt(3);
2039 break;
2040
2041 case 2:
2042 face_node0_pt = fe_pt->node_pt(0);
2043 face_node1_pt = fe_pt->node_pt(1);
2044 face_node2_pt = fe_pt->node_pt(3);
2045 break;
2046
2047 case 3:
2048 face_node0_pt = fe_pt->node_pt(0);
2049 face_node1_pt = fe_pt->node_pt(1);
2050 face_node2_pt = fe_pt->node_pt(2);
2051 break;
2052
2053 default:
2054 std::ostringstream error_message;
2055 error_message << "Wrong face number " << f << std::endl;
2056 throw OomphLibError(error_message.str(),
2057 OOMPH_CURRENT_FUNCTION,
2058 OOMPH_EXCEPTION_LOCATION);
2059 }
2060
2061 // If any of the boundary sets are empty we don't have
2062 // three nodes on the boundary...
2063 std::set<unsigned>* bc0_pt;
2064 face_node0_pt->get_boundaries_pt(bc0_pt);
2065 if (bc0_pt != 0)
2066 {
2067 std::set<unsigned>* bc1_pt;
2068 face_node1_pt->get_boundaries_pt(bc1_pt);
2069 if (bc1_pt != 0)
2070 {
2071 std::set<unsigned>* bc2_pt;
2072 face_node2_pt->get_boundaries_pt(bc2_pt);
2073 if (bc2_pt != 0)
2074 {
2075 std::set<unsigned> common_bound_0_and_1;
2076 std::set_intersection(
2077 bc0_pt->begin(),
2078 bc0_pt->end(),
2079 bc1_pt->begin(),
2080 bc1_pt->end(),
2081 std::inserter(common_bound_0_and_1,
2082 common_bound_0_and_1.begin()));
2083 std::set<unsigned> common_bound_0_and_1_and_2;
2084 std::set_intersection(
2085 common_bound_0_and_1.begin(),
2086 common_bound_0_and_1.end(),
2087 bc2_pt->begin(),
2088 bc2_pt->end(),
2089 std::inserter(common_bound_0_and_1_and_2,
2090 common_bound_0_and_1_and_2.begin()));
2091 for (std::set<unsigned>::iterator it =
2092 common_bound_0_and_1_and_2.begin();
2093 it != common_bound_0_and_1_and_2.end();
2094 it++)
2095 {
2097 .push_back(fe_pt);
2099 .push_back(f);
2100 }
2101 }
2102 }
2103 }
2104 }
2105 }
2106 }
2107
2108 // Lookup scheme has now been setup
2110
2111
2112 // At this point we've created all the elements and
2113 // created their vertex nodes. Now we need to create
2114 // the additional (midside and internal) nodes!
2115
2116 // Get number of nodes along element edge
2117 unsigned n_node_1d = finite_element_pt(0)->nnode_1d();
2118
2119 // At the moment we're only able to deal with nnode_1d=2 or 3.
2120 if ((n_node_1d != 2) && (n_node_1d != 3))
2121 {
2122 std::ostringstream error_message;
2123 error_message << "Mesh generation from gmsh currently only works\n";
2124 error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
2125 error_message << "for nnode_1d=" << n_node_1d << std::endl;
2126 throw OomphLibError(error_message.str(),
2127 OOMPH_CURRENT_FUNCTION,
2128 OOMPH_EXCEPTION_LOCATION);
2129 }
2130
2131 // Storage for the local coordinate of the new node
2132 Vector<double> s(3);
2133
2134 // Map associating edge with midside node
2135 std::map<std::pair<Node*, Node*>, Node*> midside_node_pt;
2136
2137 // Loop over all elements
2138 for (unsigned e = 0; e < nelem; e++)
2139 {
2140 // Cache pointers to the elements
2141 FiniteElement* const el_pt = this->finite_element_pt(e);
2142 FiniteElement* const simplex_el_pt =
2143 tmp_scaffold_mesh_pt->finite_element_pt(e);
2144
2145 // Loop over the edges
2146 for (unsigned j = 0; j < 6; ++j)
2147 {
2148 Node* nod0_pt = 0;
2149 Node* nod1_pt = 0;
2150 unsigned new_node_number = 0;
2151 std::pair<Node*, Node*> edge;
2152 switch (j)
2153 {
2154 case 0:
2155 nod0_pt = el_pt->node_pt(0);
2156 nod1_pt = el_pt->node_pt(1);
2157 new_node_number = 4;
2158 break;
2159
2160 case 1:
2161 nod0_pt = el_pt->node_pt(0);
2162 nod1_pt = el_pt->node_pt(2);
2163 new_node_number = 5;
2164 break;
2165
2166 case 2:
2167 nod0_pt = el_pt->node_pt(0);
2168 nod1_pt = el_pt->node_pt(3);
2169 new_node_number = 6;
2170 break;
2171
2172 case 3:
2173 nod0_pt = el_pt->node_pt(1);
2174 nod1_pt = el_pt->node_pt(2);
2175 new_node_number = 7;
2176 break;
2177
2178 case 4:
2179 nod0_pt = el_pt->node_pt(2);
2180 nod1_pt = el_pt->node_pt(3);
2181 new_node_number = 8;
2182 break;
2183
2184 case 5:
2185 nod0_pt = el_pt->node_pt(1);
2186 nod1_pt = el_pt->node_pt(3);
2187 new_node_number = 9;
2188 break;
2189
2190 default:
2191 std::ostringstream error_message;
2192 error_message << "Wrong edge number " << j << std::endl;
2193 throw OomphLibError(error_message.str(),
2194 OOMPH_CURRENT_FUNCTION,
2195 OOMPH_EXCEPTION_LOCATION);
2196 }
2197
2198 // Identify existence of node via edge
2199 edge = std::make_pair(std::min(nod0_pt, nod1_pt),
2200 std::max(nod0_pt, nod1_pt));
2201 Node* existing_node_pt = midside_node_pt[edge];
2202 if (existing_node_pt == 0)
2203 {
2204 // If any of the boundary sets are empty we don't have
2205 // two edge nodes on the boundary...
2206 std::set<unsigned> common_bound_0_and_1;
2207 std::set<unsigned>* bc0_pt;
2208 nod0_pt->get_boundaries_pt(bc0_pt);
2209 if (bc0_pt != 0)
2210 {
2211 std::set<unsigned>* bc1_pt;
2212 nod1_pt->get_boundaries_pt(bc1_pt);
2213 if (bc1_pt != 0)
2214 {
2215 std::set_intersection(
2216 bc0_pt->begin(),
2217 bc0_pt->end(),
2218 bc1_pt->begin(),
2219 bc1_pt->end(),
2220 std::inserter(common_bound_0_and_1,
2221 common_bound_0_and_1.begin()));
2222 }
2223 }
2224 // Storage for the new node
2225 Node* new_node_pt = 0;
2226
2227 // New non-boundary node:
2228 if (common_bound_0_and_1.size() == 0)
2229 {
2230 new_node_pt =
2231 el_pt->construct_node(new_node_number, time_stepper_pt);
2232 }
2233 // New boundary node
2234 else
2235 {
2236 new_node_pt = el_pt->construct_boundary_node(new_node_number,
2237 time_stepper_pt);
2238 for (std::set<unsigned>::iterator it =
2239 common_bound_0_and_1.begin();
2240 it != common_bound_0_and_1.end();
2241 it++)
2242 {
2243 this->add_boundary_node((*it), new_node_pt);
2244 }
2245 }
2246
2247 // Find the local coordinates of the node
2248 el_pt->local_coordinate_of_node(new_node_number, s);
2249
2250 // Find the coordinates of the new node from the existing
2251 // and fully-functional element in the scaffold mesh
2252 for (unsigned i = 0; i < 3; i++)
2253 {
2254 new_node_pt->x(i) = simplex_el_pt->interpolated_x(s, i);
2255 }
2256
2257 // Associate node with edge
2258 midside_node_pt[edge] = new_node_pt;
2259
2260 // Add node to mesh
2261 Node_pt.push_back(new_node_pt);
2262 }
2263 // Node already exists
2264 else
2265 {
2266 el_pt->node_pt(new_node_number) = existing_node_pt;
2267 }
2268 }
2269 }
2270 }
2271 };
2272
2273
2274 /// ///////////////////////////////////////////////////////////////////////
2275 /// ///////////////////////////////////////////////////////////////////////
2276 /// ///////////////////////////////////////////////////////////////////////
2277
2278
2279 //=========================================================================
2280 // Refineable Gmsh-based Tet Mesh
2281 //=========================================================================
2282 template<class ELEMENT>
2283 class RefineableGmshTetMesh : public virtual GmshTetMesh<ELEMENT>,
2284 public virtual RefineableTetMeshBase
2285 {
2286 public:
2287 /// Constructor. If boolean is set
2288 /// to true, the target element sizes are read from file (used during
2289 /// adaptation; otherwise uniform target size is used).
2291 GmshParameters* gmsh_parameters_pt,
2292 const bool& use_mesh_grading_from_file,
2293 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2294 : GmshTetMesh<ELEMENT>(
2295 gmsh_parameters_pt, use_mesh_grading_from_file, time_stepper_pt)
2296 {
2298 }
2299
2300 /// Constructor
2302 GmshParameters* gmsh_parameters_pt,
2303 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2304 : GmshTetMesh<ELEMENT>(gmsh_parameters_pt,
2305 false, // Don't read mesh size data from file
2306 // unless explicitly requested.
2307 time_stepper_pt)
2308 {
2310 }
2311
2312
2313 /// Adapt mesh, based on elemental error provided
2314 void adapt(const Vector<double>& elem_error);
2315
2316 /// Refine uniformly
2318 {
2319 // hierher do it!
2320 std::string error_msg("Not written yet...");
2321 throw OomphLibError(
2322 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2323 }
2324
2325 /// Unrefine uniformly
2327 {
2328 // hierher do it!
2329 std::string error_msg("Not written yet...");
2330 throw OomphLibError(
2331 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2332 }
2333
2334
2335 protected:
2336 /// Helper function to initialise data associated with adaptation
2338 {
2339 // Set max and min targets for adaptation
2344 }
2345 };
2346
2347
2348 /// //////////////////////////////////////////////////////////////////////////
2349 /// //////////////////////////////////////////////////////////////////////////
2350 /// //////////////////////////////////////////////////////////////////////////
2351
2352
2353} // namespace oomph
2354
2355#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
Information for documentation of results: Directory and file number to enable output in the form RESL...
A general Finite Element class.
Definition: elements.h:1313
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Definition: elements.h:1842
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition: elements.cc:4290
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
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
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
Class to collate parameters for Gmsh mesh generation.
Vector< TetMeshFacetedSurface * > & internal_surface_pt()
Internal boundaries.
std::string & stem_for_filename_gmsh_size_transfer()
Stem for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty str...
double Dx_for_target_size_transfer
(Isotropic) grid spacing for target size transfer
double & element_volume()
Uniform target element volume.
TetMeshFacetedClosedSurface * Outer_boundary_pt
Pointer to outer boundary.
unsigned & gmsh_onscreen_output_counter()
Counter for marker that indicates where we are in gmsh on-screen output.
std::string & geo_and_msh_file_stem()
Stem for geo and msh files (input/output to command-line gmsh invocation)
std::string Geo_and_msh_file_stem
Stem for geo and msh files (input/output to command-line gmsh invocation)
std::string & gmsh_onscreen_output_file_name()
Output filename for gmsh on-screen output.
std::string Gmsh_command_line_invocation
Gmsh command line invocation string.
double & max_permitted_edge_ratio()
Max. permitted edge ratio.
double Max_element_size
Max. element size during refinement.
GmshParameters(TetMeshFacetedClosedSurface *const &outer_boundary_pt, const std::string &gmsh_command_line_invocation)
Specify outer boundary of domain to be meshed. Other parameters get default values and can be set via...
unsigned Max_sample_points_for_limited_locate_zeta_during_target_size_transfer
Target size is transferred onto regular grid (for gmsh) by locate zeta. We try to find the exact poin...
std::string & target_size_file_name()
Filename for target volumes (for system-call based transfer to gmsh during mesh adaptation)....
int & counter_for_filename_gmsh_size_transfer()
Counter for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty ...
double Min_element_size
Min. element size during refinement.
std::string & gmsh_command_line_invocation()
String to be issued via system command to activate gmsh.
TetMeshFacetedClosedSurface *& outer_boundary_pt()
Outer boundary.
unsigned & max_sample_points_for_limited_locate_zeta_during_target_size_transfer()
Target size is transferred onto regular grid (for gmsh) by locate zeta. We try to find the exact poin...
double & max_element_size()
Max. element size during refinement.
double Element_volume
Uniform element volume.
Vector< TetMeshFacetedSurface * > Internal_surface_pt
Internal boundaries.
unsigned Gmsh_onscreen_output_counter
Counter for marker that indicates where we are in gmsh on-screen output.
double & dx_for_target_size_transfer()
(Isotropic) grid spacing for target size transfer
double Max_permitted_edge_ratio
Max edge ratio before remesh gets triggered.
bool projection_is_disabled()
Is projection of old solution onto new mesh disabled?
std::string Stem_for_filename_gmsh_size_transfer
Stem for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty str...
void disable_projection()
Disable projection of old solution onto new mesh.
std::string Target_size_file_name
Filename for target volume (for system-call based transfer to gmsh during mesh adaptation)
double & min_element_size()
Min. element size during refinement.
bool Projection_is_disabled
Is projection of old solution onto new mesh disabled?
std::string Gmsh_onscreen_output_file_name
Output filename for gmsh on-screen output.
int Counter_for_filename_gmsh_size_transfer
Counter for filename used to doc target element sizes on gmsh grid. No doc if stem is equal to empty ...
void enable_projection()
Disable projection of old solution onto new mesh.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
GmshTetMesh(GmshParameters *gmsh_parameters_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor.
GmshParameters * Gmsh_parameters_pt
Parameters.
void build_from_scaffold(GmshTetScaffoldMesh *tmp_scaffold_mesh_pt, TimeStepper *time_stepper_pt)
Build unstructured tet gmesh mesh based on output from scaffold.
GmshTetMesh(GmshParameters *gmsh_parameters_pt, const bool &use_mesh_grading_from_file, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. If boolean is set to true, the target element sizes are read from file (used during adap...
void build_it(TimeStepper *time_stepper_pt, const bool &use_mesh_grading_from_file)
void create_mesh_from_msh_file()
Create mesh from msh file (created internally via disk-based operations)
GmshTetScaffoldMesh(GmshParameters *gmsh_parameters_pt, const bool &use_mesh_grading_from_file)
Build mesh, based on specified parameters. If boolean is set to true, the target element sizes are re...
void write_geo_file(const bool &use_mesh_grading_from_file)
Write geo file for gmsh.
GmshParameters * Gmsh_parameters_pt
Parameters.
A general mesh class.
Definition: mesh.h:67
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,...
Definition: mesh.h:91
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
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition: mesh.h:95
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
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
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
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1365
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
void refine_uniformly(DocInfo &doc_info)
Unrefine uniformly.
unsigned unrefine_uniformly()
Refine uniformly.
RefineableGmshTetMesh(GmshParameters *gmsh_parameters_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor.
void adapt(const Vector< double > &elem_error)
Adapt mesh, based on elemental error provided.
RefineableGmshTetMesh(GmshParameters *gmsh_parameters_pt, const bool &use_mesh_grading_from_file, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. If boolean is set to true, the target element sizes are read from file (used during adap...
void initialise_adaptation_data()
Helper function to initialise data associated with adaptation.
DocInfo doc_info()
Access fct for DocInfo.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
double Max_element_size
Max permitted element size.
double Min_element_size
Min permitted element size.
double Max_permitted_edge_ratio
Max edge ratio before remesh gets triggered.
General TElement class.
Definition: Telements.h:1208
////////////////////////////////////////////////////////////////// //////////////////////////////////...
bool operator==(const TetEdge &tet_edge) const
Comparison operator: Edges are identical if their sorted (and therefore possibly reversed) vertex ids...
std::pair< unsigned, unsigned > Vertex_pair
The vertices (sorted by vertex ids)
unsigned second_vertex_id() const
Second vertex id.
bool is_reversed() const
Edge is reversed in the sense that vertex1 actually has a higher id than vertex2 (when specified in t...
TetEdge(const unsigned &vertex1, const unsigned &vertex2)
Constructor: Pass two vertices, identified by their indices Edge "direction" is from lower vertex to ...
bool operator<(const TetEdge &tet_edge) const
Comparison operator. Lexicographic comparison based on vertex ids.
unsigned first_vertex_id() const
First vertex id.
bool Reversed
Is it reversed? I.e. is the first input vertex stored after the second one?
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: tet_mesh.h:661
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
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
TimeStepper * Time_stepper_pt
Timestepper used to build nodes.
Definition: tet_mesh.h:1057
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
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
TetMeshFacetedClosedSurface * Outer_boundary_pt
Faceted surface that defines outer boundaries.
Definition: tet_mesh.h:1038
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
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
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:168
bool facet_is_embedded_in_a_specified_region()
Boolean indicating that facet is embedded in a specified region.
Definition: tet_mesh.h:243
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:181
TetMeshVertex * vertex_pt(const unsigned &j) const
Pointer to j-th vertex (const)
Definition: tet_mesh.h:187
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
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:538
bool faceted_volume_represents_hole_for_gmsh() const
Does closed surface represent hole for gmsh?
Definition: tet_mesh.h:562
TetMeshFacet * facet_pt(const unsigned &j) const
Pointer to j-th facet.
Definition: tet_mesh.h:373
TetMeshVertex * vertex_pt(const unsigned &j) const
Pointer to j-th vertex.
Definition: tet_mesh.h:379
unsigned nfacet() const
Number of facets.
Definition: tet_mesh.h:325
unsigned one_based_facet_boundary_id(const unsigned &j) const
One-based boundary id of j-th facet.
Definition: tet_mesh.h:331
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:319
Vertex for Tet mesh generation. Can lie on multiple boundaries (identified via one-based enumeration!...
Definition: tet_mesh.h:50
double x(const unsigned &i) const
i-th coordinate
Definition: tet_mesh.h:124
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...