diff --git a/examples/helmholtz/crossPoints/ddm2D.cpp b/examples/helmholtz/crossPoints/ddm2D.cpp
index ebaf8f8cf3604fc9f7c8100e6e5a5c2bc584a6a1..bf656f2958fcdb5bea229c4cd50b8c4f35af3e57 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 537d5a9aebf8065ab1e520097604b5f582400bea..37197a26f8f731b120da3ec7cc682c73536cdb7e 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 604f9b133051c775e20fa3e42266f37d39f62242..eb57dab2de6c2fad83fefd0d2873bd4954e03933 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);
}