OpenMEEG
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
operators.h
Go to the documentation of this file.
1 /*
2 Project Name : OpenMEEG
3 
4 © INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
5 GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
6 Emmanuel OLIVI
7 Maureen.Clerc.AT.inria.fr, keriven.AT.certis.enpc.fr,
8 kybic.AT.fel.cvut.cz, papadop.AT.inria.fr)
9 
10 The OpenMEEG software is a C++ package for solving the forward/inverse
11 problems of electroencephalography and magnetoencephalography.
12 
13 This software is governed by the CeCILL-B license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL-B
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's authors, the holders of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL-B license and that you accept its terms.
38 */
39 
43 #pragma once
44 
45 #include <iostream>
46 
47 #include <vector.h>
48 #include <matrix.h>
49 #include <symmatrix.h>
50 #include <sparse_matrix.h>
51 #include <geometry.h>
52 #include <integrator.h>
53 #include <analytics.h>
54 
55 namespace OpenMEEG {
56 
57  // #define ADAPT_LHS
58 
59  // T can be a Matrix or SymMatrix
60  void operatorSinternal(const Mesh& , Matrix& , const Vertices&, const double& );
61  void operatorDinternal(const Mesh& , Matrix& , const Vertices&, const double& );
62  void operatorFerguson(const Vect3& , const Mesh& , Matrix& , const unsigned&, const double&);
63  void operatorDipolePotDer(const Vect3& , const Vect3& , const Mesh& , Vector&, const double&, const unsigned, const bool);
64  void operatorDipolePot (const Vect3& , const Vect3& , const Mesh& , Vector&, const double&, const unsigned, const bool);
65 
66  template <typename T>
67  inline void _operatorD(const Triangle& T1,const Triangle& T2,T& mat,const double& coeff,const unsigned gauss_order) {
68  //this version of _operatorD add in the Matrix the contribution of T2 on T1
69  // for all the P1 functions it gets involved
70  // consider varying order of quadrature with the distance between T1 and T2
71  analyticD3 analyD(T2);
72 
73  #ifdef ADAPT_LHS
75  gauss.setOrder(gauss_order);
76  Vect3 total = gauss.integrate(analyD, T1);
77  #else
78  STATIC_OMP Integrator<Vect3, analyticD3> gauss(gauss_order);
79  Vect3 total = gauss.integrate(analyD, T1);
80  #endif //ADAPT_LHS
81 
82  for (unsigned i = 0; i < 3; ++i)
83  mat(T1.index(), T2(i).index()) += total(i) * coeff;
84  }
85 
86  inline void _operatorDinternal(const Triangle& T2,const Vertex& P,Matrix & mat,const double& coeff) {
87  analyticD3 analyD(T2);
88 
89  Vect3 total = analyD.f(P);
90 
91  for (unsigned i=0;i<3;++i)
92  mat(P.index(), T2(i).index()) += total(i) * coeff;
93  }
94 
95  inline double _operatorS(const Triangle& T1,const Triangle& T2,const unsigned gauss_order) {
96  STATIC_OMP Triangle *oldT = 0;
97  STATIC_OMP analyticS analyS;
98 
99  if ( oldT != &T1 ) { // a few computations are needed only when changing triangle T1
100  oldT = (Triangle*)&T1;
101  analyS.init(T1);
102  }
103  #ifdef ADAPT_LHS
105  gauss.setOrder(gauss_order);
106  return gauss.integrate(analyS, T2);
107  #else
108  STATIC_OMP Integrator<double, analyticS> gauss;
109  gauss.setOrder(gauss_order);
110  return gauss.integrate(analyS, T2);
111  #endif //ADAPT_LHS
112  }
113 
114  inline double _operatorSinternal(const Triangle& T,const Vertex& P) {
115  static analyticS analyS;
116  analyS.init(T);
117  return analyS.f(P);
118  }
119 
120  template <typename T>
121  inline double _operatorN(const Vertex& V1,const Vertex& V2,const Mesh& m1,const Mesh& m2,const T& mat) {
122 
123  const Mesh::VectPTriangle& trgs1 = m1.get_triangles_for_vertex(V1);
124  const Mesh::VectPTriangle& trgs2 = m2.get_triangles_for_vertex(V2);
125 
126  double result = 0.0;
127  for (Mesh::VectPTriangle::const_iterator tit1 = trgs1.begin(); tit1 != trgs1.end(); ++tit1 ) {
128  for (Mesh::VectPTriangle::const_iterator tit2 = trgs2.begin(); tit2 != trgs2.end(); ++tit2 ) {
129 
130  // In the second case, we here divided (precalculated) operatorS by the product of areas.
131 
132  const double Iqr = (m1.current_barrier() || m2.current_barrier()) ?
133  mat((*tit1)->index()-m1.begin()->index(),(*tit2)->index()-m2.begin()->index()) :
134  mat((*tit1)->index(),(*tit2)->index())/((*tit1)->area()*(*tit2)->area());
135 
136  Vect3 CB1 = (*tit1)->next(V1)-(*tit1)->prev(V1);
137  Vect3 CB2 = (*tit2)->next(V2)-(*tit2)->prev(V2);
138 
139  const double value = (CB1*CB2)*Iqr;
140 
141  // if it is the same shared vertex
142 
143  result -= (((&m1!=&m2) && (V1==V2)) ? 0.5 : 0.25)*value;
144  }
145  }
146  return result;
147  }
148 
149  inline double _operatorP1P0(const Triangle& T2,const Vertex& V1) {
150  double result = 0.;
151  if (T2.contains(V1))
152  result = T2.area()/3.0;
153  return result;
154  }
155 
156  template <typename T>
157  void operatorN(const Mesh& m1,const Mesh& m2,T& mat,const double& coeff,const unsigned gauss_order) {
158  // This function has the following arguments:
159  // the 2 interacting meshes
160  // the storage Matrix for the result
161  // the coefficient to be appleid to each matrix element (depending on conductivities, ...)
162  // the gauss order parameter (for adaptive integration)
163 
164  std::cout << "OPERATOR N ... (arg : mesh " << m1.name() << " , mesh " << m2.name() << " )" << std::endl;
165 
166  unsigned i = 0; // for the PROGRESSBAR
167  if (&m1==&m2) {
168  if (m1.current_barrier()) {
169  // we thus precompute operator S divided by the product of triangles area.
170 
171  SymMatrix matS(m1.nb_triangles());
172  for (Mesh::const_iterator tit1=m1.begin();tit1!=m1.end();++tit1) {
173  PROGRESSBAR(i++, m1.nb_triangles());
174  #pragma omp parallel for
175  #ifndef OPENMP_3_0
176  for (int i2=tit1-m1.begin();i2<m1.size();++i2) {
177  const Mesh::const_iterator tit2 = m1.begin()+i2;
178  #else
179  for (Mesh::const_iterator tit2=tit1;tit2<m1.end();++tit2) {
180  #endif
181  matS(tit1->index() - m1.begin()->index(), tit2->index() - m1.begin()->index()) = _operatorS(*tit1, *tit2, gauss_order) / ( tit1->area() * tit2->area());
182  }
183  }
184  i = 0 ;
185  for (Mesh::const_vertex_iterator vit1 = m1.vertex_begin(); vit1 != m1.vertex_end(); ++vit1) {
186  PROGRESSBAR(i++, m1.nb_vertices());
187  #pragma omp parallel for
188  #ifndef OPENMP_3_0
189  for (int i2=vit1-m1.vertex_begin();i2<m1.vertex_size();++i2) {
190  const Mesh::const_vertex_iterator vit2 = m1.vertex_begin()+i2;
191  #else
192  for (Mesh::const_vertex_iterator vit2=vit1;vit2<m1.vertex_end();++vit2) {
193  #endif
194  mat((*vit1)->index(), (*vit2)->index()) += _operatorN(**vit1, **vit2, m1, m1, matS) * coeff;
195  }
196  }
197  } else {
198  for (Mesh::const_vertex_iterator vit1 = m1.vertex_begin(); vit1 != m1.vertex_end(); ++vit1) {
199  PROGRESSBAR(i++, m1.nb_vertices());
200  #pragma omp parallel for
201  #ifndef OPENMP_3_0
202  for (int i2=0;i2<=vit1-m1.vertex_begin();++i2) {
203  const Mesh::const_vertex_iterator vit2 = m1.vertex_begin()+i2;
204  #else
205  for (Mesh::const_vertex_iterator vit2=m1.vertex_begin();vit2<=vit1;++vit2) {
206  #endif
207  mat((*vit1)->index(), (*vit2)->index()) += _operatorN(**vit1, **vit2, m1, m1, mat) * coeff;
208  }
209  }
210  }
211  } else {
212  if ( m1.current_barrier() || m2.current_barrier() ) {
213  // we thus precompute operator S divided by the product of triangles area.
214  Matrix matS(m1.nb_triangles(), m2.nb_triangles());
215  for (Mesh::const_iterator tit1 = m1.begin(); tit1 != m1.end(); ++tit1) {
216  PROGRESSBAR(i++, m1.nb_triangles());
217  #pragma omp parallel for
218  #ifndef OPENMP_3_0
219  for (int i2=0;i2<m2.size();++i2) {
220  const Mesh::const_iterator tit2 = m2.begin()+i2;
221  #else
222  for (Mesh::const_iterator tit2=m2.begin();tit2<m2.end();++tit2) {
223  #endif
224  matS(tit1->index()-m1.begin()->index(),tit2->index()-m2.begin()->index()) = _operatorS(*tit1,*tit2,gauss_order)/(tit1->area()*tit2->area());
225  }
226  }
227  i = 0 ;
228  for (Mesh::const_vertex_iterator vit1 = m1.vertex_begin();vit1!=m1.vertex_end();++vit1) {
229  PROGRESSBAR(i++,m1.nb_vertices());
230  #pragma omp parallel for
231  #ifndef OPENMP_3_0
232  for (int i2=0;i2<m2.vertex_size();++i2) {
233  const Mesh::const_vertex_iterator vit2 = m2.vertex_begin()+i2;
234  #else
235  for (Mesh::const_vertex_iterator vit2=m2.vertex_begin();vit2<m2.vertex_end();++vit2) {
236  #endif
237  mat((*vit1)->index(),(*vit2)->index()) += _operatorN(**vit1,**vit2,m1,m2,matS)*coeff;
238  }
239  }
240  } else {
241  // This loop is exactly the same as the one just above with just matS -> mat (factorize).
242  for (Mesh::const_vertex_iterator vit1 = m1.vertex_begin();vit1!=m1.vertex_end();++vit1) {
243  PROGRESSBAR(i++,m1.nb_vertices());
244  #pragma omp parallel for
245  #ifndef OPENMP_3_0
246  for (int i2=0;i2<m2.vertex_size();++i2) {
247  const Mesh::const_vertex_iterator vit2 = m2.vertex_begin()+i2;
248  #else
249  for (Mesh::const_vertex_iterator vit2=m2.vertex_begin();vit2<m2.vertex_end();++vit2) {
250  #endif
251  mat((*vit1)->index(),(*vit2)->index()) += _operatorN(**vit1,**vit2,m1,m2,mat)*coeff;
252  }
253  }
254  }
255  }
256  }
257 
258  template <typename T>
259  void operatorS(const Mesh& m1, const Mesh& m2, T& mat, const double& coeff, const unsigned gauss_order)
260  {
261  // This function has the following arguments:
262  // the 2 interacting meshes
263  // the storage Matrix for the result
264  // the coefficient to be appleid to each matrix element (depending on conductivities, ...)
265  // the gauss order parameter (for adaptive integration)
266 
267  std::cout << "OPERATOR S ... (arg : mesh " << m1.name() << " , mesh " << m2.name() << " )" << std::endl;
268 
269  unsigned i = 0; // for the PROGRESSBAR
270  // The operator S is given by Sij=\Int G*PSI(I, i)*Psi(J, j) with
271  // PSI(A, a) is a P0 test function on layer A and triangle a
272  if (&m1==&m2) {
273  for (Mesh::const_iterator tit1=m1.begin();tit1!=m1.end();++tit1) {
274  PROGRESSBAR(i++,m1.nb_triangles());
275  #pragma omp parallel for
276  #ifndef OPENMP_3_0
277  for (int i2=tit1-m1.begin();i2<m1.size();++i2) {
278  const Mesh::const_iterator tit2 = m1.begin()+i2;
279  #else
280  for (Mesh::const_iterator tit2=tit1;tit2<m1.end();++tit2) {
281  #endif
282  mat(tit1->index(),tit2->index()) = _operatorS(*tit1,*tit2,gauss_order)*coeff;
283  }
284  }
285  } else {
286  // TODO check the symmetry of _operatorS.
287  // if we invert tit1 with tit2: results in HeadMat differs at 4.e-5 which is too big.
288  // using ADAPT_LHS with tolerance at 0.000005 (for _opS) drops this at 6.e-6. (but increase the computation time)
289  for (Mesh::const_iterator tit1=m1.begin();tit1!=m1.end();++tit1) {
290  PROGRESSBAR(i++,m1.nb_triangles());
291  #pragma omp parallel for
292  #ifndef OPENMP_3_0
293  for (int i2=0;i2<m2.size();++i2) {
294  const Mesh::const_iterator tit2 = m2.begin()+i2;
295  #else
296  for (Mesh::const_iterator tit2=m2.begin();tit2<m2.end();++tit2) {
297  #endif
298  mat(tit1->index(),tit2->index()) = _operatorS(*tit1,*tit2,gauss_order)*coeff;
299  }
300  }
301  }
302  }
303 
304  template <typename T>
305  void operatorD(const Mesh& m1, const Mesh& m2, T& mat, const double& coeff, const unsigned gauss_order) {
306  // This function (OPTIMIZED VERSION) has the following arguments:
307  // the 2 interacting meshes
308  // the storage Matrix for the result
309  // the coefficient to be appleid to each matrix element (depending on conductivities, ...)
310  // the gauss order parameter (for adaptive integration)
311 
312  //In this version of the function, in order to skip multiple computations of the same quantities
313  // loops are run over the triangles but the Matrix cannot be filled in this function anymore
314  // That's why the filling is done is function _operatorD
315  //
316 
317  unsigned i = 0; // for the PROGRESSBAR
318  #pragma omp parallel for
319  #ifndef OPENMP_3_0
320  for (int i1=0; i1 < m1.size(); ++i1) {
321  const Mesh::const_iterator tit1 = m1.begin()+i1;
322  #else
323  for (Mesh::const_iterator tit1 = m1.begin(); tit1 < m1.end(); ++tit1) {
324  #endif
325  PROGRESSBAR(i++, m1.nb_triangles());
326  for (Mesh::const_iterator tit2=m2.begin();tit2!=m2.end();++tit2)
327  _operatorD(*tit1,*tit2,mat,coeff,gauss_order);
328  }
329  }
330 
331  template <typename T>
332  void operatorD(const Mesh& m1,const Mesh& m2, T& mat,const double& coeff,const unsigned gauss_order,const bool star) {
333  // This function (OPTIMIZED VERSION) has the following arguments:
334  // the 2 interacting meshes
335  // the storage Matrix for the result
336  // the coefficient to be appleid to each matrix element (depending on conductivities, ...)
337  // the gauss order parameter (for adaptive integration)
338  // an optional star parameter, which denotes the adjoint of the operator
339 
340  std::cout << "OPERATOR D" << ((star) ? "*" : " ") << "... (arg : mesh " << m1.name() << " , mesh " << m2.name() << " )" << std::endl;
341 
342  if (star) {
343  operatorD(m2, m1, mat, coeff, gauss_order);
344  } else {
345  operatorD(m1, m2, mat, coeff, gauss_order);
346  }
347  }
348 
349  template <typename T>
350  void operatorP1P0(const Mesh& m, T& mat,const double& coeff) {
351  // This time mat(i, j)+= ... the Matrix is incremented by the P1P0 operator
352  std::cout << "OPERATOR P1P0... (arg : mesh " << m.name() << " )" << std::endl;
353  for (Mesh::const_iterator tit = m.begin(); tit != m.end(); ++tit)
354  for (Triangle::const_iterator pit = tit->begin(); pit != tit->end(); ++pit)
355  mat(tit->index(), (*pit)->index()) += _operatorP1P0(*tit, **pit) * coeff;
356  }
357 
358  inline Vect3 _operatorFerguson(const Vect3& x,const Vertex& V1,const Mesh& m) {
359  STATIC_OMP Vect3 result;
360  STATIC_OMP analyticS analyS;
361 
362  result.x() = 0.0;
363  result.y() = 0.0;
364  result.z() = 0.0;
365 
366  //loop over triangles of which V1 is a vertex
367  const Mesh::VectPTriangle& trgs = m.get_triangles_for_vertex(V1);
368 
369  for (Mesh::VectPTriangle::const_iterator tit=trgs.begin();tit!=trgs.end();++tit) {
370 
371  const Triangle& T1 = **tit;
372 
373  // A1 , B1 are the two opposite vertices to V1 (triangle A1, B1, V1)
374  Vect3 A1 = T1.next(V1);
375  Vect3 B1 = T1.prev(V1);
376  Vect3 A1B1 = (A1 - B1) * (0.5 / T1.area());
377 
378  analyS.init(V1, A1, B1);
379  const double opS = analyS.f(x);
380 
381  result += (A1B1 * opS);
382  }
383  return result;
384  }
385 }
std::vector< Vertex > Vertices
Definition: vertex.h:73
double f(const Vect3 &x) const
Definition: analytics.h:108
double _operatorP1P0(const Triangle &T2, const Vertex &V1)
Definition: operators.h:149
size_t nb_vertices() const
Definition: mesh.h:148
bool contains(const Vertex &p) const
Definition: triangle.h:157
Vect3 f(const Vect3 &x) const
Definition: analytics.h:137
Vect3 _operatorFerguson(const Vect3 &x, const Vertex &V1, const Mesh &m)
Definition: operators.h:358
double _operatorSinternal(const Triangle &T, const Vertex &P)
Definition: operators.h:114
vertex_iterator vertex_begin()
Definition: mesh.h:133
VectPVertex::const_iterator const_vertex_iterator
Definition: mesh.h:91
const Vertex & prev(const Vertex &V) const
Definition: triangle.h:128
const bool & current_barrier() const
Definition: mesh.h:301
double _operatorS(const Triangle &T1, const Triangle &T2, const unsigned gauss_order)
Definition: operators.h:95
void operatorD(const Mesh &m1, const Mesh &m2, T &mat, const double &coeff, const unsigned gauss_order)
Definition: operators.h:305
void _operatorD(const Triangle &T1, const Triangle &T2, T &mat, const double &coeff, const unsigned gauss_order)
Definition: operators.h:67
Mesh class.
Definition: mesh.h:85
std::string & name()
Definition: mesh.h:144
void operatorP1P0(const Mesh &m, T &mat, const double &coeff)
Definition: operators.h:350
vertex_iterator vertex_end()
Definition: mesh.h:134
void setOrder(const unsigned n)
Definition: integrator.h:191
size_t vertex_size() const
Definition: mesh.h:136
size_t nb_triangles() const
Definition: mesh.h:149
Triangle.
Definition: triangle.h:54
Vertex.
Definition: vertex.h:52
void operatorN(const Mesh &m1, const Mesh &m2, T &mat, const double &coeff, const unsigned gauss_order)
Definition: operators.h:157
virtual T integrate(const I &fc, const Triangle &Trg)
Definition: integrator.h:240
void _operatorDinternal(const Triangle &T2, const Vertex &P, Matrix &mat, const double &coeff)
Definition: operators.h:86
double _operatorN(const Vertex &V1, const Vertex &V2, const Mesh &m1, const Mesh &m2, const T &mat)
Definition: operators.h:121
void init(const Triangle &T)
Definition: analytics.h:80
const Vertex & next(const Vertex &V) const
Definition: triangle.h:140
Vect3.
Definition: vect3.h:62
unsigned & index()
Definition: vertex.h:65
void operatorS(const Mesh &m1, const Mesh &m2, T &mat, const double &coeff, const unsigned gauss_order)
Definition: operators.h:259
void operatorSinternal(const Mesh &, Matrix &, const Vertices &, const double &)
void operatorDipolePot(const Vect3 &, const Vect3 &, const Mesh &, Vector &, const double &, const unsigned, const bool)
unsigned & index()
Definition: triangle.h:123
void operatorFerguson(const Vect3 &, const Mesh &, Matrix &, const unsigned &, const double &)
const VectPTriangle & get_triangles_for_vertex(const Vertex &V) const
get the triangles associated with vertex V
#define PROGRESSBAR(a, b)
Definition: om_utils.h:61
void operatorDinternal(const Mesh &, Matrix &, const Vertices &, const double &)
Matrix class.
Definition: matrix.h:61
std::vector< Triangle * > VectPTriangle
Definition: mesh.h:88
double & area()
Definition: triangle.h:120
void operatorDipolePotDer(const Vect3 &, const Vect3 &, const Mesh &, Vector &, const double &, const unsigned, const bool)