OpenMEEG
integrator.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 <cmath>
11 #include <iostream>
12 
13 #include <vertex.h>
14 #include <triangle.h>
15 #include <mesh.h>
16 
17 namespace OpenMEEG {
18 
19  class OPENMEEG_EXPORT Integrator {
20 
21  typedef Vect3 Point;
22  typedef Point TrianglePoints[3];
23 
24  static unsigned safe_order(const unsigned n) {
25  if (n>0 && n<4)
26  return n;
27  std::cout << "Unavailable Gauss order " << n << ": min is 1, max is 3" << std::endl;
28  return (n<1) ? 1 : 3;
29  }
30 
31  public:
32 
33  Integrator(const unsigned ord): Integrator(ord,0,0.0) { }
34  Integrator(const unsigned ord,const double tol): Integrator(ord,10,tol) { }
35  Integrator(const unsigned ord,const unsigned levels,const double tol=0.0001):
36  order(safe_order(ord)),tolerance(tol),max_depth(levels)
37  { }
38 
39  double norm(const double a) const { return fabs(a); }
40  double norm(const Vect3& a) const { return a.norm(); }
41 
42  // TODO: T can be deduced from Function.
43 
44  template <typename Function>
45  decltype(auto) integrate(const Function& function,const Triangle& triangle) const {
46  const TrianglePoints tripts = { triangle.vertex(0), triangle.vertex(1), triangle.vertex(2) };
47  const auto& coarse = triangle_integration(function,tripts);
48  return (max_depth==0) ? coarse : adaptive_integration(function,tripts,coarse,max_depth);
49  }
50 
51  private:
52 
53  template <typename Function>
54  decltype(auto) triangle_integration(const Function& function,const TrianglePoints& triangle) const {
55  using T = decltype(function(Vect3()));
56  T result = 0.0;
57  for (unsigned i=0;i<nbPts[order];++i) {
58  Vect3 v(0.0,0.0,0.0);
59  for (unsigned j=0; j<3; ++j)
60  v.multadd(rules[order][i].barycentric_coordinates[j],triangle[j]);
61  result += rules[order][i].weight*function(v);
62  }
63 
64  // compute double area of triangle defined by points
65 
66  const double area2 = crossprod(triangle[1]-triangle[0],triangle[2]-triangle[0]).norm();
67  return result*area2;
68  }
69 
70  template <typename T,typename Function>
71  T adaptive_integration(const Function& function,const TrianglePoints& triangle,const T& coarse,const unsigned level) const {
72  const Point midpoints[] = { 0.5*(triangle[1]+triangle[2]), 0.5*(triangle[2]+triangle[0]), 0.5*(triangle[0]+triangle[1]) };
73  const TrianglePoints new_triangles[] = {
74  { triangle[0], midpoints[1], midpoints[2] }, { midpoints[0], triangle[1], midpoints[2] },
75  { midpoints[0], midpoints[1], triangle[2] }, { midpoints[0], midpoints[1], midpoints[2] }
76  };
77 
78  T refined = 0.0;
79  T integrals[4];
80  for (unsigned i=0; i<4; ++i) {
81  integrals[i] = triangle_integration(function,new_triangles[i]);
82  refined += integrals[i];
83  }
84 
85  if (norm(coarse-refined)<=tolerance*norm(coarse) || level==0)
86  return refined;
87 
88  T sum = 0.0;
89  for (unsigned i=0; i<4; ++i)
90  sum += adaptive_integration(function,new_triangles[i],integrals[i],level-1);
91  return sum;
92  }
93 
94  static constexpr unsigned nbPts[4] = { 3, 6, 7, 16 };
95 
96  const unsigned order;
97  const double tolerance;
98  const unsigned max_depth;
99 
100  // Quadrature rules are from Marc Bonnet's book: Equations integrales..., Appendix B.3
101 
102  struct QuadratureRule {
103  double barycentric_coordinates[3];
104  double weight;
105  };
106 
107  static constexpr QuadratureRule rules[4][16] = {
108  { // Parameters for N=3
109  {{ 0.166666666666667, 0.166666666666667, 0.666666666666667 }, 0.166666666666667 },
110  {{ 0.166666666666667, 0.666666666666667, 0.166666666666667 }, 0.166666666666667 },
111  {{ 0.666666666666667, 0.166666666666667, 0.166666666666667 }, 0.166666666666667 },
112  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 },
113  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 },
114  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 },
115  {{ 0.0, 0.0, 0.0 }, 0.0 }
116  },
117  { // Parameters for N=6
118  {{ 0.445948490915965, 0.445948490915965, 0.108103018168070 }, 0.111690794839005 },
119  {{ 0.445948490915965, 0.108103018168070, 0.445948490915965 }, 0.111690794839005 },
120  {{ 0.108103018168070, 0.445948490915965, 0.445948490915965 }, 0.111690794839005 },
121  {{ 0.091576213509771, 0.091576213509771, 0.816847572980458 }, 0.054975871827661 },
122  {{ 0.091576213509771, 0.816847572980458, 0.091576213509771 }, 0.054975871827661 },
123  {{ 0.816847572980458, 0.091576213509771, 0.091576213509771 }, 0.054975871827661 },
124  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 },
125  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 },
126  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }
127  },
128  { // Parameters for N=7
129  {{ 0.333333333333333, 0.333333333333333, 0.333333333333333 }, 0.1125 },
130  {{ 0.470142064105115, 0.470142064105115, 0.059715871789770 }, 0.066197076394253 },
131  {{ 0.470142064105115, 0.059715871789770, 0.470142064105115 }, 0.066197076394253 },
132  {{ 0.059715871789770, 0.470142064105115, 0.470142064105115 }, 0.066197076394253 },
133  {{ 0.101286507323456, 0.101286507323456, 0.797426985353088 }, 0.062969590272414 },
134  {{ 0.101286507323456, 0.797426985353088, 0.101286507323456 }, 0.062969590272414 },
135  {{ 0.797426985353088, 0.101286507323456, 0.101286507323456 }, 0.062969590272414 },
136  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 },
137  {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 },
138  {{ 0.0, 0.0, 0.0 }, 0.0 }
139  },
140  { // Parameters for N=16
141  {{ 0.333333333333333, 0.333333333333333, 0.333333333333333 }, 0.072157803838893 },
142  {{ 0.081414823414554, 0.459292588292722, 0.459292588292722 }, 0.047545817133642 },
143  {{ 0.459292588292722, 0.081414823414554, 0.459292588292722 }, 0.047545817133642 },
144  {{ 0.459292588292722, 0.459292588292722, 0.081414823414554 }, 0.047545817133642 },
145  {{ 0.898905543365937, 0.050547228317031, 0.050547228317031 }, 0.016229248811599 },
146  {{ 0.050547228317031, 0.898905543365937, 0.050547228317031 }, 0.016229248811599 },
147  {{ 0.050547228317031, 0.050547228317031, 0.898905543365937 }, 0.016229248811599 },
148  {{ 0.658861384496479, 0.170569307751760, 0.170569307751761 }, 0.051608685267359 },
149  {{ 0.170569307751760, 0.658861384496479, 0.170569307751761 }, 0.051608685267359 },
150  {{ 0.170569307751760, 0.170569307751761, 0.658861384496479 }, 0.051608685267359 },
151  {{ 0.008394777409957, 0.728492392955404, 0.263112829634639 }, 0.013615157087217 },
152  {{ 0.728492392955404, 0.008394777409957, 0.263112829634639 }, 0.013615157087217 },
153  {{ 0.728492392955404, 0.263112829634639, 0.008394777409957 }, 0.013615157087217 },
154  {{ 0.008394777409957, 0.263112829634639, 0.728492392955404 }, 0.013615157087217 },
155  {{ 0.263112829634639, 0.008394777409957, 0.728492392955404 }, 0.013615157087217 },
156  {{ 0.263112829634639, 0.728492392955404, 0.008394777409957 }, 0.013615157087217 }
157  }
158  };
159  };
160 }
Integrator(const unsigned ord, const unsigned levels, const double tol=0.0001)
Definition: integrator.h:35
Integrator(const unsigned ord, const double tol)
Definition: integrator.h:34
double norm(const double a) const
Definition: integrator.h:39
Integrator(const unsigned ord)
Definition: integrator.h:33
double norm(const Vect3 &a) const
Definition: integrator.h:40
Triangle Triangle class.
Definition: triangle.h:45
Vect3.
Definition: vect3.h:28
double norm() const
Definition: vect3.h:63
Vect3 crossprod(const Vect3 &V1, const Vect3 &V2)
Definition: vect3.h:107