53 result += AB*analyS.f(x);
67 <<
"OPERATOR " << std::left << std::setw(2) << op_name
68 <<
"... (arg : mesh " << mesh.
name() <<
" )" << std::endl;
73 <<
"... (arg : mesh " << mesh1.
name() <<
" , mesh " << mesh2.
name() <<
" )"
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) {
95 for (
int i1=0; i1<static_cast<int>(triangles1.size()); ++i1) {
96 const Triangle& triangle1 = *(triangles1.begin()+i1);
99 for (
const auto& triangle2 : triangles2) {
101 const auto& Dfunc = [&analyD](
const Vect3& r) {
return analyD.
f(r); };
104 for (
unsigned i=0; i<3; ++i)
105 mat(triangle1.
index(),triangle2.vertex(i).index()) += total(i)*coeff;
115 template <
typename T>
117 return N(0.25,V1,V2,m,m,matrix);
122 template <
typename T>
124 const double coeff = (V1==V2) ? 0.5 : 0.25;
125 return N(coeff,V1,V2,m1,m2,matrix);
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) {
134 for (
const auto& tp1 : m1.
triangles(V1)) {
135 const Edge& edge1 = tp1->edge(V1);
137 for (
const auto& tp2 : m2.
triangles(V2)) {
138 const Edge& edge2 = tp2->edge(V2);
140 result -= factor*
dotprod(CB1,CB2)*matrix(tp1->index(),tp2->index())/(tp1->area()*tp2->area());
161 SymBloc(
const unsigned off,
const unsigned sz):
base(sz),offset(off) { }
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); }
168 const unsigned offset;
175 template <
typename T>
177 if (!mesh.current_barrier()) {
183 template <
typename T>
184 void set_N_block(
const double coeff,T& matrix)
const { N(coeff,matrix); }
186 template <
typename T>
188 if (!mesh.current_barrier())
192 template <
typename T>
195 template <
typename T>
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;
204 template <
typename T>
205 void S(
const double coeff,T& matrix)
const {
206 base::message(
"S",mesh,mesh);
213 const Triangles& triangles = mesh.triangles();
214 for (Triangles::const_iterator tit1=triangles.begin(); tit1!=triangles.end(); ++tit1,++pb) {
217 const auto& Sfunc = [&analyS](
const Vect3& r) {
return analyS.
f(r); };
219 #pragma omp parallel for
220 #if defined NO_OPENMP || defined OPENMP_ITERATOR
221 for (Triangles::const_iterator tit2=tit1; tit2!=triangles.end(); ++tit2) {
224 for (
int i2=tit1-triangles.begin(); i2<
static_cast<int>(triangles.size()); ++i2) {
225 const Triangle& triangle2 = *(triangles.begin()+i2);
228 matrix(triangle1.
index(),triangle2.
index()) = base::integrator.integrate(Sfunc,triangle2)*coeff;
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);
240 SymBloc Sbloc(mesh.triangles().front().index(),mesh.triangles().size());
242 N(coeff,Sbloc,matrix);
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);
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);
261 bool S_block_is_computed()
const {
return Scoeff!=0.0; }
263 static double Id(
const Triangle& T,
const Vertex& V) {
264 return (T.contains(V)) ? T.area()/3.0 : 0.0;
267 template <
typename T1,
typename T2>
268 void N(
const double coeff,
const T1& S,T2& matrix)
const {
270 base::message(
"N",mesh,mesh);
273 ProgressBar pb(mesh.vertices().size());
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) {
282 for (
int i2=0;i2<=vit1-mesh.vertices().begin();++i2) {
283 const auto vit2 = mesh.vertices().begin()+i2;
286 matrix((*vit1)->index(),(*vit2)->index()) += base::N(**vit1,**vit2,mesh,S)*coeff;
305 for (
const auto& triangle : mesh.triangles()) {
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;
317 for (
const auto& triangle : mesh.triangles()) {
319 for (
const auto& vertex : points)
320 matrix(vertex.index(),triangle.index()) = coeff*analyS.
f(vertex);
333 class Bloc:
public Matrix {
339 Bloc(
const unsigned r0,
const unsigned c0,
const unsigned n,
const unsigned m):
base(n,m),i0(r0),j0(c0) { }
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); }
359 template <
typename T>
361 if (!mesh1.current_barrier() && !mesh2.current_barrier()) {
367 template <
typename T>
368 void set_N_block(
const double coeff,T& matrix)
const { N(coeff,matrix); }
370 template <
typename T>
372 if (!mesh1.current_barrier())
376 template <
typename T>
378 if (mesh1!=mesh2 && !mesh2.current_barrier())
386 template <
typename T>
387 void S(
const double coeff,T& matrix)
const {
388 base::message(
"S",mesh1,mesh2);
398 for (
const auto& triangle1 : mesh1.triangles()) {
401 const auto& Sfunc = [&analyS](
const Vect3& r) {
return analyS.
f(r); };
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) {
411 for (
int i2=0;i2<static_cast<int>(m2_triangles.size());++i2) {
412 const Triangle& triangle2 = *(m2_triangles.begin()+i2);
415 matrix(triangle1.index(),triangle2.
index()) = base::integrator.integrate(Sfunc,triangle2)*coeff;
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);
428 Bloc Sbloc(mesh1.triangles().front().index(),mesh2.triangles().front().index(),mesh1.triangles().size(),mesh2.triangles().size());
430 N(coeff,Sbloc,matrix);
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);
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);
448 bool S_block_is_computed()
const {
return Scoeff!=0.0; }
450 template <
typename T1,
typename T2>
451 void N(
const double coeff,
const T1& S,T2& matrix)
const {
453 base::message(
"N",mesh1,mesh2);
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;
466 for (
int i2=0; i2<static_cast<int>(m2_vertices.size()); ++i2) {
467 const Vertex* vertex2 = *(m2_vertices.begin()+i2);
470 matrix(vertex1->index(),vertex2->index()) += base::N(*vertex1,*vertex2,mesh1,mesh2,S)*coeff;
483 template <
typename BlockType>
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() };
504 matrix.add_blocks(trange1,trange2);
505 matrix.add_blocks(vrange1,vrange2);
506 matrix.add_blocks(trange1,vrange2);
508 matrix.add_blocks(trange2,vrange1);
512 template <
typename T>
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);
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 addIdentity(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 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
std::ostream & log_stream(const InfoLevel level)
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)