OpenMEEG
integrator.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 
40 #pragma once
41 
42 #include <cmath>
43 #include <iostream>
44 
45 #include <vertex.h>
46 #include <triangle.h>
47 #include <mesh.h>
48 
49 namespace OpenMEEG {
50 
51  // light class containing d Vect3
52  template <int d>
53  class OPENMEEG_EXPORT Vect3array
54  {
55  Vect3 t[d];
56 
57  public:
58  Vect3array() {};
59  inline Vect3array(double x) {
60  for ( unsigned i = 0; i < d; ++i)
61  t[i] = Vect3(x);
62  }
63  inline Vect3array<d> operator*(double x) const {
64  Vect3array<d> r;
65  for ( unsigned i = 0; i < d; ++i)
66  r.t[i] = t[i]*x;
67  return r;
68  }
69  inline Vect3 operator()(int i) const { return t[i]; }
70  inline Vect3& operator()(int i) { return t[i]; }
71  };
72 
73  template <int d>
74  inline void multadd (Vect3array<d>& target, const double scale, const Vect3array<d>& incr)
75  {
76  for ( unsigned i = 0; i < d; ++i) {
77  target(i) = target(i) + scale*incr(i);
78  }
79  }
80 
81  inline void multadd (double& target, const double scale, const double incr)
82  {
83  target += scale*incr;
84  }
85 
86  inline void multadd (Vect3& target, const double scale, const Vect3& incr)
87  {
88  target = target + scale*incr;
89  }
90 
91  // Quadrature rules are from Marc Bonnet's book: Equations integrales..., Appendix B.3
92 
93  static const double cordBars[4][16][4] =
94  {
95  //parameters for N=3
96  {
97  {0.166666666666667, 0.166666666666667, 0.666666666666667, 0.166666666666667},
98  {0.166666666666667, 0.666666666666667, 0.166666666666667, 0.166666666666667},
99  {0.666666666666667, 0.166666666666667, 0.166666666666667, 0.166666666666667},
100  {0.0, 0.0, 0.0, 0.0},
101  {0.0, 0.0, 0.0, 0.0},
102  {0.0, 0.0, 0.0, 0.0},
103  {0.0, 0.0, 0.0, 0.0},
104  {0.0, 0.0, 0.0, 0.0},
105  {0.0, 0.0, 0.0, 0.0},
106  {0.0, 0.0, 0.0, 0.0},
107  {0.0, 0.0, 0.0, 0.0},
108  {0.0, 0.0, 0.0, 0.0},
109  {0.0, 0.0, 0.0, 0.0},
110  {0.0, 0.0, 0.0, 0.0},
111  {0.0, 0.0, 0.0, 0.0},
112  {0.0, 0.0, 0.0, 0.0}
113  }
114  ,
115  // parameters for N=6
116  {
117  {0.445948490915965, 0.445948490915965, 0.108103018168070, 0.111690794839005},
118  {0.445948490915965, 0.108103018168070, 0.445948490915965, 0.111690794839005},
119  {0.108103018168070, 0.445948490915965, 0.445948490915965, 0.111690794839005},
120  {0.091576213509771, 0.091576213509771, 0.816847572980458, 0.054975871827661},
121  {0.091576213509771, 0.816847572980458, 0.091576213509771, 0.054975871827661},
122  {0.816847572980458, 0.091576213509771, 0.091576213509771, 0.054975871827661},
123  {0.0, 0.0, 0.0, 0.0},
124  {0.0, 0.0, 0.0, 0.0},
125  {0.0, 0.0, 0.0, 0.0},
126  {0.0, 0.0, 0.0, 0.0},
127  {0.0, 0.0, 0.0, 0.0},
128  {0.0, 0.0, 0.0, 0.0},
129  {0.0, 0.0, 0.0, 0.0},
130  {0.0, 0.0, 0.0, 0.0},
131  {0.0, 0.0, 0.0, 0.0},
132  {0.0, 0.0, 0.0, 0.0}
133  }
134  ,
135  // parameters for N=7
136  {
137  {0.333333333333333, 0.333333333333333, 0.333333333333333, 0.1125},
138  {0.470142064105115, 0.470142064105115, 0.059715871789770, 0.066197076394253},
139  {0.470142064105115, 0.059715871789770, 0.470142064105115, 0.066197076394253},
140  {0.059715871789770, 0.470142064105115, 0.470142064105115, 0.066197076394253},
141  {0.101286507323456, 0.101286507323456, 0.797426985353088, 0.062969590272414},
142  {0.101286507323456, 0.797426985353088, 0.101286507323456, 0.062969590272414},
143  {0.797426985353088, 0.101286507323456, 0.101286507323456, 0.062969590272414},
144  {0.0, 0.0, 0.0, 0.0},
145  {0.0, 0.0, 0.0, 0.0},
146  {0.0, 0.0, 0.0, 0.0},
147  {0.0, 0.0, 0.0, 0.0},
148  {0.0, 0.0, 0.0, 0.0},
149  {0.0, 0.0, 0.0, 0.0},
150  {0.0, 0.0, 0.0, 0.0},
151  {0.0, 0.0, 0.0, 0.0},
152  {0.0, 0.0, 0.0, 0.0}
153  }
154  ,
155 
156  // parameters for N=16
157  {
158  {0.333333333333333, 0.333333333333333, 0.333333333333333, 0.072157803838893},
159  {0.081414823414554, 0.459292588292722, 0.459292588292722, 0.047545817133642},
160  {0.459292588292722, 0.081414823414554, 0.459292588292722, 0.047545817133642},
161  {0.459292588292722, 0.459292588292722, 0.081414823414554, 0.047545817133642},
162  {0.898905543365937, 0.050547228317031, 0.050547228317031, 0.016229248811599},
163  {0.050547228317031, 0.898905543365937, 0.050547228317031, 0.016229248811599},
164  {0.050547228317031, 0.050547228317031, 0.898905543365937, 0.016229248811599},
165  {0.658861384496479, 0.170569307751760, 0.170569307751761, 0.051608685267359},
166  {0.170569307751760, 0.658861384496479, 0.170569307751761, 0.051608685267359},
167  {0.170569307751760, 0.170569307751761, 0.658861384496479, 0.051608685267359},
168  {0.008394777409957, 0.728492392955404, 0.263112829634639, 0.013615157087217},
169  {0.728492392955404, 0.008394777409957, 0.263112829634639, 0.013615157087217},
170  {0.728492392955404, 0.263112829634639, 0.008394777409957, 0.013615157087217},
171  {0.008394777409957, 0.263112829634639, 0.728492392955404, 0.013615157087217},
172  {0.263112829634639, 0.008394777409957, 0.728492392955404, 0.013615157087217},
173  {0.263112829634639, 0.728492392955404, 0.008394777409957, 0.013615157087217}
174  }
175 
176  }; // end of gaussTriangleParams
177 
178  static const unsigned nbPts[4] = {3, 6, 7, 16};
179 
180  template <typename T, typename I>
181  class OPENMEEG_EXPORT Integrator
182  {
183  unsigned order;
184 
185  public:
186 
187  inline Integrator() { setOrder(3); }
188  inline Integrator(unsigned ord) { setOrder(ord); }
189  inline virtual ~Integrator() {}
190 
191  inline void setOrder(const unsigned n)
192  {
193  if ( n < 4 ) {
194  order = n;
195  } else {
196  std::cout << "Unavailable Gauss order: min is 1, max is 3" << n << std::endl;
197  order = (n < 1) ? 1 : 3;
198  }
199  }
200 
201  virtual inline T integrate(const I& fc, const Triangle& Trg)
202  {
203  const Vect3 points[3] = { Trg.s1(), Trg.s2(), Trg.s3() };
204  return triangle_integration(fc, points);
205  }
206 
207  protected:
208 
209  inline T triangle_integration(const I& fc, const Vect3 points[3])
210  {
211  // compute double area of triangle defined by points
212  Vect3 crossprod = (points[1] - points[0])^(points[2] - points[0]);
213  double S = crossprod.norm();
214  T result = 0;
215  for ( unsigned i = 0; i < nbPts[order];++i) {
216  Vect3 v(0.0, 0.0, 0.0);
217  for ( unsigned j = 0; j < 3; ++j) {
218  v.multadd(cordBars[order][i][j], points[j]);
219  }
220  multadd(result, cordBars[order][i][3], fc.f(v));
221  }
222  return result*S;
223  }
224  };
225 
226  template <typename T, typename I>
227  class OPENMEEG_EXPORT AdaptiveIntegrator: public Integrator<T, I>
228  {
230 
231  public:
232 
233  inline AdaptiveIntegrator() : tolerance(0.0001) {}
234  inline AdaptiveIntegrator(double tol) : tolerance(tol) {}
236 
237  inline double norm(const double a) { return fabs(a); }
238  inline double norm(const Vect3& a) { return a.norm(); }
239 
240  virtual inline T integrate(const I& fc, const Triangle& Trg) {
241  const Vect3 points[3] = { Trg.s1(), Trg.s2(), Trg.s3() };
242  T I0 = base::triangle_integration(fc, points);
243  return adaptive_integration(fc, points, I0, 0);
244  }
245 
246  private:
247 
248  double tolerance;
249 
250  inline T adaptive_integration(const I& fc, const Vect3 * points, T I0, unsigned n)
251  {
252  Vect3 newpoint0(0.0, 0.0, 0.0);
253  multadd(newpoint0, 0.5, points[0]);
254  multadd(newpoint0, 0.5, points[1]);
255  Vect3 newpoint1(0.0, 0.0, 0.0);
256  multadd(newpoint1, 0.5, points[1]);
257  multadd(newpoint1, 0.5, points[2]);
258  Vect3 newpoint2(0.0, 0.0, 0.0);
259  multadd(newpoint2, 0.5, points[2]);
260  multadd(newpoint2, 0.5, points[0]);
261  Vect3 points1[3] = {points[0], newpoint0, newpoint2};
262  Vect3 points2[3] = {points[1], newpoint1, newpoint0};
263  Vect3 points3[3] = {points[2], newpoint2, newpoint1};
264  Vect3 points4[3] = {newpoint0, newpoint1, newpoint2};
265  T I1 = base::triangle_integration(fc, points1);
266  T I2 = base::triangle_integration(fc, points2);
267  T I3 = base::triangle_integration(fc, points3);
268  T I4 = base::triangle_integration(fc, points4);
269  T sum = I1+I2+I3+I4;
270  if ( norm(I0-sum) > tolerance*norm(I0) ) {
271  n = n+1;
272  if ( n < 10 ) {
273  I1 = adaptive_integration(fc, points1, I1, n);
274  I2 = adaptive_integration(fc, points2, I2, n);
275  I3 = adaptive_integration(fc, points3, I3, n);
276  I4 = adaptive_integration(fc, points4, I4, n);
277  I0 = I1+I2+I3+I4;
278  }
279  }
280  return I0;
281  }
282  };
283 }
OpenMEEG::Triangle
Triangle.
Definition: triangle.h:54
OpenMEEG::Triangle::s3
const Vertex & s3() const
Definition: triangle.h:111
OpenMEEG::AdaptiveIntegrator::AdaptiveIntegrator
AdaptiveIntegrator()
Definition: integrator.h:233
OpenMEEG::AdaptiveIntegrator::integrate
virtual T integrate(const I &fc, const Triangle &Trg)
Definition: integrator.h:240
OpenMEEG::AdaptiveIntegrator::~AdaptiveIntegrator
~AdaptiveIntegrator()
Definition: integrator.h:235
triangle.h
OpenMEEG::Vect3array::operator()
Vect3 operator()(int i) const
Definition: integrator.h:69
OpenMEEG::Vect3array::t
Vect3 t[d]
Definition: integrator.h:55
OpenMEEG::Integrator::setOrder
void setOrder(const unsigned n)
Definition: integrator.h:191
OpenMEEG::AdaptiveIntegrator::norm
double norm(const Vect3 &a)
Definition: integrator.h:238
OpenMEEG::Vect3array::operator*
Vect3array< d > operator*(double x) const
Definition: integrator.h:63
OpenMEEG::AdaptiveIntegrator::tolerance
double tolerance
Definition: integrator.h:248
OpenMEEG::Integrator::Integrator
Integrator(unsigned ord)
Definition: integrator.h:188
OpenMEEG
Definition: analytics.h:45
OpenMEEG::Vect3
Vect3.
Definition: vect3.h:62
OpenMEEG::AdaptiveIntegrator::base
Integrator< T, I > base
Definition: integrator.h:229
OpenMEEG::Integrator::Integrator
Integrator()
Definition: integrator.h:187
OpenMEEG::Vect3array::operator()
Vect3 & operator()(int i)
Definition: integrator.h:70
OpenMEEG::Integrator::integrate
virtual T integrate(const I &fc, const Triangle &Trg)
Definition: integrator.h:201
OpenMEEG::Integrator::~Integrator
virtual ~Integrator()
Definition: integrator.h:189
mesh.h
OpenMEEG::Vect3array::Vect3array
Vect3array()
Definition: integrator.h:58
OpenMEEG::Triangle::s1
const Vertex & s1() const
Definition: triangle.h:109
OpenMEEG::Integrator
Definition: integrator.h:182
OpenMEEG::nbPts
static const unsigned nbPts[4]
Definition: integrator.h:178
vertex.h
OpenMEEG::AdaptiveIntegrator::norm
double norm(const double a)
Definition: integrator.h:237
OpenMEEG::AdaptiveIntegrator::AdaptiveIntegrator
AdaptiveIntegrator(double tol)
Definition: integrator.h:234
OpenMEEG::cordBars
static const double cordBars[4][16][4]
Definition: integrator.h:93
OpenMEEG::Vect3::norm
double norm() const
Definition: vect3.h:96
OpenMEEG::Vect3array
Definition: integrator.h:54
OpenMEEG::Integrator::triangle_integration
T triangle_integration(const I &fc, const Vect3 points[3])
Definition: integrator.h:209
OpenMEEG::Triangle::s2
const Vertex & s2() const
Definition: triangle.h:110
OpenMEEG::AdaptiveIntegrator
Definition: integrator.h:228
OpenMEEG::Vect3::multadd
void multadd(const double &d, const Vect3 &v)
Definition: vect3.h:107
OpenMEEG::AdaptiveIntegrator::adaptive_integration
T adaptive_integration(const I &fc, const Vect3 *points, T I0, unsigned n)
Definition: integrator.h:250
OpenMEEG::Vect3array::Vect3array
Vect3array(double x)
Definition: integrator.h:59
OpenMEEG::Integrator::order
unsigned order
Definition: integrator.h:183
OpenMEEG::multadd
void multadd(Vect3array< d > &target, const double scale, const Vect3array< d > &incr)
Definition: integrator.h:74