From 6f659a831a69e3fc4c02f83a60d055772bb0e768 Mon Sep 17 00:00:00 2001
From: Anthony Royer <anthony.royer@uliege.be>
Date: Fri, 13 May 2022 14:29:40 +0200
Subject: [PATCH] HABC and corner with different number of auxillary fields
---
.../helmholtz/crossPoints/Subproblem2D.cpp | 43 +++++++------
examples/helmholtz/crossPoints/Subproblem2D.h | 2 +-
examples/helmholtz/crossPoints/ddm2D.cpp | 63 ++++++++++++++-----
examples/helmholtz/crossPoints/monoDomain.cpp | 2 +-
4 files changed, 74 insertions(+), 36 deletions(-)
diff --git a/examples/helmholtz/crossPoints/Subproblem2D.cpp b/examples/helmholtz/crossPoints/Subproblem2D.cpp
index 0cd84471..3d847794 100755
--- a/examples/helmholtz/crossPoints/Subproblem2D.cpp
+++ b/examples/helmholtz/crossPoints/Subproblem2D.cpp
@@ -95,7 +95,7 @@ namespace D2 {
}
}
- Subproblem::Subproblem(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters, const std::string &E, const std::string &N, const std::string &W, const std::string &S) : _domains(domains), _parameters(parameters), _formulation(formulation), _fieldName(fieldName), _boundary(4, nullptr), _corner(4, nullptr)
+ Subproblem::Subproblem(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters, const std::string &E, const std::string &N, const std::string &W, const std::string &S, const std::vector< unsigned int > activatedCrossPoints) : _domains(domains), _parameters(parameters), _formulation(formulation), _fieldName(fieldName), _boundary(4, nullptr), _corner(4, nullptr)
{
std::string type[4] = {E, N, W, S};
std::string method[4];
@@ -134,6 +134,11 @@ namespace D2 {
}
for(unsigned int c = 0; c < 4; ++c) {
+ auto it = std::find(activatedCrossPoints.begin(), activatedCrossPoints.end(), c);
+ if(it == activatedCrossPoints.end()) {
+ continue;
+ }
+
if(method[c] == "sommerfeld" && method[(c+1)%4] == "sommerfeld") {
}
else if(method[c] == "habc" && method[(c+1)%4] == "habc") {
@@ -303,7 +308,6 @@ namespace D2 {
const double pi = 3.141592653589793238462643383279;
const std::complex< double > expPTheta(std::cos(_theta), std::sin(_theta));
- const std::complex< double > expMTheta(std::cos(-_theta), std::sin(-_theta));
const std::complex< double > expPTheta2(std::cos(_theta/2.), std::sin(_theta/2.));
const double M = 2 * _N + 1;
double cPade[_N];
@@ -618,11 +622,7 @@ namespace D2 {
const unsigned int N1 = static_cast< HABC * >(_bnd[1])->getN();
const double theta0 = static_cast< HABC * >(_bnd[0])->getTheta();
- // const double theta1 = static_cast< HABC * >(_bnd[1])->getTheta();
-
- if(N0 != N1) {
- throw gmshfem::common::Exception("Corner with HABC having different number of auxilary fields");
- }
+ const double theta1 = static_cast< HABC * >(_bnd[1])->getTheta();
_uHABC_C.resize(N0);
for(unsigned int m = 0; m < N0; ++m) {
@@ -642,14 +642,19 @@ namespace D2 {
}
const double pi = 3.141592653589793238462643383279;
- const std::complex< double > expPTheta(std::cos(theta0), std::sin(theta0));
- const std::complex< double > expMTheta(std::cos(-theta0), std::sin(-theta0));
- const std::complex< double > expPTheta2(std::cos(theta0/2.), std::sin(theta0/2.));
+
+ const std::complex< double > expPTheta_0(std::cos(theta0), std::sin(theta0));
+ const std::complex< double > expPTheta2_0(std::cos(theta0/2.), std::sin(theta0/2.));
+
+ const std::complex< double > expPTheta_1(std::cos(theta1), std::sin(theta1));
+ const std::complex< double > expPTheta2_1(std::cos(theta1/2.), std::sin(theta1/2.));
+
const double M0 = 2 * N0 + 1;
double cPade0[N0];
for(unsigned int m = 0; m < N0; ++m) {
cPade0[m] = std::tan((m+1) * pi / M0) * std::tan((m+1) * pi / M0);
}
+
const double M1 = 2 * N1 + 1;
double cPade1[N1];
for(unsigned int n = 0; n < N1; ++n) {
@@ -667,11 +672,11 @@ namespace D2 {
formulation.integral(-dof(*_vC[0][m]), tf(*uHABC[0][m]), corner, gauss);
formulation.integral( dof(*_vC[0][m]), tf(*_vC[0][m]), corner, gauss);
- formulation.integral(-im * kappa * expPTheta2 * dof(*uHABC[0][m]), tf(*_vC[0][m]), corner, gauss);
+ formulation.integral(-im * kappa * expPTheta2_1 * dof(*uHABC[0][m]), tf(*_vC[0][m]), corner, gauss);
for(unsigned int n = 0; n < N1; ++n) {
- formulation.integral(-im * kappa * expPTheta2 * 2./M1 * cPade1[n] * dof(*uHABC[0][m]), tf(*_vC[0][m]), corner, gauss);
- formulation.integral(-im * kappa * expPTheta2 * 2./M1 * cPade1[n] * dof(*_uHABC_C[m][n]), tf(*_vC[0][m]), corner, gauss);
+ formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*uHABC[0][m]), tf(*_vC[0][m]), corner, gauss);
+ formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*_uHABC_C[m][n]), tf(*_vC[0][m]), corner, gauss);
}
}
@@ -680,20 +685,20 @@ namespace D2 {
formulation.integral(-dof(*_vC[1][n]), tf(*uHABC[1][n]), corner, gauss);
formulation.integral( dof(*_vC[1][n]), tf(*_vC[1][n]), corner, gauss);
- formulation.integral(-im * kappa * expPTheta2 * dof(*uHABC[1][n]), tf(*_vC[1][n]), corner, gauss);
+ formulation.integral(-im * kappa * expPTheta2_0 * dof(*uHABC[1][n]), tf(*_vC[1][n]), corner, gauss);
for(unsigned int m = 0; m < N0; ++m) {
- formulation.integral(-im * kappa * expPTheta2 * 2./M0 * cPade0[m] * dof(*uHABC[1][n]), tf(*_vC[1][n]), corner, gauss);
- formulation.integral(-im * kappa * expPTheta2 * 2./M0 * cPade0[m] * dof(*_uHABC_C[m][n]), tf(*_vC[1][n]), corner, gauss);
+ formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*uHABC[1][n]), tf(*_vC[1][n]), corner, gauss);
+ formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*_uHABC_C[m][n]), tf(*_vC[1][n]), corner, gauss);
}
}
// auxiliary
for(unsigned int m = 0; m < N0; ++m) {
for(unsigned int n = 0; n < N1; ++n) {
- formulation.integral( (cPade0[m] + cPade1[n] + expMTheta) * dof(*_uHABC_C[m][n]), tf(*_uHABC_C[m][n]), corner, gauss);
- formulation.integral( (cPade1[n] + 1.) * dof(*uHABC[0][m]), tf(*_uHABC_C[m][n]), corner, gauss);
- formulation.integral( (cPade0[m] + 1.) * dof(*uHABC[1][n]), tf(*_uHABC_C[m][n]), corner, gauss);
+ formulation.integral( (expPTheta_0 * cPade0[m] + expPTheta_1 * cPade1[n] + 1.) * dof(*_uHABC_C[m][n]), tf(*_uHABC_C[m][n]), corner, gauss);
+ formulation.integral( expPTheta_1 * (cPade1[n] + 1.) * dof(*uHABC[0][m]), tf(*_uHABC_C[m][n]), corner, gauss);
+ formulation.integral( expPTheta_0 * (cPade0[m] + 1.) * dof(*uHABC[1][n]), tf(*_uHABC_C[m][n]), corner, gauss);
}
}
}
diff --git a/examples/helmholtz/crossPoints/Subproblem2D.h b/examples/helmholtz/crossPoints/Subproblem2D.h
index 57bc7692..b153d52f 100755
--- a/examples/helmholtz/crossPoints/Subproblem2D.h
+++ b/examples/helmholtz/crossPoints/Subproblem2D.h
@@ -51,7 +51,7 @@ namespace D2 {
void _parseParameters(std::string &method, std::vector< std::string > ¶meters, const std::string &str) const;
public:
- Subproblem(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters, const std::string &E, const std::string &N, const std::string &W, const std::string &S);
+ Subproblem(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters ¶meters, const std::string &E, const std::string &N, const std::string &W, const std::string &S, const std::vector< unsigned int > activatedCrossPoints);
~Subproblem();
void writeFormulation();
diff --git a/examples/helmholtz/crossPoints/ddm2D.cpp b/examples/helmholtz/crossPoints/ddm2D.cpp
index a14e1364..1e984dea 100755
--- a/examples/helmholtz/crossPoints/ddm2D.cpp
+++ b/examples/helmholtz/crossPoints/ddm2D.cpp
@@ -108,6 +108,8 @@ namespace D2 {
double thetaPade = 0.;
gmshDdm->userDefinedParameter(thetaPade, "thetaPade");
thetaPade *= pi;
+ bool switchOffInteriorCrossPoints = false;
+ gmshDdm->userDefinedParameter(switchOffInteriorCrossPoints, "switchOffInteriorCrossPoints");
// ************************
// P O S T
@@ -190,7 +192,8 @@ namespace D2 {
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::msg::info << " - interior cross-points: " << (switchOffInteriorCrossPoints ? "switch off" : "switch on") << gmshfem::msg::endl;
+ gmshfem::msg::info << " - stab: " << stab << " * lc (= " << stab * lc << ")" << gmshfem::msg::endl;
// source
@@ -378,8 +381,8 @@ namespace D2 {
gCornerPmlDiscontinuous.push_back(std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[f] + dir[(f+1)%4] + "_first", pmlBndInterface[f].first, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[f] + dir[(f+1)%4] + "_second", pmlBndInterface[f].second, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
- gCornerHABC.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) {
+ gCornerHABC.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 > * > >(std::max(N, NExt)));
+ for(unsigned int n = 0; n < std::max(N, NExt); ++n) {
gCornerHABC[f][n] = std::make_pair(
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[f] + dir[(f+1)%4] + "_" + std::to_string(n) + "_first", cornerInterface[f].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[f] + dir[(f+1)%4] + "_" + std::to_string(n) + "_second", cornerInterface[f].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
@@ -399,7 +402,7 @@ namespace D2 {
formulation.addInterfaceField(*gCornerPmlDiscontinuous[f].second);
// HABC
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < std::max(N, NExt); ++n) {
formulation.addInterfaceField(*gCornerHABC[f][n].first);
formulation.addInterfaceField(*gCornerHABC[f][n].second);
}
@@ -474,9 +477,39 @@ namespace D2 {
}
}
gmshfem::msg::info << "Subdomain (" << i << ", " << j << ") has boundaries [" << (i == nDomX-1 ? bndExt : bnd) << ", " << (j == nDomY-1 ? bndExt : bnd) << ", " << (i == 0 ? bndExt : bnd) << ", " << (j == 0 ? bndExt : bnd) << "]" << gmshfem::msg::endl;
- subproblem[i][j] = 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));
+ std::vector< unsigned int > activatedCrossPoints = {0, 1, 2, 3};
+ if(switchOffInteriorCrossPoints) {
+ activatedCrossPoints.clear();
+ if(i == nDomX-1 && j == nDomY-1) {
+ activatedCrossPoints.push_back(0);
+ }
+ if(i == 0 && j == nDomY-1) {
+ activatedCrossPoints.push_back(1);
+ }
+ if(i == 0 && j == 0) {
+ activatedCrossPoints.push_back(2);
+ }
+ if(i == nDomX-1 && j == 0) {
+ activatedCrossPoints.push_back(3);
+ }
+ }
+ subproblem[i][j] = 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), activatedCrossPoints);
subproblem[i][j]->writeFormulation();
+ unsigned int Ns[4] = {N, N, N, N};
+ if(i == nDomX-1) {
+ Ns[0] = NExt;
+ }
+ if(j == nDomY-1) {
+ Ns[1] = NExt;
+ }
+ if(i == 0) {
+ Ns[2] = NExt;
+ }
+ if(j == 0) {
+ Ns[3] = NExt;
+ }
+
// coupling
for(unsigned int b = 0; b < 4; ++b) {
for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
@@ -521,12 +554,12 @@ namespace D2 {
}
if(auto cr = dynamic_cast< HABC_HABC * >(corner)) {
if(!cornerInterface[c].first(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[c]; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gCornerHABC[(c-1)%4][n].second)(jndex, index) ), tf(*cr->firstBoundary()->getUHABC(n)), cornerInterface[c].first(index, jndex), gauss);
}
}
if(!cornerInterface[c].second(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[(c+1)%4]; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gCornerHABC[(c+1)%4][n].first)(jndex, index) ), tf(*cr->secondBoundary()->getUHABC(n)), cornerInterface[c].second(index, jndex), gauss);
}
}
@@ -582,14 +615,14 @@ namespace D2 {
}
else if(auto cr = dynamic_cast< HABC_Sommerfeld * >(corner)) {
if(!cornerInterface[c].first(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[c]; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gCornerHABC[(c-1)%4][n].second)(jndex, index) ), tf(*cr->firstBoundary()->getUHABC(n)), cornerInterface[c].first(index, jndex), gauss);
}
}
}
else if(auto cr = dynamic_cast< Sommerfeld_HABC * >(corner)) {
if(!cornerInterface[c].second(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[(c+1)%4]; ++n) {
formulation(index).integral(formulation.artificialSource( - (*gCornerHABC[(c+1)%4][n].first)(jndex, index) ), tf(*cr->secondBoundary()->getUHABC(n)), cornerInterface[c].second(index, jndex), gauss);
}
}
@@ -651,14 +684,14 @@ namespace D2 {
}
if(auto cr = dynamic_cast< HABC_HABC * >(corner)) {
if(!cornerInterface[c].first(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[c]; ++n) {
formulation(index, jndex).integral(dof((*gCornerHABC[c][n].first)(index, jndex)), tf((*gCornerHABC[c][n].first)(index, jndex)), cornerInterface[c].first(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gCornerHABC[(c-1)%4][n].second)(jndex, index) ), tf((*gCornerHABC[c][n].first)(index, jndex)), cornerInterface[c].first(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(n, 0), tf((*gCornerHABC[c][n].first)(index, jndex)), cornerInterface[c].first(index, jndex), gauss);
}
}
if(!cornerInterface[c].second(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[(c+1)%4]; ++n) {
formulation(index, jndex).integral(dof((*gCornerHABC[c][n].second)(index, jndex)), tf((*gCornerHABC[c][n].second)(index, jndex)), cornerInterface[c].second(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gCornerHABC[(c+1)%4][n].first)(jndex, index) ), tf((*gCornerHABC[c][n].second)(index, jndex)), cornerInterface[c].second(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(n, 1), tf((*gCornerHABC[c][n].second)(index, jndex)), cornerInterface[c].second(index, jndex), gauss);
@@ -731,7 +764,7 @@ namespace D2 {
}
else if(auto cr = dynamic_cast< HABC_Sommerfeld * >(corner)) {
if(!cornerInterface[c].first(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[c]; ++n) {
formulation(index, jndex).integral(dof((*gCornerHABC[c][n].first)(index, jndex)), tf((*gCornerHABC[c][n].first)(index, jndex)), cornerInterface[c].first(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gCornerHABC[(c-1)%4][n].second)(jndex, index) ), tf((*gCornerHABC[c][n].first)(index, jndex)), cornerInterface[c].first(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(n), tf((*gCornerHABC[c][n].first)(index, jndex)), cornerInterface[c].first(index, jndex), gauss);
@@ -740,7 +773,7 @@ namespace D2 {
}
else if(auto cr = dynamic_cast< Sommerfeld_HABC * >(corner)) {
if(!cornerInterface[c].second(index, jndex).isEmpty()) {
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < Ns[(c+1)%4]; ++n) {
formulation(index, jndex).integral(dof((*gCornerHABC[c][n].second)(index, jndex)), tf((*gCornerHABC[c][n].second)(index, jndex)), cornerInterface[c].second(index, jndex), gauss);
formulation(index, jndex).integral(formulation.artificialSource( (*gCornerHABC[(c+1)%4][n].first)(jndex, index) ), tf((*gCornerHABC[c][n].second)(index, jndex)), cornerInterface[c].second(index, jndex), gauss);
formulation(index, jndex).integral(2. * *cr->getV(n), tf((*gCornerHABC[c][n].second)(index, jndex)), cornerInterface[c].second(index, jndex), gauss);
@@ -851,7 +884,7 @@ namespace D2 {
}
}
gmshfem::msg::info << "Monodomain has boundaries [" << bndExt << ", " << bndExt << ", " << bndExt << ", " << bndExt << "]" << gmshfem::msg::endl;
- Subproblem subproblem(formulationMono, uMono.name(), domains, parameters, bndExt, bndExt, bndExt, bndExt);
+ Subproblem subproblem(formulationMono, uMono.name(), domains, parameters, bndExt, bndExt, bndExt, bndExt, {0,1,2,3});
subproblem.writeFormulation();
formulationMono.pre();
@@ -933,7 +966,7 @@ namespace D2 {
delete gCornerPmlContinuous[f].second;
delete gCornerPmlDiscontinuous[f].first;
delete gCornerPmlDiscontinuous[f].second;
- for(unsigned int n = 0; n < N; ++n) {
+ for(unsigned int n = 0; n < std::max(N, NExt); ++n) {
delete gCornerHABC[f][n].first;
delete gCornerHABC[f][n].second;
}
diff --git a/examples/helmholtz/crossPoints/monoDomain.cpp b/examples/helmholtz/crossPoints/monoDomain.cpp
index d5f15947..dc891331 100755
--- a/examples/helmholtz/crossPoints/monoDomain.cpp
+++ b/examples/helmholtz/crossPoints/monoDomain.cpp
@@ -187,7 +187,7 @@ namespace D2 {
bnd += "WithoutMultiplier_" + std::to_string(pmlSize) + "_" + pmlType;
}
}
- Subproblem subproblem(formulation, u.name(), domains, parameters, bnd, bnd, bnd, bnd);
+ Subproblem subproblem(formulation, u.name(), domains, parameters, bnd, bnd, bnd, bnd, {0,1,2,3});
subproblem.writeFormulation();
if(spy) {
--
GitLab