OpenMEEG
Loading...
Searching...
No Matches
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
17namespace 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 #ifndef SWIGPYTHON // SWIG sees the integrate def as a syntax error
45 template <typename Function>
46 decltype(auto) integrate(const Function& function,const Triangle& triangle) const {
47 const TrianglePoints tripts = { triangle.vertex(0), triangle.vertex(1), triangle.vertex(2) };
48 const auto& coarse = triangle_integration(function,tripts);
49 return (max_depth==0) ? coarse : adaptive_integration(function,tripts,coarse,max_depth);
50 }
51 #endif /* not SWIGPYTHON */
52
53 private:
54
55 #ifndef SWIGPYTHON // SWIG sees the triangle_integration def as a syntax error
56 template <typename Function>
57 decltype(auto) triangle_integration(const Function& function,const TrianglePoints& triangle) const {
58 using T = decltype(function(Vect3()));
59 T result = 0.0;
60 for (unsigned i=0;i<nbPts[order];++i) {
61 Vect3 v(0.0,0.0,0.0);
62 for (unsigned j=0; j<3; ++j)
63 v.multadd(rules[order][i].barycentric_coordinates[j],triangle[j]);
64 result += rules[order][i].weight*function(v);
65 }
66
67 // compute double area of triangle defined by points
68
69 const double area2 = crossprod(triangle[1]-triangle[0],triangle[2]-triangle[0]).norm();
70 return result*area2;
71 }
72 #endif /* not SWIGPYTHON */
73
74 template <typename T,typename Function>
75 T adaptive_integration(const Function& function,const TrianglePoints& triangle,const T& coarse,const unsigned level) const {
76 const Point midpoints[] = { 0.5*(triangle[1]+triangle[2]), 0.5*(triangle[2]+triangle[0]), 0.5*(triangle[0]+triangle[1]) };
77 const TrianglePoints new_triangles[] = {
78 { triangle[0], midpoints[1], midpoints[2] }, { midpoints[0], triangle[1], midpoints[2] },
79 { midpoints[0], midpoints[1], triangle[2] }, { midpoints[0], midpoints[1], midpoints[2] }
80 };
81
82 T refined = 0.0;
83 T integrals[4];
84 for (unsigned i=0; i<4; ++i) {
85 integrals[i] = triangle_integration(function,new_triangles[i]);
86 refined += integrals[i];
87 }
88
89 if (norm(coarse-refined)<=tolerance*norm(coarse) || level==0)
90 return refined;
91
92 T sum = 0.0;
93 for (unsigned i=0; i<4; ++i)
94 sum += adaptive_integration(function,new_triangles[i],integrals[i],level-1);
95 return sum;
96 }
97
98 static constexpr unsigned nbPts[4] = { 3, 6, 7, 16 };
99
100 const unsigned order;
101 const double tolerance;
102 const unsigned max_depth;
103
104 // Quadrature rules are from Marc Bonnet's book: Equations integrales..., Appendix B.3
105
106 struct QuadratureRule {
107 double barycentric_coordinates[3];
108 double weight;
109 };
110
111 static constexpr QuadratureRule rules[4][16] = {
112 { // Parameters for N=3
113 {{ 0.166666666666667, 0.166666666666667, 0.666666666666667 }, 0.166666666666667 },
114 {{ 0.166666666666667, 0.666666666666667, 0.166666666666667 }, 0.166666666666667 },
115 {{ 0.666666666666667, 0.166666666666667, 0.166666666666667 }, 0.166666666666667 },
116 {{ 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 },
117 {{ 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 },
118 {{ 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 },
119 {{ 0.0, 0.0, 0.0 }, 0.0 }
120 },
121 { // Parameters for N=6
122 {{ 0.445948490915965, 0.445948490915965, 0.108103018168070 }, 0.111690794839005 },
123 {{ 0.445948490915965, 0.108103018168070, 0.445948490915965 }, 0.111690794839005 },
124 {{ 0.108103018168070, 0.445948490915965, 0.445948490915965 }, 0.111690794839005 },
125 {{ 0.091576213509771, 0.091576213509771, 0.816847572980458 }, 0.054975871827661 },
126 {{ 0.091576213509771, 0.816847572980458, 0.091576213509771 }, 0.054975871827661 },
127 {{ 0.816847572980458, 0.091576213509771, 0.091576213509771 }, 0.054975871827661 },
128 {{ 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 },
129 {{ 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 },
130 {{ 0.0, 0.0, 0.0 }, 0.0 }, {{ 0.0, 0.0, 0.0 }, 0.0 }
131 },
132 { // Parameters for N=7
133 {{ 0.333333333333333, 0.333333333333333, 0.333333333333333 }, 0.1125 },
134 {{ 0.470142064105115, 0.470142064105115, 0.059715871789770 }, 0.066197076394253 },
135 {{ 0.470142064105115, 0.059715871789770, 0.470142064105115 }, 0.066197076394253 },
136 {{ 0.059715871789770, 0.470142064105115, 0.470142064105115 }, 0.066197076394253 },
137 {{ 0.101286507323456, 0.101286507323456, 0.797426985353088 }, 0.062969590272414 },
138 {{ 0.101286507323456, 0.797426985353088, 0.101286507323456 }, 0.062969590272414 },
139 {{ 0.797426985353088, 0.101286507323456, 0.101286507323456 }, 0.062969590272414 },
140 {{ 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 },
141 {{ 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 },
142 {{ 0.0, 0.0, 0.0 }, 0.0 }
143 },
144 { // Parameters for N=16
145 {{ 0.333333333333333, 0.333333333333333, 0.333333333333333 }, 0.072157803838893 },
146 {{ 0.081414823414554, 0.459292588292722, 0.459292588292722 }, 0.047545817133642 },
147 {{ 0.459292588292722, 0.081414823414554, 0.459292588292722 }, 0.047545817133642 },
148 {{ 0.459292588292722, 0.459292588292722, 0.081414823414554 }, 0.047545817133642 },
149 {{ 0.898905543365937, 0.050547228317031, 0.050547228317031 }, 0.016229248811599 },
150 {{ 0.050547228317031, 0.898905543365937, 0.050547228317031 }, 0.016229248811599 },
151 {{ 0.050547228317031, 0.050547228317031, 0.898905543365937 }, 0.016229248811599 },
152 {{ 0.658861384496479, 0.170569307751760, 0.170569307751761 }, 0.051608685267359 },
153 {{ 0.170569307751760, 0.658861384496479, 0.170569307751761 }, 0.051608685267359 },
154 {{ 0.170569307751760, 0.170569307751761, 0.658861384496479 }, 0.051608685267359 },
155 {{ 0.008394777409957, 0.728492392955404, 0.263112829634639 }, 0.013615157087217 },
156 {{ 0.728492392955404, 0.008394777409957, 0.263112829634639 }, 0.013615157087217 },
157 {{ 0.728492392955404, 0.263112829634639, 0.008394777409957 }, 0.013615157087217 },
158 {{ 0.008394777409957, 0.263112829634639, 0.728492392955404 }, 0.013615157087217 },
159 {{ 0.263112829634639, 0.008394777409957, 0.728492392955404 }, 0.013615157087217 },
160 {{ 0.263112829634639, 0.728492392955404, 0.008394777409957 }, 0.013615157087217 }
161 }
162 };
163 };
164}
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
decltype(auto) integrate(const Function &function, const Triangle &triangle) const
Definition integrator.h:46
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
Vertex & vertex(const unsigned &vindex)
Definition triangle.h:83
Vect3.
Definition vect3.h:28
double norm() const
Definition vect3.h:63