OpenMEEG
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 
14 namespace 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 
117  analyticD3(const Triangle& T):
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
Vertex & vertex(const unsigned &vindex)
Definition: triangle.h:83
Normal & normal()
Definition: triangle.h:95
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