From dfc6f49c413e5c4ca79ab4ee1af36075bd127a58 Mon Sep 17 00:00:00 2001
From: Anthony Royer <anthony.royer@uliege.be>
Date: Thu, 19 May 2022 09:10:06 +0200
Subject: [PATCH] Add non-right partitions
---
examples/helmholtz/crossPoints/ddm2D.cpp | 45 ++++-----
examples/helmholtz/crossPoints/mesh.cpp | 121 +++++++++++++++++++++--
examples/helmholtz/crossPoints/mesh.h | 3 +-
3 files changed, 133 insertions(+), 36 deletions(-)
diff --git a/examples/helmholtz/crossPoints/ddm2D.cpp b/examples/helmholtz/crossPoints/ddm2D.cpp
index ebaf8f8c..bf656f29 100755
--- a/examples/helmholtz/crossPoints/ddm2D.cpp
+++ b/examples/helmholtz/crossPoints/ddm2D.cpp
@@ -82,6 +82,8 @@ namespace D2 {
double sizeX = 2., sizeY = 2.;
gmshDdm->userDefinedParameter(sizeX, "sizeX");
gmshDdm->userDefinedParameter(sizeY, "sizeY");
+ double twist = 0.;
+ gmshDdm->userDefinedParameter(twist, "twist");
unsigned int N = 6;
gmshDdm->userDefinedParameter(N, "N");
@@ -140,7 +142,12 @@ namespace D2 {
gmshDdm->userDefinedParameter(saveIteration, "saveIteration");
if(benchmark == "scattering") {
- checkerboard(nDomX, nDomY, sizeX, sizeY, R, lc, (boundary == "pml"), pmlSize, (boundaryExt == "pml"), pmlSizeExt, meshOrder, 0, 0);
+ if(twist == 0.) {
+ checkerboard(nDomX, nDomY, sizeX, sizeY, R, lc, (boundary == "pml"), pmlSize, (boundaryExt == "pml"), pmlSizeExt, meshOrder, 0, 0);
+ }
+ else {
+ checkerboardTwist(nDomX, nDomY, sizeX, sizeY, R, lc, (boundary == "pml"), pmlSize, (boundaryExt == "pml"), pmlSizeExt, meshOrder, twist, 0, 0);
+ }
}
else if (benchmark == "marmousi") {
checkerboard(nDomX, nDomY, 9192./nDomX, 2904./nDomY, 0., lc, (boundary == "pml"), pmlSize, (boundaryExt == "pml"), pmlSizeExt, meshOrder, -1, -1, true);
@@ -167,6 +174,7 @@ namespace D2 {
gmshfem::msg::info << " - meshOrder: " << meshOrder << gmshfem::msg::endl;
gmshfem::msg::info << " * modeling:" << gmshfem::msg::endl;
gmshfem::msg::info << " - grid: (" << nDomX << ", " << nDomY << ")" << gmshfem::msg::endl;
+ gmshfem::msg::info << " - twist: " << twist*100. << "\%" << 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 << ")" << gmshfem::msg::endl;
@@ -364,8 +372,6 @@ namespace D2 {
else if (benchmark == " marmousi") {
u(37).domain(u(37).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 f = 0; f < 4; ++f) {
@@ -409,7 +415,14 @@ namespace D2 {
formulation.addInterfaceField(*gCornerHABC[f][n].second);
}
}
-
+
+ if(benchmark == "scattering") {
+ u(0).addConstraint(gamma, fAnalytic);
+ }
+ else if (benchmark == "marmousi") {
+ u(37).addConstraint(gamma, 1.);
+ }
+
std::vector< std::vector< Subproblem * > > subproblem(nDomX);
for(unsigned int i = 0; i < nDomX; ++i) {
subproblem[i].resize(nDomY);
@@ -417,21 +430,6 @@ namespace D2 {
unsigned int index = i * nDomY + j;
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);
- }
- }
- else if (benchmark == "marmousi") {
- if(index == 37) {
- 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(- 1.), tf(lambda), gamma, gauss);
- }
- }
SubproblemDomains domains;
domains.setOmega(omega(index));
@@ -792,7 +790,7 @@ namespace D2 {
if(saveIteration) {
gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
- file << N << thetaPade << formulation.numberOfIterations() << gmshfem::csv::endl;
+ file << nDomX << nDomY << formulation.numberOfIterations() << gmshfem::csv::endl;
file.close();
}
@@ -847,18 +845,15 @@ namespace D2 {
cornerMono[3] = corner[3]((nDomX - 1) * nDomY);
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);
+ uMono.addConstraint(gamma, fAnalytic);
}
else if (benchmark == "marmousi") {
- formulationMono.integral(- 1., tf(lambda), gamma, gauss);
+ uMono.addConstraint(gamma, 1.);
}
SubproblemDomains domains;
diff --git a/examples/helmholtz/crossPoints/mesh.cpp b/examples/helmholtz/crossPoints/mesh.cpp
index 537d5a9a..37197a26 100755
--- a/examples/helmholtz/crossPoints/mesh.cpp
+++ b/examples/helmholtz/crossPoints/mesh.cpp
@@ -124,21 +124,34 @@ static void fixOrientationOfLine(std::vector< int > &cubeLines, const std::vecto
namespace D2 {
- void subdomain(const int index, const int jndex, const double xSize, const double ySize, const double circleSize, const double lc, const double pmlSizeE, const double pmlSizeN, const double pmlSizeW, const double pmlSizeS, bool transfinite)
+ void subdomain(const int index, const int jndex, const double xSize, const double ySize, const double circleSize, const double lc, const double pmlSizeE, const double pmlSizeN, const double pmlSizeW, const double pmlSizeS, bool transfinite, const std::vector< std::pair< double, double > > &points)
{
// square
double posX[4] = {xSize, xSize, 0., 0.};
double posY[4] = {0., ySize, ySize, 0.};
if(index != -1 && jndex != -1) {
- posX[0] = (index+1) * xSize;
- posX[1] = (index+1) * xSize;
- posX[2] = index * xSize;
- posX[3] = index * xSize;
-
- posY[0] = jndex * ySize;
- posY[1] = (jndex+1) * ySize;
- posY[2] = (jndex+1) * ySize;
- posY[3] = jndex * ySize;
+ if(points.empty()) {
+ posX[0] = (index+1) * xSize;
+ posX[1] = (index+1) * xSize;
+ posX[2] = index * xSize;
+ posX[3] = index * xSize;
+
+ posY[0] = jndex * ySize;
+ posY[1] = (jndex+1) * ySize;
+ posY[2] = (jndex+1) * ySize;
+ posY[3] = jndex * ySize;
+ }
+ else {
+ posX[0] = points[0].first;
+ posX[1] = points[1].first;
+ posX[2] = points[2].first;
+ posX[3] = points[3].first;
+
+ posY[0] = points[0].second;
+ posY[1] = points[1].second;
+ posY[2] = points[2].second;
+ posY[3] = points[3].second;
+ }
}
std::vector< int > squarePoints(4);
squarePoints[0] = newPoint(posX[0], posY[0], 0., lc);
@@ -354,6 +367,94 @@ namespace D2 {
gmsh::model::mesh::generate();
gmsh::model::mesh::setOrder(order);
}
+
+ void checkerboardTwist(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, const double circleSize, const double lc, const bool pml, const double pmlSize, const bool pmlExt, const double pmlSizeExt, const int order, const double twist, const unsigned int circlePosX, const unsigned int circlePosY, bool transfinite)
+ {
+ std::vector< std::vector< std::pair< double, double > > > grid(nDomX+1);
+ for(int i = 0; i <= nDomX; ++i) {
+ grid[i].resize(nDomY+1);
+ for(int j = 0; j <= nDomY; ++j) {
+ grid[i][j] = std::make_pair(i * xSize, j * ySize);
+ }
+ }
+ int count = 0;
+ for(int j = 1; j < nDomY; ++j) {
+ for(int i = 1; i < nDomX; ++i) {
+ switch(count % 4) {
+ case 0:
+ grid[i][j].first += twist * xSize;
+ break;
+ case 1:
+ grid[i][j].second -= twist * ySize;
+ break;
+ case 2:
+ grid[i][j].second += twist * ySize;
+ break;
+ case 3:
+ grid[i][j].first -= twist * xSize;
+ break;
+ }
+ ++count;
+ }
+ }
+
+ for(int i = 0; i < nDomX; ++i) {
+ for(int j = 0; j < nDomY; ++j) {
+ double pmlSizeE = 0.;
+ if(i == nDomX - 1) {
+ if(pmlExt) {
+ pmlSizeE = pmlSizeExt;
+ }
+ }
+ else {
+ if(pml) {
+ pmlSizeE = pmlSize;
+ }
+ }
+
+ double pmlSizeN = 0.;
+ if(j == nDomY - 1) {
+ if(pmlExt) {
+ pmlSizeN = pmlSizeExt;
+ }
+ }
+ else {
+ if(pml) {
+ pmlSizeN = pmlSize;
+ }
+ }
+
+ double pmlSizeW = 0.;
+ if(i == 0) {
+ if(pmlExt) {
+ pmlSizeW = pmlSizeExt;
+ }
+ }
+ else {
+ if(pml) {
+ pmlSizeW = pmlSize;
+ }
+ }
+
+ double pmlSizeS = 0.;
+ if(j == 0) {
+ if(pmlExt) {
+ pmlSizeS = pmlSizeExt;
+ }
+ }
+ else {
+ if(pml) {
+ pmlSizeS = pmlSize;
+ }
+ }
+
+ subdomain(i, j, xSize, ySize, (circlePosX == i && circlePosY == j ? circleSize : 0.), lc, pmlSizeE, pmlSizeN, pmlSizeW, pmlSizeS, transfinite, {grid[i+1][j], grid[i+1][j+1], grid[i][j+1], grid[i][j]});
+ }
+ }
+
+ gmsh::model::mesh::generate();
+ gmsh::model::mesh::setOrder(order);
+ }
}
diff --git a/examples/helmholtz/crossPoints/mesh.h b/examples/helmholtz/crossPoints/mesh.h
index 604f9b13..eb57dab2 100755
--- a/examples/helmholtz/crossPoints/mesh.h
+++ b/examples/helmholtz/crossPoints/mesh.h
@@ -5,9 +5,10 @@
namespace D2 {
- void subdomain(const int index, const int jndex, const double xSize, const double ySize, const double circleSize, const double lc, const double pmlSizeE, const double pmlSizeN, const double pmlSizeW, const double pmlSizeS, bool transfinite);
+ void subdomain(const int index, const int jndex, const double xSize, const double ySize, const double circleSize, const double lc, const double pmlSizeE, const double pmlSizeN, const double pmlSizeW, const double pmlSizeS, bool transfinite, const std::vector< std::pair< double, double > > &points = {});
void monoDomain(const double lc, const double pmlSize, const double R, const int order);
void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, const double circleSize, const double lc, const bool pml, const double pmlSize, const bool pmlExt, const double pmlSizeExt, const int order, const unsigned int circlePosX = 0, const unsigned int circlePosY = 0, bool transfinite = false);
+ void checkerboardTwist(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, const double circleSize, const double lc, const bool pml, const double pmlSize, const bool pmlExt, const double pmlSizeExt, const int order, const double twist, const unsigned int circlePosX = 0, const unsigned int circlePosY = 0, bool transfinite = false);
}
--
GitLab