Select Git revision
monoDomain.cpp
monoDomain.cpp 31.94 KiB
#include "mesh.h"
#include "SubproblemDomains.h"
#include "Subproblem2D.h"
#include "Subproblem3D.h"
#include <gmshfem/GmshFem.h>
#include <gmshfem/FieldInterface.h>
#include <gmshfem/Formulation.h>
#include <gmshfem/AnalyticalFunction.h>
#include <gmshfem/Post.h>
#include <gmshfem/Function.h>
#include <gmshfem/io.h>
#include <gmshddm/GmshDdm.h>
#include <gmshddm/Subdomain.h>
#include <gmshddm/Interface.h>
#include <gmshddm/SubdomainField.h>
#include <gmshddm/InterfaceField.h>
#include <gmshddm/Formulation.h>
#include <gmshddm/MPIInterface.h>
#include <algorithm>
#include <fstream>
using gmshfem::equation::dof;
using gmshfem::equation::tf;
using gmshfem::function::operator-;
using gmshfem::function::operator*;
using gmshfem::function::abs;
namespace D2 {
void monoDomain()
{
gmshddm::common::GmshDdm *gmshDdm = gmshddm::common::GmshDdm::currentInstance();
// ************************
// P H Y S I C S
// ************************
double pi = 3.141592653589793238462643383279;
double k = 4. * pi;
gmshDdm->userDefinedParameter(k, "k");
double R = 0.5;
gmshDdm->userDefinedParameter(R, "R");
// ************************
// M E S H
// ************************
double lc = 0.03333333333333; // Marmousi = 20; other = 0.03333333333333
gmshDdm->userDefinedParameter(lc, "lc");
int meshOrder = 1;
gmshDdm->userDefinedParameter(meshOrder, "meshOrder");
// ************************
// M O D E L I N G
// ************************
std::string boundary = "sommerfeld"; // sommerfeld, pml, habc
gmshDdm->userDefinedParameter(boundary, "boundary");
std::string pmlMethod = "continuous"; // continuous, discontinuous, withoutMultiplier
gmshDdm->userDefinedParameter(pmlMethod, "pmlMethod");
std::string pmlType = "hs"; // hs, h, q
gmshDdm->userDefinedParameter(pmlType, "pmlType");
double N = 6;
gmshDdm->userDefinedParameter(N, "N");
double pmlSize = N * lc;
if(pmlSize == 0.) {
pmlSize = 0.3;
}
std::string gauss = "Gauss10";
gmshDdm->userDefinedParameter(gauss, "gauss");
int fieldOrder = 1;
gmshDdm->userDefinedParameter(fieldOrder, "fieldOrder");
int neumannOrder = 1;
gmshDdm->userDefinedParameter(neumannOrder, "neumannOrder");
double stab = 0.11;
gmshDdm->userDefinedParameter(stab, "stab");
double thetaPade = 0.;
gmshDdm->userDefinedParameter(thetaPade, "thetaPade");
thetaPade *= pi;
// ************************
// P O S T
// ************************
bool computeError = false;
gmshDdm->userDefinedParameter(computeError, "error");
bool computeInterfaceError = false;
gmshDdm->userDefinedParameter(computeInterfaceError, "interfaceError");
bool getMatrices = false;
gmshDdm->userDefinedParameter(getMatrices, "getMatrices");
bool convergence = false;
gmshDdm->userDefinedParameter(convergence, "convergence");
std::string fileName = "none";
gmshDdm->userDefinedParameter(fileName, "file");
bool saveNeumannTraces = false;
gmshDdm->userDefinedParameter(saveNeumannTraces, "saveNeumannTraces");
bool spy = false;
gmshDdm->userDefinedParameter(spy, "spy");
// mesh
monoDomain(lc, pmlSize, R, meshOrder);
// source
gmshfem::analytics::AnalyticalFunction< gmshfem::analytics::helmholtz2D::ScatteringByASoftCylinder< std::complex< double > > > fAnalytic(k, R, 1., 1., 2 * k);
gmshfem::msg::info << "Running 'monoDomain2D'" << gmshfem::msg::endl;
gmshfem::msg::info << "Parameters:" << gmshfem::msg::endl;
gmshfem::msg::info << " * physics:" << gmshfem::msg::endl;
gmshfem::msg::info << " - k: " << k << " (" << k/pi << "*pi" << ")" << gmshfem::msg::endl;
gmshfem::msg::info << " - R: " << R << gmshfem::msg::endl;
gmshfem::msg::info << " * mesh:" << gmshfem::msg::endl;
gmshfem::msg::info << " - lc: " << lc << " (eta_h = " << (2*pi/k)/lc << ")" << gmshfem::msg::endl;
gmshfem::msg::info << " - meshOrder: " << meshOrder << gmshfem::msg::endl;
gmshfem::msg::info << " * modeling:" << gmshfem::msg::endl;
gmshfem::msg::info << " - boundary: " << boundary << gmshfem::msg::endl;
gmshfem::msg::info << " - N: " << N << gmshfem::msg::endl;
if(boundary == "pml") {
gmshfem::msg::info << " - pmlSize: " << pmlSize << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlMethod: " << pmlMethod << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlType: " << pmlType << gmshfem::msg::endl;
}
else if(boundary == "habc") {
gmshfem::msg::info << " - thetaPade: " << thetaPade << " (" << thetaPade/pi << "*pi" << ")" << gmshfem::msg::endl;
}
gmshfem::msg::info << " - gauss: " << gauss << gmshfem::msg::endl;
gmshfem::msg::info << " - fieldOrder: " << fieldOrder << gmshfem::msg::endl;
gmshfem::msg::info << " - neumannOrder: " << neumannOrder << gmshfem::msg::endl;
gmshfem::msg::info << " - stab: " << stab << " * lc (= " << stab * lc << ")" << gmshfem::msg::endl;
gmshfem::domain::Domain omega("omega");
gmshfem::domain::Domain gamma("gamma");
gmshfem::domain::Domain sigmas;
std::vector< gmshfem::domain::Domain > pml(4);
std::vector< gmshfem::domain::Domain > pmlCorner(4);
std::vector< gmshfem::domain::Domain > sigma(4);
std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlBnd(4);
std::vector< gmshfem::domain::Domain > corner(4);
std::vector< std::string > dir {"E", "N", "W", "S"};
for(unsigned int i = 0; i < 4; ++i) {
pml[i] = gmshfem::domain::Domain("pml" + dir[i]);
pmlCorner[i] = gmshfem::domain::Domain("pmlCorner"+ dir[i] + dir[(i+1)%4]);
sigma[i] = gmshfem::domain::Domain("sigma" + dir[i]);
pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i] + "_second");
pmlBnd[i].second = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "_first");
corner[i] = gmshfem::domain::Domain("corner" + dir[i] + dir[(i+1)%4]);
sigmas |= sigma[i];
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > u("u", omega | sigmas | gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
u.addConstraint(gamma, fAnalytic);
gmshfem::problem::Formulation< std::complex< double > > formulation("HelmholtzMonoDomain");
formulation.integral(-grad(dof(u)), grad(tf(u)), omega, gauss);
formulation.integral(k * k * dof(u), tf(u), omega, gauss);
SubproblemDomains domains;
domains.setOmega(omega);
domains.setSigma(sigma);
domains.setPml(pml);
domains.setPmlBnd(pmlBnd);
domains.setPmlCorner(pmlCorner);
domains.setCorner(corner);
SubproblemParameters parameters;
parameters.setGauss(gauss);
parameters.setKappa(k);
parameters.setNeumannOrder(neumannOrder);
parameters.setFieldOrder(fieldOrder);
parameters.setStab(stab * lc);
std::string bnd = boundary;
if(boundary == "habc") {
bnd += "_" + std::to_string(N) + "_" + std::to_string(thetaPade);
}
else if(boundary == "pml") {
if(pmlMethod == "continuous") {
bnd += "Continuous_" + std::to_string(pmlSize) + "_" + pmlType;
}
else if(pmlMethod == "discontinuous") {
bnd += "Discontinuous_" + std::to_string(pmlSize) + "_" + pmlType;
}
else if(pmlMethod == "withoutMultiplier") {
bnd += "WithoutMultiplier_" + std::to_string(pmlSize) + "_" + pmlType;
}
}
Subproblem subproblem(formulation, u.name(), domains, parameters, bnd, bnd, bnd, bnd);
subproblem.writeFormulation();
if(spy) {
gmshfem::common::Options::instance()->dofsSortAlgorithm = gmshfem::problem::DofsSort::Algorithm::None;
}
formulation.pre();
gmshfem::common::Timer assembleTime = formulation.assemble();
gmshfem::common::Timer solveTime = formulation.solve();
// POST
double errorL2 = 0.;
if(computeError) {
const double num = std::real(gmshfem::post::integrate(pow(abs(u - fAnalytic), 2), omega, gauss));
const double den = std::real(gmshfem::post::integrate(pow(abs(fAnalytic), 2), omega, gauss));
errorL2 = std::sqrt(num/den);
gmshfem::msg::info << "Error L2 = " << errorL2 << gmshfem::msg::endl;
gmshfem::post::save(u - fAnalytic, omega, "e");
}
if(computeInterfaceError) {
if(boundary == "pml") {
double numInterface = 0.;
double denInterface = 0.;
for(unsigned int b = 0; b < 4; ++b) {
const Boundary *boundary = subproblem.getBoundary(b);
if(auto *bnd = dynamic_cast< const PmlContinuous * >(boundary)) {
numInterface += std::real(gmshfem::post::integrate(pow(abs(u - *bnd->getUPml()), 2), sigma[b], gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(u), 2), sigma[b], gauss));
}
else if(auto *bnd = dynamic_cast< const PmlDiscontinuous * >(boundary)) {
numInterface += std::real(gmshfem::post::integrate(pow(abs(u - *bnd->getUPml()), 2), sigma[b], gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(u), 2), sigma[b], gauss));
}
else if(auto *bnd = dynamic_cast< const PmlWithoutMultiplier * >(boundary)) {
numInterface += std::real(gmshfem::post::integrate(pow(abs(u - *bnd->getUPml()), 2), sigma[b], gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(u), 2), sigma[b], gauss));
}
}
for(unsigned int c = 0; c < 4; ++c) {
const Corner *corner = subproblem.getCorner(c);
if(auto *cr = dynamic_cast< const PmlContinuous_PmlContinuous * >(corner)) {
numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].first, gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml()), 2), pmlBnd[c].first, gauss));
numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].second, gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml()), 2), pmlBnd[c].second, gauss));
}
else if(auto *cr = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous * >(corner)) {
numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].first, gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml()), 2), pmlBnd[c].first, gauss));
numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].second, gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml()), 2), pmlBnd[c].second, gauss));
}
else if(auto *cr = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier * >(corner)) {
numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].first, gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml()), 2), pmlBnd[c].first, gauss));
numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].second, gauss));
denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml()), 2), pmlBnd[c].second, gauss));
}
}
gmshfem::msg::info << "Interface L2 error = " << std::sqrt(numInterface/denInterface) << gmshfem::msg::endl;
if(fileName != "none") {
gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
file << (2*pi/k)/lc << lc << 1./lc << std::sqrt(numInterface/denInterface) << gmshfem::csv::endl;
file.close();
}
}
else {
gmshfem::msg::error << "Unable to define a boundary error with 'boundary = " << boundary << "'" << gmshfem::msg::endl;
}
}
if(getMatrices) {
// gmshfem::algebra::MatrixCRS< std::complex< double > > U;
// gmshfem::algebra::MatrixCRS< std::complex< double > > L;
// gmshfem::algebra::MatrixCRS< std::complex< double > > M;
// gmshfem::algebra::MatrixCRS< std::complex< double > > C;
gmshfem::algebra::MatrixCRS< std::complex< double > > A;
//
// std::vector< const gmshfem::field::FieldInterface< std::complex< double > > * > fieldsU;
// fieldsU.push_back(&u);
// for(unsigned int i = 0; i < 4; ++i) {
// fieldsU.push_back(subproblem.getU(i));
// }
// for(unsigned int i = 0; i < 4; ++i) {
// fieldsU.push_back(subproblem.getUC(i));
// }
//
// std::vector< const gmshfem::field::FieldInterface< std::complex< double > > * > fieldsL;
// for(unsigned int i = 0; i < 4; ++i) {
// fieldsL.push_back(subproblem.getV(i));
// }
// if(boundary == "pml") {
// for(unsigned int i = 0; i < 4; ++i) {
// fieldsL.push_back(subproblem.getVC(i, 0));
// fieldsL.push_back(subproblem.getVC(i, 1));
// }
// }
// std::vector< const gmshfem::field::FieldInterface< std::complex< double > > * > fieldsC;
// if(boundary == "pml") {
// if(form == "form0") {
// for(unsigned int i = 0; i < 4; ++i) {
// fieldsC.push_back(static_cast< NeumannTrace< gmshfem::field::Form::Form0 > * >(neumannTrace)->get("vC_" + orientation[(i-1)%4] + "_" + orientation[i]));
// }
// }
// }
// else if (boundary == "habc") {
// for(unsigned int i = 0; i < 4; ++i) {
// for(unsigned int n = 0; n < N; ++n) {
// fieldsC.push_back(static_cast< NeumannTrace< gmshfem::field::Form::Form0 > * >(neumannTrace)->get("vC_" + orientation[i] + "_" + std::to_string(n) + "_first"));
// fieldsC.push_back(static_cast< NeumannTrace< gmshfem::field::Form::Form0 > * >(neumannTrace)->get("vC_" + orientation[(i-1)%4] + "_" + std::to_string(n) + "_second"));
// }
// }
// }
// formulation.getLHSBlock(U, fieldsU, fieldsU);
// formulation.getLHSBlock(L, fieldsL, fieldsU);
// formulation.getLHSBlock(C, fieldsC, fieldsL);
// formulation.getLHSBlock(M, fieldsL, fieldsL);
formulation.getLHS(A);
//
// U.save("U");
// L.save("L");
// C.save("C");
// M.save("M");
A.save("A");
}
gmshfem::post::save(u, omega, "u");
if(convergence) {
gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
// file << lc << (2*pi/k)/lc << errorL2 << assembleTime << solveTime << formulation.getTotalNumberOfDof() << gmshfem::csv::endl;
file << k << (2*pi/k)/lc << errorL2 << assembleTime << solveTime << gmshfem::csv::endl;
}
if(saveNeumannTraces) {
for(unsigned int b = 0; b < 4; ++b) {
const Boundary *boundary = subproblem.getBoundary(b);
if(auto *bnd = dynamic_cast< const Sommerfeld * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const HABC * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const PmlContinuous * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const PmlDiscontinuous * >(boundary)) {
auto n = gmshfem::function::normal< std::complex< double > >();
gmshfem::post::save(*(bnd->getV()) * n, bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const PmlWithoutMultiplier * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
}
for(unsigned int c = 0; c < 4; ++c) {
const Corner *corner = subproblem.getCorner(c);
if(auto *cr = dynamic_cast< const PmlContinuous_PmlContinuous * >(corner)) {
gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
}
else if(auto *cr = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous * >(corner)) {
auto n = gmshfem::function::normal< std::complex< double > >();
gmshfem::post::save(*(cr->getV(0)) * n, cr->getV(0)->domain(), cr->getV(0)->name());
gmshfem::post::save(*(cr->getV(1)) * n, cr->getV(1)->domain(), cr->getV(1)->name());
}
else if(auto *cr = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier * >(corner)) {
gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
}
}
}
if(spy) {
gmshfem::algebra::MatrixCRSFast< std::complex< double > > A;
formulation.getLHS(A);
A.saveSpyPlot("spy", 8);
}
}
}
namespace D3 {
void monoDomain()
{
gmshddm::common::GmshDdm *gmshDdm = gmshddm::common::GmshDdm::currentInstance();
// ************************
// P H Y S I C S
// ************************
double pi = 3.141592653589793238462643383279;
double k = 2. * pi;
gmshDdm->userDefinedParameter(k, "k");
double R = 0.5;
gmshDdm->userDefinedParameter(R, "R");
// ************************
// M E S H
// ************************
double lc = 0.06666666666666; // Marmousi = 20; other = 0.06666666666666
gmshDdm->userDefinedParameter(lc, "lc");
int meshOrder = 1;
gmshDdm->userDefinedParameter(meshOrder, "meshOrder");
// ************************
// M O D E L I N G
// ************************
std::string boundary = "sommerfeld"; // sommerfeld, pml, habc
gmshDdm->userDefinedParameter(boundary, "boundary");
std::string pmlMethod = "continuous"; // continuous, discontinuous
gmshDdm->userDefinedParameter(pmlMethod, "pmlMethod");
std::string pmlType = "hs"; // hs, h, q
gmshDdm->userDefinedParameter(pmlType, "pmlType");
double N = 6;
gmshDdm->userDefinedParameter(N, "N");
double pmlSize = N * lc;
if(pmlSize == 0.) {
pmlSize = 0.3;
}
std::string gauss = "Gauss10";
gmshDdm->userDefinedParameter(gauss, "gauss");
int fieldOrder = 1;
gmshDdm->userDefinedParameter(fieldOrder, "fieldOrder");
int neumannOrder = 1;
gmshDdm->userDefinedParameter(neumannOrder, "neumannOrder");
double stab = 0.11;
gmshDdm->userDefinedParameter(stab, "stab");
double thetaPade = 0.;
gmshDdm->userDefinedParameter(thetaPade, "thetaPade");
thetaPade *= pi;
// ************************
// P O S T
// ************************
bool onPlane = false;
gmshDdm->userDefinedParameter(onPlane, "onPlane");
bool computeError = false;
gmshDdm->userDefinedParameter(computeError, "error");
bool computeInterfaceError = false;
gmshDdm->userDefinedParameter(computeInterfaceError, "interfaceError");
bool getMatrices = false;
gmshDdm->userDefinedParameter(getMatrices, "getMatrices");
bool saveNeumannTraces = false;
gmshDdm->userDefinedParameter(saveNeumannTraces, "saveNeumannTraces");
// mesh
monoDomain(lc, pmlSize, R, meshOrder);
// source
gmshfem::analytics::AnalyticalFunction< gmshfem::analytics::helmholtz3D::ScatteringByASoftSphere< std::complex< double > > > fAnalytic(k, R, 1., 1., 1., 2 * k);
gmshfem::msg::info << "Running 'monoDomain3D'" << gmshfem::msg::endl;
gmshfem::msg::info << "Parameters:" << gmshfem::msg::endl;
gmshfem::msg::info << " * physics:" << gmshfem::msg::endl;
gmshfem::msg::info << " - k: " << k << " (" << k/pi << "*pi" << ")" << gmshfem::msg::endl;
gmshfem::msg::info << " - R: " << R << gmshfem::msg::endl;
gmshfem::msg::info << " * mesh:" << gmshfem::msg::endl;
gmshfem::msg::info << " - lc: " << lc << " (eta_h = " << (2*pi/k)/lc << ")" << gmshfem::msg::endl;
gmshfem::msg::info << " - meshOrder: " << meshOrder << gmshfem::msg::endl;
gmshfem::msg::info << " * modeling:" << gmshfem::msg::endl;
gmshfem::msg::info << " - boundary: " << boundary << gmshfem::msg::endl;
gmshfem::msg::info << " - N: " << N << gmshfem::msg::endl;
if(boundary == "pml") {
gmshfem::msg::info << " - pmlSize: " << pmlSize << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlMethod: " << pmlMethod << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlType: " << pmlType << gmshfem::msg::endl;
}
else if(boundary == "habc") {
gmshfem::msg::info << " - thetaPade: " << thetaPade << " (" << thetaPade/pi << "*pi" << ")" << gmshfem::msg::endl;
}
gmshfem::msg::info << " - gauss: " << gauss << gmshfem::msg::endl;
gmshfem::msg::info << " - fieldOrder: " << fieldOrder << gmshfem::msg::endl;
gmshfem::msg::info << " - neumannOrder: " << neumannOrder << gmshfem::msg::endl;
gmshfem::msg::info << " - stab: " << stab << " * lc (= " << stab * lc << ")" << gmshfem::msg::endl;
gmshfem::domain::Domain omega("omega");
gmshfem::domain::Domain gamma("gamma");
gmshfem::domain::Domain sigmas;
std::vector< gmshfem::domain::Domain > pml(6);
std::vector< gmshfem::domain::Domain > pmlEdge(12);
std::vector< gmshfem::domain::Domain > pmlCorner(8);
std::vector< gmshfem::domain::Domain > sigma(6);
std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlBnd(12);
std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlEdgeBnd(8);
std::vector< gmshfem::domain::Domain > edge(12);
std::vector< gmshfem::domain::Domain > corner(8);
std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > cornerEdge(8);
std::vector< std::string > dir {"E", "N", "W", "S", "D", "U"};
// face
for(unsigned int i = 0; i < 6; ++i) {
pml[i] = gmshfem::domain::Domain("pml" + dir[i]);
sigma[i] = gmshfem::domain::Domain("sigma" + dir[i]);
sigmas |= sigma[i];
}
// edge
for(unsigned int i = 0; i < 12; ++i) {
std::vector< std::string > count {"first", "second", "third", "fourth"};
if(i < 4) {
pmlEdge[i] = gmshfem::domain::Domain("pml" + dir[i%4] + "D");
pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "_first");
pmlBnd[i].second = gmshfem::domain::Domain("pmlBndD_" + count[i%4]);
edge[i] = gmshfem::domain::Domain("edge" + dir[i%4] + "D");
}
else if(i >= 4 && i < 8) {
pmlEdge[i] = gmshfem::domain::Domain("pml" + dir[i%4] + dir[(i+1)%4]);
pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "_second");
pmlBnd[i].second = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "_fourth");
edge[i] = gmshfem::domain::Domain("edge" + dir[i%4] + dir[(i+1)%4]);
}
else {
pmlEdge[i] = gmshfem::domain::Domain("pml" + dir[i%4] + "U");
pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "_third");
pmlBnd[i].second = gmshfem::domain::Domain("pmlBndU_" + count[i%4]);
edge[i] = gmshfem::domain::Domain("edge" + dir[i%4] + "U");
}
}
// corner
for(unsigned int i = 0; i < 8; ++i) {
if(i < 4) {
pmlCorner[i] = gmshfem::domain::Domain("pml" + dir[i%4] + dir[(i+1)%4] + "D");
std::get<0>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "D_second");
std::get<1>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "D_first");
std::get<2>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + dir[(i+1)%4] + "_first");
corner[i] = gmshfem::domain::Domain("corner" + dir[i%4] + dir[(i+1)%4] + "D");
std::get<0>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "D_first");
std::get<1>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "D_second");
std::get<2>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "D_third");
}
else {
pmlCorner[i] = gmshfem::domain::Domain("pml" + dir[i%4] + dir[(i+1)%4] + "U");
std::get<0>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "U_second");
std::get<1>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "U_first");
std::get<2>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + dir[(i+1)%4] + "_second");
corner[i] = gmshfem::domain::Domain("corner" + dir[i%4] + dir[(i+1)%4] + "U");
std::get<0>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "U_first");
std::get<1>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "U_second");
std::get<2>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "U_third");
}
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > u("u", omega | sigmas | gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
u.addConstraint(gamma, fAnalytic);
gmshfem::problem::Formulation< std::complex< double > > formulation("HelmholtzMonoDomain");
formulation.integral(-grad(dof(u)), grad(tf(u)), omega, gauss);
formulation.integral(k * k * dof(u), tf(u), omega, gauss);
SubproblemDomains domains;
domains.setOmega(omega);
domains.setSigma(sigma);
domains.setPml(pml);
domains.setPmlBnd(pmlBnd);
domains.setPmlEdge(pmlEdge);
domains.setEdge(edge);
domains.setPmlEdgeBnd(pmlEdgeBnd);
domains.setPmlCorner(pmlCorner);
domains.setCorner(corner);
domains.setCornerEdge(cornerEdge);
SubproblemParameters parameters;
parameters.setGauss(gauss);
parameters.setKappa(k);
parameters.setNeumannOrder(neumannOrder);
parameters.setFieldOrder(fieldOrder);
parameters.setStab(stab * lc);
std::string bnd = boundary;
if(boundary == "habc") {
bnd += "_" + std::to_string(N) + "_" + std::to_string(thetaPade);
}
else if(boundary == "pml") {
if(pmlMethod == "continuous") {
bnd += "Continuous_" + std::to_string(pmlSize) + "_" + pmlType;
}
else if(pmlMethod == "discontinuous") {
bnd += "Discontinuous_" + std::to_string(pmlSize) + "_" + pmlType;
}
}
Subproblem subproblem(formulation, u.name(), domains, parameters, bnd, bnd, bnd, bnd, bnd, bnd);
subproblem.writeFormulation();
formulation.pre();
formulation.assemble();
formulation.solve();
// POST
const double eps = 0.01;
gmshfem::post::Plane planeXY("planeXY", -1.+eps, -1.+eps, 0., 2.-2.*eps, 0., 0., 0., 2.-2.*eps, 0., 100., 100.);
gmshfem::post::Plane planeYZ("planeYZ", -1.+eps, -1.+eps, -1.+eps, 2.-2.*eps, 2.-2.*eps, 0., 0., 0., 2.-2.*eps, 100., 100.);
gmshfem::post::Plane planeZX("planeZX", -1.+eps, 1.-eps, -1.+eps, 2.-2.*eps, -2.+2.*eps, 0., 0., 0., 2.-2.*eps, 100., 100.);
gmshfem::post::Sphere scatter("scatter", 0., 0., 0., R+eps, 100. * R);
if(computeError) {
const double num = std::real(gmshfem::post::integrate(pow(abs(u - fAnalytic), 2), omega, gauss));
const double den = std::real(gmshfem::post::integrate(pow(abs(fAnalytic), 2), omega, gauss));
gmshfem::msg::info << "Error L2 = " << std::sqrt(num/den) << gmshfem::msg::endl;
if(onPlane) {
gmshfem::post::save(gmshfem::function::heaviside(gmshfem::function::r3d< std::complex< double > >() - R) * (u - fAnalytic), planeXY, "e_XY", "pos");
gmshfem::post::save(gmshfem::function::heaviside(gmshfem::function::r3d< std::complex< double > >() - R) * (u - fAnalytic), planeYZ, "e_YZ", "pos");
gmshfem::post::save(gmshfem::function::heaviside(gmshfem::function::r3d< std::complex< double > >() - R) * (u - fAnalytic), planeZX, "e_ZX", "pos");
gmshfem::post::save(u - fAnalytic, scatter, "e_SCAT", "pos");
}
else {
gmshfem::post::save(gmshfem::function::norm(u - fAnalytic), omega, "e");
}
}
if(onPlane) {
gmshfem::post::save(u, planeXY, "u_XY", "pos");
gmshfem::post::save(u, planeYZ, "u_YZ", "pos");
gmshfem::post::save(u, planeZX, "u_ZX", "pos");
gmshfem::post::save(u, scatter, "u_SCAT", "pos");
}
else {
gmshfem::post::save(u, omega, "u");
}
if(getMatrices) {
gmshfem::algebra::MatrixCRS< std::complex< double > > A;
formulation.getLHS(A);
A.save("A");
}
if(saveNeumannTraces) {
for(unsigned int b = 0; b < 6; ++b) {
const Boundary *boundary = subproblem.getBoundary(b);
if(auto *bnd = dynamic_cast< const Sommerfeld * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const HABC * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const PmlContinuous * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const PmlDiscontinuous * >(boundary)) {
auto n = gmshfem::function::normal< std::complex< double > >();
gmshfem::post::save(*(bnd->getV()) * n, bnd->getV()->domain(), bnd->getV()->name());
}
else if(auto *bnd = dynamic_cast< const PmlWithoutMultiplier * >(boundary)) {
gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
}
}
for(unsigned int e = 0; e < 12; ++e) {
const Edge *edge = subproblem.getEdge(e);
if(auto *eg = dynamic_cast< const PmlContinuous_PmlContinuous * >(edge)) {
gmshfem::post::save(*(eg->getV(0)), eg->getV(0)->domain(), eg->getV(0)->name());
gmshfem::post::save(*(eg->getV(1)), eg->getV(1)->domain(), eg->getV(1)->name());
}
// else if(auto *eg = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous * >(edge)) {
// auto n = gmshfem::function::normal< std::complex< double > >();
// gmshfem::post::save(*(eg->getV(0)) * n, eg->getV(0)->domain(), eg->getV(0)->name());
// gmshfem::post::save(*(eg->getV(1)) * n, eg->getV(1)->domain(), eg->getV(1)->name());
// }
// else if(auto *eg = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier * >(boundary)) {
// gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
// }
}
for(unsigned int c = 0; c < 8; ++c) {
const Corner *corner = subproblem.getCorner(c);
if(auto *cr = dynamic_cast< const PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
gmshfem::post::save(*(cr->getV(2)), cr->getV(2)->domain(), cr->getV(2)->name());
}
// else if(auto *cr = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous_PmlDiscontinuous * >(corner)) {
// auto n = gmshfem::function::normal< std::complex< double > >();
// gmshfem::post::save(*(cr->getV(0)) * n, cr->getV(0)->domain(), cr->getV(0)->name());
// gmshfem::post::save(*(cr->getV(1)) * n, cr->getV(1)->domain(), cr->getV(1)->name());
// gmshfem::post::save(*(cr->getV(2)) * n, cr->getV(2)->domain(), cr->getV(2)->name());
// }
// else if(auto *cr = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier_PmlWithoutMultiplier * >(corner)) {
// gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
// gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
// }
}
}
}
}