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;
59 Geometry(
const std::string& geomFileName,
const bool OLD_ORDERING=
false) {
60 load(geomFileName,OLD_ORDERING);
63 Geometry(
const std::string& geomFileName,
const std::string& condFileName,
const bool OLD_ORDERING=
false) {
64 load(geomFileName,condFileName,OLD_ORDERING);
69 Geometry(
const char* geomFileName,
const bool OLD_ORDERING=
false):
Geometry(std::string(geomFileName),OLD_ORDERING) { }
70 Geometry(
const char* geomFileName,
const char* condFileName,
const bool OLD_ORDERING=
false):
71 Geometry(std::string(geomFileName),std::string(condFileName),OLD_ORDERING) { }
73 void info(
const bool verbose=
false)
const;
75 for (
const auto& domain : domains())
76 if (!domain.has_conductivity())
98 const Vertices::iterator vit = std::find(vertices().begin(),vertices().end(),V);
99 if (vit!=vertices().end())
100 return vit-vertices().begin();
102 vertices().push_back(V);
103 return vertices().size()-1;
111 meshes().emplace_back(
this);
112 Mesh& mesh = meshes().back();
119 for (
unsigned i=0; i<vs.size(); ++i)
120 indmap.insert({ i, add_vertex(vs[i]) });
149 for (
const auto& domain : domains())
150 if (domain.contains(m))
151 result.push_back(&domain);
163 outer_domain = &domain;
165 boundary.interface().set_to_outermost();
177 double sigma (
const Mesh& m1,
const Mesh& m2)
const {
return eval_on_common_domains<IDENTITY>(m1,m2); }
178 double sigma_inv(
const Mesh& m1,
const Mesh& m2)
const {
return eval_on_common_domains<INVERSE>(m1,m2); }
179 double indicator(
const Mesh& m1,
const Mesh& m2)
const {
return eval_on_common_domains<INDICATOR>(m1,m2); }
186 for (
const auto& domainptr : doms)
187 res += domainptr->conductivity()*domainptr->mesh_orientation(m);
200 return (doms.size()==0) ? 0 : ((doms[0]->mesh_orientation(m1)==doms[0]->mesh_orientation(m2)) ? 1 : -1);
206 void load(
const std::string& filename,
const bool OLD_ORDERING=
false) {
208 read_geometry_file(filename);
209 finalize(OLD_ORDERING);
212 void load(
const std::string& geomFileName,
const std::string& condFileName,
const bool OLD_ORDERING=
false) {
214 read_geometry_file(geomFileName);
215 read_conductivity_file(condFileName);
216 finalize(OLD_ORDERING);
221 void save(
const std::string& filename)
const;
230 if (has_conductivities())
231 mark_current_barriers();
233 if (domains().size()!=0) {
234 set_outermost_domain(outermost_domain());
235 check_geometry_is_nested();
238 generate_indices(OLD_ORDERING);
255 geom_vertices.clear();
257 geom_domains.clear();
263 void read_geometry_file(
const std::string& filename);
264 void read_conductivity_file(
const std::string& filename);
266 void make_mesh_pairs();
274 const Domain* outer_domain = 0;
276 size_t num_params = 0;
278 void generate_indices(
const bool);
280 DomainsReference common_domains(
const Mesh& m1,
const Mesh& m2)
const {
281 const DomainsReference& doms1 = domains(m1);
282 const DomainsReference& doms2 = domains(m2);
283 DomainsReference doms;
284 std::set_intersection(doms1.begin(),doms1.end(),doms2.begin(),doms2.end(),std::back_inserter(doms));
290 static double IDENTITY(
const Domain& domain) {
return domain.conductivity(); }
291 static double INVERSE(
const Domain& domain) {
return 1.0/domain.conductivity(); }
292 static double INDICATOR(
const Domain&) {
return 1.0; }
294 template <
double Function(const Domain&)>
295 double eval_on_common_domains(
const Mesh& m1,
const Mesh& m2)
const {
296 const DomainsReference& doms = common_domains(m1,m2);
298 for (
const auto& domainptr : doms)
299 result += Function(*domainptr);
305 std::set<Vertex> invalid_vertices_;
306 size_t nb_current_barrier_triangles_ = 0;
308 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.
std::vector< const Domain * > DomainsReference
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:
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()
std::vector< MeshPair > MeshPairs
Domain & outermost_domain()
Returns the outermost domain.
Mesh & add_mesh(const std::string &name="")
double sigma_inv(const Mesh &m1, const Mesh &m2) const
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
std::vector< std::pair< std::string, std::string > > MeshList
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.
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::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