OpenMEEG
mesh.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 <iostream>
11 
12 #include <vector>
13 #include <map>
14 #include <string>
15 #include <memory>
16 
17 #include <om_common.h>
18 #include <triangle.h>
19 #include <om_utils.h>
20 
21 #include <symmatrix.h>
22 #include <block_matrix.h>
23 
24 namespace OpenMEEG {
25 
26  class Geometry;
27  using maths::Range;
28  using maths::Ranges;
29 
30  // Mesh class
31  // \brief Mesh is a collection of triangles associated to a geometry containing the points
32  // on which triangles are based.
33 
34  class OPENMEEG_EXPORT Mesh {
35 
36  static Geometry* create_geometry(Geometry* geom);
37 
38  public:
39 
40  friend class Geometry;
41  friend class MeshIO;
42 
43  typedef std::map<const Vertex*,TrianglesRefs> VertexTriangles;
44 
47 
48  Mesh(Geometry* geometry=nullptr): geom(create_geometry(geometry)) { }
49 
54  // Do we need this ? TODO
55 
56  Mesh(const unsigned nv,const unsigned nt,Geometry* geometry=nullptr);
57 
58  Mesh(const Mesh&) = default;
59  Mesh(Mesh&& m) = default;
60  // MSVC 2019 https://stackoverflow.com/questions/31264984/c-compiler-error-c2280-attempting-to-reference-a-deleted-function-in-visual
61  Mesh& operator=(const Mesh&) = default;
62 
67 
68  Mesh(const std::string& filename,const bool verbose,Geometry* geometry=nullptr): Mesh(geometry) {
69  load(filename,verbose);
70  }
71 
74 
75  Mesh(const std::string& filename,Geometry* geometry=nullptr): Mesh(filename,false,geometry) { }
76 
78 
79  ~Mesh() { clear(); }
80 
81  std::string& name() { return mesh_name; }
82  const std::string& name() const { return mesh_name; }
83 
84  VerticesRefs& vertices() { return mesh_vertices; }
85  const VerticesRefs& vertices() const { return mesh_vertices; }
86 
87  Geometry& geometry() const { return *geom; }
88 
89  Triangles& triangles() { return mesh_triangles; }
90  const Triangles& triangles() const { return mesh_triangles; }
91 
93 
94  bool current_barrier() const { return current_barrier_; }
95  bool& current_barrier() { return current_barrier_; }
96  bool isolated() const { return isolated_; }
97  bool& isolated() { return isolated_; }
98 
100 
102 
103  Triangle& add_triangle(const TriangleIndices inds,const IndexMap& indmap) {
104  const TriangleIndices t = { indmap.at(inds[0]), indmap.at(inds[1]), indmap.at(inds[2])};
105  return add_triangle(t);
106  }
107 
108  void add(const std::vector<TriangleIndices>& trgs) {
109  for (const auto& triangle : trgs)
110  add_triangle(triangle);
111  }
112 
113  void add(const std::vector<TriangleIndices>& trgs,const IndexMap& indmap) {
114  for (const auto& triangle : trgs)
115  add_triangle(triangle,indmap);
116  }
117 
118  bool operator==(const Mesh& m) const { return triangles()==m.triangles(); }
119  bool operator!=(const Mesh& m) const { return triangles()!=m.triangles(); }
120 
124 
125  void info(const bool verbose=false) const;
126  bool has_self_intersection() const;
127  bool intersection(const Mesh&) const;
128  bool has_correct_orientation() const;
130  void update(const bool topology_changed);
131  void merge(const Mesh&,const Mesh&);
132 
135 
137  std::vector<size_t> indices;
138  for (const auto& vertex : vertices())
139  indices.push_back(vertex->index());
140  std::sort(indices.begin(),indices.end());
141  Ranges result;
142  for (auto it=indices.begin(); it!=indices.end();) {
143  auto it1 = it;
144  for (auto it2=it1+1; it2!=indices.end() && *it2==*it1+1; it1=it2++);
145  result.push_back(Range(*it,*it1));
146  it = it1+1;
147  }
148  return result;
149  }
150 
151  // Triangles always have a contiguous range as they are never shared between meshes.
152 
153  Range triangles_range() const { return Range(triangles().front().index(),triangles().back().index()); }
154 
156 
157  TrianglesRefs triangles(const Vertex& V) const { return vertex_triangles.at(&V); }
158 
160 
161  TrianglesRefs adjacent_triangles(const Triangle& triangle) const {
162  std::map<Triangle*,unsigned> mapt;
163  TrianglesRefs result;
164  for (auto& vertex : triangle)
165  for (const auto& t2 : triangles(*vertex))
166  if (++mapt[t2]==2)
167  result.push_back(t2);
168  return result;
169  }
170 
172 
174  for (auto& triangle : triangles())
175  triangle.change_orientation();
176  }
177 
180  double solid_angle(const Vect3& p) const;
181  Normal normal(const Vertex& v) const;
182  void laplacian(SymMatrix &A) const;
183 
184  bool& outermost() { return outermost_; }
185  bool outermost() const { return outermost_; }
186 
191 
192  void smooth(const double& smoothing_intensity, const unsigned& niter);
193 
195 
196  void gradient_norm2(SymMatrix &A) const;
197 
204 
205  void load(const std::string& filename,const bool verbose=true);
206 
209 
210  void save(const std::string& filename) const ;
211 
212  #if !defined(SWIGPYTHON) && !defined(_MSC_VER)
213  private:
214  #endif
215 
216  // This private method must be accessible from swig.
217 
218  void reference_vertices(const IndexMap& indmap);
219 
220  private:
221 
223 
224  typedef std::map<std::pair<const Vertex*,const Vertex*>,int> EdgeMap;
225 
226  void clear();
227 
230 
231  void add_mesh(const Mesh& m);
232 
233  // regarding mesh orientation
234 
235  const EdgeMap compute_edge_map() const;
236  bool triangle_intersection(const Triangle&,const Triangle&) const;
237 
239 
240  Vect3 P1gradient(const Vect3& p0,const Vect3& p1,const Vect3& p2) const { return crossprod(p1,p2)/det(p0,p1,p2); }
241 
243 
244  double P0gradient_norm2(const Triangle& t1,const Triangle& t2) const {
245  return sqr(dotprod(t1.normal(),t2.normal()))/(t1.center()-t2.center()).norm2();
246  }
247 
248  // Create the map that for each vertex gives the triangles containing it.
249 
250  void make_adjacencies() {
251  vertex_triangles.clear();
252  for (auto& triangle : triangles())
253  for (const auto& vertex : triangle)
254  vertex_triangles[vertex].push_back(&triangle);
255  }
256 
257  typedef std::shared_ptr<Geometry> Geom;
258 
259  std::string mesh_name = "";
260  VertexTriangles vertex_triangles;
261  Geometry* geom;
262  VerticesRefs mesh_vertices;
263  Triangles mesh_triangles;
264  bool outermost_ = false;
265 
267 
268  bool current_barrier_ = false;
269  bool isolated_ = false;
270  };
271 
272  typedef std::vector<Mesh> Meshes;
273 }
Geometry contains the electrophysiological model Vertices, meshes and domains are stored in this geom...
Definition: geometry.h:31
const Triangles & triangles() const
Definition: mesh.h:90
void laplacian(SymMatrix &A) const
Compute mesh laplacian.
bool operator!=(const Mesh &m) const
Definition: mesh.h:119
const VerticesRefs & vertices() const
Definition: mesh.h:85
TriangleIndices triangle(const Triangle &t) const
bool outermost() const
Returns True if it is an outermost mesh.
Definition: mesh.h:185
bool has_correct_orientation() const
Check local orientation of mesh triangles.
std::string & name()
Definition: mesh.h:81
VerticesRefs & vertices()
Definition: mesh.h:84
void correct_local_orientation()
Correct the local orientation of the mesh triangles.
Mesh(const unsigned nv, const unsigned nt, Geometry *geometry=nullptr)
Constructor from scratch (vertices/triangles t be added)
bool & isolated()
Definition: mesh.h:97
Mesh(Mesh &&m)=default
bool intersection(const Mesh &) const
Check whether the mesh intersects another mesh.
Triangles & triangles()
Definition: mesh.h:89
void smooth(const double &smoothing_intensity, const unsigned &niter)
Smooth Mesh.
Normal normal(const Vertex &v) const
Get normal at vertex.`.
void correct_global_orientation()
Correct the global orientation (if there is one).
void merge(const Mesh &, const Mesh &)
Merge two meshes.
bool isolated() const
Definition: mesh.h:96
void add(const std::vector< TriangleIndices > &trgs, const IndexMap &indmap)
Definition: mesh.h:113
Mesh(const std::string &filename, Geometry *geometry=nullptr)
Definition: mesh.h:75
void add(const std::vector< TriangleIndices > &trgs)
Definition: mesh.h:108
~Mesh()
Destructor.
Definition: mesh.h:79
bool has_self_intersection() const
Check whether the mesh self-intersects.
void generate_indices()
Generate indices (if allocate).
TrianglesRefs triangles(const Vertex &V) const
Get the triangles adjacent to vertex.
Definition: mesh.h:157
const std::string & name() const
Definition: mesh.h:82
Mesh & operator=(const Mesh &)=default
std::map< const Vertex *, TrianglesRefs > VertexTriangles
Definition: mesh.h:43
void info(const bool verbose=false) const
Print info Print to std::cout some info about the mesh.
bool operator==(const Mesh &m) const
Definition: mesh.h:118
bool & outermost()
Definition: mesh.h:184
void save(const std::string &filename) const
Save mesh to file.
Triangle & add_triangle(const TriangleIndices inds)
Add a triangle specified by its indices in the geometry.
Triangle & add_triangle(const TriangleIndices inds, const IndexMap &indmap)
Definition: mesh.h:103
Mesh(const std::string &filename, const bool verbose, Geometry *geometry=nullptr)
Constructors.
Definition: mesh.h:68
void change_orientation()
Change mesh orientation.
Definition: mesh.h:173
void update(const bool topology_changed)
Recompute triangles normals, area, and vertex triangles.
bool & current_barrier()
Definition: mesh.h:95
Mesh(Geometry *geometry=nullptr)
Default constructor or constructor using a provided geometry.
Definition: mesh.h:48
bool current_barrier() const
Definition: mesh.h:94
TrianglesRefs adjacent_triangles(const Triangle &triangle) const
Get the triangles adjacent to.
Definition: mesh.h:161
Geometry & geometry() const
Definition: mesh.h:87
void gradient_norm2(SymMatrix &A) const
Compute the square norm of the surfacic gradient.
Ranges vertices_ranges() const
Get the ranges of the specific mesh in the global matrix.
Definition: mesh.h:136
Mesh(const Mesh &)=default
double solid_angle(const Vect3 &p) const
Given a point p, computes the solid angle of the mesh seen from.
void load(const std::string &filename, const bool verbose=true)
Read mesh from file.
Range triangles_range() const
Definition: mesh.h:153
Triangle Triangle class.
Definition: triangle.h:45
Vect3 center() const
Definition: triangle.h:104
Normal & normal()
Definition: triangle.h:95
Vect3.
Definition: vect3.h:28
Vertex.
Definition: vertex.h:20
std::vector< Triangle > Triangles
Definition: triangle.h:138
Vect3 crossprod(const Vect3 &V1, const Vect3 &V2)
Definition: vect3.h:107
std::vector< Mesh > Meshes
Definition: mesh.h:272
std::vector< Triangle * > TrianglesRefs
Definition: triangle.h:139
double det(const Vect3 &V1, const Vect3 &V2, const Vect3 &V3)
Definition: vect3.h:108
std::vector< Vertex * > VerticesRefs
Definition: vertex.h:38
std::map< unsigned, unsigned > IndexMap
Definition: triangle.h:141
double sqr(const double x)
Definition: vect3.h:24
double dotprod(const Vect3 &V1, const Vect3 &V2)
Definition: vect3.h:106