Skip to content
Snippets Groups Projects
Commit 06a5981f authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

remove unused + bcgs by default

parent 1c0d7c25
No related branches found
No related tags found
No related merge requests found
...@@ -10,7 +10,7 @@ using gmshfem::equation::tf; ...@@ -10,7 +10,7 @@ using gmshfem::equation::tf;
using gmshfem::function::operator-; using gmshfem::function::operator-;
using gmshfem::function::operator*; using gmshfem::function::operator*;
void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, double pmlSize, const double lc, bool pml, const int order, const unsigned int sourcePosX, const unsigned int sourcePosY) void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, const double lc, const int order, const unsigned int sourcePosX, const unsigned int sourcePosY)
{ {
// *********** // ***********
// j * * // j * *
...@@ -41,58 +41,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub ...@@ -41,58 +41,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub
} }
} }
// ******************
// pmls
// ******************
// 1 0
// 0 ****(1)**** 1
// * *
// (2) (0)
// * *
// 1 ****(3)**** 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::pair< int, 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);
}
}
// ******************
// pmlsSquare
// ******************
// (1)*******(0)
// * *
// * *
// * *
// (2)*******(3)
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< std::vector< std::vector< int > > > corners(nDomX);
// for(unsigned int i = 0; i < nDomX; ++i) {
// corners[i].resize(nDomY);
// for(unsigned int j = 0; j < nDomY; ++j) {
// corners[i][j].resize(4);
// }
// }
std::vector< int > left(nDomY); std::vector< int > left(nDomY);
int bottom = 0; int bottom = 0;
for(unsigned int i = 0; i < nDomX; ++i) { for(unsigned int i = 0; i < nDomX; ++i) {
...@@ -108,11 +56,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub ...@@ -108,11 +56,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub
p[2] = gmsh::model::geo::addPoint((i + 1) * xSize, (j + 1) * 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);
// corners[i][j][0] = p[2];
// corners[i][j][1] = p[3];
// corners[i][j][2] = p[0];
// corners[i][j][3] = p[1];
int l[4]; int l[4];
borders[i][j][0] = l[0] = gmsh::model::geo::addLine(p[1], p[2]); borders[i][j][0] = l[0] = gmsh::model::geo::addLine(p[1], p[2]);
borders[i][j][1] = l[1] = gmsh::model::geo::addLine(p[2], p[3]); borders[i][j][1] = l[1] = gmsh::model::geo::addLine(p[2], p[3]);
...@@ -151,45 +94,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub ...@@ -151,45 +94,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub
gmsh::model::setPhysicalName(0, 1, "source"); gmsh::model::setPhysicalName(0, 1, "source");
} }
if(pml) {
// pmls
const double extrudeX[4] = {pmlSize, 0., -pmlSize, 0.};
const double extrudeY[4] = {0., pmlSize, 0., -pmlSize};
for(unsigned int k = 0; k < 4; ++k) {
std::vector< std::pair< int, int > > extrudeEntities;
gmsh::model::geo::extrude({std::make_pair(1, -borders[i][j][k])}, extrudeX[k], extrudeY[k], 0., extrudeEntities, {static_cast< int >(/*(k % 2 ? 1. : 2.) */ pmlSize / lc)}, std::vector< double >(), true);
pmlsInf[i][j][k] = extrudeEntities[0].second;
pmls[i][j][k] = extrudeEntities[1].second;
if(-borders[i][j][k] == extrudeEntities[2].second) {
pmlsBnd[i][j][k].first = extrudeEntities[3].second;
pmlsBnd[i][j][k].second = extrudeEntities[4].second;
}
else {
pmlsBnd[i][j][k].first = extrudeEntities[2].second;
pmlsBnd[i][j][k].second = extrudeEntities[3].second;
}
}
// pmlsSquare
std::vector< std::pair< int, int > > extrudeEntities;
gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][0].second)}, 0., pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
pmlsSquare[i][j][0] = extrudeEntities[1].second;
extrudeEntities.clear();
gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][2].first)}, 0., pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
pmlsSquare[i][j][1] = extrudeEntities[1].second;
extrudeEntities.clear();
gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][2].second)}, 0., -pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
pmlsSquare[i][j][2] = extrudeEntities[1].second;
extrudeEntities.clear();
gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][0].first)}, 0., -pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
pmlsSquare[i][j][3] = extrudeEntities[1].second;
extrudeEntities.clear();
}
} }
} }
gmsh::model::geo::removeAllDuplicates(); gmsh::model::geo::removeAllDuplicates();
...@@ -216,68 +120,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub ...@@ -216,68 +120,6 @@ void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const doub
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)); gmsh::model::setPhysicalName(1, 40000 + i * nDomY + j + 1, "sigmaS_" + std::to_string(i) + "_" + std::to_string(j));
// gmsh::model::addPhysicalGroup(0, {corners[i][j][0]}, 10000 + i * nDomY + j + 1);
// gmsh::model::setPhysicalName(0, 10000 + i * nDomY + j + 1, "crossPointNE_" + std::to_string(i) + "_" + std::to_string(j));
//
// gmsh::model::addPhysicalGroup(0, {corners[i][j][1]}, 20000 + i * nDomY + j + 1);
// gmsh::model::setPhysicalName(0, 20000 + i * nDomY + j + 1, "crossPointNW_" + std::to_string(i) + "_" + std::to_string(j));
//
// gmsh::model::addPhysicalGroup(0, {corners[i][j][2]}, 30000 + i * nDomY + j + 1);
// gmsh::model::setPhysicalName(0, 30000 + i * nDomY + j + 1, "crossPointSW_" + std::to_string(i) + "_" + std::to_string(j));
//
// gmsh::model::addPhysicalGroup(0, {corners[i][j][3]}, 40000 + i * nDomY + j + 1);
// gmsh::model::setPhysicalName(0, 40000 + i * nDomY + j + 1, "crossPointSE_" + 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, "pmlE_" + 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, "pmlN_" + 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, "pmlW_" + 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, "pmlS_" + 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, "pmlNE_" + 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, "pmlNW_" + 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, "pmlSW_" + 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, "pmlSE_" + std::to_string(i) + "_" + std::to_string(j));
// Sigmas in PMLs
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0].first}, 100000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 100000 + i * nDomY + j + 1, "sigmaE_S_" + std::to_string(i) + "_" + std::to_string(j));
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0].second}, 150000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 150000 + i * nDomY + j + 1, "sigmaE_N_" + std::to_string(i) + "_" + std::to_string(j));
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1].first}, 200000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 200000 + i * nDomY + j + 1, "sigmaN_E_" + std::to_string(i) + "_" + std::to_string(j));
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1].second}, 250000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 250000 + i * nDomY + j + 1, "sigmaN_W_" + std::to_string(i) + "_" + std::to_string(j));
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2].first}, 300000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 300000 + i * nDomY + j + 1, "sigmaW_N_" + std::to_string(i) + "_" + std::to_string(j));
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2].second}, 350000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 350000 + i * nDomY + j + 1, "sigmaW_S_" + std::to_string(i) + "_" + std::to_string(j));
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3].first}, 400000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 400000 + i * nDomY + j + 1, "sigmaS_W_" + std::to_string(i) + "_" + std::to_string(j));
gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3].second}, 450000 + i * nDomY + j + 1);
gmsh::model::setPhysicalName(1, 450000 + i * nDomY + j + 1, "sigmaS_E_" + std::to_string(i) + "_" + std::to_string(j));
}
} }
} }
...@@ -296,10 +138,6 @@ int main(int argc, char **argv) ...@@ -296,10 +138,6 @@ int main(int argc, char **argv)
double xSize = 0.1, ySize = 0.1; double xSize = 0.1, ySize = 0.1;
gmshDdm.userDefinedParameter(xSize, "xSize"); gmshDdm.userDefinedParameter(xSize, "xSize");
gmshDdm.userDefinedParameter(ySize, "ySize"); gmshDdm.userDefinedParameter(ySize, "ySize");
double pmlSize = (xSize + ySize) / 10.;
gmshDdm.userDefinedParameter(pmlSize, "pmlSize");
bool withPml = false;
gmshDdm.userDefinedParameter(withPml, "pml");
double lc = 0.1; double lc = 0.1;
gmshDdm.userDefinedParameter(lc, "lc"); gmshDdm.userDefinedParameter(lc, "lc");
int meshOrder = 1; int meshOrder = 1;
...@@ -309,7 +147,7 @@ int main(int argc, char **argv) ...@@ -309,7 +147,7 @@ int main(int argc, char **argv)
gmshDdm.userDefinedParameter(sourcePosY, "sourcePosY"); gmshDdm.userDefinedParameter(sourcePosY, "sourcePosY");
// Build geometry and mesh using gmsh API // Build geometry and mesh using gmsh API
checkerboard(nDomX, nDomY, xSize, ySize, pmlSize, lc, withPml, meshOrder, sourcePosX, sourcePosY); checkerboard(nDomX, nDomY, xSize, ySize, lc, meshOrder, sourcePosX, sourcePosY);
const double pi = 3.14159265359; const double pi = 3.14159265359;
const std::complex< double > im(0., 1.); const std::complex< double > im(0., 1.);
...@@ -328,7 +166,6 @@ int main(int argc, char **argv) ...@@ -328,7 +166,6 @@ int main(int argc, char **argv)
std::vector< gmshddm::domain::Subdomain > sigma(4, nDom); std::vector< gmshddm::domain::Subdomain > sigma(4, nDom);
std::vector< gmshddm::domain::Interface > sigmaInterface(4, nDom); std::vector< gmshddm::domain::Interface > sigmaInterface(4, nDom);
std::vector< std::pair< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > corner(4, std::make_pair(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)) );
// Define topology // Define topology
std::vector< std::vector< unsigned int > > topology(nDom); std::vector< std::vector< unsigned int > > topology(nDom);
for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) { for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) {
...@@ -339,8 +176,6 @@ int main(int argc, char **argv) ...@@ -339,8 +176,6 @@ int main(int argc, char **argv)
for(unsigned int f = 0; f < 4; ++f) { 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);
// corner[f].first(index) = gmshfem::domain::Domain(0, ((f-1)%4 + 1) * 10000 + i * nDomY + j + 1);
// corner[f].second(index) = gmshfem::domain::Domain(0, (f + 1) * 10000 + i * nDomY + j + 1);
} }
if(i != static_cast< unsigned int >(nDomX) - 1) { // E if(i != static_cast< unsigned int >(nDomX) - 1) { // E
...@@ -432,7 +267,7 @@ int main(int argc, char **argv) ...@@ -432,7 +267,7 @@ int main(int argc, char **argv)
double tol = 1e-4; double tol = 1e-4;
gmshDdm.userDefinedParameter(tol, "tol"); gmshDdm.userDefinedParameter(tol, "tol");
std::string solver = "gmres"; std::string solver = "bcgs";
gmshDdm.userDefinedParameter(solver, "solver"); gmshDdm.userDefinedParameter(solver, "solver");
formulation.pre(); formulation.pre();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment