Select Git revision
delaunay_refinement.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
ddm3D.cpp 102.54 KiB
#include "ddm3D.h"
#include "mesh.h"
#include "SubproblemDomains.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>
#include <vector>
using gmshfem::equation::dof;
using gmshfem::equation::tf;
using gmshfem::function::operator-;
using gmshfem::function::operator*;
using gmshfem::function::abs;
namespace D3 {
void ddm()
{
gmshddm::common::GmshDdm *gmshDdm = gmshddm::common::GmshDdm::currentInstance();
// ************************
// P H Y S I C S
// ************************
double pi = 3.14159265359;
double k = 2. * pi;
gmshDdm->userDefinedParameter(k, "k");
double R = 0.5;
gmshDdm->userDefinedParameter(R, "R");
std::string benchmark = "scattering"; // scattering
gmshDdm->userDefinedParameter(benchmark, "benchmark");
// ************************
// M E S H
// ************************
double lc = 0.1; // 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 boundaryExt = "sommerfeld"; // sommerfeld, pml, habc
gmshDdm->userDefinedParameter(boundaryExt, "boundaryExt");
std::string pmlMethod = "continuous"; // continuous, discontinuous
gmshDdm->userDefinedParameter(pmlMethod, "pmlMethod");
std::string pmlMethodExt = "continuous"; // continuous, discontinuous
gmshDdm->userDefinedParameter(pmlMethodExt, "pmlMethodExt");
std::string pmlType = "hs"; // hs, h, q
gmshDdm->userDefinedParameter(pmlType, "pmlType");
std::string pmlTypeExt = "hs"; // hs, h, q
gmshDdm->userDefinedParameter(pmlTypeExt, "pmlTypeExt");
unsigned int nDomX = 2, nDomY = 2, nDomZ = 2;
gmshDdm->userDefinedParameter(nDomX, "nDomX");
gmshDdm->userDefinedParameter(nDomY, "nDomY");
gmshDdm->userDefinedParameter(nDomZ, "nDomZ");
unsigned int iterMax = 1000;
gmshDdm->userDefinedParameter(iterMax, "iterMax");
double res = 1e-6;
gmshDdm->userDefinedParameter(res, "res");
double sizeX = 2., sizeY = 2., sizeZ = 2.;
gmshDdm->userDefinedParameter(sizeX, "sizeX");
gmshDdm->userDefinedParameter(sizeY, "sizeY");
gmshDdm->userDefinedParameter(sizeZ, "sizeZ");
unsigned int N = 6;
gmshDdm->userDefinedParameter(N, "N");
double pmlSize = N * lc;
if(pmlSize == 0.) {
pmlSize = 0.3;
}
unsigned int NExt = N;
gmshDdm->userDefinedParameter(NExt, "NExt");
double pmlSizeExt = NExt * lc;
if(pmlSizeExt == 0.) {
pmlSizeExt = 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 monoDomainError = false;
gmshDdm->userDefinedParameter(monoDomainError, "monoDomainError");
std::string fileName = "none";
gmshDdm->userDefinedParameter(fileName, "file");
bool saveEMono = false;
gmshDdm->userDefinedParameter(saveEMono, "saveEMono");
bool saveUMono = false;
gmshDdm->userDefinedParameter(saveUMono, "saveUMono");
bool saveU = false;
gmshDdm->userDefinedParameter(saveU, "saveU");
bool wavenumberPlot = false;
gmshDdm->userDefinedParameter(wavenumberPlot, "wavenumber");
bool meshPlot = false;
gmshDdm->userDefinedParameter(meshPlot, "mesh");
bool saveMesh = false;
gmshDdm->userDefinedParameter(saveMesh, "saveMesh");
if(benchmark == "scattering") {
checkerboard(nDomX, nDomY, nDomZ, sizeX, sizeY, sizeZ, R, lc, (boundary == "pml"), pmlSize, (boundaryExt == "pml"), pmlSizeExt, meshOrder, 0, 0, 0);
}
if(saveMesh) {
gmsh::write("mesh.msh");
}
gmshfem::msg::info << "Running 'ddm3D'" << 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 << " - grid: (" << nDomX << ", " << nDomY << ", " << nDomZ << ")" << gmshfem::msg::endl;
gmshfem::msg::info << " - iterMax: " << iterMax << gmshfem::msg::endl;
gmshfem::msg::info << " - res: " << res << gmshfem::msg::endl;
gmshfem::msg::info << " - subdomain size: (" << sizeX << ", " << sizeY << ", " << sizeZ << ")" << gmshfem::msg::endl;
gmshfem::msg::info << " - boundary: " << boundary << gmshfem::msg::endl;
gmshfem::msg::info << " - boundaryExt: " << boundaryExt << gmshfem::msg::endl;
if(boundary == "pml") {
gmshfem::msg::info << " - N: " << N << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlSize: " << pmlSize << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlType: " << pmlType << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlMethod: " << pmlMethod << gmshfem::msg::endl;
}
else if(boundary == "habc") {
gmshfem::msg::info << " - N: " << N << gmshfem::msg::endl;
gmshfem::msg::info << " - thetaPade: " << thetaPade << " (" << thetaPade/pi << "*pi" << ")" << gmshfem::msg::endl;
}
if(boundaryExt == "pml") {
gmshfem::msg::info << " - NExt: " << NExt << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlSizeExt: " << pmlSizeExt << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlTypeExt: " << pmlTypeExt << gmshfem::msg::endl;
gmshfem::msg::info << " - pmlMethodExt: " << pmlMethodExt << gmshfem::msg::endl;
}
else if(boundaryExt == "habc") {
gmshfem::msg::info << " - NExt: " << NExt << 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 << gmshfem::msg::endl;
// source
gmshfem::analytics::AnalyticalFunction< gmshfem::analytics::helmholtz3D::ScatteringByASoftSphere< std::complex< double > > > fAnalytic(k, R, sizeX/2., sizeY/2., sizeZ/2., 2 * k);
const unsigned int nDom = nDomX * nDomY * nDomZ;
// Define domain
gmshddm::domain::Subdomain omega(nDom);
gmshfem::domain::Domain gamma;
if(benchmark == "scattering") {
gamma = gmshfem::domain::Domain("gamma");
}
std::vector< gmshddm::domain::Subdomain > pml(6, nDom);
std::vector< gmshddm::domain::Subdomain > pmlEdge(12, nDom);
std::vector< gmshddm::domain::Subdomain > pmlCorner(8, nDom);
std::vector< gmshddm::domain::Subdomain > sigma(6, nDom);
std::vector< gmshddm::domain::Interface > sigmaInterface(6, nDom);
std::vector< std::pair< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > pmlBnd(12, std::make_pair(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)) );
std::vector< std::pair< gmshddm::domain::Interface, gmshddm::domain::Interface > > pmlBndInterface(12, std::make_pair(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
std::vector< std::tuple< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > pmlEdgeBnd(8, std::make_tuple(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)) );
std::vector< std::tuple< gmshddm::domain::Interface, gmshddm::domain::Interface, gmshddm::domain::Interface > > pmlEdgeBndInterface(8, std::make_tuple(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
std::vector< gmshddm::domain::Subdomain > edge(12, nDom);
std::vector< std::pair< gmshddm::domain::Interface, gmshddm::domain::Interface > > edgeInterface(12, std::make_pair(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
std::vector< gmshddm::domain::Subdomain > corner(8, nDom);
std::vector< std::tuple< gmshddm::domain::Interface, gmshddm::domain::Interface, gmshddm::domain::Interface > > cornerInterface(8, std::make_tuple(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
std::vector< std::tuple< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > cornerEdge(8, std::make_tuple(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)));
// Define topology
std::vector< std::vector< unsigned int > > topology(nDom);
std::vector< std::string > dir {"E", "N", "W", "S", "D", "U"};
for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) {
for(unsigned int j = 0; j < static_cast< unsigned int >(nDomY); ++j) {
for(unsigned int k = 0; k < static_cast< unsigned int >(nDomZ); ++k) {
unsigned int index = i * nDomZ * nDomY + j * nDomZ + k;
const std::string subdomainTag = "_" + std::to_string(i) + "_" + std::to_string(j) + "_" + std::to_string(k);
omega(index) = gmshfem::domain::Domain("omega" + subdomainTag);
for(unsigned int b = 0; b < 6; ++b) {
if(boundary == "pml" || boundaryExt == "pml") {
pml[b](index) = gmshfem::domain::Domain("pml" + dir[b] + subdomainTag);
}
sigma[b](index) = gmshfem::domain::Domain("sigma" + dir[b] + subdomainTag);
}
for(unsigned int e = 0; e < 12; ++e) {
std::vector< std::string > count {"first", "second", "third", "fourth"};
if(e < 4) {
if(boundary == "pml" || boundaryExt == "pml") {
pmlEdge[e](index) = gmshfem::domain::Domain("pml" + dir[e%4] + "D" + subdomainTag);
pmlBnd[e].first(index) = gmshfem::domain::Domain("pmlBnd" + dir[e%4] + "_first" + subdomainTag);
pmlBnd[e].second(index) = gmshfem::domain::Domain("pmlBndD_" + count[e%4] + subdomainTag);
}
edge[e](index) = gmshfem::domain::Domain("edge" + dir[e%4] + "D" + subdomainTag);
}
else if(e >= 4 && e < 8) {
if(boundary == "pml" || boundaryExt == "pml") {
pmlEdge[e](index) = gmshfem::domain::Domain("pml" + dir[e%4] + dir[(e+1)%4] + subdomainTag);
pmlBnd[e].first(index) = gmshfem::domain::Domain("pmlBnd" + dir[e%4] + "_second" + subdomainTag);
pmlBnd[e].second(index) = gmshfem::domain::Domain("pmlBnd" + dir[(e+1)%4] + "_fourth" + subdomainTag);
}
edge[e](index) = gmshfem::domain::Domain("edge" + dir[e%4] + dir[(e+1)%4] + subdomainTag);
}
else {
if(boundary == "pml" || boundaryExt == "pml") {
pmlEdge[e](index) = gmshfem::domain::Domain("pml" + dir[e%4] + "U" + subdomainTag);
pmlBnd[e].first(index) = gmshfem::domain::Domain("pmlBnd" + dir[e%4] + "_third" + subdomainTag);
pmlBnd[e].second(index) = gmshfem::domain::Domain("pmlBndU_" + count[e%4] + subdomainTag);
}
edge[e](index) = gmshfem::domain::Domain("edge" + dir[e%4] + "U" + subdomainTag);
}
}
for(unsigned int c = 0; c < 8; ++c) {
if(c < 4) {
if(boundary == "pml" || boundaryExt == "pml") {
pmlCorner[c](index) = gmshfem::domain::Domain("pml" + dir[c%4] + dir[(c+1)%4] + "D" + subdomainTag);
std::get<0>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + "D_second" + subdomainTag);
std::get<1>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[(c+1)%4] + "D_first" + subdomainTag);
std::get<2>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + dir[(c+1)%4] + "_first" + subdomainTag);
std::get<0>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "D_first" + subdomainTag);
std::get<1>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "D_second" + subdomainTag);
std::get<2>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "D_third" + subdomainTag);
}
corner[c](index) = gmshfem::domain::Domain("corner" + dir[c%4] + dir[(c+1)%4] + "D" + subdomainTag);
}
else {
if(boundary == "pml" || boundaryExt == "pml") {
pmlCorner[c](index) = gmshfem::domain::Domain("pml" + dir[c%4] + dir[(c+1)%4] + "U" + subdomainTag);
std::get<0>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + "U_second" + subdomainTag);
std::get<1>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[(c+1)%4] + "U_first" + subdomainTag);
std::get<2>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + dir[(c+1)%4] + "_second" + subdomainTag);
std::get<0>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "U_first" + subdomainTag);
std::get<1>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "U_second" + subdomainTag);
std::get<2>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "U_third" + subdomainTag);
}
corner[c](index) = gmshfem::domain::Domain("corner" + dir[c%4] + dir[(c+1)%4] + "U" + subdomainTag);
}
}
if(i != static_cast< unsigned int >(nDomX) - 1) { // E
if(boundary == "pml" || boundaryExt == "pml") {
pmlBndInterface[0].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[0].second(index);
pmlBndInterface[4].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[4].second(index);
pmlBndInterface[8].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[8].second(index);
pmlBndInterface[7].first(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[7].first(index);
std::get<1>(pmlEdgeBndInterface[0])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[0])(index);
std::get<1>(pmlEdgeBndInterface[4])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[4])(index);
std::get<0>(pmlEdgeBndInterface[7])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[7])(index);
std::get<0>(pmlEdgeBndInterface[3])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[3])(index);
}
sigmaInterface[0](index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = sigma[0](index);
edgeInterface[0].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[0](index);
edgeInterface[4].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[4](index);
edgeInterface[8].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[8](index);
edgeInterface[7].first(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[7](index);
std::get<1>(cornerInterface[0])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[0](index);
std::get<1>(cornerInterface[4])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[4](index);
std::get<0>(cornerInterface[7])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[7](index);
std::get<0>(cornerInterface[3])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[3](index);
topology[index].push_back((i+1) * nDomY * nDomZ + j * nDomZ + k);
}
if(j != static_cast< unsigned int >(nDomY) - 1) { // N
if(boundary == "pml" || boundaryExt == "pml") {
pmlBndInterface[1].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[1].second(index);
pmlBndInterface[5].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[5].second(index);
pmlBndInterface[9].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[9].second(index);
pmlBndInterface[4].first(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[4].first(index);
std::get<1>(pmlEdgeBndInterface[1])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[1])(index);
std::get<1>(pmlEdgeBndInterface[5])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[5])(index);
std::get<0>(pmlEdgeBndInterface[4])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[4])(index);
std::get<0>(pmlEdgeBndInterface[0])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[0])(index);
}
sigmaInterface[1](index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = sigma[1](index);
edgeInterface[1].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[1](index);
edgeInterface[5].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[5](index);
edgeInterface[9].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[9](index);
edgeInterface[4].first(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[4](index);
std::get<1>(cornerInterface[1])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[1](index);
std::get<1>(cornerInterface[5])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[5](index);
std::get<0>(cornerInterface[4])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[4](index);
std::get<0>(cornerInterface[0])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[0](index);
topology[index].push_back(i * nDomY * nDomZ + (j+1) * nDomZ + k);
}
if(i != 0) { // W
if(boundary == "pml" || boundaryExt == "pml") {
pmlBndInterface[2].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[2].second(index);
pmlBndInterface[6].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[6].second(index);
pmlBndInterface[10].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[10].second(index);
pmlBndInterface[5].first(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[5].first(index);
std::get<1>(pmlEdgeBndInterface[2])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[2])(index);
std::get<1>(pmlEdgeBndInterface[6])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[6])(index);
std::get<0>(pmlEdgeBndInterface[5])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[5])(index);
std::get<0>(pmlEdgeBndInterface[1])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[1])(index);
}
sigmaInterface[2](index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = sigma[2](index);
edgeInterface[2].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[2](index);
edgeInterface[6].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[6](index);
edgeInterface[10].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[10](index);
edgeInterface[5].first(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[5](index);
std::get<1>(cornerInterface[2])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[2](index);
std::get<1>(cornerInterface[6])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[6](index);
std::get<0>(cornerInterface[5])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[5](index);
std::get<0>(cornerInterface[1])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[1](index);
topology[index].push_back((i-1) * nDomY * nDomZ + j * nDomZ + k);
}
if(j != 0) { // S
if(boundary == "pml" || boundaryExt == "pml") {
pmlBndInterface[3].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[3].second(index);
pmlBndInterface[7].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[7].second(index);
pmlBndInterface[11].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[11].second(index);
pmlBndInterface[6].first(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[6].first(index);
std::get<1>(pmlEdgeBndInterface[3])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[3])(index);
std::get<1>(pmlEdgeBndInterface[7])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[7])(index);
std::get<0>(pmlEdgeBndInterface[6])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[6])(index);
std::get<0>(pmlEdgeBndInterface[2])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[2])(index);
}
sigmaInterface[3](index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = sigma[3](index);
edgeInterface[3].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[3](index);
edgeInterface[7].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[7](index);
edgeInterface[11].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[11](index);
edgeInterface[6].first(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[6](index);
std::get<1>(cornerInterface[3])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[3](index);
std::get<1>(cornerInterface[7])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[7](index);
std::get<0>(cornerInterface[6])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[6](index);
std::get<0>(cornerInterface[2])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[2](index);
topology[index].push_back(i * nDomY * nDomZ + (j-1) * nDomZ + k);
}
if(k != 0) { // D
if(boundary == "pml" || boundaryExt == "pml") {
pmlBndInterface[0].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[0].first(index);
pmlBndInterface[1].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[1].first(index);
pmlBndInterface[2].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[2].first(index);
pmlBndInterface[3].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[3].first(index);
std::get<2>(pmlEdgeBndInterface[0])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[0])(index);
std::get<2>(pmlEdgeBndInterface[1])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[1])(index);
std::get<2>(pmlEdgeBndInterface[2])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[2])(index);
std::get<2>(pmlEdgeBndInterface[3])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[3])(index);
}
sigmaInterface[4](index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = sigma[4](index);
edgeInterface[0].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[0](index);
edgeInterface[1].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[1](index);
edgeInterface[2].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[2](index);
edgeInterface[3].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[3](index);
std::get<2>(cornerInterface[0])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[0](index);
std::get<2>(cornerInterface[1])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[1](index);
std::get<2>(cornerInterface[2])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[2](index);
std::get<2>(cornerInterface[3])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[3](index);
topology[index].push_back(i * nDomY * nDomZ + j * nDomZ + (k-1));
}
if(k != static_cast< unsigned int >(nDomZ) - 1) { // U
if(boundary == "pml" || boundaryExt == "pml") {
pmlBndInterface[8].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[8].first(index);
pmlBndInterface[9].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[9].first(index);
pmlBndInterface[10].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[10].first(index);
pmlBndInterface[11].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[11].first(index);
std::get<2>(pmlEdgeBndInterface[4])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[4])(index);
std::get<2>(pmlEdgeBndInterface[5])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[5])(index);
std::get<2>(pmlEdgeBndInterface[6])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[6])(index);
std::get<2>(pmlEdgeBndInterface[7])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[7])(index);
}
sigmaInterface[5](index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = sigma[5](index);
edgeInterface[8].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[8](index);
edgeInterface[9].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[9](index);
edgeInterface[10].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[10](index);
edgeInterface[11].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[11](index);
std::get<2>(cornerInterface[4])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[4](index);
std::get<2>(cornerInterface[5])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[5](index);
std::get<2>(cornerInterface[6])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[6](index);
std::get<2>(cornerInterface[7])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[7](index);
topology[index].push_back(i * nDomY * nDomZ + j * nDomZ + (k+1));
}
}
}
}
// Kappa definition
gmshfem::function::ScalarFunction< std::complex< double > > kappa;
if(benchmark == "scattering") {
kappa = k;
}
// Fields definition
gmshddm::field::SubdomainField< std::complex< double >, gmshfem::field::Form::Form0 > u("u", omega | sigma[0] | sigma[1] | sigma[2] | sigma[3], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
if(benchmark == "scattering") {
u(0).domain(u(0).domain() | gamma);
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambda("lambda", gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
std::vector< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > gContinuous;
// std::vector< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > * > gDiscontinuous;
for(unsigned int b = 0; b < 6; ++b) {
gContinuous.push_back(new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gC_" + dir[b], sigmaInterface[b], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
// gDiscontinuous.push_back(new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gD_" + dir[b], sigmaInterface[b], gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder));
}
std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > gEdgePmlContinuous;
// std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > * > > gEdgePmlDiscontinuous;
std::vector< std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > > gEdgeHABC;
for(unsigned int e = 0; e < 12; ++e) {
if(e < 4) {
gEdgePmlContinuous.push_back(std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "D_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "D_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
// gEdgePmlDiscontinuous.push_back(std::make_pair(
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "D_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "D_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
gEdgeHABC.push_back(std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N));
for(unsigned int n = 0; n < N; ++n) {
gEdgeHABC[e][n] = std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "D_" + std::to_string(n) + "_first", edgeInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "D_" + std::to_string(n) + "_second", edgeInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
}
}
else if(e >= 4 && e < 8) {
gEdgePmlContinuous.push_back(std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + dir[(e+1)%4] + "_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + dir[(e+1)%4] + "_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
// gEdgePmlDiscontinuous.push_back(std::make_pair(
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + dir[(e+1)%4] + "_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + dir[(e+1)%4] + "_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
gEdgeHABC.push_back(std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N));
for(unsigned int n = 0; n < N; ++n) {
gEdgeHABC[e][n] = std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + dir[(e+1)%4] + "_" + std::to_string(n) + "_first", edgeInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + dir[(e+1)%4] + "_" + std::to_string(n) + "_second", edgeInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
}
}
else {
gEdgePmlContinuous.push_back(std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "U_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "U_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
// gEdgePmlDiscontinuous.push_back(std::make_pair(
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "U_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "U_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
gEdgeHABC.push_back(std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N));
for(unsigned int n = 0; n < N; ++n) {
gEdgeHABC[e][n] = std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "U_" + std::to_string(n) + "_first", edgeInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "U_" + std::to_string(n) + "_second", edgeInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
}
}
}
std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > gCornerPmlContinuous;
// std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > * > > gCornerPmlDiscontinuous;
std::vector< std::vector< std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > > > gCornerHABC;
for(unsigned int c = 0; c < 8; ++c) {
if(c < 4) {
gCornerPmlContinuous.push_back(std::make_tuple(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "D_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "D_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "D_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
// gCornerPmlDiscontinuous.push_back(std::make_tuple(
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "D_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "D_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "D_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
gCornerHABC.push_back(std::vector< std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > >(N));
for(unsigned int n = 0; n < N; ++n) {
gCornerHABC[c][n] = std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N);
for(unsigned int m = 0; m < N; ++m) {
gCornerHABC[c][n][m] = std::make_tuple(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_first", std::get<0>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_second", std::get<1>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_third", std::get<2>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<0>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_first");
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<1>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_second");
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<2>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_third");
}
}
}
else {
gCornerPmlContinuous.push_back(std::make_tuple(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "U_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "U_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "U_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
// gCornerPmlDiscontinuous.push_back(std::make_tuple(
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "U_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "U_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
// new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "U_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
gCornerHABC.push_back(std::vector< std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > >(N));
for(unsigned int n = 0; n < N; ++n) {
gCornerHABC[c][n] = std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N);
for(unsigned int m = 0; m < N; ++m) {
gCornerHABC[c][n][m] = std::make_tuple(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_first", std::get<0>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_second", std::get<1>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_third", std::get<2>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<0>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_first");
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<1>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_second");
// gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<2>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_third");
}
}
}
}
gmshddm::problem::Formulation< std::complex< double > > formulation("HelmholtzDDMScattering", topology);
for(unsigned int b = 0; b < 6; ++b) {
formulation.addInterfaceField(*gContinuous[b]);
// formulation.addInterfaceField(*gDiscontinuous[b]);
}
for(unsigned int e = 0; e < 12; ++e) {
formulation.addInterfaceField(*gEdgePmlContinuous[e].first);
formulation.addInterfaceField(*gEdgePmlContinuous[e].second);
// formulation.addInterfaceField(*gEdgePmlDiscontinuous[e].first);
// formulation.addInterfaceField(*gEdgePmlDiscontinuous[e].second);
for(unsigned int n = 0; n < N; ++n) {
formulation.addInterfaceField(*gEdgeHABC[e][n].first);
formulation.addInterfaceField(*gEdgeHABC[e][n].second);
}
}
for(unsigned int c = 0; c < 8; ++c) {
formulation.addInterfaceField(*std::get<0>(gCornerPmlContinuous[c]));
formulation.addInterfaceField(*std::get<1>(gCornerPmlContinuous[c]));
formulation.addInterfaceField(*std::get<2>(gCornerPmlContinuous[c]));
// formulation.addInterfaceField(*std::get<0>(gCornerPmlDiscontinuous[c]));
// formulation.addInterfaceField(*std::get<1>(gCornerPmlDiscontinuous[c]));
// formulation.addInterfaceField(*std::get<2>(gCornerPmlDiscontinuous[c]));
for(unsigned int n = 0; n < N; ++n) {
for(unsigned int m = 0; m < N; ++m) {
formulation.addInterfaceField(*std::get<0>(gCornerHABC[c][n][m]));
formulation.addInterfaceField(*std::get<1>(gCornerHABC[c][n][m]));
formulation.addInterfaceField(*std::get<2>(gCornerHABC[c][n][m]));
}
}
}
std::vector< std::vector< std::vector< Subproblem * > > > subproblem(nDomX);
for(unsigned int i = 0; i < nDomX; ++i) {
subproblem[i].resize(nDomY);
for(unsigned int j = 0; j < nDomY; ++j) {
subproblem[i][j].resize(nDomZ);
for(unsigned int k = 0; k < nDomZ; ++k) {
unsigned int index = i * nDomY * nDomZ + j * nDomZ + k;
formulation(index).integral(-grad(dof(u(index))), grad(tf(u(index))), omega(index), gauss);
formulation(index).integral(kappa * kappa * dof(u(index)), tf(u(index)), omega(index), gauss);
if(benchmark == "scattering") {
if(index == 0) {
formulation(index).integral(dof(lambda), tf(u(index)), gamma, gauss);
formulation(index).integral(dof(u(index)), tf(lambda), gamma, gauss);
formulation(index).integral(formulation.physicalSource(- fAnalytic), tf(lambda), gamma, gauss);
}
}
SubproblemDomains domains;
domains.setOmega(omega(index));
domains.setSigma({ sigma[0](index), sigma[1](index), sigma[2](index), sigma[3](index), sigma[4](index), sigma[5](index) });
domains.setPml({ pml[0](index), pml[1](index), pml[2](index), pml[3](index), pml[4](index), pml[5](index) });
domains.setPmlBnd({ std::make_pair(pmlBnd[0].first(index), pmlBnd[0].second(index)), std::make_pair(pmlBnd[1].first(index), pmlBnd[1].second(index)), std::make_pair(pmlBnd[2].first(index), pmlBnd[2].second(index)), std::make_pair(pmlBnd[3].first(index), pmlBnd[3].second(index)), std::make_pair(pmlBnd[4].first(index), pmlBnd[4].second(index)), std::make_pair(pmlBnd[5].first(index), pmlBnd[5].second(index)), std::make_pair(pmlBnd[6].first(index), pmlBnd[6].second(index)), std::make_pair(pmlBnd[7].first(index), pmlBnd[7].second(index)), std::make_pair(pmlBnd[8].first(index), pmlBnd[8].second(index)), std::make_pair(pmlBnd[9].first(index), pmlBnd[9].second(index)), std::make_pair(pmlBnd[10].first(index), pmlBnd[10].second(index)), std::make_pair(pmlBnd[11].first(index), pmlBnd[11].second(index)) });
domains.setPmlEdge({ pmlEdge[0](index), pmlEdge[1](index), pmlEdge[2](index), pmlEdge[3](index), pmlEdge[4](index), pmlEdge[5](index), pmlEdge[6](index), pmlEdge[7](index), pmlEdge[8](index), pmlEdge[9](index), pmlEdge[10](index), pmlEdge[11](index) });
domains.setEdge({ edge[0](index), edge[1](index), edge[2](index), edge[3](index), edge[4](index), edge[5](index), edge[6](index), edge[7](index), edge[8](index), edge[9](index), edge[10](index), edge[11](index) });
domains.setPmlEdgeBnd({ std::make_tuple(std::get<0>(pmlEdgeBnd[0])(index), std::get<1>(pmlEdgeBnd[0])(index), std::get<2>(pmlEdgeBnd[0])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[1])(index), std::get<1>(pmlEdgeBnd[1])(index), std::get<2>(pmlEdgeBnd[1])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[2])(index), std::get<1>(pmlEdgeBnd[2])(index), std::get<2>(pmlEdgeBnd[2])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[3])(index), std::get<1>(pmlEdgeBnd[3])(index), std::get<2>(pmlEdgeBnd[3])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[4])(index), std::get<1>(pmlEdgeBnd[4])(index), std::get<2>(pmlEdgeBnd[4])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[5])(index), std::get<1>(pmlEdgeBnd[5])(index), std::get<2>(pmlEdgeBnd[5])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[6])(index), std::get<1>(pmlEdgeBnd[6])(index), std::get<2>(pmlEdgeBnd[6])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[7])(index), std::get<1>(pmlEdgeBnd[7])(index), std::get<2>(pmlEdgeBnd[7])(index)) });
domains.setPmlCorner({ pmlCorner[0](index), pmlCorner[1](index), pmlCorner[2](index), pmlCorner[3](index), pmlCorner[4](index), pmlCorner[5](index), pmlCorner[6](index), pmlCorner[7](index) });
domains.setCorner({ corner[0](index), corner[1](index), corner[2](index), corner[3](index), corner[4](index), corner[5](index), corner[6](index), corner[7](index) });
domains.setCornerEdge({ std::make_tuple(std::get<0>(cornerEdge[0])(index), std::get<1>(cornerEdge[0])(index), std::get<2>(cornerEdge[0])(index)), std::make_tuple(std::get<0>(cornerEdge[1])(index), std::get<1>(cornerEdge[1])(index), std::get<2>(cornerEdge[1])(index)), std::make_tuple(std::get<0>(cornerEdge[2])(index), std::get<1>(cornerEdge[2])(index), std::get<2>(cornerEdge[2])(index)), std::make_tuple(std::get<0>(cornerEdge[3])(index), std::get<1>(cornerEdge[3])(index), std::get<2>(cornerEdge[3])(index)), std::make_tuple(std::get<0>(cornerEdge[4])(index), std::get<1>(cornerEdge[4])(index), std::get<2>(cornerEdge[4])(index)), std::make_tuple(std::get<0>(cornerEdge[5])(index), std::get<1>(cornerEdge[5])(index), std::get<2>(cornerEdge[5])(index)), std::make_tuple(std::get<0>(cornerEdge[6])(index), std::get<1>(cornerEdge[6])(index), std::get<2>(cornerEdge[6])(index)), std::make_tuple(std::get<0>(cornerEdge[7])(index), std::get<1>(cornerEdge[7])(index), std::get<2>(cornerEdge[7])(index)) });
SubproblemParameters parameters;
parameters.setGauss(gauss);
parameters.setKappa(kappa);
parameters.setNeumannOrder(neumannOrder);
parameters.setFieldOrder(fieldOrder);
parameters.setStab(stab);
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;
}
}
std::string bndExt = boundaryExt;
if(boundaryExt == "habc") {
bndExt += "_" + std::to_string(NExt) + "_" + std::to_string(thetaPade);
}
else if(boundaryExt == "pml") {
if(pmlMethodExt == "continuous") {
bndExt += "Continuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
}
else if(pmlMethodExt == "discontinuous") {
bndExt += "Discontinuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
}
}
gmshfem::msg::info << "Subdomain (" << i << ", " << j << ", " << k << ") has boundaries [" << (i == nDomX-1 ? bndExt : bnd) << ", " << (j == nDomY-1 ? bndExt : bnd) << ", " << (i == 0 ? bndExt : bnd) << ", " << (j == 0 ? bndExt : bnd) << ", " << (k == 0 ? bndExt : bnd) << ", " << (k == nDomZ-1 ? bndExt : bnd) << "]" << gmshfem::msg::endl;
subproblem[i][j][k] = new Subproblem(formulation(index), u(index).name(), domains, parameters, (i == nDomX-1 ? bndExt : bnd), (j == nDomY-1 ? bndExt : bnd), (i == 0 ? bndExt : bnd), (j == 0 ? bndExt : bnd), (k == 0 ? bndExt : bnd), (k == nDomZ-1 ? bndExt : bnd));
subproblem[i][j][k]->writeFormulation();
// coupling
const int faceNeighbor[6] = {
2, 3, 0, 1, 5, 4
};
for(unsigned int b = 0; b < 6; ++b) {
for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
const unsigned int jndex = topology[index][jj];
Boundary *boundary = subproblem[i][j][k]->getBoundary(b);
if(dynamic_cast< Sommerfeld * >(boundary)) {
if(!sigmaInterface[b](index, jndex).isEmpty()) {
formulation(index).integral(formulation.artificialSource( - (*gContinuous[faceNeighbor[b]])(jndex, index)), tf(u(index)), sigmaInterface[b](index, jndex), gauss);
}
}
else if(dynamic_cast< HABC * >(boundary)) {
if(!sigmaInterface[b](index, jndex).isEmpty()) {
formulation(index).integral(formulation.artificialSource( - (*gContinuous[faceNeighbor[b]])(jndex, index)), tf(u(index)), sigmaInterface[b](index, jndex), gauss);
}
}
else if(dynamic_cast< PmlContinuous * >(boundary)) {
if(!sigmaInterface[b](index, jndex).isEmpty()) {
formulation(index).integral(formulation.artificialSource( - (*gContinuous[faceNeighbor[b]])(jndex, index)), tf(u(index)), sigmaInterface[b](index, jndex), gauss);
}
}
}
}
const int edgeNeighbor[12][2] = {
// Down
{8, 2},
{9, 3},
{10, 0},
{11, 1},
// Rot
{7, 5},
{4, 6},
{5, 7},
{6, 4},
// Up
{0, 10},
{1, 11},
{2, 8},
{3, 9}
};
for(unsigned int e = 0; e < 12; ++e) {
for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
const unsigned int jndex = topology[index][jj];
Edge *edge = subproblem[i][j][k]->getEdge(e);
if(edge == nullptr) {
continue;
}
if(auto eg = dynamic_cast< HABC_HABC * >(edge)) {
if(!edgeInterface[e].first(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
for(unsigned int n = 0; n < N; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index) ), tf(*eg->firstBoundary()->getUHABC(n)), edgeInterface[e].first(index, jndex), gauss);
}
}
else {
for(unsigned int n = 0; n < N; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index) ), tf(*eg->firstBoundary()->getUHABC(n)), edgeInterface[e].first(index, jndex), gauss);
}
}
}
if(!edgeInterface[e].second(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
for(unsigned int n = 0; n < N; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index) ), tf(*eg->secondBoundary()->getUHABC(n)), edgeInterface[e].second(index, jndex), gauss);
}
}
else {
for(unsigned int n = 0; n < N; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index) ), tf(*eg->secondBoundary()->getUHABC(n)), edgeInterface[e].second(index, jndex), gauss);
}
}
}
}
else if(auto eg = dynamic_cast< PmlContinuous_PmlContinuous * >(edge)) {
if(!pmlBndInterface[e].first(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index) ), tf(*eg->firstBoundary()->getUPml()), pmlBndInterface[e].first(index, jndex), gauss);
}
else {
formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index) ), tf(*eg->firstBoundary()->getUPml()), pmlBndInterface[e].first(index, jndex), gauss);
}
}
if(!pmlBndInterface[e].second(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index) ), tf(*eg->secondBoundary()->getUPml()), pmlBndInterface[e].second(index, jndex), gauss);
}
else {
formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index) ), tf(*eg->secondBoundary()->getUPml()), pmlBndInterface[e].second(index, jndex), gauss);
}
}
}
}
}
const int cornerNeighbor[8][3] = {
// Down
{3, 1, 4},
{0, 2, 5},
{1, 3, 6},
{2, 0, 7},
// Up
{7, 5, 0},
{4, 6, 1},
{5, 7, 2},
{6, 4, 3}
};
for(unsigned int c = 0; c < 8; ++c) {
for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
const unsigned int jndex = topology[index][jj];
Corner *corner = subproblem[i][j][k]->getCorner(c);
if(corner == nullptr) {
continue;
}
if(auto cr = dynamic_cast< HABC_HABC_HABC * >(corner)) {
if(!std::get<0>(cornerInterface[c])(index, jndex).isEmpty()) {
for(unsigned int n = 0; n < N; ++n) {
for(unsigned int m = 0; m < N; ++m) {
formulation(index).integral(formulation.artificialSource( - (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index) ), tf(*cr->firstEdge()->getUHABC(n, m)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
}
}
}
if(!std::get<1>(cornerInterface[c])(index, jndex).isEmpty()) {
for(unsigned int n = 0; n < N; ++n) {
for(unsigned int m = 0; m < N; ++m) {
formulation(index).integral(formulation.artificialSource( - (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index) ), tf(*cr->secondEdge()->getUHABC(n, m)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
}
}
}
if(!std::get<2>(cornerInterface[c])(index, jndex).isEmpty()) {
for(unsigned int n = 0; n < N; ++n) {
for(unsigned int m = 0; m < N; ++m) {
formulation(index).integral(formulation.artificialSource( - (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index) ), tf(*cr->thirdEdge()->getUHABC(n, m)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
}
}
}
}
else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
formulation(index).integral(formulation.artificialSource( - (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index) ), tf(*cr->firstEdge()->getUPml()), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
}
if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
formulation(index).integral(formulation.artificialSource( - (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index) ), tf(*cr->secondEdge()->getUPml()), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
}
if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
formulation(index).integral(formulation.artificialSource( - (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index) ), tf(*cr->thirdEdge()->getUPml()), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
}
}
}
}
// interface
for(unsigned int b = 0; b < 6; ++b) {
for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
const unsigned int jndex = topology[index][jj];
Boundary *boundary = subproblem[i][j][k]->getBoundary(b);
if(auto bnd = dynamic_cast< Sommerfeld * >(boundary)) {
if(!sigmaInterface[b](index, jndex).isEmpty()) {
formulation(index, jndex).integral(dof((*gContinuous[b])(index, jndex)), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gContinuous[faceNeighbor[b]])(jndex, index) ), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
formulation(index, jndex).integral(2. * *bnd->getV(), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
}
}
else if(auto bnd = dynamic_cast< HABC * >(boundary)) {
if(!sigmaInterface[b](index, jndex).isEmpty()) {
formulation(index, jndex).integral(dof((*gContinuous[b])(index, jndex)), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gContinuous[faceNeighbor[b]])(jndex, index) ), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
formulation(index, jndex).integral(2. * *bnd->getV(), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
}
}
else if(auto bnd = dynamic_cast< PmlContinuous * >(boundary)) {
if(!sigmaInterface[b](index, jndex).isEmpty()) {
formulation(index, jndex).integral(dof((*gContinuous[b])(index, jndex)), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gContinuous[faceNeighbor[b]])(jndex, index) ), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
formulation(index, jndex).integral(2. * *bnd->getV(), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
}
}
}
}
for(unsigned int e = 0; e < 12; ++e) {
for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
const unsigned int jndex = topology[index][jj];
Edge *edge = subproblem[i][j][k]->getEdge(e);
if(edge == nullptr) {
continue;
}
if(auto eg = dynamic_cast< HABC_HABC * >(edge)) {
if(!edgeInterface[e].first(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
for(unsigned int n = 0; n < N; ++n) {
formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].first)(index, jndex)), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index) ), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(n, 0), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
}
}
else {
for(unsigned int n = 0; n < N; ++n) {
formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].first)(index, jndex)), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index) ), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(n, 0), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
}
}
}
if(!edgeInterface[e].second(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
for(unsigned int n = 0; n < N; ++n) {
formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].second)(index, jndex)), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index) ), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(n, 1), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
}
}
else {
for(unsigned int n = 0; n < N; ++n) {
formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].second)(index, jndex)), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index) ), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(n, 1), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
}
}
}
}
else if(auto eg = dynamic_cast< PmlContinuous_PmlContinuous * >(edge)) {
if(!pmlBndInterface[e].first(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].first)(index, jndex)), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index) ), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(0), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
}
else {
formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].first)(index, jndex)), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index) ), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(0), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
}
}
if(!pmlBndInterface[e].second(index, jndex).isEmpty()) {
if(e >= 4 && e < 8) {
formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].second)(index, jndex)), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index) ), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(1), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
}
else {
formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].second)(index, jndex)), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index) ), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
formulation(index, jndex).integral(2. * *eg->getV(1), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
}
}
}
}
}
for(unsigned int c = 0; c < 8; ++c) {
for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
const unsigned int jndex = topology[index][jj];
Corner *corner = subproblem[i][j][k]->getCorner(c);
if(corner == nullptr) {
continue;
}
if(auto cr = dynamic_cast< HABC_HABC_HABC * >(corner)) {
if(!std::get<0>(cornerInterface[c])(index, jndex).isEmpty()) {
for(unsigned int n = 0; n < N; ++n) {
for(unsigned int m = 0; m < N; ++m) {
formulation(index, jndex).integral(dof((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), tf((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index) ), tf((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(n, m, 0), tf((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
}
}
}
if(!std::get<1>(cornerInterface[c])(index, jndex).isEmpty()) {
for(unsigned int n = 0; n < N; ++n) {
for(unsigned int m = 0; m < N; ++m) {
formulation(index, jndex).integral(dof((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), tf((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index) ), tf((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(n, m, 1), tf((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
}
}
}
if(!std::get<2>(cornerInterface[c])(index, jndex).isEmpty()) {
for(unsigned int n = 0; n < N; ++n) {
for(unsigned int m = 0; m < N; ++m) {
formulation(index, jndex).integral(dof((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), tf((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index) ), tf((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(n, m, 2), tf((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
}
}
}
}
else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
formulation(index, jndex).integral(dof((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index) ), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(0), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
}
if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
formulation(index, jndex).integral(dof((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index) ), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(1), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
}
if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
formulation(index, jndex).integral(dof((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index) ), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(2), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
}
}
}
}
// for(unsigned int e = 0; e < 12; ++e) {
// for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
// const unsigned int jndex = topology[index][jj];
//
// Edge *edge = subproblem[i][j][k]->getEdge(e);
// if(edge == nullptr) {
// continue;
// }
// if(auto eg = dynamic_cast< HABC_HABC * >(edge)) {
// if(!edgeInterface[e].first(index, jndex).isEmpty()) {
// if(e >= 4 && e < 8) {
// for(unsigned int n = 0; n < N; ++n) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].first)(index, jndex), (*gEdgeHABC[e][n].first)(index, jndex).domain(), (*gEdgeHABC[e][n].first)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index).name() + "_pair");
// }
// }
// else {
// for(unsigned int n = 0; n < N; ++n) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].first)(index, jndex), (*gEdgeHABC[e][n].first)(index, jndex).domain(), (*gEdgeHABC[e][n].first)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index).name() + "_pair");
// }
// }
// }
// if(!edgeInterface[e].second(index, jndex).isEmpty()) {
// if(e >= 4 && e < 8) {
// for(unsigned int n = 0; n < N; ++n) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].second)(index, jndex), (*gEdgeHABC[e][n].second)(index, jndex).domain(), (*gEdgeHABC[e][n].second)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index).name() + "_pair");
// }
// }
// else {
// for(unsigned int n = 0; n < N; ++n) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].second)(index, jndex), (*gEdgeHABC[e][n].second)(index, jndex).domain(), (*gEdgeHABC[e][n].second)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index).name() + "_pair");
// }
// }
// }
// }
// else if(auto eg = dynamic_cast< PmlContinuous_PmlContinuous * >(edge)) {
// if(!pmlBndInterface[e].first(index, jndex).isEmpty()) {
// if(e >= 4 && e < 8) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].first)(index, jndex), (*gEdgePmlContinuous[e].first)(index, jndex).domain(), (*gEdgePmlContinuous[e].first)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index).name() + "_pair");
// }
// else {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].first)(index, jndex), (*gEdgePmlContinuous[e].first)(index, jndex).domain(), (*gEdgePmlContinuous[e].first)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index).name() + "_pair");
// }
// }
// if(!pmlBndInterface[e].second(index, jndex).isEmpty()) {
// if(e >= 4 && e < 8) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].second)(index, jndex), (*gEdgePmlContinuous[e].second)(index, jndex).domain(), (*gEdgePmlContinuous[e].second)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index).name() + "_pair");
// }
// else {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].second)(index, jndex), (*gEdgePmlContinuous[e].second)(index, jndex).domain(), (*gEdgePmlContinuous[e].second)(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index).name() + "_pair");
// }
// }
// }
// }
// }
// for(unsigned int c = 0; c < 8; ++c) {
// for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
// const unsigned int jndex = topology[index][jj];
//
// Corner *corner = subproblem[i][j][k]->getCorner(c);
// if(corner == nullptr) {
// continue;
// }
// if(auto cr = dynamic_cast< HABC_HABC_HABC * >(corner)) {
// if(!std::get<0>(cornerInterface[c])(index, jndex).isEmpty()) {
// for(unsigned int n = 0; n < N; ++n) {
// for(unsigned int m = 0; m < N; ++m) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerHABC[c][n][m]))(index, jndex), (*std::get<0>(gCornerHABC[c][n][m]))(index, jndex).domain(), (*std::get<0>(gCornerHABC[c][n][m]))(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index), (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index).domain(), (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index).name() + "_pair");
// }
// }
// }
// if(!std::get<1>(cornerInterface[c])(index, jndex).isEmpty()) {
// for(unsigned int n = 0; n < N; ++n) {
// for(unsigned int m = 0; m < N; ++m) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerHABC[c][n][m]))(index, jndex), (*std::get<1>(gCornerHABC[c][n][m]))(index, jndex).domain(), (*std::get<1>(gCornerHABC[c][n][m]))(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index), (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index).domain(), (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index).name() + "_pair");
// }
// }
// }
// if(!std::get<2>(cornerInterface[c])(index, jndex).isEmpty()) {
// for(unsigned int n = 0; n < N; ++n) {
// for(unsigned int m = 0; m < N; ++m) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerHABC[c][n][m]))(index, jndex), (*std::get<2>(gCornerHABC[c][n][m]))(index, jndex).domain(), (*std::get<2>(gCornerHABC[c][n][m]))(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index), (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index).domain(), (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index).name() + "_pair");
// }
// }
// }
// }
// else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
// if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerPmlContinuous[c]))(index, jndex), (*std::get<0>(gCornerPmlContinuous[c]))(index, jndex).domain(), (*std::get<0>(gCornerPmlContinuous[c]))(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index), (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index).domain(), (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index).name() + "_pair");
// }
// if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerPmlContinuous[c]))(index, jndex), (*std::get<1>(gCornerPmlContinuous[c]))(index, jndex).domain(), (*std::get<1>(gCornerPmlContinuous[c]))(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index), (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index).domain(), (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index).name() + "_pair");
// }
// if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerPmlContinuous[c]))(index, jndex), (*std::get<2>(gCornerPmlContinuous[c]))(index, jndex).domain(), (*std::get<2>(gCornerPmlContinuous[c]))(index, jndex).name());
// gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index), (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index).domain(), (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index).name() + "_pair");
// }
// }
// }
// }
// exit(0);
}
}
}
formulation.pre();
formulation.solve("gmres", res, iterMax);
for(unsigned int i = 0; i < nDomX; ++i) {
for(unsigned int j = 0; j < nDomY; ++j) {
for(unsigned int k = 0; k < nDomZ; ++k) {
unsigned int index = i * nDomY * nDomZ + j * nDomZ + k;
if(gmshddm::mpi::isItMySubdomain(index)) {
if(saveU) {
gmshfem::post::save(u(index), omega(index), "u_" + std::to_string(index), "msh");
}
}
}
}
}
if(monoDomainError) {
gmshfem::problem::Formulation< std::complex< double > > formulationMono("HelmholtzMonoDomain");
std::vector< gmshfem::domain::Domain > sigmaMono(6);
std::vector< gmshfem::domain::Domain > pmlMono(6);
std::vector< gmshfem::domain::Domain > pmlEdgeMono(12);
std::vector< gmshfem::domain::Domain > pmlCornerMono(8);
std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlBndMono(12);
std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlEdgeBndMono(8);
std::vector< gmshfem::domain::Domain > edgeMono(12);
std::vector< gmshfem::domain::Domain > cornerMono(8);
std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > cornerEdgeMono(8);
for(unsigned int i = 0; i < nDomX; ++i) {
for(unsigned int j = 0; j < nDomY; ++j) {
for(unsigned int k = 0; k < nDomZ; ++k) {
unsigned int index = i * nDomY * nDomZ + j * nDomZ + k;
if(i == 0) {
sigmaMono[2] |= sigma[2](index);
pmlMono[2] |= pml[2](index);
if(j == 0) {
pmlEdgeMono[6] |= pmlEdge[6](index);
edgeMono[6] |= edge[6](index);
pmlBndMono[6].first |= pmlBnd[6].first(index);
pmlBndMono[6].second |= pmlBnd[6].second(index);
if(k == 0) {
pmlCornerMono[2] |= pmlCorner[2](index);
cornerMono[2] |= corner[2](index);
std::get<0>(pmlEdgeBndMono[2]) |= std::get<0>(pmlEdgeBnd[2])(index);
std::get<1>(pmlEdgeBndMono[2]) |= std::get<1>(pmlEdgeBnd[2])(index);
std::get<2>(pmlEdgeBndMono[2]) |= std::get<2>(pmlEdgeBnd[2])(index);
std::get<0>(cornerEdgeMono[2]) |= std::get<0>(cornerEdge[2])(index);
std::get<1>(cornerEdgeMono[2]) |= std::get<1>(cornerEdge[2])(index);
std::get<2>(cornerEdgeMono[2]) |= std::get<2>(cornerEdge[2])(index);
}
if(k == nDomZ-1) {
pmlCornerMono[6] |= pmlCorner[6](index);
cornerMono[6] |= corner[6](index);
std::get<0>(pmlEdgeBndMono[6]) |= std::get<0>(pmlEdgeBnd[6])(index);
std::get<1>(pmlEdgeBndMono[6]) |= std::get<1>(pmlEdgeBnd[6])(index);
std::get<2>(pmlEdgeBndMono[6]) |= std::get<2>(pmlEdgeBnd[6])(index);
std::get<0>(cornerEdgeMono[6]) |= std::get<0>(cornerEdge[6])(index);
std::get<1>(cornerEdgeMono[6]) |= std::get<1>(cornerEdge[6])(index);
std::get<2>(cornerEdgeMono[6]) |= std::get<2>(cornerEdge[6])(index);
}
}
if(k == 0) {
pmlEdgeMono[2] |= pmlEdge[2](index);
edgeMono[2] |= edge[2](index);
pmlBndMono[2].first |= pmlBnd[2].first(index);
pmlBndMono[2].second |= pmlBnd[2].second(index);
}
if(j == nDomY-1) {
pmlEdgeMono[5] |= pmlEdge[5](index);
edgeMono[5] |= edge[5](index);
pmlBndMono[5].first |= pmlBnd[5].first(index);
pmlBndMono[5].second |= pmlBnd[5].second(index);
if(k == 0) {
pmlCornerMono[1] |= pmlCorner[1](index);
cornerMono[1] |= corner[1](index);
std::get<0>(pmlEdgeBndMono[1]) |= std::get<0>(pmlEdgeBnd[1])(index);
std::get<1>(pmlEdgeBndMono[1]) |= std::get<1>(pmlEdgeBnd[1])(index);
std::get<2>(pmlEdgeBndMono[1]) |= std::get<2>(pmlEdgeBnd[1])(index);
std::get<0>(cornerEdgeMono[1]) |= std::get<0>(cornerEdge[1])(index);
std::get<1>(cornerEdgeMono[1]) |= std::get<1>(cornerEdge[1])(index);
std::get<2>(cornerEdgeMono[1]) |= std::get<2>(cornerEdge[1])(index);
}
if(k == nDomZ-1) {
pmlCornerMono[5] |= pmlCorner[5](index);
cornerMono[5] |= corner[5](index);
std::get<0>(pmlEdgeBndMono[5]) |= std::get<0>(pmlEdgeBnd[5])(index);
std::get<1>(pmlEdgeBndMono[5]) |= std::get<1>(pmlEdgeBnd[5])(index);
std::get<2>(pmlEdgeBndMono[5]) |= std::get<2>(pmlEdgeBnd[5])(index);
std::get<0>(cornerEdgeMono[5]) |= std::get<0>(cornerEdge[5])(index);
std::get<1>(cornerEdgeMono[5]) |= std::get<1>(cornerEdge[5])(index);
std::get<2>(cornerEdgeMono[5]) |= std::get<2>(cornerEdge[5])(index);
}
}
if(k == nDomZ-1) {
pmlEdgeMono[10] |= pmlEdge[10](index);
edgeMono[10] |= edge[10](index);
pmlBndMono[10].first |= pmlBnd[10].first(index);
pmlBndMono[10].second |= pmlBnd[10].second(index);
}
}
if(j == 0) {
sigmaMono[3] |= sigma[3](index);
pmlMono[3] |= pml[3](index);
if(k == 0) {
pmlEdgeMono[3] |= pmlEdge[3](index);
edgeMono[3] |= edge[3](index);
pmlBndMono[3].first |= pmlBnd[3].first(index);
pmlBndMono[3].second |= pmlBnd[3].second(index);
}
if(k == nDomZ-1) {
pmlEdgeMono[11] |= pmlEdge[11](index);
edgeMono[11] |= edge[11](index);
pmlBndMono[11].first |= pmlBnd[11].first(index);
pmlBndMono[11].second |= pmlBnd[11].second(index);
}
}
if(k == 0) {
sigmaMono[4] |= sigma[4](index);
pmlMono[4] |= pml[4](index);
}
if(i == nDomX-1) {
sigmaMono[0] |= sigma[0](index);
pmlMono[0] |= pml[0](index);
if(j == 0) {
pmlEdgeMono[7] |= pmlEdge[7](index);
edgeMono[7] |= edge[7](index);
pmlBndMono[7].first |= pmlBnd[7].first(index);
pmlBndMono[7].second |= pmlBnd[7].second(index);
if(k == 0) {
pmlCornerMono[3] |= pmlCorner[3](index);
cornerMono[3] |= corner[3](index);
std::get<0>(pmlEdgeBndMono[3]) |= std::get<0>(pmlEdgeBnd[3])(index);
std::get<1>(pmlEdgeBndMono[3]) |= std::get<1>(pmlEdgeBnd[3])(index);
std::get<2>(pmlEdgeBndMono[3]) |= std::get<2>(pmlEdgeBnd[3])(index);
std::get<0>(cornerEdgeMono[3]) |= std::get<0>(cornerEdge[3])(index);
std::get<1>(cornerEdgeMono[3]) |= std::get<1>(cornerEdge[3])(index);
std::get<2>(cornerEdgeMono[3]) |= std::get<2>(cornerEdge[3])(index);
}
if(k == nDomZ-1) {
pmlCornerMono[7] |= pmlCorner[7](index);
cornerMono[7] |= corner[7](index);
std::get<0>(pmlEdgeBndMono[7]) |= std::get<0>(pmlEdgeBnd[7])(index);
std::get<1>(pmlEdgeBndMono[7]) |= std::get<1>(pmlEdgeBnd[7])(index);
std::get<2>(pmlEdgeBndMono[7]) |= std::get<2>(pmlEdgeBnd[7])(index);
std::get<0>(cornerEdgeMono[7]) |= std::get<0>(cornerEdge[7])(index);
std::get<1>(cornerEdgeMono[7]) |= std::get<1>(cornerEdge[7])(index);
std::get<2>(cornerEdgeMono[7]) |= std::get<2>(cornerEdge[7])(index);
}
}
if(k == 0) {
pmlEdgeMono[0] |= pmlEdge[0](index);
edgeMono[0] |= edge[0](index);
pmlBndMono[0].first |= pmlBnd[0].first(index);
pmlBndMono[0].second |= pmlBnd[0].second(index);
}
if(j == nDomY-1) {
pmlEdgeMono[4] |= pmlEdge[4](index);
edgeMono[4] |= edge[4](index);
pmlBndMono[4].first |= pmlBnd[4].first(index);
pmlBndMono[4].second |= pmlBnd[4].second(index);
if(k == 0) {
pmlCornerMono[0] |= pmlCorner[0](index);
cornerMono[0] |= corner[0](index);
std::get<0>(pmlEdgeBndMono[0]) |= std::get<0>(pmlEdgeBnd[0])(index);
std::get<1>(pmlEdgeBndMono[0]) |= std::get<1>(pmlEdgeBnd[0])(index);
std::get<2>(pmlEdgeBndMono[0]) |= std::get<2>(pmlEdgeBnd[0])(index);
std::get<0>(cornerEdgeMono[0]) |= std::get<0>(cornerEdge[0])(index);
std::get<1>(cornerEdgeMono[0]) |= std::get<1>(cornerEdge[0])(index);
std::get<2>(cornerEdgeMono[0]) |= std::get<2>(cornerEdge[0])(index);
}
if(k == nDomZ-1) {
pmlCornerMono[4] |= pmlCorner[4](index);
cornerMono[4] |= corner[4](index);
std::get<0>(pmlEdgeBndMono[4]) |= std::get<0>(pmlEdgeBnd[4])(index);
std::get<1>(pmlEdgeBndMono[4]) |= std::get<1>(pmlEdgeBnd[4])(index);
std::get<2>(pmlEdgeBndMono[4]) |= std::get<2>(pmlEdgeBnd[4])(index);
std::get<0>(cornerEdgeMono[4]) |= std::get<0>(cornerEdge[4])(index);
std::get<1>(cornerEdgeMono[4]) |= std::get<1>(cornerEdge[4])(index);
std::get<2>(cornerEdgeMono[4]) |= std::get<2>(cornerEdge[4])(index);
}
}
if(k == nDomZ-1) {
pmlEdgeMono[8] |= pmlEdge[8](index);
edgeMono[8] |= edge[8](index);
pmlBndMono[8].first |= pmlBnd[8].first(index);
pmlBndMono[8].second |= pmlBnd[8].second(index);
}
}
if(j == nDomY-1) {
sigmaMono[1] |= sigma[1](index);
pmlMono[1] |= pml[1](index);
if(k == 0) {
pmlEdgeMono[1] |= pmlEdge[1](index);
edgeMono[1] |= edge[1](index);
pmlBndMono[1].first |= pmlBnd[1].first(index);
pmlBndMono[1].second |= pmlBnd[1].second(index);
}
if(k == nDomZ-1) {
pmlEdgeMono[9] |= pmlEdge[9](index);
edgeMono[9] |= edge[9](index);
pmlBndMono[9].first |= pmlBnd[9].first(index);
pmlBndMono[9].second |= pmlBnd[9].second(index);
}
}
if(k == nDomZ-1) {
sigmaMono[5] |= sigma[5](index);
pmlMono[5] |= pml[5](index);
}
}
}
}
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > uMono("uMono", omega.getUnion(), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambdaMono("lamdbaMono", gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
formulationMono.integral(-grad(dof(uMono)), grad(tf(uMono)), omega.getUnion(), gauss);
formulationMono.integral(kappa * kappa * dof(uMono), tf(uMono), omega.getUnion(), gauss);
formulationMono.integral(dof(lambda), tf(uMono), gamma, gauss);
formulationMono.integral(dof(uMono), tf(lambda), gamma, gauss);
if(benchmark == "scattering") {
formulationMono.integral(- fAnalytic, tf(lambda), gamma, gauss);
}
SubproblemDomains domains;
domains.setOmega(omega.getUnion());
domains.setSigma(sigmaMono);
domains.setPml(pmlMono);
domains.setPmlBnd(pmlBndMono);
domains.setPmlEdge(pmlEdgeMono);
domains.setEdge(edgeMono);
domains.setPmlEdgeBnd(pmlEdgeBndMono);
domains.setPmlCorner(pmlCornerMono);
domains.setCorner(cornerMono);
domains.setCornerEdge(cornerEdgeMono);
SubproblemParameters parameters;
parameters.setGauss(gauss);
parameters.setKappa(kappa);
parameters.setNeumannOrder(neumannOrder);
parameters.setFieldOrder(fieldOrder);
parameters.setStab(stab);
std::string bndExt = boundaryExt;
if(boundaryExt == "habc") {
bndExt += "_" + std::to_string(NExt) + "_" + std::to_string(thetaPade);
}
else if(boundaryExt == "pml") {
if(pmlMethodExt == "continuous") {
bndExt += "Continuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
}
else if(pmlMethodExt == "discontinuous") {
bndExt += "Discontinuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
}
}
gmshfem::msg::info << "Monodomain has boundaries [" << bndExt << ", " << bndExt << ", " << bndExt << ", " << bndExt << ", " << bndExt << ", " << bndExt << "]" << gmshfem::msg::endl;
Subproblem subproblem(formulationMono, uMono.name(), domains, parameters, bndExt, bndExt, bndExt, bndExt, bndExt, bndExt);
subproblem.writeFormulation();
formulationMono.pre();
formulationMono.assemble();
formulationMono.solve();
double num = 0.;
double den = 0.;
for(unsigned int i = 0; i < nDomX * nDomY * nDomZ; ++i) {
num += std::real(gmshfem::post::integrate(pow(abs(u(i) - uMono), 2), omega(i), gauss));
den += std::real(gmshfem::post::integrate(pow(abs(uMono), 2), omega(i), gauss));
}
gmshfem::msg::info << "Error mono L2 = " << std::sqrt(num/den) << gmshfem::msg::endl;
for(unsigned int i = 0; i < nDomX * nDomY * nDomZ; ++i) {
if(gmshddm::mpi::isItMySubdomain(i)) {
if(saveEMono) {
gmshfem::post::save(u(i) - uMono, omega(i), "eMono_" + std::to_string(i));
}
if(saveUMono) {
gmshfem::post::save(uMono, omega(i), "uMono_" + std::to_string(i));
}
}
}
if(fileName != "none") {
gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
file << iterMax << formulation.relativeResidual().back() << std::sqrt(num/den) << gmshfem::csv::endl;
file.close();
}
}
if(wavenumberPlot) {
gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
file << k << formulation.relativeResidual().size()-1 << gmshfem::csv::endl;
file.close();
}
if(meshPlot) {
gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
file << (2*pi/k)/lc << lc << 1./lc << formulation.relativeResidual().size()-1 << gmshfem::csv::endl;
file.close();
}
}
}