OpenMEEG
Loading...
Searching...
No Matches
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
27namespace 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
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
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 &)
void operatorFerguson(const Vect3 &, const Mesh &, Matrix &, const unsigned &, const double)
double dotprod(const Vect3 &V1, const Vect3 &V2)
Definition: vect3.h:106