OpenMEEG
Loading...
Searching...
No Matches
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
24namespace 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
58
59 Geometry(const unsigned n) {
60 meshes().reserve(n);
61 }
62
63 Geometry(const std::string& geomFileName,const bool OLD_ORDERING=false) {
64 load(geomFileName,OLD_ORDERING);
65 }
66
67 Geometry(const std::string& geomFileName,const std::string& condFileName,const bool OLD_ORDERING=false) {
68 load(geomFileName,condFileName,OLD_ORDERING);
69 }
70
71 // Absolutely necessary or wrong constructor is called because of conversion of char* to bool.
72
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) { }
76
77 void info(const bool verbose=false) const;
78 bool has_conductivities() const {
79 for (const auto& domain : domains())
80 if (!domain.has_conductivity())
81 return false;
82 return true;
83 }
84
85 bool selfCheck() const;
86 bool check(const Mesh& m) const;
87 bool check_inner(const Matrix& m) const;
88
90
91 bool is_nested() const { return nested; }
92
94
95 Vertices& vertices() { return geom_vertices; }
96 const Vertices& vertices() const { return geom_vertices; }
97
99
100 unsigned add_vertex(const Vertex& V) {
101 // Insert the vertex in the set of vertices if it is not already in.
102
103 const Vertices::iterator vit = std::find(vertices().begin(),vertices().end(),V);
104 if (vit!=vertices().end())
105 return vit-vertices().begin();
106
107 vertices().push_back(V);
108 return vertices().size()-1;
109 }
110
111 Mesh& add_mesh(const std::string& name="") {
112
113 // It is dangerous to store the returned mesh because the vector can be reallocated.
114 // Use mesh(name) after all meshes have been added....
115
116 meshes().emplace_back(this);
117 Mesh& mesh = meshes().back();
118 mesh.name() = name;
119 return mesh;
120 }
121
123 IndexMap indmap;
124 for (unsigned i=0; i<vs.size(); ++i)
125 indmap.insert({ i, add_vertex(vs[i]) });
126 return indmap;
127 }
128
130
131 Meshes& meshes() { return geom_meshes; }
132 const Meshes& meshes() const { return geom_meshes; }
133
134 const MeshPairs& communicating_mesh_pairs() const { return meshpairs; }
135
137
138 Mesh& mesh(const std::string& name);
139
141
142 Domains& domains() { return geom_domains; }
143 const Domains& domains() const { return geom_domains; }
144
146
147 const Domain& domain(const std::string& name) const;
148 const Domain& domain(const Vect3& p) const;
149
151
152 DomainsReference domains(const Mesh& m) const {
153 DomainsReference result;
154 for (const auto& domain : domains())
155 if (domain.contains(m))
156 result.push_back(&domain);
157 return result;
158 }
159
160 size_t nb_parameters() const { return num_params; }
161
163 // It is unclear whether outermost_domain and set_outermost_domain need to be in the public interface.
164
166
168 outer_domain = &domain;
169 for (auto& boundary : domain.boundaries())
170 boundary.interface().set_to_outermost();
171 }
172
173 bool is_outermost(const Domain& domain) const { return outer_domain==&domain; }
174
177
178 const Interface& interface(const std::string& name) const;
179
180 // TODO: Find better names for the next two methods.
181
182 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).
183 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).
184 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).
185
187
188 double conductivity_jump(const Mesh& m) const {
189 const DomainsReference& doms = domains(m);
190 double res = 0.0;
191 for (const auto& domainptr : doms)
192 res += domainptr->conductivity()*domainptr->mesh_orientation(m);
193 return res;
194 }
195
200
201 int relative_orientation(const Mesh& m1,const Mesh& m2) const {
202 if (&m1==&m2) // Fast path for identical meshes.
203 return 1;
204 const DomainsReference& doms = common_domains(m1,m2); // 2 meshes have either 0, 1 or 2 domains in common
205 return (doms.size()==0) ? 0 : ((doms[0]->mesh_orientation(m1)==doms[0]->mesh_orientation(m2)) ? 1 : -1);
206 }
207
208
209 // Calling this method read induces failures due do wrong conversions when read is passed with one or two arguments...
210
211 void load(const std::string& filename,const bool OLD_ORDERING=false) {
212 clear();
213 read_geometry_file(filename);
214 finalize(OLD_ORDERING);
215 }
216
217 void load(const std::string& geomFileName,const std::string& condFileName,const bool OLD_ORDERING=false) {
218 clear();
219 read_geometry_file(geomFileName);
220 read_conductivity_file(condFileName);
221 finalize(OLD_ORDERING);
222 }
223
224 void import(const MeshList& meshes);
225
226 void save(const std::string& filename) const;
227
228 void finalize(const bool OLD_ORDERING=false) {
229 // TODO: We should check the correct decomposition of the geometry into domains here.
230 // In a correct decomposition, each interface is used exactly once ?? Unsure...
231 // Search for the outermost domain and set boolean OUTERMOST on the domain in the vector domains.
232 // An outermost domain is defined as the only domain which has no inside. It is supposed to be
233 // unique.
234
235 if (has_conductivities())
236 mark_current_barriers(); // mark meshes that touch the domains of null conductivity.
237
238 if (domains().size()!=0) {
239 set_outermost_domain(outermost_domain());
240 check_geometry_is_nested();
241 }
242
243 generate_indices(OLD_ORDERING);
244 make_mesh_pairs();
245 #ifdef DEBUG
246 for (const auto& mesh : meshes())
247 mesh.check_consistency("geometry finalize step");
248 #endif
249 }
250
252
253 size_t nb_current_barrier_triangles() const { return nb_current_barrier_triangles_; }
254 size_t& nb_current_barrier_triangles() { return nb_current_barrier_triangles_; }
255 size_t nb_invalid_vertices() { return invalid_vertices_.size(); }
256
257 const MeshParts& isolated_parts() const { return independant_parts; }
259 const Mesh& mesh(const std::string& id) const; // Is this useful ?? TODO.
260
261 private:
262
263 void clear() {
264 geom_vertices.clear();
265 geom_meshes.clear();
266 geom_domains.clear();
267 nested = false;
268 outer_domain = 0;
269 num_params = 0;
270 }
271
272 void read_geometry_file(const std::string& filename);
273 void read_conductivity_file(const std::string& filename);
274
275 void make_mesh_pairs();
276
278
279 Vertices geom_vertices;
280 Meshes geom_meshes;
281 Domains geom_domains;
282
283 const Domain* outer_domain = 0;
284 bool nested = false;
285 size_t num_params = 0; // total number = nb of vertices + nb of triangles
286
287 void generate_indices(const bool);
288
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));
294 return doms;
295 }
296
297 // Accumulate a function over the domain common to two meshes.
298
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; }
302
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);
306 double result = 0.0;
307 for (const auto& domainptr : doms)
308 result += Function(*domainptr);
309 return result;
310 }
311
313
314 std::set<Vertex> invalid_vertices_;
315 size_t nb_current_barrier_triangles_ = 0;
316
317 MeshParts independant_parts;
318 MeshPairs meshpairs;
319 };
320}
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
Domains & domains()
Return the list of domains.
Definition geometry.h:142
size_t nb_invalid_vertices()
Definition geometry.h:255
Meshes & meshes()
Return the list of meshes involved in the geometry.
Definition geometry.h:131
const Interface & interface(const std::string &name) const
returns the Interface called
void finalize(const bool OLD_ORDERING=false)
Definition geometry.h:228
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:201
std::vector< const Domain * > DomainsReference
Definition geometry.h:50
bool check(const Mesh &m) const
check if m intersect geometry meshes
Vertices & vertices()
Return the list of vertices involved in the geometry.
Definition geometry.h:95
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="")
Definition geometry.h:111
double sigma_inv(const Mesh &m1, const Mesh &m2) const
Definition geometry.h:183
Geometry(const unsigned n)
Definition geometry.h:59
unsigned add_vertex(const Vertex &V)
Add a vertex.
Definition geometry.h:100
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
Definition geometry.h:132
const Vertices & vertices() const
Definition geometry.h:96
Geometry(const std::string &geomFileName, const bool OLD_ORDERING=false)
Definition geometry.h:63
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:188
Geometry(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
Definition geometry.h:67
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)
Definition geometry.h:74
bool is_outermost(const Domain &domain) const
Definition geometry.h:173
Geometry()
Constructors.
Definition geometry.h:57
size_t & nb_current_barrier_triangles()
Definition geometry.h:254
const MeshParts & isolated_parts() const
Definition geometry.h:257
Geometry(const char *geomFileName, const bool OLD_ORDERING=false)
Definition geometry.h:73
const MeshPairs & communicating_mesh_pairs() const
Definition geometry.h:134
void set_outermost_domain(Domain &domain)
Definition geometry.h:167
bool is_nested() const
Definition geometry.h:91
size_t nb_parameters() const
the total number of vertices + triangles
Definition geometry.h:160
size_t nb_current_barrier_triangles() const
Handle multiple isolated domains.
Definition geometry.h:253
std::vector< MeshPair > MeshPairs
Definition geometry.h:48
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.
Definition geometry.h:152
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)
Definition geometry.h:122
void load(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
Definition geometry.h:217
double indicator(const Mesh &m1, const Mesh &m2) const
Definition geometry.h:184
const Mesh & mesh(const std::string &id) const
void load(const std::string &filename, const bool OLD_ORDERING=false)
Definition geometry.h:211
const Domains & domains() const
Definition geometry.h:143
std::vector< std::pair< std::string, std::string > > MeshList
Definition geometry.h:53
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:182
bool has_conductivities() const
Definition geometry.h:78
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:276
std::vector< Domain > Domains
A vector of Domain is called Domains.
Definition domain.h:112
std::map< unsigned, unsigned > IndexMap
Definition triangle.h:148
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