add edge and corner connections in domain decomposition

This commit is contained in:
Arne Morten Kvarving 2016-08-29 15:55:24 +02:00
parent 8d5dddf382
commit b486773a67
17 changed files with 479 additions and 23 deletions

View File

@ -14,12 +14,15 @@
#include "DomainDecomposition.h"
#include "ASMbase.h"
#include "ASMstruct.h"
#include "ASM2D.h"
#include "ASM3D.h"
#include "LinSolParams.h"
#include "ProcessAdm.h"
#include "SAMpatch.h"
#include "SIMbase.h"
#include "IFEM.h"
#include "Utilities.h"
#include "Vec3.h"
#include <cassert>
#include <numeric>
#include <iostream>
@ -33,12 +36,34 @@ public:
//! \param orient Orientation of boundary on slave
//! \param lIdx Index of boundary on slave
//! \param basis Basis to use for boundary on slave
NodeIterator(const ASMstruct* pch, int orient, int lIdx, int basis)
//! \param dim Dimension to iterate over
NodeIterator(const ASMstruct* pch, int orient, int lIdx, int basis, int dim = 2)
{
if (dim == 0) {
nodes.resize(1);
return;
}
int n1, n2, n3;
pch->getSize(n1,n2,n3,basis);
int nsd = pch->getNoSpaceDim();
if (nsd == 3) {
if (dim == 1) {
if (lIdx <= 4)
nodes.resize(n1);
else if (lIdx >= 5 && lIdx <= 8)
nodes.resize(n2);
else
nodes.resize(n3);
if (orient == 0)
std::iota(nodes.begin(), nodes.end(), 0);
else
std::iota(nodes.rbegin(), nodes.rend(), 0);
return;
}
int dim1, dim2;
if (lIdx == 1 || lIdx == 2)
dim1 = n2, dim2 = n3;
@ -270,6 +295,58 @@ std::vector<IntVec> DomainDecomposition::calcSubdomains3D(size_t nel1, size_t ne
}
void DomainDecomposition::setupNodeNumbers(int basis, IntVec& lNodes,
std::set<int>& cbasis,
const ASMbase* pch,
int dim, int lidx)
{
if (basis != 0) // specified base
cbasis = utl::getDigits(basis);
else if (dim == 0 || (dim == 1 && pch->getNoSpaceDim() == 3)) {
// need to expand to all bases for corners and edges
for (size_t b = 1; b <= pch->getNoBasis(); ++b)
cbasis.insert(b);
} else // directly add nodes, cbasis remains empty
pch->getBoundaryNodes(lidx, lNodes);
const ASM2D* pch2D = dynamic_cast<const ASM2D*>(pch);
const ASM3D* pch3D = dynamic_cast<const ASM3D*>(pch);
for (auto& it2 : cbasis) {
if (dim == 0) {
int node = 0;
if (pch2D) {
static std::map<int, std::array<int,2>> vmap2D
{{1, {-1, -1}},
{2, { 1, -1}},
{3, {-1, 1}},
{4, { 1, 1}}};
auto itc = vmap2D.find(lidx);
node = pch2D->getCorner(itc->second[0], itc->second[1], it2);
} else if (pch3D) {
static std::map<int, std::array<int,3>> vmap3D
{{1, {-1, -1, -1}},
{2, { 1, -1, -1}},
{3, {-1, 1, -1}},
{4, { 1, 1, -1}},
{5, {-1, -1, 1}},
{6, { 1, -1, 1}},
{7, {-1, 1, 1}},
{8, { 1, 1, 1}}};
auto itc = vmap3D.find(lidx);
node = pch3D->getCorner(itc->second[0], itc->second[1],
itc->second[2], it2);
}
lNodes.push_back(pch->getNodeID(node));
} else if (dim == 1 && pch3D) {
std::vector<int> eNodes = pch3D->getEdge(lidx, false, it2);
for (const int& it : eNodes)
lNodes.push_back(pch->getNodeID(it));
} else
pch->getBoundaryNodes(lidx, lNodes, it2);
}
}
bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
const SIMbase& sim)
{
@ -323,12 +400,8 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
std::set<int> cbasis;
IntVec lNodes;
if (it.basis != 0) {
cbasis = utl::getDigits(it.basis);
for (auto& it2 : cbasis)
sim.getPatch(sidx)->getBoundaryNodes(it.sidx, lNodes, it2);
} else
sim.getPatch(sidx)->getBoundaryNodes(it.sidx, lNodes);
setupNodeNumbers(it.basis, lNodes, cbasis, sim.getPatch(sidx), it.dim, it.sidx);
int nRecv;
adm.receive(nRecv, getPatchOwner(it.master));
@ -344,7 +417,7 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
for (size_t b = 1; b <= sim.getPatch(sidx)->getNoBasis(); ++b) {
if (cbasis.empty() || cbasis.find(b) != cbasis.end()) {
NodeIterator iter(dynamic_cast<const ASMstruct*>(sim.getPatch(sidx)),
it.orient, it.sidx, b);
it.orient, it.sidx, b, it.dim);
auto it_n = iter.begin();
for (size_t i = 0; i < iter.size(); ++i, ++it_n) {
int node = MLGN[lNodes[i+ofs]-1];
@ -389,13 +462,8 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
continue;
std::set<int> cbasis;
if (it.basis != 0)
cbasis = utl::getDigits(it.basis);
std::vector<int> glbNodes;
for (size_t b = 1; b <= sim.getPatch(midx)->getNoBasis(); ++b)
if (cbasis.empty() || cbasis.find(b) != cbasis.end())
sim.getPatch(midx)->getBoundaryNodes(it.midx, glbNodes, b);
setupNodeNumbers(it.basis, glbNodes, cbasis, sim.getPatch(midx), it.dim, it.midx);
for (size_t i = 0; i < glbNodes.size(); ++i)
glbNodes[i] = MLGN[glbNodes[i]-1];
@ -410,7 +478,9 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
std::vector<int> DomainDecomposition::setupEquationNumbers(const SIMbase& sim,
int pidx, int lidx, const std::set<int>& cbasis)
int pidx, int lidx,
const std::set<int>& cbasis,
int dim)
{
std::vector<IntVec> lNodes(sim.getPatch(pidx)->getNoBasis());
std::vector<int> result;
@ -427,8 +497,10 @@ std::vector<int> DomainDecomposition::setupEquationNumbers(const SIMbase& sim,
if (!cbasis.empty() && cbasis.find(basis) == cbasis.end())
continue;
if (lNodes[basis-1].empty())
sim.getPatch(pidx)->getBoundaryNodes(lidx, lNodes[basis-1], basis);
if (lNodes[basis-1].empty()) {
IntSet dummy;
setupNodeNumbers(basis, lNodes[basis-1], dummy, sim.getPatch(pidx), dim, lidx);
}
std::set<int> components;
if (block == 0 || sim.getSolParams()->getBlock(block-1).comps == 0) {
@ -552,7 +624,7 @@ bool DomainDecomposition::calcGlobalEqNumbers(const ProcessAdm& adm,
if (it.basis != 0)
cbasis = utl::getDigits(it.basis);
IntVec locEqs = setupEquationNumbers(sim, sidx, it.sidx, cbasis);
IntVec locEqs = setupEquationNumbers(sim, sidx, it.sidx, cbasis, it.dim);
int nRecv;
adm.receive(nRecv, getPatchOwner(it.master));
@ -587,7 +659,7 @@ bool DomainDecomposition::calcGlobalEqNumbers(const ProcessAdm& adm,
components = utl::getDigits(sim.getSolParams()->getBlock(block-1).comps);
NodeIterator iter(dynamic_cast<const ASMstruct*>(sim.getPatch(sidx)),
it.orient, it.sidx, basis);
it.orient, it.sidx, basis, it.dim);
auto it_n = iter.begin();
int nodeDofs = components.size();
for (size_t i = 0; i < iter.size(); ++i, ++it_n) {
@ -667,7 +739,7 @@ bool DomainDecomposition::calcGlobalEqNumbers(const ProcessAdm& adm,
if (it.basis != 0)
cbasis = utl::getDigits(it.basis);
IntVec glbEqs = setupEquationNumbers(sim, midx, it.midx, cbasis);
IntVec glbEqs = setupEquationNumbers(sim, midx, it.midx, cbasis, it.dim);
adm.send(int(glbEqs.size()), getPatchOwner(it.slave));
adm.send(glbEqs, getPatchOwner(it.slave));
}
@ -696,10 +768,96 @@ int DomainDecomposition::getGlobalEq(int lEq, size_t idx) const
}
bool DomainDecomposition::sanityCheckCorners(const SIMbase& sim)
{
#ifdef HAVE_MPI
const ProcessAdm& adm = sim.getProcessAdm();
if (!adm.isParallel())
return true;
std::vector<int> sizes(adm.getNoProcs());
for (int i = 1; i <= sim.getNoPatches(); ++i)
sizes[adm.dd.getPatchOwner(i)] += 5*sim.getPatch(1)->getNoBasis() *
pow(2,sim.getNoSpaceDim());
std::vector<int> displ(adm.getNoProcs());
for (int i = 1; i < adm.getNoProcs(); ++i)
displ[i] = displ[i-1] + sizes[i-1];
std::vector<double> loc_data;
for (int i = 1; i <= sim.getNoPatches(); ++i) {
int pIdx;
if ((pIdx=sim.getLocalPatchIndex(i)) > 0) {
const ASM2D* pch2D = nullptr;
if (sim.getNoSpaceDim() == 2)
pch2D = dynamic_cast<const ASM2D*>(sim.getPatch(pIdx));
const ASM3D* pch3D = nullptr;
if (sim.getNoSpaceDim() == 3)
pch3D = dynamic_cast<const ASM3D*>(sim.getPatch(pIdx));
if (pch2D || pch3D) {
for (size_t c = 0; c < pow(2,sim.getPatch(pIdx)->getNoSpaceDim()); ++c) {
for (size_t b = 1; b <= sim.getPatch(pIdx)->getNoBasis(); ++b) {
int node;
if (pch2D)
node = pch2D->getCorner((c == 1 || c == 3) ? 1 : -1,
c >= 2 ? 1 : -1, b);
else if (pch3D)
node = pch3D->getCorner((c == 1 || c == 3 || c == 5 || c == 7) ? 1 : -1,
(c == 2 || c == 3 || c == 6 || c == 7) ? 1 : -1,
c >= 4 ? 1 : -1, b);
Vec3 pos = sim.getPatch(pIdx)->getCoord(node);
loc_data.push_back(pos[0]);
loc_data.push_back(pos[1]);
loc_data.push_back(pos[2]);
loc_data.push_back(b);
loc_data.push_back(adm.dd.getMLGN()[sim.getPatch(pIdx)->getNodeID(node)-1]);
}
}
}
}
}
std::vector<double> glob_data(displ.back()+sizes.back());
MPI_Allgatherv(loc_data.data(), loc_data.size(), MPI_DOUBLE,
glob_data.data(), sizes.data(), displ.data(), MPI_DOUBLE,
*adm.getCommunicator());
// unpack data
std::vector<std::array<double,5>> corners;
for (size_t i = 0; i < glob_data.size(); i += 5)
corners.push_back({glob_data[i], glob_data[i+1],
glob_data[i+2], glob_data[i+3], glob_data[i+4]});
for (const auto& c : corners) {
auto fail = std::find_if(corners.begin(), corners.end(),
[c](const std::array<double,5>& C)
{
return std::fabs(C[0]-c[0]) < 1e-6 &&
std::fabs(C[1]-c[1]) < 1e-6 &&
std::fabs(C[2]-c[2]) < 1e-6 &&
(int)c[3] == (int)C[3] &&
(int)c[4] != (int)C[4];
});
if (fail != corners.end()) {
std::cerr << "** DomainDecomposition::setup ** Corner node " << (int)c[4]
<< " with coordinates (" << c[0] << "," << c[1] << "," << c[2] << ")"
<< " found with different global ID " << (int)(*fail)[4]
<< " for basis " << (int)c[3] << "."
<< " You have to add vertex / edge connections." << std::endl;
return false;
}
}
#endif
return true;
}
bool DomainDecomposition::setup(const ProcessAdm& adm, const SIMbase& sim)
{
#ifdef HAVE_MPI
if (adm.isParallel()) {
if (!adm.isParallel())
IFEM::cout << "Establishing domain decomposition" << std::endl;
#if SP_DEBUG > 1
@ -709,8 +867,8 @@ bool DomainDecomposition::setup(const ProcessAdm& adm, const SIMbase& sim)
", slave/idx=" << it.slave << "/" << it.sidx <<
", orient=" << it.orient << ", dim=" << it.dim <<
", basis=" << it.basis << std::endl;
#endif
}
#endif
#endif
sam = dynamic_cast<const SAMpatch*>(sim.getSAM());
@ -736,7 +894,12 @@ bool DomainDecomposition::setup(const ProcessAdm& adm, const SIMbase& sim)
if (ok < adm.getNoProcs())
return false;
// sanity check all corners of the patches
if (!sanityCheckCorners(sim))
return false;
ok = 1;
// Establish local equation mappings for each block.
if (sim.getSolParams() && sim.getSolParams()->getNoBlocks() > 1) {
IFEM::cout << "\tEstablishing local block equation numbers" << std::endl;

View File

@ -19,6 +19,7 @@
#include <vector>
#include <cstddef>
class ASMbase;
class ProcessAdm;
class SAMpatch;
class SIMbase;
@ -64,6 +65,10 @@ public:
if (A.slave != B.slave)
return A.slave < B.slave;
// then by dim
if (A.dim != B.dim)
return A.dim < B.dim;
// finally by index on master
return A.midx < B.midx;
}
@ -181,7 +186,20 @@ private:
//! \param cbasis If non-empty, bases to connect
std::vector<int> setupEquationNumbers(const SIMbase& sim,
int pidx, int lidx,
const std::set<int>& cbasis);
const std::set<int>& cbasis,
int dim);
//! \brief Setup node numbers for all bases on a boundary.
//! \param basis Bases to grab nodes for
//! \param lNodes Resulting nodes
//! \param cbasis Bases to connect
//! \param pch Patch to obtain nodes for
//! \param dim Dimension of boundary to obtain nodes for
//! \param lidx Local index of boundary to obtain nodes for
void setupNodeNumbers(int basis, std::vector<int>& lNodes,
std::set<int>& cbasis,
const ASMbase* pch,
int dim, int lidx);
//! \brief Calculate the global node numbers for given finite element model.
bool calcGlobalNodeNumbers(const ProcessAdm& adm, const SIMbase& sim);
@ -189,6 +207,11 @@ private:
//! \brief Calculate the global equation numbers for given finite element model.
bool calcGlobalEqNumbers(const ProcessAdm& adm, const SIMbase& sim);
//! \brief Sanity check model.
//! \details Collects the corners of all patches in the model and make sure nodes
//! with matching coordinates have the same global node ID.
bool sanityCheckCorners(const SIMbase& sim);
std::map<int,int> patchOwner; //!< Process that owns a particular patch
//! \brief Struct with information per matrix block.

View File

@ -13,6 +13,7 @@
#include "DomainDecomposition.h"
#include "SAM.h"
#include "ASMbase.h"
#include "IntegrandBase.h"
#include "SIM2D.h"
#include "SIM3D.h"
#include "IFEM.h"
@ -44,6 +45,14 @@ protected:
};
template<class Dim>
class DummySIM : public Dim {
public:
class DummyIntegrand : public IntegrandBase {};
DummySIM() : Dim(1) { Dim::myProblem = new DummyIntegrand; }
};
typedef std::vector<int> IntVec;
static IntVec readIntVector(const std::string& file)
@ -84,6 +93,29 @@ class TestDomainDecomposition3D : public testing::Test,
};
TEST_P(TestDomainDecomposition2D, Corner)
{
DummySIM<SIM2D> sim;
std::stringstream str;
str << "src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner";
if (GetParam() == 1)
str << "_fail";
str << ".xinp";
sim.read(str.str().c_str());
if (GetParam() == 0) {
ASSERT_TRUE(sim.preprocess());
const ProcessAdm& adm = sim.getProcessAdm();
str.str("");
str << "src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes";
str << adm.getProcId() << ".ref";
IntVec B = readIntVector(str.str());
check_intvectors_equal(adm.dd.getMLGN(), B);
} else
ASSERT_FALSE(sim.preprocess());
}
TEST_P(TestDomainDecomposition2D, SetupSingleBasis)
{
SIM2D sim(2);
@ -468,6 +500,32 @@ TEST_P(TestDomainDecomposition3D, SetupMixedBasisPeriodicLM)
}
TEST_P(TestDomainDecomposition3D, Corner)
{
if (GetParam() > 1)
return;
DummySIM<SIM3D> sim;
std::stringstream str;
str << "src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner";
if (GetParam() == 1)
str << "_fail";
str << ".xinp";
sim.read(str.str().c_str());
if (GetParam() == 0) {
ASSERT_TRUE(sim.preprocess());
const ProcessAdm& adm = sim.getProcessAdm();
str.str("");
str << "src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes";
str << adm.getProcId() << ".ref";
IntVec B = readIntVector(str.str());
check_intvectors_equal(adm.dd.getMLGN(), B);
} else
ASSERT_FALSE(sim.preprocess());
}
const std::vector<int> orientations2D = {0,1};
INSTANTIATE_TEST_CASE_P(TestDomainDecomposition2D, TestDomainDecomposition2D, testing::ValuesIn(orientations2D));
const std::vector<int> orientations3D = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};

View File

@ -0,0 +1,2 @@
9
1 2 3 4 5 6 7 8 9

View File

@ -0,0 +1,2 @@
9
3 10 11 6 12 13 9 14 15

View File

@ -0,0 +1,2 @@
9
16 17 18 19 20 21 11 22 23

View File

@ -0,0 +1,2 @@
9
11 22 23 13 24 25 15 26 27

View File

@ -0,0 +1,22 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<simulation>
<console>
<logging output_prefix="foo"/>
</console>
<!-- General - geometry definitions !-->
<geometry>
<partitioning procs="4" nperproc="1"/>
<patchfile>src/ASM/Test/refdata/backstep2D-4.g2</patchfile>
<refine lowerpatch="1" upperpatch="4" u="1" v="1"/>
<topology>
<connection master="1" medge="2" slave="2" sedge="1"/>
<connection master="2" medge="2" slave="4" sedge="1"/>
<connection master="3" medge="4" slave="4" sedge="3"/>
<connection master="2" midx="2" slave="3" sidx="3" dim="0"/>
</topology>
</geometry>
</simulation>

View File

@ -0,0 +1,21 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<simulation>
<console>
<logging output_prefix="foo"/>
</console>
<!-- General - geometry definitions !-->
<geometry>
<partitioning procs="4" nperproc="1"/>
<patchfile>src/ASM/Test/refdata/backstep2D-4.g2</patchfile>
<refine lowerpatch="1" upperpatch="4" u="1" v="1"/>
<topology>
<connection master="1" medge="2" slave="2" sedge="1"/>
<connection master="2" medge="2" slave="4" sedge="1"/>
<connection master="3" medge="4" slave="4" sedge="3"/>
</topology>
</geometry>
</simulation>

View File

@ -0,0 +1,2 @@
27
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

View File

@ -0,0 +1,2 @@
27
3 28 29 6 30 31 9 32 33 12 34 35 15 36 37 18 38 39 21 40 41 24 42 43 27 44 45

View File

@ -0,0 +1,2 @@
27
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 29 64 65 31 66 67 33 68 69

View File

@ -0,0 +1,2 @@
27
29 64 65 31 66 67 33 68 69 35 70 71 37 72 73 39 74 75 41 76 77 43 78 79 45 80 81

View File

@ -0,0 +1,22 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<simulation>
<console>
<logging output_prefix="foo"/>
</console>
<!-- General - geometry definitions !-->
<geometry>
<partitioning procs="4" nperproc="1"/>
<patchfile>src/ASM/Test/refdata/backstep-4.g2</patchfile>
<refine lowerpatch="1" upperpatch="4" u="1" v="1" w="1"/>
<topology>
<connection master="3" mface="6" slave="4" sface="5"/>
<connection master="2" mface="2" slave="4" sface="1"/>
<connection master="1" mface="2" slave="2" sface="1"/>
<connection master="2" midx="6" slave="3" sidx="7" dim="1"/>
</topology>
</geometry>
</simulation>

View File

@ -0,0 +1,21 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<simulation>
<console>
<logging output_prefix="foo"/>
</console>
<!-- General - geometry definitions !-->
<geometry>
<partitioning procs="4" nperproc="1"/>
<patchfile>src/ASM/Test/refdata/backstep-4.g2</patchfile>
<refine lowerpatch="1" upperpatch="4" u="1" v="1" w="1"/>
<topology>
<connection master="3" mface="6" slave="4" sface="5"/>
<connection master="2" mface="2" slave="4" sface="1"/>
<connection master="1" mface="2" slave="2" sface="1"/>
</topology>
</geometry>
</simulation>

View File

@ -0,0 +1,67 @@
700 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
2 2
0 0 1 1
0 0 0
1 0 0
0 1 0
1 1 0
0 0 1
1 0 1
0 1 1
1 1 1
700 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
2 2
0 0 1 1
1 0 0
2 0 0
1 1 0
2 1 0
1 0 1
2 0 1
1 1 1
2 1 1
700 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
2 2
0 0 1 1
2 0 -1
3 0 -1
2 1 -1
3 1 -1
2 0 0
3 0 0
2 1 0
3 1 0
700 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
2 2
0 0 1 1
2 0 0
3 0 0
2 1 0
3 1 0
2 0 1
3 0 1
2 1 1
3 1 1

View File

@ -0,0 +1,43 @@
200 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
0 0 0
1 0 0
0 1 0
1 1 0
200 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
1 0 0
2 0 0
1 1 0
2 1 0
200 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
2 -1 0
3 -1 0
2 0 0
3 0 0
200 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
2 0 0
3 0 0
2 1 0
3 1 0