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