42 for (
const auto&
tp : m.triangles(
V)) {
44 const Edge& edge =
T.edge(
V);
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
91 #elif defined OPENMP_ITERATOR
104 for (
unsigned i=0;
i<3; ++
i)
115 template <
typename T>
122 template <
typename T>
124 const double coeff = (
V1==
V2) ? 0.5 : 0.25;
130 template <
typename T>
134 for (
const auto&
tp1 :
m1.triangles(
V1)) {
137 for (
const auto&
tp2 :
m2.triangles(
V2)) {
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>
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>
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) {
219 #pragma omp parallel for
220 #if defined NO_OPENMP || defined OPENMP_ITERATOR
235 template <
typename T>
237 if (S_block_is_computed()) {
240 SymBloc
Sbloc(mesh.triangles().front().index(),mesh.triangles().size());
246 template <
typename T>
248 base::message(
"D",mesh,mesh);
249 base::D(mesh.triangles(),mesh.triangles(),
coeff,
matrix);
252 template <
typename T>
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) {
309 for (
unsigned i=0;
i<3;++
i)
317 for (
const auto& triangle : mesh.triangles()) {
319 for (
const auto& vertex :
points)
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>
370 template <
typename T>
372 if (!mesh1.current_barrier())
376 template <
typename T>
378 if (mesh1!=mesh2 && !mesh2.current_barrier())
386 template <
typename T>
388 base::message(
"S",mesh1,mesh2);
398 for (
const auto&
triangle1 : mesh1.triangles()) {
404 #pragma omp parallel for
405 #if defined NO_OPENMP || defined OPENMP_RANGEFOR
407 #elif defined OPENMP_ITERATOR
423 template <
typename T>
425 if (S_block_is_computed()) {
428 Bloc
Sbloc(mesh1.triangles().front().index(),mesh2.triangles().front().index(),mesh1.triangles().size(),mesh2.triangles().size());
434 template <
typename T>
436 base::message(
"D",mesh1,mesh2);
437 base::D(mesh1.triangles(),mesh2.triangles(),
coeff,
matrix);
440 template <
typename T>
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>
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>
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)
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)
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)