22#include <OMExceptions.H>
36 MeshPair(
const Mesh& m1,
const Mesh& m2,
const int o): meshes{&m1,&m2},orientation(o) { }
44 const Mesh* meshes[2];
51 typedef std::vector<std::vector<const Mesh*>>
MeshParts;
53 typedef std::vector<std::pair<std::string,std::string>>
MeshList;
63 Geometry(
const std::string& geomFileName,
const bool OLD_ORDERING=
false) {
64 load(geomFileName,OLD_ORDERING);
67 Geometry(
const std::string& geomFileName,
const std::string& condFileName,
const bool OLD_ORDERING=
false) {
68 load(geomFileName,condFileName,OLD_ORDERING);
73 Geometry(
const char* geomFileName,
const bool OLD_ORDERING=
false):
Geometry(std::string(geomFileName),OLD_ORDERING) { }
74 Geometry(
const char* geomFileName,
const char* condFileName,
const bool OLD_ORDERING=
false):
75 Geometry(std::string(geomFileName),std::string(condFileName),OLD_ORDERING) { }
77 void info(
const bool verbose=
false)
const;
79 for (
const auto& domain : domains())
80 if (!domain.has_conductivity())
103 const Vertices::iterator vit = std::find(vertices().begin(),vertices().end(),V);
104 if (vit!=vertices().end())
105 return vit-vertices().begin();
107 vertices().push_back(V);
108 return vertices().size()-1;
116 meshes().emplace_back(
this);
117 Mesh& mesh = meshes().back();
124 for (
unsigned i=0; i<vs.size(); ++i)
125 indmap.insert({ i, add_vertex(vs[i]) });
154 for (
const auto& domain : domains())
155 if (domain.contains(m))
156 result.push_back(&domain);
168 outer_domain = &domain;
170 boundary.interface().set_to_outermost();
182 double sigma (
const Mesh& m1,
const Mesh& m2)
const {
return eval_on_common_domains<IDENTITY>(m1,m2); }
183 double sigma_inv(
const Mesh& m1,
const Mesh& m2)
const {
return eval_on_common_domains<INVERSE>(m1,m2); }
184 double indicator(
const Mesh& m1,
const Mesh& m2)
const {
return eval_on_common_domains<INDICATOR>(m1,m2); }
191 for (
const auto& domainptr : doms)
192 res += domainptr->conductivity()*domainptr->mesh_orientation(m);
205 return (doms.size()==0) ? 0 : ((doms[0]->mesh_orientation(m1)==doms[0]->mesh_orientation(m2)) ? 1 : -1);
211 void load(
const std::string& filename,
const bool OLD_ORDERING=
false) {
213 read_geometry_file(filename);
214 finalize(OLD_ORDERING);
217 void load(
const std::string& geomFileName,
const std::string& condFileName,
const bool OLD_ORDERING=
false) {
219 read_geometry_file(geomFileName);
220 read_conductivity_file(condFileName);
221 finalize(OLD_ORDERING);
226 void save(
const std::string& filename)
const;
235 if (has_conductivities())
236 mark_current_barriers();
238 if (domains().size()!=0) {
239 set_outermost_domain(outermost_domain());
240 check_geometry_is_nested();
243 generate_indices(OLD_ORDERING);
246 for (
const auto& mesh : meshes())
247 mesh.check_consistency(
"geometry finalize step");
264 geom_vertices.clear();
266 geom_domains.clear();
272 void read_geometry_file(
const std::string& filename);
273 void read_conductivity_file(
const std::string& filename);
275 void make_mesh_pairs();
283 const Domain* outer_domain = 0;
285 size_t num_params = 0;
287 void generate_indices(
const bool);
289 DomainsReference common_domains(
const Mesh& m1,
const Mesh& m2)
const {
290 const DomainsReference& doms1 = domains(m1);
291 const DomainsReference& doms2 = domains(m2);
292 DomainsReference doms;
293 std::set_intersection(doms1.begin(),doms1.end(),doms2.begin(),doms2.end(),std::back_inserter(doms));
299 static double IDENTITY(
const Domain& domain) {
return domain.conductivity(); }
300 static double INVERSE(
const Domain& domain) {
return 1.0/domain.conductivity(); }
301 static double INDICATOR(
const Domain&) {
return 1.0; }
303 template <
double Function(const Domain&)>
304 double eval_on_common_domains(
const Mesh& m1,
const Mesh& m2)
const {
305 const DomainsReference& doms = common_domains(m1,m2);
307 for (
const auto& domainptr : doms)
308 result += Function(*domainptr);
314 std::set<Vertex> invalid_vertices_;
315 size_t nb_current_barrier_triangles_ = 0;
317 MeshParts independant_parts;
a Domain is a vector of SimpleDomain A Domain is the intersection of simple domains (of type SimpleDo...
Boundaries & boundaries()
Boundaries of the domain.
Geometry contains the electrophysiological model Vertices, meshes and domains are stored in this geom...
Domains & domains()
Return the list of domains.
size_t nb_invalid_vertices()
Meshes & meshes()
Return the list of meshes involved in the geometry.
const Interface & interface(const std::string &name) const
returns the Interface called
void finalize(const bool OLD_ORDERING=false)
void info(const bool verbose=false) const
Print information on the geometry.
int relative_orientation(const Mesh &m1, const Mesh &m2) const
Give the relative orientation of two meshes:
std::vector< const Domain * > DomainsReference
bool check(const Mesh &m) const
check if m intersect geometry meshes
Vertices & vertices()
Return the list of vertices involved in the geometry.
const Domain & domain(const std::string &name) const
Get specific domains.
void check_geometry_is_nested()
Domain & outermost_domain()
Returns the outermost domain.
Mesh & add_mesh(const std::string &name="")
double sigma_inv(const Mesh &m1, const Mesh &m2) const
Geometry(const unsigned n)
unsigned add_vertex(const Vertex &V)
Add a vertex.
const Interface & outermost_interface() const
returns the outermost interface (only valid for nested geometries).
bool selfCheck() const
the geometry meshes intersect each other
void mark_current_barriers()
const Meshes & meshes() const
const Vertices & vertices() const
Geometry(const std::string &geomFileName, const bool OLD_ORDERING=false)
double conductivity_jump(const Mesh &m) const
Return the conductivity jump across a mesh (i.e. between the 2 domains it separates).
Geometry(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
const Domain & domain(const Vect3 &p) const
returns the Domain containing the point p
Geometry(const char *geomFileName, const char *condFileName, const bool OLD_ORDERING=false)
bool is_outermost(const Domain &domain) const
size_t & nb_current_barrier_triangles()
const MeshParts & isolated_parts() const
Geometry(const char *geomFileName, const bool OLD_ORDERING=false)
const MeshPairs & communicating_mesh_pairs() const
void set_outermost_domain(Domain &domain)
size_t nb_parameters() const
the total number of vertices + triangles
size_t nb_current_barrier_triangles() const
Handle multiple isolated domains.
std::vector< MeshPair > MeshPairs
const Interface & innermost_interface() const
returns the innermost interface (only valid for nested geometries).
DomainsReference domains(const Mesh &m) const
Return the list of domains containing a mesh.
bool check_inner(const Matrix &m) const
check if dipoles are outside of geometry meshes
Mesh & mesh(const std::string &name)
returns the Mesh called
IndexMap add_vertices(const Vertices &vs)
void load(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
double indicator(const Mesh &m1, const Mesh &m2) const
const Mesh & mesh(const std::string &id) const
void load(const std::string &filename, const bool OLD_ORDERING=false)
const Domains & domains() const
std::vector< std::pair< std::string, std::string > > MeshList
std::vector< std::vector< const Mesh * > > MeshParts
void save(const std::string &filename) const
double sigma(const Mesh &m1, const Mesh &m2) const
bool has_conductivities() const
Interface class An interface is a closed-shape composed of oriented meshes (vector of oriented meshes...
Matrix class Matrix class.
file containing the definition of a Domain.
std::vector< Vertex > Vertices
std::vector< Mesh > Meshes
std::vector< Domain > Domains
A vector of Domain is called Domains.
std::map< unsigned, unsigned > IndexMap
MeshPair(const Mesh &m1, const Mesh &m2, const int o)
const Mesh & operator()(const unsigned i) const
int relative_orientation() const