Select Git revision
-
Christophe Geuzaine authored- removed old mesh-based Curvature calculation class (using RBFs) - removed multiscale laplace partitioner
Christophe Geuzaine authored- removed old mesh-based Curvature calculation class (using RBFs) - removed multiscale laplace partitioner
Subproblem3D.cpp 57.43 KiB
#include "Subproblem3D.h"
using gmshfem::equation::dof;
using gmshfem::equation::tf;
using gmshfem::function::operator-;
using gmshfem::function::operator+;
using gmshfem::function::operator*;
namespace D3 {
// **********************************
// Boundary
// **********************************
Boundary::Boundary(const unsigned int boundary) : _boundary(boundary)
{
}
Boundary::~Boundary()
{
}
std::string Boundary::orientation() const
{
switch (_boundary) {
case 0:
return "E";
break;
case 1:
return "N";
break;
case 2:
return "W";
break;
case 3:
return "S";
break;
case 4:
return "D";
break;
case 5:
return "U";
break;
default:
break;
}
return "null";
}
// **********************************
// Edge
// **********************************
Edge::Edge(const unsigned int edge, Boundary *first, Boundary *second) : _edge(edge), _bnd {first, second}
{
}
Edge::~Edge()
{
}
std::string Edge::orientation() const
{
switch (_edge) {
case 0:
return "ED";
break;
case 1:
return "ND"; break;
case 2:
return "WD";
break;
case 3:
return "SD";
break;
case 4:
return "EN";
break;
case 5:
return "NW";
break;
case 6:
return "WS";
break;
case 7:
return "SE";
break;
case 8:
return "EU";
break;
case 9:
return "NU";
break;
case 10:
return "WU";
break;
case 11:
return "SU";
break;
default:
break;
}
return "null";
}
// **********************************
// Corner
// **********************************
Corner::Corner(const unsigned int corner, Edge *first, Edge *second, Edge *third) : _corner(corner), _bnd {first, second, third}
{
}
Corner::~Corner()
{
}
std::string Corner::orientation() const
{
switch (_corner) {
case 0:
return "END";
break;
case 1:
return "NWD";
break;
case 2:
return "WSD";
break;
case 3:
return "SED";
break;
case 4:
return "ENU";
break; case 5:
return "NWU";
break;
case 6:
return "WSU";
break;
case 7:
return "SEU";
break;
default:
break;
}
return "null";
}
// **********************************
// Subproblem
// **********************************
void Subproblem::_parseParameters(std::string &method, std::vector< std::string > ¶meters, const std::string &str) const
{
method.clear();
parameters.clear();
std::string *current = &method;
for(unsigned int i = 0; i < str.size(); ++i) {
if(str[i] == '_') {
parameters.push_back(std::string());
current = ¶meters.back();
}
else {
(*current) += str[i];
}
}
}
Subproblem::Subproblem(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters, const std::string &E, const std::string &N, const std::string &W, const std::string &S, const std::string &D, const std::string &U) : _domains(domains), _parameters(parameters), _formulation(formulation), _fieldName(fieldName), _boundary(6, nullptr), _edge(12, nullptr), _corner(8, nullptr)
{
std::string type[6] = {E, N, W, S, D, U};
std::string method[6];
std::vector< std::string > bndParameters[6];
for(unsigned int i = 0; i < 6; ++i) {
_parseParameters(method[i], bndParameters[i], type[i]);
}
for(unsigned int b = 0; b < 6; ++b) {
if(method[b] == "sommerfeld") { // str: sommerfeld
_boundary[b] = new Sommerfeld(b);
}
else if(method[b] == "habc") { // str: habc_[N]_[theta]
const unsigned int N = std::stoi(bndParameters[b][0]);
const double theta = std::stof(bndParameters[b][1]);
_boundary[b] = new HABC(b, N, theta);
}
else if(method[b] == "pmlContinuous") { // str: pmlContinuous_[N_pml]_[Type]
const double N_pml = std::stof(bndParameters[b][0]);
const std::string type = bndParameters[b][1];
_boundary[b] = new PmlContinuous(b, N_pml, type);
}
// marche pas, il faudrait un vrai H(div)
else if(method[b] == "pmlDiscontinuous") { // str: pmlDiscontinuous_[N_pml]_[Type]
const double N_pml = std::stof(bndParameters[b][0]);
const std::string type = bndParameters[b][1];
_boundary[b] = new PmlDiscontinuous(b, N_pml, type);
}
else {
gmshfem::msg::error << "Unknown method '" << method[b] << "'" << gmshfem::msg::endl;
}
}
for(unsigned int e = 0; e < 12; ++e) { if(e < 4) {
if(method[e%4] == "habc" && method[4] == "habc") {
_edge[e] = new HABC_HABC(e, static_cast< HABC * >(_boundary[e%4]), static_cast< HABC * >(_boundary[4]));
}
else if(method[e%4] == "pmlContinuous" && method[4] == "pmlContinuous") {
_edge[e] = new PmlContinuous_PmlContinuous(e, static_cast< PmlContinuous * >(_boundary[e%4]), static_cast< PmlContinuous * >(_boundary[4]));
}
}
else if(e >= 4 && e < 8) {
if(method[e%4] == "habc" && method[(e+1)%4] == "habc") {
_edge[e] = new HABC_HABC(e, static_cast< HABC * >(_boundary[e%4]), static_cast< HABC * >(_boundary[(e+1)%4]));
}
else if(method[e%4] == "pmlContinuous" && method[(e+1)%4] == "pmlContinuous") {
_edge[e] = new PmlContinuous_PmlContinuous(e, static_cast< PmlContinuous * >(_boundary[e%4]), static_cast< PmlContinuous * >(_boundary[(e+1)%4]));
}
}
else {
if(method[e%4] == "habc" && method[5] == "habc") {
_edge[e] = new HABC_HABC(e, static_cast< HABC * >(_boundary[e%4]), static_cast< HABC * >(_boundary[5]));
}
else if(method[e%4] == "pmlContinuous" && method[5] == "pmlContinuous") {
_edge[e] = new PmlContinuous_PmlContinuous(e, static_cast< PmlContinuous * >(_boundary[e%4]), static_cast< PmlContinuous * >(_boundary[5]));
}
}
}
for(unsigned int c = 0; c < 8; ++c) {
if(c < 4) {
if(method[c%4] == "habc" && method[(c+1)%4] == "habc" && method[4] == "habc") {
_corner[c] = new HABC_HABC_HABC(c, static_cast< HABC_HABC * >(_edge[c%4]), static_cast< HABC_HABC * >(_edge[(c+1)%4]), static_cast< HABC_HABC * >(_edge[c+4]));
}
else if(method[c%4] == "pmlContinuous" && method[(c+1)%4] == "pmlContinuous" && method[4] == "pmlContinuous") {
_corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c+4]));
}
}
else {
if(method[c%4] == "habc" && method[(c+1)%4] == "habc" && method[5] == "habc") {
_corner[c] = new HABC_HABC_HABC(c, static_cast< HABC_HABC * >(_edge[c%4+8]), static_cast< HABC_HABC * >(_edge[(c+1)%4+8]), static_cast< HABC_HABC * >(_edge[c]));
}
else if(method[c%4] == "pmlContinuous" && method[(c+1)%4] == "pmlContinuous" && method[5] == "pmlContinuous") {
_corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c]));
}
}
}
}
Subproblem::~Subproblem()
{
for(unsigned int c = 0; c < 8; ++c) {
if(_corner[c] != nullptr) {
delete _corner[c];
}
}
for(unsigned int e = 0; e < 12; ++e) {
if(_edge[e] != nullptr) {
delete _edge[e];
}
}
for(unsigned int b = 0; b < 6; ++b) {
if(_boundary[b] != nullptr) {
delete _boundary[b];
}
}
}
void Subproblem::writeFormulation()
{
for(unsigned int b = 0; b < 6; ++b) {
_boundary[b]->writeFormulation(_formulation, _fieldName, _domains, _parameters);
}
for(unsigned int e = 0; e < 12; ++e) {
if(_edge[e] != nullptr) {
_edge[e]->writeFormulation(_formulation, _domains, _parameters);
}
}
for(unsigned int c = 0; c < 8; ++c) {
if(_corner[c] != nullptr) {
_corner[c]->writeFormulation(_formulation, _domains, _parameters);
}
}
}
Boundary *Subproblem::getBoundary(const unsigned int b) const
{
return _boundary[b];
}
Edge *Subproblem::getEdge(const unsigned int e) const
{
return _edge[e];
}
Corner *Subproblem::getCorner(const unsigned int c) const
{
return _corner[c];
}
// **********************************
// Sommerfeld
// **********************************
Sommerfeld::Sommerfeld(const unsigned int boundary) : Boundary(boundary)
{
}
Sommerfeld::~Sommerfeld()
{
delete _v;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* Sommerfeld::getV() const
{
return _v;
}
void Sommerfeld::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
_v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
const std::complex< double > im(0., 1.);
formulation.integral(-dof(*_v), tf(*u), sigma, gauss);
formulation.integral(dof(*_v), tf(*_v), sigma, gauss);
formulation.integral(im * kappa * dof(*u), tf(*_v), sigma, gauss);
}
// **********************************
// HABC
// **********************************
HABC::HABC(const unsigned int boundary, const unsigned int N, const double theta) : Boundary(boundary), _N(N), _theta(theta)
{
}
HABC::~HABC()
{
delete _v;
for(unsigned int m = 0; m < _N; ++m) {
delete _uHABC[m];
}
}
unsigned int HABC::getN() const
{
return _N;
}
double HABC::getTheta() const
{
return _theta;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * HABC::getUHABC(const unsigned int m) const
{
if(m > _N) throw gmshfem::common::Exception("Try to get uHABC field " + std::to_string(m));
return _uHABC[m];
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* HABC::getV() const
{
return _v;
}
void HABC::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const unsigned int fieldOrder = parameters.getFieldOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
_uHABC.resize(_N);
for(unsigned int m = 0; m < _N; ++m) {
_uHABC[m] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uHABC_" + orientation() + "_" + std::to_string(m), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
}
_v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
const double pi = 3.141592653589793238462643383279;
const std::complex< double > expPTheta(std::cos(_theta), std::sin(_theta));
const std::complex< double > expPTheta2(std::cos(_theta/2.), std::sin(_theta/2.));
const double M = 2 * _N + 1;
double cPade[_N];
for(unsigned int m = 0; m < _N; ++m) {
cPade[m] = std::tan((m+1) * pi / M) * std::tan((m+1) * pi / M);
}
const std::complex< double > im(0., 1.);
formulation.integral(-dof(*_v), tf(*u), sigma, gauss);
formulation.integral( dof(*_v), tf(*_v), sigma, gauss);
formulation.integral( im * kappa * expPTheta2 * dof(*u), tf(*_v), sigma, gauss);
for(unsigned int m = 0; m < _N; ++m) {
formulation.integral( im * kappa * expPTheta2 * 2./M * cPade[m] * dof(*u), tf(*_v), sigma, gauss);
formulation.integral( im * kappa * expPTheta2 * 2./M * cPade[m] * dof(*_uHABC[m]), tf(*_v), sigma, gauss); }
// auxiliary
for(unsigned int m = 0; m < _N; ++m) {
formulation.integral(grad(dof(*_uHABC[m])), grad(tf(*_uHABC[m])), sigma, gauss);
formulation.integral(- kappa * kappa * (expPTheta * cPade[m] + 1.) * dof(*_uHABC[m]), tf(*_uHABC[m]), sigma, gauss);
formulation.integral(- kappa * kappa * expPTheta * (cPade[m] + 1.) * dof(*u), tf(*_uHABC[m]), sigma, gauss);
}
}
// **********************************
// PML
// **********************************
Pml::Pml(const unsigned int boundary, const double pmlSize, const std::string type) : Boundary(boundary), _size(pmlSize), _type(type)
{
}
Pml::~Pml()
{
}
double Pml::getSize() const
{
return _size;
}
gmshfem::function::ScalarFunction< std::complex< double > > Pml::pmlCoefficients(gmshfem::function::TensorFunction< std::complex< double > > &D, gmshfem::function::ScalarFunction< std::complex< double > > &E, gmshfem::domain::Domain &bnd, const gmshfem::function::ScalarFunction< std::complex< double > > &kappa) const
{
gmshfem::function::ScalarFunction< std::complex< double > > kappaMod;
const std::complex< double > im(0., 1.);
gmshfem::function::ScalarFunction< std::complex< double > > distSigma;
gmshfem::function::ScalarFunction< std::complex< double > > sigma;
if(orientation() == "E" || orientation() == "W") {
distSigma = abs(gmshfem::function::x< std::complex< double > >() - gmshfem::function::x< std::complex< double > >(bnd));
if(bnd.maxDim() == 0 || bnd.maxDim() == 1) {
kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >(bnd));
}
else {
kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(), gmshfem::function::z< double >());
}
}
else if(orientation() == "N" || orientation() == "S") {
distSigma = abs(gmshfem::function::y< std::complex< double > >() - gmshfem::function::y< std::complex< double > >(bnd));
if(bnd.maxDim() == 0 || bnd.maxDim() == 1) {
kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >(bnd));
}
else {
kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >());
}
}
else if(orientation() == "D" || orientation() == "U") {
distSigma = abs(gmshfem::function::z< std::complex< double > >() - gmshfem::function::z< std::complex< double > >(bnd));
if(bnd.maxDim() == 0 || bnd.maxDim() == 1) {
kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >(bnd));
}
else {
kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(), gmshfem::function::y< double >(), gmshfem::function::z< double >(bnd));
}
}
if(_type == "hs") {
sigma = 1. / (_size - distSigma) - 1./ _size;
}
else if(_type == "h") {
sigma = 1. / (_size - distSigma);
}
else if(_type == "q") {
const double sigmaStar = 664.55; sigma = sigmaStar * distSigma * distSigma / (_size * _size);
}
gmshfem::function::ScalarFunction< std::complex< double > > kMod = 1. + im * sigma / kappaMod;
if(orientation() == "E" || orientation() == "W") {
D = gmshfem::function::tensorDiag< std::complex< double > >(1./kMod, kMod, kMod);
E = kMod;
}
else if(orientation() == "N" || orientation() == "S") {
D = gmshfem::function::tensorDiag< std::complex< double > >(kMod, 1./kMod, kMod);
E = kMod;
}
else if(orientation() == "D" || orientation() == "U") {
D = gmshfem::function::tensorDiag< std::complex< double > >(kMod, kMod, 1./kMod);
E = kMod;
}
return kappaMod;
}
// **********************************
// PML continuous
// **********************************
PmlContinuous::PmlContinuous(const unsigned int boundary, const double pmlSize, const std::string type) : Pml(boundary, pmlSize, type)
{
}
PmlContinuous::~PmlContinuous()
{
delete _v;
delete _uPml;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlContinuous::getUPml() const
{
return _uPml;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlContinuous::getV() const
{
return _v;
}
void PmlContinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
gmshfem::domain::Domain pml = domains.getPml(_boundary);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const unsigned int fieldOrder = parameters.getFieldOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
_uPml = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPml_" + orientation(), pml, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
_v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
gmshfem::function::TensorFunction< std::complex< double > > D;
gmshfem::function::ScalarFunction< std::complex< double > > E;
gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = pmlCoefficients(D, E, sigma, kappa);
formulation.integral(- D * grad(dof(*_uPml)), grad(tf(*_uPml)), pml, gauss);
formulation.integral(kappaMod * kappaMod * E * dof(*_uPml), tf(*_uPml), pml, gauss);
formulation.integral(dof(*_v), tf(*_uPml), sigma, gauss);
formulation.integral(- dof(*_v), tf(*u), sigma, gauss);
formulation.integral(dof(*_uPml), tf(*_v), sigma, gauss); formulation.integral(- dof(*u), tf(*_v), sigma, gauss);
}
// **********************************
// PML discontinuous
// **********************************
PmlDiscontinuous::PmlDiscontinuous(const unsigned int boundary, const double pmlSize, const std::string type) : Pml(boundary, pmlSize, type)
{
}
PmlDiscontinuous::~PmlDiscontinuous()
{
delete _v;
delete _uPml;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlDiscontinuous::getUPml() const
{
return _uPml;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 > * PmlDiscontinuous::getV() const
{
return _v;
}
void PmlDiscontinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
// Marche pas, il faudrait du vrai H(div) ...
// gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
// gmshfem::domain::Domain pml = domains.getPml(_boundary);
//
// const unsigned int neumannOrder = parameters.getNeumannOrder();
// const unsigned int fieldOrder = parameters.getFieldOrder();
// const std::string gauss = parameters.getGauss();
// const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
//
// _uPml = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPml_" + orientation(), pml, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
//
// _v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
//
// gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
//
// gmshfem::function::TensorFunction< std::complex< double > > D;
// gmshfem::function::ScalarFunction< std::complex< double > > E;
// gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = pmlCoefficients(D, E, sigma, kappa);
//
// formulation.integral(- D * grad(dof(*_uPml)), grad(tf(*_uPml)), pml, gauss);
// formulation.integral(kappaMod * kappaMod * E * dof(*_uPml), tf(*_uPml), pml, gauss);
//
// formulation.integral(dof(*_v), tf(*_uPml), sigma, gauss);
// formulation.integral(- dof(*_v), tf(*u), sigma, gauss);
//
// formulation.integral(dof(*_uPml), tf(*_v), sigma, gauss);
// formulation.integral(- dof(*u), tf(*_v), sigma, gauss);
}
// **********************************
// PML without multiplier
// **********************************
PmlWithoutMultiplier::PmlWithoutMultiplier(const unsigned int boundary, const double pmlSize, const std::string type) : Pml(boundary, pmlSize, type)
{
}
PmlWithoutMultiplier::~PmlWithoutMultiplier()
{
delete _v;
delete _uPml; }
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlWithoutMultiplier::getUPml() const
{
return _uPml;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlWithoutMultiplier::getV() const
{
return _v;
}
void PmlWithoutMultiplier::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
gmshfem::domain::Domain pml = domains.getPml(_boundary);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const unsigned int fieldOrder = parameters.getFieldOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
_uPml = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPml_" + orientation(), pml, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
_v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
gmshfem::function::TensorFunction< std::complex< double > > D;
gmshfem::function::ScalarFunction< std::complex< double > > E;
gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = pmlCoefficients(D, E, sigma, kappa);
formulation.integral(- D * grad(dof(*_uPml)), grad(tf(*_uPml)), pml, gauss);
formulation.integral(kappaMod * kappaMod * E * dof(*_uPml), tf(*_uPml), pml, gauss);
formulation.integral(dof(*_v), tf(*_uPml), sigma, gauss);
formulation.integral(- dof(*_v), tf(*u), sigma, gauss);
formulation.integral(dof(*_uPml), tf(*_v), sigma, gauss);
formulation.integral(- dof(*u), tf(*_v), sigma, gauss);
}
// **********************************
// HABC_HABC
// **********************************
HABC_HABC::HABC_HABC(const unsigned int edge, HABC *first, HABC *second) : Edge(edge, first, second)
{
}
HABC_HABC::~HABC_HABC()
{
const unsigned int N0 = static_cast< HABC * >(_bnd[0])->getN();
const unsigned int N1 = static_cast< HABC * >(_bnd[1])->getN();
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N1; ++n) {
delete _uHABC_E[m][n];
}
}
for(unsigned int m = 0; m < N0; ++m) {
delete _vE[0][m];
}
for(unsigned int n = 0; n < N1; ++n) {
delete _vE[1][n];
}
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC::getUHABC(const unsigned int m, const unsigned int n) const {
return _uHABC_E[m][n];
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC::getV(const unsigned int m, const unsigned int i) const
{
return _vE[i][m];
}
HABC *HABC_HABC::firstBoundary() const
{
return static_cast< HABC * >(_bnd[0]);
}
HABC *HABC_HABC::secondBoundary() const
{
return static_cast< HABC * >(_bnd[1]);
}
void HABC_HABC::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain edge = domains.getEdge(_edge);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const unsigned int fieldOrder = parameters.getFieldOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
const unsigned int N0 = static_cast< HABC * >(_bnd[0])->getN();
const unsigned int N1 = static_cast< HABC * >(_bnd[1])->getN();
const double theta0 = static_cast< HABC * >(_bnd[0])->getTheta();
const double theta1 = static_cast< HABC * >(_bnd[1])->getTheta();
_uHABC_E.resize(N0);
for(unsigned int m = 0; m < N0; ++m) {
_uHABC_E[m].resize(N1);
for(unsigned int n = 0; n < N1; ++n) {
_uHABC_E[m][n] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uHABC_Edge_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(n), edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
}
}
_vE[0].resize(N0);
for(unsigned int m = 0; m < N0; ++m) {
_vE[0][m] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vE_" + orientation() + "_" + std::to_string(m) + "_first", edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
}
_vE[1].resize(N1);
for(unsigned int n = 0; n < N1; ++n) {
_vE[1][n] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vE_" + orientation() + "_" + std::to_string(n) + "_second", edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
}
const double pi = 3.141592653589793238462643383279;
const std::complex< double > expPTheta_0(std::cos(theta0), std::sin(theta0));
const std::complex< double > expPTheta2_0(std::cos(theta0/2.), std::sin(theta0/2.));
const std::complex< double > expPTheta_1(std::cos(theta1), std::sin(theta1));
const std::complex< double > expPTheta2_1(std::cos(theta1/2.), std::sin(theta1/2.));
const double M0 = 2 * N0 + 1;
double cPade0[N0];
for(unsigned int m = 0; m < N0; ++m) {
cPade0[m] = std::tan((m+1) * pi / M0) * std::tan((m+1) * pi / M0);
}
const double M1 = 2 * N1 + 1;
double cPade1[N1];
for(unsigned int n = 0; n < N1; ++n) {
cPade1[n] = std::tan((n+1) * pi / M1) * std::tan((n+1) * pi / M1);
}
const std::complex< double > im(0., 1.);
std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > uHABC[2];
uHABC[0].resize(N0);
uHABC[1].resize(N1);
for(unsigned int m = 0; m < N0; ++m) {
uHABC[0][m] = static_cast< HABC * >(_bnd[0])->getUHABC(m);
formulation.integral(-dof(*_vE[0][m]), tf(*uHABC[0][m]), edge, gauss);
formulation.integral( dof(*_vE[0][m]), tf(*_vE[0][m]), edge, gauss);
formulation.integral(-im * kappa * expPTheta2_1 * dof(*uHABC[0][m]), tf(*_vE[0][m]), edge, gauss);
for(unsigned int n = 0; n < N1; ++n) {
formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*uHABC[0][m]), tf(*_vE[0][m]), edge, gauss);
formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*_uHABC_E[m][n]), tf(*_vE[0][m]), edge, gauss);
}
}
for(unsigned int n = 0; n < N1; ++n) {
uHABC[1][n] = static_cast< HABC * >(_bnd[1])->getUHABC(n);
formulation.integral(-dof(*_vE[1][n]), tf(*uHABC[1][n]), edge, gauss);
formulation.integral( dof(*_vE[1][n]), tf(*_vE[1][n]), edge, gauss);
formulation.integral(-im * kappa * expPTheta2_0 * dof(*uHABC[1][n]), tf(*_vE[1][n]), edge, gauss);
for(unsigned int m = 0; m < N0; ++m) {
formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*uHABC[1][n]), tf(*_vE[1][n]), edge, gauss);
formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*_uHABC_E[m][n]), tf(*_vE[1][n]), edge, gauss);
}
}
// auxiliary
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N1; ++n) {
formulation.integral(grad(dof(*_uHABC_E[m][n])), grad(tf(*_uHABC_E[m][n])), edge, gauss);
formulation.integral(- kappa * kappa * (expPTheta_0 * cPade0[m] + expPTheta_1 * cPade1[n] + 1.) * dof(*_uHABC_E[m][n]), tf(*_uHABC_E[m][n]), edge, gauss);
formulation.integral(- kappa * kappa * expPTheta_1 * (cPade1[n] + 1.) * dof(*uHABC[0][m]), tf(*_uHABC_E[m][n]), edge, gauss);
formulation.integral(- kappa * kappa * expPTheta_0 * (cPade0[m] + 1.) * dof(*uHABC[1][n]), tf(*_uHABC_E[m][n]), edge, gauss);
}
}
}
// **********************************
// PML continuous _ PML continuous
// **********************************
PmlContinuous_PmlContinuous::PmlContinuous_PmlContinuous(const unsigned int edge, PmlContinuous *first, PmlContinuous *second) : Edge(edge, first, second)
{
}
PmlContinuous_PmlContinuous::~PmlContinuous_PmlContinuous()
{
delete _uPmlEdge;
delete _vC[0];
delete _vC[1];
delete _vEdge;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous::getUPml() const
{
return _uPmlEdge;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous::getV(const unsigned int i) const
{
return _vC[i];
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous::getVEdge() const
{
return _vEdge;
}
PmlContinuous *PmlContinuous_PmlContinuous::firstBoundary() const
{
return static_cast< PmlContinuous * >(_bnd[0]);
}
PmlContinuous *PmlContinuous_PmlContinuous::secondBoundary() const
{
return static_cast< PmlContinuous * >(_bnd[1]);
}
void PmlContinuous_PmlContinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain edge = domains.getEdge(_edge);
gmshfem::domain::Domain pmlEdge = domains.getPmlEdge(_edge);
std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > pmlBnd = domains.getPmlBnd(_edge);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const unsigned int fieldOrder = parameters.getFieldOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
_uPmlEdge = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPmlEdge_" + orientation(), pmlEdge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
_vC[0] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_first", pmlBnd.first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
_vC[1] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_second", pmlBnd.second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *uPml[2];
uPml[0] = static_cast< PmlContinuous * >(_bnd[0])->getUPml();
uPml[1] = static_cast< PmlContinuous * >(_bnd[1])->getUPml();
gmshfem::function::TensorFunction< std::complex< double > > D[2];
gmshfem::function::ScalarFunction< std::complex< double > > E[2];
gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = static_cast< PmlContinuous * >(_bnd[0])->pmlCoefficients(D[0], E[0], edge, kappa);
static_cast< PmlContinuous * >(_bnd[1])->pmlCoefficients(D[1], E[1], edge, kappa);
formulation.integral(- D[0] * D[1] * grad(dof(*_uPmlEdge)), grad(tf(*_uPmlEdge)), pmlEdge, gauss);
formulation.integral(kappaMod * kappaMod * E[0] * E[1] * dof(*_uPmlEdge), tf(*_uPmlEdge), pmlEdge, gauss);
formulation.integral(dof(*_vC[0]), tf(*_uPmlEdge), pmlBnd.first, gauss);
formulation.integral(- dof(*_vC[0]), tf(*uPml[0]), pmlBnd.first, gauss);
formulation.integral(dof(*_uPmlEdge), tf(*_vC[0]), pmlBnd.first, gauss);
formulation.integral(- dof(*uPml[0]), tf(*_vC[0]), pmlBnd.first, gauss);
formulation.integral(dof(*_vC[1]), tf(*_uPmlEdge), pmlBnd.second, gauss);
formulation.integral(- dof(*_vC[1]), tf(*uPml[1]), pmlBnd.second, gauss);
formulation.integral(dof(*_uPmlEdge), tf(*_vC[1]), pmlBnd.second, gauss);
formulation.integral(- dof(*uPml[1]), tf(*_vC[1]), pmlBnd.second, gauss);
// edge equation
_vEdge = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vEdge_" + orientation(), edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
formulation.integral(dof(*_vC[0]), tf(*_vEdge), edge, gauss);
formulation.integral(dof(*_vEdge), tf(*_vC[0]), edge, gauss);
formulation.integral(-dof(*_vC[1]), tf(*_vEdge), edge, gauss);
formulation.integral(-dof(*_vEdge), tf(*_vC[1]), edge, gauss);
formulation.integral(dof(*static_cast< PmlContinuous * >(_bnd[0])->getV()), tf(*_vEdge), edge, gauss);
formulation.integral(dof(*_vEdge), tf(*static_cast< PmlContinuous * >(_bnd[0])->getV()), edge, gauss);
formulation.integral(-dof(*static_cast< PmlContinuous * >(_bnd[1])->getV()), tf(*_vEdge), edge, gauss);
formulation.integral(-dof(*_vEdge), tf(*static_cast< PmlContinuous * >(_bnd[1])->getV()), edge, gauss);
}
// **********************************
// HABC_HABC_HABC
// **********************************
HABC_HABC_HABC::HABC_HABC_HABC(const unsigned int corner, HABC_HABC *first, HABC_HABC *second, HABC_HABC *third) : Corner(corner, first, second, third)
{
}
HABC_HABC_HABC::~HABC_HABC_HABC()
{
std::vector< HABC * > boundaries(3, nullptr);
std::map< std::string, HABC * > boundaryMap;
for(unsigned int i = 0; i < 3; ++i) {
boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()));
boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()));
}
for(auto it = boundaryMap.begin(); it != boundaryMap.end(); ++it) {
if(it->first == "D" || it->first == "U") {
boundaries[2] = it->second;
}
else {
if(boundaries[0] == nullptr) {
boundaries[0] = it->second;
}
else {
boundaries[1] = it->second;
}
}
}
if((boundaries[0]->orientation() == "N" && boundaries[1]->orientation() == "E") ||
(boundaries[0]->orientation() == "W" && boundaries[1]->orientation() == "N") ||
(boundaries[0]->orientation() == "S" && boundaries[1]->orientation() == "W") ||
(boundaries[0]->orientation() == "E" && boundaries[1]->orientation() == "S")) {
std::swap(boundaries[0], boundaries[1]);
}
const unsigned int N0 = static_cast< HABC * >(boundaries[0])->getN();
const unsigned int N1 = static_cast< HABC * >(boundaries[1])->getN();
const unsigned int N2 = static_cast< HABC * >(boundaries[2])->getN();
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N1; ++n) {
for(unsigned int o = 0; o < N2; ++o) {
delete _uHABC_C[m][n][o];
}
}
}
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N2; ++n) {
delete _vC[0][m][n];
}
}
for(unsigned int n = 0; n < N1; ++n) {
for(unsigned int o = 0; o < N2; ++o) {
delete _vC[1][n][o];
}
}
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N1; ++n) {
delete _vC[2][m][n];
}
}
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC_HABC::getUHABC(const unsigned int m, const unsigned int n, const unsigned int o) const
{
return _uHABC_C[m][n][o];
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC_HABC::getV(const unsigned int m, const unsigned int n, const unsigned int i) const
{
return _vC[i][m][n];
}
HABC_HABC *HABC_HABC_HABC::firstEdge() const
{
return static_cast< HABC_HABC * >(_bnd[0]);
}
HABC_HABC *HABC_HABC_HABC::secondEdge() const
{
return static_cast< HABC_HABC * >(_bnd[1]);
}
HABC_HABC *HABC_HABC_HABC::thirdEdge() const
{
return static_cast< HABC_HABC * >(_bnd[2]);
}
void HABC_HABC_HABC::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain corner = domains.getCorner(_corner);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const unsigned int fieldOrder = parameters.getFieldOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
std::vector< HABC * > boundaries(3, nullptr);
std::map< std::string, HABC * > boundaryMap;
for(unsigned int i = 0; i < 3; ++i) {
boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()));
boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()));
}
for(auto it = boundaryMap.begin(); it != boundaryMap.end(); ++it) {
if(it->first == "D" || it->first == "U") {
boundaries[2] = it->second;
}
else {
if(boundaries[0] == nullptr) {
boundaries[0] = it->second;
}
else {
boundaries[1] = it->second;
}
}
}
if((boundaries[0]->orientation() == "N" && boundaries[1]->orientation() == "E") ||
(boundaries[0]->orientation() == "W" && boundaries[1]->orientation() == "N") ||
(boundaries[0]->orientation() == "S" && boundaries[1]->orientation() == "W") ||
(boundaries[0]->orientation() == "E" && boundaries[1]->orientation() == "S")) {
std::swap(boundaries[0], boundaries[1]);
}
const unsigned int N0 = static_cast< HABC * >(boundaries[0])->getN();
const unsigned int N1 = static_cast< HABC * >(boundaries[1])->getN();
const unsigned int N2 = static_cast< HABC * >(boundaries[2])->getN();
const double theta0 = static_cast< HABC * >(boundaries[0])->getTheta();
const double theta1 = static_cast< HABC * >(boundaries[1])->getTheta();
const double theta2 = static_cast< HABC * >(boundaries[2])->getTheta();
_uHABC_C.resize(N0);
for(unsigned int m = 0; m < N0; ++m) {
_uHABC_C[m].resize(N1);
for(unsigned int n = 0; n < N1; ++n) {
_uHABC_C[m][n].resize(N2);
for(unsigned int o = 0; o < N2; ++o) {
_uHABC_C[m][n][o] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uHABC_Corner_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(n) + "_" + std::to_string(o), corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
}
}
}
_vC[0].resize(N0);
for(unsigned int m = 0; m < N0; ++m) {
_vC[0][m].resize(N2);
for(unsigned int o = 0; o < N2; ++o) {
_vC[0][m][o] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(o) + "_first", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
}
}
_vC[1].resize(N1);
for(unsigned int n = 0; n < N1; ++n) {
_vC[1][n].resize(N2);
for(unsigned int o = 0; o < N2; ++o) {
_vC[1][n][o] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(n) + "_" + std::to_string(o) + "_second", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
}
}
_vC[2].resize(N0);
for(unsigned int m = 0; m < N0; ++m) {
_vC[2][m].resize(N1);
for(unsigned int n = 0; n < N1; ++n) {
_vC[2][m][n] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(n) + "_third", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
}
}
const double pi = 3.141592653589793238462643383279;
const std::complex< double > expPTheta_0(std::cos(theta0), std::sin(theta0));
const std::complex< double > expPTheta2_0(std::cos(theta0/2.), std::sin(theta0/2.));
const std::complex< double > expPTheta_1(std::cos(theta1), std::sin(theta1));
const std::complex< double > expPTheta2_1(std::cos(theta1/2.), std::sin(theta1/2.));
const std::complex< double > expPTheta_2(std::cos(theta2), std::sin(theta2));
const std::complex< double > expPTheta2_2(std::cos(theta2/2.), std::sin(theta2/2.));
const double M0 = 2 * N0 + 1;
double cPade0[N0];
for(unsigned int m = 0; m < N0; ++m) {
cPade0[m] = std::tan((m+1) * pi / M0) * std::tan((m+1) * pi / M0);
}
const double M1 = 2 * N1 + 1;
double cPade1[N1];
for(unsigned int n = 0; n < N1; ++n) {
cPade1[n] = std::tan((n+1) * pi / M1) * std::tan((n+1) * pi / M1);
}
const double M2 = 2 * N2 + 1;
double cPade2[N2];
for(unsigned int o = 0; o < N2; ++o) {
cPade2[o] = std::tan((o+1) * pi / M2) * std::tan((o+1) * pi / M2);
}
const std::complex< double > im(0., 1.);
std::vector< std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > > uHABC[3];
uHABC[0].resize(N0, std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(N2, nullptr));
uHABC[1].resize(N1, std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(N2, nullptr));
uHABC[2].resize(N0, std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(N1, nullptr));
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int o = 0; o < N2; ++o) {
uHABC[0][m][o] = static_cast< HABC_HABC * >(_bnd[0])->getUHABC(m,o);
formulation.integral(-dof(*_vC[0][m][o]), tf(*uHABC[0][m][o]), corner, gauss);
formulation.integral( dof(*_vC[0][m][o]), tf(*_vC[0][m][o]), corner, gauss);
formulation.integral(-im * kappa * expPTheta2_1 * dof(*uHABC[0][m][o]), tf(*_vC[0][m][o]), corner, gauss);
for(unsigned int n = 0; n < N1; ++n) {
formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*uHABC[0][m][o]), tf(*_vC[0][m][o]), corner, gauss);
formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*_uHABC_C[m][n][o]), tf(*_vC[0][m][o]), corner, gauss);
}
}
}
for(unsigned int n = 0; n < N1; ++n) {
for(unsigned int o = 0; o < N2; ++o) {
uHABC[1][n][o] = static_cast< HABC_HABC * >(_bnd[1])->getUHABC(n,o);
formulation.integral(-dof(*_vC[1][n][o]), tf(*uHABC[1][n][o]), corner, gauss);
formulation.integral( dof(*_vC[1][n][o]), tf(*_vC[1][n][o]), corner, gauss);
formulation.integral(-im * kappa * expPTheta2_0 * dof(*uHABC[1][n][o]), tf(*_vC[1][n][o]), corner, gauss);
for(unsigned int m = 0; m < N0; ++m) {
formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*uHABC[1][n][o]), tf(*_vC[1][n][o]), corner, gauss);
formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*_uHABC_C[m][n][o]), tf(*_vC[1][n][o]), corner, gauss);
}
}
}
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N1; ++n) {
uHABC[2][m][n] = static_cast< HABC_HABC * >(_bnd[2])->getUHABC(m,n);
formulation.integral(-dof(*_vC[2][m][n]), tf(*uHABC[2][m][n]), corner, gauss);
formulation.integral( dof(*_vC[2][m][n]), tf(*_vC[2][m][n]), corner, gauss);
formulation.integral(-im * kappa * expPTheta2_2 * dof(*uHABC[2][m][n]), tf(*_vC[2][m][n]), corner, gauss);
for(unsigned int o = 0; o < N2; ++o) {
formulation.integral(-im * kappa * expPTheta2_2 * 2./M2 * cPade2[o] * dof(*uHABC[2][m][n]), tf(*_vC[2][m][n]), corner, gauss);
formulation.integral(-im * kappa * expPTheta2_2 * 2./M2 * cPade2[o] * dof(*_uHABC_C[m][n][o]), tf(*_vC[2][m][n]), corner, gauss);
}
}
}
// auxiliary
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N1; ++n) {
for(unsigned int o = 0; o < N2; ++o) {
formulation.integral( (expPTheta_0 * cPade0[m] + expPTheta_1 * cPade1[n] + expPTheta_2 * cPade2[o] + 1.) * dof(*_uHABC_C[m][n][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
formulation.integral( expPTheta_1 * (cPade1[n] + 1.) * dof(*uHABC[0][m][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
formulation.integral( expPTheta_0 * (cPade0[m] + 1.) * dof(*uHABC[1][n][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
formulation.integral( expPTheta_2 * (cPade2[o] + 1.) * dof(*uHABC[2][m][n]), tf(*_uHABC_C[m][n][o]), corner, gauss);
}
}
}
}
// **********************************
// PML continuous _ PML continuous _ PML continuous
// **********************************
PmlContinuous_PmlContinuous_PmlContinuous::PmlContinuous_PmlContinuous_PmlContinuous(const unsigned int corner, PmlContinuous_PmlContinuous *first, PmlContinuous_PmlContinuous *second, PmlContinuous_PmlContinuous *third) : Corner(corner, first, second, third)
{
}
PmlContinuous_PmlContinuous_PmlContinuous::~PmlContinuous_PmlContinuous_PmlContinuous()
{
delete _uPmlCorner;
delete _vC[0]; delete _vC[1];
delete _vC[2];
delete _vEdge[0];
delete _vEdge[1];
delete _vEdge[2];
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous_PmlContinuous::getUPml() const
{
return _uPmlCorner;
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous_PmlContinuous::getV(const unsigned int i) const
{
return _vC[i];
}
PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous::firstEdge() const
{
return static_cast< PmlContinuous_PmlContinuous * >(_bnd[0]);
}
PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous::secondEdge() const
{
return static_cast< PmlContinuous_PmlContinuous * >(_bnd[1]);
}
PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous::thirdEdge() const
{
return static_cast< PmlContinuous_PmlContinuous * >(_bnd[2]);
}
void PmlContinuous_PmlContinuous_PmlContinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters ¶meters)
{
gmshfem::domain::Domain corner = domains.getCorner(_corner);
gmshfem::domain::Domain pmlCorner = domains.getPmlCorner(_corner);
std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > pmlBnd = domains.getPmlEdgeBnd(_corner);
std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > cornerEdge = domains.getCornerEdge(_corner);
const unsigned int neumannOrder = parameters.getNeumannOrder();
const unsigned int fieldOrder = parameters.getFieldOrder();
const std::string gauss = parameters.getGauss();
const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
_uPmlCorner = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPmlCorner_" + orientation(), pmlCorner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
_vC[0] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_first", std::get<0>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
_vC[1] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_second", std::get<1>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
_vC[2] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_third", std::get<2>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<0>(pmlBnd), "vC_" + orientation() + "_first");
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<1>(pmlBnd), "vC_" + orientation() + "_second");
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<2>(pmlBnd), "vC_" + orientation() + "_third");
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *uPml[3];
uPml[0] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getUPml();
uPml[1] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getUPml();
uPml[2] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getUPml();
std::vector< PmlContinuous * > boundaries(6);
for(unsigned int i = 0; i < 3; ++i) {
boundaries[2 * i] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[i])->firstBoundary();
boundaries[2 * i + 1] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[i])->secondBoundary();
}
std::sort(boundaries.begin(), boundaries.end());
gmshfem::function::TensorFunction< std::complex< double > > D[3];
gmshfem::function::ScalarFunction< std::complex< double > > E[3];
gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = boundaries[0]->pmlCoefficients(D[0], E[0], corner, kappa);
boundaries[2]->pmlCoefficients(D[1], E[1], corner, kappa);
boundaries[4]->pmlCoefficients(D[2], E[2], corner, kappa);
formulation.integral(- D[0] * D[1] * D[2] * grad(dof(*_uPmlCorner)), grad(tf(*_uPmlCorner)), pmlCorner, gauss);
formulation.integral(kappaMod * kappaMod * E[0] * E[1] * E[2] * dof(*_uPmlCorner), tf(*_uPmlCorner), pmlCorner, gauss);
formulation.integral(dof(*_vC[0]), tf(*_uPmlCorner), std::get<0>(pmlBnd), gauss);
formulation.integral(- dof(*_vC[0]), tf(*uPml[0]), std::get<0>(pmlBnd), gauss);
formulation.integral(dof(*_uPmlCorner), tf(*_vC[0]), std::get<0>(pmlBnd), gauss);
formulation.integral(- dof(*uPml[0]), tf(*_vC[0]), std::get<0>(pmlBnd), gauss);
formulation.integral(dof(*_vC[1]), tf(*_uPmlCorner), std::get<1>(pmlBnd), gauss);
formulation.integral(- dof(*_vC[1]), tf(*uPml[1]), std::get<1>(pmlBnd), gauss);
formulation.integral(dof(*_uPmlCorner), tf(*_vC[1]), std::get<1>(pmlBnd), gauss);
formulation.integral(- dof(*uPml[1]), tf(*_vC[1]), std::get<1>(pmlBnd), gauss);
formulation.integral(dof(*_vC[2]), tf(*_uPmlCorner), std::get<2>(pmlBnd), gauss);
formulation.integral(- dof(*_vC[2]), tf(*uPml[2]), std::get<2>(pmlBnd), gauss);
formulation.integral(dof(*_uPmlCorner), tf(*_vC[2]), std::get<2>(pmlBnd), gauss);
formulation.integral(- dof(*uPml[2]), tf(*_vC[2]), std::get<2>(pmlBnd), gauss);
// corner equation
_vEdge[0] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCornerEdge_" + orientation() + "_first", std::get<0>(cornerEdge), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
formulation.integral(dof(*_vC[0]), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
formulation.integral(dof(*_vEdge[0]), tf(*_vC[0]), std::get<0>(cornerEdge), gauss);
formulation.integral(-dof(*_vC[2]), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
formulation.integral(-dof(*_vEdge[0]), tf(*_vC[2]), std::get<0>(cornerEdge), gauss);
formulation.integral(dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(0)), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
formulation.integral(dof(*_vEdge[0]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(0)), std::get<0>(cornerEdge), gauss);
formulation.integral(-dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(0)), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
formulation.integral(-dof(*_vEdge[0]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(0)), std::get<0>(cornerEdge), gauss);
_vEdge[1] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCornerEdge_" + orientation() + "_second", std::get<1>(cornerEdge), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
formulation.integral(dof(*_vC[1]), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
formulation.integral(dof(*_vEdge[1]), tf(*_vC[1]), std::get<1>(cornerEdge), gauss);
formulation.integral(-dof(*_vC[2]), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
formulation.integral(-dof(*_vEdge[1]), tf(*_vC[2]), std::get<1>(cornerEdge), gauss);
formulation.integral(dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(0)), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
formulation.integral(dof(*_vEdge[1]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(0)), std::get<1>(cornerEdge), gauss);
formulation.integral(-dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(1)), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
formulation.integral(-dof(*_vEdge[1]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(1)), std::get<1>(cornerEdge), gauss);
_vEdge[2] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCornerEdge_" + orientation() + "_third", std::get<2>(cornerEdge), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
formulation.integral(dof(*_vC[0]), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
formulation.integral(dof(*_vEdge[2]), tf(*_vC[0]), std::get<2>(cornerEdge), gauss);
formulation.integral(-dof(*_vC[1]), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
formulation.integral(-dof(*_vEdge[2]), tf(*_vC[1]), std::get<2>(cornerEdge), gauss);
formulation.integral(dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(1)), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
formulation.integral(dof(*_vEdge[2]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(1)), std::get<2>(cornerEdge), gauss);
formulation.integral(-dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(1)), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
formulation.integral(-dof(*_vEdge[2]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(1)), std::get<2>(cornerEdge), gauss);
_vCorner = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCorner_" + orientation(), corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
formulation.integral(dof(*_vEdge[0]), tf(*_vCorner), corner, gauss);
formulation.integral(dof(*_vCorner), tf(*_vEdge[0]), corner, gauss);
formulation.integral(dof(*_vEdge[1]), tf(*_vCorner), corner, gauss);
formulation.integral(dof(*_vCorner), tf(*_vEdge[1]), corner, gauss);
formulation.integral(dof(*_vEdge[2]), tf(*_vCorner), corner, gauss);
formulation.integral(dof(*_vCorner), tf(*_vEdge[2]), corner, gauss);
}
}