62 void operatorFerguson(
const Vect3& ,
const Mesh& , Matrix& ,
const unsigned&,
const double&);
63 void operatorDipolePotDer(
const Vect3& ,
const Vect3& ,
const Mesh& , Vector&,
const double&,
const unsigned,
const bool);
64 void operatorDipolePot (
const Vect3& ,
const Vect3& ,
const Mesh& , Vector&,
const double&,
const unsigned,
const bool);
79 Vect3 total = gauss.integrate(analyD, T1);
82 for (
unsigned i = 0; i < 3; ++i)
83 mat(T1.
index(), T2(i).
index()) += total(i) * coeff;
91 for (
unsigned i=0;i<3;++i)
92 mat(P.
index(), T2(i).
index()) += total(i) * coeff;
110 return gauss.integrate(analyS, T2);
120 template <
typename T>
127 for (Mesh::VectPTriangle::const_iterator tit1 = trgs1.begin(); tit1 != trgs1.end(); ++tit1 ) {
128 for (Mesh::VectPTriangle::const_iterator tit2 = trgs2.begin(); tit2 != trgs2.end(); ++tit2 ) {
133 mat((*tit1)->index()-m1.begin()->index(),(*tit2)->index()-m2.begin()->index()) :
134 mat((*tit1)->index(),(*tit2)->index())/((*tit1)->area()*(*tit2)->area());
136 Vect3 CB1 = (*tit1)->next(V1)-(*tit1)->prev(V1);
137 Vect3 CB2 = (*tit2)->next(V2)-(*tit2)->prev(V2);
139 const double value = (CB1*CB2)*Iqr;
143 result -= (((&m1!=&m2) && (V1==V2)) ? 0.5 : 0.25)*value;
152 result = T2.
area()/3.0;
156 template <
typename T>
157 void operatorN(
const Mesh& m1,
const Mesh& m2,T& mat,
const double& coeff,
const unsigned gauss_order) {
164 std::cout <<
"OPERATOR N ... (arg : mesh " << m1.
name() <<
" , mesh " << m2.
name() <<
" )" << std::endl;
172 for (Mesh::const_iterator tit1=m1.begin();tit1!=m1.end();++tit1) {
174 #pragma omp parallel for
176 for (
int i2=tit1-m1.begin();i2<m1.size();++i2) {
177 const Mesh::const_iterator tit2 = m1.begin()+i2;
179 for (Mesh::const_iterator tit2=tit1;tit2<m1.end();++tit2) {
181 matS(tit1->index() - m1.begin()->index(), tit2->index() - m1.begin()->index()) =
_operatorS(*tit1, *tit2, gauss_order) / ( tit1->area() * tit2->area());
187 #pragma omp parallel for
194 mat((*vit1)->index(), (*vit2)->index()) +=
_operatorN(**vit1, **vit2, m1, m1, matS) * coeff;
200 #pragma omp parallel for
207 mat((*vit1)->index(), (*vit2)->index()) +=
_operatorN(**vit1, **vit2, m1, m1, mat) * coeff;
215 for (Mesh::const_iterator tit1 = m1.begin(); tit1 != m1.end(); ++tit1) {
217 #pragma omp parallel for
219 for (
int i2=0;i2<m2.size();++i2) {
220 const Mesh::const_iterator tit2 = m2.begin()+i2;
222 for (Mesh::const_iterator tit2=m2.begin();tit2<m2.end();++tit2) {
224 matS(tit1->index()-m1.begin()->index(),tit2->index()-m2.begin()->index()) =
_operatorS(*tit1,*tit2,gauss_order)/(tit1->area()*tit2->area());
230 #pragma omp parallel for
237 mat((*vit1)->index(),(*vit2)->index()) +=
_operatorN(**vit1,**vit2,m1,m2,matS)*coeff;
244 #pragma omp parallel for
251 mat((*vit1)->index(),(*vit2)->index()) +=
_operatorN(**vit1,**vit2,m1,m2,mat)*coeff;
258 template <
typename T>
259 void operatorS(
const Mesh& m1,
const Mesh& m2, T& mat,
const double& coeff,
const unsigned gauss_order)
267 std::cout <<
"OPERATOR S ... (arg : mesh " << m1.
name() <<
" , mesh " << m2.
name() <<
" )" << std::endl;
273 for (Mesh::const_iterator tit1=m1.begin();tit1!=m1.end();++tit1) {
275 #pragma omp parallel for
277 for (
int i2=tit1-m1.begin();i2<m1.size();++i2) {
278 const Mesh::const_iterator tit2 = m1.begin()+i2;
280 for (Mesh::const_iterator tit2=tit1;tit2<m1.end();++tit2) {
282 mat(tit1->index(),tit2->index()) =
_operatorS(*tit1,*tit2,gauss_order)*coeff;
289 for (Mesh::const_iterator tit1=m1.begin();tit1!=m1.end();++tit1) {
291 #pragma omp parallel for
293 for (
int i2=0;i2<m2.size();++i2) {
294 const Mesh::const_iterator tit2 = m2.begin()+i2;
296 for (Mesh::const_iterator tit2=m2.begin();tit2<m2.end();++tit2) {
298 mat(tit1->index(),tit2->index()) =
_operatorS(*tit1,*tit2,gauss_order)*coeff;
304 template <
typename T>
305 void operatorD(
const Mesh& m1,
const Mesh& m2, T& mat,
const double& coeff,
const unsigned gauss_order) {
318 #pragma omp parallel for
320 for (
int i1=0; i1 < m1.size(); ++i1) {
321 const Mesh::const_iterator tit1 = m1.begin()+i1;
323 for (Mesh::const_iterator tit1 = m1.begin(); tit1 < m1.end(); ++tit1) {
326 for (Mesh::const_iterator tit2=m2.begin();tit2!=m2.end();++tit2)
327 _operatorD(*tit1,*tit2,mat,coeff,gauss_order);
331 template <
typename T>
332 void operatorD(
const Mesh& m1,
const Mesh& m2, T& mat,
const double& coeff,
const unsigned gauss_order,
const bool star) {
340 std::cout <<
"OPERATOR D" << ((star) ?
"*" :
" ") <<
"... (arg : mesh " << m1.
name() <<
" , mesh " << m2.
name() <<
" )" << std::endl;
343 operatorD(m2, m1, mat, coeff, gauss_order);
345 operatorD(m1, m2, mat, coeff, gauss_order);
349 template <
typename T>
352 std::cout <<
"OPERATOR P1P0... (arg : mesh " << m.
name() <<
" )" << std::endl;
353 for (Mesh::const_iterator tit = m.begin(); tit != m.end(); ++tit)
355 mat(tit->index(), (*pit)->index()) +=
_operatorP1P0(*tit, **pit) * coeff;
359 STATIC_OMP
Vect3 result;
369 for (Mesh::VectPTriangle::const_iterator tit=trgs.begin();tit!=trgs.end();++tit) {
376 Vect3 A1B1 = (A1 - B1) * (0.5 / T1.
area());
378 analyS.init(V1, A1, B1);
379 const double opS = analyS.f(x);
381 result += (A1B1 * opS);
std::vector< Vertex > Vertices
double f(const Vect3 &x) const
double _operatorP1P0(const Triangle &T2, const Vertex &V1)
size_t nb_vertices() const
bool contains(const Vertex &p) const
Vect3 f(const Vect3 &x) const
Vect3 _operatorFerguson(const Vect3 &x, const Vertex &V1, const Mesh &m)
double _operatorSinternal(const Triangle &T, const Vertex &P)
vertex_iterator vertex_begin()
VectPVertex::const_iterator const_vertex_iterator
const Vertex & prev(const Vertex &V) const
const bool & current_barrier() const
double _operatorS(const Triangle &T1, const Triangle &T2, const unsigned gauss_order)
void operatorD(const Mesh &m1, const Mesh &m2, T &mat, const double &coeff, const unsigned gauss_order)
void _operatorD(const Triangle &T1, const Triangle &T2, T &mat, const double &coeff, const unsigned gauss_order)
void operatorP1P0(const Mesh &m, T &mat, const double &coeff)
vertex_iterator vertex_end()
void setOrder(const unsigned n)
size_t vertex_size() const
size_t nb_triangles() const
void operatorN(const Mesh &m1, const Mesh &m2, T &mat, const double &coeff, const unsigned gauss_order)
virtual T integrate(const I &fc, const Triangle &Trg)
void _operatorDinternal(const Triangle &T2, const Vertex &P, Matrix &mat, const double &coeff)
double _operatorN(const Vertex &V1, const Vertex &V2, const Mesh &m1, const Mesh &m2, const T &mat)
void init(const Triangle &T)
const Vertex & next(const Vertex &V) const
void operatorS(const Mesh &m1, const Mesh &m2, T &mat, const double &coeff, const unsigned gauss_order)
void operatorSinternal(const Mesh &, Matrix &, const Vertices &, const double &)
void operatorDipolePot(const Vect3 &, const Vect3 &, const Mesh &, Vector &, const double &, const unsigned, const bool)
void operatorFerguson(const Vect3 &, const Mesh &, Matrix &, const unsigned &, const double &)
const VectPTriangle & get_triangles_for_vertex(const Vertex &V) const
get the triangles associated with vertex V
#define PROGRESSBAR(a, b)
void operatorDinternal(const Mesh &, Matrix &, const Vertices &, const double &)
std::vector< Triangle * > VectPTriangle
void operatorDipolePotDer(const Vect3 &, const Vect3 &, const Mesh &, Vector &, const double &, const unsigned, const bool)