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 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
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
Domains & domains()
Return the list of domains.
Definition: geometry.h:137
std::vector< const Domain * > DomainsReference
Definition: geometry.h:50
size_t nb_invalid_vertices()
Definition: geometry.h:246
Meshes & meshes()
Return the list of meshes involved in the geometry.
Definition: geometry.h:126
const Interface & interface(const std::string &name) const
returns the Interface called
void finalize(const bool OLD_ORDERING=false)
Definition: geometry.h:223
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
Vertices & vertices()
Return the list of vertices involved in the geometry.
Definition: geometry.h:90
const Domain & domain(const std::string &name) const
Get specific domains.
void check_geometry_is_nested()
std::vector< MeshPair > MeshPairs
Definition: geometry.h:48
Domain & outermost_domain()
Returns the outermost domain.
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
unsigned add_vertex(const Vertex &V)
Add a vertex.
Definition: geometry.h:95
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:127
std::vector< std::pair< std::string, std::string > > MeshList
Definition: geometry.h:53
const Vertices & vertices() const
Definition: geometry.h:91
Geometry(const std::string &geomFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:59
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
Geometry(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:63
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:70
bool is_outermost(const Domain &domain) const
Definition: geometry.h:168
Geometry()
Constructors.
Definition: geometry.h:57
size_t & nb_current_barrier_triangles()
Definition: geometry.h:245
const MeshParts & isolated_parts() const
Definition: geometry.h:248
Geometry(const char *geomFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:69
const MeshPairs & communicating_mesh_pairs() const
Definition: geometry.h:129
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
size_t nb_current_barrier_triangles() const
Handle multiple isolated domains.
Definition: geometry.h:244
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:147
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:117
void load(const std::string &geomFileName, const std::string &condFileName, const bool OLD_ORDERING=false)
Definition: geometry.h:212
double indicator(const Mesh &m1, const Mesh &m2) const
Definition: geometry.h:179
const Mesh & mesh(const std::string &id) const
void load(const std::string &filename, const bool OLD_ORDERING=false)
Definition: geometry.h:206
const Domains & domains() const
Definition: geometry.h:138
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
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