Select Git revision
main.cpp 25.45 KiB
#include "AnalyticalFunction.h"
#include <gmshddm/Formulation.h>
#include <gmshddm/GmshDdm.h>
#include <pybind11/complex.h>
#include <pybind11/embed.h>
#include <pybind11/stl.h>
using namespace gmshfem;
using namespace gmshfem::field;
using namespace gmshfem::function;
using namespace gmshfem::analytics;
using namespace gmshfem::post;
using namespace gmshfem::equation;
//***************************************
// Create surface mesh
//***************************************
class BemppMesh
{
private:
int _tagSurface;
int _dim;
std::unordered_map< unsigned int, unsigned int >
_nodeConnectBemppMeshToGmshMesh;
std::vector< unsigned int > _elementsNode0;
std::vector< unsigned int > _elementsNode1;
std::vector< unsigned int > _elementsNode2;
std::vector< double > _nodesCoordX;
std::vector< double > _nodesCoordY;
std::vector< double > _nodesCoordZ;
unsigned int _renameNode(unsigned int nodeNum,
std::unordered_map< unsigned int, unsigned int >
&nodeConnectGmshMeshToBemppMesh)
{
auto it = nodeConnectGmshMeshToBemppMesh.find(nodeNum);
if(it != nodeConnectGmshMeshToBemppMesh.end()) {
return it->second;
}
else {
unsigned int size = nodeConnectGmshMeshToBemppMesh.size();
nodeConnectGmshMeshToBemppMesh.insert(
std::pair< unsigned int, unsigned int >(nodeNum, size));
_nodeConnectBemppMeshToGmshMesh.insert(
std::pair< unsigned int, unsigned int >(size, nodeNum));
return size;
}
}
void _createSurfaceMesh()
{
gmshfem::common::Timer time;
time.tick();
std::vector< int > vEntities;
gmsh::model::getEntitiesForPhysicalGroup(_dim, _tagSurface, vEntities);
std::unordered_map< unsigned int, unsigned int >
nodeConnectGmshMeshToBemppMesh;
for(unsigned int ent = 0; ent < vEntities.size(); ent++) {
std::vector< int > elementTypes;
std::vector< std::vector< std::size_t > > elementTags;
std::vector< std::vector< std::size_t > > nodeTags;
gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTags, _dim,
vEntities[ent]);
std::vector< std::size_t > nodeElem = nodeTags[0];
int numElem = elementTags[0].size();
_elementsNode0.reserve(_elementsNode0.size() + numElem);
_elementsNode1.reserve(_elementsNode1.size() + numElem);
_elementsNode2.reserve(_elementsNode2.size() + numElem);
unsigned int numNode0;
unsigned int numNode1;
unsigned int numNode2;
for(int i = 0; i < numElem; i++) { numNode0 = _renameNode(nodeElem[3 * i], nodeConnectGmshMeshToBemppMesh);
numNode1 =
_renameNode(nodeElem[3 * i + 1], nodeConnectGmshMeshToBemppMesh);
numNode2 =
_renameNode(nodeElem[3 * i + 2], nodeConnectGmshMeshToBemppMesh);
_elementsNode0.push_back(numNode0);
_elementsNode1.push_back(numNode1);
_elementsNode2.push_back(numNode2);
}
}
_nodesCoordX.reserve(nodeConnectGmshMeshToBemppMesh.size());
_nodesCoordY.reserve(nodeConnectGmshMeshToBemppMesh.size());
_nodesCoordZ.reserve(nodeConnectGmshMeshToBemppMesh.size());
for(unsigned int i = 0; i < nodeConnectGmshMeshToBemppMesh.size(); i++) {
std::vector< double > coord;
std::vector< double > parametricCoord;
int doo;
int dii;
gmsh::model::mesh::getNode(_nodeConnectBemppMeshToGmshMesh[i], coord,
parametricCoord, doo, dii);
_nodesCoordX.push_back(coord[0]);
_nodesCoordY.push_back(coord[1]);
_nodesCoordZ.push_back(coord[2]);
}
time.tock();
msg::info << "bem mesh created in " << time << msg::endl;
}
public:
BemppMesh(const int &tagSurface) :
_tagSurface(tagSurface), _dim(2)
{
_createSurfaceMesh();
}
~BemppMesh() {}
const std::vector< unsigned int > &getElementsNode0() const
{
return _elementsNode0;
}
const std::vector< unsigned int > &getElementsNode1() const
{
return _elementsNode1;
}
const std::vector< unsigned int > &getElementsNode2() const
{
return _elementsNode2;
}
const std::unordered_map< unsigned int, unsigned int > &
getnodeConnectBemppMeshToGmshMesh() const
{
return _nodeConnectBemppMeshToGmshMesh;
}
const std::vector< double > &getnodesCoordX() const
{
return _nodesCoordX;
}
const std::vector< double > &getnodesCoordY() const
{
return _nodesCoordY;
}
const std::vector< double > &getnodesCoordZ() const
{
return _nodesCoordZ;
}
};
//***************************************
// BemFormulation
//***************************************
#pragma GCC visibility push(hidden)class BemFormulation
: public gmshfem::problem::Formulation< std::complex< double > >
{
private:
bool _firstIt;
pybind11::module _bem;
gmshfem::field::FieldInterface< std::complex< double > > *_E;
gmshfem::field::FieldInterface< std::complex< double > > *_H;
BemppMesh _surfaceMesh;
std::vector< int > _typeKeys;
std::vector< unsigned long long > _entityKeys;
std::complex< double > _Zi;
std::complex< double > _Ze;
double _ki;
double _ke;
Field< std::complex< double >, Form::Form2 > *_interfaceFieldG;
std::vector< double > _ChangeOfBasisVectorBemToGmsh;
double *_bembemL2errorM;
double *_bembemL2errorJ;
double _opT;
public:
BemFormulation(const std::string &name,
gmshfem::field::FieldInterface< std::complex< double > > *E,
gmshfem::field::FieldInterface< std::complex< double > > *H,
const BemppMesh &surfaceMesh, const std::complex< double > &Zi,
const std::complex< double > &Ze, const double &ki,
const double &ke, Field< std::complex< double >, Form::Form2 > *g,
double *bembemL2errorM, double *bembemL2errorJ,
const double &opT) :
gmshfem::problem::Formulation< std::complex< double > >(name),
_firstIt(true), _E(E), _H(H), _surfaceMesh(surfaceMesh), _Zi(Zi),
_Ze(Ze), _ki(ki), _ke(ke), _interfaceFieldG(g),
_ChangeOfBasisVectorBemToGmsh(), _bembemL2errorM(bembemL2errorM),
_bembemL2errorJ(bembemL2errorJ), _opT(opT)
{
// load our bem.py module
_bem = pybind11::module::import("bem");
}
gmshfem::common::Timer
pre(const gmshfem::problem::DofsSort::Algorithm algo =
gmshfem::common::Options::instance()->dofsSortAlgorithm)
{
msg::info << "Bem pre-process" << _name << "..." << msg::endl;
gmshfem::common::Timer time;
time.tick();
if(_firstIt) {
_bem.attr("pre")(
_ke, _surfaceMesh.getElementsNode0(), _surfaceMesh.getElementsNode1(),
_surfaceMesh.getElementsNode2(), _surfaceMesh.getnodesCoordX(),
_surfaceMesh.getnodesCoordY(), _surfaceMesh.getnodesCoordZ(), _opT,
_Ze, _surfaceMesh.getnodeConnectBemppMeshToGmshMesh());
pybind11::object edgeNumFollowingBemDofOrderPython =
_bem.attr("getEdgeNumInBemDofOrder")();
_entityKeys = pybind11::cast< std::vector< unsigned long long > >(
edgeNumFollowingBemDofOrderPython);
_typeKeys.resize(_entityKeys.size(), 1);
pybind11::object ChangeOfBasisVectorPython =
_bem.attr("getChangeOfBasisVector")();
_ChangeOfBasisVectorBemToGmsh =
pybind11::cast< std::vector< double > >(ChangeOfBasisVectorPython);
_bem.attr("clearGlobalVariable")();
for(unsigned int i = 0; i < _entityKeys.size(); ++i) {
dofs::UnknownDof *dof = new dofs::UnknownDof(
_typeKeys[i] + GMSHFEM_DOF_FIELD_OFFSET * _E->tag(),
_entityKeys[i]); dof->numDof(i + 1);
_E->setDof(dof);
}
for(unsigned int i = 0; i < _entityKeys.size(); ++i) {
dofs::UnknownDof *dof = new dofs::UnknownDof(
_typeKeys[i] + GMSHFEM_DOF_FIELD_OFFSET * _H->tag(),
_entityKeys[i]);
dof->numDof(i + 1);
_H->setDof(dof);
}
_firstIt = false;
time.tock();
msg::info << _entityKeys.size() << " unknown dofs created in " << time
<< "s:" << msg::endl;
}
else {
time.tock();
msg::info << "Bem pre-process" << _name << "...Done" << msg::endl;
}
return time;
}
gmshfem::common::Timer solve(const bool reusePreconditioner = false)
{
msg::info << "bem solve" << msg::endl;
gmshfem::common::Timer time;
time.tick();
unsigned int nDof = _entityKeys.size();
std::vector< std::complex< double > > g_gmshfem(nDof);
std::vector< std::complex< double > > g_bem;
g_bem.resize(nDof);
_interfaceFieldG->getValues(_typeKeys, _entityKeys, g_gmshfem, 0, nDof);
// low order basis Hdiv in gmshfem (ei low order hcurl basis): n^ei
// low order basis Hdiv in bempp (ei low order hcurl basis): ei^n
for(unsigned int i = 0; i < nDof; i++) {
g_bem[i] = -(g_gmshfem[i]) * _ChangeOfBasisVectorBemToGmsh[i];
}
bool physicalSource;
bool artificialSource;
getAttribute("ddm::physicalCommutator", physicalSource);
getAttribute("ddm::artificialCommutator", artificialSource);
pybind11::object pythonVect =
_bem.attr("solveIE")(g_bem, physicalSource, artificialSource, _opT);
std::vector< std::complex< double > > EH_bem =
pybind11::cast< std::vector< std::complex< double > > >(pythonVect);
std::vector< std::complex< double > > E_gmshfem(nDof);
std::vector< std::complex< double > > H_gmshfem(nDof);
for(unsigned int i = 0; i < nDof; i++) {
E_gmshfem[i] = EH_bem[i] / _ChangeOfBasisVectorBemToGmsh[i];
}
for(unsigned int i = 0; i < nDof; i++) {
H_gmshfem[i] = EH_bem[i + nDof] / _ChangeOfBasisVectorBemToGmsh[i];
}
_E->assignValues(E_gmshfem);
_H->assignValues(H_gmshfem);
if(artificialSource && physicalSource) {
pybind11::object errorPython =
_bem.attr("getL2errorJM")(g_bem, _Zi, _ki, _ke, _Ze, _opT);
std::vector< double > error(2);
error = pybind11::cast< std::vector< double > >(errorPython);
*_bembemL2errorM = error[0];
*_bembemL2errorJ = error[1];
}
time.tock();
msg::info << "bem solve in " << time << "s" << msg::endl;
return time;
}
};
#pragma GCC visibility pop
//***************************************
// Geometry
//***************************************
void meshSphere(double const &lc, const int &tagSurface, const int &tagVolume,
const double &r)
{
gmsh::model::add("sphere");
gmsh::model::occ::addPoint(0.5, 0, 0, 1);
gmsh::model::occ::addSphere(0, 0, 0, r, 1);
gmsh::model::occ::synchronize();
gmsh::model::addPhysicalGroup(3, {1}, tagVolume);
gmsh::model::setPhysicalName(3, 1, "volume");
gmsh::model::addPhysicalGroup(2, {1}, tagSurface);
gmsh::model::setPhysicalName(2, tagSurface, "surface");
gmsh::option::setNumber("Mesh.CharacteristicLengthMin", lc);
gmsh::option::setNumber("Mesh.CharacteristicLengthMax", lc);
gmsh::model::mesh::generate(3);
}
int main(int argc, char **argv)
{
gmshddm::common::GmshDdm gmshDdm(argc, argv);
/****
* start python interpreter
****/
pybind11::scoped_interpreter guard{};
/****
* input Data
****/
double pi = 3.14159265359;
double epsrext = 1;
double murext = 1;
double epsrint = 4;
double murint = 1;
double mu0 = 4 * pi * 1e-7;
double eps0 = 8.854187817e-12;
double muext = murext * mu0;
double epsext = eps0 * epsrext;
double muint = murint * mu0;
double epsint = eps0 * epsrint;
double frequency =
4 / (2. * pi * std::sqrt(epsint * muint));
double w = 2 * pi * frequency;
std::complex< double > im(0., 1.);
std::complex< double > k1e = im * w * epsext;
std::complex< double > k2e = im * w * muext;
double ke = std::sqrt(real(-k1e * k2e));
std::complex< double > k1i = im * w * epsint;
std::complex< double > k2i = im * w * muint;
double ki = std::sqrt(real(-k1i * k2i));
double Zi = std::sqrt(muint / epsint);
double Ze = std::sqrt(muext / epsext);
int tagSurface = 1;
int tagVolume = 3;
double r = 1;
int FEorderAlpha = 2;
int FEorderBeta = 3;
double GMRESTol = 4;
double pointsByWl = 10;
int pade = 1;
double angleDivisor = 2.;
int padeOrder = 2;
gmshDdm.userDefinedParameter(FEorderAlpha, "FEorderAlpha");
gmshDdm.userDefinedParameter(FEorderBeta, "FEorderBeta");
gmshDdm.userDefinedParameter(pointsByWl, "pointsByWl");
gmshDdm.userDefinedParameter(pade, "pade");
double lc = 2 * pi / (pointsByWl * std::max(ki, ke));
std::string gauss = "Gauss" + std::to_string(2 * (std::max(FEorderAlpha, FEorderBeta) + 1));
gmshDdm.userDefinedParameter(gauss, "gauss");
std::vector< unsigned int > degree = {0, unsigned(FEorderBeta)};
double mL2Error = 0;
double jL2Error = 0;
double opT = Ze;
double Tp = 1 / opT;
double Tm = 1 / Ze;
/****
* Generate mesh
****/
meshSphere(lc, tagSurface, tagVolume, r);
/****
* create surface mesh
****/
BemppMesh surfaceMesh(tagSurface);
/****
* ddm part
******/
msg::info << "********************************" << msg::endl;
msg::info << " Wave properties:" << msg::endl;
msg::info << " - kint = " << ki << " [1/m]." << msg::endl;
msg::info << " - kout = " << ke << " [1/m]." << msg::endl;
msg::info << " Problem properties:" << msg::endl;
msg::info << " - Ze : " << Ze << msg::endl;
msg::info << " - Zi : " << Zi << msg::endl;
msg::info << " - FEorderAlpha = " << FEorderAlpha << "" << msg::endl;
msg::info << " - FEorderBeta = " << FEorderBeta << "" << msg::endl;
msg::info << " Mesh properties:" << msg::endl;
msg::info << " - Approximate number of points by wavelength = " << pointsByWl
<< "" << msg::endl;
msg::info << " GMRES tolerance:" << msg::endl;
msg::info << " - " << GMRESTol << msg::endl;
msg::info << " Use Pade:" << msg::endl;
msg::info << " - " << pade << msg::endl;
msg::info << "********************************" << msg::endl;
// create domains
gmshddm::domain::Subdomain omega(2);
gmshddm::domain::Interface sigma(2);
omega(0) = gmshfem::domain::Domain(3, tagVolume);
omega(1) = gmshfem::domain::Domain(2, tagSurface);
sigma(0, 1) = gmshfem::domain::Domain(2, tagSurface);
sigma(1, 0) = gmshfem::domain::Domain(2, tagSurface);
// Create fields
gmshddm::field::SubdomainField< std::complex< double >,
gmshfem::field::Form::Form1 >
E("E", omega | sigma,
gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl, FEorderAlpha);
Field< std::complex< double >, Form::Form1 > Hbem(
"Hbem", omega(1), FunctionSpaceTypeForm1::HierarchicalHCurl, 0);
gmshddm::field::InterfaceField< std::complex< double >,
gmshfem::field::Form::Form2 >
g("g", sigma, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl,
degree);
std::vector< gmshfem::problem::Formulation< std::complex< double > > * > vol(2);
vol[0] = new gmshfem::problem::Formulation< std::complex< double > >("int");
vol[1] = new BemFormulation("ext", &E(1), &Hbem, surfaceMesh, Zi, Ze, ki, ke,
&g(0, 1), &mL2Error, &jL2Error, opT);
// Define interface formulations
std::vector< std::unordered_map<
unsigned int, gmshfem::problem::Formulation< std::complex< double > > * > >
sur(2);
sur[0][1] =
new gmshfem::problem::Formulation< std::complex< double > >("interface_01");
sur[1][0] =
new gmshfem::problem::Formulation< std::complex< double > >("interface_10");
// Define ddm formulation
gmshddm::problem::Formulation< std::complex< double > > formulation(
"weak coupling", vol, sur); VectorFunction< std::complex< double > > n = normal< std::complex< double > >();
// Tell to the formulation that g is field that have to be exchanged
// between subdomains
formulation.addInterfaceField(g);
// interior problem formulation
formulation(0).integral(curl(dof(E(0))), curl(tf(E(0))), omega(0), gauss);
formulation(0).integral(-ki * ki * dof(E(0)), tf(E(0)), omega(0), gauss);
// Boundary conditions
formulation(0).integral(-im * ki * Zi * formulation.artificialSource(g(1, 0)),
tf(E(0)), sigma(0, 1), gauss);
formulation(0, 1).integral(dof(g(0, 1)), tf(g(0, 1)), sigma(0, 1), gauss);
formulation(0, 1).integral(-Tp * n % (E(0) % n), tf(g(0, 1)), sigma(0, 1),
gauss);
formulation(0, 1).integral(-formulation.artificialSource(g(1, 0)),
tf(g(0, 1)), sigma(0, 1), gauss);
formulation(1, 0).integral(dof(g(1, 0)), tf(g(1, 0)), sigma(1, 0), gauss);
formulation(1, 0).integral(Tp * E(1), tf(g(1, 0)), sigma(1, 0), gauss);
formulation(1, 0).integral(-formulation.artificialSource(g(0, 1)), tf(g(1, 0)), sigma(1, 0), gauss);
std::vector< FieldInterface< std::complex< double > > * > fieldBucket;
if(pade == 0) {
formulation(0).integral(-im * ki * Zi * Tm * dof(E(0)), tf(E(0)),
sigma(0, 1), gauss);
formulation(1, 0).integral(Tm * E(1), tf(g(1, 0)), sigma(1, 0), gauss);
formulation(0, 1).integral(-Tm * n % (E(0) % n), tf(g(0, 1)), sigma(0, 1),
gauss);
}
else {
// all what we need for pade approximations
const double angle = pi / angleDivisor;
const double M = 2. * padeOrder + 1.;
std::complex< double > exp1 =
std::complex< double >(std::cos(-angle / 2), std::sin(-angle / 2));
std::complex< double > expC =
std::complex< double >(std::cos(angle / 2.), std::sin(angle / 2.));
std::complex< double > exp3 =
std::complex< double >(std::cos(-angle), std::sin(-angle));
std::complex< double > C0 = expC;
std::vector< std::complex< double > > Aj(padeOrder, 0.);
std::vector< std::complex< double > > Bj(padeOrder, 0.);
std::complex< double > kepsilonsquare = std::complex< double >(
ke, 0.4 * pow(ke, 1. / 3.) * pow((1 / r), 2. / 3.));
kepsilonsquare *= kepsilonsquare;
for(int i = 0; i < padeOrder; i++) {
std::complex< double > aj = std::sin((i + 1) * pi / M);
aj = aj * aj * 2. / M;
std::complex< double > bj = std::cos((i + 1) * pi / M);
bj *= bj;
std::complex< double > den = 1. + bj * (exp3 - 1.);
Aj[i] = exp1 * aj / (den * den);
Bj[i] = exp3 * bj / den;
C0 += expC * aj * (exp3 - 1.) / den;
}
// Declare Auxiliary fields
Field< std::complex< double >, Form::Form1 > *E_auxi0 =
new Field< std::complex< double >, Form::Form1 >(
"E_auxi0", sigma(0, 1),
gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
FEorderBeta);
Field< std::complex< double >, Form::Form1 > *E_auxi1 =
new Field< std::complex< double >, Form::Form1 >(
"E_auxi1", sigma(1, 0),
gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
0);
fieldBucket.push_back(E_auxi0);
fieldBucket.push_back(E_auxi1);
std::vector< Field< std::complex< double >, Form::Form0 > * > zeta01;
std::vector< Field< std::complex< double >, Form::Form1 > * > phi01;
std::vector< Field< std::complex< double >, Form::Form1 > * > phi10;
std::vector< Field< std::complex< double >, Form::Form0 > * > zeta10;
for(int i = 0; i < padeOrder; ++i) { phi01.push_back(new Field< std::complex< double >, Form::Form1 >(
"phi01_" + std::to_string(i), sigma(0, 1),
gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
0));
zeta01.push_back(new Field< std::complex< double >, Form::Form0 >(
"zeta01_" + std::to_string(i), sigma(0, 1),
gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, 1));
phi10.push_back(new Field< std::complex< double >, Form::Form1 >(
"phi10_" + std::to_string(i), sigma(1, 0),
gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
0));
zeta10.push_back(new Field< std::complex< double >, Form::Form0 >(
"zeta10_" + std::to_string(i), sigma(1, 0),
gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, 1));
fieldBucket.push_back(zeta10.back());
fieldBucket.push_back(phi10.back());
fieldBucket.push_back(zeta01.back());
fieldBucket.push_back(phi01.back());
}
// Boundary conditions
formulation(0).integral(im * ki * Zi / Ze * dof(*E_auxi0), tf(E(0)),
sigma(0, 1), gauss);
// Auxiliary equations E_auxi0
formulation(0).integral(C0 * dof(*E_auxi0), tf(*E_auxi0), sigma(0, 1),
gauss);
formulation(0).integral(dof(E(0)), tf(*E_auxi0), sigma(0, 1), gauss);
formulation(0).integral(-1. / kepsilonsquare * curl(dof(E(0))),
curl(tf(*E_auxi0)), sigma(0, 1), gauss);
for(int i = 0; i < padeOrder; ++i) {
formulation(0).integral(Aj[i] * grad(dof(*zeta01[i])), tf(*E_auxi0),
sigma(0, 1), gauss);
formulation(0).integral(-Aj[i] / kepsilonsquare * curl(dof(*phi01[i])),
curl(tf(*E_auxi0)), sigma(0, 1), gauss);
// Auxiliary equations phi01
formulation(0).integral(dof(*phi01[i]), tf(*phi01[i]), sigma(0, 1),
gauss);
formulation(0).integral(Bj[i] * grad(dof(*zeta01[i])), tf(*phi01[i]),
sigma(0, 1), gauss);
formulation(0).integral(-Bj[i] / kepsilonsquare * curl(dof(*phi01[i])),
curl(tf(*phi01[i])), sigma(0, 1), gauss);
formulation(0).integral(-dof(*E_auxi0), tf(*phi01[i]), sigma(0, 1),
gauss);
// Auxiliary equations zeta01
formulation(0).integral(dof(*zeta01[i]), tf(*zeta01[i]), sigma(0, 1),
gauss);
formulation(0).integral(1. / kepsilonsquare * dof(*phi01[i]),
grad(tf(*zeta01[i])), sigma(0, 1), gauss);
}
formulation(0, 1).integral(1. / Ze * (*E_auxi0), tf(g(0, 1)), sigma(0, 1),
gauss);
formulation(1, 0).integral(-1. / Ze * dof(*E_auxi1), tf(g(1, 0)),
sigma(1, 0), gauss);
// Auxiliary equations E_auxi1
formulation(1, 0).integral(C0 * dof(*E_auxi1), tf(*E_auxi1), sigma(1, 0),
gauss);
formulation(1, 0).integral(-(E(1) % n) % n, tf(*E_auxi1), sigma(1, 0),
gauss);
formulation(1, 0).integral(-1. / kepsilonsquare * curl(E(1)),
curl(tf(*E_auxi1)), sigma(1, 0), gauss);
for(int i = 0; i < padeOrder; ++i) {
formulation(1, 0).integral(Aj[i] * grad(dof(*zeta10[i])), tf(*E_auxi1),
sigma(1, 0), gauss);
formulation(1, 0).integral(-Aj[i] / kepsilonsquare * curl(dof(*phi10[i])),
curl(tf(*E_auxi1)), sigma(1, 0), gauss);
// Auxiliary equations phi10 formulation(1, 0).integral(dof(*phi10[i]), tf(*phi10[i]), sigma(1, 0),
gauss);
formulation(1, 0).integral(Bj[i] * grad(dof(*zeta10[i])), tf(*phi10[i]),
sigma(1, 0), gauss);
formulation(1, 0).integral(-Bj[i] / kepsilonsquare * curl(dof(*phi10[i])),
curl(tf(*phi10[i])), sigma(1, 0), gauss);
formulation(1, 0).integral(-dof(*E_auxi1), tf(*phi10[i]), sigma(1, 0),
gauss);
// Auxiliary equations zeta10
formulation(1, 0).integral(dof(*zeta10[i]), tf(*zeta10[i]), sigma(1, 0),
gauss);
formulation(1, 0).integral(1. / kepsilonsquare * dof(*phi10[i]),
grad(tf(*zeta10[i])), sigma(1, 0), gauss);
}
}
// Solve the DDM formulation
formulation.pre();
formulation.solve("gmres", pow(10, -GMRESTol), 100, true);
for(unsigned int i = 0; i < fieldBucket.size(); ++i) {
delete fieldBucket[i];
}
msg::info << "***************************************************************************************" <<msg::endl;
msg::info << "Differences in L2-norm between the strong BEM-BEM coupling and the weak FEM-BEM coupling"<<msg::endl;
msg::info << "***************************************************************************************" <<msg::endl;
msg::info <<"e_bembem_M: " << mL2Error << msg::endl;
msg::info << "e_bembem_J: " << jL2Error << msg::endl;
msg::info << "***************************************************************************************" <<msg::endl;
msg::info << "Errors in L2-norm between the analytial solution and the weak FEM-BEM coupling"<<msg::endl;
msg::info << "***************************************************************************************" <<msg::endl;
AnalyticalFunction< maxwell3D::MMagneticSurfaceCurrent< std::complex< double > > >
Mexact(ki, ke, k2i, k2e, r);
msg::info << "e_analy_M: "
<< real(sqrt(
integrate(pow(abs(xComp(Mexact) - xComp(E(1) % n)), 2) +
pow(abs(yComp(Mexact) - yComp(E(1) % n)), 2) +
pow(abs(zComp(Mexact) - zComp(E(1) % n)), 2),
omega(1), gauss) /
integrate(pow(abs(xComp(Mexact)), 2) +
pow(abs(yComp(Mexact)), 2) +
pow(abs(zComp(Mexact)), 2),
omega(1), gauss)))
<< msg::endl;
AnalyticalFunction<
maxwell3D::JElectricSurfaceCurrent< std::complex< double > > >
Jexact(ki, ke, k2i, k2e, r);
msg::info << "e_analy_J: "<<real(
sqrt(integrate(pow(abs(xComp(Jexact) - xComp(Hbem % n)), 2) +
pow(abs(yComp(Jexact) - yComp(Hbem % n)), 2) +
pow(abs(zComp(Jexact) - zComp(Hbem % n)), 2),
omega(1), gauss) /
integrate(pow(abs(xComp(Jexact)), 2) +
pow(abs(yComp(Jexact)), 2) +
pow(abs(zComp(Jexact)), 2),
omega(1), gauss)))<< msg::endl;;
return 0;
}