OpenMEEG
operators.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 
12 
13 #include <iostream>
14 #include <iomanip>
15 
16 #include <vector.h>
17 #include <matrix.h>
18 #include <symmatrix.h>
19 #include <symm_block_matrix.h>
20 #include <sparse_matrix.h>
21 #include <geometry.h>
22 #include <integrator.h>
23 #include <analytics.h>
24 
25 #include <progressbar.h>
26 
27 namespace OpenMEEG {
28 
29  void operatorFerguson(const Vect3&,const Mesh&,Matrix&,const unsigned&,const double);
30  void operatorDipolePotDer(const Dipole&,const Mesh&,Vector&,const double,const Integrator&);
31  void operatorDipolePot(const Dipole&,const Mesh&,Vector&,const double,const Integrator&);
32 
33  namespace Details {
34 
35  inline Vect3 operatorFerguson(const Vect3& x,const Vertex& V,const Mesh& m) {
36  Vect3 result;
37 
38  // Loop over triangles of which V is a vertex
39 
40  result = 0.0;
41  for (const auto& tp : m.triangles(V)) {
42  const Triangle& T = *tp;
43  const Edge& edge = T.edge(V);
44 
45  // A, B are the two opposite vertices to V (triangle A, B, V)
46  const Vertex& A = edge.vertex(0);
47  const Vertex& B = edge.vertex(1);
48  const Vect3& AB = (A-B)/(2*T.area());
49 
50  const analyticS analyS(V,A,B);
51 
52  result += AB*analyS.f(x);
53  }
54 
55  return result;
56  }
57  }
58 
59  class BlocksBase {
60  public:
61 
62  BlocksBase(const Integrator& intg): integrator(intg) { }
63 
64  void message(const char* op_name,const Mesh& mesh) const {
65  if (verbose)
66  std::cout << "OPERATOR " << std::left << std::setw(2) << op_name << "... (arg : mesh " << mesh.name() << " )" << std::endl;
67  }
68 
69  void message(const char* op_name,const Mesh& mesh1,const Mesh& mesh2) const {
70  if (verbose)
71  std::cout << "OPERATOR " << std::left << std::setw(2) << op_name << "... (arg : mesh " << mesh1.name() << " , mesh " << mesh2.name() << " )" << std::endl;
72  }
73 
74  protected:
75 
76  template <typename T>
77  void D(const Triangles& triangles1,const Triangles& triangles2,const double coeff,T& mat) const {
78  // This function (OPTIMIZED VERSION) has the following arguments:
79  // - 2 interacting meshes (as sets of triangles).
80  // - coefficient to be applied to each matrix element (depending on conductivities, ...)
81  // - storage Matrix for the result.
82 
83  ProgressBar pb(triangles1.size());
84  #pragma omp parallel for
85  #if defined NO_OPENMP || defined OPENMP_RANGEFOR
86  for (const auto& triangle1 : triangles1) {
87  #elif defined OPENMP_ITERATOR
88  for (Triangles::const_iterator tit1=triangles1.begin(); tit1<triangles1.end(); ++tit1) {
89  const Triangle& triangle1 = *tit1;
90  #else
91  for (int i1=0; i1<static_cast<int>(triangles1.size()); ++i1) {
92  const Triangle& triangle1 = *(triangles1.begin()+i1);
93  #endif
94  for (const auto& triangle2 : triangles2) {
95  const analyticD3 analyD(triangle2);
96  const auto& Dfunc = [&analyD](const Vect3& r) { return analyD.f(r); };
97  const Vect3& total = integrator.integrate(Dfunc,triangle1);
98 
99  for (unsigned i=0; i<3; ++i)
100  mat(triangle1.index(),triangle2.vertex(i).index()) += total(i)*coeff;
101  }
102  ++pb;
103  }
104  }
105 
106  // Operator N for two vertices of the same mesh.
107 
108  template <typename T>
109  static double N(const Vertex& V1,const Vertex& V2,const Mesh& m,const T& matrix) {
110  return N(0.25,V1,V2,m,m,matrix);
111  }
112 
113  // This function assumes that the two meshes are different.
114 
115  template <typename T>
116  static double N(const Vertex& V1,const Vertex& V2,const Mesh& m1,const Mesh& m2,const T& matrix) {
117  const double coeff = (V1==V2) ? 0.5 : 0.25;
118  return N(coeff,V1,V2,m1,m2,matrix);
119  }
120 
121  private:
122 
123  template <typename T>
124  static double N(const double factor,const Vertex& V1,const Vertex& V2,const Mesh& m1,const Mesh& m2,const T& matrix) {
125 
126  double result = 0.0;
127  for (const auto& tp1 : m1.triangles(V1)) {
128  const Edge& edge1 = tp1->edge(V1);
129  const Vect3& CB1 = edge1.vertex(0)-edge1.vertex(1);
130  for (const auto& tp2 : m2.triangles(V2)) {
131 
132  const Edge& edge2 = tp2->edge(V2);
133  const Vect3& CB2 = edge2.vertex(0)-edge2.vertex(1);
134 
135  result -= factor*dotprod(CB1,CB2)*matrix(tp1->index(),tp2->index())/(tp1->area()*tp2->area());
136  }
137  }
138  return result;
139  }
140 
141  protected:
142 
144  bool verbose = true;
145  };
146 
147  class DiagonalBlock: public BlocksBase {
148 
149  typedef BlocksBase base;
150 
151  class SymBloc: public SymMatrix {
152 
153  typedef SymMatrix base;
154 
155  public:
156 
157  SymBloc(const unsigned off,const unsigned sz): base(sz),offset(off) { }
158 
159  double& operator()(const unsigned i,const unsigned j) { return base::operator()(i-offset,j-offset); }
160  double operator()(const unsigned i,const unsigned j) const { return base::operator()(i-offset,j-offset); }
161 
162  private:
163 
164  const unsigned offset;
165  };
166 
167  public:
168 
169  DiagonalBlock(const Mesh& m,const Integrator& intg): base(intg),mesh(m) { }
170 
171  template <typename T>
172  void set_S_block(const double coeff,T& matrix) {
173  if (!mesh.current_barrier()) {
174  S(coeff,matrix);
175  Scoeff = coeff;
176  }
177  }
178 
179  template <typename T>
180  void set_N_block(const double coeff,T& matrix) const { N(coeff,matrix); }
181 
182  template <typename T>
183  void set_D_block(const double coeff,T& matrix) const {
184  if (!mesh.current_barrier())
185  D(coeff,matrix);
186  }
187 
188  template <typename T>
189  void set_Dstar_block(const double /* coeff */,T& /* matrix */) const { }
190 
191  template <typename T>
192  void addId(const double coeff,T& matrix) const {
193  // The Matrix is incremented by the identity P1P0 operator
194  base::message("Id",mesh);
195  for (const auto& triangle : mesh.triangles())
196  for (const auto& vertex : triangle)
197  matrix(triangle.index(),vertex->index()) += Id(triangle,*vertex)*coeff;
198  }
199 
200  template <typename T>
201  void S(const double coeff,T& matrix) const {
202  base::message("S",mesh,mesh);
203  ProgressBar pb(mesh.triangles().size());
204 
205  // Operator S is given by Sij=\Int G*PSI(I,i)*Psi(J,j) with PSI(l,t) a P0 test function on layer l and triangle t.
206  // When meshes are equal, optimized computation for a symmetric matrix.
207 
208  const Triangles& triangles = mesh.triangles();
209  for (Triangles::const_iterator tit1=triangles.begin(); tit1!=triangles.end(); ++tit1,++pb) {
210  const Triangle& triangle1 = *tit1;
211  const analyticS analyS(triangle1);
212  const auto& Sfunc = [&analyS](const Vect3& r) { return analyS.f(r); };
213 
214  #pragma omp parallel for
215  #if defined NO_OPENMP || defined OPENMP_ITERATOR
216  for (Triangles::const_iterator tit2=tit1; tit2!=triangles.end(); ++tit2) {
217  const Triangle& triangle2 = *tit2;
218  #else
219  for (int i2=tit1-triangles.begin(); i2<static_cast<int>(triangles.size()); ++i2) {
220  const Triangle& triangle2 = *(triangles.begin()+i2);
221  #endif
222  matrix(triangle1.index(),triangle2.index()) = base::integrator.integrate(Sfunc,triangle2)*coeff;
223  }
224  }
225  }
226 
227  template <typename T>
228  void N(const double coeff,T& matrix) const {
229  if (S_block_is_computed()) {
230  N(coeff/Scoeff,matrix,matrix);
231  } else {
232  SymBloc Sbloc(mesh.triangles().front().index(),mesh.triangles().size());
233  S(1.0,Sbloc);
234  N(coeff,Sbloc,matrix);
235  }
236  }
237 
238  template <typename T>
239  void D(const double coeff,T& matrix) const {
240  base::message("D",mesh,mesh);
241  base::D(mesh.triangles(),mesh.triangles(),coeff,matrix);
242  }
243 
244  template <typename T>
245  void Dstar(const double coeff,T& matrix) const {
246  base::message("D*",mesh,mesh);
247  base::D(mesh.triangles(),mesh.triangles(),coeff,matrix);
248  }
249 
250 
251  private:
252 
253  bool S_block_is_computed() const { return Scoeff!=0.0; }
254 
255  static double Id(const Triangle& T,const Vertex& V) {
256  return (T.contains(V)) ? T.area()/3.0 : 0.0;
257  }
258 
259  template <typename T1,typename T2>
260  void N(const double coeff,const T1& S,T2& matrix) const {
261 
262  base::message("N",mesh,mesh);
263 
264  ProgressBar pb(mesh.vertices().size());
265 
266  // When meshes are equal, optimized computation for a symmetric matrix.
267 
268  for (auto vit1=mesh.vertices().begin(); vit1!=mesh.vertices().end(); ++vit1) {
269  #pragma omp parallel for
270  #if defined NO_OPENMP || defined OPENMP_ITERATOR
271  for (auto vit2=vit1; vit2<mesh.vertices().end(); ++vit2) {
272  #else
273  for (int i2=0;i2<=vit1-mesh.vertices().begin();++i2) {
274  const auto vit2 = mesh.vertices().begin()+i2;
275  #endif
276  matrix((*vit1)->index(),(*vit2)->index()) += base::N(**vit1,**vit2,mesh,S)*coeff;
277  }
278  ++pb;
279  }
280  }
281 
282  const Mesh& mesh;
283  double Scoeff = 0.0;
284  };
285 
286  class PartialBlock {
287  public:
288 
289  PartialBlock(const Mesh& m): mesh(m) { }
290 
291  void addD(const double coeff,const Vertices& points,Matrix& matrix) const {
292  std::cout << "PARTAL OPERATOR D..." << std::endl;
293  for (const auto& triangle : mesh.triangles()) {
294  const analyticD3 analyD(triangle);
295  for (const auto& vertex : points) {
296  const Vect3& integrals = analyD.f(vertex);
297  for (unsigned i=0;i<3;++i)
298  matrix(vertex.index(),triangle.vertex(i).index()) += integrals(i)*coeff;
299  }
300  }
301  }
302 
303  void S(const double coeff,const Vertices& points,Matrix& matrix) const {
304  std::cout << "PARTIAL OPERATOR S..." << std::endl;
305  for (const auto& triangle : mesh.triangles()) {
306  const analyticS analyS(triangle);
307  for (const auto& vertex : points)
308  matrix(vertex.index(),triangle.index()) = coeff*analyS.f(vertex);
309  }
310  }
311 
312 
313  private:
314 
315  const Mesh& mesh;
316  };
317 
318  class NonDiagonalBlock: public BlocksBase {
319 
320  typedef BlocksBase base;
321 
322  class Bloc: public Matrix {
323 
324  typedef Matrix base;
325 
326  public:
327 
328  Bloc(const unsigned r0,const unsigned c0,const unsigned n,const unsigned m): base(n,m),i0(r0),j0(c0) { }
329 
330  double& operator()(const unsigned i,const unsigned j) { return base::operator()(i-i0,j-j0); }
331  double operator()(const unsigned i,const unsigned j) const { return base::operator()(i-i0,j-j0); }
332 
333  private:
334 
335  const unsigned i0;
336  const unsigned j0;
337  };
338 
339  public:
340 
341  // This constructor takes the following arguments:
342  // - The 2 interacting meshes.
343  // - The gauss order parameter (for adaptive integration).
344  // - A verbosity parameters (for printing the action on the terminal).
345 
346  NonDiagonalBlock(const Mesh& m1,const Mesh& m2,const Integrator& intg): base(intg),mesh1(m1),mesh2(m2) { }
347 
348  template <typename T>
349  void set_S_block(const double coeff,T& matrix) {
350  if (!mesh1.current_barrier() && !mesh2.current_barrier()) {
351  S(coeff,matrix);
352  Scoeff = coeff;
353  }
354  }
355 
356  template <typename T>
357  void set_N_block(const double coeff,T& matrix) const { N(coeff,matrix); }
358 
359  template <typename T>
360  void set_D_block(const double coeff,T& matrix) const {
361  if (!mesh1.current_barrier())
362  D(coeff,matrix);
363  }
364 
365  template <typename T>
366  void set_Dstar_block(const double coeff,T& matrix) const {
367  if (mesh1!=mesh2 && !mesh2.current_barrier())
368  Dstar(coeff,matrix);
369  }
370 
371  // Various operators take the following arguments:
372  // - The coefficient to be applied to each matrix element (depending on conductivities, ...)
373  // - The storage Matrix for the result
374 
375  template <typename T>
376  void S(const double coeff,T& matrix) const {
377  base::message("S",mesh1,mesh2);
378  ProgressBar pb(mesh1.triangles().size());
379 
380  // Operator S is given by Sij=\Int G*PSI(I,i)*Psi(J,j) with PSI(l,t) a P0 test function on layer l and triangle t.
381 
382  // TODO check the symmetry of S.
383  // if we invert tit1 with tit2: results in HeadMat differs at 4.e-5 which is too big.
384  // using ADAPT_LHS with tolerance at 0.000005 (for S) drops this at 6.e-6 (but increase the computation time).
385 
386  for (const auto& triangle1 : mesh1.triangles()) {
387 
388  const analyticS analyS(triangle1);
389  const auto& Sfunc = [&analyS](const Vect3& r) { return analyS.f(r); };
390 
391  const Triangles& m2_triangles = mesh2.triangles();
392  #pragma omp parallel for
393  #if defined NO_OPENMP || defined OPENMP_RANGEFOR
394  for (const auto& triangle2 : m2_triangles) {
395  #elif defined OPENMP_ITERATOR
396  for (Triangles::const_iterator tit2=m2_triangles.begin();tit2<m2_triangles.end();++tit2) {
397  const Triangle& triangle2 = *tit2;
398  #else
399  for (int i2=0;i2<static_cast<int>(m2_triangles.size());++i2) {
400  const Triangle& triangle2 = *(m2_triangles.begin()+i2);
401  #endif
402  matrix(triangle1.index(),triangle2.index()) = base::integrator.integrate(Sfunc,triangle2)*coeff;
403  }
404  ++pb;
405  }
406  }
407 
408  template <typename T>
409  void N(const double coeff,T& matrix) const {
410  if (S_block_is_computed()) {
411  N(coeff/Scoeff,matrix,matrix);
412  } else {
413  Bloc Sbloc(mesh1.triangles().front().index(),mesh2.triangles().front().index(),mesh1.triangles().size(),mesh2.triangles().size());
414  S(1.0,Sbloc);
415  N(coeff,Sbloc,matrix);
416  }
417  }
418 
419  template <typename T>
420  void D(const double coeff,T& matrix) const {
421  base::message("D",mesh1,mesh2);
422  base::D(mesh1.triangles(),mesh2.triangles(),coeff,matrix);
423  }
424 
425  template <typename T>
426  void Dstar(const double coeff,T& matrix) const {
427  base::message("D*",mesh1,mesh2);
428  base::D(mesh2.triangles(),mesh1.triangles(),coeff,matrix);
429  }
430 
431  private:
432 
433  bool S_block_is_computed() const { return Scoeff!=0.0; }
434 
435  template <typename T1,typename T2>
436  void N(const double coeff,const T1& S,T2& matrix) const {
437 
438  base::message("N",mesh1,mesh2);
439 
440  ProgressBar pb(mesh1.vertices().size());
441  const VerticesRefs& m2_vertices = mesh2.vertices();
442  for (const auto& vertex1 : mesh1.vertices()) {
443  #pragma omp parallel for
444  #if defined NO_OPENMP || defined OPENMP_RANGEFOR
445  for (const auto& vertex2 : m2_vertices) {
446  #elif defined OPENMP_ITERATOR
447  for (auto vit2=m2_vertices.begin(); vit2<m2_vertices.end(); ++vit2) {
448  const Vertex* vertex2 = *vit2;
449  #else
450  for (int i2=0; i2<static_cast<int>(m2_vertices.size()); ++i2) {
451  const Vertex* vertex2 = *(m2_vertices.begin()+i2);
452  #endif
453  matrix(vertex1->index(),vertex2->index()) += base::N(*vertex1,*vertex2,mesh1,mesh2,S)*coeff;
454  }
455  ++pb;
456  }
457  }
458 
459  const Mesh& mesh1;
460  const Mesh& mesh2;
461  double Scoeff = 0.0;
462  };
463 
464  template <typename BlockType>
466  public:
467 
468  HeadMatrixBlocks(const BlockType& blk): block(blk) { }
469 
470  // SymMatrix is initialized at once, and there is nothing to for blockwise matrix.
471 
472  static void init(SymMatrix& matrix) { matrix.set(0.0); }
473  void set_blocks(SymMatrix&) const { }
474 
475  // SymmetricBlockMatrix is initialized blockwise.
476 
478  #if 0
479  void set_blocks(maths::SymmetricBlockMatrix& matrix) const { // Integrate this in blocks.
480  const Ranges& vrange1 = mesh1.vertices_ranges();
481  const Ranges& vrange2 = mesh2.vertices_ranges();
482  const Ranges& trange1 = { mesh1.triangles_range() };
483  const Ranges& trange2 = { mesh2.triangles_range() };
484 
485  matrix.add_blocks(trange1,trange2); // S blocks.
486  matrix.add_blocks(vrange1,vrange2); // N blocks.
487  matrix.add_blocks(trange1,vrange2); // D blocks.
488  if (!same_mesh)
489  matrix.add_blocks(trange2,vrange1); // D* blocks when they are not the transpose of D blocks.
490  }
491  #endif
492 
493  template <typename T>
494  void set_blocks(const double coeffs[3],T& matrix) {
495  const double SCondCoeff = coeffs[0];
496  const double NCondCoeff = coeffs[1];
497  const double DCondCoeff = coeffs[2];
498  block.set_S_block(SCondCoeff,matrix);
499  block.set_N_block(NCondCoeff,matrix);
500  block.set_D_block(DCondCoeff,matrix);
501  block.set_Dstar_block(DCondCoeff,matrix);
502  }
503 
504  private:
505 
506  BlockType block;
507  };
508 }
static double N(const Vertex &V1, const Vertex &V2, const Mesh &m1, const Mesh &m2, const T &matrix)
Definition: operators.h:116
BlocksBase(const Integrator &intg)
Definition: operators.h:62
static double N(const Vertex &V1, const Vertex &V2, const Mesh &m, const T &matrix)
Definition: operators.h:109
void message(const char *op_name, const Mesh &mesh) const
Definition: operators.h:64
void message(const char *op_name, const Mesh &mesh1, const Mesh &mesh2) const
Definition: operators.h:69
void D(const Triangles &triangles1, const Triangles &triangles2, const double coeff, T &mat) const
Definition: operators.h:77
const Integrator integrator
Definition: operators.h:143
void set_D_block(const double coeff, T &matrix) const
Definition: operators.h:183
void set_N_block(const double coeff, T &matrix) const
Definition: operators.h:180
void S(const double coeff, T &matrix) const
Definition: operators.h:201
void Dstar(const double coeff, T &matrix) const
Definition: operators.h:245
DiagonalBlock(const Mesh &m, const Integrator &intg)
Definition: operators.h:169
void N(const double coeff, T &matrix) const
Definition: operators.h:228
void D(const double coeff, T &matrix) const
Definition: operators.h:239
void set_S_block(const double coeff, T &matrix)
Definition: operators.h:172
void addId(const double coeff, T &matrix) const
Definition: operators.h:192
void set_Dstar_block(const double, T &) const
Definition: operators.h:189
Edge Edge class.
Definition: edge.h:18
const Vertex & vertex(const unsigned i) const
Definition: edge.h:31
static void init(maths::SymmetricBlockMatrix &)
Definition: operators.h:477
static void init(SymMatrix &matrix)
Definition: operators.h:472
HeadMatrixBlocks(const BlockType &blk)
Definition: operators.h:468
void set_blocks(SymMatrix &) const
Definition: operators.h:473
void set_blocks(const double coeffs[3], T &matrix)
Definition: operators.h:494
decltype(auto) integrate(const Function &function, const Triangle &triangle) const
Definition: integrator.h:45
Matrix class Matrix class.
Definition: matrix.h:28
std::string & name()
Definition: mesh.h:81
Triangles & triangles()
Definition: mesh.h:89
void Dstar(const double coeff, T &matrix) const
Definition: operators.h:426
void set_S_block(const double coeff, T &matrix)
Definition: operators.h:349
NonDiagonalBlock(const Mesh &m1, const Mesh &m2, const Integrator &intg)
Definition: operators.h:346
void S(const double coeff, T &matrix) const
Definition: operators.h:376
void D(const double coeff, T &matrix) const
Definition: operators.h:420
void set_N_block(const double coeff, T &matrix) const
Definition: operators.h:357
void set_D_block(const double coeff, T &matrix) const
Definition: operators.h:360
void set_Dstar_block(const double coeff, T &matrix) const
Definition: operators.h:366
void N(const double coeff, T &matrix) const
Definition: operators.h:409
void S(const double coeff, const Vertices &points, Matrix &matrix) const
Definition: operators.h:303
void addD(const double coeff, const Vertices &points, Matrix &matrix) const
Definition: operators.h:291
PartialBlock(const Mesh &m)
Definition: operators.h:289
void set(double x)
Triangle Triangle class.
Definition: triangle.h:45
Edge edge(const Vertex &V) const
Definition: triangle.h:86
unsigned & index()
Definition: triangle.h:101
double area() const
Definition: triangle.h:98
Vect3.
Definition: vect3.h:28
Vertex.
Definition: vertex.h:20
Vect3 f(const Vect3 &x) const
Definition: analytics.h:121
double f(const Vect3 &x) const
Definition: analytics.h:70
Block symmetric matrix class Block symmetric matrix class.
std::vector< Vertex > Vertices
Definition: vertex.h:37
void operatorDipolePot(const Dipole &, const Mesh &, Vector &, const double, const Integrator &)
std::vector< Triangle > Triangles
Definition: triangle.h:138
void operatorDipolePotDer(const Dipole &, const Mesh &, Vector &, const double, const Integrator &)
std::vector< Vertex * > VerticesRefs
Definition: vertex.h:38
void operatorFerguson(const Vect3 &, const Mesh &, Matrix &, const unsigned &, const double)
double dotprod(const Vect3 &V1, const Vect3 &V2)
Definition: vect3.h:106