OpenMEEG
Loading...
Searching...
No Matches
analytics.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 <triangle.h>
12#include <dipole.h>
13
14namespace OpenMEEG {
15
16 inline double integral_simplified_green(const Vect3& p0x, const double norm2p0x,
17 const Vect3& p1x, const double norm2p1x,
18 const Vect3& p1p0,const double norm2p1p0)
19 {
20 // The quantity arg is normally >= 1, verifying this relates to a triangular inequality
21 // between p0, p1 and x.
22 // Consequently, there is no need of an absolute value in the first case.
23
24 const double arg = (norm2p0x*norm2p1p0-dotprod(p0x,p1p0))/(norm2p1x*norm2p1p0-dotprod(p1x,p1p0));
25 return (std::isnormal(arg) && arg>0.0) ? log(arg) : fabs(log(norm2p1x/norm2p0x));
26 }
27
28 class OPENMEEG_EXPORT analyticS {
29
30 void initialize(const Vect3& v0,const Vect3& v1,const Vect3& v2) {
31 // All computations needed when the first triangle of integration is changed
32
33 p0 = v0;
34 p1 = v1;
35 p2 = v2;
36
37 p1p0 = p1-p0;
38 p2p1 = p2-p1;
39 p0p2 = p0-p2;
40
41 norm2p1p0 = p1p0.norm();
42 norm2p2p1 = p2p1.norm();
43 norm2p0p2 = p0p2.norm();
44 }
45
46 void finish_intialization() {
47 nu0 = (p1p0^n);
48 nu1 = (p2p1^n);
49 nu2 = (p0p2^n);
50 nu0.normalize();
51 nu1.normalize();
52 nu2.normalize();
53 }
54
55 public:
56
57 analyticS(const Triangle& T) {
58 initialize(T.vertex(0),T.vertex(1),T.vertex(2));
59 n = T.normal();
60 finish_intialization();
61 }
62
63 analyticS(const Vect3& v0,const Vect3& v1,const Vect3& v2) {
64 initialize(v0,v1,v2);
65 n = p1p0^p0p2;
66 n /= n.norm();
67 finish_intialization();
68 }
69
70 double f(const Vect3& x) const {
71 // analytical value of the internal integral of S operator at point X
72 const Vect3& p0x = p0-x;
73 const Vect3& p1x = p1-x;
74 const Vect3& p2x = p2-x;
75 const double norm2p0x = p0x.norm();
76 const double norm2p1x = p1x.norm();
77 const double norm2p2x = p2x.norm();
78
79 const double g0 = integral_simplified_green(p0x,norm2p0x,p1x,norm2p1x,p1p0,norm2p1p0);
80 const double g1 = integral_simplified_green(p1x,norm2p1x,p2x,norm2p2x,p2p1,norm2p2p1);
81 const double g2 = integral_simplified_green(p2x,norm2p2x,p0x,norm2p0x,p0p2,norm2p0p2);
82
83 const double alpha = dotprod(p0x,n);
84
85 return ((dotprod(p0x,nu0)*g0+dotprod(p1x,nu1)*g1+dotprod(p2x,nu2)*g2)-alpha*x.solid_angle(p0,p1,p2));
86 }
87
88 private:
89
90 Vect3 p0, p1, p2;
91 Vect3 p2p1, p1p0, p0p2;
92 Vect3 nu0, nu1, nu2;
93 Vect3 n;
94 double norm2p2p1, norm2p1p0, norm2p0p2;
95 };
96
97 class OPENMEEG_EXPORT analyticD3 {
98
99 static Vect3 unit_vector(const Vect3& V) { return V/V.norm(); }
100
101 Vect3 diff(const unsigned i,const unsigned j) const { return triangle.vertex(i)-triangle.vertex(j); }
102
103 // TODO: Introduce a D matrix, and a dipole....
104
105 #if 0
106 Matrix initD() const {
107 Matrix res(3,3);
108 res.setlin(1,D2);
109 res.setlin(2,D3);
110 res.setlin(3,D1);
111 return res;
112 }
113 #endif
114
115 public:
116
118 triangle(T),D1(diff(1,0)),D2(diff(2,1)),D3(diff(0,2)),U1(unit_vector(D1)),U2(unit_vector(D2)),U3(unit_vector(D3))
119 { }
120
121 inline Vect3 f(const Vect3& x) const {
122 // Analytical value of the inner integral in operator D. See DeMunck article for further details.
123 // Used in non-optimized version of operator D.
124 // Returns a vector of the inner integrals of operator D on a triangle wrt its three P1 functions.
125
126 // First part omega is just x.solid_angle(triangle.vertex(0),triangle.vertex(1),triangle.vertex(2))
127
128 const Vect3& Y1 = triangle.vertex(0)-x;
129 const Vect3& Y2 = triangle.vertex(1)-x;
130 const Vect3& Y3 = triangle.vertex(2)-x;
131 const double y1 = Y1.norm();
132 const double y2 = Y2.norm();
133 const double y3 = Y3.norm();
134 const double d = det(Y1,Y2,Y3);
135
136 if (fabs(d)<1e-10)
137 return 0.0;
138
139 const double omega = 2*atan2(d,(y1*y2*y3+y1*dotprod(Y2,Y3)+y2*dotprod(Y3,Y1)+y3*dotprod(Y1,Y2)));
140
141 const Vect3& Z1 = crossprod(Y2,Y3);
142 const Vect3& Z2 = crossprod(Y3,Y1);
143 const Vect3& Z3 = crossprod(Y1,Y2);
144 const double g1 = log((y2+dotprod(Y2,U1))/(y1+dotprod(Y1,U1)));
145 const double g2 = log((y3+dotprod(Y3,U2))/(y2+dotprod(Y2,U2)));
146 const double g3 = log((y1+dotprod(Y1,U3))/(y3+dotprod(Y3,U3)));
147 const Vect3& N = Z1+Z2+Z3;
148 const Vect3& S = U1*g1+U2*g2+U3*g3;
149
150 return (omega*Vect3(dotprod(Z1,N),dotprod(Z2,N),dotprod(Z3,N))+d*Vect3(dotprod(D2,S),dotprod(D3,S),dotprod(D1,S)))/N.norm2();
151 }
152
153 private:
154
155 const Triangle& triangle;
156 //const Matrix D;
157 const Vect3 D1;
158 const Vect3 D2;
159 const Vect3 D3;
160 const Vect3 U1;
161 const Vect3 U2;
162 const Vect3 U3;
163 };
164
165 class OPENMEEG_EXPORT analyticDipPotDer {
166 public:
167
168 analyticDipPotDer(const Dipole& dip,const Triangle& T): dipole(dip) {
169
170 const Vect3& p0 = T.vertex(0);
171 const Vect3& p1 = T.vertex(1);
172 const Vect3& p2 = T.vertex(2);
173
174 const Vect3& p1p0 = p0-p1;
175 const Vect3& p2p1 = p1-p2;
176 const Vect3& p0p2 = p2-p0;
177 const Vect3& p1p0n = p1p0/p1p0.norm();
178 const Vect3& p2p1n = p2p1/p2p1.norm();
179 const Vect3& p0p2n = p0p2/p0p2.norm();
180
181 const Vect3& p1H0 = dotprod(p1p0,p2p1n)*p2p1n;
182 H0 = p1H0+p1;
183 H0p0DivNorm2 = p0-H0;
184 H0p0DivNorm2 = H0p0DivNorm2/H0p0DivNorm2.norm2();
185 const Vect3& p2H1 = dotprod(p2p1,p0p2n)*p0p2n;
186 H1 = p2H1+p2;
187 H1p1DivNorm2 = p1-H1;
188 H1p1DivNorm2 = H1p1DivNorm2/H1p1DivNorm2.norm2();
189 const Vect3& p0H2 = dotprod(p0p2,p1p0n)*p1p0n;
190 H2 = p0H2+p0;
191 H2p2DivNorm2 = p2-H2;
192 H2p2DivNorm2 = H2p2DivNorm2/H2p2DivNorm2.norm2();
193
194 n = -crossprod(p1p0,p0p2);
195 n.normalize();
196 }
197
198 Vect3 f(const Vect3& r) const {
199 Vect3 P1part(dotprod(H0p0DivNorm2,r-H0),dotprod(H1p1DivNorm2,r-H1),dotprod(H2p2DivNorm2,r-H2));
200
201 // B = n.grad_x(A) with grad_x(A)= q/||^3 - 3r(q.r)/||^5
202
203 const Vect3& x = r-dipole.position();
204 const double inv_xnrm2 = 1.0/x.norm2();
205 const double EMpart = dotprod(n,dipole.moment()-3*dotprod(dipole.moment(),x)*x*inv_xnrm2)*(inv_xnrm2*sqrt(inv_xnrm2));
206
207 return -EMpart*P1part; // RK: why - sign ?
208 }
209
210 private:
211
212 const Dipole& dipole;
213
214 Vect3 H0, H1, H2;
215 Vect3 H0p0DivNorm2, H1p1DivNorm2, H2p2DivNorm2, n;
216 };
217}
Matrix class Matrix class.
Definition: matrix.h:28
void setlin(const Index i, const Vector &v)
Definition: matrix.h:261
Triangle Triangle class.
Definition: triangle.h:45
Normal & normal()
Definition: triangle.h:95
Vertex & vertex(const unsigned &vindex)
Definition: triangle.h:83
Vect3.
Definition: vect3.h:28
double norm() const
Definition: vect3.h:63
double solid_angle(const Vect3 &v1, const Vect3 &v2, const Vect3 &v3) const
Definition: vect3.h:110
double norm2() const
Definition: vect3.h:64
analyticD3(const Triangle &T)
Definition: analytics.h:117
Vect3 f(const Vect3 &x) const
Definition: analytics.h:121
Vect3 f(const Vect3 &r) const
Definition: analytics.h:198
analyticDipPotDer(const Dipole &dip, const Triangle &T)
Definition: analytics.h:168
analyticS(const Triangle &T)
Definition: analytics.h:57
analyticS(const Vect3 &v0, const Vect3 &v1, const Vect3 &v2)
Definition: analytics.h:63
double f(const Vect3 &x) const
Definition: analytics.h:70
double integral_simplified_green(const Vect3 &p0x, const double norm2p0x, const Vect3 &p1x, const double norm2p1x, const Vect3 &p1p0, const double norm2p1p0)
Definition: analytics.h:16
Vect3 crossprod(const Vect3 &V1, const Vect3 &V2)
Definition: vect3.h:107
double det(const Vect3 &V1, const Vect3 &V2, const Vect3 &V3)
Definition: vect3.h:108
double dotprod(const Vect3 &V1, const Vect3 &V2)
Definition: vect3.h:106