From b5960dda1b9b8f43e1acc89bd3b24221f19968d4 Mon Sep 17 00:00:00 2001
From: Matthieu Aussal <aussal@mbp-de-matthieu.home>
Date: Fri, 20 Nov 2020 11:46:32 +0100
Subject: [PATCH] pre-release
---
demo/dirichlet.cpp | 15 +-
demo/dirichletBW.cpp | 11 +-
demo/head.cpp | 19 +-
demo/neumann.cpp | 13 +-
demo/neumannBW.cpp | 13 +-
demo/overview_bem.cpp | 124 ++-----
demo/overview_fem.cpp | 212 ++---------
include/bem.hpp | 737 +-------------------------------------
include/fem.hpp | 807 ++----------------------------------------
9 files changed, 121 insertions(+), 1830 deletions(-)
diff --git a/demo/dirichlet.cpp b/demo/dirichlet.cpp
index ba19399..5b119a8 100644
--- a/demo/dirichlet.cpp
+++ b/demo/dirichlet.cpp
@@ -8,8 +8,8 @@
| _#_ | AUTHOR(S) : Matthieu Aussal |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Acoustical scattering using BEM with Dirichlet|
- | `---' | boundary condition |
+ | ( === ) | SYNOPSIS : Acoustical scattering using H-Matrix BEM |
+ | `---' | solver for Dirichlet boundary condition |
+========================================================================+
*/
@@ -18,17 +18,18 @@
#include <castor/hmatrix.hpp>
#include <castor/linalg.hpp>
#include <castor/graphics.hpp>
-#include "bembuilder.hpp"
+#include "fem.hpp"
+#include "bem.hpp"
using namespace castor;
-int main (int argc, char* argv[])
+int main(int argc, char* argv[])
{
// Load meshes
matrix<double> Svtx, Pvtx;
matrix<size_t> Stri, Ptri;
- std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","unitSphere_1e3.ply");
- std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","Oxyz_1e4.ply");
+ std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","unitSphere_1e3.ply");
+ std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","Oxyz_1e4.ply");
// Graphical representation
figure fig;
@@ -83,7 +84,7 @@ int main (int argc, char* argv[])
// Solve linear system : S lambda = -PO
hmatrix<std::complex<double>> Lh,Uh;
tic();
- std::tie(Lh,Uh) = lu(Sbnd,tol);
+ std::tie(Lh,Uh) = lu(Sbnd);
toc();
disp(Lh);
disp(Uh);
diff --git a/demo/dirichletBW.cpp b/demo/dirichletBW.cpp
index 23ad4dd..c95e45f 100644
--- a/demo/dirichletBW.cpp
+++ b/demo/dirichletBW.cpp
@@ -8,8 +8,8 @@
| _#_ | AUTHOR(S) : Matthieu Aussal |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Acoustical scattering using BEM with Dirichlet|
- | `---' | boundary condition |
+ | ( === ) | SYNOPSIS : Acoustical scattering using H-Matrix BEM |
+ | `---' | solver for Dirichlet boundary condition |
+========================================================================+
*/
@@ -18,7 +18,8 @@
#include <castor/hmatrix.hpp>
#include <castor/linalg.hpp>
#include <castor/graphics.hpp>
-#include "bembuilder.hpp"
+#include "fem.hpp"
+#include "bem.hpp"
using namespace castor;
@@ -27,8 +28,8 @@ int main (int argc, char* argv[])
// Load meshes
matrix<double> Svtx, Pvtx;
matrix<size_t> Stri, Ptri;
- std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","unitSphere_1e3.ply");
- std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","Oxyz_1e4.ply");
+ std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","unitSphere_1e3.ply");
+ std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","Oxyz_1e4.ply");
// Graphical representation
figure fig;
diff --git a/demo/head.cpp b/demo/head.cpp
index 7ba825a..fe18bfa 100644
--- a/demo/head.cpp
+++ b/demo/head.cpp
@@ -8,9 +8,9 @@
| _#_ | AUTHOR(S) : Matthieu Aussal |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Acoustical scattering using BEM with Neumann |
- | `---' | boundary condition |
- +========================================================================+
+ | ( === ) | SYNOPSIS : Acoustical scattering using H-Matrix BEM |
+ | `---' | solver for Neumann boundary condition |
+ +========================================================================+=+
*/
#include <castor/matrix.hpp>
@@ -18,7 +18,8 @@
#include <castor/hmatrix.hpp>
#include <castor/linalg.hpp>
#include <castor/graphics.hpp>
-#include "bembuilder.hpp"
+#include "fem.hpp"
+#include "bem.hpp"
using namespace castor;
@@ -27,7 +28,7 @@ int main (int argc, char* argv[])
// Load meshes
matrix<double> Svtx;
matrix<size_t> Stri;
- std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","Head03_04kHz.ply");
+ std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","Head03_04kHz.ply");
// Graphical representation
figure fig;
@@ -35,18 +36,14 @@ int main (int argc, char* argv[])
// Parameters
matrix<double> U = {0,0,-1};
- double f = 4000;
+ double f = 2000;
double k = 2*M_PI*f/340;
float tol = 1e-3;
- // FEM
+ // FEM and identity matrix
tic();
femdata<double> v(Stri,Svtx,lagrangeP1,3);
femdata<double> u(Stri,Svtx,lagrangeP1,3);
- toc();
-
- // Identity matrix
- tic();
auto Id = mass<std::complex<double>>(v);
toc();
diff --git a/demo/neumann.cpp b/demo/neumann.cpp
index 4d8433e..b11ab93 100644
--- a/demo/neumann.cpp
+++ b/demo/neumann.cpp
@@ -8,9 +8,9 @@
| _#_ | AUTHOR(S) : Matthieu Aussal |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Acoustical scattering using BEM with Neumann |
- | `---' | boundary condition |
- +========================================================================+
+ | ( === ) | SYNOPSIS : Acoustical scattering using H-Matrix BEM |
+ | `---' | solver for Neumann boundary condition |
+ +========================================================================+=+
*/
#include <castor/matrix.hpp>
@@ -18,7 +18,8 @@
#include <castor/hmatrix.hpp>
#include <castor/linalg.hpp>
#include <castor/graphics.hpp>
-#include "bembuilder.hpp"
+#include "fem.hpp"
+#include "bem.hpp"
using namespace castor;
@@ -27,8 +28,8 @@ int main (int argc, char* argv[])
// Load meshes
matrix<double> Svtx, Pvtx;
matrix<size_t> Stri, Ptri;
- std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","unitSphere_1e3.ply");
- std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","Oxyz_1e4.ply");
+ std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","unitSphere_1e3.ply");
+ std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","Oxyz_1e4.ply");
// Graphical representation
figure fig;
diff --git a/demo/neumannBW.cpp b/demo/neumannBW.cpp
index 56b52fc..ac94c77 100644
--- a/demo/neumannBW.cpp
+++ b/demo/neumannBW.cpp
@@ -8,9 +8,9 @@
| _#_ | AUTHOR(S) : Matthieu Aussal |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Acoustical scattering using BEM with Neumann |
- | `---' | boundary condition |
- +========================================================================+
+ | ( === ) | SYNOPSIS : Acoustical scattering using H-Matrix BEM |
+ | `---' | solver for Neumann boundary condition |
+ +========================================================================+=+
*/
#include <castor/matrix.hpp>
@@ -18,7 +18,8 @@
#include <castor/hmatrix.hpp>
#include <castor/linalg.hpp>
#include <castor/graphics.hpp>
-#include "bembuilder.hpp"
+#include "fem.hpp"
+#include "bem.hpp"
using namespace castor;
@@ -27,8 +28,8 @@ int main (int argc, char* argv[])
// Load meshes
matrix<double> Svtx, Pvtx;
matrix<size_t> Stri, Ptri;
- std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","unitSphere_1e3.ply");
- std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/bembuilder/data/","Oxyz_1e4.ply");
+ std::tie(Stri,Svtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","unitSphere_1e3.ply");
+ std::tie(Ptri,Pvtx) = triread("/Users/aussal/Git/leprojetcastor/fembem/demo/data/","Oxyz_1e4.ply");
// Graphical representation
figure fig;
diff --git a/demo/overview_bem.cpp b/demo/overview_bem.cpp
index 1e87992..3f852bc 100644
--- a/demo/overview_bem.cpp
+++ b/demo/overview_bem.cpp
@@ -3,13 +3,13 @@
| (c) 2020 - PROPERTY OF ECOLE POLYTECHNIQUE - LGPL 3.0 |
|________________________________________________________________________|
| '&` | |
- | # | FILE : overview_bembuilder.cpp |
+ | # | FILE : overview_bem.cpp |
| # | VERSION : 0.1.0 |
| _#_ | AUTHOR(S) : Matthieu Aussal & Marc Bakry |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Boundary Element Methods with Finite elements |
- | `---' | representation |
+ | ( === ) | SYNOPSIS : Boundary Elements Methods for standard |
+ | `---' | operators S, D, Dt, H in dense format |
+========================================================================+
*/
@@ -17,66 +17,13 @@
#include <castor/smatrix.hpp>
#include <castor/linalg.hpp>
#include <castor/graphics.hpp>
-#include "bembuilder.hpp"
+#include "fem.hpp"
+#include "bem.hpp"
using namespace castor;
-int main (int argc, char* argv[])
+int main(int argc, char* argv[])
{
- //===============================================================
- std::cout << "+=====================+" << std::endl;
- std::cout << "| INITIALIZE |" << std::endl;
- std::cout << "+=====================+" << std::endl;
-
- // Documentation
- documentationFiles =
- {
- "/Users/aussal/Git/leprojetcastor/matrix/include/matrix/matrix.hpp",
- "/Users/aussal/Git/leprojetcastor/graphics/graphics.hpp"
- };
- help("");
-
- //===============================================================
- std::cout << "+====================+" << std::endl;
- std::cout << "| TRIANGLEDATA |" << std::endl;
- std::cout << "+====================+" << std::endl;
-
- // Unitary triangle mesh
- matrix<double> Tvtx = {{0,0,0},{2,0,0},{0,2,0}};
- matrix<std::size_t> Ttri = {0,1,2};
- matrix<double> ctredg = {{1,1,0},{0,1,0},{1,0,0}};
-
- // Triangle data
- triangledata<double> triangle;
- triangle.update(Tvtx);
- disp(triangle.isfar(0,0,3));
-
- // Graphical representation
- figure fig;
- caxis(fig,{0,1});
- trimesh(fig,Ttri,Tvtx);
- quiver(fig,triangle.ctr(),triangle.nrm());
- quiver(fig,ctredg,triangle.tau());
- quiver(fig,ctredg,triangle.nu());
-
- //===============================================================
- std::cout << "+====================+" << std::endl;
- std::cout << "| DOF/QUADRATURE |" << std::endl;
- std::cout << "+====================+" << std::endl;
-
- // Finite element data
- femdata<double> fem(Ttri,Tvtx,lagrangeP1,3);
- disp(sum(fem.wgt()) - 2);
- disp(fem.base(0,0)+fem.base(0,1)+fem.base(0,2));
- disp(fem.grad(0,0,0)+fem.grad(0,1,0)+fem.grad(0,2,0));
-
- // Graphical representation
- caxis(fig,{0,1});
- std::size_t nd = size(fem.dof(),1), ng = size(fem.gss(),1);
- vermesh(fig,range(0,nd),fem.dof(),zeros(nd,1));
- vermesh(fig,range(0,ng),fem.gss(),ones(ng,1));
- quiver(fig,fem.gss(),0.2*fem.gssnrm());
-
//===============================================================
std::cout << "+===================+" << std::endl;
std::cout << "| MESHES |" << std::endl;
@@ -97,38 +44,9 @@ int main (int argc, char* argv[])
std::tie(Ptri,Pvtx) = tridelaunay(X,Y,zeros(size(X)));
// Graphical representation
- figure fig2;
- trimesh(fig2,Stri,Svtx);
- trimesh(fig2,Ptri,Pvtx);
-
- //===============================================================
- std::cout << "+======================+" << std::endl;
- std::cout << "| TRIANGLE INTEGRATION |" << std::endl;
- std::cout << "+======================+" << std::endl;
-
- // Exact integration VS singular
- X = 3*Svtx;
- matrix<double> ref(size(X,1),2), sol(size(X,1),2);
- double rxy, rdotn;
- for (std::size_t l=0; l<size(X,1); ++l)
- {
- // Analytic
- std::tie(ref(l,0),ref(l,1)) = triangle.intsing(X(l,0),X(l,1),X(l,2));
-
- // Numeric (gauss)
- for (std::size_t g=0; g<size(fem.gss(),1); ++g)
- {
- rxy = std::sqrt(std::pow(X(l,0)-fem.gss(g,0),2) +
- std::pow(X(l,1)-fem.gss(g,1),2) +
- std::pow(X(l,2)-fem.gss(g,2),2) );
- rdotn = ((X(l,0)-fem.gss(g,0))*fem.gssnrm(g,0) +
- (X(l,1)-fem.gss(g,1))*fem.gssnrm(g,1) +
- (X(l,2)-fem.gss(g,2))*fem.gssnrm(g,2) );
- sol(l,0) += 1./rxy * fem.wgt(g);
- sol(l,1) += rdotn/std::pow(rxy,3) * fem.wgt(g);
- }
- }
- disp(norm(ref-sol,"inf")/norm(ref,"inf"));
+ figure fig;
+ trimesh(fig,Stri,Svtx);
+ trimesh(fig,Ptri,Pvtx);
//===============================================================
std::cout << "+=======================+" << std::endl;
@@ -175,9 +93,9 @@ int main (int argc, char* argv[])
toc();
// Graphical representation
- figure fig5;
- trimesh(fig5,Ptri,Pvtx,abs(Pdom));
- trimesh(fig5,Stri,Svtx,abs(Pbnd));
+ figure fig2;
+ trimesh(fig2,Ptri,Pvtx,abs(Pdom));
+ trimesh(fig2,Stri,Svtx,abs(Pbnd));
//===============================================================
std::cout << "+=======================+" << std::endl;
@@ -208,9 +126,9 @@ int main (int argc, char* argv[])
}
// Graphical representation
- figure fig6;
- trimesh(fig6,Ptri,Pvtx,abs(Pdom));
- trimesh(fig6,Stri,Svtx,abs(Pbnd));
+ figure fig3;
+ trimesh(fig3,Ptri,Pvtx,abs(Pdom));
+ trimesh(fig3,Stri,Svtx,abs(Pbnd));
//===============================================================
std::cout << "+========================+" << std::endl;
@@ -242,9 +160,9 @@ int main (int argc, char* argv[])
}
// Graphical representation
- figure fig7;
- trimesh(fig7,Stri,Svtx,abs(Pbnd));
- trimesh(fig7,Ptri,Pvtx,abs(Pdom));
+ figure fig4;
+ trimesh(fig4,Stri,Svtx,abs(Pbnd));
+ trimesh(fig4,Ptri,Pvtx,abs(Pdom));
//===============================================================
std::cout << "+===================+" << std::endl;
@@ -276,9 +194,9 @@ int main (int argc, char* argv[])
}
// Graphical representation
- figure fig8;
- trimesh(fig8,Stri,Svtx,abs(Pbnd));
- trimesh(fig8,Ptri,Pvtx,abs(Pdom));
+ figure fig5;
+ trimesh(fig5,Stri,Svtx,abs(Pbnd));
+ trimesh(fig5,Ptri,Pvtx,abs(Pdom));
//===============================================================
std::cout << "+====================+" << std::endl;
diff --git a/demo/overview_fem.cpp b/demo/overview_fem.cpp
index c9da151..160fa26 100644
--- a/demo/overview_fem.cpp
+++ b/demo/overview_fem.cpp
@@ -3,13 +3,13 @@
| (c) 2020 - PROPERTY OF ECOLE POLYTECHNIQUE - LGPL 3.0 |
|________________________________________________________________________|
| '&` | |
- | # | FILE : overview_bembuilder.cpp |
+ | # | FILE : overview_fem.cpp |
| # | VERSION : 0.1.0 |
| _#_ | AUTHOR(S) : Matthieu Aussal & Marc Bakry |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Boundary Element Methods with Finite elements |
- | `---' | representation |
+ | ( === ) | SYNOPSIS : Mesh, quadrature and finite elements tools |
+ | `---' | |
+========================================================================+
*/
@@ -18,25 +18,11 @@
#include <castor/linalg.hpp>
#include <castor/graphics.hpp>
#include "fem.hpp"
-//#include "bem.hpp"
using namespace castor;
-int main (int argc, char* argv[])
+int main(int argc, char* argv[])
{
- //===============================================================
- std::cout << "+=====================+" << std::endl;
- std::cout << "| INITIALIZE |" << std::endl;
- std::cout << "+=====================+" << std::endl;
-
- // Documentation
- documentationFiles =
- {
- "/Users/aussal/Git/leprojetcastor/matrix/include/matrix/matrix.hpp",
- "/Users/aussal/Git/leprojetcastor/graphics/graphics.hpp"
- };
- help("");
-
//===============================================================
std::cout << "+====================+" << std::endl;
std::cout << "| TRIANGLEDATA |" << std::endl;
@@ -92,15 +78,10 @@ int main (int argc, char* argv[])
Z = horzcat(0,reshape(Z,1,numel(Z)));
std::tie(Stri,Svtx) = tetdelaunay(X,Y,Z);
std::tie(Stri,Svtx) = tetboundary(Stri,Svtx);
-
- // Planar mesh
- std::tie(X,Y) = meshgrid(linspace(-4,4,50));
- std::tie(Ptri,Pvtx) = tridelaunay(X,Y,zeros(size(X)));
// Graphical representation
figure fig2;
trimesh(fig2,Stri,Svtx);
- trimesh(fig2,Ptri,Pvtx);
//===============================================================
std::cout << "+======================+" << std::endl;
@@ -130,181 +111,34 @@ int main (int argc, char* argv[])
}
}
disp(norm(ref-sol,"inf")/norm(ref,"inf"));
-
- //===============================================================
- std::cout << "+=======================+" << std::endl;
- std::cout << "| SINGLE LAYER |" << std::endl;
- std::cout << "+=======================+" << std::endl;
-
- // FEM
- femdata<double> v(Stri,Svtx,lagrangeP1,3);
- femdata<double> u(Stri,Svtx,lagrangeP1,3);
- // Incident wave : PW(x,u) = exp(i*k*x.u)
- matrix<double> U = {-1,0,0};
- double k = 5;
- auto Pinc = planeWave(Pvtx,U,k);
-
- // Identity matrix
- tic();
- auto Id = full(mass<std::complex<double>>(v));
- toc();
-
- // Left hand side : \int_sigma_x \int_sigma_y v(x) G(x,y) u(y) dx dy
- tic();
- auto LHS = singleLayer<std::complex<double>>(v,u,k);
- toc();
- disp(sum(LHS)); // -0.6614 + 2.3185i
-
- // Right hand side : - \int_sigma_x v(x) PW(x,u) dx
- tic();
- auto RHS = - rightHandSide<std::complex<double>>(v,PWsource,U,k);
- toc();
-
- // Solve linear system : S lambda = -PO
- tic();
- auto lambda = linsolve(LHS,RHS);
- toc();
-
- // Radiation
- tic();
- auto Pbnd = linsolve(Id,mtimes(LHS,lambda)) + planeWave(v.dof(),U,k);
- toc();
- tic();
- auto Sdom = singleLayer<std::complex<double>>(Pvtx,u,k);
- auto Pdom = mtimes(Sdom,lambda) + Pinc;
- toc();
-
- // Graphical representation
- figure fig5;
- trimesh(fig5,Ptri,Pvtx,abs(Pdom));
- trimesh(fig5,Stri,Svtx,abs(Pbnd));
-
- //===============================================================
- std::cout << "+=======================+" << std::endl;
- std::cout << "| DOUBLE LAYER |" << std::endl;
- std::cout << "+=======================+" << std::endl;
-
- // Left hand side : \int_sigma_x \int_sigma_y v(x) dnyG(x,y) u(y) dx dy
- tic();
- LHS = -(0.5*Id + doubleLayer<std::complex<double>>(v,u,k));
- toc();
- disp(sum(LHS)); // -1.6068 + 5.6271i
-
- // Right hand side : - \int_sigma_x v(x) PW(x,u) dx
- RHS = - rightHandSide<std::complex<double>>(v,PWsource,U,k);
-
- // Solve -(Id/2 - D) = -P0
- auto mu = linsolve(LHS,RHS);
-
- // Boundary radiation
- Pbnd = linsolve(Id,mtimes(LHS,mu)) + planeWave(v.dof(),U,k);
-
- // Domain radiation
- auto Ddom = doubleLayer<std::complex<double>>(Pvtx,u,k);
- Pdom = - mtimes(Ddom,mu) + Pinc;
- for (std::size_t l=0; l<size(Pvtx,1); ++l)
- {
- if (std::pow(Pvtx(l,0),2)+std::pow(Pvtx(l,1),2)+std::pow(Pvtx(l,2),2)<1.1) {Pdom(l)=Pinc(l);}
- }
-
- // Graphical representation
- figure fig6;
- trimesh(fig6,Ptri,Pvtx,abs(Pdom));
- trimesh(fig6,Stri,Svtx,abs(Pbnd));
-
//===============================================================
- std::cout << "+========================+" << std::endl;
- std::cout << "| DOUBLE LAYER TRANSPOSE |" << std::endl;
- std::cout << "+========================+" << std::endl;
-
- // Left hand side : \int_sigma_x \int_sigma_y v(x) dnyG(x,y) u(y) dx dy
+ std::cout << "+=========================+" << std::endl;
+ std::cout << "| VARIATIONAL FORMULATION |" << std::endl;
+ std::cout << "+=========================+" << std::endl;
+
+ // Finite element space
+ fem = femdata<double>(Stri,Svtx,lagrangeP1,3);
+
+ // Mass matrix
tic();
- LHS = -0.5*Id + transpose(doubleLayer<std::complex<double>>(v,u,k));
+ auto Ms = mass<double>(fem);
toc();
- disp(sum(LHS)); //
-
- // Right hand side : - \int_sigma_x v(x) PW(x,u) dx
- RHS = - rightHandSide<std::complex<double>>(v,dnPWsource,U,k);
-
- // Solve (-Id/2 + Dt) = -Pdn0
- lambda = linsolve(LHS,RHS);
-
- // Boundary radiation
- auto Sbnd = singleLayer<std::complex<double>>(v,u,k);
- Pbnd = linsolve(Id,mtimes(Sbnd,lambda)) + planeWave(v.dof(),U,k);
-
- // Domain radiation
- Sdom = singleLayer<std::complex<double>>(Pvtx,u,k);
- Pdom = mtimes(Sdom,lambda) + Pinc;
- for (std::size_t l=0; l<size(Pvtx,1); ++l)
- {
- if (std::pow(Pvtx(l,0),2)+std::pow(Pvtx(l,1),2)+std::pow(Pvtx(l,2),2)<1.1) {Pdom(l)=Pinc(l);}
- }
-
- // Graphical representation
- figure fig7;
- trimesh(fig7,Stri,Svtx,abs(Pbnd));
- trimesh(fig7,Ptri,Pvtx,abs(Pdom));
-
- //===============================================================
- std::cout << "+===================+" << std::endl;
- std::cout << "| HYPERSINGULAR |" << std::endl;
- std::cout << "+===================+" << std::endl;
-
- // Left hand side : k^2 * \int_Sx \int_Sy n.psi(x) G(x,y) n.psi(y) dx dy - \int_Sx \int_Sy nxgrad(psi(x)) G(x,y) nxgrad(psi(y)) dx dy
+ disp(4*M_PI - sum(full(Ms)));
+
+ // Rigidity matrix
tic();
- LHS = - hypersingular<std::complex<double>>(v,u,k);
+ auto Ks = rigidity<double>(fem);
toc();
- disp(sum(LHS)); // -22.4022 -13.6573i
-
- // Right hand side : - \int_sigma_x v(x) PW(x,u) dx
- RHS = - rightHandSide<std::complex<double>>(v,dnPWsource,U,k);
-
- // Solve -H = -dnP0
- mu = linsolve(LHS,RHS);
-
- // Boundary radiation
- auto Dbnd = 0.5*Id + doubleLayer<std::complex<double>>(v,u,k);
- Pbnd = -linsolve(Id,mtimes(Dbnd,mu)) + planeWave(v.dof(),U,k);
-
- // Domain radiation
- Ddom = doubleLayer<std::complex<double>>(Pvtx,u,k);
- Pdom = - mtimes(Ddom,mu) + Pinc;
- for (std::size_t l=0; l<size(Pvtx,1); ++l)
- {
- if (std::pow(Pvtx(l,0),2)+std::pow(Pvtx(l,1),2)+std::pow(Pvtx(l,2),2)<1.1) {Pdom(l)=Pinc(l);}
- }
-
- // Graphical representation
- figure fig8;
- trimesh(fig8,Stri,Svtx,abs(Pbnd));
- trimesh(fig8,Ptri,Pvtx,abs(Pdom));
+ disp(max(full(Ks)));
+ disp(sum(full(Ks)));
- //===============================================================
- std::cout << "+====================+" << std::endl;
- std::cout << "| SUB-MATRIX |" << std::endl;
- std::cout << "+====================+" << std::endl;
-
- // Collocation
- matrix<std::size_t> I = colon(30,2,900);
- matrix<std::size_t> J = colon(900,-3,7);
- matrix<double> S = singleLayer<double>(Svtx,u,0.);
- matrix<double> Sb = singleLayer<double>(Svtx,u,0.,I,J);
- disp(norm(eval(S(I,J))-Sb,"inf"));
-
- // Galerkin
- S = singleLayer<double>(v,u,0.);
- Sb = singleLayer<double>(v,u,0.,I,J);
- disp(norm(eval(S(I,J))-Sb,"inf"));
-
- //===============================================================
- std::cout << "+===================+" << std::endl;
- std::cout << "| DRAWNOW |" << std::endl;
- std::cout << "+===================+" << std::endl;
+ // Graphical representation
+ figure fig3;
+ spy(fig3,Ms,"b","Mass matrix");
- drawnow(fig);
+ drawnow(fig);
disp("done !");
return 0;
}
diff --git a/include/bem.hpp b/include/bem.hpp
index 8802bde..1a4a70d 100644
--- a/include/bem.hpp
+++ b/include/bem.hpp
@@ -3,721 +3,35 @@
| (c) 2020 - PROPERTY OF ECOLE POLYTECHNIQUE - LGPL 3.0 |
|________________________________________________________________________|
| '&` | |
- | # | FILE : bembuilder.hpp |
+ | # | FILE : bem.hpp |
| # | VERSION : 0.1.0 |
| _#_ | AUTHOR(S) : Matthieu Aussal & Marc Bakry |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Boundary Element Methods with Finite elements |
- | `---' | representation |
+ | ( === ) | SYNOPSIS : Boundary Elements Methods for standard |
+ | `---' | operators S, D, Dt, H in dense format |
+========================================================================+
*/
#pragma once
+
+#define CASTOR_BEM_HPP
+
#include <castor/matrix.hpp>
#include <castor/smatrix.hpp>
+#include "fem.hpp"
-using namespace castor;
+namespace castor
+{
//==========================================================================//
// ENVIRONMENT DATA //
//==========================================================================//
-enum femoperator{none, gradient};
-enum femtype{vertices, lagrangeP0, lagrangeP1};
enum bemsource{PWsource, dnPWsource, SWsource, dnSWsource};
//==========================================================================//
-// TRIANGLEDATA CLASS //
-//==========================================================================//
-//==========================================================================
-// [triangledata]
-///
-template<typename T>
-class triangledata
-{
-public:
- auto intsing(T x0, T y0, T z0) const;
- bool isfar(T x0, T y0, T z0) const;
- void projectP1(T x0, T y0, T z0, matrix<T>& P) const;
- void update(matrix<T>const& S);
- void update(matrix<T>const& A, matrix<T>const& B, matrix<T>const& C);
-
- matrix<T>const& ctr() const {return m_ctr;};
- matrix<T>const& nrm() const {return m_nrm;};
- matrix<T>const& nu() const {return m_nu;};
- matrix<T>const& tau() const {return m_tau;};
- T const& A(std::size_t i) const {return m_A(i);};
- T const& B(std::size_t i) const {return m_B(i);};
- T const& C(std::size_t i) const {return m_C(i);};
- T const& ctr(std::size_t i) const {return m_ctr(i);};
- T const& hgt(std::size_t i) const {return m_hgt(i);};
- T const& nrm(std::size_t i) const {return m_nrm(i);};
- T const& nu(std::size_t i, std::size_t j) const {return m_nu(i,j);};
- T const& srf() const {return m_srf;}
- T const& tau(std::size_t i, std::size_t j) const {return m_tau(i,j);};
-
-private:
- matrix<T> m_A = matrix<T>(1,3); // FIRST NODE
- matrix<T> m_B = matrix<T>(1,3); // SECOND NODE
- matrix<T> m_C = matrix<T>(1,3); // THIRD NODE
- matrix<T> m_ctr = matrix<T>(1,3); // CENTER
- matrix<T> m_hgt = matrix<T>(1,3); // HEIGHT
- matrix<T> m_lgt = matrix<T>(1,3); // EDGE LENGTH
- matrix<T> m_nrm = matrix<T>(1,3); // TRIANGLE NORMAL
- matrix<T> m_nu = matrix<T>(3,3); // EDGE NORMALS
- T m_srf = 0; // SURFACE
- T m_stp = 0; // MAX EDGE LENGTH
- matrix<T> m_tau = matrix<T>(3,3); // EDGE TANGENT
-};
-
-//==========================================================================
-//
-///
-template<typename T>
-auto triangledata<T>::intsing(T x0, T y0, T z0) const
-{
- // initialize
- std::size_t ip, ipp;
- static matrix<T> xX0S(1,3), yX0S(1,3), zX0S(1,3);
- static matrix<T> nrmX0S(1,3), ca(1,3), sa(1,3), xnu(1,3);
- T omega, h, ps, psp1, ps2, ah, ar, intaRm1;
- T Rm1, dnRm1;
- T tol = 1e-7;
-
- // Point to vertex
- xX0S(0) = m_A(0) - x0;
- xX0S(1) = m_B(0) - x0;
- xX0S(2) = m_C(0) - x0;
-
- yX0S(0) = m_A(1) - y0;
- yX0S(1) = m_B(1) - y0;
- yX0S(2) = m_C(1) - y0;
-
- zX0S(0) = m_A(2) - z0;
- zX0S(1) = m_B(2) - z0;
- zX0S(2) = m_C(2) - z0;
-
- for(std::size_t i=0; i<3; ++i)
- {
- nrmX0S(i) = std::sqrt(xX0S(i)*xX0S(i)+yX0S(i)*yX0S(i)+zX0S(i)*zX0S(i));
- }
-
- // Point height to triangle
- h = xX0S(0)*m_nrm(0) + yX0S(0)*m_nrm(1) + zX0S(0)*m_nrm(2);
-
- // solid angle
- if (std::abs(h)<tol) {omega = 0;}
- else
- {
- for(std::size_t i=0; i<3; ++i)
- {
- ip = (i+1)%3;
- ca(i) = (xX0S(i)*xX0S(ip)+yX0S(i)*yX0S(ip)+zX0S(i)*zX0S(ip))/(nrmX0S(i)*nrmX0S(ip));
- sa(i) = std::sqrt(1 - ca(i)*ca(i));
- }
- omega = -M_PI;
- for(std::size_t i=0; i<3; ++i)
- {
- ip = (i+1)%3; ipp = (ip+1)%3;
- omega += std::acos((ca(i)-ca(ip)*ca(ipp))/(sa(ip)*sa(ipp)));
- }
- if (h < 0) {omega = -omega;}
- }
- if (std::isnan(omega))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Point is close to triangle, change h tolerance.");
- }
-
- // initialization of integration
- Rm1 = -h*omega;
- dnRm1 = -omega;
-
- // edge integration
- for(std::size_t i=0; i<3; ++i)
- {
- ip = (i+1)%3;
- ipp = (ip+1)%3;
- // projections
- ps = xX0S(ip)*m_tau(i,0) + yX0S(ip)*m_tau(i,1) + zX0S(ip)*m_tau(i,2);
- psp1 = xX0S(ipp)*m_tau(i,0) + yX0S(ipp)*m_tau(i,1) + zX0S(ipp)*m_tau(i,2);
- ps2 = ps*(nrmX0S(ipp) - nrmX0S(ip));
- // Subtil
- ar = std::abs(psp1 - ps);
- xnu(0) = xX0S(ip) - ps*m_tau(i,0);
- xnu(1) = yX0S(ip) - ps*m_tau(i,1);
- xnu(2) = zX0S(ip) - ps*m_tau(i,2);
- ah = std::sqrt(xnu(0)*xnu(0) + xnu(1)*xnu(1) + xnu(2)*xnu(2));
- // integration of 1/R on the edge - singular case
- if( nrmX0S(ip)-std::abs(ps) < tol)
- {
- intaRm1 = 0;
- // Particle on the left of the edge
- if (ps > tol && psp1 > tol) {intaRm1 = std::log(nrmX0S(ipp)/nrmX0S(ip));}
- // Particle on the right of the edge
- else if (ps < -tol && psp1 < -tol) {intaRm1 = -std::log(nrmX0S(ipp)/nrmX0S(ip));}
- }
- // integration of 1/R on the edge - general case h>tol
- else
- {
- intaRm1 = std::asinh(psp1/ah) - std::asinh(ps/ah);
- }
- // 1/r
- Rm1 += intaRm1*(xX0S(ip)*m_nu(i,0) + yX0S(ip)*m_nu(i,1) + zX0S(ip)*m_nu(i,2));
- // grad[1/r]
- dnRm1 += intaRm1*(m_nu(i,0)*m_nrm(0) + m_nu(i,1)*m_nrm(1) + m_nu(i,2)*m_nrm(2));
- }
- return std::make_tuple(Rm1, dnRm1);
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-bool triangledata<T>::isfar(T x0, T y0, T z0) const
-{
- return ((x0-m_ctr(0))*(x0-m_ctr(0)) +
- (y0-m_ctr(1))*(y0-m_ctr(1)) +
- (z0-m_ctr(2))*(z0-m_ctr(2)) ) > 1.1*m_stp;
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-void triangledata<T>::projectP1(T x0, T y0, T z0, matrix<T>& P) const
-{
- // ((Sk-X) dot NUj) / Hj
- P(0) = ((m_B(0)-x0)*m_nu(0,0) +
- (m_B(1)-y0)*m_nu(0,1) +
- (m_B(2)-z0)*m_nu(0,2) ) / m_hgt(0);
- P(1) = ((m_C(0)-x0)*m_nu(1,0) +
- (m_C(1)-y0)*m_nu(1,1) +
- (m_C(2)-z0)*m_nu(1,2) ) / m_hgt(0);
- P(2) = ((m_A(0)-x0)*m_nu(2,0) +
- (m_A(1)-y0)*m_nu(2,1) +
- (m_A(2)-z0)*m_nu(2,2) ) / m_hgt(2);
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-inline void triangledata<T>::update(matrix<T>const& S)
-{
- update({S(0,0),S(0,1),S(0,2)},{S(1,0),S(1,1),S(1,2)},{S(2,0),S(2,1),S(2,2)});
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-inline void triangledata<T>::update(matrix<T>const& A, matrix<T>const& B, matrix<T>const& C)
-{
- // Nodes and center
- for(std::size_t j=0; j<3; ++j)
- {
- m_A(j) = A(j);
- m_B(j) = B(j);
- m_C(j) = C(j);
- m_ctr(j) = (m_A(j)+m_B(j)+m_C(j))/3.;
- }
-
- // Edges
- matrix<T> AB = B-A;
- matrix<T> BC = C-B;
- matrix<T> CA = A-C;
-
- // Normal to ABC and surface
- m_nrm(0) = AB(1)*BC(2) - AB(2)*BC(1);
- m_nrm(1) = AB(2)*BC(0) - AB(0)*BC(2);
- m_nrm(2) = AB(0)*BC(1) - AB(1)*BC(0);
- m_srf = 0.5*std::sqrt(m_nrm(0)*m_nrm(0) + m_nrm(1)*m_nrm(1) + m_nrm(2)*m_nrm(2));
- m_nrm /= 2*m_srf;
-
- // length edges
- m_lgt(2) = std::sqrt(AB(0)*AB(0)+AB(1)*AB(1)+AB(2)*AB(2));
- m_lgt(0) = std::sqrt(BC(0)*BC(0)+BC(1)*BC(1)+BC(2)*BC(2));
- m_lgt(1) = std::sqrt(CA(0)*CA(0)+CA(1)*CA(1)+CA(2)*CA(2));
- m_stp = std::max(std::max(m_lgt(0),m_lgt(1)),m_lgt(2));
-
- // tangent to edges
- for(std::size_t j=0; j<3; ++j)
- {
- m_tau(2,j) = AB(j)/m_lgt(2);
- m_tau(0,j) = BC(j)/m_lgt(0);
- m_tau(1,j) = CA(j)/m_lgt(1);
- }
-
- // normal to edges
- for(std::size_t i=0; i<3; ++i)
- {
- m_nu(i,0) = m_tau(i,1)*m_nrm(2)-m_tau(i,2)*m_nrm(1);
- m_nu(i,1) = m_tau(i,2)*m_nrm(0)-m_tau(i,0)*m_nrm(2);
- m_nu(i,2) = m_tau(i,0)*m_nrm(1)-m_tau(i,1)*m_nrm(0);
- }
-
- // Height
- m_hgt(0) = AB(0)*m_nu(0,0) + AB(1)*m_nu(0,1) + AB(2)*m_nu(0,2);
- m_hgt(1) = BC(0)*m_nu(1,0) + BC(1)*m_nu(1,1) + BC(2)*m_nu(1,2);
- m_hgt(2) = CA(0)*m_nu(2,0) + CA(1)*m_nu(2,1) + CA(2)*m_nu(2,2);
-};
-
-
-//==========================================================================//
-// FEMDATA CLASS //
-//==========================================================================//
-//==========================================================================
-// [femdata]
-///
-template<typename T>
-class femdata
-{
-public:
- femdata(matrix<std::size_t>const& elt, matrix<T>const& vtx, femtype fet, std::size_t ngs=0);
- auto dofsubdata(matrix<std::size_t>const& Idof) const;
- auto refQuadrature() const;
-
- matrix<T>const& dof() const {return m_dof;};
- matrix<T>const& gss() const {return m_gss;};
- matrix<T>const& gssnrm() const {return m_gssnrm;};
- matrix<T>const& wgt() const {return m_wgt;};
-
- matrix<std::size_t>const& elt2dof() const {return m_elt2dof;};
- matrix<std::size_t>const& elt2gss() const {return m_elt2gss;};
-
- T const& base(std::size_t i, std::size_t j) const {return m_base(i,j);};
- T const& dof(std::size_t i, std::size_t j) const {return m_dof(i,j);};
- T const& grad(std::size_t i, std::size_t j, std::size_t k) const {return m_grad[k](i,j);};
- T const& gss(std::size_t i, std::size_t j) const {return m_gss(i,j);};
- T const& gssnrm(std::size_t i, std::size_t j) const {return m_gssnrm(i,j);};
- T const& nxgrad(std::size_t i, std::size_t j, std::size_t k) const {return m_nxgrad[k](i,j);};
- T const& wgt(std::size_t i) const {return m_wgt(i);};
-
- std::size_t const& elt2dof(std::size_t i, std::size_t j) const {return m_elt2dof(i,j);};
- std::size_t const& elt2gss(std::size_t i, std::size_t j) const {return m_elt2gss(i,j);};
-
- triangledata<T>const& tridata(std::size_t i) const {return m_tridata[i];};
-
-private:
- // PRIVATE BUILDER
- void m_buildBase();
- void m_buildDof();
- void m_buildGeometry();
- void m_buildGrad();
- void m_buildQuadrature();
-
- // ATTRIBUTS
- femtype m_fet; // ENUM : TYPE OF FINITE ELEMENT
- std::size_t m_ngs; // DIMENSION : SIZE OF ELEMENTARY QUADRATURE
- matrix<T> m_dof; // COORDINATES: DEGREES OF FREEDOM
- matrix<T> m_gss; // COORDINATES: GAUSSIAN QUADRATURE
- matrix<T> m_gssnrm; // COORDINATES: NORMALS AT GAUSS POINTS
- matrix<T> m_vtx; // COORDINATES: NODES
- matrix<T> m_wgt; // VALUES : GAUSSIAN WEIGHTS
- std::vector<matrix<std::size_t>> m_dof2elt; // TABLE : DOF TO ELEMENT
- matrix<std::size_t> m_elt2dof; // TABLE : ELEMENTS TO DOF
- matrix<std::size_t> m_elt2gss; // TABLE : ELEMENTS TO GAUSSIAN QUADRATURE
- matrix<std::size_t> m_elt2vtx; // TABLE : ELEMENTS TO NODES
- matrix<T> m_base; // OPERATOR : BASIS
- std::vector<matrix<T>> m_grad; // OPERATOR : GRADIENT
- std::vector<matrix<T>> m_nxgrad; // OPERATOR : NORMALxGRADIENT
- std::vector<triangledata<T>> m_tridata; // GEOMETRY : TRIANGLE DATA
-};
-
-//==========================================================================
-//
-///
-template<typename T>
-femdata<T>::femdata(matrix<std::size_t>const& elt, matrix<T>const& vtx, femtype fet, std::size_t ngs)
-{
- // Check input
- if (size(vtx,2)!=3 || size(elt,2)!=3)
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Vertices should be a Nvtx-by-3 real matrix and elements a Nelt-by-3 integer matrix.");
- }
-
- // Input data
- m_vtx = vtx;
- m_elt2vtx = elt;
- m_ngs = ngs;
- m_fet = fet;
-
- // Computed data
- m_buildGeometry();
- m_buildQuadrature();
- m_buildDof();
- m_buildBase();
- m_buildGrad();
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-auto femdata<T>::dofsubdata(matrix<std::size_t>const& Idof) const
-{
- // Elements associated to dof indices
- std::vector<std::size_t> v;
- for (std::size_t l=0; l<numel(Idof); ++l)
- {
- v.insert(v.end(),std::begin(m_dof2elt[Idof(l)].val()),std::end(m_dof2elt[Idof(l)].val()));
- }
- matrix<std::size_t> Ielt = unique<std::size_t>(v);
-
- // Local renumerotation of dof indices
- matrix<std::size_t> e2d = eval(m_elt2dof(Ielt,col(m_elt2dof)));
- matrix<std::size_t> i1=argsort(Idof), i2=argsort(e2d);
- matrix<long> e2dl(size(e2d,1),size(e2d,2),-1);
- std::size_t k=0, l=0;
- while (l<numel(e2d))
- {
- if (e2d(i2(l))<Idof(i1(k))) {++l;}
- else if (e2d(i2(l))==Idof(i1(k)))
- {
- e2dl(i2(l)) = i1(k);
- ++l;
- }
- else if (e2d(i2(l))>Idof(i1(k)))
- {
- ++k;
- if (k==numel(Idof)){break;}
- }
- }
- return std::make_tuple(Ielt,e2dl);
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-auto femdata<T>::refQuadrature() const
-{
- matrix<T> gssloc(m_ngs,2), wgtloc(m_ngs,1);
- double a, b, w;
- switch(m_ngs)
- {
- case 0:
- break;
-
- case 1:
- gssloc(0,0) = 1/3.; gssloc(0,1) = 1/3.; wgtloc(0) = 0.5;
- break;
-
- case 3:
- a = 2/3.; b = 1/6.; w = 1/6.;
- wgtloc(0) = w; wgtloc(1) = w; wgtloc(2) = w;
- gssloc(0,0) = a; gssloc(0,1) = b;
- gssloc(1,0) = b; gssloc(1,1) = a;
- gssloc(2,0) = b; gssloc(2,1) = b;
- break;
-
- case 6:
- a = 0.816847572980459; b = 0.091576213509771; w = 0.109951743655322/2.;
- wgtloc(0) = w; wgtloc(1) = w; wgtloc(2) = w;
- gssloc(0,0) = a; gssloc(0,1) = b;
- gssloc(1,0) = b; gssloc(1,1) = a;
- gssloc(2,0) = b; gssloc(2,1) = b;
- a = 0.108103018168070; b = 0.445948490915965; w = 0.223381589678011/2.;
- wgtloc(3) = w; wgtloc(4) = w; wgtloc(5) = w;
- gssloc(3,0) = a; gssloc(3,1) = b;
- gssloc(4,0) = b; gssloc(4,1) = a;
- gssloc(5,0) = b; gssloc(5,1) = b;
- break;
-
- default:
- error(__FILE__, __LINE__, __FUNCTION__,"Unvalid quadrature size (available are 0, 1, 3 and 6).");
- break;
- }
- return std::make_tuple(gssloc,wgtloc);
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-void femdata<T>::m_buildBase()
-{
- // Initialize
- matrix<T> gssloc,wgtloc;
- std::vector<matrix<T>> fct(size(m_elt2dof,2),matrix<T>(1,m_ngs));
- m_base = zeros<T>(size(m_gss,1),size(m_elt2dof,2));
-
- // Reference quadrature
- std::tie(gssloc,wgtloc) = refQuadrature();
-
- // Basis functions
- switch(m_fet)
- {
- case vertices :
- break;
-
- case lagrangeP0:
- fct[0] = ones<T>(1,m_ngs);
- break;
-
- case lagrangeP1:
- for(std::size_t g=0; g<m_ngs; ++g)
- {
- fct[0](g) = 1 - gssloc(g,0) - gssloc(g,1);
- fct[1](g) = gssloc(g,0);
- fct[2](g) = gssloc(g,1);
- }
- break;
-
- default:
- error(__FILE__, __LINE__, __FUNCTION__,"Invalid basis function type");
- break;
- }
-
- // Assembling
- for (std::size_t e=0; e<size(m_elt2vtx,1); ++e)
- {
- for (std::size_t d=0; d<size(m_elt2dof,2); ++d)
- {
- for (std::size_t g=0; g<m_ngs; ++g)
- {
- m_base(m_elt2gss(e,g),d) = fct[d](g);
- }
- }
- }
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-void femdata<T>::m_buildDof()
-{
- // Dof, elt2dof
- switch(m_fet)
- {
- case vertices :
- break;
-
- case lagrangeP0 :
- m_dof = zeros<T>(size(m_elt2vtx,1),size(m_vtx,2));
- for (std::size_t i=0; i<size(m_elt2vtx,1); ++i)
- {
- for (std::size_t j=0; j<size(m_elt2vtx,2); ++j)
- {
- m_dof(i,0) += m_vtx(m_elt2vtx(i,j),0);
- m_dof(i,1) += m_vtx(m_elt2vtx(i,j),1);
- m_dof(i,2) += m_vtx(m_elt2vtx(i,j),2);
- }
- }
- m_dof /= size(m_elt2vtx,2);
- m_elt2dof = transpose(range(0,size(m_elt2vtx,1)));
- break;
-
- case lagrangeP1 :
- m_dof = m_vtx;
- m_elt2dof = m_elt2vtx;
- break;
-
- default:
- error(__FILE__, __LINE__, __FUNCTION__,"Invalid basis function type");
- break;
- }
-
- // Number of elements close to dof
- matrix<std::size_t> ind(1,size(m_dof,1));
- for (std::size_t l=0; l<numel(m_elt2dof); ++l)
- {
- ++ind(m_elt2dof(l));
- }
-
- // Initialize tables
- m_dof2elt.resize(size(m_dof,1));
- for (std::size_t d=0; d<size(m_dof,1); ++d)
- {
- m_dof2elt[d].resize(1,ind(d));
- ind(d)=0;
- }
-
- // Fill tables
- std::size_t d;
- for (std::size_t e=0; e<size(m_elt2dof,1); ++e)
- {
- for (std::size_t j=0; j<size(m_elt2dof,2); ++j)
- {
- d = m_elt2dof(e,j);
- m_dof2elt[d](ind(d)) = e;
- ++ind(d);
- }
- }
-
- // Check error
- if (sum(ind)!=numel(m_elt2dof))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Unavailable case.");
- }
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-void femdata<T>::m_buildGrad()
-{
- // Initialize
- matrix<T> gssloc,wgtloc;
- std::vector<matrix<T>> fctX(size(m_elt2dof,2),zeros<T>(1,m_ngs));
- std::vector<matrix<T>> fctY(size(m_elt2dof,2),zeros<T>(1,m_ngs));
- m_grad = std::vector<matrix<T>>(3, matrix<T>(size(m_gss,1),size(m_elt2dof,2)));
- m_nxgrad = std::vector<matrix<T>>(3, matrix<T>(size(m_gss,1),size(m_elt2dof,2)));
- matrix<T> AB(1,3), AC(1,3), G(2,2), Gm1(2,2), P(2,3);
- T detG;
-
- // Reference quadrature
- std::tie(gssloc,wgtloc) = refQuadrature();
-
- // Basis functions
- switch (m_fet)
- {
- case vertices :
- break;
-
- case lagrangeP0:
- break;
-
- case lagrangeP1:
- for(std::size_t g=0; g<m_ngs; ++g)
- {
- fctX[0](g) = -1; fctY[0](g) = -1;
- fctX[1](g) = 1; fctY[1](g) = 0;
- fctX[2](g) = 0; fctY[2](g) = 1;
- }
- break;
-
- default:
- error(__FILE__, __LINE__, __FUNCTION__,"Invalid basis function type");
- break;
- }
-
- // Build grad operator
- for(std::size_t e=0; e<size(m_elt2vtx,1); ++e)
- {
- // For each coordinate
- G = zeros<T>(2,2);
- for(std::size_t j=0; j<3; ++j)
- {
- // Vector basis
- AB(j) = m_vtx(m_elt2vtx(e,1),j) - m_vtx(m_elt2vtx(e,0),j);
- AC(j) = m_vtx(m_elt2vtx(e,2),j) - m_vtx(m_elt2vtx(e,0),j);
-
- // Gramm matrix (symetric)
- G(0,0) += AB(j)*AB(j);
- G(0,1) += AB(j)*AC(j);
- G(1,1) += AC(j)*AC(j);
- }
- G(1,0) = G(0,1);
-
- // inverse of Gramm matrix
- detG = G(0,0)*G(1,1) - G(1,0)*G(0,1);
- Gm1(0,0) = G(1,1)/detG; Gm1(0,1) = -G(0,1)/detG;
- Gm1(1,0) = -G(1,0)/detG; Gm1(1,1) = G(0,0)/detG;
-
- // P = (G^{-1})^t {AB;AC}
- for(std::size_t j=0; j<3; ++j)
- {
- P(0,j) = Gm1(0,0)*AB(j) + Gm1(1,0)*AC(j);
- P(1,j) = Gm1(0,1)*AB(j) + Gm1(1,1)*AC(j);
- }
-
- // grad(PHI) = P * grad(PHI_ref)
- for(std::size_t j=0; j<3; ++j)
- {
- for(std::size_t d=0; d<size(m_elt2dof,2); ++d)
- {
- for(std::size_t g=0; g<m_ngs; ++g)
- {
- m_grad[j](m_elt2gss(e,g),d) = P(0,j)*fctX[d](g) + P(1,j)*fctY[d](g);
- }
- }
- }
- }
-
- // Build nxgrad operator
- for(std::size_t k=0; k<m_grad.size(); ++k)
- {
- for(std::size_t i=0; i<size(m_grad[k],1); ++i)
- {
- for(std::size_t j=0; j<size(m_grad[k],2); ++j)
- {
- m_nxgrad[k](i,j) = (m_gssnrm(i,(k+1)%3)*m_grad[(k+2)%3](i,j) -
- m_gssnrm(i,(k+2)%3)*m_grad[(k+1)%3](i,j) );
- }
- }
- }
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-void femdata<T>::m_buildGeometry()
-{
- m_tridata = std::vector<triangledata<T>>(size(m_elt2vtx,1),triangledata<T>());
- matrix<T> A(1,3), B(1,3), C(1,3);
- for (std::size_t e=0; e<size(m_elt2vtx,1); ++e)
- {
- for (std::size_t j=0; j<3; ++j)
- {
- A(j) = m_vtx(m_elt2vtx(e,0),j);
- B(j) = m_vtx(m_elt2vtx(e,1),j);
- C(j) = m_vtx(m_elt2vtx(e,2),j);
- }
- m_tridata[e].update(A,B,C);
- }
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-void femdata<T>::m_buildQuadrature()
-{
- // Initialize data
- std::size_t ne = size(m_elt2vtx,1);
- m_gss = zeros<T>(m_ngs*ne,3);
- m_wgt = zeros<T>(m_ngs*ne,1);
- m_gssnrm = zeros<T>(m_ngs*ne,3);
- m_elt2gss = zeros<std::size_t>(ne,m_ngs);
- matrix<T> AB(1,3), AC(1,3), gssloc, wgtloc;
- std::size_t ind;
-
- // Quadrature on reference element {(0,0), (1,0), (0,1)}
- std::tie(gssloc,wgtloc) = refQuadrature();
-
- // For each element
- for(std::size_t e=0; e<ne; ++e)
- {
- // Basis vectors (AB and AC)
- for(std::size_t j=0; j<3; ++j)
- {
- AB(j) = m_tridata[e].B(j) - m_tridata[e].A(j);
- AC(j) = m_tridata[e].C(j) - m_tridata[e].A(j);
- }
-
- // Local gauss nodes: XG = A + [B-A C-A]*XGloc, weights and normals
- for(std::size_t g=0; g<m_ngs; ++g)
- {
- ind = e*m_ngs + g;
- m_elt2gss(e,g) = ind;
- m_wgt(ind) = 2*m_tridata[e].srf()*wgtloc(g);
- for(std::size_t j=0; j<3; ++j)
- {
- m_gss(ind,j) = m_tridata[e].A(j) + AB(j)*gssloc(g,0) + AC(j)*gssloc(g,1);
- m_gssnrm(ind,j) = m_tridata[e].nrm(j);
- }
- }
- }
-}
-
-
-//==========================================================================//
-// INTEGRAL EQUATION //
+// EXCITATIONS //
//==========================================================================//
//==========================================================================
//
@@ -793,33 +107,10 @@ matrix<T> rightHandSide(femdata<S>const& v, bemsource source, matrix<S>const& Y,
return B;
}
-//==========================================================================
-//
-///
-template<typename T, typename S>
-smatrix<T> mass(femdata<S>const& v)
-{
- std::vector<std::size_t> ind;
- std::vector<T> val;
- std::size_t ndof = size(v.dof(),1);
- for (std::size_t e=0; e<size(v.elt2dof(),1); ++e)
- {
- for (std::size_t g=0; g<size(v.elt2gss(),2); ++g)
- {
- std::size_t ig = v.elt2gss(e,g);
- for (std::size_t i=0; i<size(v.elt2dof(),2); ++i)
- {
- for (std::size_t j=0; j<size(v.elt2dof(),2); ++j)
- {
- ind.push_back(v.elt2dof(e,i)*ndof+v.elt2dof(e,j));
- val.push_back(v.base(ig,i) * v.wgt(ig) * v.base(ig,j));
- }
- }
- }
- }
- return smatrix<T>(ndof,ndof,ind,val);
-}
+//==========================================================================//
+// INTEGRAL OPERATORS //
+//==========================================================================//
//==========================================================================
//
///
@@ -1537,3 +828,5 @@ matrix<T> hypersingular(femdata<S>const& v, femdata<S>const& u, S wavenumber=0,
M /= (T)(4*M_PI);
return M;
}
+
+}
diff --git a/include/fem.hpp b/include/fem.hpp
index 8802bde..3cf1cbc 100644
--- a/include/fem.hpp
+++ b/include/fem.hpp
@@ -3,28 +3,30 @@
| (c) 2020 - PROPERTY OF ECOLE POLYTECHNIQUE - LGPL 3.0 |
|________________________________________________________________________|
| '&` | |
- | # | FILE : bembuilder.hpp |
+ | # | FILE : fem.hpp |
| # | VERSION : 0.1.0 |
| _#_ | AUTHOR(S) : Matthieu Aussal & Marc Bakry |
| ( # ) | CREATION : 01.04.2020 |
| / 0 \ | LAST MODIF : 31.10.2020 |
- | ( === ) | SYNOPSIS : Boundary Element Methods with Finite elements |
- | `---' | representation |
+ | ( === ) | SYNOPSIS : Quadrature and Finite Element basic tools |
+ | `---' | |
+========================================================================+
*/
#pragma once
+
+#define CASTOR_FEM_HPP
+
#include <castor/matrix.hpp>
#include <castor/smatrix.hpp>
-using namespace castor;
+namespace castor
+{
//==========================================================================//
// ENVIRONMENT DATA //
//==========================================================================//
-enum femoperator{none, gradient};
enum femtype{vertices, lagrangeP0, lagrangeP1};
-enum bemsource{PWsource, dnPWsource, SWsource, dnSWsource};
//==========================================================================//
@@ -717,82 +719,8 @@ void femdata<T>::m_buildQuadrature()
//==========================================================================//
-// INTEGRAL EQUATION //
+// VARIATIONAL FORMULATIONS //
//==========================================================================//
-//==========================================================================
-//
-///
-template<typename T>
-matrix<std::complex<T>> planeWave(matrix<T>const& X, matrix<T>const& U, T wavenumber)
-{
- if (size(X,2)!=size(U,2))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Vertices X and direction U should have same number of columns.");
- }
- if (norm(U,"2") != 1.)
- {
- warning(__FILE__, __LINE__, __FUNCTION__,"Plane wave direction is not normalized.");
- }
- std::complex<T> imk(0,wavenumber);
- return exp(imk*mtimes(X,transpose(U)));
-}
-
-//==========================================================================
-//
-///
-template<typename T>
-matrix<std::complex<T>> dnplaneWave(matrix<T>const& X, matrix<T>const& N, matrix<T>const& U, T wavenumber)
-{
- if (size(X,2)!=size(U,2) || size(N,2)!=size(U,2))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Vertices X, normals N and direction U should have same number of columns.");
- }
- if (norm(U,"2") != 1.)
- {
- warning(__FILE__, __LINE__, __FUNCTION__,"Plane wave direction is not normalized.");
- }
- std::complex<T> imk(0,wavenumber);
- return imk * mtimes(N,transpose(U)) * exp(imk*mtimes(X,transpose(U)));
-}
-
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> rightHandSide(femdata<S>const& v, bemsource source, matrix<S>const& Y, S wavenumber=0)
-{
- // Initialization
- matrix<T> B(size(v.dof(),1),1);
- matrix<T> Sxy;
- switch (source)
- {
- case PWsource:
- Sxy = planeWave<S>(v.gss(),Y,wavenumber);
- break;
- case dnPWsource:
- Sxy = dnplaneWave<S>(v.gss(),v.gssnrm(),Y,wavenumber);
- break;
- default:
- error(__FILE__, __LINE__, __FUNCTION__,"Source term not defined.");
- break;
- }
-
- // Builder
- std::size_t ig;
- for (std::size_t e=0; e<size(v.elt2dof(),1); ++e)
- {
- for (std::size_t g=0; g<size(v.elt2gss(),2); ++g)
- {
- ig = v.elt2gss(e,g);
- for (std::size_t d=0; d<size(v.elt2dof(),2); ++d)
- {
- B(v.elt2dof(e,d)) += Sxy(ig) * v.wgt(ig) * v.base(ig,d);
- }
- }
- }
- return B;
-}
-
//==========================================================================
//
///
@@ -824,716 +752,33 @@ smatrix<T> mass(femdata<S>const& v)
//
///
template<typename T, typename S>
-matrix<T> singleLayer(matrix<S>const& X, matrix<S>const& Y, S wavenumber=0,
- matrix<std::size_t> Ix={}, matrix<std::size_t> Iy={})
-{
- // Initialize dof data
- if (numel(Ix)==0) {Ix = range(0,size(X,1));}
- if (numel(Iy)==0) {Iy = range(0,size(Y,1));}
- if (max(Ix)>size(X,1) || max(Iy)>size(Y,1))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Indices are not in data bound.");
- }
- matrix<T> M(numel(Ix),numel(Iy));
-
- // Initialize local data
- std::size_t ix, iy;
- S rxy;
- T imk=std::sqrt(T(-1))*wavenumber;
-
- // Build
- for (std::size_t i=0; i<numel(Ix); ++i)
- {
- ix = Ix(i);
- for (std::size_t j=0; j<numel(Iy); ++j)
- {
- iy = Iy(j);
- rxy = std::sqrt((X(ix,0)-Y(iy,0))*(X(ix,0)-Y(iy,0)) +
- (X(ix,1)-Y(iy,1))*(X(ix,1)-Y(iy,1)) +
- (X(ix,2)-Y(iy,2))*(X(ix,2)-Y(iy,2)) );
- if (wavenumber==0)
- {
- M(i,j) = (rxy>1e-12)? (1./rxy) : (0);
- }
- else
- {
- M(i,j) = (rxy>1e-12)? (std::exp(imk*rxy)/rxy) : (imk);
- }
- }
- }
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
-}
-
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> singleLayer(matrix<S>const& X, femdata<S>const& u, S wavenumber=0,
- matrix<std::size_t> Ix={}, matrix<std::size_t> Jd={})
-{
- // Initialize dof data
- if (numel(Ix)==0) {Ix = range(0,size(X,1));}
- if (numel(Jd)==0) {Jd = range(0,size(u.dof(),1));}
- if (max(Ix)>size(X,1) || max(Jd)>size(u.dof(),1))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Indices are not in data bound.");
- }
- matrix<T> M(numel(Ix),numel(Jd));
-
- // Initialize element data
- matrix<std::size_t> Je;
- matrix<long> e2du;
- std::tie(Je,e2du) = u.dofsubdata(Jd);
-
- // Initialize local data
- std::size_t ix, jg;
- matrix<S> projectP1(1,3);
- S rxy, Rm1, dnyRm1;
- T Gxy, imk=std::sqrt(T(-1))*wavenumber;
-
- // Collocation points
- for (std::size_t i=0; i<numel(Ix); ++i)
- {
- ix = Ix(i);
- // Emmiters elements
- for (std::size_t j=0; j<length(Je); ++j)
- {
- // Analytic integration
- if (u.tridata(Je(j)).isfar(X(ix,0),X(ix,1),X(ix,2))) {Rm1 = 0;}
- else
- {
- std::tie(Rm1,dnyRm1) = u.tridata(Je(j)).intsing(X(ix,0),X(ix,1),X(ix,2));
- }
-
- // Gaussian quadrature
- for (std::size_t k=0; k<size(u.elt2gss(),2); ++k)
- {
- jg = u.elt2gss(Je(j),k);
-
- // Green kernel
- rxy = std::sqrt((X(ix,0)-u.gss(jg,0))*(X(ix,0)-u.gss(jg,0)) +
- (X(ix,1)-u.gss(jg,1))*(X(ix,1)-u.gss(jg,1)) +
- (X(ix,2)-u.gss(jg,2))*(X(ix,2)-u.gss(jg,2)));
- if (wavenumber==0)
- {
- Gxy = (rxy>1e-12)? (1./rxy) : (0);
- }
- else
- {
- Gxy = (rxy>1e-12)? (std::exp(imk*rxy)/rxy) : (imk);
- }
-
- // Integration
- Gxy *= u.wgt(jg);
-
- // Dof repartition
- for (std::size_t l=0; l<size(e2du,2); ++l)
- {
- if (e2du(j,l)!=-1)
- {
- M(i,e2du(j,l)) += Gxy * u.base(jg,l);
- }
- }
-
- // Gaussian integration
- if (Rm1!=0) {Rm1 -= ((rxy>1e-12)?(1./rxy):(0)) * u.wgt(jg);}
- }
-
- // Regularization
- if (Rm1!=0)
- {
- // Lagrange P0
- if (size(e2du,2)==1) {M(i,e2du(j)) += Rm1;}
- // Lagrange P1
- else if (size(e2du,2)==3)
- {
- u.tridata(Je(j)).projectP1(X(ix,0),X(ix,1),X(ix,2),projectP1);
- for (std::size_t l=0; l<size(e2du,2); ++l)
- {
- if (e2du(j,l)!=-1)
- {
- M(i,e2du(j,l)) += Rm1 * projectP1(l);
- }
- }
- }
- }
- }
- }
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
-}
-
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> singleLayer(femdata<S>const& v, femdata<S>const& u, S wavenumber=0,
- matrix<std::size_t> Id={}, matrix<std::size_t> Jd={})
-{
- // Initialize dof data
- if (numel(Id)==0) {Id = range(0,size(v.dof(),1));}
- if (numel(Jd)==0) {Jd = range(0,size(u.dof(),1));}
- matrix<T> M(numel(Id),numel(Jd));
-
- // Initialize element data
- matrix<std::size_t> Ie, Je;
- matrix<long> e2dv, e2du;
- std::tie(Ie,e2dv) = v.dofsubdata(Id);
- std::tie(Je,e2du) = u.dofsubdata(Jd);
-
- // Initialize local data
- std::size_t ig, jg;
- matrix<S> projectP1(1,3), xg(1,3);
- S rxy, Rm1, dnyRm1;
- T tmp, Gxy, imk=std::sqrt(T(-1))*wavenumber;
-
- // Receiver elements
- for (std::size_t i=0; i<length(Ie); ++i)
- {
- // Emitters elements
- for (std::size_t j=0; j<length(Je); ++j)
- {
- // Receiver gauss points
- for (std::size_t k=0; k<size(v.elt2gss(),2); ++k)
- {
- // Emmiters analytic integration for receiver gauss point
- ig = v.elt2gss(Ie(i),k);
- xg(0)=v.gss(ig,0); xg(1)=v.gss(ig,1); xg(2)=v.gss(ig,2);
- if (u.tridata(Je(j)).isfar(xg(0),xg(1),xg(2))) {Rm1 = 0;}
- else
- {
- std::tie(Rm1,dnyRm1) = u.tridata(Je(j)).intsing(xg(0),xg(1),xg(2));
- }
-
- // Emmiter gauss points
- for (std::size_t l=0; l<size(u.elt2gss(),2); ++l)
- {
- jg = u.elt2gss(Je(j),l);
-
- // Green kernel
- rxy = std::sqrt((xg(0)-u.gss(jg,0))*(xg(0)-u.gss(jg,0))+
- (xg(1)-u.gss(jg,1))*(xg(1)-u.gss(jg,1))+
- (xg(2)-u.gss(jg,2))*(xg(2)-u.gss(jg,2)));
- if (wavenumber==0)
- {
- Gxy = (rxy>1e-12)? (1./rxy) : (0);
- }
- else
- {
- Gxy = (rxy>1e-12)? (std::exp(imk*rxy)/rxy) : (imk);
- }
-
- // Integration
- Gxy *= v.wgt(ig) * u.wgt(jg);
-
- // Dof repartition
- for (std::size_t m=0; m<size(e2dv,2); ++m)
- {
- tmp = v.base(ig,m) * Gxy;
- for (std::size_t n=0; n<size(e2du,2); ++n)
- {
- if (e2dv(i,m)!=-1 && e2du(j,n)!=-1)
- {
- M(e2dv(i,m),e2du(j,n)) += tmp * u.base(jg,n);
- }
- }
- }
-
- // Gaussian integration
- if (Rm1!=0) {Rm1 -= ((rxy>1e-12)?(1./rxy):(0)) * u.wgt(jg);}
- }
-
- // Regularization
- if (Rm1!=0)
- {
- for (std::size_t m=0; m<size(e2dv,2); ++m)
- {
- if (e2dv(i,m)!=-1)
- {
- tmp = v.base(ig,m) * v.wgt(ig) * Rm1;
- if (size(e2du,2)==1) {M(e2dv(i,m),e2du(j)) += tmp;}
- else if (size(e2du,2)==3)
- {
- u.tridata(Je(j)).projectP1(xg(0),xg(1),xg(2),projectP1);
- for (std::size_t n=0; n<size(e2du,2); ++n)
- {
- if (e2du(j,n)!=-1)
- {
- M(e2dv(i,m),e2du(j,n)) += tmp * projectP1(n);
- }
- }
- }
- }
- }
- }
- }
- }
- }
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
-}
-
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> singleLayerInfinite(matrix<S>const& X, femdata<S>const& u, S wavenumber=0,
- matrix<std::size_t> Ix={}, matrix<std::size_t> Jd={})
-{
- // Initialize dof data
- if (numel(Ix)==0) {Ix = range(0,size(X,1));}
- if (numel(Jd)==0) {Jd = range(0,size(u.dof(),1));}
- if (max(Ix)>size(X,1) || max(Jd)>size(u.dof(),1))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Indices are not in data bound.");
- }
- matrix<T> M(numel(Ix),numel(Jd));
-
- // Initialize element data
- matrix<std::size_t> Je;
- matrix<long> e2du;
- std::tie(Je,e2du) = u.dofsubdata(Jd);
-
- // Initialize local data
- std::size_t ix, jg;
- S xdoty;
- T Gxy, imk=std::sqrt(T(-1))*wavenumber;
-
- // Collocation points
- for (std::size_t i=0; i<numel(Ix); ++i)
- {
- ix = Ix(i);
- // Emmiters elements
- for (std::size_t j=0; j<length(Je); ++j)
- {
- // Gaussian quadrature
- for (std::size_t k=0; k<size(u.elt2gss(),2); ++k)
- {
- jg = u.elt2gss(Je(j),k);
-
- // Green kernel
- xdoty = X(ix,0)*u.gss(jg,0) + X(ix,1)*u.gss(jg,1) + X(ix,2)*u.gss(jg,2);
- Gxy = std::exp(-imk*xdoty) * u.wgt(jg);
-
- // Dof repartition
- for (std::size_t l=0; l<size(e2du,2); ++l)
- {
- if (e2du(j,l)!=-1)
- {
- M(i,e2du(j,l)) += Gxy * u.base(jg,l);
- }
- }
- }
- }
- }
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
-}
-
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> doubleLayer(matrix<S>const& X, femdata<S>const& u, S wavenumber=0,
- matrix<std::size_t> Ix={}, matrix<std::size_t> Jd={})
-{
- // Initialize dof data
- if (numel(Ix)==0) {Ix = range(0,size(X,1));}
- if (numel(Jd)==0) {Jd = range(0,size(u.dof(),1));}
- if (max(Ix)>size(X,1) || max(Jd)>size(u.dof(),1))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Indices are not in data bound.");
- }
- matrix<T> M(numel(Ix),numel(Jd));
-
- // Initialize element data
- matrix<std::size_t> Je;
- matrix<long> e2du;
- std::tie(Je,e2du) = u.dofsubdata(Jd);
-
- // Initialize local data
- std::size_t ix, jg;
- matrix<S> projectP1(1,3);
- S rxy, rdotny, Rm1, dnyRm1;
- T Gxy, imk=std::sqrt(T(-1))*wavenumber;
-
- // Collocation points
- for (std::size_t i=0; i<numel(Ix); ++i)
- {
- ix = Ix(i);
- // Emmiters elements
- for (std::size_t j=0; j<length(Je); ++j)
- {
- // Analytic integration
- if (u.tridata(Je(j)).isfar(X(ix,0),X(ix,1),X(ix,2))) {dnyRm1 = 0;}
- else
- {
- std::tie(Rm1,dnyRm1) = u.tridata(Je(j)).intsing(X(ix,0),X(ix,1),X(ix,2));
- }
-
- // Gaussian quadrature
- for (std::size_t k=0; k<size(u.elt2gss(),2); ++k)
- {
- jg = u.elt2gss(Je(j),k);
-
- // Green kernel
- rxy = std::sqrt((X(ix,0)-u.gss(jg,0))*(X(ix,0)-u.gss(jg,0)) +
- (X(ix,1)-u.gss(jg,1))*(X(ix,1)-u.gss(jg,1)) +
- (X(ix,2)-u.gss(jg,2))*(X(ix,2)-u.gss(jg,2)) );
- rdotny = ((X(ix,0)-u.gss(jg,0))*u.gssnrm(jg,0) +
- (X(ix,1)-u.gss(jg,1))*u.gssnrm(jg,1) +
- (X(ix,2)-u.gss(jg,2))*u.gssnrm(jg,2) );
- if (wavenumber==0)
- {
- Gxy = (rxy>1e-12)? (rdotny/(rxy*rxy*rxy)) : (0);
- }
- else
- {
- Gxy = (rxy>1e-12)? (rdotny*(1./rxy-imk)*std::exp(imk*rxy)/(rxy*rxy)) : (0);
- }
-
- // Integration
- Gxy *= u.wgt(jg);
-
- // Dof repartition
- for (std::size_t l=0; l<size(e2du,2); ++l)
- {
- if (e2du(j,l)!=-1)
- {
- M(i,e2du(j,l)) += Gxy * u.base(jg,l);
- }
- }
-
- // Gaussian integration
- if (dnyRm1!=0) {dnyRm1 -= ((rxy>1e-12)?(rdotny/(rxy*rxy*rxy)):(0)) * u.wgt(jg);}
- }
-
- // Regularization
- if (dnyRm1!=0)
- {
- // Lagrange P0
- if (size(e2du,2)==1) {M(i,e2du(j)) += dnyRm1;}
- // Lagrange P1
- else if (size(e2du,2)==3)
- {
- u.tridata(Je(j)).projectP1(X(ix,0),X(ix,1),X(ix,2),projectP1);
- for (std::size_t l=0; l<size(e2du,2); ++l)
- {
- if (e2du(j,l)!=-1)
- {
- M(i,e2du(j,l)) += dnyRm1 * projectP1(l);
- }
- }
- }
- }
- }
- }
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
-}
-
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> doubleLayer(femdata<S>const& v, femdata<S>const& u, S wavenumber=0,
- matrix<std::size_t> Id={}, matrix<std::size_t> Jd={})
+smatrix<T> rigidity(femdata<S>const& v)
{
- // Initialize dof data
- if (numel(Id)==0) {Id = range(0,size(v.dof(),1));}
- if (numel(Jd)==0) {Jd = range(0,size(u.dof(),1));}
- matrix<T> M(numel(Id),numel(Jd));
-
- // Initialize element data
- matrix<std::size_t> Ie, Je;
- matrix<long> e2dv, e2du;
- std::tie(Ie,e2dv) = v.dofsubdata(Id);
- std::tie(Je,e2du) = u.dofsubdata(Jd);
-
- // Initialize local data
- std::size_t ig, jg;
- matrix<S> projectP1(1,3), xg(1,3);
- S rxy, rdotny, Rm1, dnyRm1;
- T tmp, Gxy, imk=std::sqrt(T(-1))*wavenumber;
-
- // Receiver elements
- for (std::size_t i=0; i<length(Ie); ++i)
+ std::vector<std::size_t> ind;
+ std::vector<T> val;
+ std::size_t ndof = size(v.dof(),1);
+ T tmp;
+ for (std::size_t e=0; e<size(v.elt2dof(),1); ++e)
{
- // Emitters elements
- for (std::size_t j=0; j<length(Je); ++j)
+ for (std::size_t g=0; g<size(v.elt2gss(),2); ++g)
{
- // Receiver gauss points
- for (std::size_t k=0; k<size(v.elt2gss(),2); ++k)
+ std::size_t ig = v.elt2gss(e,g);
+ for (std::size_t i=0; i<size(v.elt2dof(),2); ++i)
{
- // Emmiters analytic integration for receiver gauss point
- ig = v.elt2gss(Ie(i),k);
- xg(0)=v.gss(ig,0); xg(1)=v.gss(ig,1); xg(2)=v.gss(ig,2);
- if (u.tridata(Je(j)).isfar(xg(0),xg(1),xg(2))) {dnyRm1 = 0;}
- else
- {
- std::tie(Rm1,dnyRm1) = u.tridata(Je(j)).intsing(xg(0),xg(1),xg(2));
- }
-
- // Emmiter gauss points
- for (std::size_t l=0; l<size(u.elt2gss(),2); ++l)
- {
- jg = u.elt2gss(Je(j),l);
-
- // Green kernel
- rxy = std::sqrt((xg(0)-u.gss(jg,0))*(xg(0)-u.gss(jg,0)) +
- (xg(1)-u.gss(jg,1))*(xg(1)-u.gss(jg,1)) +
- (xg(2)-u.gss(jg,2))*(xg(2)-u.gss(jg,2)) );
- rdotny = ((xg(0)-u.gss(jg,0))*u.gssnrm(jg,0) +
- (xg(1)-u.gss(jg,1))*u.gssnrm(jg,1) +
- (xg(2)-u.gss(jg,2))*u.gssnrm(jg,2) );
- if (wavenumber==0)
- {
- Gxy = (rxy>1e-12)? (rdotny/(rxy*rxy*rxy)) : (0);
- }
- else
- {
- Gxy = (rxy>1e-12)? (rdotny*(1./rxy-imk)*std::exp(imk*rxy)/(rxy*rxy)) : (0);
- }
-
- // Integration
- Gxy *= v.wgt(ig) * u.wgt(jg);
-
- // Dof repartition
- for (std::size_t m=0; m<size(e2dv,2); ++m)
- {
- tmp = v.base(ig,m) * Gxy;
- for (std::size_t n=0; n<size(e2du,2); ++n)
- {
- if (e2dv(i,m)!=-1 && e2du(j,n)!=-1)
- {
- M(e2dv(i,m),e2du(j,n)) += tmp * u.base(jg,n);
- }
- }
- }
-
- // Gaussian integration
- if (dnyRm1!=0) {dnyRm1 -= ((rxy>1e-12)?(rdotny/(rxy*rxy*rxy)):(0)) * u.wgt(jg);}
- }
-
- // Regularization
- if (dnyRm1!=0)
+ for (std::size_t j=0; j<size(v.elt2dof(),2); ++j)
{
- for (std::size_t m=0; m<size(e2dv,2); ++m)
- {
- if (e2dv(i,m)!=-1)
- {
- tmp = v.base(ig,m) * v.wgt(ig) * dnyRm1;
- if (size(e2du,2)==1) {M(e2dv(i,m),e2du(j)) += tmp;}
- else if (size(e2du,2)==3)
- {
- u.tridata(Je(j)).projectP1(xg(0),xg(1),xg(2),projectP1);
- for (std::size_t n=0; n<size(e2du,2); ++n)
- {
- if (e2du(j,n)!=-1)
- {
- M(e2dv(i,m),e2du(j,n)) += tmp * projectP1(n);
- }
- }
- }
- }
- }
+ ind.push_back(v.elt2dof(e,i)*ndof+v.elt2dof(e,j));
+ tmp = v.wgt(ig) *
+ ( v.grad(ig,i,0) * v.grad(ig,j,0)
+ + v.grad(ig,i,1) * v.grad(ig,j,1)
+ + v.grad(ig,i,2) * v.grad(ig,j,2) );
+ val.push_back(tmp);
}
}
}
}
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
+ return smatrix<T>(ndof,ndof,ind,val);
}
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> doubleLayerInfinite(matrix<S>const& X, femdata<S>const& u, S wavenumber=0,
- matrix<std::size_t> Ix={}, matrix<std::size_t> Jd={})
-{
- // Initialize dof data
- if (numel(Ix)==0) {Ix = range(0,size(X,1));}
- if (numel(Jd)==0) {Jd = range(0,size(u.dof(),1));}
- if (max(Ix)>size(X,1) || max(Jd)>size(u.dof(),1))
- {
- error(__FILE__, __LINE__, __FUNCTION__,"Indices are not in data bound.");
- }
- matrix<T> M(numel(Ix),numel(Jd));
-
- // Initialize element data
- matrix<std::size_t> Je;
- matrix<long> e2du;
- std::tie(Je,e2du) = u.dofsubdata(Jd);
-
- // Initialize local data
- std::size_t ix, jg;
- S xdoty, xdotny;
- T Gxy, imk=std::sqrt(T(-1))*wavenumber;
-
- // Collocation points
- for (std::size_t i=0; i<numel(Ix); ++i)
- {
- ix = Ix(i);
- // Emmiters elements
- for (std::size_t j=0; j<length(Je); ++j)
- {
- // Gaussian quadrature
- for (std::size_t k=0; k<size(u.elt2gss(),2); ++k)
- {
- jg = u.elt2gss(Je(j),k);
-
- // Green kernel
- xdoty = X(ix,0)*u.gss(jg,0) + X(ix,1)*u.gss(jg,1) + X(ix,2)*u.gss(jg,2);
- xdotny = X(ix,0)*u.gssnrm(jg,0) + X(ix,1)*u.gssnrm(jg,1) + X(ix,2)*u.gssnrm(jg,2);
- Gxy = (-imk*xdotny) * std::exp(-imk*xdoty) * u.wgt(jg);
-
- // Dof repartition
- for (std::size_t l=0; l<size(e2du,2); ++l)
- {
- if (e2du(j,l)!=-1)
- {
- M(i,e2du(j,l)) += Gxy * u.base(jg,l);
- }
- }
- }
- }
- }
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
}
-//==========================================================================
-//
-///
-template<typename T, typename S>
-matrix<T> hypersingular(femdata<S>const& v, femdata<S>const& u, S wavenumber=0,
- matrix<std::size_t> Id={}, matrix<std::size_t> Jd={})
-{
- // Initialize dof data
- if (numel(Id)==0) {Id = range(0,size(v.dof(),1));}
- if (numel(Jd)==0) {Jd = range(0,size(u.dof(),1));}
- matrix<T> M(numel(Id),numel(Jd));
-
- // Initialize element data
- matrix<std::size_t> Ie, Je;
- matrix<long> e2dv, e2du;
- std::tie(Ie,e2dv) = v.dofsubdata(Id);
- std::tie(Je,e2du) = u.dofsubdata(Jd);
-
- // Initialize local data
- std::size_t ig, jg;
- matrix<S> projectP1(1,3), xg(1,3);
- S rxy, Rm1, dnyRm1, k2=std::pow(wavenumber,2);
- T Gxy, k2ndotnG, imk=std::sqrt(T(-1))*wavenumber;
-
- // Receiver elements
- for (std::size_t i=0; i<length(Ie); ++i)
- {
- // Emitters elements
- for (std::size_t j=0; j<length(Je); ++j)
- {
- // Receiver gauss points
- for (std::size_t k=0; k<size(v.elt2gss(),2); ++k)
- {
- // Emmiters analytic integration for receiver gauss point
- ig = v.elt2gss(Ie(i),k);
- xg(0)=v.gss(ig,0); xg(1)=v.gss(ig,1); xg(2)=v.gss(ig,2);
- if (u.tridata(Je(j)).isfar(xg(0),xg(1),xg(2))) {Rm1 = 0;}
- else
- {
- std::tie(Rm1,dnyRm1) = u.tridata(Je(j)).intsing(xg(0),xg(1),xg(2));
- }
-
- // Emmiter gauss points
- for (std::size_t l=0; l<size(u.elt2gss(),2); ++l)
- {
- jg = u.elt2gss(Je(j),l);
-
- // Green kernel
- rxy = std::sqrt((xg(0)-u.gss(jg,0))*(xg(0)-u.gss(jg,0))+
- (xg(1)-u.gss(jg,1))*(xg(1)-u.gss(jg,1))+
- (xg(2)-u.gss(jg,2))*(xg(2)-u.gss(jg,2)));
- if (wavenumber==0)
- {
- Gxy = (rxy>1e-12)? (1./rxy) : (0);
- }
- else
- {
- Gxy = (rxy>1e-12)? (std::exp(imk*rxy)/rxy) : (imk);
- }
-
- // Integration
- Gxy *= v.wgt(ig) * u.wgt(jg);
- k2ndotnG = Gxy * k2 * (v.gssnrm(ig,0)*u.gssnrm(jg,0) +
- v.gssnrm(ig,1)*u.gssnrm(jg,1) +
- v.gssnrm(ig,2)*u.gssnrm(jg,2) );
-
- // Dof repartition
- for (std::size_t m=0; m<size(e2dv,2); ++m)
- {
- for (std::size_t n=0; n<size(e2du,2); ++n)
- {
- if (e2dv(i,m)!=-1 && e2du(j,n)!=-1)
- {
- M(e2dv(i,m),e2du(j,n)) += (k2ndotnG * v.base(ig,m) * u.base(jg,n) -
- Gxy * (
- v.nxgrad(ig,m,0) * u.nxgrad(jg,n,0) +
- v.nxgrad(ig,m,1) * u.nxgrad(jg,n,1) +
- v.nxgrad(ig,m,2) * u.nxgrad(jg,n,2) ));
- }
- }
- }
-
- // Gaussian integration
- if (Rm1!=0) {Rm1 -= ((rxy>1e-12)?(1./rxy):(0)) * u.wgt(jg);}
- }
-
- // Regularization
- if (Rm1!=0)
- {
- for (std::size_t m=0; m<size(e2dv,2); ++m)
- {
- if (e2dv(i,m)!=-1)
- {
- jg = u.elt2gss(Je(j),0);
- k2ndotnG = v.base(ig,m) * v.wgt(ig) * Rm1 * k2 * (v.gssnrm(ig,0)*u.gssnrm(jg,0) +
- v.gssnrm(ig,1)*u.gssnrm(jg,1) +
- v.gssnrm(ig,2)*u.gssnrm(jg,2) );
- u.tridata(Je(j)).projectP1(xg(0),xg(1),xg(2),projectP1);
- for (std::size_t n=0; n<size(e2du,2); ++n)
- {
- if (e2du(j,n)!=-1)
- {
- M(e2dv(i,m),e2du(j,n)) += k2ndotnG * projectP1(n) - v.wgt(ig) * Rm1 *
- (- v.nxgrad(ig,m,0) * u.tridata(Je(j)).tau(n,0)
- - v.nxgrad(ig,m,1) * u.tridata(Je(j)).tau(n,1)
- - v.nxgrad(ig,m,2) * u.tridata(Je(j)).tau(n,2)) / u.tridata(Je(j)).hgt(n);
- }
- }
- }
- }
- }
- }
- }
- }
-
- // Return with 4*pi normalization
- M /= (T)(4*M_PI);
- return M;
-}
--
GitLab