OpenMEEG
geometry.h
Go to the documentation of this file.
1 // Project Name: OpenMEEG (http://openmeeg.github.io)
2 // © INRIA and ENPC under the French open source license CeCILL-B.
3 // See full copyright notice in the file LICENSE.txt
4 // If you make a copy of this file, you must either:
5 // - provide also LICENSE.txt and modify this header to refer to it.
6 // - replace this header by the LICENSE.txt content.
7 
8 #pragma once
9 
10 #include <iterator>
11 #include <string>
12 #include <vector>
13 #include <set>
14 
15 #include <om_common.h>
16 #include <vertex.h>
17 #include <triangle.h>
18 #include <interface.h>
19 #include <domain.h>
20 #include <matrix.h>
21 
22 #include <OMExceptions.H>
23 
24 namespace OpenMEEG {
25 
26  class Mesh;
27 
30 
31  class OPENMEEG_EXPORT Geometry {
32  public:
33 
34  struct MeshPair {
35 
36  MeshPair(const Mesh& m1,const Mesh& m2,const int o): meshes{&m1,&m2},orientation(o) { }
37 
38  const Mesh& operator()(const unsigned i) const { return *meshes[i]; }
39 
40  int relative_orientation() const { return orientation; }
41 
42  private:
43 
44  const Mesh* meshes[2];
45  int orientation;
46  };
47 
48  typedef std::vector<MeshPair> MeshPairs;
49 
50  typedef std::vector<const Domain*> DomainsReference;
51  typedef std::vector<std::vector<const Mesh*>> MeshParts;
52 
53  typedef std::vector<std::pair<std::string,std::string>> MeshList;
54 
56 
57  Geometry() {}
58 
59  Geometry(const std::string& geomFileName,const bool OLD_ORDERING=false) {
60  load(geomFileName,OLD_ORDERING);
61  }
62 
63  Geometry(const std::string& geomFileName,const std::string& condFileName,const bool OLD_ORDERING=false) {
64  load(geomFileName,condFileName,OLD_ORDERING);
65  }
66 
67  // Absolutely necessary or wrong constructor is called because of conversion of char* to bool.
68 
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) { }
72 
73  void info(const bool verbose=false) const;
74  bool has_conductivities() const {
75  for (const auto& domain : domains())
76  if (!domain.has_conductivity())
77  return false;
78  return true;
79  }
80  bool selfCheck() const;
81  bool check(const Mesh& m) const;
82  bool check_inner(const Matrix& m) const;
83 
85 
86  bool is_nested() const { return nested; }
87 
89 
90  Vertices& vertices() { return geom_vertices; }
91  const Vertices& vertices() const { return geom_vertices; }
92 
94 
95  unsigned add_vertex(const Vertex& V) {
96  // Insert the vertex in the set of vertices if it is not already in.
97 
98  const Vertices::iterator vit = std::find(vertices().begin(),vertices().end(),V);
99  if (vit!=vertices().end())
100  return vit-vertices().begin();
101 
102  vertices().push_back(V);
103  return vertices().size()-1;
104  }
105 
106  Mesh& add_mesh(const std::string& name="") {
107 
108  // It is dangerous to store the returned mesh because the vector can be reallocated.
109  // Use mesh(name) after all meshes have been added....
110 
111  meshes().emplace_back(this);
112  Mesh& mesh = meshes().back();
113  mesh.name() = name;
114  return mesh;
115  }
116 
118  IndexMap indmap;
119  for (unsigned i=0; i<vs.size(); ++i)
120  indmap.insert({ i, add_vertex(vs[i]) });
121  return indmap;
122  }
123 
125 
126  Meshes& meshes() { return geom_meshes; }
127  const Meshes& meshes() const { return geom_meshes; }
128 
129  const MeshPairs& communicating_mesh_pairs() const { return meshpairs; }
130 
132 
133  Mesh& mesh(const std::string& name);
134 
136 
137  Domains& domains() { return geom_domains; }
138  const Domains& domains() const { return geom_domains; }
139 
141 
142  const Domain& domain(const std::string& name) const;
143  const Domain& domain(const Vect3& p) const;
144 
146 
147  DomainsReference domains(const Mesh& m) const {
148  DomainsReference result;
149  for (const auto& domain : domains())
150  if (domain.contains(m))
151  result.push_back(&domain);
152  return result;
153  }
154 
155  size_t nb_parameters() const { return num_params; }
156 
158  // It is unclear whether outermost_domain and set_outermost_domain need to be in the public interface.
159 
161 
162  void set_outermost_domain(Domain& domain) {
163  outer_domain = &domain;
164  for (auto& boundary : domain.boundaries())
165  boundary.interface().set_to_outermost();
166  }
167 
168  bool is_outermost(const Domain& domain) const { return outer_domain==&domain; }
169 
172 
173  const Interface& interface(const std::string& name) const;
174 
175  // TODO: Find better names for the next two methods.
176 
177  double sigma (const Mesh& m1,const Mesh& m2) const { return eval_on_common_domains<IDENTITY>(m1,m2); } // return the (sum) conductivity(ies) of the shared domain(s).
178  double sigma_inv(const Mesh& m1,const Mesh& m2) const { return eval_on_common_domains<INVERSE>(m1,m2); } // return the (sum) inverse of conductivity(ies) of the shared domain(s).
179  double indicator(const Mesh& m1,const Mesh& m2) const { return eval_on_common_domains<INDICATOR>(m1,m2); } // return the (sum) indicator function of the shared domain(s).
180 
182 
183  double conductivity_jump(const Mesh& m) const {
184  const DomainsReference& doms = domains(m);
185  double res = 0.0;
186  for (const auto& domainptr : doms)
187  res += domainptr->conductivity()*domainptr->mesh_orientation(m);
188  return res;
189  }
190 
195 
196  int relative_orientation(const Mesh& m1,const Mesh& m2) const {
197  if (&m1==&m2) // Fast path for identical meshes.
198  return 1;
199  const DomainsReference& doms = common_domains(m1,m2); // 2 meshes have either 0, 1 or 2 domains in common
200  return (doms.size()==0) ? 0 : ((doms[0]->mesh_orientation(m1)==doms[0]->mesh_orientation(m2)) ? 1 : -1);
201  }
202 
203 
204  // Calling this method read induces failures due do wrong conversions when read is passed with one or two arguments...
205 
206  void load(const std::string& filename,const bool OLD_ORDERING=false) {
207  clear();
208  read_geometry_file(filename);
209  finalize(OLD_ORDERING);
210  }
211 
212  void load(const std::string& geomFileName,const std::string& condFileName,const bool OLD_ORDERING=false) {
213  clear();
214  read_geometry_file(geomFileName);
215  read_conductivity_file(condFileName);
216  finalize(OLD_ORDERING);
217  }
218 
219  void import(const MeshList& meshes);
220 
221  void save(const std::string& filename) const;
222 
223  void finalize(const bool OLD_ORDERING=false) {
224  // TODO: We should check the correct decomposition of the geometry into domains here.
225  // In a correct decomposition, each interface is used exactly once ?? Unsure...
226  // Search for the outermost domain and set boolean OUTERMOST on the domain in the vector domains.
227  // An outermost domain is defined as the only domain which has no inside. It is supposed to be
228  // unique.
229 
230  if (has_conductivities())
231  mark_current_barriers(); // mark meshes that touch the domains of null conductivity.
232 
233  if (domains().size()!=0) {
234  set_outermost_domain(outermost_domain());
235  check_geometry_is_nested();
236  }
237 
238  generate_indices(OLD_ORDERING);
239  make_mesh_pairs();
240  }
241 
243 
244  size_t nb_current_barrier_triangles() const { return nb_current_barrier_triangles_; }
245  size_t& nb_current_barrier_triangles() { return nb_current_barrier_triangles_; }
246  size_t nb_invalid_vertices() { return invalid_vertices_.size(); }
247 
248  const MeshParts& isolated_parts() const { return independant_parts; }
250  const Mesh& mesh(const std::string& id) const; // Is this useful ?? TODO.
251 
252  private:
253 
254  void clear() {
255  geom_vertices.clear();
256  geom_meshes.clear();
257  geom_domains.clear();
258  nested = false;
259  outer_domain = 0;
260  num_params = 0;
261  }
262 
263  void read_geometry_file(const std::string& filename);
264  void read_conductivity_file(const std::string& filename);
265 
266  void make_mesh_pairs();
267 
269 
270  Vertices geom_vertices;
271  Meshes geom_meshes;
272  Domains geom_domains;
273 
274  const Domain* outer_domain = 0;
275  bool nested = false;
276  size_t num_params = 0; // total number = nb of vertices + nb of triangles
277 
278  void generate_indices(const bool);
279 
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));
285  return doms;
286  }
287 
288  // Accumulate a function over the domain common to two meshes.
289 
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; }
293 
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);
297  double result = 0.0;
298  for (const auto& domainptr : doms)
299  result += Function(*domainptr);
300  return result;
301  }
302 
304 
305  std::set<Vertex> invalid_vertices_;
306  size_t nb_current_barrier_triangles_ = 0;
307 
308  MeshParts independant_parts;
309  MeshPairs meshpairs;
310  };
311 }
a Domain is a vector of SimpleDomain A Domain is the intersection of simple domains (of type SimpleDo...
Definition: domain.h:57
Boundaries & boundaries()
Boundaries of the domain.
Definition: domain.h:67
Geometry contains the electrophysiological model Vertices, meshes and domains are stored in this geom...
Definition: geometry.h:31
std::vector< const Domain * > DomainsReference
Definition: geometry.h:50
size_t nb_invalid_vertices()
Definition: geometry.h:246
void finalize(const bool OLD_ORDERING=false)
Definition: geometry.h:223
Mesh & mesh(const std::string &name)
returns the Mesh called
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:
Definition: geometry.h:196
bool check(const Mesh &m) const
check if m intersect geometry meshes
const Domain & domain(const std::string &name) const
Get specific domains.
const Interface & innermost_interface() const
returns the innermost interface (only valid for nested geometries).
const Interface & interface(const std::string &name) const
returns the Interface called
void check_geometry_is_nested()
std::vector< MeshPair > MeshPairs
Definition: geometry.h:48
Mesh & add_mesh(const std::string &name="")
Definition: geometry.h:106
double sigma_inv(const Mesh &m1, const Mesh &m2) const
Definition: geometry.h:178
Meshes & meshes()
Return the list of meshes involved in the geometry.
Definition: geometry.h:126
unsigned add_vertex(const Vertex &V)
Add a vertex.
Definition: geometry.h:95
bool selfCheck() const
the geometry meshes intersect each other
void mark_current_barriers()
std::vector< std::pair< std::string, std::string > > MeshList
Definition: geometry.h:53
const Domains & domains() const
Definition: geometry.h:138
Geometry(const std::string &geomFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:59
Domains & domains()
Return the list of domains.
Definition: geometry.h:137
Vertices & vertices()
Return the list of vertices involved in the geometry.
Definition: geometry.h:90
double conductivity_jump(const Mesh &m) const
Return the conductivity jump across a mesh (i.e. between the 2 domains it separates).
Definition: geometry.h:183
const Domain & domain(const Vect3 &p) const
returns the Domain containing the point p
Geometry(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:63
Geometry(const char *geomFileName, const char *condFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:70
const Vertices & vertices() const
Definition: geometry.h:91
bool is_outermost(const Domain &domain) const
Definition: geometry.h:168
Geometry()
Constructors.
Definition: geometry.h:57
const Mesh & mesh(const std::string &id) const
const MeshParts & isolated_parts() const
Definition: geometry.h:248
Geometry(const char *geomFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:69
void set_outermost_domain(Domain &domain)
Definition: geometry.h:162
bool is_nested() const
Definition: geometry.h:86
size_t nb_parameters() const
the total number of vertices + triangles
Definition: geometry.h:155
const MeshPairs & communicating_mesh_pairs() const
Definition: geometry.h:129
const Meshes & meshes() const
Definition: geometry.h:127
const Interface & outermost_interface() const
returns the outermost interface (only valid for nested geometries).
size_t nb_current_barrier_triangles() const
Handle multiple isolated domains.
Definition: geometry.h:244
DomainsReference domains(const Mesh &m) const
Return the list of domains containing a mesh.
Definition: geometry.h:147
bool check_inner(const Matrix &m) const
check if dipoles are outside of geometry meshes
IndexMap add_vertices(const Vertices &vs)
Definition: geometry.h:117
void load(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:212
size_t & nb_current_barrier_triangles()
Definition: geometry.h:245
double indicator(const Mesh &m1, const Mesh &m2) const
Definition: geometry.h:179
void load(const std::string &filename, const bool OLD_ORDERING=false)
Definition: geometry.h:206
std::vector< std::vector< const Mesh * > > MeshParts
Definition: geometry.h:51
void save(const std::string &filename) const
double sigma(const Mesh &m1, const Mesh &m2) const
Definition: geometry.h:177
Domain & outermost_domain()
Returns the outermost domain.
bool has_conductivities() const
Definition: geometry.h:74
Interface class An interface is a closed-shape composed of oriented meshes (vector of oriented meshes...
Definition: interface.h:49
Matrix class Matrix class.
Definition: matrix.h:28
std::string & name()
Definition: mesh.h:81
Vect3.
Definition: vect3.h:28
Vertex.
Definition: vertex.h:20
file containing the definition of a Domain.
std::vector< Vertex > Vertices
Definition: vertex.h:37
std::vector< Mesh > Meshes
Definition: mesh.h:272
std::vector< Domain > Domains
A vector of Domain is called Domains.
Definition: domain.h:112
std::map< unsigned, unsigned > IndexMap
Definition: triangle.h:141
MeshPair(const Mesh &m1, const Mesh &m2, const int o)
Definition: geometry.h:36
const Mesh & operator()(const unsigned i) const
Definition: geometry.h:38
int relative_orientation() const
Definition: geometry.h:40