52 result += AB*analyS.f(x);
66 std::cout <<
"OPERATOR " << std::left << std::setw(2) << op_name <<
"... (arg : mesh " << mesh.
name() <<
" )" << std::endl;
71 std::cout <<
"OPERATOR " << std::left << std::setw(2) << op_name <<
"... (arg : mesh " << mesh1.
name() <<
" , mesh " << mesh2.
name() <<
" )" << std::endl;
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) {
91 for (
int i1=0; i1<static_cast<int>(triangles1.size()); ++i1) {
92 const Triangle& triangle1 = *(triangles1.begin()+i1);
94 for (
const auto& triangle2 : triangles2) {
96 const auto& Dfunc = [&analyD](
const Vect3& r) {
return analyD.
f(r); };
99 for (
unsigned i=0; i<3; ++i)
100 mat(triangle1.
index(),triangle2.vertex(i).index()) += total(i)*coeff;
108 template <
typename T>
110 return N(0.25,V1,V2,m,m,matrix);
115 template <
typename T>
117 const double coeff = (V1==V2) ? 0.5 : 0.25;
118 return N(coeff,V1,V2,m1,m2,matrix);
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) {
127 for (
const auto& tp1 : m1.
triangles(V1)) {
128 const Edge& edge1 = tp1->edge(V1);
130 for (
const auto& tp2 : m2.
triangles(V2)) {
132 const Edge& edge2 = tp2->edge(V2);
135 result -= factor*
dotprod(CB1,CB2)*matrix(tp1->index(),tp2->index())/(tp1->area()*tp2->area());
157 SymBloc(
const unsigned off,
const unsigned sz):
base(sz),offset(off) { }
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); }
164 const unsigned offset;
171 template <
typename T>
173 if (!mesh.current_barrier()) {
179 template <
typename T>
180 void set_N_block(
const double coeff,T& matrix)
const { N(coeff,matrix); }
182 template <
typename T>
184 if (!mesh.current_barrier())
188 template <
typename T>
191 template <
typename T>
192 void addId(
const double coeff,T& matrix)
const {
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;
200 template <
typename T>
201 void S(
const double coeff,T& matrix)
const {
202 base::message(
"S",mesh,mesh);
208 const Triangles& triangles = mesh.triangles();
209 for (Triangles::const_iterator tit1=triangles.begin(); tit1!=triangles.end(); ++tit1,++pb) {
212 const auto& Sfunc = [&analyS](
const Vect3& r) {
return analyS.
f(r); };
214 #pragma omp parallel for
215 #if defined NO_OPENMP || defined OPENMP_ITERATOR
216 for (Triangles::const_iterator tit2=tit1; tit2!=triangles.end(); ++tit2) {
219 for (
int i2=tit1-triangles.begin(); i2<static_cast<int>(triangles.size()); ++i2) {
220 const Triangle& triangle2 = *(triangles.begin()+i2);
222 matrix(triangle1.
index(),triangle2.
index()) = base::integrator.integrate(Sfunc,triangle2)*coeff;
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);
232 SymBloc Sbloc(mesh.triangles().front().index(),mesh.triangles().size());
234 N(coeff,Sbloc,matrix);
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);
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);
253 bool S_block_is_computed()
const {
return Scoeff!=0.0; }
255 static double Id(
const Triangle& T,
const Vertex& V) {
256 return (T.contains(V)) ? T.area()/3.0 : 0.0;
259 template <
typename T1,
typename T2>
260 void N(
const double coeff,
const T1& S,T2& matrix)
const {
262 base::message(
"N",mesh,mesh);
264 ProgressBar pb(mesh.vertices().size());
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) {
273 for (
int i2=0;i2<=vit1-mesh.vertices().begin();++i2) {
274 const auto vit2 = mesh.vertices().begin()+i2;
276 matrix((*vit1)->index(),(*vit2)->index()) += base::N(**vit1,**vit2,mesh,S)*coeff;
292 std::cout <<
"PARTAL OPERATOR D..." << std::endl;
293 for (
const auto& triangle : mesh.triangles()) {
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;
304 std::cout <<
"PARTIAL OPERATOR S..." << std::endl;
305 for (
const auto& triangle : mesh.triangles()) {
307 for (
const auto& vertex : points)
308 matrix(vertex.index(),triangle.index()) = coeff*analyS.
f(vertex);
322 class Bloc:
public Matrix {
328 Bloc(
const unsigned r0,
const unsigned c0,
const unsigned n,
const unsigned m):
base(n,m),i0(r0),j0(c0) { }
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); }
348 template <
typename T>
350 if (!mesh1.current_barrier() && !mesh2.current_barrier()) {
356 template <
typename T>
357 void set_N_block(
const double coeff,T& matrix)
const { N(coeff,matrix); }
359 template <
typename T>
361 if (!mesh1.current_barrier())
365 template <
typename T>
367 if (mesh1!=mesh2 && !mesh2.current_barrier())
375 template <
typename T>
376 void S(
const double coeff,T& matrix)
const {
377 base::message(
"S",mesh1,mesh2);
386 for (
const auto& triangle1 : mesh1.triangles()) {
389 const auto& Sfunc = [&analyS](
const Vect3& r) {
return analyS.
f(r); };
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) {
399 for (
int i2=0;i2<static_cast<int>(m2_triangles.size());++i2) {
400 const Triangle& triangle2 = *(m2_triangles.begin()+i2);
402 matrix(triangle1.index(),triangle2.
index()) = base::integrator.integrate(Sfunc,triangle2)*coeff;
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);
413 Bloc Sbloc(mesh1.triangles().front().index(),mesh2.triangles().front().index(),mesh1.triangles().size(),mesh2.triangles().size());
415 N(coeff,Sbloc,matrix);
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);
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);
433 bool S_block_is_computed()
const {
return Scoeff!=0.0; }
435 template <
typename T1,
typename T2>
436 void N(
const double coeff,
const T1& S,T2& matrix)
const {
438 base::message(
"N",mesh1,mesh2);
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;
450 for (
int i2=0; i2<static_cast<int>(m2_vertices.size()); ++i2) {
451 const Vertex* vertex2 = *(m2_vertices.begin()+i2);
453 matrix(vertex1->index(),vertex2->index()) += base::N(*vertex1,*vertex2,mesh1,mesh2,S)*coeff;
464 template <
typename BlockType>
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() };
485 matrix.add_blocks(trange1,trange2);
486 matrix.add_blocks(vrange1,vrange2);
487 matrix.add_blocks(trange1,vrange2);
489 matrix.add_blocks(trange2,vrange1);
493 template <
typename T>
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);
static double N(const Vertex &V1, const Vertex &V2, const Mesh &m1, const Mesh &m2, const T &matrix)
BlocksBase(const Integrator &intg)
static double N(const Vertex &V1, const Vertex &V2, const Mesh &m, const T &matrix)
void message(const char *op_name, const Mesh &mesh) const
void message(const char *op_name, const Mesh &mesh1, const Mesh &mesh2) const
void D(const Triangles &triangles1, const Triangles &triangles2, const double coeff, T &mat) const
const Integrator integrator
void set_D_block(const double coeff, T &matrix) const
void set_N_block(const double coeff, T &matrix) const
void S(const double coeff, T &matrix) const
void Dstar(const double coeff, T &matrix) const
DiagonalBlock(const Mesh &m, const Integrator &intg)
void N(const double coeff, T &matrix) const
void D(const double coeff, T &matrix) const
void set_S_block(const double coeff, T &matrix)
void addId(const double coeff, T &matrix) const
void set_Dstar_block(const double, T &) const
const Vertex & vertex(const unsigned i) const
static void init(maths::SymmetricBlockMatrix &)
static void init(SymMatrix &matrix)
HeadMatrixBlocks(const BlockType &blk)
void set_blocks(SymMatrix &) const
void set_blocks(const double coeffs[3], T &matrix)
decltype(auto) integrate(const Function &function, const Triangle &triangle) const
Matrix class Matrix class.
void Dstar(const double coeff, T &matrix) const
void set_S_block(const double coeff, T &matrix)
NonDiagonalBlock(const Mesh &m1, const Mesh &m2, const Integrator &intg)
void S(const double coeff, T &matrix) const
void D(const double coeff, T &matrix) const
void set_N_block(const double coeff, T &matrix) const
void set_D_block(const double coeff, T &matrix) const
void set_Dstar_block(const double coeff, T &matrix) const
void N(const double coeff, T &matrix) const
void S(const double coeff, const Vertices &points, Matrix &matrix) const
void addD(const double coeff, const Vertices &points, Matrix &matrix) const
PartialBlock(const Mesh &m)
Edge edge(const Vertex &V) const
Vect3 f(const Vect3 &x) const
double f(const Vect3 &x) const
Block symmetric matrix class Block symmetric matrix class.
std::vector< Vertex > Vertices
void operatorDipolePot(const Dipole &, const Mesh &, Vector &, const double, const Integrator &)
std::vector< Triangle > Triangles
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)