27#ifndef OOMPH_GMSH_TET_MESH_HEADER
28#define OOMPH_GMSH_TET_MESH_HEADER
32#include <oomph-lib-config.h>
40#include "../generic/sample_point_parameters.h"
41#include "../generic/mesh_as_geometric_object.h"
42#include "../generic/projection.h"
286 "Muppet! You can't build an edge from one vertex to the same vertex!",
360 template<
class ELEMENT>
371 template<
class ELEMENT>
381 double t_start = TimingHelpers::timer();
387 << TimingHelpers::timer() -
t_start <<
" sec " << std::endl;
388 t_start = TimingHelpers::timer();
393 std::string gmsh_onscreen_output_file_name =
397 if (gmsh_onscreen_output_file_name !=
"")
399 marker <<
"\n\n====================================================\n"
400 <<
" gmsh invocation: "
403 <<
"====================================================\n\n\n"
416 if (gmsh_onscreen_output_file_name !=
"")
427 << TimingHelpers::timer() -
t_start <<
" sec " << std::endl;
428 t_start = TimingHelpers::timer();
434 oomph_info <<
"Time for creating mesh from msh file: "
435 << TimingHelpers::timer() -
t_start <<
" sec " << std::endl;
469 std::string
error_msg(
"Failed to open mesh file: ");
480 if (
line !=
"$MeshFormat")
483 "First line has to contain the string \"$MeshFormat\"; ");
502 if (
line ==
"$Nodes")
510 for (
unsigned j = 0;
j <
nnod;
j++)
518 std::istringstream
iss(
s);
527 if (
line !=
"$EndNodes")
529 std::string
error_msg(
"Line has to contain the string \"$EndNodes\"; ");
540 if (
line !=
"$MeshFormat")
543 "First line has to contain the string \"$MeshFormat\"; ");
570 if (
line ==
"$Elements")
583 for (
unsigned e = 0;
e <
nel;
e++)
624 std::string
error_msg(
"Can't handle element type: ");
644 for (
unsigned i = 2;
i <
ntags;
i++)
656 std::istringstream
iss(
s);
697 if (
line !=
"$EndElements")
700 "Line has to contain the string \"$EndElements\"; ");
711 if (
line !=
"$MeshFormat")
714 "First line has to contain the string \"$MeshFormat\"; ");
737 unsigned nnod = ((*it).second).
size();
741 for (
unsigned j = 0;
j <
nnod;
j++)
744 std::map<unsigned, std::set<unsigned>>::iterator
itt =
748 for (std::set<unsigned>::iterator
ittt = (
itt->second).begin();
757 for (std::map<unsigned, unsigned>::iterator
itt =
762 if ((*itt).second == 3)
796 unsigned el_nnod = (*it).second.size();
840 for (std::set<unsigned>::iterator
it =
851 for (
unsigned i = 0;
i < 3;
i++)
868 for (std::map<
unsigned, std::set<unsigned>>::iterator
it =
877 for (std::set<unsigned>::iterator
itt = (*it).second.begin();
878 itt != (*it).second.end();
908 unsigned el_id = (*it).first;
911 itt != (*it).second.end();
917 unsigned n = x.size();
918 for (
unsigned i = 0;
i <
n;
i++)
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)
928 prefix =
", N=3, E=1, F=FEPOINT, ET=TRIANGLE";
942 for (std::map<
unsigned, std::set<unsigned>>::iterator
it =
949 <<
"\"" << std::endl;
950 for (std::set<unsigned>::iterator
itt = (*it).second.begin();
951 itt != (*it).second.end();
956 unsigned n = x.size();
957 for (
unsigned i = 0;
i <
n;
i++)
969 for (std::map<
unsigned, std::set<unsigned>>::iterator
it =
979 for (std::set<unsigned>::iterator
itt = (*it).second.begin();
980 itt != (*it).second.end();
984 <<
"\", N=4, E=1, F=FEPOINT, ET=TETRAHEDRON" << std::endl;
987 for (
unsigned j = 0;
j <
nnod;
j++)
991 unsigned n = x.size();
992 for (
unsigned i = 0;
i <
n;
i++)
998 outfile <<
"1 2 3 4" << std::endl;
1008 for (
unsigned e = 0;
e <
nel;
e++)
1012 oomph_info <<
"Total volume of all elements in scaffold mesh: " <<
vol
1031 std::string& target_size_file_name =
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;
1049 geo_file <<
"// Outer box" << std::endl;
1050 geo_file <<
"//==========" << std::endl;
1052 geo_file <<
"// Vertices" << std::endl;
1053 geo_file <<
"//---------" << std::endl;
1055 unsigned nv = outer_boundary_pt->nvertex();
1056 for (
unsigned j = 0;
j <
nv;
j++)
1061 for (
unsigned i = 0;
i < 3;
i++)
1071 unsigned nfacet = outer_boundary_pt->nfacet();
1076 for (
unsigned j = 0;
j <
nv - 1;
j++)
1093 geo_file <<
"// Edge of outer box" << std::endl;
1094 geo_file <<
"//------------------" << std::endl;
1096 std::map<TetEdge, unsigned>
tet_edge;
1102 geo_file <<
"Line(" <<
count + 1 <<
")={" << (*it).first_vertex_id()
1103 <<
"," << (*it).second_vertex_id() <<
"};" << std::endl;
1110 geo_file <<
"// Faces of outer box" << std::endl;
1111 geo_file <<
"//-------------------" << std::endl;
1114 geo_file <<
"Line Loop(" <<
f + 1 <<
")={";
1118 for (
unsigned j = 0;
j <
nv - 1;
j++)
1148 geo_file << ((
it->second) + 1) <<
"};" << std::endl;
1150 geo_file <<
"Plane Surface(" <<
f + 1 <<
")={" <<
f + 1 <<
"};"
1156 geo_file <<
"// Define Plane Surfaces bounding the volume" << std::endl;
1157 geo_file <<
"//------------------------------------------" << std::endl;
1159 for (
unsigned f = 0;
f <
nfacet - 1;
f++)
1167 geo_file <<
"// Define one-based boundary IDs" << std::endl;
1168 geo_file <<
"//------------------------------" << std::endl;
1172 outer_boundary_pt->one_based_facet_boundary_id(
f);
1174 <<
f + 1 <<
"};" << std::endl;
1189 std::map<unsigned, Vector<unsigned>>
1194 unsigned n_internal = internal_surface_pt.size();
1214 geo_file <<
"//==========================" << std::endl;
1216 geo_file <<
"// Vertices" << std::endl;
1217 geo_file <<
"//---------" << std::endl;
1219 unsigned nv = internal_surface_pt[
i_internal]->nvertex();
1220 for (
unsigned j = 0;
j <
nv;
j++)
1226 for (
unsigned i = 0;
i < 3;
i++)
1240 for (
unsigned j = 0;
j <
nv - 1;
j++)
1256 geo_file <<
"// Edge of inner faceted surface" << std::endl;
1257 geo_file <<
"//------------------------------" << std::endl;
1259 std::map<TetEdge, unsigned>
tet_edge;
1266 << (*it).first_vertex_id() <<
"," << (*it).second_vertex_id()
1267 <<
"};" << std::endl;
1274 geo_file <<
"// Faces of inner faceted surface" << std::endl;
1275 geo_file <<
"//-------------------------------" << std::endl;
1282 for (
unsigned j = 0;
j <
nv - 1;
j++)
1289 std::map<TetEdge, unsigned>::iterator
it =
1313 geo_file << ((
it->second) + 1) <<
"};" << std::endl;
1321 facet_pt->facet_is_embedded_in_a_specified_region();
1325 facet_pt->one_based_region_that_facet_is_embedded_in();
1367 facet_pt->one_based_adjacent_region_id());
1377 for (std::set<unsigned>::iterator
it =
region_id.begin();
1396 <<
"Something fishy going on! "
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";
1414 geo_file <<
"// Define Plane Surfaces bounding compound volume"
1416 <<
"//-----------------------------------------------"
1418 geo_file <<
"// that will have to be treated as hole in main volume"
1420 <<
"//-----------------------------------------------"
1425 for (
unsigned f = 0;
f <
n - 1;
f++)
1446 geo_file <<
"// Define Plane Surfaces bounding the region volume "
1447 << (*it).first << std::endl;
1448 geo_file <<
"//----------------------------------------------------"
1453 unsigned n = (*it).second.size();
1454 for (
unsigned f = 0;
f <
n - 1;
f++)
1458 geo_file << ((*it).second)[
n - 1] <<
"};" << std::endl;
1463 <<
" as the volume bounded by Surface Loop "
1467 <<
"//--------------------------------------------------------"
1473 <<
"};" << std::endl;
1484 if (
closed_srf_pt->faceted_volume_represents_hole_for_gmsh())
1491 geo_file <<
"// Define one-based region IDs" << std::endl;
1492 geo_file <<
"//----------------------------" << std::endl;
1493 geo_file <<
"Physical Volume(" << (*it).first <<
")={"
1495 <<
"};" << std::endl;
1504 <<
" embedded surfaces\n";
1527 geo_file <<
"// Define one-based boundary IDs" << std::endl;
1528 geo_file <<
"//------------------------------" << std::endl;
1532 internal_surface_pt[
i_internal]->one_based_facet_boundary_id(
f);
1548 geo_file <<
"// Define volume 1 as the volume bounded by Surface Loop 1"
1550 geo_file <<
"//--------------------------------------------------------"
1556 for (
unsigned i = 0;
i <
n;
i++)
1560 geo_file <<
"removed." << std::endl;
1561 geo_file <<
"//--------------------------------------------------------"
1569 for (
unsigned i = 0;
i <
n;
i++)
1577 geo_file <<
"// Define one-based region IDs" << std::endl;
1578 geo_file <<
"//----------------------------" << std::endl;
1580 <<
")={1};" << std::endl;
1588 geo_file <<
"// Embed Plane Surfaces in main volume (volume 1)"
1590 geo_file <<
"//-----------------------------------------------"
1593 for (
unsigned s = 0;
s <
ns - 1;
s++)
1598 <<
"} In Volume{1};" << std::endl;
1610 std::ifstream
file(target_size_file_name.c_str(), std::ios_base::in);
1613 if (!
file.is_open())
1615 std::string
error_msg(
"Failed to open target volume file: ");
1616 error_msg +=
"\"" + target_size_file_name;
1622 geo_file <<
"Field[1]=Structured;" << std::endl;
1623 geo_file <<
"Field[1].FileName=\"" << target_size_file_name <<
"\";"
1625 geo_file <<
"Field[1].TextFormat=1;" << std::endl;
1626 geo_file <<
"Background Field = 1;" << std::endl;
1632 geo_file <<
"Mesh 3;" << std::endl;
1660 template<
class ELEMENT>
1691 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
1704 Outer_boundary_pt = outer_boundary_pt;
1708 unsigned n_facet = Outer_boundary_pt->nfacet();
1711 unsigned b = Outer_boundary_pt->one_based_facet_boundary_id(
f);
1720 error_message <<
"Boundary IDs have to be one-based. Yours is " <<
b
1731 Internal_surface_pt = internal_surface_pt;
1735 unsigned n = Internal_surface_pt.size();
1736 for (
unsigned i = 0;
i <
n;
i++)
1738 unsigned n_facet = Internal_surface_pt[
i]->nfacet();
1741 unsigned b = Internal_surface_pt[
i]->one_based_facet_boundary_id(
f);
1750 error_message <<
"Boundary IDs have to be one-based. Yours is "
1773 for (
unsigned b = 0;
b <
nb;
b++)
1796 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
1820 for (
unsigned e = 0;
e <
nelem;
e++)
1835 for (
unsigned e = 0;
e <
nelem;
e++)
1916 for (
unsigned f = 0;
f < 4;
f++)
1959 std::set<unsigned>*
bc0_pt;
1963 std::set<unsigned>*
bc1_pt;
1967 std::set<unsigned>*
bc2_pt;
1972 std::set_intersection(
1980 std::set_intersection(
1987 for (std::set<unsigned>::iterator
it =
2005 for (
unsigned i = 0;
i <
nr;
i++)
2010 for (
unsigned e = 0;
e <
nel;
e++)
2021 for (
unsigned f = 0;
f < 4;
f++)
2063 std::set<unsigned>*
bc0_pt;
2067 std::set<unsigned>*
bc1_pt;
2071 std::set<unsigned>*
bc2_pt;
2076 std::set_intersection(
2084 std::set_intersection(
2091 for (std::set<unsigned>::iterator
it =
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";
2138 for (
unsigned e = 0;
e <
nelem;
e++)
2146 for (
unsigned j = 0;
j < 6; ++
j)
2151 std::pair<Node*, Node*>
edge;
2207 std::set<unsigned>*
bc0_pt;
2211 std::set<unsigned>*
bc1_pt;
2215 std::set_intersection(
2238 for (std::set<unsigned>::iterator
it =
2252 for (
unsigned i = 0;
i < 3;
i++)
2282 template<
class ELEMENT>
2284 public virtual RefineableTetMeshBase
2320 std::string
error_msg(
"Not written yet...");
2329 std::string
error_msg(
"Not written yet...");
2342 this->Max_permitted_edge_ratio =
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.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
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.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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?
////////////////////////////////////////////////////////////////////// //////////////////////////////...