diff --git a/CMakeLists.txt b/CMakeLists.txt index 75934f7c0f3025853c4c0bd44c1e358abb8f52ad..c598830685321a7c8a8a2587e25f819a1567d2d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,9 @@ cmake_minimum_required(VERSION 3.0 FATAL_ERROR) project(gmshddm CXX) +include(CTest) +enable_testing() + set(GMSHDDM_MAJOR_VERSION 1) set(GMSHDDM_MINOR_VERSION 0) set(GMSHDDM_PATCH_VERSION 0) @@ -142,6 +145,9 @@ add_subdirectory(src/problem) include_directories(src/common src/domain src/field src/problem) include_directories(contrib/eigen) +add_subdirectory(examples) +add_subdirectory(tutorials) + #---------------------------------- # # O T H E R L I B S diff --git a/demos/helmholtz/scattering/pml/main.cpp b/demos/helmholtz/scattering/pml/main.cpp deleted file mode 100644 index 33d688d38b64e131fcb9a93a596c0116ef8257ca..0000000000000000000000000000000000000000 --- a/demos/helmholtz/scattering/pml/main.cpp +++ /dev/null @@ -1,853 +0,0 @@ -#include "GmshDdm.h" - -#include "FormulationDDM.h" - -#include "Message.h" -#include "AnalyticalFunction.h" - -#include "gmsh.h" - -using namespace gmshddm; -using namespace gmshddm::common; -using namespace gmshddm::domain; -using namespace gmshddm::problem; -using namespace gmshddm::field; - -using namespace gmshfem; -using namespace gmshfem::analytics; -using namespace gmshfem::domain; -using namespace gmshfem::problem; -using namespace gmshfem::field; -using namespace gmshfem::function; -using namespace gmshfem::post; - -gmshfem::function::TensorFunction< std::complex< double > > getPmlDParameter(const double k, const double pmlSize, const double pos, const bool onX); -gmshfem::function::ScalarFunction< std::complex< double > > getPmlEParameter(const double k, const double pmlSize, const double pos, const bool onX); - -void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, const double circleSize, const double pmlSize, const double lc, const bool pml = true, const unsigned int circlePosX = 0, const unsigned int circlePosY = 0) -{ - std::vector< std::vector< int > > omega(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - omega[i].resize(nDomY); - } - - // ****************** - // borders - // ****************** - - // ****(2)**** - // * * - // (3) (1) - // * * - // ****(0)**** - std::vector< std::vector< std::vector< int > > > borders(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - borders[i].resize(nDomY); - for(unsigned int j = 0; j < nDomY; ++j) { - borders[i][j].resize(4); - } - } - - // ****************** - // pmls - // ****************** - - // 1 0 - // 0 ****(2)**** 1 - // * * - // (3) (1) - // * * - // 1 ****(0)**** 0 - // 0 1 - - std::vector< std::vector< std::vector< int > > > pmls(nDomX); - std::vector< std::vector< std::vector< int > > > pmlsInf(nDomX); - std::vector< std::vector< std::vector< std::vector< int > > > > pmlsBnd(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - pmls[i].resize(nDomY); - pmlsInf[i].resize(nDomY); - pmlsBnd[i].resize(nDomY); - for(unsigned int j = 0; j < nDomY; ++j) { - pmls[i][j].resize(4); - pmlsInf[i][j].resize(4); - pmlsBnd[i][j].resize(4); - for(unsigned int k = 0; k < 4; ++k) { - pmlsBnd[i][j][k].resize(2); - } - } - } - - // ****************** - // pmlsSquare - // ****************** - - // (2)*******(1) - // * * - // * * - // * * - // (3)*******(0) - - std::vector< std::vector< std::vector< int > > > pmlsSquare(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - pmlsSquare[i].resize(nDomY); - for(unsigned int j = 0; j < nDomY; ++j) { - pmlsSquare[i][j].resize(4); - } - } - - std::vector< int > left(nDomY); - int bottom = 0; - for(unsigned int i = 0; i < nDomX; ++i) { - for(unsigned int j = 0; j < nDomY; ++j) { - int p[4]; - p[0] = gmsh::model::geo::addPoint(i * xSize, j * ySize, 0., lc); - p[1] = gmsh::model::geo::addPoint((i + 1) * xSize, j * ySize, 0., lc); - p[2] = gmsh::model::geo::addPoint((i + 1) * xSize, (j + 1) * ySize, 0., lc); - p[3] = gmsh::model::geo::addPoint(i * xSize, (j + 1) * ySize, 0., lc); - - int l[4]; - borders[i][j][0] = l[0] = gmsh::model::geo::addLine(p[0], p[1]); - borders[i][j][1] = l[1] = gmsh::model::geo::addLine(p[1], p[2]); - borders[i][j][2] = l[2] = gmsh::model::geo::addLine(p[2], p[3]); - borders[i][j][3] = l[3] = gmsh::model::geo::addLine(p[3], p[0]); - - if(i == 0) { - left[j] = l[1]; - } - else { - borders[i][j][3] = left[j]; - left[j] = l[1]; - } - - if(j == 0) { - bottom = l[2]; - } - else { - borders[i][j][0] = bottom; - bottom = l[2]; - } - - int ll = gmsh::model::geo::addCurveLoop({l[0], l[1], l[2], l[3]}); - - int surface = 0; - if(i == circlePosX && j == circlePosY) { - // circle - int pC[5]; - pC[0] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2., circlePosY * j + ySize / 2., 0., lc); - pC[1] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2. + circleSize, circlePosY * j + ySize / 2., 0., lc); - pC[2] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2., circlePosY * j + ySize / 2. + circleSize, 0., lc); - pC[3] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2. - circleSize, circlePosY * j + ySize / 2., 0., lc); - pC[4] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2., circlePosY * j + ySize / 2. - circleSize, 0., lc); - - int lC[4]; - lC[0] = gmsh::model::geo::addCircleArc(pC[1], pC[0], pC[2]); - lC[1] = gmsh::model::geo::addCircleArc(pC[2], pC[0], pC[3]); - lC[2] = gmsh::model::geo::addCircleArc(pC[3], pC[0], pC[4]); - lC[3] = gmsh::model::geo::addCircleArc(pC[4], pC[0], pC[1]); - - // surface - int llC; - llC = gmsh::model::geo::addCurveLoop({lC[0], lC[1], lC[2], lC[3]}); - - surface = gmsh::model::geo::addPlaneSurface({ll, llC}); - - gmsh::model::addPhysicalGroup(1, {lC[0], lC[1], lC[2], lC[3]}, 1); - gmsh::model::setPhysicalName(1, 1, "gammaScat"); - } - else { - surface = gmsh::model::geo::addPlaneSurface({ll}); - } - omega[i][j] = surface; - - if(pml) { - // pmls - double extrudeX[4] = {0., pmlSize, 0., -pmlSize}; - double extrudeY[4] = {-pmlSize, 0., pmlSize, 0.}; - for(unsigned int k = 0; k < 4; ++k) { - std::vector< std::pair< int, int > > extrudeEntities; - gmsh::model::geo::extrude({std::make_pair(1, ((k == 3 && i != 0) || (k == 0 && j != 0) ? 1 : -1) * borders[i][j][k])}, extrudeX[k], extrudeY[k], 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true); - pmlsInf[i][j][k] = extrudeEntities[0].second; - pmls[i][j][k] = extrudeEntities[1].second; - if((k == 3 && i != 0) || (k == 0 && j != 0)) { - pmlsBnd[i][j][k][0] = extrudeEntities[2].second; - pmlsBnd[i][j][k][1] = extrudeEntities[3].second; - } - else { - pmlsBnd[i][j][k][0] = extrudeEntities[3].second; - pmlsBnd[i][j][k][1] = extrudeEntities[4].second; - } - } - - // pmlsSquare - double extrudeSquareX[4] = {pmlSize, pmlSize, -pmlSize, -pmlSize}; - int side[4] = {0, 2, 2, 0}; - int sideBnd[4] = {1, 0, 1, 0}; - for(unsigned int k = 0; k < 4; ++k) { - std::vector< std::pair< int, int > > extrudeEntities; - gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][side[k]][sideBnd[k]])}, extrudeSquareX[k], 0., 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true); - pmlsSquare[i][j][k] = extrudeEntities[1].second; - } - } - } - } - - gmsh::model::geo::synchronize(); - - // ****************** - // physicals - // ****************** - - for(unsigned int i = 0; i < nDomX; ++i) { - for(unsigned int j = 0; j < nDomY; ++j) { - gmsh::model::addPhysicalGroup(2, {omega[i][j]}, i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, i * nDomY + j + 1, "omega_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][0]}, 10000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 10000 + i * nDomY + j + 1, "sigmaS_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][1]}, 20000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 20000 + i * nDomY + j + 1, "sigmaE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][2]}, 30000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 30000 + i * nDomY + j + 1, "sigmaN_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][3]}, 40000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 40000 + i * nDomY + j + 1, "sigmaW_" + std::to_string(i) + "_" + std::to_string(j)); - - if(pml) { - // PMLs - gmsh::model::addPhysicalGroup(2, {pmls[i][j][0]}, 10000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 10000 + i * nDomY + j + 1, "pmlS_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmls[i][j][1]}, 20000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 20000 + i * nDomY + j + 1, "pmlE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmls[i][j][2]}, 30000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 30000 + i * nDomY + j + 1, "pmlN_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmls[i][j][3]}, 40000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 40000 + i * nDomY + j + 1, "pmlW_" + std::to_string(i) + "_" + std::to_string(j)); - - // PMLs square - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][0]}, 100000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 100000 + i * nDomY + j + 1, "pmlSE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][1]}, 200000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 200000 + i * nDomY + j + 1, "pmlNE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][2]}, 300000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 300000 + i * nDomY + j + 1, "pmlNW_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][3]}, 400000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 400000 + i * nDomY + j + 1, "pmlSW_" + std::to_string(i) + "_" + std::to_string(j)); - - // Sigmas in PMLs - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3][1]}, 100000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 100000 + i * nDomY + j + 1, "sigmaS_W_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1][0]}, 150000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 150000 + i * nDomY + j + 1, "sigmaS_E_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0][1]}, 200000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 200000 + i * nDomY + j + 1, "sigmaE_S_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2][0]}, 250000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 250000 + i * nDomY + j + 1, "sigmaE_N_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1][1]}, 300000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 300000 + i * nDomY + j + 1, "sigmaN_E_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3][0]}, 350000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 350000 + i * nDomY + j + 1, "sigmaN_W_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2][1]}, 400000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 400000 + i * nDomY + j + 1, "sigmaW_N_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0][0]}, 450000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 450000 + i * nDomY + j + 1, "sigmaW_S_" + std::to_string(i) + "_" + std::to_string(j)); - - } - } - } - - gmsh::model::addPhysicalGroup(0, {3}, 1); - gmsh::model::setPhysicalName(0, 1, "crossPoint"); - - gmsh::model::geo::synchronize(); - gmsh::model::mesh::generate(); -} - -int main(int argc, char **argv) -{ - GmshDdm gmshDdm(argc, argv); - - double pi = 3.14159265359; - - int nDomX = 2, nDomY = 2; - gmshDdm.userDefinedParameter(nDomX, "nDomX"); - gmshDdm.userDefinedParameter(nDomY, "nDomY"); - double xSize = 1., ySize = 1.; - gmshDdm.userDefinedParameter(xSize, "xSize"); - gmshDdm.userDefinedParameter(ySize, "ySize"); - double circleSize = 0.1; - gmshDdm.userDefinedParameter(circleSize, "circleSize"); - double pmlSize = 0.1; - gmshDdm.userDefinedParameter(pmlSize, "pmlSize"); - double k = 10. * pi + 0.3; - gmshDdm.userDefinedParameter(k, "k"); - double lc = 0.02; - gmshDdm.userDefinedParameter(lc, "lc"); - std::string gauss = "Gauss3"; - gmshDdm.userDefinedParameter(gauss, "gauss"); - int circlePosX = 0, circlePosY = 0; - gmshDdm.userDefinedParameter(circlePosX, "circlePosX"); - gmshDdm.userDefinedParameter(circlePosY, "circlePosY"); - bool errorMono = false; - gmshDdm.userDefinedParameter(errorMono, "errorMono"); - bool withSquarePML = true; - gmshDdm.userDefinedParameter(withSquarePML, "withSquarePML"); - - // Build geometry and mesh using gmsh API - checkerboard(nDomX, nDomY, xSize, ySize, circleSize, pmlSize, lc, true, circlePosX, circlePosY); - - unsigned int nDom = nDomX * nDomY; - - // | | - // NW | N | NE - // ----**************---- - // * * - // * * - // W * * E - // * * - // * * - // ----**************---- - // SW | S | SE - // | | - - // Define domain - Subdomain omega(nDom); - Subdomain gammaDir(nDom); - - Subdomain pmlS(nDom); - Subdomain pmlE(nDom); - Subdomain pmlN(nDom); - Subdomain pmlW(nDom); - - Subdomain pmlSE(nDom); - Subdomain pmlNE(nDom); - Subdomain pmlNW(nDom); - Subdomain pmlSW(nDom); - - InterfaceDomain sigmaS(nDom); - InterfaceDomain sigmaE(nDom); - InterfaceDomain sigmaN(nDom); - InterfaceDomain sigmaW(nDom); - - InterfaceDomain sigmaS_W(nDom); - InterfaceDomain sigmaS_E(nDom); - InterfaceDomain sigmaE_S(nDom); - InterfaceDomain sigmaE_N(nDom); - InterfaceDomain sigmaN_E(nDom); - InterfaceDomain sigmaN_W(nDom); - InterfaceDomain sigmaW_N(nDom); - InterfaceDomain sigmaW_S(nDom); - - InterfaceDomain pmlS_W(nDom); - InterfaceDomain pmlS_E(nDom); - InterfaceDomain pmlE_S(nDom); - InterfaceDomain pmlE_N(nDom); - InterfaceDomain pmlN_E(nDom); - InterfaceDomain pmlN_W(nDom); - InterfaceDomain pmlW_N(nDom); - InterfaceDomain pmlW_S(nDom); - - InterfaceDomain crossPoint(nDom); - - // Define topology - std::vector< std::vector< unsigned int > > topology(nDom); - for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) { - for(unsigned int j = 0; j < static_cast< unsigned int >(nDomY); ++j) { - unsigned int index = i * nDomY + j; - - omega(index) = Domain(2, index + 1); - - pmlS(index) = Domain(2, 10000 + index + 1); - pmlE(index) = Domain(2, 20000 + index + 1); - pmlN(index) = Domain(2, 30000 + index + 1); - pmlW(index) = Domain(2, 40000 + index + 1); - - if(withSquarePML || i == static_cast< unsigned int >(nDomX) - 1 || j == 0) { - pmlSE(index) = Domain(2, 100000 + index + 1); - } - - if(withSquarePML || i == static_cast< unsigned int >(nDomX) - 1 || j == static_cast< unsigned int >(nDomY) - 1) { - pmlNE(index) = Domain(2, 200000 + index + 1); - } - - if(withSquarePML || i == 0 || j == static_cast< unsigned int >(nDomY) - 1) { - pmlNW(index) = Domain(2, 300000 + index + 1); - } - - if(withSquarePML || i == 0 || j == 0) { - pmlSW(index) = Domain(2, 400000 + index + 1); - } - - if(i == circlePosX && j == circlePosY) { - gammaDir(index) = Domain(1, 1); - } - - if(i != 0) { - sigmaW(index, (i-1) * nDomY + j) = Domain(1, 40000 + i * nDomY + j + 1); - - crossPoint(index, (i-1) * nDomY + j) = Domain(0, 1); - - if(withSquarePML || j == static_cast< unsigned int >(nDomY) - 1) sigmaW_N(index, (i-1) * nDomY + j) = Domain(1, 400000 + i * nDomY + j + 1); - if(withSquarePML || j == 0) sigmaW_S(index, (i-1) * nDomY + j) = Domain(1, 450000 + i * nDomY + j + 1); - - if(j != static_cast< unsigned int >(nDomY) - 1) pmlW_N(index, (i-1) * nDomY + j) = Domain(1, 350000 + i * nDomY + j + 1); - if(j != 0) pmlW_S(index, (i-1) * nDomY + j) = Domain(1, 100000 + i * nDomY + j + 1); - - topology[index].push_back((i-1) * nDomY + j); - } - if(i != static_cast< unsigned int >(nDomX) - 1) { - sigmaE(index, (i+1) * nDomY + j) = Domain(1, 20000 + i * nDomY + j + 1); - - crossPoint(index, (i+1) * nDomY + j) = Domain(0, 1); - - if(withSquarePML || j == 0) sigmaE_S(index, (i+1) * nDomY + j) = Domain(1, 200000 + i * nDomY + j + 1); - if(withSquarePML || j == static_cast< unsigned int >(nDomY) - 1) sigmaE_N(index, (i+1) * nDomY + j) = Domain(1, 250000 + i * nDomY + j + 1); - - if(j != 0) pmlE_S(index, (i+1) * nDomY + j) = Domain(1, 150000 + i * nDomY + j + 1); - if(j != static_cast< unsigned int >(nDomY) - 1) pmlE_N(index, (i+1) * nDomY + j) = Domain(1, 300000 + i * nDomY + j + 1); - - topology[index].push_back((i+1) * nDomY + j); - } - - if(j != 0) { - sigmaS(index, i * nDomY + (j-1)) = Domain(1, 10000 + i * nDomY + j + 1); - - crossPoint(index, i * nDomY + (j-1)) = Domain(0, 1); - - if(withSquarePML || i == 0) sigmaS_W(index, i * nDomY + (j-1)) = Domain(1, 100000 + i * nDomY + j + 1); - if(withSquarePML || i == static_cast< unsigned int >(nDomX) - 1) sigmaS_E(index, i * nDomY + (j-1)) = Domain(1, 150000 + i * nDomY + j + 1); - - if(i != 0) pmlS_W(index, i * nDomY + (j-1)) = Domain(1, 450000 + i * nDomY + j + 1); - if(i != static_cast< unsigned int >(nDomX) - 1) pmlS_E(index, i * nDomY + (j-1)) = Domain(1, 200000 + i * nDomY + j + 1); - - topology[index].push_back(i * nDomY + (j-1)); - } - if(j != static_cast< unsigned int >(nDomY) - 1) { - sigmaN(index, i * nDomY + (j+1)) = Domain(1, 30000 + i * nDomY + j + 1); - - crossPoint(index, i * nDomY + (j+1)) = Domain(0, 1); - - if(withSquarePML || i == static_cast< unsigned int >(nDomX) - 1) sigmaN_E(index, i * nDomY + (j+1)) = Domain(1, 300000 + i * nDomY + j + 1); - if(withSquarePML || i == 0) sigmaN_W(index, i * nDomY + (j+1)) = Domain(1, 350000 + i * nDomY + j + 1); - - if(i != static_cast< unsigned int >(nDomX) - 1) pmlN_E(index, i * nDomY + (j+1)) = Domain(1, 250000 + i * nDomY + j + 1); - if(i != 0) pmlN_W(index, i * nDomY + (j+1)) = Domain(1, 400000 + i * nDomY + j + 1); - - topology[index].push_back(i * nDomY + (j+1)); - } - } - } - - // Create resolution - Resolution< std::complex< double > > resolution("DDM"); - - // Create DDM formulation - FormulationDDM< std::complex< double > > formulation("Helmholtz", resolution, topology); - - // Create fields - SubdomainField< std::complex< double >, Form::Form0 > u(resolution, "u", omega | gammaDir | pmlS | pmlE | pmlN | pmlW | pmlSE | pmlNE | pmlNW | pmlSW, FunctionSpaceType::Lagrange); - SubdomainField< std::complex< double >, Form::Form0 > lambda(resolution, "lambda", gammaDir, FunctionSpaceType::Lagrange); - - InterfaceField< std::complex< double >, Form::Form0 > gS(resolution, "gS", sigmaS | sigmaS_W | sigmaS_E, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > &gS_W = gS; - InterfaceField< std::complex< double >, Form::Form0 > &gS_E = gS; -// InterfaceField< std::complex< double >, Form::Form0 > gS_W(resolution, "gS_W", sigmaS_W, FunctionSpaceType::Lagrange); -// InterfaceField< std::complex< double >, Form::Form0 > gS_E(resolution, "gS_E", sigmaS_E, FunctionSpaceType::Lagrange); - - InterfaceField< std::complex< double >, Form::Form0 > gE(resolution, "gE", sigmaE | sigmaE_S | sigmaE_N, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > &gE_S = gE; - InterfaceField< std::complex< double >, Form::Form0 > &gE_N = gE; -// InterfaceField< std::complex< double >, Form::Form0 > gE_S(resolution, "gE_S", sigmaE_S, FunctionSpaceType::Lagrange); -// InterfaceField< std::complex< double >, Form::Form0 > gE_N(resolution, "gE_N", sigmaE_N, FunctionSpaceType::Lagrange); - - InterfaceField< std::complex< double >, Form::Form0 > gN(resolution, "gN", sigmaN | sigmaN_E | sigmaN_W, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > &gN_E = gN; - InterfaceField< std::complex< double >, Form::Form0 > &gN_W = gN; -// InterfaceField< std::complex< double >, Form::Form0 > gN_E(resolution, "gN_E", sigmaN_E, FunctionSpaceType::Lagrange); -// InterfaceField< std::complex< double >, Form::Form0 > gN_W(resolution, "gN_W", sigmaN_W, FunctionSpaceType::Lagrange); - - InterfaceField< std::complex< double >, Form::Form0 > gW(resolution, "gW", sigmaW | sigmaW_N | sigmaW_S, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > &gW_N = gW; - InterfaceField< std::complex< double >, Form::Form0 > &gW_S = gW; -// InterfaceField< std::complex< double >, Form::Form0 > gW_N(resolution, "gW_N", sigmaW_N, FunctionSpaceType::Lagrange); -// InterfaceField< std::complex< double >, Form::Form0 > gW_S(resolution, "gW_S", sigmaW_S, FunctionSpaceType::Lagrange); - - std::complex< double > im(0., 1.); - AnalyticalFunction< helmholtz2D::SoftScatteringByACylinder< std::complex< double > > > solution(k, circleSize, xSize / 2., ySize / 2., 2 * k); - - // Tell to the formulation that g is field that have to be exchanged between subdomains - formulation.addInterfaceField(gS); -// formulation.addInterfaceField(gS_W); -// formulation.addInterfaceField(gS_E); - - formulation.addInterfaceField(gE); -// formulation.addInterfaceField(gE_S); -// formulation.addInterfaceField(gE_N); - - formulation.addInterfaceField(gN); -// formulation.addInterfaceField(gN_E); -// formulation.addInterfaceField(gN_W); - - formulation.addInterfaceField(gW); -// formulation.addInterfaceField(gW_N); -// formulation.addInterfaceField(gW_S); - - // Add terms to the formulation - for(unsigned int i = 0; i < nDom; ++i) { - // Compute PML parameters - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaS = abs(y< std::complex< double > >() - (y< std::complex< double > >(omega(i)) - ySize / 2.)); - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaE = abs(x< std::complex< double > >() - (x< std::complex< double > >(omega(i)) + xSize / 2.)); - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaN = abs(y< std::complex< double > >() - (y< std::complex< double > >(omega(i)) + ySize / 2.)); - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaW = abs(x< std::complex< double > >() - (x< std::complex< double > >(omega(i)) - xSize / 2.)); - - gmshfem::function::ScalarFunction< std::complex< double > > sigS = 1. / (pmlSize - distSigmaS) - 1./ pmlSize; - gmshfem::function::ScalarFunction< std::complex< double > > sigE = 1. / (pmlSize - distSigmaE) - 1./ pmlSize; - gmshfem::function::ScalarFunction< std::complex< double > > sigN = 1. / (pmlSize - distSigmaN) - 1./ pmlSize; - gmshfem::function::ScalarFunction< std::complex< double > > sigW = 1. / (pmlSize - distSigmaW) - 1./ pmlSize; - - gmshfem::function::ScalarFunction< std::complex< double > > kyS = 1. + im * sigS / k; - gmshfem::function::ScalarFunction< std::complex< double > > kxE = 1. + im * sigE / k; - gmshfem::function::ScalarFunction< std::complex< double > > kyN = 1. + im * sigN / k; - gmshfem::function::ScalarFunction< std::complex< double > > kxW = 1. + im * sigW / k; - - gmshfem::function::TensorFunction< std::complex< double > > DS = tensorDiag< std::complex< double > >(kyS, 1. / kyS, kyS); - gmshfem::function::ScalarFunction< std::complex< double > > ES = kyS; - - gmshfem::function::TensorFunction< std::complex< double > > DE = tensorDiag< std::complex< double > >(1. / kxE, kxE, kxE); - gmshfem::function::ScalarFunction< std::complex< double > > EE = kxE; - - gmshfem::function::TensorFunction< std::complex< double > > DN = tensorDiag< std::complex< double > >(kyN, 1. / kyN, kyN); - gmshfem::function::ScalarFunction< std::complex< double > > EN = kyN; - - gmshfem::function::TensorFunction< std::complex< double > > DW = tensorDiag< std::complex< double > >(1. / kxW, kxW, kxW); - gmshfem::function::ScalarFunction< std::complex< double > > EW = kxW; - - // VOLUME TERMS - - formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(- k * k * dof(u(i)), tf(u(i)), omega(i), gauss); - - // Physical source - formulation(i).integral(dof(lambda(i)), tf(u(i)), gammaDir(i), gauss); - formulation(i).integral(dof(u(i)), tf(lambda(i)), gammaDir(i), gauss); - formulation(i).integral(formulation.physicalSource(-solution), tf(lambda(i)), gammaDir(i), gauss); - - // S - formulation(i).integral(DS * grad(dof(u(i))), grad(tf(u(i))), pmlS(i), gauss); - formulation(i).integral(- ES * k * k * dof(u(i)), tf(u(i)), pmlS(i), gauss); - - // E - formulation(i).integral(DE * grad(dof(u(i))), grad(tf(u(i))), pmlE(i), gauss); - formulation(i).integral(- EE * k * k * dof(u(i)), tf(u(i)), pmlE(i), gauss); - - // N - formulation(i).integral(DN * grad(dof(u(i))), grad(tf(u(i))), pmlN(i), gauss); - formulation(i).integral(- EN * k * k * dof(u(i)), tf(u(i)), pmlN(i), gauss); - - // W - formulation(i).integral(DW * grad(dof(u(i))), grad(tf(u(i))), pmlW(i), gauss); - formulation(i).integral(- EW * k * k * dof(u(i)), tf(u(i)), pmlW(i), gauss); - - // SE - formulation(i).integral(DS * DE * grad(dof(u(i))), grad(tf(u(i))), pmlSE(i), gauss); - formulation(i).integral(- ES * EE * k * k * dof(u(i)), tf(u(i)), pmlSE(i), gauss); - - // NE - formulation(i).integral(DN * DE * grad(dof(u(i))), grad(tf(u(i))), pmlNE(i), gauss); - formulation(i).integral(- EN * EE * k * k * dof(u(i)), tf(u(i)), pmlNE(i), gauss); - - // NW - formulation(i).integral(DN * DW * grad(dof(u(i))), grad(tf(u(i))), pmlNW(i), gauss); - formulation(i).integral(- EN * EW * k * k * dof(u(i)), tf(u(i)), pmlNW(i), gauss); - - // SW - formulation(i).integral(DS * DW * grad(dof(u(i))), grad(tf(u(i))), pmlSW(i), gauss); - formulation(i).integral(- ES * EW * k * k * dof(u(i)), tf(u(i)), pmlSW(i), gauss); - - // Artificial sources - for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { - const unsigned int j = topology[i][jj]; - if(!sigmaS(i,j).isEmpty()) { - formulation(i).integral(formulation.artificialSource( -gN(j,i) ), tf(u(i)), sigmaS(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gN_W(j,i) ), tf(u(i)), sigmaS_W(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gN_E(j,i) ), tf(u(i)), sigmaS_E(i,j), gauss); - } - if(!sigmaE(i,j).isEmpty()) { - formulation(i).integral(formulation.artificialSource( -gW(j,i) ), tf(u(i)), sigmaE(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gW_S(j,i) ), tf(u(i)), sigmaE_S(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gW_N(j,i) ), tf(u(i)), sigmaE_N(i,j), gauss); - } - if(!sigmaN(i,j).isEmpty()) { - formulation(i).integral(formulation.artificialSource( -gS(j,i) ), tf(u(i)), sigmaN(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gS_E(j,i) ), tf(u(i)), sigmaN_E(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gS_W(j,i) ), tf(u(i)), sigmaN_W(i,j), gauss); - } - if(!sigmaW(i,j).isEmpty()) { - formulation(i).integral(formulation.artificialSource( -gE(j,i) ), tf(u(i)), sigmaW(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gE_N(j,i) ), tf(u(i)), sigmaW_N(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gE_S(j,i) ), tf(u(i)), sigmaW_S(i,j), gauss); - } - } - - // INTERFACE TERMS - for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { - const unsigned int j = topology[i][jj]; - // S - if(!sigmaS(i,j).isEmpty()) { - formulation(i,j).integral(dof(gS(i,j)), tf(gS(i,j)), sigmaS(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gN(j,i) ), tf(gS(i,j)), sigmaS(i,j), gauss); - - formulation(i,j).integral(dof(gS_W(i,j)), tf(gS_W(i,j)), sigmaS_W(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gN_W(j,i) ), tf(gS_W(i,j)), sigmaS_W(i,j), gauss); - - formulation(i,j).integral(dof(gS_E(i,j)), tf(gS_E(i,j)), sigmaS_E(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gN_E(j,i) ), tf(gS_E(i,j)), sigmaS_E(i,j), gauss); - - if(!pmlS_W(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gE_S(i - nDomY,i) ), tf(gS(i,j)), pmlS_W(i,j), gauss); - if(!pmlS_E(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gW_S(i + nDomY,i) ), tf(gS(i,j)), pmlS_E(i,j), gauss); - - formulation(i,j).integral(- 2. * DS * grad(u(i)), grad(tf(gS(i,j))), pmlS(i), gauss); - formulation(i,j).integral(2. * ES * k * k * u(i), tf(gS(i,j)), pmlS(i), gauss); - - formulation(i,j).integral(- 2. * DS * DW * grad(u(i)), grad(tf(gS_W(i,j))), pmlSW(i), gauss); - formulation(i,j).integral(2. * ES * EW * k * k * u(i), tf(gS_W(i,j)), pmlSW(i), gauss); - - formulation(i,j).integral(- 2. * DS * DE * grad(u(i)), grad(tf(gS_E(i,j))), pmlSE(i), gauss); - formulation(i,j).integral(2. * ES * EE * k * k * u(i), tf(gS_E(i,j)), pmlSE(i), gauss); - } - - // E - if(!sigmaE(i,j).isEmpty()) { - formulation(i,j).integral(dof(gE(i,j)), tf(gE(i,j)), sigmaE(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gW(j,i) ), tf(gE(i,j)), sigmaE(i,j), gauss); - - formulation(i,j).integral(dof(gE_S(i,j)), tf(gE_S(i,j)), sigmaE_S(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gW_S(j,i) ), tf(gE_S(i,j)), sigmaE_S(i,j), gauss); - - formulation(i,j).integral(dof(gE_N(i,j)), tf(gE_N(i,j)), sigmaE_N(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gW_N(j,i) ), tf(gE_N(i,j)), sigmaE_N(i,j), gauss); - - if(!pmlE_S(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gN_E(i - 1,i) ), tf(gE(i,j)), pmlE_S(i,j), gauss); - if(!pmlE_N(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gS_E(i + 1,i) ), tf(gE(i,j)), pmlE_N(i,j), gauss); - - formulation(i,j).integral(- 2. * DE * grad(u(i)), grad(tf(gE(i,j))), pmlE(i), gauss); - formulation(i,j).integral(2. * EE * k * k * u(i), tf(gE(i,j)), pmlE(i), gauss); - - formulation(i,j).integral(- 2. * DS * DE * grad(u(i)), grad(tf(gE_S(i,j))), pmlSE(i), gauss); - formulation(i,j).integral(2. * ES * EE * k * k * u(i), tf(gE_S(i,j)), pmlSE(i), gauss); - - formulation(i,j).integral(- 2. * DN * DE * grad(u(i)), grad(tf(gE_N(i,j))), pmlNE(i), gauss); - formulation(i,j).integral(2. * EN * EE * k * k * u(i), tf(gE_N(i,j)), pmlNE(i), gauss); - } - - // N - if(!sigmaN(i,j).isEmpty()) { - formulation(i,j).integral(dof(gN(i,j)), tf(gN(i,j)), sigmaN(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gS(j,i) ), tf(gN(i,j)), sigmaN(i,j), gauss); - - formulation(i,j).integral(dof(gN_E(i,j)), tf(gN_E(i,j)), sigmaN_E(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gS_E(j,i) ), tf(gN_E(i,j)), sigmaN_E(i,j), gauss); - - formulation(i,j).integral(dof(gN_W(i,j)), tf(gN_W(i,j)), sigmaN_W(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gS_W(j,i) ), tf(gN_W(i,j)), sigmaN_W(i,j), gauss); - - if(!pmlN_E(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gW_N(i + nDomY,i) ), tf(gN(i,j)), pmlN_E(i,j), gauss); - if(!pmlN_W(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gE_N(i - nDomY,i) ), tf(gN(i,j)), pmlN_W(i,j), gauss); - - formulation(i,j).integral(- 2. * DN * grad(u(i)), grad(tf(gN(i,j))), pmlN(i), gauss); - formulation(i,j).integral(2. * EN * k * k * u(i), tf(gN(i,j)), pmlN(i), gauss); - - formulation(i,j).integral(- 2. * DN * DE * grad(u(i)), grad(tf(gN_E(i,j))), pmlNE(i), gauss); - formulation(i,j).integral(2. * EN * EE * k * k * u(i), tf(gN_E(i,j)), pmlNE(i), gauss); - - formulation(i,j).integral(- 2. * DN * DW * grad(u(i)), grad(tf(gN_W(i,j))), pmlNW(i), gauss); - formulation(i,j).integral(2. * EN * EW * k * k * u(i), tf(gN_W(i,j)), pmlNW(i), gauss); - } - - // W - if(!sigmaW(i,j).isEmpty()) { - formulation(i,j).integral(dof(gW(i,j)), tf(gW(i,j)), sigmaW(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gE(j,i) ), tf(gW(i,j)), sigmaW(i,j), gauss); - - formulation(i,j).integral(dof(gW_N(i,j)), tf(gW_N(i,j)), sigmaW_N(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gE_N(j,i) ), tf(gW_N(i,j)), sigmaW_N(i,j), gauss); - - formulation(i,j).integral(dof(gW_S(i,j)), tf(gW_S(i,j)), sigmaW_S(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gE_S(j,i) ), tf(gW_S(i,j)), sigmaW_S(i,j), gauss); - - if(!pmlW_N(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gS_W(i + 1,i) ), tf(gW(i,j)), pmlW_N(i,j), gauss); - if(!pmlW_S(i,j).isEmpty() && withSquarePML) formulation(i,j).integral(formulation.artificialSource( 2. * gN_W(i - 1,i) ), tf(gW(i,j)), pmlW_S(i,j), gauss); - - formulation(i,j).integral(- 2. * DW * grad(u(i)), grad(tf(gW(i,j))), pmlW(i), gauss); - formulation(i,j).integral(2. * EW * k * k * u(i), tf(gW(i,j)), pmlW(i), gauss); - - formulation(i,j).integral(- 2. * DN * DW * grad(u(i)), grad(tf(gW_N(i,j))), pmlNW(i), gauss); - formulation(i,j).integral(2. * EN * EW * k * k * u(i), tf(gW_N(i,j)), pmlNW(i), gauss); - - formulation(i,j).integral(- 2. * DS * DW * grad(u(i)), grad(tf(gW_S(i,j))), pmlSW(i), gauss); - formulation(i,j).integral(2. * ES * EW * k * k * u(i), tf(gW_S(i,j)), pmlSW(i), gauss); - } - } - } - - // Solve the DDM formulation - formulation.pre(); - formulation.solve("gmres", 1e-6, 1000); - - if(errorMono) { - gmshfem::domain::Domain omegaMono, pmlSMono, pmlEMono, pmlNMono, pmlWMono, pmlSEMono, pmlNEMono, pmlNWMono, pmlSWMono; - for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { - omegaMono |= omega(i); - } - - for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) { - pmlSMono |= pmlS(i * nDomY + 0); - pmlNMono |= pmlN(i * nDomY + nDomY - 1); - } - for(unsigned int j = 0; j < static_cast< unsigned int >(nDomY); ++j) { - pmlEMono |= pmlE((nDomX - 1) * nDomY + j); - pmlWMono |= pmlW(0 * nDomY + j); - } - pmlSEMono = pmlSE((nDomX - 1) * nDomY + 0); - pmlNEMono = pmlNE((nDomX - 1) * nDomY + nDomY - 1); - pmlNWMono = pmlNW(0 * nDomY + nDomY - 1); - pmlSWMono = pmlSW(0); - - gmshfem::function::TensorFunction< std::complex< double > > DSmono = getPmlDParameter(k, pmlSize, 0., false); - gmshfem::function::TensorFunction< std::complex< double > > DEmono = getPmlDParameter(k, pmlSize, nDomX * xSize, true); - gmshfem::function::TensorFunction< std::complex< double > > DNmono = getPmlDParameter(k, pmlSize, nDomY * ySize, false); - gmshfem::function::TensorFunction< std::complex< double > > DWmono = getPmlDParameter(k, pmlSize, 0., true); - - gmshfem::function::ScalarFunction< std::complex< double > > ESmono = getPmlEParameter(k, pmlSize, 0., false); - gmshfem::function::ScalarFunction< std::complex< double > > EEmono = getPmlEParameter(k, pmlSize, nDomX * xSize, true); - gmshfem::function::ScalarFunction< std::complex< double > > ENmono = getPmlEParameter(k, pmlSize, nDomY * ySize, false); - gmshfem::function::ScalarFunction< std::complex< double > > EWmono = getPmlEParameter(k, pmlSize, 0., true); - - gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > uMono(resolution, "uMono", gammaDir(circlePosX * nDomY + circlePosY) | omegaMono | pmlSMono | pmlEMono | pmlNMono | pmlWMono | pmlSEMono | pmlNEMono | pmlNWMono | pmlSWMono, FunctionSpaceType::Lagrange); - gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambdaMono(resolution, "lambdaMono", gammaDir(circlePosX * nDomY + circlePosY), FunctionSpaceType::Lagrange); - - Formulation< std::complex< double > > monodomain("HelmholtzMono", resolution); - - monodomain.integral(grad(dof(uMono)), grad(tf(uMono)), omegaMono, gauss); - monodomain.integral(- k * k * dof(uMono), tf(uMono), omegaMono, gauss); - - monodomain.integral(dof(lambdaMono), tf(uMono), gammaDir(circlePosX * nDomY + circlePosY), gauss); - monodomain.integral(dof(uMono), tf(lambdaMono), gammaDir(circlePosX * nDomY + circlePosY), gauss); - monodomain.integral(- solution, tf(lambdaMono), gammaDir(circlePosX * nDomY + circlePosY), gauss); - - monodomain.integral(DSmono * grad(dof(uMono)), grad(tf(uMono)), pmlSMono, gauss); - monodomain.integral(- k * k * ESmono * dof(uMono), tf(uMono), pmlSMono, gauss); - - monodomain.integral(DEmono * grad(dof(uMono)), grad(tf(uMono)), pmlEMono, gauss); - monodomain.integral(- k * k * EEmono * dof(uMono), tf(uMono), pmlEMono, gauss); - - monodomain.integral(DNmono * grad(dof(uMono)), grad(tf(uMono)), pmlNMono, gauss); - monodomain.integral(- k * k * ENmono * dof(uMono), tf(uMono), pmlNMono, gauss); - - monodomain.integral(DWmono * grad(dof(uMono)), grad(tf(uMono)), pmlWMono, gauss); - monodomain.integral(- k * k * EWmono * dof(uMono), tf(uMono), pmlWMono, gauss); - - monodomain.integral(DSmono * DEmono * grad(dof(uMono)), grad(tf(uMono)), pmlSEMono, gauss); - monodomain.integral(- k * k * ESmono * EEmono * dof(uMono), tf(uMono), pmlSEMono, gauss); - - monodomain.integral(DNmono * DEmono * grad(dof(uMono)), grad(tf(uMono)), pmlNEMono, gauss); - monodomain.integral(- k * k * ENmono * EEmono * dof(uMono), tf(uMono), pmlNEMono, gauss); - - monodomain.integral(DNmono * DWmono * grad(dof(uMono)), grad(tf(uMono)), pmlNWMono, gauss); - monodomain.integral(- k * k * ENmono * EWmono * dof(uMono), tf(uMono), pmlNWMono, gauss); - - monodomain.integral(DSmono * DWmono * grad(dof(uMono)), grad(tf(uMono)), pmlSWMono, gauss); - monodomain.integral(- k * k * ESmono * EWmono * dof(uMono), tf(uMono), pmlSWMono, gauss); - - monodomain.pre(); - monodomain.assemble(); - monodomain.solve(); - - save(uMono); - - for(unsigned int i = 0; i < nDom; ++i) { - save(uMono - u(i), omega(i), "eMono_" + std::to_string(i)); - } - - // L2 error - double num = 0., den = 0.; - for(unsigned int i = 0; i < nDom; ++i) { - std::complex< double > denLocal = integrate(pow(abs(uMono), 2), omega(i), gauss); - std::complex< double > numLocal = integrate(pow(abs(uMono - u(i)), 2), omega(i), gauss); - num += numLocal.real(); - den += denLocal.real(); - } - msg::info << "L2 mono domain error = " << std::sqrt(num / den) << msg::endl; - } - - for(unsigned int i = 0; i < nDom; ++i) { - save(+u(i), omega(i), "u_" + std::to_string(i)); - //save(u(i)); - save(solution - u(i), omega(i), "e_" + std::to_string(i)); - for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { - const unsigned int j = topology[i][jj]; - if(!sigmaS(i,j).isEmpty()) save(gS(i,j)); - if(!sigmaE(i,j).isEmpty()) save(gE(i,j)); - if(!sigmaN(i,j).isEmpty()) save(gN(i,j)); - if(!sigmaW(i,j).isEmpty()) save(gW(i,j)); - } - } - - return 0; -} - -gmshfem::function::TensorFunction< std::complex< double > > getPmlDParameter(const double k, const double pmlSize, const double pos, const bool onX) -{ - gmshfem::function::ScalarFunction< std::complex< double > > kx = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > ky = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > kz = 1.; - - std::complex< double > I(0., 1.); - if(onX) { - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(x< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaX = 1. / (pmlSize - distSigma) - 1./ pmlSize; - kx = 1. + I * sigmaX / k; - } - else { // onY - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(y< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaY = 1. / (pmlSize - distSigma) - 1./ pmlSize; - ky = 1. + I * sigmaY / k; - } - - return tensorDiag< std::complex< double > >(ky * kz / kx, kx * kz / ky, kx * ky / kz); -} - -gmshfem::function::ScalarFunction< std::complex< double > > getPmlEParameter(const double k, const double pmlSize, const double pos, const bool onX) -{ - gmshfem::function::ScalarFunction< std::complex< double > > kx = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > ky = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > kz = 1.; - - std::complex< double > I(0., 1.); - if(onX) { - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(x< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaX = 1. / (pmlSize - distSigma) - 1./ pmlSize; - kx = 1. + I * sigmaX / k; - } - else { // onY - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(y< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaY = 1. / (pmlSize - distSigma) - 1./ pmlSize; - ky = 1. + I * sigmaY / k; - } - - return kx * ky * kz; -} diff --git a/demos/helmholtz/scattering/testWithoutGInPML/main.cpp b/demos/helmholtz/scattering/testWithoutGInPML/main.cpp deleted file mode 100644 index 8a35635126531b1be69c9b5a12122c9bf0455616..0000000000000000000000000000000000000000 --- a/demos/helmholtz/scattering/testWithoutGInPML/main.cpp +++ /dev/null @@ -1,767 +0,0 @@ -#include "GmshDdm.h" - -#include "FormulationDDM.h" - -#include "Message.h" -#include "AnalyticalFunction.h" - -#include "gmsh.h" - -using namespace gmshddm; -using namespace gmshddm::common; -using namespace gmshddm::domain; -using namespace gmshddm::problem; -using namespace gmshddm::field; - -using namespace gmshfem; -using namespace gmshfem::analytics; -using namespace gmshfem::domain; -using namespace gmshfem::problem; -using namespace gmshfem::field; -using namespace gmshfem::function; -using namespace gmshfem::post; - -gmshfem::function::TensorFunction< std::complex< double > > getPmlDParameter(const double k, const double pmlSize, const double pos, const bool onX); -gmshfem::function::ScalarFunction< std::complex< double > > getPmlEParameter(const double k, const double pmlSize, const double pos, const bool onX); - -void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, const double circleSize, const double pmlSize, const double lc, const bool pml = true, const unsigned int circlePosX = 0, const unsigned int circlePosY = 0) -{ - std::vector< std::vector< int > > omega(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - omega[i].resize(nDomY); - } - - // ****************** - // borders - // ****************** - - // ****(2)**** - // * * - // (3) (1) - // * * - // ****(0)**** - std::vector< std::vector< std::vector< int > > > borders(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - borders[i].resize(nDomY); - for(unsigned int j = 0; j < nDomY; ++j) { - borders[i][j].resize(4); - } - } - - // ****************** - // pmls - // ****************** - - // 1 0 - // 0 ****(2)**** 1 - // * * - // (3) (1) - // * * - // 1 ****(0)**** 0 - // 0 1 - - std::vector< std::vector< std::vector< int > > > pmls(nDomX); - std::vector< std::vector< std::vector< int > > > pmlsInf(nDomX); - std::vector< std::vector< std::vector< std::vector< int > > > > pmlsBnd(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - pmls[i].resize(nDomY); - pmlsInf[i].resize(nDomY); - pmlsBnd[i].resize(nDomY); - for(unsigned int j = 0; j < nDomY; ++j) { - pmls[i][j].resize(4); - pmlsInf[i][j].resize(4); - pmlsBnd[i][j].resize(4); - for(unsigned int k = 0; k < 4; ++k) { - pmlsBnd[i][j][k].resize(2); - } - } - } - - // ****************** - // pmlsSquare - // ****************** - - // (2)*******(1) - // * * - // * * - // * * - // (3)*******(0) - - std::vector< std::vector< std::vector< int > > > pmlsSquare(nDomX); - for(unsigned int i = 0; i < nDomX; ++i) { - pmlsSquare[i].resize(nDomY); - for(unsigned int j = 0; j < nDomY; ++j) { - pmlsSquare[i][j].resize(4); - } - } - - std::vector< int > left(nDomY); - int bottom = 0; - for(unsigned int i = 0; i < nDomX; ++i) { - for(unsigned int j = 0; j < nDomY; ++j) { - int p[4]; - p[0] = gmsh::model::geo::addPoint(i * xSize, j * ySize, 0., lc); - p[1] = gmsh::model::geo::addPoint((i + 1) * xSize, j * ySize, 0., lc); - p[2] = gmsh::model::geo::addPoint((i + 1) * xSize, (j + 1) * ySize, 0., lc); - p[3] = gmsh::model::geo::addPoint(i * xSize, (j + 1) * ySize, 0., lc); - - int l[4]; - borders[i][j][0] = l[0] = gmsh::model::geo::addLine(p[0], p[1]); - borders[i][j][1] = l[1] = gmsh::model::geo::addLine(p[1], p[2]); - borders[i][j][2] = l[2] = gmsh::model::geo::addLine(p[2], p[3]); - borders[i][j][3] = l[3] = gmsh::model::geo::addLine(p[3], p[0]); - - if(i == 0) { - left[j] = l[1]; - } - else { - borders[i][j][3] = left[j]; - left[j] = l[1]; - } - - if(j == 0) { - bottom = l[2]; - } - else { - borders[i][j][0] = bottom; - bottom = l[2]; - } - - int ll = gmsh::model::geo::addCurveLoop({l[0], l[1], l[2], l[3]}); - - int surface = 0; - if(i == circlePosX && j == circlePosY) { - // circle - int pC[5]; - pC[0] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2., circlePosY * j + ySize / 2., 0., lc); - pC[1] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2. + circleSize, circlePosY * j + ySize / 2., 0., lc); - pC[2] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2., circlePosY * j + ySize / 2. + circleSize, 0., lc); - pC[3] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2. - circleSize, circlePosY * j + ySize / 2., 0., lc); - pC[4] = gmsh::model::geo::addPoint(circlePosX * i + xSize / 2., circlePosY * j + ySize / 2. - circleSize, 0., lc); - - int lC[4]; - lC[0] = gmsh::model::geo::addCircleArc(pC[1], pC[0], pC[2]); - lC[1] = gmsh::model::geo::addCircleArc(pC[2], pC[0], pC[3]); - lC[2] = gmsh::model::geo::addCircleArc(pC[3], pC[0], pC[4]); - lC[3] = gmsh::model::geo::addCircleArc(pC[4], pC[0], pC[1]); - - // surface - int llC; - llC = gmsh::model::geo::addCurveLoop({lC[0], lC[1], lC[2], lC[3]}); - - surface = gmsh::model::geo::addPlaneSurface({ll, llC}); - - gmsh::model::addPhysicalGroup(1, {lC[0], lC[1], lC[2], lC[3]}, 1); - gmsh::model::setPhysicalName(1, 1, "gammaScat"); - } - else { - surface = gmsh::model::geo::addPlaneSurface({ll}); - } - omega[i][j] = surface; - - if(pml) { - // pmls - double extrudeX[4] = {0., pmlSize, 0., -pmlSize}; - double extrudeY[4] = {-pmlSize, 0., pmlSize, 0.}; - for(unsigned int k = 0; k < 4; ++k) { - std::vector< std::pair< int, int > > extrudeEntities; - gmsh::model::geo::extrude({std::make_pair(1, ((k == 3 && i != 0) || (k == 0 && j != 0) ? 1 : -1) * borders[i][j][k])}, extrudeX[k], extrudeY[k], 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true); - pmlsInf[i][j][k] = extrudeEntities[0].second; - pmls[i][j][k] = extrudeEntities[1].second; - if((k == 3 && i != 0) || (k == 0 && j != 0)) { - pmlsBnd[i][j][k][0] = extrudeEntities[2].second; - pmlsBnd[i][j][k][1] = extrudeEntities[3].second; - } - else { - pmlsBnd[i][j][k][0] = extrudeEntities[3].second; - pmlsBnd[i][j][k][1] = extrudeEntities[4].second; - } - } - - // pmlsSquare - double extrudeSquareX[4] = {pmlSize, pmlSize, -pmlSize, -pmlSize}; - int side[4] = {0, 2, 2, 0}; - int sideBnd[4] = {1, 0, 1, 0}; - for(unsigned int k = 0; k < 4; ++k) { - std::vector< std::pair< int, int > > extrudeEntities; - gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][side[k]][sideBnd[k]])}, extrudeSquareX[k], 0., 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true); - pmlsSquare[i][j][k] = extrudeEntities[1].second; - } - } - } - } - - gmsh::model::geo::synchronize(); - - // ****************** - // physicals - // ****************** - - for(unsigned int i = 0; i < nDomX; ++i) { - for(unsigned int j = 0; j < nDomY; ++j) { - gmsh::model::addPhysicalGroup(2, {omega[i][j]}, i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, i * nDomY + j + 1, "omega_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][0]}, 10000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 10000 + i * nDomY + j + 1, "sigmaS_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][1]}, 20000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 20000 + i * nDomY + j + 1, "sigmaE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][2]}, 30000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 30000 + i * nDomY + j + 1, "sigmaN_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {borders[i][j][3]}, 40000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 40000 + i * nDomY + j + 1, "sigmaW_" + std::to_string(i) + "_" + std::to_string(j)); - - if(pml) { - // PMLs - gmsh::model::addPhysicalGroup(2, {pmls[i][j][0]}, 10000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 10000 + i * nDomY + j + 1, "pmlS_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmls[i][j][1]}, 20000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 20000 + i * nDomY + j + 1, "pmlE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmls[i][j][2]}, 30000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 30000 + i * nDomY + j + 1, "pmlN_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmls[i][j][3]}, 40000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 40000 + i * nDomY + j + 1, "pmlW_" + std::to_string(i) + "_" + std::to_string(j)); - - // PMLs square - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][0]}, 100000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 100000 + i * nDomY + j + 1, "pmlSE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][1]}, 200000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 200000 + i * nDomY + j + 1, "pmlNE_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][2]}, 300000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 300000 + i * nDomY + j + 1, "pmlNW_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][3]}, 400000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(2, 400000 + i * nDomY + j + 1, "pmlSW_" + std::to_string(i) + "_" + std::to_string(j)); - - // Sigmas in PMLs - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3][1]}, 100000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 100000 + i * nDomY + j + 1, "sigmaS_W_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1][0]}, 150000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 150000 + i * nDomY + j + 1, "sigmaS_E_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0][1]}, 200000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 200000 + i * nDomY + j + 1, "sigmaE_S_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2][0]}, 250000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 250000 + i * nDomY + j + 1, "sigmaE_N_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1][1]}, 300000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 300000 + i * nDomY + j + 1, "sigmaN_E_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3][0]}, 350000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 350000 + i * nDomY + j + 1, "sigmaN_E_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2][1]}, 400000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 400000 + i * nDomY + j + 1, "sigmaW_N_" + std::to_string(i) + "_" + std::to_string(j)); - - gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0][0]}, 450000 + i * nDomY + j + 1); - gmsh::model::setPhysicalName(1, 450000 + i * nDomY + j + 1, "sigmaW_S_" + std::to_string(i) + "_" + std::to_string(j)); - - } - } - } - - gmsh::model::geo::synchronize(); - gmsh::model::mesh::generate(); -} - -int main(int argc, char **argv) -{ - GmshDdm gmshDdm(argc, argv); - - double pi = 3.14159265359; - - int nDomX = 2, nDomY = 2; - gmshDdm.userDefinedParameter(nDomX, "nDomX"); - gmshDdm.userDefinedParameter(nDomY, "nDomY"); - double xSize = 1., ySize = 1.; - gmshDdm.userDefinedParameter(xSize, "xSize"); - gmshDdm.userDefinedParameter(ySize, "ySize"); - double circleSize = 0.1; - gmshDdm.userDefinedParameter(circleSize, "circleSize"); - double pmlSize = 0.1; - gmshDdm.userDefinedParameter(pmlSize, "pmlSize"); - double k = 10. * pi; - gmshDdm.userDefinedParameter(k, "k"); - double lc = 0.02; - gmshDdm.userDefinedParameter(lc, "lc"); - std::string gauss = "Gauss3"; - gmshDdm.userDefinedParameter(gauss, "gauss"); - int circlePosX = 0, circlePosY = 0; - gmshDdm.userDefinedParameter(circlePosX, "circlePosX"); - gmshDdm.userDefinedParameter(circlePosY, "circlePosY"); - bool errorMono = false; - gmshDdm.userDefinedParameter(errorMono, "errorMono"); - - // Build geometry and mesh using gmsh API - checkerboard(nDomX, nDomY, xSize, ySize, circleSize, pmlSize, lc, true, circlePosX, circlePosY); - - unsigned int nDom = nDomX * nDomY; - - // | | - // NW | N | NE - // ----**************---- - // * * - // * * - // W * * E - // * * - // * * - // ----**************---- - // SW | S | SE - // | | - - // Define domain - Subdomain omega(nDom); - Subdomain gammaDir(nDom); - - Subdomain pmlS(nDom); - Subdomain pmlE(nDom); - Subdomain pmlN(nDom); - Subdomain pmlW(nDom); - - Subdomain pmlSE(nDom); - Subdomain pmlNE(nDom); - Subdomain pmlNW(nDom); - Subdomain pmlSW(nDom); - - Interface sigmaS(nDom); - Interface sigmaE(nDom); - Interface sigmaN(nDom); - Interface sigmaW(nDom); - - Interface sigmaS_W(nDom); - Interface sigmaS_E(nDom); - Interface sigmaE_S(nDom); - Interface sigmaE_N(nDom); - Interface sigmaN_E(nDom); - Interface sigmaN_W(nDom); - Interface sigmaW_N(nDom); - Interface sigmaW_S(nDom); - - Interface pmlS_W(nDom); - Interface pmlS_E(nDom); - Interface pmlE_S(nDom); - Interface pmlE_N(nDom); - Interface pmlN_E(nDom); - Interface pmlN_W(nDom); - Interface pmlW_N(nDom); - Interface pmlW_S(nDom); - - // Define topology - std::vector< std::vector< unsigned int > > topology(nDom); - for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) { - for(unsigned int j = 0; j < static_cast< unsigned int >(nDomY); ++j) { - unsigned int index = i * nDomY + j; - - omega(index) = Domain(2, index + 1); - - pmlS(index) = Domain(2, 10000 + index + 1); - pmlE(index) = Domain(2, 20000 + index + 1); - pmlN(index) = Domain(2, 30000 + index + 1); - pmlW(index) = Domain(2, 40000 + index + 1); - - if(i == static_cast< unsigned int >(nDomX) - 1 || j == 0) pmlSE(index) = Domain(2, 100000 + index + 1); - if(i == static_cast< unsigned int >(nDomX) - 1 || j == static_cast< unsigned int >(nDomY) - 1) pmlNE(index) = Domain(2, 200000 + index + 1); - if(i == 0 || j == static_cast< unsigned int >(nDomY) - 1) pmlNW(index) = Domain(2, 300000 + index + 1); - if(i == 0 || j == 0) pmlSW(index) = Domain(2, 400000 + index + 1); - - if(i == circlePosX && j == circlePosY) { - gammaDir(index) = Domain(1, 1); - } - - if(i != 0) { - sigmaW(index, (i-1) * nDomY + j) = Domain(1, 40000 + i * nDomY + j + 1); - - if(j == static_cast< unsigned int >(nDomY) - 1) sigmaW_N(index, (i-1) * nDomY + j) = Domain(1, 400000 + i * nDomY + j + 1); - if(j == 0) sigmaW_S(index, (i-1) * nDomY + j) = Domain(1, 450000 + i * nDomY + j + 1); - - if(j != static_cast< unsigned int >(nDomY) - 1) pmlW_N(index, (i-1) * nDomY + j) = Domain(1, 350000 + i * nDomY + j + 1); - if(j != 0) pmlW_S(index, (i-1) * nDomY + j) = Domain(1, 100000 + i * nDomY + j + 1); - - topology[index].push_back((i-1) * nDomY + j); - } - if(i != static_cast< unsigned int >(nDomX) - 1) { - sigmaE(index, (i+1) * nDomY + j) = Domain(1, 20000 + i * nDomY + j + 1); - - if(j == 0) sigmaE_S(index, (i+1) * nDomY + j) = Domain(1, 200000 + i * nDomY + j + 1); - if(j == static_cast< unsigned int >(nDomY) - 1) sigmaE_N(index, (i+1) * nDomY + j) = Domain(1, 250000 + i * nDomY + j + 1); - - if(j != 0) pmlE_S(index, (i+1) * nDomY + j) = Domain(1, 150000 + i * nDomY + j + 1); - if(j != static_cast< unsigned int >(nDomY) - 1) pmlE_N(index, (i+1) * nDomY + j) = Domain(1, 300000 + i * nDomY + j + 1); - - topology[index].push_back((i+1) * nDomY + j); - } - - if(j != 0) { - sigmaS(index, i * nDomY + (j-1)) = Domain(1, 10000 + i * nDomY + j + 1); - - if(i == 0) sigmaS_W(index, i * nDomY + (j-1)) = Domain(1, 100000 + i * nDomY + j + 1); - if(i == static_cast< unsigned int >(nDomX) - 1) sigmaS_E(index, i * nDomY + (j-1)) = Domain(1, 150000 + i * nDomY + j + 1); - - if(i != 0) pmlS_W(index, i * nDomY + (j-1)) = Domain(1, 450000 + i * nDomY + j + 1); - if(i != static_cast< unsigned int >(nDomX) - 1) pmlS_E(index, i * nDomY + (j-1)) = Domain(1, 200000 + i * nDomY + j + 1); - - topology[index].push_back(i * nDomY + (j-1)); - } - if(j != static_cast< unsigned int >(nDomY) - 1) { - sigmaN(index, i * nDomY + (j+1)) = Domain(1, 30000 + i * nDomY + j + 1); - - if(i == static_cast< unsigned int >(nDomX) - 1) sigmaN_E(index, i * nDomY + (j+1)) = Domain(1, 300000 + i * nDomY + j + 1); - if(i == 0) sigmaN_W(index, i * nDomY + (j+1)) = Domain(1, 350000 + i * nDomY + j + 1); - - if(i != static_cast< unsigned int >(nDomX) - 1) pmlN_E(index, i * nDomY + (j+1)) = Domain(1, 250000 + i * nDomY + j + 1); - if(i != 0) pmlN_W(index, i * nDomY + (j+1)) = Domain(1, 400000 + i * nDomY + j + 1); - - topology[index].push_back(i * nDomY + (j+1)); - } - } - } - - // Create resolution - Resolution< std::complex< double > > resolution("DDM"); - - // Create DDM formulation - FormulationDDM< std::complex< double > > formulation("Helmholtz", resolution, topology); - - // Create fields - SubdomainField< std::complex< double >, Form::Form0 > u(resolution, "u", omega | gammaDir | pmlS | pmlE | pmlN | pmlW | pmlSE | pmlNE | pmlNW | pmlSW, FunctionSpaceType::Lagrange); - SubdomainField< std::complex< double >, Form::Form0 > lambda(resolution, "lambda", gammaDir, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > gS(resolution, "gS", sigmaS | sigmaS_W | sigmaS_E, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > gE(resolution, "gE", sigmaE | sigmaE_S | sigmaE_N, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > gN(resolution, "gN", sigmaN | sigmaN_E | sigmaN_W, FunctionSpaceType::Lagrange); - InterfaceField< std::complex< double >, Form::Form0 > gW(resolution, "gW", sigmaW | sigmaW_N | sigmaW_S, FunctionSpaceType::Lagrange); - - std::complex< double > im(0., 1.); - AnalyticalFunction< helmholtz2D::SoftScatteringByACylinder< std::complex< double > > > solution(k, circleSize, xSize / 2., ySize / 2., 2 * k); - - // Tell to the formulation that g is field that have to be exchanged between subdomains - formulation.addInterfaceField(gS); - formulation.addInterfaceField(gE); - formulation.addInterfaceField(gN); - formulation.addInterfaceField(gW); - - // Add terms to the formulation - for(unsigned int i = 0; i < nDom; ++i) { - // Compute PML parameters - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaS = abs(y< std::complex< double > >() - (y< std::complex< double > >(omega(i)) - ySize / 2.)); - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaE = abs(x< std::complex< double > >() - (x< std::complex< double > >(omega(i)) + xSize / 2.)); - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaN = abs(y< std::complex< double > >() - (y< std::complex< double > >(omega(i)) + ySize / 2.)); - gmshfem::function::ScalarFunction< std::complex< double > > distSigmaW = abs(x< std::complex< double > >() - (x< std::complex< double > >(omega(i)) - xSize / 2.)); - - gmshfem::function::ScalarFunction< std::complex< double > > sigS = 1. / (pmlSize - distSigmaS) - 1./ pmlSize; - gmshfem::function::ScalarFunction< std::complex< double > > sigE = 1. / (pmlSize - distSigmaE) - 1./ pmlSize; - gmshfem::function::ScalarFunction< std::complex< double > > sigN = 1. / (pmlSize - distSigmaN) - 1./ pmlSize; - gmshfem::function::ScalarFunction< std::complex< double > > sigW = 1. / (pmlSize - distSigmaW) - 1./ pmlSize; - - gmshfem::function::ScalarFunction< std::complex< double > > kyS = 1. + im * sigS / k; - gmshfem::function::ScalarFunction< std::complex< double > > kxE = 1. + im * sigE / k; - gmshfem::function::ScalarFunction< std::complex< double > > kyN = 1. + im * sigN / k; - gmshfem::function::ScalarFunction< std::complex< double > > kxW = 1. + im * sigW / k; - - gmshfem::function::TensorFunction< std::complex< double > > DS = tensorDiag< std::complex< double > >(kyS, 1. / kyS, kyS); - gmshfem::function::ScalarFunction< std::complex< double > > ES = kyS; - - gmshfem::function::TensorFunction< std::complex< double > > DE = tensorDiag< std::complex< double > >(1. / kxE, kxE, kxE); - gmshfem::function::ScalarFunction< std::complex< double > > EE = kxE; - - gmshfem::function::TensorFunction< std::complex< double > > DN = tensorDiag< std::complex< double > >(kyN, 1. / kyN, kyN); - gmshfem::function::ScalarFunction< std::complex< double > > EN = kyN; - - gmshfem::function::TensorFunction< std::complex< double > > DW = tensorDiag< std::complex< double > >(1. / kxW, kxW, kxW); - gmshfem::function::ScalarFunction< std::complex< double > > EW = kxW; - - // VOLUME TERMS - - formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(- k * k * dof(u(i)), tf(u(i)), omega(i), gauss); - - // Physical source - formulation(i).integral(dof(lambda(i)), tf(u(i)), gammaDir(i), gauss); - formulation(i).integral(dof(u(i)), tf(lambda(i)), gammaDir(i), gauss); - formulation(i).integral(formulation.physicalSource(-solution), tf(lambda(i)), gammaDir(i), gauss); - - // S - formulation(i).integral(DS * grad(dof(u(i))), grad(tf(u(i))), pmlS(i), gauss); - formulation(i).integral(- ES * k * k * dof(u(i)), tf(u(i)), pmlS(i), gauss); - - // E - formulation(i).integral(DE * grad(dof(u(i))), grad(tf(u(i))), pmlE(i), gauss); - formulation(i).integral(- EE * k * k * dof(u(i)), tf(u(i)), pmlE(i), gauss); - - // N - formulation(i).integral(DN * grad(dof(u(i))), grad(tf(u(i))), pmlN(i), gauss); - formulation(i).integral(- EN * k * k * dof(u(i)), tf(u(i)), pmlN(i), gauss); - - // W - formulation(i).integral(DW * grad(dof(u(i))), grad(tf(u(i))), pmlW(i), gauss); - formulation(i).integral(- EW * k * k * dof(u(i)), tf(u(i)), pmlW(i), gauss); - - // SE - formulation(i).integral(DS * DE * grad(dof(u(i))), grad(tf(u(i))), pmlSE(i), gauss); - formulation(i).integral(- ES * EE * k * k * dof(u(i)), tf(u(i)), pmlSE(i), gauss); - - // NE - formulation(i).integral(DN * DE * grad(dof(u(i))), grad(tf(u(i))), pmlNE(i), gauss); - formulation(i).integral(- EN * EE * k * k * dof(u(i)), tf(u(i)), pmlNE(i), gauss); - - // NW - formulation(i).integral(DN * DW * grad(dof(u(i))), grad(tf(u(i))), pmlNW(i), gauss); - formulation(i).integral(- EN * EW * k * k * dof(u(i)), tf(u(i)), pmlNW(i), gauss); - - // SW - formulation(i).integral(DS * DW * grad(dof(u(i))), grad(tf(u(i))), pmlSW(i), gauss); - formulation(i).integral(- ES * EW * k * k * dof(u(i)), tf(u(i)), pmlSW(i), gauss); - - // Artificial sources - for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { - const unsigned int j = topology[i][jj]; - if(!sigmaS(i,j).isEmpty()) formulation(i).integral(formulation.artificialSource( -gN(j,i) ), tf(u(i)), sigmaS(i,j) | sigmaS_W(i,j) | sigmaS_E(i,j), gauss); - if(!sigmaE(i,j).isEmpty()) formulation(i).integral(formulation.artificialSource( -gW(j,i) ), tf(u(i)), sigmaE(i,j) | sigmaE_S(i,j) | sigmaE_N(i,j), gauss); - if(!sigmaN(i,j).isEmpty()) formulation(i).integral(formulation.artificialSource( -gS(j,i) ), tf(u(i)), sigmaN(i,j) | sigmaN_E(i,j) | sigmaN_W(i,j), gauss); - if(!sigmaW(i,j).isEmpty()) formulation(i).integral(formulation.artificialSource( -gE(j,i) ), tf(u(i)), sigmaW(i,j) | sigmaW_N(i,j) | sigmaW_S(i,j), gauss); - } - - // INTERFACE TERMS - for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { - const unsigned int j = topology[i][jj]; - // S - if(!sigmaS(i,j).isEmpty()) { - formulation(i,j).integral(dof(gS(i,j)), tf(gS(i,j)), sigmaS(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gN(j,i) ), tf(gS(i,j)), sigmaS(i,j), gauss); - - formulation(i,j).integral(dof(gS(i,j)), tf(gS(i,j)), sigmaS_W(i,j) | sigmaS_E(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gN(j,i) ), tf(gS(i,j)), sigmaS_W(i,j) | sigmaS_E(i,j), gauss); - - if(!pmlS_W(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gE(i - nDomY,i) ), tf(gS(i,j)), pmlS_W(i,j), gauss); - if(!pmlS_E(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gW(i + nDomY,i) ), tf(gS(i,j)), pmlS_E(i,j), gauss); - - formulation(i,j).integral(- 2. * DS * grad(u(i)), grad(tf(gS(i,j))), pmlS(i), gauss); - formulation(i,j).integral(2. * ES * k * k * u(i), tf(gS(i,j)), pmlS(i), gauss); - - formulation(i,j).integral(- 2. * DS * DW * grad(u(i)), grad(tf(gS(i,j))), pmlSW(i), gauss); - formulation(i,j).integral(2. * ES * EW * k * k * u(i), tf(gS(i,j)), pmlSW(i), gauss); - - formulation(i,j).integral(- 2. * DS * DE * grad(u(i)), grad(tf(gS(i,j))), pmlSE(i), gauss); - formulation(i,j).integral(2. * ES * EE * k * k * u(i), tf(gS(i,j)), pmlSE(i), gauss); - } - - // E - if(!sigmaE(i,j).isEmpty()) { - formulation(i,j).integral(dof(gE(i,j)), tf(gE(i,j)), sigmaE(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gW(j,i) ), tf(gE(i,j)), sigmaE(i,j), gauss); - - formulation(i,j).integral(dof(gE(i,j)), tf(gE(i,j)), sigmaE_S(i,j) | sigmaE_N(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gW(j,i) ), tf(gE(i,j)), sigmaE_S(i,j) | sigmaE_N(i,j), gauss); - - if(!pmlE_S(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gN(i - 1,i) ), tf(gE(i,j)), pmlE_S(i,j), gauss); - if(!pmlE_N(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gS(i + 1,i) ), tf(gE(i,j)), pmlE_N(i,j), gauss); - - formulation(i,j).integral(- 2. * DE * grad(u(i)), grad(tf(gE(i,j))), pmlE(i), gauss); - formulation(i,j).integral(2. * EE * k * k * u(i), tf(gE(i,j)), pmlE(i), gauss); - - formulation(i,j).integral(- 2. * DS * DE * grad(u(i)), grad(tf(gE(i,j))), pmlSE(i), gauss); - formulation(i,j).integral(2. * ES * EE * k * k * u(i), tf(gE(i,j)), pmlSE(i), gauss); - - formulation(i,j).integral(- 2. * DN * DE * grad(u(i)), grad(tf(gE(i,j))), pmlNE(i), gauss); - formulation(i,j).integral(2. * EN * EE * k * k * u(i), tf(gE(i,j)), pmlNE(i), gauss); - } - - // N - if(!sigmaN(i,j).isEmpty()) { - formulation(i,j).integral(dof(gN(i,j)), tf(gN(i,j)), sigmaN(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gS(j,i) ), tf(gN(i,j)), sigmaN(i,j), gauss); - - formulation(i,j).integral(dof(gN(i,j)), tf(gN(i,j)), sigmaN_E(i,j) | sigmaN_W(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gS(j,i) ), tf(gN(i,j)), sigmaN_E(i,j) | sigmaN_W(i,j), gauss); - - if(!pmlN_E(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gW(i + nDomY,i) ), tf(gN(i,j)), pmlN_E(i,j), gauss); - if(!pmlN_W(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gE(i - nDomY,i) ), tf(gN(i,j)), pmlN_W(i,j), gauss); - - formulation(i,j).integral(- 2. * DN * grad(u(i)), grad(tf(gN(i,j))), pmlN(i), gauss); - formulation(i,j).integral(2. * EN * k * k * u(i), tf(gN(i,j)), pmlN(i), gauss); - - formulation(i,j).integral(- 2. * DN * DE * grad(u(i)), grad(tf(gN(i,j))), pmlNE(i), gauss); - formulation(i,j).integral(2. * EN * EE * k * k * u(i), tf(gN(i,j)), pmlNE(i), gauss); - - formulation(i,j).integral(- 2. * DN * DW * grad(u(i)), grad(tf(gN(i,j))), pmlNW(i), gauss); - formulation(i,j).integral(2. * EN * EW * k * k * u(i), tf(gN(i,j)), pmlNW(i), gauss); - } - - // W - if(!sigmaW(i,j).isEmpty()) { - formulation(i,j).integral(dof(gW(i,j)), tf(gW(i,j)), sigmaW(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gE(j,i) ), tf(gW(i,j)), sigmaW(i,j), gauss); - - formulation(i,j).integral(dof(gW(i,j)), tf(gW(i,j)), sigmaW_N(i,j) | sigmaW_S(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gE(j,i) ), tf(gW(i,j)), sigmaW_N(i,j) | sigmaW_S(i,j), gauss); - - if(!pmlW_N(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gS(i + 1,i) ), tf(gW(i,j)), pmlW_N(i,j), gauss); - if(!pmlW_S(i,j).isEmpty()) formulation(i,j).integral(formulation.artificialSource( 2. * gN(i - 1,i) ), tf(gW(i,j)), pmlW_S(i,j), gauss); - - formulation(i,j).integral(- 2. * DW * grad(u(i)), grad(tf(gW(i,j))), pmlW(i), gauss); - formulation(i,j).integral(2. * EW * k * k * u(i), tf(gW(i,j)), pmlW(i), gauss); - - formulation(i,j).integral(- 2. * DN * DW * grad(u(i)), grad(tf(gW(i,j))), pmlNW(i), gauss); - formulation(i,j).integral(2. * EN * EW * k * k * u(i), tf(gW(i,j)), pmlNW(i), gauss); - - formulation(i,j).integral(- 2. * DS * DW * grad(u(i)), grad(tf(gW(i,j))), pmlSW(i), gauss); - formulation(i,j).integral(2. * ES * EW * k * k * u(i), tf(gW(i,j)), pmlSW(i), gauss); - } - } - } - - // Solve the DDM formulation - formulation.pre(); - formulation.solve("gmres"); - - if(errorMono) { - gmshfem::domain::Domain omegaMono, pmlSMono, pmlEMono, pmlNMono, pmlWMono, pmlSEMono, pmlNEMono, pmlNWMono, pmlSWMono; - for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { - omegaMono |= omega(i); - } - - for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) { - pmlSMono |= pmlS(i * nDomY + 0); - pmlNMono |= pmlN(i * nDomY + nDomY - 1); - } - for(unsigned int j = 0; j < static_cast< unsigned int >(nDomY); ++j) { - pmlEMono |= pmlE((nDomX - 1) * nDomY + j); - pmlWMono |= pmlW(0 * nDomY + j); - } - pmlSEMono = pmlSE((nDomX - 1) * nDomY + 0); - pmlNEMono = pmlNE((nDomX - 1) * nDomY + nDomY - 1); - pmlNWMono = pmlNW(0 * nDomY + nDomY - 1); - pmlSWMono = pmlSW(0); - - gmshfem::function::TensorFunction< std::complex< double > > DSmono = getPmlDParameter(k, pmlSize, 0., false); - gmshfem::function::TensorFunction< std::complex< double > > DEmono = getPmlDParameter(k, pmlSize, nDomX * xSize, true); - gmshfem::function::TensorFunction< std::complex< double > > DNmono = getPmlDParameter(k, pmlSize, nDomY * ySize, false); - gmshfem::function::TensorFunction< std::complex< double > > DWmono = getPmlDParameter(k, pmlSize, 0., true); - - gmshfem::function::ScalarFunction< std::complex< double > > ESmono = getPmlEParameter(k, pmlSize, 0., false); - gmshfem::function::ScalarFunction< std::complex< double > > EEmono = getPmlEParameter(k, pmlSize, nDomX * xSize, true); - gmshfem::function::ScalarFunction< std::complex< double > > ENmono = getPmlEParameter(k, pmlSize, nDomY * ySize, false); - gmshfem::function::ScalarFunction< std::complex< double > > EWmono = getPmlEParameter(k, pmlSize, 0., true); - - gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > uMono(resolution, "uMono", gammaDir(circlePosX * nDomY + circlePosY) | omegaMono | pmlSMono | pmlEMono | pmlNMono | pmlWMono | pmlSEMono | pmlNEMono | pmlNWMono | pmlSWMono, FunctionSpaceType::Lagrange); - gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambdaMono(resolution, "lambdaMono", gammaDir(circlePosX * nDomY + circlePosY), FunctionSpaceType::Lagrange); - - Formulation< std::complex< double > > monodomain("HelmholtzMono", resolution); - - monodomain.integral(grad(dof(uMono)), grad(tf(uMono)), omegaMono, gauss); - monodomain.integral(- k * k * dof(uMono), tf(uMono), omegaMono, gauss); - - monodomain.integral(dof(lambdaMono), tf(uMono), gammaDir(circlePosX * nDomY + circlePosY), gauss); - monodomain.integral(dof(uMono), tf(lambdaMono), gammaDir(circlePosX * nDomY + circlePosY), gauss); - monodomain.integral(- solution, tf(lambdaMono), gammaDir(circlePosX * nDomY + circlePosY), gauss); - - monodomain.integral(DSmono * grad(dof(uMono)), grad(tf(uMono)), pmlSMono, gauss); - monodomain.integral(- k * k * ESmono * dof(uMono), tf(uMono), pmlSMono, gauss); - - monodomain.integral(DEmono * grad(dof(uMono)), grad(tf(uMono)), pmlEMono, gauss); - monodomain.integral(- k * k * EEmono * dof(uMono), tf(uMono), pmlEMono, gauss); - - monodomain.integral(DNmono * grad(dof(uMono)), grad(tf(uMono)), pmlNMono, gauss); - monodomain.integral(- k * k * ENmono * dof(uMono), tf(uMono), pmlNMono, gauss); - - monodomain.integral(DWmono * grad(dof(uMono)), grad(tf(uMono)), pmlWMono, gauss); - monodomain.integral(- k * k * EWmono * dof(uMono), tf(uMono), pmlWMono, gauss); - - monodomain.integral(DSmono * DEmono * grad(dof(uMono)), grad(tf(uMono)), pmlSEMono, gauss); - monodomain.integral(- k * k * ESmono * EEmono * dof(uMono), tf(uMono), pmlSEMono, gauss); - - monodomain.integral(DNmono * DEmono * grad(dof(uMono)), grad(tf(uMono)), pmlNEMono, gauss); - monodomain.integral(- k * k * ENmono * EEmono * dof(uMono), tf(uMono), pmlNEMono, gauss); - - monodomain.integral(DNmono * DWmono * grad(dof(uMono)), grad(tf(uMono)), pmlNWMono, gauss); - monodomain.integral(- k * k * ENmono * EWmono * dof(uMono), tf(uMono), pmlNWMono, gauss); - - monodomain.integral(DSmono * DWmono * grad(dof(uMono)), grad(tf(uMono)), pmlSWMono, gauss); - monodomain.integral(- k * k * ESmono * EWmono * dof(uMono), tf(uMono), pmlSWMono, gauss); - - monodomain.pre(); - monodomain.assemble(); - monodomain.solve(); - - save(uMono); - - for(unsigned int i = 0; i < nDom; ++i) { - save(uMono - u(i), omega(i), "eMono_" + std::to_string(i)); - } - - // L2 error - double num = 0., den = 0.; - for(unsigned int i = 0; i < nDom; ++i) { - std::complex< double > denLocal = integrate(pow(abs(uMono), 2), omega(i), gauss); - std::complex< double > numLocal = integrate(pow(abs(uMono - u(i)), 2), omega(i), gauss); - num += numLocal.real(); - den += denLocal.real(); - } - msg::info << "L2 mono domain error = " << std::sqrt(num / den) << msg::endl; - } - - for(unsigned int i = 0; i < nDom; ++i) { - save(+u(i), omega(i), "u_" + std::to_string(i)); - save(solution - u(i), omega(i), "e_" + std::to_string(i)); - for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { - const unsigned int j = topology[i][jj]; - if(!sigmaS(i,j).isEmpty()) save(gS(i,j)); - if(!sigmaE(i,j).isEmpty()) save(gE(i,j)); - if(!sigmaN(i,j).isEmpty()) save(gN(i,j)); - if(!sigmaW(i,j).isEmpty()) save(gW(i,j)); - } - } - - return 0; -} - -gmshfem::function::TensorFunction< std::complex< double > > getPmlDParameter(const double k, const double pmlSize, const double pos, const bool onX) -{ - gmshfem::function::ScalarFunction< std::complex< double > > kx = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > ky = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > kz = 1.; - - std::complex< double > I(0., 1.); - if(onX) { - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(x< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaX = 1. / (pmlSize - distSigma) - 1./ pmlSize; - kx = 1. + I * sigmaX / k; - } - else { // onY - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(y< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaY = 1. / (pmlSize - distSigma) - 1./ pmlSize; - ky = 1. + I * sigmaY / k; - } - - return tensorDiag< std::complex< double > >(ky * kz / kx, kx * kz / ky, kx * ky / kz); -} - -gmshfem::function::ScalarFunction< std::complex< double > > getPmlEParameter(const double k, const double pmlSize, const double pos, const bool onX) -{ - gmshfem::function::ScalarFunction< std::complex< double > > kx = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > ky = 1.; - gmshfem::function::ScalarFunction< std::complex< double > > kz = 1.; - - std::complex< double > I(0., 1.); - if(onX) { - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(x< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaX = 1. / (pmlSize - distSigma) - 1./ pmlSize; - kx = 1. + I * sigmaX / k; - } - else { // onY - gmshfem::function::ScalarFunction< std::complex< double > > distSigma = abs(y< std::complex< double > >() - pos); - gmshfem::function::ScalarFunction< std::complex< double > > sigmaY = 1. / (pmlSize - distSigma) - 1./ pmlSize; - ky = 1. + I * sigmaY / k; - } - - return kx * ky * kz; -} diff --git a/demos/helmholtzflow/freefield/runTests.sh b/demos/helmholtzflow/freefield/runTests.sh deleted file mode 100644 index 8b4741f05eb05a01dd98340e9114503a4e962603..0000000000000000000000000000000000000000 --- a/demos/helmholtzflow/freefield/runTests.sh +++ /dev/null @@ -1,19 +0,0 @@ -cd build -cmake .. && make - -echo "****************************" -echo " zeroth order TC " -echo "****************************" -./demo -transmission T0 -alpha 0 -M 0.8 -echo "****************************" -echo " zeroth order TC with rotation " -echo "****************************" -./demo -transmission T0 -alpha -1 -M 0.8 -echo "****************************" -echo " second order TC with rotation " -echo "****************************" -./demo -transmission T2 -alpha -1 -M 0.8 -echo "****************************" -echo " Pade TC with rotation " -echo "****************************" -./demo -transmission Pade -padeOrder 4 -alpha -1 -M 0.8 diff --git a/demos/navier/scattering/CMakeLists.txt b/demos/navier/scattering/CMakeLists.txt deleted file mode 100755 index 4d2e4b0c91d2a7e0da5a224e581c56eca24b2f71..0000000000000000000000000000000000000000 --- a/demos/navier/scattering/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -cmake_minimum_required(VERSION 3.0 FATAL_ERROR) - -project(demo CXX) - -include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake) - -add_executable(demo main.cpp ${EXTRA_INCS}) -target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/examples/.clang-format b/examples/.clang-format new file mode 100644 index 0000000000000000000000000000000000000000..d902ab979dca0361bb1b063642225163bda4eb7d --- /dev/null +++ b/examples/.clang-format @@ -0,0 +1,112 @@ +--- +Language: Cpp +AccessModifierOffset: -1 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +AlignEscapedNewlines: Right +AlignOperands: true +AlignTrailingComments: false +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: false +AllowShortCaseLabelsOnASingleLine: true +AllowShortFunctionsOnASingleLine: Empty +AllowShortIfStatementsOnASingleLine: true +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: true +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: true + AfterControlStatement: false + AfterEnum: false + AfterFunction: true + AfterNamespace: true + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + AfterExternBlock: false + BeforeCatch: true + BeforeElse: true + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Custom +BreakBeforeInheritanceComma: false +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: AfterColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 0 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 2 +ContinuationIndentWidth: 2 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IncludeBlocks: Regroup +IncludeCategories: + - Regex: '^"(?!gmsh\.h"|petsc\.h"|Eigen")' + Priority: 2 + - Regex: '"gmsh.h"|"petsc.h"|"Eigen"' + Priority: 3 + - Regex: '^<' + Priority: 4 + - Regex: '.*' + Priority: 1 +IncludeIsMainRegex: '$' +IndentCaseLabels: false +IndentPPDirectives: None +IndentWidth: 2 +IndentWrappedFunctionNames: false +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 2 +NamespaceIndentation: All +ObjCBlockIndentWidth: 2 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 60 +PointerAlignment: Right +ReflowComments: false +SortIncludes: true +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterTemplateKeyword: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: Never +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: true +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Cpp11 +TabWidth: 4 +UseTab: Never +... + diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..a3b6e79a56fd112f02ac9a3fa15103ccec678000 --- /dev/null +++ b/examples/CMakeLists.txt @@ -0,0 +1,14 @@ +set(EXAMPLES + helmholtz/checkerboard_bcgs + helmholtz/scattering/planeWaveDisk + helmholtzflow/freefield + helmholtzflow/waveguide + maxwell + navier/scattering +) + +foreach(EXAMPLE ${EXAMPLES}) + add_test(NAME "example:${EXAMPLE}:mkdir" COMMAND "${CMAKE_COMMAND}" -E make_directory build WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/${EXAMPLE}) + add_test(NAME "example:${EXAMPLE}:cmake" COMMAND "${CMAKE_COMMAND}" -DCMAKE_PREFIX_PATH=${CMAKE_INSTALL_PREFIX} .. WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/${EXAMPLE}/build) + add_test(NAME "example:${EXAMPLE}:build" COMMAND "${CMAKE_MAKE_PROGRAM}" WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/${EXAMPLE}/build) +endforeach() diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000000000000000000000000000000000000..82741b66daaabf64ca87e4b08fc3b59078bc0483 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,4 @@ +# Examples +> This directory contains the GmsmDDM research codes. + +This directory contains research codes based on GmshDDM. diff --git a/examples/README_template.md b/examples/README_template.md new file mode 100644 index 0000000000000000000000000000000000000000..a4bf0e06dc73e5da7d4111f01de47222af3005b6 --- /dev/null +++ b/examples/README_template.md @@ -0,0 +1,23 @@ +# Project name +> author 1, author 2, ... + +Dependency(ies): +* GmshDDM : [commit hash] +* GmshFEM : [commit hash] +* ... + +## Project description + +Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum. + +## Installation + +Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum. + +## Usage + +Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum. + +## References + +Papers, proceedings, ... diff --git a/demos/demos.cmake b/examples/example.cmake similarity index 97% rename from demos/demos.cmake rename to examples/example.cmake index 38ebee16abc036fb6360ea366bb3cc8d4fb54c25..3e2c7ef61683c3ee4d8e3652f5103ed81c1c80c3 100644 --- a/demos/demos.cmake +++ b/examples/example.cmake @@ -50,7 +50,8 @@ if(NOT GMSHFEM_INC) message(FATAL_ERROR "Could not find GmshFEM header 'gmshfem/GmshFem.h'") endif(NOT GMSHFEM_INC) list(APPEND EXTRA_LIBS ${GMSHFEM_LIB}) -list(APPEND EXTRA_INCS ${GMSHFEM_INC}/gmshfem/) +list(APPEND EXTRA_INCS ${GMSHFEM_INC}) +list(APPEND EXTRA_INCS ${GMSHFEM_INC}/gmshfem/contrib) # Find GmshDDM find_library(GMSHDDM_LIB gmshddm) diff --git a/demos/helmholtz/waveguide/CMakeLists.txt b/examples/helmholtz/checkerboard_bcgs/CMakeLists.txt similarity index 73% rename from demos/helmholtz/waveguide/CMakeLists.txt rename to examples/helmholtz/checkerboard_bcgs/CMakeLists.txt index 4d2e4b0c91d2a7e0da5a224e581c56eca24b2f71..c368ce59f63aed4e569c1b3e727df6dd0b1d6a45 100644 --- a/demos/helmholtz/waveguide/CMakeLists.txt +++ b/examples/helmholtz/checkerboard_bcgs/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.0 FATAL_ERROR) project(demo CXX) -include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake) +include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake) add_executable(demo main.cpp ${EXTRA_INCS}) target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/demos/helmholtz/checkerboard_bcgs/conv_test.txt b/examples/helmholtz/checkerboard_bcgs/conv_test.txt similarity index 100% rename from demos/helmholtz/checkerboard_bcgs/conv_test.txt rename to examples/helmholtz/checkerboard_bcgs/conv_test.txt diff --git a/demos/helmholtz/checkerboard_bcgs/main.cpp b/examples/helmholtz/checkerboard_bcgs/main.cpp similarity index 82% rename from demos/helmholtz/checkerboard_bcgs/main.cpp rename to examples/helmholtz/checkerboard_bcgs/main.cpp index 69213f7a905984ab1a6ff5f465eb588f6f527d5c..df7859e44cbccd77b4c8994433ec6d1650e27eba 100644 --- a/demos/helmholtz/checkerboard_bcgs/main.cpp +++ b/examples/helmholtz/checkerboard_bcgs/main.cpp @@ -1,10 +1,8 @@ -#include <gmshddm/GmshDdm.h> +#include <gmsh.h> #include <gmshddm/Formulation.h> - +#include <gmshddm/GmshDdm.h> #include <gmshfem/Message.h> -#include <gmsh.h> - using gmshfem::equation::dof; using gmshfem::equation::tf; using gmshfem::function::operator-; @@ -51,10 +49,10 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub // * * // (p0)*******(p1) int p[4]; - p[0] = gmsh::model::geo::addPoint(i * xSize, j * ySize, 0., lc); - p[1] = gmsh::model::geo::addPoint((i + 1) * xSize, j * ySize, 0., lc); + p[0] = gmsh::model::geo::addPoint(i * xSize, j * ySize, 0., lc); + p[1] = gmsh::model::geo::addPoint((i + 1) * xSize, j * ySize, 0., lc); p[2] = gmsh::model::geo::addPoint((i + 1) * xSize, (j + 1) * ySize, 0., lc); - p[3] = gmsh::model::geo::addPoint(i * xSize, (j + 1) * ySize, 0., lc); + p[3] = gmsh::model::geo::addPoint(i * xSize, (j + 1) * ySize, 0., lc); int l[4]; borders[i][j][0] = l[0] = gmsh::model::geo::addLine(p[1], p[2]); @@ -93,7 +91,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub gmsh::model::addPhysicalGroup(0, {pS}, 1); gmsh::model::setPhysicalName(0, 1, "source"); } - } } gmsh::model::geo::removeAllDuplicates(); @@ -109,16 +106,16 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub gmsh::model::addPhysicalGroup(2, {omega[i][j]}, i * nDomY + j + 1); gmsh::model::setPhysicalName(2, i * nDomY + j + 1, "omega_" + std::to_string(i) + "_" + std::to_string(j)); - gmsh::model::addPhysicalGroup(1, {borders[i][j][0]}, 10000 + i * nDomY + j + 1); + gmsh::model::addPhysicalGroup(1, {borders[i][j][0]}, 10000 + i * nDomY + j + 1); gmsh::model::setPhysicalName(1, 10000 + i * nDomY + j + 1, "sigmaE_" + std::to_string(i) + "_" + std::to_string(j)); - gmsh::model::addPhysicalGroup(1, {borders[i][j][1]}, 20000 + i * nDomY + j + 1); + gmsh::model::addPhysicalGroup(1, {borders[i][j][1]}, 20000 + i * nDomY + j + 1); gmsh::model::setPhysicalName(1, 20000 + i * nDomY + j + 1, "sigmaN_" + std::to_string(i) + "_" + std::to_string(j)); - gmsh::model::addPhysicalGroup(1, {borders[i][j][2]}, 30000 + i * nDomY + j + 1); + gmsh::model::addPhysicalGroup(1, {borders[i][j][2]}, 30000 + i * nDomY + j + 1); gmsh::model::setPhysicalName(1, 30000 + i * nDomY + j + 1, "sigmaW_" + std::to_string(i) + "_" + std::to_string(j)); - gmsh::model::addPhysicalGroup(1, {borders[i][j][3]}, 40000 + i * nDomY + j + 1); + gmsh::model::addPhysicalGroup(1, {borders[i][j][3]}, 40000 + i * nDomY + j + 1); gmsh::model::setPhysicalName(1, 40000 + i * nDomY + j + 1, "sigmaS_" + std::to_string(i) + "_" + std::to_string(j)); } } @@ -175,27 +172,27 @@ int main(int argc, char **argv) omega(index) = gmshfem::domain::Domain(2, index + 1); for(unsigned int f = 0; f < 4; ++f) { - sigma[f](index) = gmshfem::domain::Domain(1, (f+1) * 10000 + i * nDomY + j + 1); + sigma[f](index) = gmshfem::domain::Domain(1, (f + 1) * 10000 + i * nDomY + j + 1); } if(i != static_cast< unsigned int >(nDomX) - 1) { // E - sigmaInterface[0](index, (i+1) * nDomY + j) = sigma[0](index); - topology[index].push_back((i+1) * nDomY + j); + sigmaInterface[0](index, (i + 1) * nDomY + j) = sigma[0](index); + topology[index].push_back((i + 1) * nDomY + j); } if(j != static_cast< unsigned int >(nDomY) - 1) { // N - sigmaInterface[1](index, i * nDomY + (j+1)) = sigma[1](index); - topology[index].push_back(i * nDomY + (j+1)); + sigmaInterface[1](index, i * nDomY + (j + 1)) = sigma[1](index); + topology[index].push_back(i * nDomY + (j + 1)); } if(i != 0) { // W - sigmaInterface[2](index, (i-1) * nDomY + j) = sigma[2](index); - topology[index].push_back((i-1) * nDomY + j); + sigmaInterface[2](index, (i - 1) * nDomY + j) = sigma[2](index); + topology[index].push_back((i - 1) * nDomY + j); } if(j != 0) { // S - sigmaInterface[3](index, i * nDomY + (j-1)) = sigma[3](index); - topology[index].push_back(i * nDomY + (j-1)); + sigmaInterface[3](index, i * nDomY + (j - 1)) = sigma[3](index); + topology[index].push_back(i * nDomY + (j - 1)); } } } @@ -208,7 +205,7 @@ int main(int argc, char **argv) gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambda("lambda", gamma, gmshfem::field::FunctionSpaceTypeForm0::Lagrange); - const std::vector< std::string > orientation {"E", "N", "W", "S"}; + const std::vector< std::string > orientation{"E", "N", "W", "S"}; std::vector< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > > g; g.reserve(4); for(unsigned int f = 0; f < 4; ++f) { @@ -229,7 +226,7 @@ int main(int argc, char **argv) if(i == sourcePosX * nDomY + sourcePosY) { formulation(i).integral(dof(lambda), tf(u(i)), gamma, gauss); formulation(i).integral(dof(u(i)), tf(lambda), gamma, gauss); - formulation(i).integral(formulation.physicalSource(- 1.), tf(lambda), gamma, gauss); + formulation(i).integral(formulation.physicalSource(-1.), tf(lambda), gamma, gauss); } // boundary @@ -245,8 +242,8 @@ int main(int argc, char **argv) for(unsigned int f = 0; f < 4; ++f) { for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - if(!sigmaInterface[f](i,j).isEmpty()) { - formulation(i).integral(formulation.artificialSource( - g[(f+2)%4](j,i) ), tf(u(i)), sigmaInterface[f](i,j), gauss); + if(!sigmaInterface[f](i, j).isEmpty()) { + formulation(i).integral(formulation.artificialSource(-g[(f + 2) % 4](j, i)), tf(u(i)), sigmaInterface[f](i, j), gauss); } } } @@ -255,11 +252,10 @@ int main(int argc, char **argv) for(unsigned int f = 0; f < 4; ++f) { for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - if(!sigmaInterface[f](i,j).isEmpty()) { - formulation(i,j).integral(dof(g[f](i,j)), tf(g[f](i,j)), sigmaInterface[f](i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( g[(f+2)%4](j,i) ), tf(g[f](i,j)), sigmaInterface[f](i,j), gauss); - formulation(i,j).integral(2. * v[f](i), tf(g[f](i,j)), sigmaInterface[f](i,j), gauss); - + if(!sigmaInterface[f](i, j).isEmpty()) { + formulation(i, j).integral(dof(g[f](i, j)), tf(g[f](i, j)), sigmaInterface[f](i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(g[(f + 2) % 4](j, i)), tf(g[f](i, j)), sigmaInterface[f](i, j), gauss); + formulation(i, j).integral(2. * v[f](i), tf(g[f](i, j)), sigmaInterface[f](i, j), gauss); } } } diff --git a/demos/helmholtz/scattering/planeWaveDisk/CMakeLists.txt b/examples/helmholtz/scattering/planeWaveDisk/CMakeLists.txt similarity index 72% rename from demos/helmholtz/scattering/planeWaveDisk/CMakeLists.txt rename to examples/helmholtz/scattering/planeWaveDisk/CMakeLists.txt index f01d5714177b61745622f922015d7abd78c1cdcf..78dcc127355d6e26101a160e271312b86d9565ab 100644 --- a/demos/helmholtz/scattering/planeWaveDisk/CMakeLists.txt +++ b/examples/helmholtz/scattering/planeWaveDisk/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.0 FATAL_ERROR) project(demo CXX) -include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake) +include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake) add_executable(demo main.cpp ${EXTRA_INCS}) target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/demos/helmholtz/scattering/planeWaveDisk/main.cpp b/examples/helmholtz/scattering/planeWaveDisk/main.cpp similarity index 70% rename from demos/helmholtz/scattering/planeWaveDisk/main.cpp rename to examples/helmholtz/scattering/planeWaveDisk/main.cpp index c58ab82948217de38e5090b2145ddda19fa432c2..22f639c6e792c03b96693e727eeebcae6c58104d 100644 --- a/demos/helmholtz/scattering/planeWaveDisk/main.cpp +++ b/examples/helmholtz/scattering/planeWaveDisk/main.cpp @@ -1,8 +1,10 @@ -#include <gmshddm/GmshDdm.h> +#include "gmsh.h" + #include <gmshddm/Formulation.h> -#include <gmshfem/Message.h> +#include <gmshddm/GmshDdm.h> +#include <gmshddm/MPIInterface.h> #include <gmshfem/AnalyticalFunction.h> -#include "gmsh.h" +#include <gmshfem/Message.h> using namespace gmshddm; using namespace gmshddm::common; @@ -27,20 +29,20 @@ void circlesConcentric(const unsigned int nDom, const double r0, const double rI } gmsh::model::add("circlesConcentric"); - std::vector< int > tagsCircle(nDom+2); + std::vector< int > tagsCircle(nDom + 2); for(unsigned int i = 0; i <= nDom; ++i) { tagsCircle[i] = gmsh::model::occ::addCircle(0., 0., 0., std::sqrt(double(i * (rInf * rInf - r0 * r0)) / nDom + r0 * r0)); } - tagsCircle[nDom+1] = gmsh::model::occ::addCircle(0., 0., 0., std::sqrt(double( nDom * (rPml * rPml - r0 * r0)) / nDom + r0 * r0)); + tagsCircle[nDom + 1] = gmsh::model::occ::addCircle(0., 0., 0., std::sqrt(double(nDom * (rPml * rPml - r0 * r0)) / nDom + r0 * r0)); std::vector< int > tagsCurve(nDom + 2); - for(unsigned int i = 0; i <= (nDom+1); ++i) { + for(unsigned int i = 0; i <= (nDom + 1); ++i) { tagsCurve[i] = gmsh::model::occ::addCurveLoop({tagsCircle[i]}); } - - std::vector< int > tagsSurface(nDom+1); + + std::vector< int > tagsSurface(nDom + 1); for(unsigned int i = 0; i <= nDom; ++i) { - tagsSurface[i] = gmsh::model::occ::addPlaneSurface({tagsCurve[i], tagsCurve[i+1]}); + tagsSurface[i] = gmsh::model::occ::addPlaneSurface({tagsCurve[i], tagsCurve[i + 1]}); } gmsh::model::occ::synchronize(); @@ -50,24 +52,24 @@ void circlesConcentric(const unsigned int nDom, const double r0, const double rI gmsh::model::addPhysicalGroup(1, {tagsCircle[nDom]}, 92); gmsh::model::setPhysicalName(1, 92, "gammaInf"); - - gmsh::model::addPhysicalGroup(1, {tagsCircle[nDom+1]}, 93); + + gmsh::model::addPhysicalGroup(1, {tagsCircle[nDom + 1]}, 93); gmsh::model::setPhysicalName(1, 93, "gammaPml"); - + for(unsigned int i = 0; i < nDom; ++i) { - gmsh::model::addPhysicalGroup(2, {tagsCircle[i]}, i+1); - gmsh::model::setPhysicalName(2, i+1, "omega_" + std::to_string(i)); + gmsh::model::addPhysicalGroup(2, {tagsCircle[i]}, i + 1); + gmsh::model::setPhysicalName(2, i + 1, "omega_" + std::to_string(i)); if(i != nDom - 1) { - gmsh::model::addPhysicalGroup(1, {tagsCircle[i+1]}, 200 + i); + gmsh::model::addPhysicalGroup(1, {tagsCircle[i + 1]}, 200 + i); gmsh::model::setPhysicalName(1, 200 + i, "sigma_" + std::to_string(i)); } } gmsh::model::addPhysicalGroup(2, {tagsCircle[nDom]}, 1 + nDom); gmsh::model::setPhysicalName(2, 1 + nDom, "omega_pml"); - + std::vector< std::pair< int, int > > out; gmsh::model::getEntities(out, 0); gmsh::model::mesh::setSize(out, lc); @@ -81,21 +83,19 @@ int main(int argc, char **argv) unsigned int nDom = 3; //gmshDdm.userDefinedParameter(nDom, "nDom"); - unsigned int MPI_Size = getSize_MPI(); - unsigned int MPI_Rank = getRank_MPI(); - + double rInf = 5.; double r0 = 1; double lc = 0.1; int Npml = 4; - double rPml = rInf+Npml*lc; - + double rPml = rInf + Npml * lc; + gmshDdm.userDefinedParameter(rInf, "rInf"); gmshDdm.userDefinedParameter(r0, "r0"); gmshDdm.userDefinedParameter(lc, "lc"); - - int FEMorder = 3; + + int FEMorder = 3; std::string gauss = "Gauss10"; gmshDdm.userDefinedParameter(gauss, "gauss"); @@ -103,7 +103,7 @@ int main(int argc, char **argv) // transmission conditions std::string Transmission = "Taylor2"; gmshDdm.userDefinedParameter(Transmission, "transmission"); - double alpha = pi/4.; + double alpha = pi / 4.; gmshDdm.userDefinedParameter(alpha, "alpha"); msg::info << " Transmission condition " << Transmission << " with rotation alpha = " << alpha << msg::endl; // Build geometry and mesh using gmsh API @@ -112,16 +112,16 @@ int main(int argc, char **argv) // Define domain Subdomain omega(nDom); Subdomain gammaPml(nDom); - Subdomain gammaInf(nDom); + Subdomain gammaInf(nDom); Subdomain gammaScat(nDom); Interface sigma(nDom); for(unsigned int i = 0; i < nDom; ++i) { // ddm entities - - if(i == (nDom - 1) ) { + + if(i == (nDom - 1)) { gammaPml(i) = Domain(1, 93); gammaInf(i) = Domain(1, 92); - omega(i) = ( Domain(2, i+1) | Domain (2, i+2) ); + omega(i) = (Domain(2, i + 1) | Domain(2, i + 2)); } else { omega(i) = Domain(2, i + 1); @@ -154,15 +154,15 @@ int main(int argc, char **argv) gmshddm::problem::Formulation< std::complex< double > > formulation("helmholtz", topology); double w = 2 * pi; - double theta = pi/3.; + double theta = pi / 3.; // Create fields SubdomainField< std::complex< double >, Form::Form0 > u("u", omega | gammaScat | gammaPml | sigma, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder); InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder); Function< std::complex< double >, Degree::Degree0 > *solution = nullptr; - - solution = new AnalyticalFunction< helmholtz2D::ScatteringByASoftCylinder< std::complex< double > > >(w, r0, 0., 0., 2*w, theta); + + solution = new AnalyticalFunction< helmholtz2D::ScatteringByASoftCylinder< std::complex< double > > >(w, r0, 0., 0., 2 * w, theta); for(unsigned int i = 0; i < nDom; ++i) { // physical source constraint u(i).addConstraint(gammaScat(i), *solution); @@ -170,15 +170,15 @@ int main(int argc, char **argv) solution->activateMemorization(); ScalarFunction< std::complex< double > > S, dS; - if ( Transmission == "Taylor0" ) { - S = -im*w*exp(im*alpha/2.); + if(Transmission == "Taylor0") { + S = -im * w * exp(im * alpha / 2.); dS = 0.; } - else if ( Transmission == "Taylor2" ) { - S = -im*w*cos(alpha/2.); - dS = -im*exp(-im*alpha/2.0)/(2.0*w); + else if(Transmission == "Taylor2") { + S = -im * w * cos(alpha / 2.); + dS = -im * exp(-im * alpha / 2.0) / (2.0 * w); } - + // add interface field formulation.addInterfaceField(g); @@ -186,25 +186,25 @@ int main(int argc, char **argv) const double WPml = rPml - rInf; ScalarFunction< std::complex< double > > cosT = x< std::complex< double > >() / r2d< std::complex< double > >(); ScalarFunction< std::complex< double > > sinT = y< std::complex< double > >() / r2d< std::complex< double > >(); - ScalarFunction< std::complex< double > > dampingProfileR = heaviside( abs(r2d< std::complex< double > >()) - rInf ) / (WPml - (r2d< std::complex< double > >() - rInf)); - ScalarFunction< std::complex< double > > dampingProfileInt = - heaviside( abs(r2d< std::complex< double > >()) - rInf ) * ln((WPml - (r2d< std::complex< double > >() - rInf)) / WPml); + ScalarFunction< std::complex< double > > dampingProfileR = heaviside(abs(r2d< std::complex< double > >()) - rInf) / (WPml - (r2d< std::complex< double > >() - rInf)); + ScalarFunction< std::complex< double > > dampingProfileInt = -heaviside(abs(r2d< std::complex< double > >()) - rInf) * ln((WPml - (r2d< std::complex< double > >() - rInf)) / WPml); ScalarFunction< std::complex< double > > cR = 1. + im * dampingProfileR / w; ScalarFunction< std::complex< double > > cStretch = 1. + im * (1. / r2d< std::complex< double > >()) * dampingProfileInt / w; ScalarFunction< std::complex< double > > S_PML = cR * cStretch; - TensorFunction< std::complex< double > > D = tensor< std::complex <double > >(cStretch / cR * cosT * cosT + cR / cStretch * sinT * sinT, - cStretch / cR * cosT * sinT - cR / cStretch * cosT * sinT, - 0., - cStretch / cR * cosT * sinT - cR / cStretch * cosT * sinT, - cStretch / cR * sinT * sinT + cR / cStretch * cosT * cosT, - 0., - 0., 0., 0.); + TensorFunction< std::complex< double > > D = tensor< std::complex< double > >(cStretch / cR * cosT * cosT + cR / cStretch * sinT * sinT, + cStretch / cR * cosT * sinT - cR / cStretch * cosT * sinT, + 0., + cStretch / cR * cosT * sinT - cR / cStretch * cosT * sinT, + cStretch / cR * sinT * sinT + cR / cStretch * cosT * cosT, + 0., + 0., 0., 0.); // Add terms to the formulation for(unsigned int i = 0; i < nDom; ++i) { // VOLUME TERMS - if ( i == (nDom -1) ) { // pre-processing is slower with the PML formulation + if(i == (nDom - 1)) { // pre-processing is slower with the PML formulation formulation(i).integral(D * grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(- w * w * S_PML * dof(u(i)), tf(u(i)), omega(i), gauss); + formulation(i).integral(-w * w * S_PML * dof(u(i)), tf(u(i)), omega(i), gauss); } else { formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss); @@ -213,18 +213,18 @@ int main(int argc, char **argv) // Artificial source for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(S * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(- dS * grad( dof(u(i)) ), grad( tf(u(i)) ) , sigma(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -g(j,i) ), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(S * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(-dS * grad(dof(u(i))), grad(tf(u(i))), sigma(i, j), gauss); + formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss); } // INTERFACE TERMS for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i,j).integral(dof(g(i,j)), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( g(j,i) ), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(-2.0 * S * u(i), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(2.0 * dS * grad( u(i) ), grad ( tf(g(i,j)) ), sigma(i,j), gauss); + formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(g(j, i)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(-2.0 * S * u(i), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(2.0 * dS * grad(u(i)), grad(tf(g(i, j))), sigma(i, j), gauss); } } @@ -235,14 +235,14 @@ int main(int argc, char **argv) gmshDdm.userDefinedParameter(maxIter, "maxIter"); formulation.pre(); formulation.solve("gmres", tolerence, maxIter); - + // L2 error double num = 0., den = 0.; Subdomain omega_err(nDom); for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { - if ( i % MPI_Size == MPI_Rank ) { - omega_err(i) = Domain (2, i + 1); - save(+u(i), omega(i),"u_"+std::to_string(i)); + if(gmshddm::mpi::isItMySubdomain(i)) { + omega_err(i) = Domain(2, i + 1); + save(+u(i), omega(i), "u_" + std::to_string(i)); save(*solution, omega_err(i), "uex_" + std::to_string(i)); save(*solution - u(i), omega_err(i), "err_" + std::to_string(i)); diff --git a/demos/helmholtz/checkerboard_bcgs/CMakeLists.txt b/examples/helmholtzflow/freefield/CMakeLists.txt similarity index 73% rename from demos/helmholtz/checkerboard_bcgs/CMakeLists.txt rename to examples/helmholtzflow/freefield/CMakeLists.txt index 4d2e4b0c91d2a7e0da5a224e581c56eca24b2f71..c368ce59f63aed4e569c1b3e727df6dd0b1d6a45 100644 --- a/demos/helmholtz/checkerboard_bcgs/CMakeLists.txt +++ b/examples/helmholtzflow/freefield/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.0 FATAL_ERROR) project(demo CXX) -include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake) +include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake) add_executable(demo main.cpp ${EXTRA_INCS}) target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/demos/helmholtzflow/freefield/main.cpp b/examples/helmholtzflow/freefield/main.cpp similarity index 63% rename from demos/helmholtzflow/freefield/main.cpp rename to examples/helmholtzflow/freefield/main.cpp index 373429a682c851d3b2f20721cda2fb4ee0f946c3..36a41a07dd2d2fcd3e0f7f2e5c33fc69d58b09a6 100644 --- a/demos/helmholtzflow/freefield/main.cpp +++ b/examples/helmholtzflow/freefield/main.cpp @@ -1,8 +1,9 @@ -#include <gmshddm/GmshDdm.h> -#include <gmshddm/Formulation.h> -#include <gmshfem/Message.h> #include "gmsh.h" + +#include <gmshddm/Formulation.h> +#include <gmshddm/GmshDdm.h> #include <gmshfem/AnalyticalFunction.h> +#include <gmshfem/Message.h> using namespace gmshddm; using namespace gmshddm::common; @@ -37,9 +38,9 @@ void circlesConcentric(const unsigned int nDom, const double rInf, const double std::vector< int > tagsCircle(nDom + 1); for(unsigned int i = 0; i < nDom; ++i) { - tagsCircle[i] = gmsh::model::occ::addCircle(0., 0., 0., std::sqrt(double( (i+1) * (rInf * rInf)) / nDom )); + tagsCircle[i] = gmsh::model::occ::addCircle(0., 0., 0., std::sqrt(double((i + 1) * (rInf * rInf)) / nDom)); } - tagsCircle[nDom] = gmsh::model::occ::addCircle(0., 0., 0., std::sqrt( double(rPml * rPml) ) ); + tagsCircle[nDom] = gmsh::model::occ::addCircle(0., 0., 0., std::sqrt(double(rPml * rPml))); std::vector< int > tagsCurve(nDom + 1); for(unsigned int i = 0; i < (nDom + 1); ++i) { @@ -48,12 +49,12 @@ void circlesConcentric(const unsigned int nDom, const double rInf, const double std::vector< int > tagsSurface(nDom + 1); tagsSurface[0] = gmsh::model::occ::addPlaneSurface({tagsCurve[0]}); - for(unsigned int i = 1; i < (nDom+1); ++i) { - tagsSurface[i] = gmsh::model::occ::addPlaneSurface({tagsCurve[i-1], tagsCurve[i]}); + for(unsigned int i = 1; i < (nDom + 1); ++i) { + tagsSurface[i] = gmsh::model::occ::addPlaneSurface({tagsCurve[i - 1], tagsCurve[i]}); } gmsh::model::occ::synchronize(); - gmsh::model::mesh::embed(0,{tagPoint},2,1); + gmsh::model::mesh::embed(0, {tagPoint}, 2, 1); for(unsigned int i = 0; i < nDom; ++i) { gmsh::model::addPhysicalGroup(2, {tagsCircle[i]}, 100 + i); @@ -65,7 +66,7 @@ void circlesConcentric(const unsigned int nDom, const double rInf, const double } } - gmsh::model::addPhysicalGroup(1, {tagsCircle[nDom-1]}, 92); + gmsh::model::addPhysicalGroup(1, {tagsCircle[nDom - 1]}, 92); gmsh::model::setPhysicalName(1, 92, "gammaInf"); gmsh::model::addPhysicalGroup(1, {tagsCircle[nDom]}, 93); @@ -95,18 +96,18 @@ int main(int argc, char **argv) unsigned int nDom = 5; // gmshDdm.userDefinedParameter(nDom, "nDom"); const double pi = 3.14159265358979323846264338327950288; - + // geometry parameters double r0 = 2.; double lc = 0.07; - double rPML = r0+4.*lc; + double rPML = r0 + 4. * lc; gmshDdm.userDefinedParameter(rPML, "rPML"); gmshDdm.userDefinedParameter(r0, "r0"); gmshDdm.userDefinedParameter(lc, "lc"); // souce point location double xs = -0.1; double ys = 0.2; - + // numerical parameters int FEMorder = 4; gmshDdm.userDefinedParameter(FEMorder, "FEMorder"); @@ -133,13 +134,13 @@ int main(int argc, char **argv) if(i == (nDom - 1)) { gammaPml(i) = Domain(1, 93); gammaInf(i) = Domain(1, 92); - omega(i) = ( Domain(2, i + 100) | Domain (2, 100 + i + 1) ); + omega(i) = (Domain(2, i + 100) | Domain(2, 100 + i + 1)); } else { omega(i) = Domain(2, i + 100); } - if ( i == 0 ) { + if(i == 0) { source(i) = Domain(0, 999); } @@ -167,30 +168,30 @@ int main(int argc, char **argv) const std::complex< double > im(0., 1.); double M = 0.5; // Mach number gmshDdm.userDefinedParameter(M, "M"); - double theta = pi/3.; // mean flow orientation + double theta = pi / 3.; // mean flow orientation msg::info << "Mach number " << M << " theta " << theta << msg::endl; - double Mx = M*std::cos(theta); - double My = M*std::sin(theta); + double Mx = M * std::cos(theta); + double My = M * std::sin(theta); // Parameters for Inverse Lorentz transformation - double beta = sqrt(1-M*M); // Jacobian - double Alphax = 1 + Mx*Mx/(beta*(1+beta)); - double Alphay = 1 + My*My/(beta*(1+beta)); - double K = Mx*My/(beta*(1+beta)); + double beta = sqrt(1 - M * M); // Jacobian + double Alphax = 1 + Mx * Mx / (beta * (1 + beta)); + double Alphay = 1 + My * My / (beta * (1 + beta)); + double K = Mx * My / (beta * (1 + beta)); // Inverse Lorentz transformation - TensorFunction< std::complex< double > > Linv = tensor< std::complex <double > > (beta*Alphay, -beta*K, 0., -beta*K, beta*Alphax,0.,0.,0.,0.); + TensorFunction< std::complex< double > > Linv = tensor< std::complex< double > >(beta * Alphay, -beta * K, 0., -beta * K, beta * Alphax, 0., 0., 0., 0.); // Mach-velocity vector - VectorFunction< std::complex< double > > MM = vector< std::complex< double > > (Mx,My,0.); - ScalarFunction< std::complex< double > > Mn = MM*normal< std::complex< double > >(); - ScalarFunction< std::complex< double > > Mt = MM*tangent< std::complex< double > >(); - ScalarFunction< std::complex< double > > beta_n = sqrt(1-Mn*Mn); // Jacobian + VectorFunction< std::complex< double > > MM = vector< std::complex< double > >(Mx, My, 0.); + ScalarFunction< std::complex< double > > Mn = MM * normal< std::complex< double > >(); + ScalarFunction< std::complex< double > > Mt = MM * tangent< std::complex< double > >(); + ScalarFunction< std::complex< double > > beta_n = sqrt(1 - Mn * Mn); // Jacobian // Transmission condition ScalarFunction< std::complex< double > > Si, Sj, Sti, Stj, dS; - float pointsByWl = ( 2*pi*FEMorder*( 1 - abs(M) ) ) / (lc*k); + float pointsByWl = (2 * pi * FEMorder * (1 - abs(M))) / (lc * k); msg::info << " - FEM basis order = " << FEMorder << "" << msg::endl; msg::info << " - Minimum number of dofs by wavelength = " << pointsByWl << msg::endl; @@ -207,7 +208,7 @@ int main(int argc, char **argv) Function< std::complex< double >, Degree::Degree0 > *solution = nullptr; solution = new AnalyticalFunction< helmholtz2D::MonopoleFreeField< std::complex< double > > >(k, M, theta, 1., xs, ys, 0., 0.); solution->activateMemorization(); - + // PML parameters double Sigma0 = beta; ScalarFunction< std::complex< double > > det_J; @@ -216,104 +217,104 @@ int main(int argc, char **argv) ScalarFunction< std::complex< double > > cosT = x< std::complex< double > >() / r2d< std::complex< double > >(); ScalarFunction< std::complex< double > > sinT = y< std::complex< double > >() / r2d< std::complex< double > >(); - ScalarFunction< std::complex< double > > dampingProfileR = Sigma0*heaviside( abs(r2d< std::complex< double > >()) - r0 ) / (Wpml - (r2d< std::complex< double > >() - r0)); - ScalarFunction< std::complex< double > > dampingProfileInt = -Sigma0*heaviside( abs(r2d< std::complex< double > >()) - r0 ) * ln((Wpml - (r2d< std::complex< double > >() - r0)) / Wpml); + ScalarFunction< std::complex< double > > dampingProfileR = Sigma0 * heaviside(abs(r2d< std::complex< double > >()) - r0) / (Wpml - (r2d< std::complex< double > >() - r0)); + ScalarFunction< std::complex< double > > dampingProfileInt = -Sigma0 * heaviside(abs(r2d< std::complex< double > >()) - r0) * ln((Wpml - (r2d< std::complex< double > >() - r0)) / Wpml); ScalarFunction< std::complex< double > > gamma = 1. - im * dampingProfileR / k; ScalarFunction< std::complex< double > > gamma_hat = 1. - im * (1. / r2d< std::complex< double > >()) * dampingProfileInt / k; det_J = gamma * gamma_hat; - J_PML_inv_T = tensor< std::complex <double > >(cosT/gamma,sinT/gamma,0.,-sinT/gamma_hat,cosT/gamma_hat,0.,0.,0.,0.); + J_PML_inv_T = tensor< std::complex< double > >(cosT / gamma, sinT / gamma, 0., -sinT / gamma_hat, cosT / gamma_hat, 0., 0., 0., 0.); // build vector and matrix for the weak form VectorFunction< std::complex< double > > J_PML_inv_T_M = J_PML_inv_T * MM; - TensorFunction< std::complex< double > > J_PML_Linv = J_PML_inv_T * Linv; + TensorFunction< std::complex< double > > J_PML_Linv = J_PML_inv_T * Linv; msg::info << "Use a PML as radiation condition - circular stable formulation" << msg::endl; for(unsigned int i = 0; i < nDom; ++i) { // VOLUME TERMS - if ( i != (nDom - 1) ) { + if(i != (nDom - 1)) { // convected Helmholz weak form - formulation(i).integral(vector< std::complex< double > >(1-Mx*Mx, -Mx*My, 0.) * grad(dof(u(i))), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(vector< std::complex< double > >(-Mx*My, 1-My*My, 0.) * grad(dof(u(i))), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(dof(u(i)), vector< std::complex< double > >(-im*k*Mx,-im*k*My, 0.) * grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(vector< std::complex< double > >(im*k*Mx, im*k*My, 0.) * grad(dof(u(i))), tf(u(i)), omega(i), gauss); - formulation(i).integral(- k * k * dof(u(i)), tf(u(i)), omega(i), gauss); + formulation(i).integral(vector< std::complex< double > >(1 - Mx * Mx, -Mx * My, 0.) * grad(dof(u(i))), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u(i))), omega(i), gauss); + formulation(i).integral(vector< std::complex< double > >(-Mx * My, 1 - My * My, 0.) * grad(dof(u(i))), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u(i))), omega(i), gauss); + formulation(i).integral(dof(u(i)), vector< std::complex< double > >(-im * k * Mx, -im * k * My, 0.) * grad(tf(u(i))), omega(i), gauss); + formulation(i).integral(vector< std::complex< double > >(im * k * Mx, im * k * My, 0.) * grad(dof(u(i))), tf(u(i)), omega(i), gauss); + formulation(i).integral(-k * k * dof(u(i)), tf(u(i)), omega(i), gauss); } else { // convected Helmholz weak form with PML - formulation(i).integral( det_J * J_PML_Linv*grad(dof(u(i))) , J_PML_Linv*grad(tf(u(i))) , omega(i), gauss); + formulation(i).integral(det_J * J_PML_Linv * grad(dof(u(i))), J_PML_Linv * grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(+ k*k/(beta*beta) * det_J * J_PML_inv_T_M * dof(u(i)) , J_PML_inv_T_M * tf(u(i)), omega(i), gauss); - formulation(i).integral(+ im*k/beta* det_J * J_PML_Linv * grad(dof(u(i))), J_PML_inv_T_M * tf(u(i)), omega(i), gauss); + formulation(i).integral(+k * k / (beta * beta) * det_J * J_PML_inv_T_M * dof(u(i)), J_PML_inv_T_M * tf(u(i)), omega(i), gauss); + formulation(i).integral(+im * k / beta * det_J * J_PML_Linv * grad(dof(u(i))), J_PML_inv_T_M * tf(u(i)), omega(i), gauss); - formulation(i).integral(- im*k/beta* det_J * J_PML_inv_T_M * dof(u(i)), J_PML_Linv * grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(- k*k/(beta*beta) * det_J * dof(u(i)), tf(u(i)), omega(i), gauss); + formulation(i).integral(-im * k / beta * det_J * J_PML_inv_T_M * dof(u(i)), J_PML_Linv * grad(tf(u(i))), omega(i), gauss); + formulation(i).integral(-k * k / (beta * beta) * det_J * dof(u(i)), tf(u(i)), omega(i), gauss); } // point source formulation(i).integral(formulation.physicalSource(-1.), tf(u(i)), source(i), gauss); - + // Artificial source - Zeroth order transmission condition for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(formulation.artificialSource( -g(j,i) ), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss); // Contribution from the boundary terms - formulation(i).integral( (im*k*Mn) * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(Mn * Mt * tangent< std::complex< double > >() * grad( dof(u(i)) ), dof(u(i)), sigma(i,j), gauss); + formulation(i).integral((im * k * Mn) * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(Mn * Mt * tangent< std::complex< double > >() * grad(dof(u(i))), dof(u(i)), sigma(i, j), gauss); } // Interface weak form for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i,j).integral(dof(g(i,j)), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( g(j,i) ), tf(g(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(g(j, i)), tf(g(i, j)), sigma(i, j), gauss); } - if ( transmission == "T0" ) { - Si = im*k*(exp(im*alpha/2.)-Mn); // alpha=0, beta*beta*Si = i*k*(1-Mn) - Sti = - Mn * Mt; - Sj = im*k*(exp(im*alpha/2.)+Mn); + if(transmission == "T0") { + Si = im * k * (exp(im * alpha / 2.) - Mn); // alpha=0, beta*beta*Si = i*k*(1-Mn) + Sti = -Mn * Mt; + Sj = im * k * (exp(im * alpha / 2.) + Mn); for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral( Si * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral( Sti * tangent< std::complex< double > >() * grad( dof(u(i)) ), dof(u(i)), sigma(i,j), gauss); - formulation(i,j).integral(- (Si+Sj) * u(i) , tf(g(i,j)), sigma(i,j), gauss); - } + formulation(i).integral(Si * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(Sti * tangent< std::complex< double > >() * grad(dof(u(i))), dof(u(i)), sigma(i, j), gauss); + formulation(i, j).integral(-(Si + Sj) * u(i), tf(g(i, j)), sigma(i, j), gauss); + } } - else if ( transmission == "T2" ) { - Si = im*k*(cos(alpha/2.)-Mn) ; - Sj = im*k*(cos(alpha/2.)+Mn) ; - Sti = ( exp(-im*alpha/2.) - Mn ) * Mt ; - Stj = ( exp(-im*alpha/2.) + Mn ) * Mt ; - dS = im*(beta*beta)*exp(-im*alpha/2.) / (2.0*k); - + else if(transmission == "T2") { + Si = im * k * (cos(alpha / 2.) - Mn); + Sj = im * k * (cos(alpha / 2.) + Mn); + Sti = (exp(-im * alpha / 2.) - Mn) * Mt; + Stj = (exp(-im * alpha / 2.) + Mn) * Mt; + dS = im * (beta * beta) * exp(-im * alpha / 2.) / (2.0 * k); + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral( Si * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral( Sti * tangent< std::complex< double > >() * grad( dof(u(i)) ), dof(u(i)), sigma(i,j), gauss); - formulation(i).integral(- dS * grad( dof(u(i)) ), grad( tf(u(i)) ), sigma(i,j), gauss); - - formulation(i,j).integral(- (Si+Sj) * u(i) , tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(- (Sti+Stj) * tangent< std::complex< double > >() * grad( u(i) ) , tf(g(i,j)), sigma(i,j), gauss); - - formulation(i,j).integral( 2. * dS * grad( u(i) ), grad( tf(g(i,j)) ), sigma(i,j), gauss); + formulation(i).integral(Si * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(Sti * tangent< std::complex< double > >() * grad(dof(u(i))), dof(u(i)), sigma(i, j), gauss); + formulation(i).integral(-dS * grad(dof(u(i))), grad(tf(u(i))), sigma(i, j), gauss); + + formulation(i, j).integral(-(Si + Sj) * u(i), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(-(Sti + Stj) * tangent< std::complex< double > >() * grad(u(i)), tf(g(i, j)), sigma(i, j), gauss); + + formulation(i, j).integral(2. * dS * grad(u(i)), grad(tf(g(i, j))), sigma(i, j), gauss); } } - else if ( transmission == "Pade" ) { + else if(transmission == "Pade") { int padeOrder = 2; gmshDdm.userDefinedParameter(padeOrder, "padeOrder"); const double Np = 2. * padeOrder + 1.; - const std::complex< double > exp1 = std::complex<double>(std::cos(alpha),std::sin(alpha)); - const std::complex< double > exp2 = std::complex<double>(std::cos(alpha/2.),std::sin(alpha/2.)); - const std::complex< double > coef = 2./Np; - + const std::complex< double > exp1 = std::complex< double >(std::cos(alpha), std::sin(alpha)); + const std::complex< double > exp2 = std::complex< double >(std::cos(alpha / 2.), std::sin(alpha / 2.)); + const std::complex< double > coef = 2. / Np; + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral( im * k * exp2 * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(- (im*Mn*k) * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(- Mn*Mt * tangent< std::complex< double > >() * grad( dof(u(i)) ), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(im * k * exp2 * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(-(im * Mn * k) * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(-Mn * Mt * tangent< std::complex< double > >() * grad(dof(u(i))), tf(u(i)), sigma(i, j), gauss); // interface unknowns - formulation(i,j).integral(-2.*(im * k * exp2) * u(i), tf(g(i,j)), sigma(i,j), gauss); - + formulation(i, j).integral(-2. * (im * k * exp2) * u(i), tf(g(i, j)), sigma(i, j), gauss); + // define Pade coefficients std::vector< std::complex< double > > c(padeOrder, 0.); for(int l = 0; l < padeOrder; ++l) { @@ -322,27 +323,27 @@ int main(int argc, char **argv) } // define the auxiliary fields - std::vector< Field< std::complex< double >, Form::Form0 >* > phi; + std::vector< Field< std::complex< double >, Form::Form0 > * > phi; for(int l = 0; l < padeOrder; ++l) { - phi.push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(l), sigma(i,j), FunctionSpaceTypeForm0::HierarchicalH1, FEMorder)); + phi.push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(l), sigma(i, j), FunctionSpaceTypeForm0::HierarchicalH1, FEMorder)); fieldBucket.push_back(phi.back()); } - + // write the augmented weak form for(int l = 0; l < padeOrder; ++l) { // boundary integral terms relating the auxiliary fields - formulation(i).integral(im * k * exp2 * coef * c[l] * dof(*phi[l]), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(im * k * exp2 * coef * c[l] * dof(u(i)), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(im * k * exp2 * coef * c[l] * dof(*phi[l]), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(im * k * exp2 * coef * c[l] * dof(u(i)), tf(u(i)), sigma(i, j), gauss); // coupling of the auxiliary equations - formulation(i).integral(- (beta*beta) * grad(dof(*phi[l])), grad(tf(*phi[l])), sigma(i,j), gauss); - formulation(i).integral(- 2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[l])), tf(*phi[l]), sigma(i,j), gauss); - formulation(i).integral( (k*k) * (exp1 * c[l] + 1.) * dof(*phi[l]), tf(*phi[l]), sigma(i,j), gauss); - formulation(i).integral( (k*k) * exp1 * (c[l] + 1.) * dof(u(i)), tf(*phi[l]), sigma(i,j), gauss); + formulation(i).integral(-(beta * beta) * grad(dof(*phi[l])), grad(tf(*phi[l])), sigma(i, j), gauss); + formulation(i).integral(-2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[l])), tf(*phi[l]), sigma(i, j), gauss); + formulation(i).integral((k * k) * (exp1 * c[l] + 1.) * dof(*phi[l]), tf(*phi[l]), sigma(i, j), gauss); + formulation(i).integral((k * k) * exp1 * (c[l] + 1.) * dof(u(i)), tf(*phi[l]), sigma(i, j), gauss); } // interfaces unknowns for(int l = 0; l < padeOrder; ++l) { - formulation(i,j).integral(- 2. * (im * k * exp2) * coef * c[l] * (*phi[l]), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(- 2. * (im * k * exp2) * coef * c[l] * u(i), tf(g(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(-2. * (im * k * exp2) * coef * c[l] * (*phi[l]), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(-2. * (im * k * exp2) * coef * c[l] * u(i), tf(g(i, j)), sigma(i, j), gauss); } } } @@ -355,8 +356,8 @@ int main(int argc, char **argv) double num = 0., den = 0.; // L2 error Subdomain omega_err(nDom); for(unsigned int i = 0; i < nDom; ++i) { - omega_err(i) = Domain(2, i+100); - save(+u(i), omega(i),"u_"+std::to_string(i)); + omega_err(i) = Domain(2, i + 100); + save(+u(i), omega(i), "u_" + std::to_string(i)); save(*solution, omega_err(i), "uex_" + std::to_string(i)); save(*solution - u(i), omega_err(i), "err_" + std::to_string(i)); diff --git a/examples/helmholtzflow/freefield/runTests.sh b/examples/helmholtzflow/freefield/runTests.sh new file mode 100644 index 0000000000000000000000000000000000000000..5486cd76fcd8942657652e9ad212faaeed70004a --- /dev/null +++ b/examples/helmholtzflow/freefield/runTests.sh @@ -0,0 +1,7 @@ +cd build + cmake..&&make + + echo "****************************" echo " zeroth order TC " echo "****************************" + ./ + demo - + transmission T0 - alpha 0 - M 0.8 echo "****************************" echo " zeroth order TC with rotation " echo "****************************"./ demo - transmission T0 - alpha - 1 - M 0.8 echo "****************************" echo " second order TC with rotation " echo "****************************"./ demo - transmission T2 - alpha - 1 - M 0.8 echo "****************************" echo " Pade TC with rotation " echo "****************************"./ demo - transmission Pade - padeOrder 4 - alpha - 1 - M 0.8 diff --git a/demos/helmholtzflow/freefield/CMakeLists.txt b/examples/helmholtzflow/waveguide/CMakeLists.txt similarity index 73% rename from demos/helmholtzflow/freefield/CMakeLists.txt rename to examples/helmholtzflow/waveguide/CMakeLists.txt index 4d2e4b0c91d2a7e0da5a224e581c56eca24b2f71..c368ce59f63aed4e569c1b3e727df6dd0b1d6a45 100644 --- a/demos/helmholtzflow/freefield/CMakeLists.txt +++ b/examples/helmholtzflow/waveguide/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.0 FATAL_ERROR) project(demo CXX) -include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake) +include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake) add_executable(demo main.cpp ${EXTRA_INCS}) target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/demos/helmholtzflow/waveguide/main.cpp b/examples/helmholtzflow/waveguide/main.cpp similarity index 64% rename from demos/helmholtzflow/waveguide/main.cpp rename to examples/helmholtzflow/waveguide/main.cpp index 34e495ea5f12cdfb7913f576b210682c7944c1e6..454510b8861f93faecba63a33f80eafd979ffa62 100644 --- a/demos/helmholtzflow/waveguide/main.cpp +++ b/examples/helmholtzflow/waveguide/main.cpp @@ -1,10 +1,11 @@ -#include <gmshddm/GmshDdm.h> #include <gmshddm/Formulation.h> +#include <gmshddm/GmshDdm.h> +#include <gmshddm/MPIInterface.h> #include <gmshfem/Message.h> // #include "mpi.h" // #include "petscsys.h" -#include "gmsh.h" -#include "AnalyticalFunction.h" +#include <gmsh.h> +#include <gmshfem/AnalyticalFunction.h> using namespace gmshddm; using namespace gmshddm::common; @@ -81,21 +82,16 @@ int main(int argc, char **argv) { GmshDdm gmshDdm(argc, argv); - // rank and size - unsigned int MPI_Rank = getRank_MPI(); - unsigned int MPI_Size = getSize_MPI(); - printf("I am processor of rank %d out of %d processes \n", MPI_Rank, MPI_Size); - double pi = 3.14159265359; int nDom = 4; gmshDdm.userDefinedParameter(nDom, "nDom"); - + // duct geometry parameters double L = 0.5; double H = 0.25; - double lc = sqrt(L*H)/(40.0); + double lc = sqrt(L * H) / (40.0); gmshDdm.userDefinedParameter(lc, "lc"); - + // upper wall boundary condition std::string wallType = "hard"; gmshDdm.userDefinedParameter(wallType, "wallType"); @@ -104,33 +100,31 @@ int main(int argc, char **argv) double k = 70; gmshDdm.userDefinedParameter(k, "k"); int mode = 3; - double ky = mode*pi/H; - + double ky = mode * pi / H; + // relative Mach number double M = 0.8; gmshDdm.userDefinedParameter(M, "M"); - double beta = sqrt(1-M*M); - double beta2 = beta*beta; - + double beta = sqrt(1 - M * M); + double beta2 = beta * beta; + // transmission conditions std::string Transmission = "DtN"; gmshDdm.userDefinedParameter(Transmission, "transmission"); // init transmission operators ScalarFunction< std::complex< double > > Si, Sj, dS; - double alpha=0.; // rotation branch-cut, in [0, -pi] + double alpha = 0.; // rotation branch-cut, in [0, -pi] gmshDdm.userDefinedParameter(alpha, "alpha"); - if ( MPI_Rank == 0 ) { - if (Transmission != "DtN") { - msg::info << " use " << Transmission << " as transmission condition with alpha = " << alpha << msg::endl; - } - else { - msg::info << " use analytical DtN as transmission condition " << msg::endl; - } - if ( wallType == "soft" && Transmission == "Pade" ) { - msg::error << "Pade condition not available yet for Dirichlet wall " << msg::endl; - } + if(Transmission != "DtN") { + msg::info << " use " << Transmission << " as transmission condition with alpha = " << alpha << msg::endl; + } + else { + msg::info << " use analytical DtN as transmission condition " << msg::endl; + } + if(wallType == "soft" && Transmission == "Pade") { + msg::error << "Pade condition not available yet for Dirichlet wall " << msg::endl; } - + // numerical parameters int FEMorder = 4; gmshDdm.userDefinedParameter(FEMorder, "FEMorder"); @@ -139,7 +133,7 @@ int main(int argc, char **argv) // Build geometry and mesh using gmsh API // for the moment all the procs do the meshing - waveguide(nDom, H, L, lc); + waveguide(nDom, H, L, lc); // Define domain Subdomain omega(nDom); @@ -150,7 +144,7 @@ int main(int argc, char **argv) for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { omega(i) = Domain(2, i + 1); - if(i == (static_cast< unsigned int >(nDom-1)) ) { + if(i == (static_cast< unsigned int >(nDom - 1))) { gammaInf(i) = Domain(1, 2); // outgoing boundary } @@ -162,7 +156,7 @@ int main(int argc, char **argv) if(i != 0) { sigma(i, i - 1) = Domain(1, 200 + i - 1); } - if(i != (static_cast< unsigned int >(nDom-1)) ) { + if(i != (static_cast< unsigned int >(nDom - 1))) { sigma(i, i + 1) = Domain(1, 200 + i); } } @@ -174,9 +168,9 @@ int main(int argc, char **argv) if(i != 0) { topology[i].push_back(i - 1); } - if(i != (static_cast< unsigned int >(nDom-1)) ) { + if(i != (static_cast< unsigned int >(nDom - 1))) { topology[i].push_back(i + 1); - } + } } // Create DDM formulation @@ -186,15 +180,15 @@ int main(int argc, char **argv) // Create fields SubdomainField< std::complex< double >, Form::Form0 > u("u", omega | gammaDir | gammaInf | border | sigma, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder); InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder); - + // define analytical solution Function< std::complex< double >, Degree::Degree0 > *solution = nullptr; - if ( wallType == "hard" ) { + if(wallType == "hard") { solution = new AnalyticalFunction< helmholtz2D::DuctModeSolution< std::complex< double > > >(k, M, H, 0., 0., mode, 0); } - else if ( wallType == "soft" ) { + else if(wallType == "soft") { solution = new AnalyticalFunction< helmholtz2D::DuctModeSolution< std::complex< double > > >(k, M, H, 0., 0., mode, 1); - for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { + for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { u(i).addConstraint(border(i), 0.); } } @@ -202,124 +196,118 @@ int main(int argc, char **argv) std::complex< double > im(0., 1.); std::complex< double > kx; // propagative wavenumber - double Argsqrt = k*k-beta2*ky*ky; - if (Argsqrt >= 0) { - kx = (1/(beta2))*(-M*k + sqrt(Argsqrt)); - if ( MPI_Rank == 0 ) { - msg::info << " propagative mode " << msg::endl; - } + double Argsqrt = k * k - beta2 * ky * ky; + if(Argsqrt >= 0) { + kx = (1 / (beta2)) * (-M * k + sqrt(Argsqrt)); + msg::info << " propagative mode " << msg::endl; } else { - kx = (1/(beta2))*(-M*k - im*sqrt(abs(Argsqrt))); - if ( MPI_Rank == 0 ) { - msg::info << " evanescent mode " << msg::endl; - } + kx = (1 / (beta2)) * (-M * k - im * sqrt(abs(Argsqrt))); + msg::info << " evanescent mode " << msg::endl; } - if ( k>(beta*ky) && k<ky ) { + if(k > (beta * ky) && k < ky) { msg::info << " - inverse upstream mode ! " << msg::endl; } - float pointsByWl = (2*pi*FEMorder*(1+M)/k) / lc; - if ( MPI_Rank == 0 ) { - msg::info << " - Approximate dofs by wavelength = " << pointsByWl << "" << msg::endl; - msg::info << " - Mach number = " << M << "" << msg::endl; - } - + float pointsByWl = (2 * pi * FEMorder * (1 + M) / k) / lc; + msg::info << " - Approximate dofs by wavelength = " << pointsByWl << "" << msg::endl; + msg::info << " - Mach number = " << M << "" << msg::endl; + // define output absorbing boundary condition - std::complex< double > ABC = im*kx; - + std::complex< double > ABC = im * kx; + // Tell to the formulation that g is field that have to be exchanged between subdomains formulation.addInterfaceField(g); - for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { + for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { // VOLUME TERMS formulation(i).integral(vector< std::complex< double > >(beta2, 0., 0.) * grad(dof(u(i))), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u(i))), omega(i), gauss); formulation(i).integral(vector< std::complex< double > >(0., 1., 0.) * grad(dof(u(i))), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(vector< std::complex< double > >(im*k*M, 0., 0.) * grad(dof(u(i))), tf(u(i)), omega(i), gauss); - formulation(i).integral(dof(u(i)), vector< std::complex< double > >(-im*k*M, 0., 0.) * grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(- k * k * dof(u(i)), tf(u(i)), omega(i), gauss); + formulation(i).integral(vector< std::complex< double > >(im * k * M, 0., 0.) * grad(dof(u(i))), tf(u(i)), omega(i), gauss); + formulation(i).integral(dof(u(i)), vector< std::complex< double > >(-im * k * M, 0., 0.) * grad(tf(u(i))), omega(i), gauss); + formulation(i).integral(-k * k * dof(u(i)), tf(u(i)), omega(i), gauss); // SURFACE TERMS // artificial boundary condition, use analytical DtN formulation(i).integral(beta2 * ABC * dof(u(i)), tf(u(i)), gammaInf(i), gauss); // Contribution from the mass boundary matrix - formulation(i).integral( (im*k*M) * dof(u(i)), tf(u(i)), gammaInf(i), gauss); + formulation(i).integral((im * k * M) * dof(u(i)), tf(u(i)), gammaInf(i), gauss); // Physical source - if ( wallType == "soft" ) { - formulation(i).integral(formulation.physicalSource( -im*kx*beta2*sin< std::complex< double > >(ky * y< std::complex< double > >() ) ), tf(u(i)), gammaDir(i), gauss); + if(wallType == "soft") { + formulation(i).integral(formulation.physicalSource(-im * kx * beta2 * sin< std::complex< double > >(ky * y< std::complex< double > >())), tf(u(i)), gammaDir(i), gauss); } - else if ( wallType == "hard" ) { - formulation(i).integral(formulation.physicalSource( -im*kx*beta2*cos< std::complex< double > >(ky * y< std::complex< double > >() ) ), tf(u(i)), gammaDir(i), gauss); + else if(wallType == "hard") { + formulation(i).integral(formulation.physicalSource(-im * kx * beta2 * cos< std::complex< double > >(ky * y< std::complex< double > >())), tf(u(i)), gammaDir(i), gauss); } // Contribution from the mass boundary matrix - formulation(i).integral( (-im*k*M) * dof(u(i)), tf(u(i)), gammaDir(i), gauss); + formulation(i).integral((-im * k * M) * dof(u(i)), tf(u(i)), gammaDir(i), gauss); // Artificial source for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(formulation.artificialSource( -g(j,i) ), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss); // Contribution from the mass boundary matrix - formulation(i).integral( (im*k*M) * dof(u(i)), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral((im * k * M) * dof(u(i)), tf(u(i)), sigma(i, j), gauss); } // INTERFACE TERMS for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i,j).integral(dof(g(i,j)), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( g(j,i) ), tf(g(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(g(j, i)), tf(g(i, j)), sigma(i, j), gauss); } // Transmission condition - if ( Transmission == "Taylor0" ) { - Si = im*k*(exp(im*alpha/2.)-M)/beta2; - Sj = im*k*(exp(im*alpha/2.)+M)/beta2; + if(Transmission == "Taylor0") { + Si = im * k * (exp(im * alpha / 2.) - M) / beta2; + Sj = im * k * (exp(im * alpha / 2.) + M) / beta2; for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(beta2 * Si * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i,j).integral(- beta2 * (Si+Sj) * u(i), tf(g(i,j)), sigma(i,j), gauss); + formulation(i).integral(beta2 * Si * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i, j).integral(-beta2 * (Si + Sj) * u(i), tf(g(i, j)), sigma(i, j), gauss); } } - else if ( Transmission == "Taylor2" ) { - Si = im*k*(cos(alpha/2.)-M)/beta2; - Sj = im*k*(cos(alpha/2.)+M)/beta2; - dS = im*exp(-im*alpha/2.0)/(2.0*k); + else if(Transmission == "Taylor2") { + Si = im * k * (cos(alpha / 2.) - M) / beta2; + Sj = im * k * (cos(alpha / 2.) + M) / beta2; + dS = im * exp(-im * alpha / 2.0) / (2.0 * k); for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(beta2 * Si * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(- beta2 * dS * grad( dof(u(i)) ), grad( tf(u(i)) ), sigma(i,j), gauss); - formulation(i,j).integral(- beta2 * (Si+Sj) * u(i), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(+ beta2 * (2.*dS) * grad( u(i) ), grad( tf(g(i,j)) ), sigma(i,j), gauss); + formulation(i).integral(beta2 * Si * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(-beta2 * dS * grad(dof(u(i))), grad(tf(u(i))), sigma(i, j), gauss); + formulation(i, j).integral(-beta2 * (Si + Sj) * u(i), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(+beta2 * (2. * dS) * grad(u(i)), grad(tf(g(i, j))), sigma(i, j), gauss); } } - else if ( Transmission == "DtN" ) { // analytic DtN relative the normal of the artificial boundary, for a given mode - Si = im*kx; - Sj = im*( kx + 2*M*k/(beta2) ); // or Sj = im*( (M*k + sqrt(Argsqrt)) / beta2 ); + else if(Transmission == "DtN") { // analytic DtN relative the normal of the artificial boundary, for a given mode + Si = im * kx; + Sj = im * (kx + 2 * M * k / (beta2)); // or Sj = im*( (M*k + sqrt(Argsqrt)) / beta2 ); for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(beta2 * Si * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i,j).integral(- beta2 * (Si+Sj) * u(i), tf(g(i,j)), sigma(i,j), gauss); + formulation(i).integral(beta2 * Si * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i, j).integral(-beta2 * (Si + Sj) * u(i), tf(g(i, j)), sigma(i, j), gauss); } } - else if ( Transmission == "Pade" ) { + else if(Transmission == "Pade") { int padeOrder = 4; gmshDdm.userDefinedParameter(padeOrder, "padeOrder"); - if ( i == 0 && MPI_Rank == 0 ) { + if(i == 0) { msg::info << "Use Pade transmission condition of order " << padeOrder << msg::endl; } const double Np = 2. * padeOrder + 1.; - const std::complex< double > exp1 = std::complex<double>(std::cos(alpha),std::sin(alpha)); - const std::complex< double > exp2 = std::complex<double>(std::cos(alpha/2.),std::sin(alpha/2.)); - const std::complex< double > coef = 2./Np; - + const std::complex< double > exp1 = std::complex< double >(std::cos(alpha), std::sin(alpha)); + const std::complex< double > exp2 = std::complex< double >(std::cos(alpha / 2.), std::sin(alpha / 2.)); + const std::complex< double > coef = 2. / Np; + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral( im * k * exp2 * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(- (im*M*k) * dof(u(i)), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(im * k * exp2 * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(-(im * M * k) * dof(u(i)), tf(u(i)), sigma(i, j), gauss); // interface unknowns - formulation(i,j).integral(- 2. * (im * k * exp2) * u(i), tf(g(i,j)), sigma(i,j), gauss); - + formulation(i, j).integral(-2. * (im * k * exp2) * u(i), tf(g(i, j)), sigma(i, j), gauss); + // define Pade coefficients std::vector< std::complex< double > > c(padeOrder, 0.); for(int l = 0; l < padeOrder; ++l) { @@ -328,47 +316,46 @@ int main(int argc, char **argv) } // define the auxiliary fields - std::vector< Field< std::complex< double >, Form::Form0 >* > phi; + std::vector< Field< std::complex< double >, Form::Form0 > * > phi; for(int l = 0; l < padeOrder; ++l) { - phi.push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(l), sigma(i,j), FunctionSpaceTypeForm0::HierarchicalH1, FEMorder)); + phi.push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(l), sigma(i, j), FunctionSpaceTypeForm0::HierarchicalH1, FEMorder)); fieldBucket.push_back(phi.back()); } - + // write the augmented weak form for(int l = 0; l < padeOrder; ++l) { // boundary integral terms relating the auxiliary fields - formulation(i).integral(im * k * exp2 * coef * c[l] * dof(*phi[l]), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(im * k * exp2 * coef * c[l] * dof(u(i)), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(im * k * exp2 * coef * c[l] * dof(*phi[l]), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(im * k * exp2 * coef * c[l] * dof(u(i)), tf(u(i)), sigma(i, j), gauss); // coupling of the auxiliary equations - formulation(i).integral(grad(dof(*phi[l])), grad(tf(*phi[l])), sigma(i,j), gauss); - formulation(i).integral(- (k*k)/ beta2 * (exp1 * c[l] + 1.) * dof(*phi[l]), tf(*phi[l]), sigma(i,j), gauss); - formulation(i).integral(- (k*k)/ beta2 * exp1 * (c[l] + 1.) * dof(u(i)), tf(*phi[l]), sigma(i,j), gauss); + formulation(i).integral(grad(dof(*phi[l])), grad(tf(*phi[l])), sigma(i, j), gauss); + formulation(i).integral(-(k * k) / beta2 * (exp1 * c[l] + 1.) * dof(*phi[l]), tf(*phi[l]), sigma(i, j), gauss); + formulation(i).integral(-(k * k) / beta2 * exp1 * (c[l] + 1.) * dof(u(i)), tf(*phi[l]), sigma(i, j), gauss); } - + // interfaces unknowns for(int l = 0; l < padeOrder; ++l) { - formulation(i,j).integral(- 2. * (im * k * exp2) * coef * c[l] * (*phi[l]), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(- 2. * (im * k * exp2) * coef * c[l] * u(i), tf(g(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(-2. * (im * k * exp2) * coef * c[l] * (*phi[l]), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(-2. * (im * k * exp2) * coef * c[l] * u(i), tf(g(i, j)), sigma(i, j), gauss); } } } - else { + else { msg::error << "Transmission condition not available !" << msg::endl; exit(0); } - } // Solve the DDM formulation formulation.pre(); - formulation.solve("jacobi",1e-6); + formulation.solve("jacobi", 1e-6); //formulation.solve("gmres"); double num = 0., den = 0.; double local_err; // L2 local error for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) { - if ( i % MPI_Size == MPI_Rank ) { - save(+u(i),omega(i),"u_"+std::to_string(i)); + if(gmshddm::mpi::isItMySubdomain(i)) { + save(+u(i), omega(i), "u_" + std::to_string(i)); save(*solution, omega(i), "uex_" + std::to_string(i)); save(*solution - u(i), omega(i), "err_" + std::to_string(i)); @@ -390,7 +377,7 @@ int main(int argc, char **argv) } */ - msg::info << "Local L2 error = " << local_err << " on rank " << MPI_Rank <<msg::endl; + msg::info << "Local L2 error = " << local_err << " on rank " << gmshddm::mpi::getMPIRank() << msg::endl; for(unsigned int i = 0; i < fieldBucket.size(); ++i) { delete fieldBucket[i]; diff --git a/demos/maxwell/CMakeLists.txt b/examples/maxwell/CMakeLists.txt similarity index 74% rename from demos/maxwell/CMakeLists.txt rename to examples/maxwell/CMakeLists.txt index c78119517ddd374d9663ff309fd250afbcd67ab9..14a0aa41dcb8ea385d1342e3e16ff57ca12f1ea2 100644 --- a/demos/maxwell/CMakeLists.txt +++ b/examples/maxwell/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.0 FATAL_ERROR) project(demo CXX) -include(${CMAKE_CURRENT_SOURCE_DIR}/../demos.cmake) +include(${CMAKE_CURRENT_SOURCE_DIR}/../example.cmake) add_executable(demo main.cpp ${EXTRA_INCS}) target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/demos/maxwell/main.cpp b/examples/maxwell/main.cpp similarity index 70% rename from demos/maxwell/main.cpp rename to examples/maxwell/main.cpp index 7bfc36611a4a49a5dcd07d061949d1e1b2b2ee65..9b3e8c39a5c056b8942cea491a8f8bd5073dc43d 100644 --- a/demos/maxwell/main.cpp +++ b/examples/maxwell/main.cpp @@ -1,10 +1,9 @@ +#include "gmsh.h" + #include <gmshddm/Formulation.h> #include <gmshddm/GmshDdm.h> - #include <gmshfem/Message.h> -#include "gmsh.h" - using namespace gmshddm; using namespace gmshddm::common; using namespace gmshddm::domain; @@ -20,8 +19,9 @@ using namespace gmshfem::function; using namespace gmshfem::post; void waveguide(const unsigned int nDom, const double a, const double b, - const double L, const double lc) { - if (nDom < 2) { + const double L, const double lc) +{ + if(nDom < 2) { gmshfem::msg::error << "nDom should be at least equal to 2." << gmshfem::msg::endl; exit(0); @@ -30,9 +30,9 @@ void waveguide(const unsigned int nDom, const double a, const double b, const double domainSize = L / nDom; - std::vector<std::vector<int>> borders(nDom); - std::vector<int> omegas(nDom); - std::vector<int> sigmas(nDom - 1); + std::vector< std::vector< int > > borders(nDom); + std::vector< int > omegas(nDom); + std::vector< int > sigmas(nDom - 1); int p[4]; @@ -48,7 +48,7 @@ void waveguide(const unsigned int nDom, const double a, const double b, line[3] = gmsh::model::geo::addLine(p[3], p[0]); int curveLoop = - gmsh::model::geo::addCurveLoop({line[3], line[0], -line[1], line[2]}); + gmsh::model::geo::addCurveLoop({line[3], line[0], -line[1], line[2]}); int tagSurface = gmsh::model::geo::addPlaneSurface({curveLoop}); gmsh::model::geo::mesh::setTransfiniteSurface(tagSurface); @@ -56,11 +56,11 @@ void waveguide(const unsigned int nDom, const double a, const double b, gmsh::model::addPhysicalGroup(2, {tagSurface}, 1); gmsh::model::setPhysicalName(2, 1, "gammaDir"); - for (unsigned int i = 0; i < nDom; ++i) { - std::vector<std::pair<int, int>> extrudeEntities; + for(unsigned int i = 0; i < nDom; ++i) { + std::vector< std::pair< int, int > > extrudeEntities; gmsh::model::geo::extrude( - {std::make_pair(2, tagSurface)}, domainSize, 0., 0., extrudeEntities, - {static_cast<int>((domainSize) / lc)}, std::vector<double>(), true); + {std::make_pair(2, tagSurface)}, domainSize, 0., 0., extrudeEntities, + {static_cast< int >((domainSize) / lc)}, std::vector< double >(), true); tagSurface = extrudeEntities[0].second; omegas[i] = extrudeEntities[1].second; borders[i].push_back(extrudeEntities[2].second); @@ -73,26 +73,27 @@ void waveguide(const unsigned int nDom, const double a, const double b, gmsh::model::addPhysicalGroup(2, {tagSurface}, 2); gmsh::model::setPhysicalName(2, 2, "gammaInf"); - for (unsigned int i = 0; i < borders.size(); ++i) { + for(unsigned int i = 0; i < borders.size(); ++i) { gmsh::model::addPhysicalGroup(2, borders[i], 100 + i); gmsh::model::setPhysicalName(2, 100 + i, "border_" + std::to_string(i)); } - for (unsigned int i = 0; i < omegas.size(); ++i) { + for(unsigned int i = 0; i < omegas.size(); ++i) { gmsh::model::addPhysicalGroup(3, {omegas[i]}, i + 1); gmsh::model::setPhysicalName(3, i + 1, "omega_" + std::to_string(i)); } - for (unsigned int i = 0; i < sigmas.size(); ++i) { + for(unsigned int i = 0; i < sigmas.size(); ++i) { gmsh::model::addPhysicalGroup(2, {sigmas[i]}, 200 + i); gmsh::model::setPhysicalName( - 2, 200 + i, "sigma_" + std::to_string(i) + "_" + std::to_string(i + 1)); + 2, 200 + i, "sigma_" + std::to_string(i) + "_" + std::to_string(i + 1)); } gmsh::model::geo::synchronize(); gmsh::model::mesh::generate(); } -int main(int argc, char **argv) { +int main(int argc, char **argv) +{ GmshDdm gmshDdm(argc, argv); /**** * input Data @@ -121,34 +122,35 @@ int main(int argc, char **argv) { double ky = mode_m * pi / a; double kz = mode_n * pi / b; double kc = sqrt(ky * ky + kz * kz); - std::complex<double> im(0., 1.); - std::complex<double> beta; - if (-kc * kc + k * k >= 0) { + std::complex< double > im(0., 1.); + std::complex< double > beta; + if(-kc * kc + k * k >= 0) { beta = sqrt(-kc * kc + k * k); - } else { + } + else { beta = -im * sqrt(kc * kc - k * k); } // TM - VectorFunction<std::complex<double>> Einc = vector<std::complex<double>>( - sin<std::complex<double>>(ky * y<std::complex<double>>()) * - sin<std::complex<double>>(kz * z<std::complex<double>>()), - im * beta * ky / (kc * kc) * - cos<std::complex<double>>(ky * y<std::complex<double>>()) * - sin<std::complex<double>>(kz * z<std::complex<double>>()), - im * beta * kz / (kc * kc) * - cos<std::complex<double>>(kz * z<std::complex<double>>()) * - sin<std::complex<double>>(ky * y<std::complex<double>>())); - /**** + VectorFunction< std::complex< double > > Einc = vector< std::complex< double > >( + sin< std::complex< double > >(ky * y< std::complex< double > >()) * + sin< std::complex< double > >(kz * z< std::complex< double > >()), + im * beta * ky / (kc * kc) * + cos< std::complex< double > >(ky * y< std::complex< double > >()) * + sin< std::complex< double > >(kz * z< std::complex< double > >()), + im * beta * kz / (kc * kc) * + cos< std::complex< double > >(kz * z< std::complex< double > >()) * + sin< std::complex< double > >(ky * y< std::complex< double > >())); + /**** * FEM ****/ - if (!ddm) { - /**** + if(!ddm) { + /**** * Gmsh part ****/ gmsh::model::add("waveguide"); const double domainSize = L; - std::vector<std::vector<int>> borders(1); - std::vector<int> omegas(1); + std::vector< std::vector< int > > borders(1); + std::vector< int > omegas(1); int p[4]; p[0] = gmsh::model::geo::addPoint(0., 0, 0., lc); p[1] = gmsh::model::geo::addPoint(0, a, 0., lc); @@ -162,17 +164,17 @@ int main(int argc, char **argv) { line[3] = gmsh::model::geo::addLine(p[3], p[0]); int curveLoop = - gmsh::model::geo::addCurveLoop({line[3], line[0], -line[1], line[2]}); + gmsh::model::geo::addCurveLoop({line[3], line[0], -line[1], line[2]}); int tagSurface = gmsh::model::geo::addPlaneSurface({curveLoop}); gmsh::model::geo::mesh::setTransfiniteSurface(tagSurface); gmsh::model::geo::mesh::setRecombine(2, tagSurface); gmsh::model::addPhysicalGroup(2, {tagSurface}, 1); gmsh::model::setPhysicalName(2, 1, "gammaDir"); - std::vector<std::pair<int, int>> extrudeEntities; + std::vector< std::pair< int, int > > extrudeEntities; gmsh::model::geo::extrude( - {std::make_pair(2, tagSurface)}, domainSize, 0., 0., extrudeEntities, - {static_cast<int>((domainSize) / lc)}, std::vector<double>(), true); + {std::make_pair(2, tagSurface)}, domainSize, 0., 0., extrudeEntities, + {static_cast< int >((domainSize) / lc)}, std::vector< double >(), true); tagSurface = extrudeEntities[0].second; omegas[0] = extrudeEntities[1].second; borders[0].push_back(extrudeEntities[2].second); @@ -193,22 +195,22 @@ int main(int argc, char **argv) { Domain gammaInf(2, 2); Domain gammaDir(2, 1); Domain border(2, 100); - gmshfem::problem::Formulation<std::complex<double>> formulation("maxwell"); + gmshfem::problem::Formulation< std::complex< double > > formulation("maxwell"); /**** * Create fields ****/ - Field<std::complex<double>, Form::Form1> E( - "E", omega | gammaDir | gammaInf | border, - FunctionSpaceTypeForm1::HierarchicalHCurl, FEorder); - Field<std::complex<double>, Form::Form1> lambda( - "lambda", gammaDir , FunctionSpaceTypeForm1::HierarchicalHCurl, - FEorder); + Field< std::complex< double >, Form::Form1 > E( + "E", omega | gammaDir | gammaInf | border, + FunctionSpaceTypeForm1::HierarchicalHCurl, FEorder); + Field< std::complex< double >, Form::Form1 > lambda( + "lambda", gammaDir, FunctionSpaceTypeForm1::HierarchicalHCurl, + FEorder); // VectorFunction<std::complex<double>> n = normal<std::complex<double>>(); // VOLUME TERMS - E.addConstraint(border, vector<std::complex<double>>(0.,0.,0.)); - lambda.addConstraint(border, vector<std::complex<double>>(0.,0.,0.)); + E.addConstraint(border, vector< std::complex< double > >(0., 0., 0.)); + lambda.addConstraint(border, vector< std::complex< double > >(0., 0., 0.)); formulation.integral(curl(dof(E)), curl(tf(E)), omega, gauss); formulation.integral(-k * k * dof(E), tf(E), omega, gauss); @@ -217,7 +219,7 @@ int main(int argc, char **argv) { // gauss); formulation.integral(-im * k * dof(E), tf(E), gammaInf, gauss); // Physical source using Lagrange mutiplier - formulation.integral(dof(lambda), tf(E), gammaDir , gauss); + formulation.integral(dof(lambda), tf(E), gammaDir, gauss); formulation.integral(dof(E), tf(lambda), gammaDir, gauss); formulation.integral(Einc, tf(lambda), gammaDir, gauss); formulation.pre(); @@ -242,69 +244,69 @@ int main(int argc, char **argv) { Subdomain gammaDir(nDom); Subdomain border(nDom); Interface sigma(nDom); - for (unsigned int i = 0; i < unsigned(nDom); ++i) { + for(unsigned int i = 0; i < unsigned(nDom); ++i) { omega(i) = Domain(3, i + 1); - if (i == unsigned(nDom) - 1) { + if(i == unsigned(nDom) - 1) { gammaInf(i) = Domain(2, 2); } - if (i == 0) { + if(i == 0) { gammaDir(i) = Domain(2, 1); } border(i) = Domain(2, 100 + i); - if (i != 0) { + if(i != 0) { sigma(i, i - 1) = Domain(2, 200 + i - 1); } - if (i != unsigned(nDom) - 1) { + if(i != unsigned(nDom) - 1) { sigma(i, i + 1) = Domain(2, 200 + i); } } /**** * Define topology ****/ - std::vector<std::vector<unsigned int>> topology(nDom); - for (unsigned int i = 0; i < unsigned(nDom); ++i) { - if (i != 0) { + std::vector< std::vector< unsigned int > > topology(nDom); + for(unsigned int i = 0; i < unsigned(nDom); ++i) { + if(i != 0) { topology[i].push_back(i - 1); } - if (i != unsigned(nDom) - 1) { + if(i != unsigned(nDom) - 1) { topology[i].push_back(i + 1); } } // Create DDM formulation - gmshddm::problem::Formulation<std::complex<double>> formulation("Maxwell", - topology); + gmshddm::problem::Formulation< std::complex< double > > formulation("Maxwell", + topology); // Create fields - SubdomainField<std::complex<double>, Form::Form1> E( - "E", omega | gammaDir | gammaInf | sigma, - FunctionSpaceTypeForm1::HierarchicalHCurl, FEorder); - SubdomainField<std::complex<double>, Form::Form1> lambda( - "lambda", gammaDir, FunctionSpaceTypeForm1::HierarchicalHCurl, - FEorder); - InterfaceField<std::complex<double>, Form::Form1> g( - "g", sigma, FunctionSpaceTypeForm1::HierarchicalHCurl, 0); + SubdomainField< std::complex< double >, Form::Form1 > E( + "E", omega | gammaDir | gammaInf | sigma, + FunctionSpaceTypeForm1::HierarchicalHCurl, FEorder); + SubdomainField< std::complex< double >, Form::Form1 > lambda( + "lambda", gammaDir, FunctionSpaceTypeForm1::HierarchicalHCurl, + FEorder); + InterfaceField< std::complex< double >, Form::Form1 > g( + "g", sigma, FunctionSpaceTypeForm1::HierarchicalHCurl, 0); // Tell to the formulation that g is field that have to be exchanged between // subdomains formulation.addInterfaceField(g); - VectorFunction<std::complex<double>> n = normal<std::complex<double>>(); - VectorFunction<std::complex<double>> functionZero = - vector<std::complex<double>>(0., 0., 0.); - for (unsigned int i = 0; i < unsigned(nDom); ++i){ + VectorFunction< std::complex< double > > n = normal< std::complex< double > >(); + VectorFunction< std::complex< double > > functionZero = + vector< std::complex< double > >(0., 0., 0.); + for(unsigned int i = 0; i < unsigned(nDom); ++i) { E(i).addConstraint(border(i), functionZero); lambda(i).addConstraint(border(i), functionZero); - for (unsigned int jj = 0; jj < topology[i].size(); ++jj) { + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; g(i, j).addConstraint(border(i), functionZero); } } // Add terms to the formulation - for (unsigned int i = 0; i < unsigned(nDom); ++i) { + for(unsigned int i = 0; i < unsigned(nDom); ++i) { // VOLUME TERMS formulation(i).integral(curl(dof(E(i))), curl(tf(E(i))), omega(i), gauss); formulation(i).integral(-k * k * dof(E(i)), tf(E(i)), omega(i), gauss); @@ -316,13 +318,13 @@ int main(int argc, char **argv) { // Physical source using Lagrange mutiplier formulation(i).integral(dof(lambda(i)), tf(E(i)), gammaDir(i), gauss); - formulation(i).integral(dof(E(i)), tf(lambda(i)), gammaDir(i) , + formulation(i).integral(dof(E(i)), tf(lambda(i)), gammaDir(i), gauss); formulation(i).integral(formulation.physicalSource(Einc), tf(lambda(i)), gammaDir(i), gauss); // Artificial source on transmission boundaries - for (unsigned int jj = 0; jj < topology[i].size(); ++jj) { + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; /* formulation(i).integral(-im * k * n % (dof(E(i)) % n), tf(E(i)), sigma(i, j), gauss);*/ @@ -332,7 +334,7 @@ int main(int argc, char **argv) { sigma(i, j), gauss); } // INTERFACE TERMS - for (unsigned int jj = 0; jj < topology[i].size(); ++jj) { + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss); @@ -348,7 +350,7 @@ int main(int argc, char **argv) { // Solve the DDM formulation formulation.solve("gmres"); //save field E - for (int i = 0; i < nDom; ++i) { + for(int i = 0; i < nDom; ++i) { save(E(i), omega(i), "E" + std::to_string(i)); } } diff --git a/demos/helmholtzflow/waveguide/CMakeLists.txt b/examples/navier/scattering/CMakeLists.txt old mode 100644 new mode 100755 similarity index 73% rename from demos/helmholtzflow/waveguide/CMakeLists.txt rename to examples/navier/scattering/CMakeLists.txt index 4d2e4b0c91d2a7e0da5a224e581c56eca24b2f71..c368ce59f63aed4e569c1b3e727df6dd0b1d6a45 --- a/demos/helmholtzflow/waveguide/CMakeLists.txt +++ b/examples/navier/scattering/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.0 FATAL_ERROR) project(demo CXX) -include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake) +include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake) add_executable(demo main.cpp ${EXTRA_INCS}) target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/demos/navier/scattering/main.cpp b/examples/navier/scattering/main.cpp similarity index 68% rename from demos/navier/scattering/main.cpp rename to examples/navier/scattering/main.cpp index d94b532a9331cc24d5eda9c7b2611dbf542a3b70..b717993b57c63bbc36ccf4636bcbe6f59d1cee42 100644 --- a/demos/navier/scattering/main.cpp +++ b/examples/navier/scattering/main.cpp @@ -1,10 +1,9 @@ -#include <gmshddm/GmshDdm.h> -#include <gmshddm/Formulation.h> +#include "gmsh.h" +#include <gmshddm/Formulation.h> +#include <gmshddm/GmshDdm.h> #include <gmshfem/Message.h> -#include "gmsh.h" - using namespace gmshddm; using namespace gmshddm::common; using namespace gmshddm::domain; @@ -39,7 +38,7 @@ void circlesConcentric(const unsigned int nDom, const double r0, const double rI std::vector< int > tagsSurface(nDom); for(unsigned int i = 0; i < nDom; ++i) { - tagsSurface[i] = gmsh::model::occ::addPlaneSurface({tagsCurve[i], tagsCurve[i+1]}); + tagsSurface[i] = gmsh::model::occ::addPlaneSurface({tagsCurve[i], tagsCurve[i + 1]}); } gmsh::model::occ::synchronize(); @@ -55,7 +54,7 @@ void circlesConcentric(const unsigned int nDom, const double r0, const double rI gmsh::model::setPhysicalName(2, 1 + i, "omega_" + std::to_string(i)); if(i != nDom - 1) { - gmsh::model::addPhysicalGroup(1, {tagsCircle[i+1]}, 200 + i); + gmsh::model::addPhysicalGroup(1, {tagsCircle[i + 1]}, 200 + i); gmsh::model::setPhysicalName(1, 200 + i, "sigma_" + std::to_string(i)); } } @@ -127,20 +126,20 @@ int main(int argc, char **argv) gmshddm::problem::Formulation< std::complex< double > > formulation("Navier", topology); const double pi = 3.14159265358979323846264338327950288; - const double rho = 1.; //.............[kg/m3] masse volumique - double f = 1.; //.............[Hz] fréquence - double w = 2 * pi * f; //.............[rad/s] pulsation - const double lambda = 1.; //.............[Pa] premier coefficient de Lamé - const double mu = 1.; //.............[Pa] deuxième coefficient de Lamé - const double cP = sqrt((lambda+2.0*mu)/rho); //.............[m/s] celerite de l'onde P - const double cS = sqrt(mu/rho); //.............[m/s] celerite de l'onde S - double kP = w/cP; //.............[1/m] nombre d'onde p - double kS = w/cS; //.............[1/m] nombre d'onde s - - TensorFunction< std::complex< double > > C_xx = tensor< std::complex< double > >((lambda+2.*mu), 0., 0., 0., mu, 0., 0., 0., 0.); + const double rho = 1.; //.............[kg/m3] masse volumique + double f = 1.; //.............[Hz] fréquence + double w = 2 * pi * f; //.............[rad/s] pulsation + const double lambda = 1.; //.............[Pa] premier coefficient de Lamé + const double mu = 1.; //.............[Pa] deuxième coefficient de Lamé + const double cP = sqrt((lambda + 2.0 * mu) / rho); //.............[m/s] celerite de l'onde P + const double cS = sqrt(mu / rho); //.............[m/s] celerite de l'onde S + double kP = w / cP; //.............[1/m] nombre d'onde p + double kS = w / cS; //.............[1/m] nombre d'onde s + + TensorFunction< std::complex< double > > C_xx = tensor< std::complex< double > >((lambda + 2. * mu), 0., 0., 0., mu, 0., 0., 0., 0.); TensorFunction< std::complex< double > > C_xy = tensor< std::complex< double > >(0., lambda, 0., mu, 0., 0., 0., 0., 0.); TensorFunction< std::complex< double > > C_yx = tensor< std::complex< double > >(0., mu, 0., lambda, 0., 0., 0., 0., 0.); - TensorFunction< std::complex< double > > C_yy = tensor< std::complex< double > >(mu, 0., 0., 0., (lambda+2.*mu), 0., 0., 0., 0.); + TensorFunction< std::complex< double > > C_yy = tensor< std::complex< double > >(mu, 0., 0., 0., (lambda + 2. * mu), 0., 0., 0., 0.); // Create fields SubdomainField< std::complex< double >, Form::Form0 > ux("ux", omega | gammaScat | gammaInf | sigma, FunctionSpaceTypeForm0::Lagrange); @@ -172,7 +171,7 @@ int main(int argc, char **argv) // Physical source formulation(i).integral(dof(lambdax(i)), tf(ux(i)), gammaScat(i), gauss); formulation(i).integral(dof(ux(i)), tf(lambdax(i)), gammaScat(i), gauss); - formulation(i).integral(formulation.physicalSource(- (cos(kP * x< std::complex< double > >()) + im * sin(kP * x< std::complex< double > >()))), tf(lambdax(i)), gammaScat(i), gauss); + formulation(i).integral(formulation.physicalSource(-(cos(kP * x< std::complex< double > >()) + im * sin(kP * x< std::complex< double > >()))), tf(lambdax(i)), gammaScat(i), gauss); formulation(i).integral(dof(lambday(i)), tf(uy(i)), gammaScat(i), gauss); formulation(i).integral(dof(uy(i)), tf(lambday(i)), gammaScat(i), gauss); @@ -185,45 +184,45 @@ int main(int argc, char **argv) formulation(i).integral(-im * w * cP * rho * nY * nY * dof(uy(i)), tf(uy(i)), gammaInf(i), gauss); formulation(i).integral(-im * w * cS * rho * nY * nY * dof(ux(i)), tf(ux(i)), gammaInf(i), gauss); - formulation(i).integral( im * w * cS * rho * nX * nY * dof(uy(i)), tf(ux(i)), gammaInf(i), gauss); + formulation(i).integral(im * w * cS * rho * nX * nY * dof(uy(i)), tf(ux(i)), gammaInf(i), gauss); formulation(i).integral(-im * w * cS * rho * nX * nX * dof(uy(i)), tf(uy(i)), gammaInf(i), gauss); - formulation(i).integral( im * w * cS * rho * nX * nY * dof(ux(i)), tf(uy(i)), gammaInf(i), gauss); + formulation(i).integral(im * w * cS * rho * nX * nY * dof(ux(i)), tf(uy(i)), gammaInf(i), gauss); // Artificial source for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(-im * w * cP * rho * nX * nX * dof(ux(i)), tf(ux(i)), sigma(i,j), gauss); - formulation(i).integral(-im * w * cP * rho * nX * nY * dof(uy(i)), tf(ux(i)), sigma(i,j), gauss); - formulation(i).integral(-im * w * cP * rho * nX * nY * dof(ux(i)), tf(uy(i)), sigma(i,j), gauss); - formulation(i).integral(-im * w * cP * rho * nY * nY * dof(uy(i)), tf(uy(i)), sigma(i,j), gauss); - - formulation(i).integral(-im * w * cS * rho * nY * nY * dof(ux(i)), tf(ux(i)), sigma(i,j), gauss); - formulation(i).integral( im * w * cS * rho * nX * nY * dof(uy(i)), tf(ux(i)), sigma(i,j), gauss); - formulation(i).integral(-im * w * cS * rho * nX * nX * dof(uy(i)), tf(uy(i)), sigma(i,j), gauss); - formulation(i).integral( im * w * cS * rho * nX * nY * dof(ux(i)), tf(uy(i)), sigma(i,j), gauss); - - formulation(i).integral(formulation.artificialSource( -gx(j,i) ), tf(ux(i)), sigma(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -gy(j,i) ), tf(uy(i)), sigma(i,j), gauss); + formulation(i).integral(-im * w * cP * rho * nX * nX * dof(ux(i)), tf(ux(i)), sigma(i, j), gauss); + formulation(i).integral(-im * w * cP * rho * nX * nY * dof(uy(i)), tf(ux(i)), sigma(i, j), gauss); + formulation(i).integral(-im * w * cP * rho * nX * nY * dof(ux(i)), tf(uy(i)), sigma(i, j), gauss); + formulation(i).integral(-im * w * cP * rho * nY * nY * dof(uy(i)), tf(uy(i)), sigma(i, j), gauss); + + formulation(i).integral(-im * w * cS * rho * nY * nY * dof(ux(i)), tf(ux(i)), sigma(i, j), gauss); + formulation(i).integral(im * w * cS * rho * nX * nY * dof(uy(i)), tf(ux(i)), sigma(i, j), gauss); + formulation(i).integral(-im * w * cS * rho * nX * nX * dof(uy(i)), tf(uy(i)), sigma(i, j), gauss); + formulation(i).integral(im * w * cS * rho * nX * nY * dof(ux(i)), tf(uy(i)), sigma(i, j), gauss); + + formulation(i).integral(formulation.artificialSource(-gx(j, i)), tf(ux(i)), sigma(i, j), gauss); + formulation(i).integral(formulation.artificialSource(-gy(j, i)), tf(uy(i)), sigma(i, j), gauss); } // INTERFACE TERMS for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i,j).integral(dof(gx(i,j)), tf(gx(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gx(j,i) ), tf(gx(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(dof(gx(i, j)), tf(gx(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(gx(j, i)), tf(gx(i, j)), sigma(i, j), gauss); - formulation(i,j).integral(dof(gy(i,j)), tf(gy(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( gy(j,i) ), tf(gy(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(dof(gy(i, j)), tf(gy(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(gy(j, i)), tf(gy(i, j)), sigma(i, j), gauss); - formulation(i,j).integral(2. * im * w * cP * rho * nX * nX * ux(i), tf(gx(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(2. * im * w * cP * rho * nX * nY * uy(i), tf(gx(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(2. * im * w * cP * rho * nX * nY * ux(i), tf(gy(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(2. * im * w * cP * rho * nY * nY * uy(i), tf(gy(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(2. * im * w * cP * rho * nX * nX * ux(i), tf(gx(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(2. * im * w * cP * rho * nX * nY * uy(i), tf(gx(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(2. * im * w * cP * rho * nX * nY * ux(i), tf(gy(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(2. * im * w * cP * rho * nY * nY * uy(i), tf(gy(i, j)), sigma(i, j), gauss); - formulation(i,j).integral(2. * im * w * cS * rho * nY * nY * ux(i), tf(gx(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(-2. * im * w * cS * rho * nX * nY * uy(i), tf(gx(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(2. * im * w * cS * rho * nX * nX * uy(i), tf(gy(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(-2. * im * w * cS * rho * nX * nY * ux(i), tf(gy(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(2. * im * w * cS * rho * nY * nY * ux(i), tf(gx(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(-2. * im * w * cS * rho * nX * nY * uy(i), tf(gx(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(2. * im * w * cS * rho * nX * nX * uy(i), tf(gy(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(-2. * im * w * cS * rho * nX * nY * ux(i), tf(gy(i, j)), sigma(i, j), gauss); } } diff --git a/tutorials/.clang-format b/tutorials/.clang-format new file mode 100644 index 0000000000000000000000000000000000000000..d902ab979dca0361bb1b063642225163bda4eb7d --- /dev/null +++ b/tutorials/.clang-format @@ -0,0 +1,112 @@ +--- +Language: Cpp +AccessModifierOffset: -1 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +AlignEscapedNewlines: Right +AlignOperands: true +AlignTrailingComments: false +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: false +AllowShortCaseLabelsOnASingleLine: true +AllowShortFunctionsOnASingleLine: Empty +AllowShortIfStatementsOnASingleLine: true +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: true +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: true + AfterControlStatement: false + AfterEnum: false + AfterFunction: true + AfterNamespace: true + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + AfterExternBlock: false + BeforeCatch: true + BeforeElse: true + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Custom +BreakBeforeInheritanceComma: false +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: AfterColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 0 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 2 +ContinuationIndentWidth: 2 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IncludeBlocks: Regroup +IncludeCategories: + - Regex: '^"(?!gmsh\.h"|petsc\.h"|Eigen")' + Priority: 2 + - Regex: '"gmsh.h"|"petsc.h"|"Eigen"' + Priority: 3 + - Regex: '^<' + Priority: 4 + - Regex: '.*' + Priority: 1 +IncludeIsMainRegex: '$' +IndentCaseLabels: false +IndentPPDirectives: None +IndentWidth: 2 +IndentWrappedFunctionNames: false +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 2 +NamespaceIndentation: All +ObjCBlockIndentWidth: 2 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 60 +PointerAlignment: Right +ReflowComments: false +SortIncludes: true +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterTemplateKeyword: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: Never +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: true +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Cpp11 +TabWidth: 4 +UseTab: Never +... + diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..210187a9649003866f9d32cff8264416605deb81 --- /dev/null +++ b/tutorials/CMakeLists.txt @@ -0,0 +1,11 @@ +set(TUTOS + helmholtz/waveguide + maxwell/waveguide +) + +foreach(TUTO ${TUTOS}) + add_test(NAME "tuto:${TUTO}:mkdir" COMMAND "${CMAKE_COMMAND}" -E make_directory build WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/${TUTO}) + add_test(NAME "tuto:${TUTO}:cmake" COMMAND "${CMAKE_COMMAND}" -DCMAKE_PREFIX_PATH=${CMAKE_INSTALL_PREFIX} .. WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/${TUTO}/build) + add_test(NAME "tuto:${TUTO}:build" COMMAND "${CMAKE_MAKE_PROGRAM}" WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/${TUTO}/build) + add_test(NAME "tuto:${TUTO}:run" COMMAND ./tuto WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/${TUTO}/build) +endforeach() diff --git a/tutorials/README.md b/tutorials/README.md new file mode 100644 index 0000000000000000000000000000000000000000..fc6fff349f9a001d05bf1af2d7fe1e7ec88bf1c9 --- /dev/null +++ b/tutorials/README.md @@ -0,0 +1,3 @@ +# GmshDDM tutorials + +Todo: Create tutorials based on problems of increasing difficulties that shows first the basic features of GmshDDM ('basic tutorial'), then that presents some more advanced uses ('advanced tutorial') @CGeuzaine ? diff --git a/tutorials/helmholtz/waveguide/CMakeLists.txt b/tutorials/helmholtz/waveguide/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..7481f3ad58dcd475585824c7f40fddf332e191dc --- /dev/null +++ b/tutorials/helmholtz/waveguide/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(tuto CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../tutorial.cmake) + +add_executable(tuto main.cpp ${EXTRA_INCS}) +target_link_libraries(tuto ${EXTRA_LIBS}) diff --git a/demos/helmholtz/waveguide/main.cpp b/tutorials/helmholtz/waveguide/main.cpp similarity index 78% rename from demos/helmholtz/waveguide/main.cpp rename to tutorials/helmholtz/waveguide/main.cpp index a19176604364f6c8aa79f9e501051c45ef92468f..0a9b5b028798234adc46d604fa99b4669227eee2 100644 --- a/demos/helmholtz/waveguide/main.cpp +++ b/tutorials/helmholtz/waveguide/main.cpp @@ -1,11 +1,10 @@ -#include <gmshddm/GmshDdm.h> +#include "gmsh.h" + #include <gmshddm/Formulation.h> +#include <gmshddm/GmshDdm.h> #include <gmshddm/MPIInterface.h> - #include <gmshfem/Message.h> -#include "gmsh.h" - using namespace gmshddm; using namespace gmshddm::common; using namespace gmshddm::domain; @@ -44,6 +43,7 @@ void waveguide(const unsigned int nDom, const double l, const double L, const do gmsh::model::addPhysicalGroup(1, {line}, 1); gmsh::model::setPhysicalName(1, 1, "gammaDir"); + std::vector< int > cornerTags; for(unsigned int i = 0; i < nDom; ++i) { std::vector< std::pair< int, int > > extrudeEntities; gmsh::model::geo::extrude({std::make_pair(1, line)}, domainSize, 0., 0., extrudeEntities, {static_cast< int >((domainSize) / lc)}, std::vector< double >(), true); @@ -52,6 +52,10 @@ void waveguide(const unsigned int nDom, const double l, const double L, const do borders[i].push_back(extrudeEntities[2].second); borders[i].push_back(extrudeEntities[3].second); sigmas[i] = line; + if(i != nDom - 1) { + cornerTags.push_back(3 + cornerTags.size()); + cornerTags.push_back(3 + cornerTags.size()); + } } gmsh::model::addPhysicalGroup(1, {line}, 2); @@ -72,6 +76,9 @@ void waveguide(const unsigned int nDom, const double l, const double L, const do gmsh::model::setPhysicalName(1, 200 + i, "sigma_" + std::to_string(i) + "_" + std::to_string(i + 1)); } + gmsh::model::addPhysicalGroup(0, cornerTags, 1000); + gmsh::model::setPhysicalName(0, 1000, "corner"); + gmsh::model::geo::synchronize(); gmsh::model::mesh::generate(); } @@ -79,10 +86,7 @@ void waveguide(const unsigned int nDom, const double l, const double L, const do int main(int argc, char **argv) { GmshDdm gmshDdm(argc, argv); - - const unsigned int MPI_Size = getMPISize(); - const unsigned int MPI_Rank = getMPIRank(); - + double pi = 3.14159265359; int nDom = 3; @@ -111,6 +115,7 @@ int main(int argc, char **argv) Subdomain gammaDir(nDom); Subdomain border(nDom); Interface sigma(nDom); + Domain corner = Domain("corner"); for(unsigned int i = 0; i < nDom; ++i) { omega(i) = Domain(2, i + 1); @@ -150,8 +155,17 @@ int main(int argc, char **argv) // Create fields SubdomainField< std::complex< double >, Form::Form0 > u("u", omega | gammaDir | gammaInf | border | sigma, FunctionSpaceTypeForm0::HierarchicalH1, order); - SubdomainField< std::complex< double >, Form::Form0 > lambda("lambda", gammaDir | border, FunctionSpaceTypeForm0::HierarchicalH1, order); + for(unsigned int i = 0; i < nDom; ++i) { + u(i).addConstraint(border(i), 0.); + } + u(0).addConstraint(gammaDir(0), sin< std::complex< double > >(pi * y< std::complex< double > >() / l)); InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, order); + for(unsigned int i = 0; i < nDom; ++i) { + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { + const unsigned int j = topology[i][jj]; + g(i, j).addConstraint(corner, 0.); + } + } std::complex< double > im(0., 1.); @@ -161,29 +175,23 @@ int main(int argc, char **argv) // Add terms to the formulation for(unsigned int i = 0; i < nDom; ++i) { // VOLUME TERMS - formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss); - formulation(i).integral(- k * k * dof(u(i)), tf(u(i)), omega(i), gauss); - formulation(i).integral(- im * k * dof(u(i)), tf(u(i)), gammaInf(i), gauss); - - // Physical source - formulation(i).integral(dof(lambda(i)), tf(u(i)), gammaDir(i) | border(i), gauss); - formulation(i).integral(dof(u(i)), tf(lambda(i)), gammaDir(i) | border(i), gauss); - formulation(i).integral(formulation.physicalSource( -sin< std::complex< double > >((kIn * pi / l) * y< std::complex< double > >() ) ), tf(lambda(i)), gammaDir(i), gauss); + formulation(i).integral(-k * k * dof(u(i)), tf(u(i)), omega(i), gauss); + formulation(i).integral(-im * k * dof(u(i)), tf(u(i)), gammaInf(i), gauss); // Artificial source for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i).integral(- im * k * dof(u(i)), tf(u(i)), sigma(i,j), gauss); - formulation(i).integral(formulation.artificialSource( -g(j,i) ), tf(u(i)), sigma(i,j), gauss); + formulation(i).integral(-im * k * dof(u(i)), tf(u(i)), sigma(i, j), gauss); + formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss); } // INTERFACE TERMS for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { const unsigned int j = topology[i][jj]; - formulation(i,j).integral(dof(g(i,j)), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(formulation.artificialSource( g(j,i) ), tf(g(i,j)), sigma(i,j), gauss); - formulation(i,j).integral(2. * im * k * u(i), tf(g(i,j)), sigma(i,j), gauss); + formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(g(j, i)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(2. * im * k * u(i), tf(g(i, j)), sigma(i, j), gauss); } } @@ -194,8 +202,8 @@ int main(int argc, char **argv) formulation.solve("gmres", 1e-6, iterMax); for(unsigned int i = 0; i < nDom; ++i) { - if(i % MPI_Size == MPI_Rank) { - save(u(i), omega(i), "u_" + std::to_string(i)); + if(isItMySubdomain(i)) { + save(u(i), omega(i), "u_" + std::to_string(i)); } } diff --git a/tutorials/maxwell/main.cpp b/tutorials/maxwell/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..9b3e8c39a5c056b8942cea491a8f8bd5073dc43d --- /dev/null +++ b/tutorials/maxwell/main.cpp @@ -0,0 +1,359 @@ +#include "gmsh.h" + +#include <gmshddm/Formulation.h> +#include <gmshddm/GmshDdm.h> +#include <gmshfem/Message.h> + +using namespace gmshddm; +using namespace gmshddm::common; +using namespace gmshddm::domain; +using namespace gmshddm::problem; +using namespace gmshddm::field; + +using namespace gmshfem; +using namespace gmshfem::domain; +using namespace gmshfem::equation; +using namespace gmshfem::problem; +using namespace gmshfem::field; +using namespace gmshfem::function; +using namespace gmshfem::post; + +void waveguide(const unsigned int nDom, const double a, const double b, + const double L, const double lc) +{ + if(nDom < 2) { + gmshfem::msg::error << "nDom should be at least equal to 2." + << gmshfem::msg::endl; + exit(0); + } + gmsh::model::add("waveguide"); + + const double domainSize = L / nDom; + + std::vector< std::vector< int > > borders(nDom); + std::vector< int > omegas(nDom); + std::vector< int > sigmas(nDom - 1); + + int p[4]; + + p[0] = gmsh::model::geo::addPoint(0., 0, 0., lc); + p[1] = gmsh::model::geo::addPoint(0, a, 0., lc); + p[2] = gmsh::model::geo::addPoint(0, a, b, lc); + p[3] = gmsh::model::geo::addPoint(0, 0., b, lc); + int line[4]; + + line[0] = gmsh::model::geo::addLine(p[0], p[1]); + line[1] = gmsh::model::geo::addLine(p[2], p[1]); + line[2] = gmsh::model::geo::addLine(p[2], p[3]); + line[3] = gmsh::model::geo::addLine(p[3], p[0]); + + int curveLoop = + gmsh::model::geo::addCurveLoop({line[3], line[0], -line[1], line[2]}); + + int tagSurface = gmsh::model::geo::addPlaneSurface({curveLoop}); + gmsh::model::geo::mesh::setTransfiniteSurface(tagSurface); + gmsh::model::geo::mesh::setRecombine(2, tagSurface); + gmsh::model::addPhysicalGroup(2, {tagSurface}, 1); + gmsh::model::setPhysicalName(2, 1, "gammaDir"); + + for(unsigned int i = 0; i < nDom; ++i) { + std::vector< std::pair< int, int > > extrudeEntities; + gmsh::model::geo::extrude( + {std::make_pair(2, tagSurface)}, domainSize, 0., 0., extrudeEntities, + {static_cast< int >((domainSize) / lc)}, std::vector< double >(), true); + tagSurface = extrudeEntities[0].second; + omegas[i] = extrudeEntities[1].second; + borders[i].push_back(extrudeEntities[2].second); + borders[i].push_back(extrudeEntities[3].second); + borders[i].push_back(extrudeEntities[4].second); + borders[i].push_back(extrudeEntities[5].second); + sigmas[i] = tagSurface; + } + + gmsh::model::addPhysicalGroup(2, {tagSurface}, 2); + gmsh::model::setPhysicalName(2, 2, "gammaInf"); + + for(unsigned int i = 0; i < borders.size(); ++i) { + gmsh::model::addPhysicalGroup(2, borders[i], 100 + i); + gmsh::model::setPhysicalName(2, 100 + i, "border_" + std::to_string(i)); + } + + for(unsigned int i = 0; i < omegas.size(); ++i) { + gmsh::model::addPhysicalGroup(3, {omegas[i]}, i + 1); + gmsh::model::setPhysicalName(3, i + 1, "omega_" + std::to_string(i)); + } + + for(unsigned int i = 0; i < sigmas.size(); ++i) { + gmsh::model::addPhysicalGroup(2, {sigmas[i]}, 200 + i); + gmsh::model::setPhysicalName( + 2, 200 + i, "sigma_" + std::to_string(i) + "_" + std::to_string(i + 1)); + } + gmsh::model::geo::synchronize(); + gmsh::model::mesh::generate(); +} + +int main(int argc, char **argv) +{ + GmshDdm gmshDdm(argc, argv); + /**** + * input Data + ****/ + double pi = 3.14159265359; + int nDom = 2; + double L = 2; + double a = 1; + double b = 1; + double k = 10; + double mode_m = 2.; + double mode_n = 1.; + int FEorder = 0; + double lc = 0.5; + bool ddm = true; + std::string gauss = "Gauss4"; + gmshDdm.userDefinedParameter(ddm, "ddm"); + gmshDdm.userDefinedParameter(nDom, "nDom"); + gmshDdm.userDefinedParameter(L, "L"); + gmshDdm.userDefinedParameter(a, "a"); + gmshDdm.userDefinedParameter(b, "b"); + gmshDdm.userDefinedParameter(lc, "lc"); + gmshDdm.userDefinedParameter(k, "k"); + gmshDdm.userDefinedParameter(FEorder, "FEorder"); + gmshDdm.userDefinedParameter(gauss, "gauss"); + double ky = mode_m * pi / a; + double kz = mode_n * pi / b; + double kc = sqrt(ky * ky + kz * kz); + std::complex< double > im(0., 1.); + std::complex< double > beta; + if(-kc * kc + k * k >= 0) { + beta = sqrt(-kc * kc + k * k); + } + else { + beta = -im * sqrt(kc * kc - k * k); + } + // TM + VectorFunction< std::complex< double > > Einc = vector< std::complex< double > >( + sin< std::complex< double > >(ky * y< std::complex< double > >()) * + sin< std::complex< double > >(kz * z< std::complex< double > >()), + im * beta * ky / (kc * kc) * + cos< std::complex< double > >(ky * y< std::complex< double > >()) * + sin< std::complex< double > >(kz * z< std::complex< double > >()), + im * beta * kz / (kc * kc) * + cos< std::complex< double > >(kz * z< std::complex< double > >()) * + sin< std::complex< double > >(ky * y< std::complex< double > >())); + /**** + * FEM + ****/ + if(!ddm) { + /**** + * Gmsh part + ****/ + gmsh::model::add("waveguide"); + const double domainSize = L; + std::vector< std::vector< int > > borders(1); + std::vector< int > omegas(1); + int p[4]; + p[0] = gmsh::model::geo::addPoint(0., 0, 0., lc); + p[1] = gmsh::model::geo::addPoint(0, a, 0., lc); + p[2] = gmsh::model::geo::addPoint(0, a, b, lc); + p[3] = gmsh::model::geo::addPoint(0, 0., b, lc); + int line[4]; + + line[0] = gmsh::model::geo::addLine(p[0], p[1]); + line[1] = gmsh::model::geo::addLine(p[2], p[1]); + line[2] = gmsh::model::geo::addLine(p[2], p[3]); + line[3] = gmsh::model::geo::addLine(p[3], p[0]); + + int curveLoop = + gmsh::model::geo::addCurveLoop({line[3], line[0], -line[1], line[2]}); + + int tagSurface = gmsh::model::geo::addPlaneSurface({curveLoop}); + gmsh::model::geo::mesh::setTransfiniteSurface(tagSurface); + gmsh::model::geo::mesh::setRecombine(2, tagSurface); + gmsh::model::addPhysicalGroup(2, {tagSurface}, 1); + gmsh::model::setPhysicalName(2, 1, "gammaDir"); + std::vector< std::pair< int, int > > extrudeEntities; + gmsh::model::geo::extrude( + {std::make_pair(2, tagSurface)}, domainSize, 0., 0., extrudeEntities, + {static_cast< int >((domainSize) / lc)}, std::vector< double >(), true); + tagSurface = extrudeEntities[0].second; + omegas[0] = extrudeEntities[1].second; + borders[0].push_back(extrudeEntities[2].second); + borders[0].push_back(extrudeEntities[3].second); + borders[0].push_back(extrudeEntities[4].second); + borders[0].push_back(extrudeEntities[5].second); + + gmsh::model::addPhysicalGroup(2, {tagSurface}, 2); + gmsh::model::setPhysicalName(2, 2, "gammaInf"); + gmsh::model::addPhysicalGroup(2, borders[0], 100); + gmsh::model::setPhysicalName(2, 100, "border_" + std::to_string(0)); + gmsh::model::addPhysicalGroup(3, {omegas[0]}, 1); + gmsh::model::setPhysicalName(3, 1, "omega_" + std::to_string(0)); + + gmsh::model::geo::synchronize(); + gmsh::model::mesh::generate(); + Domain omega(3, 1); + Domain gammaInf(2, 2); + Domain gammaDir(2, 1); + Domain border(2, 100); + gmshfem::problem::Formulation< std::complex< double > > formulation("maxwell"); + + /**** + * Create fields + ****/ + Field< std::complex< double >, Form::Form1 > E( + "E", omega | gammaDir | gammaInf | border, + FunctionSpaceTypeForm1::HierarchicalHCurl, FEorder); + Field< std::complex< double >, Form::Form1 > lambda( + "lambda", gammaDir, FunctionSpaceTypeForm1::HierarchicalHCurl, + FEorder); + // VectorFunction<std::complex<double>> n = normal<std::complex<double>>(); + // VOLUME TERMS + + E.addConstraint(border, vector< std::complex< double > >(0., 0., 0.)); + lambda.addConstraint(border, vector< std::complex< double > >(0., 0., 0.)); + + formulation.integral(curl(dof(E)), curl(tf(E)), omega, gauss); + formulation.integral(-k * k * dof(E), tf(E), omega, gauss); + // Silver Muller + // formulation.integral(- im * k * n % (dof(E) % n), tf(E), gammaInf, + // gauss); + formulation.integral(-im * k * dof(E), tf(E), gammaInf, gauss); + // Physical source using Lagrange mutiplier + formulation.integral(dof(lambda), tf(E), gammaDir, gauss); + formulation.integral(dof(E), tf(lambda), gammaDir, gauss); + formulation.integral(Einc, tf(lambda), gammaDir, gauss); + formulation.pre(); + formulation.assemble(); + formulation.solve(); + save(E, omega, "E"); + } + /**** + * DDM + ****/ + else { + /**** + * Build geometry and mesh using gmsh API + ****/ + waveguide(nDom, a, b, L, lc); + + /**** + * Define domains + ****/ + Subdomain omega(nDom); + Subdomain gammaInf(nDom); + Subdomain gammaDir(nDom); + Subdomain border(nDom); + Interface sigma(nDom); + for(unsigned int i = 0; i < unsigned(nDom); ++i) { + omega(i) = Domain(3, i + 1); + + if(i == unsigned(nDom) - 1) { + gammaInf(i) = Domain(2, 2); + } + + if(i == 0) { + gammaDir(i) = Domain(2, 1); + } + + border(i) = Domain(2, 100 + i); + + if(i != 0) { + sigma(i, i - 1) = Domain(2, 200 + i - 1); + } + if(i != unsigned(nDom) - 1) { + sigma(i, i + 1) = Domain(2, 200 + i); + } + } + /**** + * Define topology + ****/ + std::vector< std::vector< unsigned int > > topology(nDom); + for(unsigned int i = 0; i < unsigned(nDom); ++i) { + if(i != 0) { + topology[i].push_back(i - 1); + } + if(i != unsigned(nDom) - 1) { + topology[i].push_back(i + 1); + } + } + // Create DDM formulation + gmshddm::problem::Formulation< std::complex< double > > formulation("Maxwell", + topology); + // Create fields + SubdomainField< std::complex< double >, Form::Form1 > E( + "E", omega | gammaDir | gammaInf | sigma, + FunctionSpaceTypeForm1::HierarchicalHCurl, FEorder); + SubdomainField< std::complex< double >, Form::Form1 > lambda( + "lambda", gammaDir, FunctionSpaceTypeForm1::HierarchicalHCurl, + FEorder); + InterfaceField< std::complex< double >, Form::Form1 > g( + "g", sigma, FunctionSpaceTypeForm1::HierarchicalHCurl, 0); + + + // Tell to the formulation that g is field that have to be exchanged between + // subdomains + formulation.addInterfaceField(g); + VectorFunction< std::complex< double > > n = normal< std::complex< double > >(); + VectorFunction< std::complex< double > > functionZero = + vector< std::complex< double > >(0., 0., 0.); + for(unsigned int i = 0; i < unsigned(nDom); ++i) { + E(i).addConstraint(border(i), functionZero); + lambda(i).addConstraint(border(i), functionZero); + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { + const unsigned int j = topology[i][jj]; + g(i, j).addConstraint(border(i), functionZero); + } + } + + // Add terms to the formulation + for(unsigned int i = 0; i < unsigned(nDom); ++i) { + // VOLUME TERMS + formulation(i).integral(curl(dof(E(i))), curl(tf(E(i))), omega(i), gauss); + formulation(i).integral(-k * k * dof(E(i)), tf(E(i)), omega(i), gauss); + // Silver Muller + /* formulation(i).integral(-im * k * n % (dof(E(i)) % n), tf(E(i)), + gammaInf(i), gauss);*/ + formulation(i).integral(-im * k * dof(E(i)), tf(E(i)), gammaInf(i), + gauss); + // Physical source using Lagrange mutiplier + formulation(i).integral(dof(lambda(i)), tf(E(i)), gammaDir(i), + gauss); + formulation(i).integral(dof(E(i)), tf(lambda(i)), gammaDir(i), + gauss); + formulation(i).integral(formulation.physicalSource(Einc), tf(lambda(i)), + gammaDir(i), gauss); + + // Artificial source on transmission boundaries + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { + const unsigned int j = topology[i][jj]; + /* formulation(i).integral(-im * k * n % (dof(E(i)) % n), tf(E(i)), + sigma(i, j), gauss);*/ + formulation(i).integral(-im * k * dof(E(i)), tf(E(i)), sigma(i, j), + gauss); + formulation(i).integral(formulation.artificialSource(g(j, i)), tf(E(i)), + sigma(i, j), gauss); + } + // INTERFACE TERMS + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { + const unsigned int j = topology[i][jj]; + formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), + gauss); + formulation(i, j).integral(formulation.artificialSource(g(j, i)), + tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(-2. * im * k * E(i), tf(g(i, j)), + sigma(i, j), gauss); + /*formulation(i, j).integral(-2. * im * k * n % (E(i) % n), tf(g(i, j)), + sigma(i, j), gauss);*/ + } + } + formulation.pre(); + // Solve the DDM formulation + formulation.solve("gmres"); + //save field E + for(int i = 0; i < nDom; ++i) { + save(E(i), omega(i), "E" + std::to_string(i)); + } + } + + return 0; +} diff --git a/tutorials/maxwell/waveguide/CMakeLists.txt b/tutorials/maxwell/waveguide/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..7481f3ad58dcd475585824c7f40fddf332e191dc --- /dev/null +++ b/tutorials/maxwell/waveguide/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(tuto CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../tutorial.cmake) + +add_executable(tuto main.cpp ${EXTRA_INCS}) +target_link_libraries(tuto ${EXTRA_LIBS}) diff --git a/tutorials/maxwell/waveguide/main.cpp b/tutorials/maxwell/waveguide/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3d9c6310dfd15d8b37b661c0f4e74a89b765be62 --- /dev/null +++ b/tutorials/maxwell/waveguide/main.cpp @@ -0,0 +1,213 @@ +#include "gmsh.h" + +#include <gmshddm/Formulation.h> +#include <gmshddm/GmshDdm.h> +#include <gmshddm/MPIInterface.h> +#include <gmshfem/Message.h> + +using namespace gmshddm; +using namespace gmshddm::common; +using namespace gmshddm::domain; +using namespace gmshddm::problem; +using namespace gmshddm::field; +using namespace gmshddm::mpi; + +using namespace gmshfem; +using namespace gmshfem::domain; +using namespace gmshfem::equation; +using namespace gmshfem::problem; +using namespace gmshfem::field; +using namespace gmshfem::function; +using namespace gmshfem::post; + +void waveguide(const unsigned int nDom, const double l, const double L, const double lc) +{ + if(nDom < 2) { + gmshfem::msg::error << "nDom should be at least equal to 2." << gmshfem::msg::endl; + exit(0); + } + gmsh::model::add("waveguide"); + + const double domainSize = L / nDom; + + std::vector< std::vector< int > > borders(nDom); + std::vector< int > omegas(nDom); + std::vector< int > sigmas(nDom - 1); + + int p[2]; + + p[0] = gmsh::model::geo::addPoint(0., l, 0., lc); + p[1] = gmsh::model::geo::addPoint(0., 0., 0., lc); + + int line = gmsh::model::geo::addLine(p[0], p[1]); + gmsh::model::addPhysicalGroup(1, {line}, 1); + gmsh::model::setPhysicalName(1, 1, "gammaDir"); + + std::vector< int > cornerTags; + for(unsigned int i = 0; i < nDom; ++i) { + std::vector< std::pair< int, int > > extrudeEntities; + gmsh::model::geo::extrude({std::make_pair(1, line)}, domainSize, 0., 0., extrudeEntities, {static_cast< int >((domainSize) / lc)}, std::vector< double >(), true); + line = extrudeEntities[0].second; + omegas[i] = extrudeEntities[1].second; + borders[i].push_back(extrudeEntities[2].second); + borders[i].push_back(extrudeEntities[3].second); + sigmas[i] = line; + if(i != nDom - 1) { + cornerTags.push_back(3 + cornerTags.size()); + cornerTags.push_back(3 + cornerTags.size()); + } + } + + gmsh::model::addPhysicalGroup(1, {line}, 2); + gmsh::model::setPhysicalName(1, 2, "gammaInf"); + + for(unsigned int i = 0; i < borders.size(); ++i) { + gmsh::model::addPhysicalGroup(1, borders[i], 100 + i); + gmsh::model::setPhysicalName(1, 100 + i, "border_" + std::to_string(i)); + } + + for(unsigned int i = 0; i < omegas.size(); ++i) { + gmsh::model::addPhysicalGroup(2, {omegas[i]}, i + 1); + gmsh::model::setPhysicalName(2, i + 1, "omega_" + std::to_string(i)); + } + + for(unsigned int i = 0; i < sigmas.size(); ++i) { + gmsh::model::addPhysicalGroup(1, {sigmas[i]}, 200 + i); + gmsh::model::setPhysicalName(1, 200 + i, "sigma_" + std::to_string(i) + "_" + std::to_string(i + 1)); + } + + gmsh::model::addPhysicalGroup(0, cornerTags, 1000); + gmsh::model::setPhysicalName(0, 1000, "corner"); + + gmsh::model::geo::synchronize(); + gmsh::model::mesh::generate(); +} + +int main(int argc, char **argv) +{ + GmshDdm gmshDdm(argc, argv); + + double pi = 3.14159265359; + + int nDom = 3; + gmshDdm.userDefinedParameter(nDom, "nDom"); + double L = 5; + double l = 1; + double lc = 0.02; + double k = 2 * pi; + double kIn = 1.; + int order = 1; + gmshDdm.userDefinedParameter(L, "L"); + gmshDdm.userDefinedParameter(l, "l"); + gmshDdm.userDefinedParameter(lc, "lc"); + gmshDdm.userDefinedParameter(k, "k"); + gmshDdm.userDefinedParameter(kIn, "kIn"); + gmshDdm.userDefinedParameter(order, "order"); + std::string gauss = "Gauss10"; + gmshDdm.userDefinedParameter(gauss, "gauss"); + + // Build geometry and mesh using gmsh API + waveguide(nDom, l, L, lc); + + // Define domain + Subdomain omega(nDom); + Subdomain gammaInf(nDom); + Subdomain gammaDir(nDom); + Subdomain border(nDom); + Interface sigma(nDom); + Domain corner = Domain("corner"); + + for(unsigned int i = 0; i < nDom; ++i) { + omega(i) = Domain(2, i + 1); + + if(i == nDom - 1) { + gammaInf(i) = Domain(1, 2); + } + + if(i == 0) { + gammaDir(i) = Domain(1, 1); + } + + border(i) = Domain(1, 100 + i); + + if(i != 0) { + sigma(i, i - 1) = Domain(1, 200 + i - 1); + } + if(i != nDom - 1) { + sigma(i, i + 1) = Domain(1, 200 + i); + } + } + + // Define topology + std::vector< std::vector< unsigned int > > topology(nDom); + for(unsigned int i = 0; i < nDom; ++i) { + if(i != 0) { + topology[i].push_back(i - 1); + } + if(i != nDom - 1) { + topology[i].push_back(i + 1); + } + } + + + // Create DDM formulation + gmshddm::problem::Formulation< std::complex< double > > formulation("Maxwell", topology); + + // Create fields + SubdomainField< std::complex< double >, Form::Form1 > e("e", omega | gammaDir | gammaInf | border | sigma, FunctionSpaceTypeForm1::HierarchicalHCurl, order); + for(unsigned int i = 0; i < nDom; ++i) { + e(i).addConstraint(border(i), vector< std::complex< double > >(0., 0., 0.)); + } + e(0).addConstraint(gammaDir(0), sin< std::complex< double > >(2. * pi * y< std::complex< double > >() / l) * vector< std::complex< double > >(0., 1., 0.)); + InterfaceField< std::complex< double >, Form::Form1 > g("g", sigma, FunctionSpaceTypeForm1::HierarchicalHCurl, order); + for(unsigned int i = 0; i < nDom; ++i) { + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { + const unsigned int j = topology[i][jj]; + g(i, j).addConstraint(border(i), vector< std::complex< double > >(0., 0., 0.)); + } + } + + std::complex< double > im(0., 1.); + + // Tell to the formulation that g is field that have to be exchanged between subdomains + formulation.addInterfaceField(g); + + VectorFunction< std::complex< double > > n = normal< std::complex< double > >(); + + // Add terms to the formulation + for(unsigned int i = 0; i < nDom; ++i) { + // VOLUME TERMS + formulation(i).integral(curl(dof(e(i))), curl(tf(e(i))), omega(i), gauss); + formulation(i).integral(-k * k * dof(e(i)), tf(e(i)), omega(i), gauss); + formulation(i).integral(-im * k * n % (dof(e(i)) % n), tf(e(i)), gammaInf(i), gauss); + + // Artificial source + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { + const unsigned int j = topology[i][jj]; + formulation(i).integral(-im * k * n % (dof(e(i)) % n), tf(e(i)), sigma(i, j), gauss); + formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(e(i)), sigma(i, j), gauss); + } + + // INTERFACE TERMS + for(unsigned int jj = 0; jj < topology[i].size(); ++jj) { + const unsigned int j = topology[i][jj]; + formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(formulation.artificialSource(g(j, i)), tf(g(i, j)), sigma(i, j), gauss); + formulation(i, j).integral(2. * im * k * n % (e(i) % n), tf(g(i, j)), sigma(i, j), gauss); + } + } + + // Solve the DDM formulation + formulation.pre(); + int iterMax = 1000; + gmshDdm.userDefinedParameter(iterMax, "iter"); + formulation.solve("gmres", 1e-6, iterMax); + + for(unsigned int i = 0; i < nDom; ++i) { + if(isItMySubdomain(i)) { + save(e(i), omega(i), "e_" + std::to_string(i)); + } + } + + return 0; +} diff --git a/tutorials/tutorial.cmake b/tutorials/tutorial.cmake new file mode 100644 index 0000000000000000000000000000000000000000..3e2c7ef61683c3ee4d8e3652f5103ed81c1c80c3 --- /dev/null +++ b/tutorials/tutorial.cmake @@ -0,0 +1,98 @@ +# The following commands are useful for the demos, and can be included in each +# demos's CMakeLists.txt; they fill the EXTRA_INCS and EXTRA_LIBS variables. + +# Detect if we use the MinGW compilers on Cygwin - if we do, handle the build as +# a pure Windows build +if(CYGWIN OR MSYS) + if(CMAKE_CXX_COMPILER_ID MATCHES "GNU" OR + CMAKE_CXX_COMPILER_ID MATCHES "Clang") + execute_process(COMMAND ${CMAKE_CXX_COMPILER} -dumpmachine + OUTPUT_VARIABLE CXX_COMPILER_MACHINE + OUTPUT_STRIP_TRAILING_WHITESPACE) + if(CXX_COMPILER_MACHINE MATCHES "mingw") + set(WIN32 1) + add_definitions(-DWIN32) + set(CMAKE_FIND_LIBRARY_PREFIXES "lib" "") + set(CMAKE_FIND_LIBRARY_SUFFIXES ".a" ".so" ".lib" ".LIB" ".dll" ".DLL") + endif() + endif() +endif() + +# Make sure we have C++11 support +set(CMAKE_CXX_STANDARD 11) + +# Enable all warnings +include(CheckCXXCompilerFlag) +check_cxx_compiler_flag("-Wall" HAVE_WALL) +if(HAVE_WALL) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall") +endif() + +# Find Gmsh +find_library(GMSH_LIB gmsh) +if(NOT GMSH_LIB) + message(FATAL_ERROR "Could not find Gmsh library") +endif(NOT GMSH_LIB) +find_path(GMSH_INC gmsh.h) +if(NOT GMSH_INC) + message(FATAL_ERROR "Could not find Gmsh header 'gmsh.h'") +endif(NOT GMSH_INC) +list(APPEND EXTRA_INCS ${GMSH_INC}) +list(APPEND EXTRA_LIBS ${GMSH_LIB}) + +# Find GmshFEM +find_library(GMSHFEM_LIB gmshfem) +if(NOT GMSHFEM_LIB) + message(FATAL_ERROR "Could not find GmshFEM library") +endif(NOT GMSHFEM_LIB) +find_path(GMSHFEM_INC gmshfem/GmshFem.h) +if(NOT GMSHFEM_INC) + message(FATAL_ERROR "Could not find GmshFEM header 'gmshfem/GmshFem.h'") +endif(NOT GMSHFEM_INC) +list(APPEND EXTRA_LIBS ${GMSHFEM_LIB}) +list(APPEND EXTRA_INCS ${GMSHFEM_INC}) +list(APPEND EXTRA_INCS ${GMSHFEM_INC}/gmshfem/contrib) + +# Find GmshDDM +find_library(GMSHDDM_LIB gmshddm) +if(NOT GMSHDDM_LIB) + message(FATAL_ERROR "Could not find GmshDDM library") +endif(NOT GMSHDDM_LIB) + +find_path(GMSHDDM_INC gmshddm/GmshDdm.h) +if(NOT GMSHDDM_INC) + message(FATAL_ERROR "Could not find GmshDDM header 'gmshddm/GmshDdm.h'") +endif(NOT GMSHDDM_INC) +list(APPEND EXTRA_LIBS ${GMSHDDM_LIB}) +list(APPEND EXTRA_INCS ${GMSHDDM_INC}) + +# OpenMP +find_package(OpenMP) +if(OPENMP_FOUND) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +endif() + +if(APPLE AND EXISTS "/usr/local/opt/libomp") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xpreprocessor -fopenmp -I/usr/local/opt/libomp/include") + list(APPEND EXTRA_LIBS "-L/usr/local/opt/libomp/lib -lomp") + message(STATUS "Found OpenMP[Homebrew]") +elseif(APPLE AND EXISTS "/opt/local/lib/libomp") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xpreprocessor -fopenmp -I/opt/local/include/libomp") + list(APPEND EXTRA_LIBS "-L/opt/local/lib/libomp -lomp") + message(STATUS "Found OpenMP[MacPorts]") +endif() + +# MPI +find_package(MPI) +if(MPI_FOUND) + list(APPEND EXTRA_INCS ${MPI_CXX_INCLUDE_PATH}) + list(APPEND EXTRA_LIBS ${MPI_CXX_LIBRARIES}) + set(CMAKE_C_COMPILER ${MPI_C_COMPILER}) + set(CMAKE_CXX_COMPILER ${MPI_CXX_COMPILER}) + set(CMAKE_Fortran_COMPILER ${MPI_Fortran_COMPILER}) +endif() + +# Eigen +list(APPEND EXTRA_INCS ${CMAKE_CURRENT_LIST_DIR}/../contrib/eigen) + +include_directories(${EXTRA_INCS})