From aded054030b05da705c265862409b62e6e781331 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 29 Jan 2010 12:29:11 +0000
Subject: [PATCH] Revision of parallel domain decomposition code.
---
AMDiS/libtool | 24 +-
AMDiS/src/Element.cc | 67 ------
AMDiS/src/Element.h | 12 -
AMDiS/src/InteriorBoundary.cc | 34 +++
AMDiS/src/InteriorBoundary.h | 21 +-
AMDiS/src/ParMetisPartitioner.cc | 6 +-
AMDiS/src/ParallelDomainBase.cc | 365 ++++++++++++++++++++-----------
AMDiS/src/ParallelDomainBase.h | 37 +++-
AMDiS/src/StdMpi.h | 7 +
AMDiS/src/Tetrahedron.cc | 12 +-
10 files changed, 338 insertions(+), 247 deletions(-)
diff --git a/AMDiS/libtool b/AMDiS/libtool
index e6b2b7be..7a44f979 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host deimos103:
+# Libtool was configured on host p2q071:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -82,13 +82,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="gcc"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -171,7 +171,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host deimos103:
+# Libtool was configured on host p2q071:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6798,13 +6798,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="g++"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpiCC"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -6887,7 +6887,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6954,11 +6954,11 @@ predeps=""
# Dependencies to place after the objects being linked to create a
# shared library.
-postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
+postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -libverbs -lrt -lnuma -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path="-L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
+compiler_lib_search_path="-L/usr/lib64 -L/licsoft/libraries/openmpi/1.2.6/64bit/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host deimos103:
+# Libtool was configured on host p2q071:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7103,7 +7103,7 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc
index e251c974..62c5618f 100644
--- a/AMDiS/src/Element.cc
+++ b/AMDiS/src/Element.cc
@@ -566,71 +566,4 @@ namespace AMDiS {
return result;
}
- bool fitElementToMeshCode(RefinementManager *refineManager,
- MeshStructure &code,
- Element *el,
- int ithSide,
- int elType)
- {
- FUNCNAME("fitElementToMeshCode()");
-
- if (code.empty())
- return false;
-
- int s1 = el->getSideOfChild(0, ithSide, elType);
- int s2 = el->getSideOfChild(1, ithSide, elType);
-
- TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");
-
- if (s1 != -1 && s2 != -1)
- return fitElementToMeshCode2(refineManager, code, el, ithSide, elType);
-
- if (el->isLeaf()) {
- if (code.getNumElements() == 1 && code.isLeafElement())
- return false;
-
- el->setMark(1);
- refineManager->refineElement(el->getMesh(), el);
- }
-
- if (s1 != -1)
- return fitElementToMeshCode2(refineManager, code,
- el->getFirstChild(), s1, el->getChildType(elType));
- else
- return fitElementToMeshCode2(refineManager, code,
- el->getSecondChild(), s2, el->getChildType(elType));
- }
-
- bool fitElementToMeshCode2(RefinementManager *refineManager, MeshStructure &code,
- Element *el, int ithSide, int elType)
- {
- if (code.isLeafElement())
- return false;
-
- bool value = false;
-
- if (el->isLeaf()) {
- el->setMark(1);
- refineManager->refineElement(el->getMesh(), el);
- value = true;
- }
-
- int s1 = el->getSideOfChild(0, ithSide, elType);
- int s2 = el->getSideOfChild(1, ithSide, elType);
-
- if (s1 != -1) {
- code.nextElement();
- value |= fitElementToMeshCode2(refineManager, code, el->getFirstChild(),
- s1, el->getChildType(elType));
- }
-
- if (s2 != -1) {
- code.nextElement();
- value |= fitElementToMeshCode2(refineManager, code, el->getSecondChild(),
- s2, el->getChildType(elType));
- }
-
- return value;
- }
-
}
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index 3955d629..4cfc027e 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -574,18 +574,6 @@ namespace AMDiS {
friend class Mesh;
};
-
- bool fitElementToMeshCode(RefinementManager *refineManager,
- MeshStructure &code,
- Element *el,
- int ithSide,
- int elType);
-
- bool fitElementToMeshCode2(RefinementManager *refineManager,
- MeshStructure &code,
- Element *el,
- int ithSide,
- int elType);
}
#endif // AMDIS_ELEMENT_H
diff --git a/AMDiS/src/InteriorBoundary.cc b/AMDiS/src/InteriorBoundary.cc
index 2683b314..869f3705 100644
--- a/AMDiS/src/InteriorBoundary.cc
+++ b/AMDiS/src/InteriorBoundary.cc
@@ -76,12 +76,14 @@ namespace AMDiS {
SerUtil::serialize(out, bound.rankObj.subObj);
SerUtil::serialize(out, bound.rankObj.ithObj);
SerUtil::serialize(out, bound.rankObj.reverseMode);
+ serializeExcludeList(out, bound.rankObj.excludedSubstructures);
SerUtil::serialize(out, bound.neighObj.elIndex);
SerUtil::serialize(out, bound.neighObj.elType);
SerUtil::serialize(out, bound.neighObj.subObj);
SerUtil::serialize(out, bound.neighObj.ithObj);
SerUtil::serialize(out, bound.neighObj.reverseMode);
+ serializeExcludeList(out, bound.neighObj.excludedSubstructures);
}
}
}
@@ -107,16 +109,48 @@ namespace AMDiS {
SerUtil::deserialize(in, bound.rankObj.subObj);
SerUtil::deserialize(in, bound.rankObj.ithObj);
SerUtil::deserialize(in, bound.rankObj.reverseMode);
+ deserializeExcludeList(in, bound.rankObj.excludedSubstructures);
SerUtil::deserialize(in, bound.neighObj.elIndex);
SerUtil::deserialize(in, bound.neighObj.elType);
SerUtil::deserialize(in, bound.neighObj.subObj);
SerUtil::deserialize(in, bound.neighObj.ithObj);
SerUtil::deserialize(in, bound.neighObj.reverseMode);
+ deserializeExcludeList(in, bound.neighObj.excludedSubstructures);
bound.rankObj.el = elIndexMap[bound.rankObj.elIndex];
bound.neighObj.el = NULL;
}
}
}
+
+
+ void InteriorBoundary::serializeExcludeList(std::ostream &out, ExcludeList &list)
+ {
+ int size = list.size();
+ SerUtil::serialize(out, size);
+ for (int i = 0; i < size; i++) {
+ SerUtil::serialize(out, list[i].first);
+ SerUtil::serialize(out, list[i].second);
+ }
+ }
+
+
+ void InteriorBoundary::deserializeExcludeList(std::istream &in, ExcludeList &list)
+ {
+ int size = 0;
+ SerUtil::deserialize(in, size);
+ list.resize(0);
+ list.reserve(size);
+
+ for (int i = 0; i < size; i++) {
+ GeoIndex a;
+ int b;
+
+ SerUtil::deserialize(in, a);
+ SerUtil::deserialize(in, b);
+ list.push_back(std::make_pair(a, b));
+ }
+ }
+
}
diff --git a/AMDiS/src/InteriorBoundary.h b/AMDiS/src/InteriorBoundary.h
index 60078920..ae85f77b 100644
--- a/AMDiS/src/InteriorBoundary.h
+++ b/AMDiS/src/InteriorBoundary.h
@@ -31,12 +31,14 @@
namespace AMDiS {
+ typedef std::vector<std::pair<GeoIndex, int> > ExcludeList;
+
/// Defines the geometrical objects that forms the boundary;
struct BoundaryObject {
BoundaryObject()
: reverseMode(false),
- notIncludedSubStructures(0)
+ excludedSubstructures(0)
{}
BoundaryObject(Element *e, int eType, GeoIndex sObj, int iObj)
@@ -46,7 +48,7 @@ namespace AMDiS {
subObj(sObj),
ithObj(iObj),
reverseMode(false),
- notIncludedSubStructures(0)
+ excludedSubstructures(0)
{}
void setReverseMode(BoundaryObject &otherBound, FiniteElemSpace *feSpace);
@@ -83,7 +85,15 @@ namespace AMDiS {
bool reverseMode;
- std::vector<std::pair<GeoIndex, int> > notIncludedSubStructures;
+ /** \brief
+ * In many situations it may be necessary to exclude some parts of the element to
+ * be part of the boundary. In 3d, when a face is part of the boundary, an edge or
+ * an vertex may be exludeded. In 2d only vertices may be exluded to be part of
+ * an edge boundary. This list contains pairs of exludeded structures. The first
+ * component of every pair denotes if it is a vertex or an edge, and the second
+ * component denotes the local index of the structure.
+ */
+ ExcludeList excludedSubstructures;
};
/** \brief
@@ -117,6 +127,11 @@ namespace AMDiS {
/// Reads the state of an interior boundary from a file.
void deserialize(std::istream &in, std::map<int, Element*> &elIndexMap);
+ protected:
+ void serializeExcludeList(std::ostream &out, ExcludeList &list);
+
+ void deserializeExcludeList(std::istream &in, ExcludeList &list);
+
public:
RankToBoundMap boundary;
};
diff --git a/AMDiS/src/ParMetisPartitioner.cc b/AMDiS/src/ParMetisPartitioner.cc
index 57b98cb0..e9084ae4 100644
--- a/AMDiS/src/ParMetisPartitioner.cc
+++ b/AMDiS/src/ParMetisPartitioner.cc
@@ -287,14 +287,12 @@ namespace AMDiS {
switch(mode) {
case INITIAL:
- // TODO: Removed weight because it does not work correctly for very large
- // macro meshes.
ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
parMetisGraph.getXAdj(),
parMetisGraph.getAdjncy(),
+ wgts,
NULL,
- NULL,
- NULL,
+ &wgtflag,
&numflag,
&ncon,
&nparts,
diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc
index 23aaf48b..13bb00e4 100644
--- a/AMDiS/src/ParallelDomainBase.cc
+++ b/AMDiS/src/ParallelDomainBase.cc
@@ -1,4 +1,6 @@
#include <algorithm>
+#include <iostream>
+#include <fstream>
#include "ParallelDomainBase.h"
#include "ParMetisPartitioner.h"
@@ -49,6 +51,8 @@ namespace AMDiS {
mesh(fe->getMesh()),
refineManager(refinementManager),
initialPartitionMesh(true),
+ d_nnz(NULL),
+ o_nnz(NULL),
nRankDofs(0),
rstart(0),
nComponents(1),
@@ -364,38 +368,25 @@ namespace AMDiS {
}
- void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
+ void ParallelDomainBase::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
{
- FUNCNAME("ParallelDomainBase::fillPetscMatrix()");
-
- clock_t first = clock();
-
- // === Create PETSc vector (rhs, solution and a temporary vector). ===
+ FUNCNAME("ParallelDomainBase::createPetscNnzStructure()");
- VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
- VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
- VecSetType(petscRhsVec, VECMPI);
-
- VecCreate(PETSC_COMM_WORLD, &petscSolVec);
- VecSetSizes(petscSolVec, nRankRows, nOverallRows);
- VecSetType(petscSolVec, VECMPI);
+ TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
+ TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");
- VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
- VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
- VecSetType(petscTmpVec, VECMPI);
+ d_nnz = new int[nRankRows];
+ o_nnz = new int[nRankRows];
+ for (int i = 0; i < nRankRows; i++) {
+ d_nnz[i] = 0;
+ o_nnz[i] = 0;
+ }
using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits = mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
- int d_nnz[nRankRows];
- int o_nnz[nRankRows];
- for (int i = 0; i < nRankRows; i++) {
- d_nnz[i] = 0;
- o_nnz[i] = 0;
- }
-
// Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
// that this rank will send to. This nnz entries will be assembled on this rank,
// but because the row DOFs are not DOFs of this rank they will be send to the
@@ -531,12 +522,39 @@ namespace AMDiS {
}
}
}
+ }
+
+
+ void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
+ {
+ FUNCNAME("ParallelDomainBase::fillPetscMatrix()");
+
+ clock_t first = clock();
+
+ // === Create PETSc vector (rhs, solution and a temporary vector). ===
+
+ VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
+ VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
+ VecSetType(petscRhsVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscSolVec);
+ VecSetSizes(petscSolVec, nRankRows, nOverallRows);
+ VecSetType(petscSolVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
+ VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
+ VecSetType(petscTmpVec, VECMPI);
+
+ if (!d_nnz)
+ createPetscNnzStructure(mat);
// === Create PETSc matrix with the computed nnz data structure. ===
MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
0, d_nnz, 0, o_nnz, &petscMatrix);
+ INFO(info, 8)("Fill petsc matrix 1 needed %.5f seconds\n", TIME_USED(first, clock()));
+
#if (DEBUG != 0)
int a, b;
MatGetOwnershipRange(petscMatrix, &a, &b);
@@ -551,6 +569,8 @@ namespace AMDiS {
if ((*mat)[i][j])
setDofMatrix((*mat)[i][j], nComponents, i, j);
+ INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", TIME_USED(first, clock()));
+
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
@@ -727,10 +747,20 @@ namespace AMDiS {
do {
// To check the interior boundaries, the ownership of the boundaries is not
// important. Therefore, we add all boundaries to one boundary container.
- RankToBoundMap allBound = myIntBoundary.boundary;
+ RankToBoundMap allBound;
+ for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
+ it != myIntBoundary.boundary.end(); ++it)
+ for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
+ bit != it->second.end(); ++bit)
+ if (bit->rankObj.subObj == FACE)
+ allBound[it->first].push_back(*bit);
+
for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
it != otherIntBoundary.boundary.end(); ++it)
- allBound[it->first] = it->second;
+ for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
+ bit != it->second.end(); ++bit)
+ if (bit->rankObj.subObj == FACE)
+ allBound[it->first].push_back(*bit);
// === Check the boundaries and adapt mesh if necessary. ===
@@ -741,6 +771,7 @@ namespace AMDiS {
int sendValue = static_cast<int>(!meshChanged);
recvAllValues = 0;
mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
+ MSG("LOOP\n");
} while (recvAllValues != 0);
@@ -754,6 +785,19 @@ namespace AMDiS {
// === Because the mesh has been changed, update the DOF numbering and mappings. ===
updateLocalGlobalNumbering();
+
+ // === If there is a non zero pattern computed for Petsc matrix, delete it. So ===
+ // === it will be recomputed after next assembling. ===
+
+ if (d_nnz) {
+ delete [] d_nnz;
+ d_nnz = NULL;
+ }
+
+ if (o_nnz) {
+ delete [] o_nnz;
+ o_nnz = NULL;
+ }
}
@@ -761,6 +805,8 @@ namespace AMDiS {
{
FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()");
+ MSG("CHECK_AND_ADAPT 1!\n");
+
// === Create mesh structure codes for all ranks boundary elements. ===
std::map<int, MeshCodeVec> sendCodes;
@@ -774,14 +820,21 @@ namespace AMDiS {
}
}
+ MSG("CHECK_AND_ADAPT 2!\n");
+
StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
stdMpi.send(sendCodes);
stdMpi.recv(allBound);
stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
+ MSG("CHECK_AND_ADAPT 3!\n");
+
+ clock_t first = clock();
+
// === Compare received mesh structure codes. ===
bool meshFitTogether = true;
+ int superCounter = 0;
for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
@@ -790,18 +843,27 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
-
+
MeshStructure elCode;
elCode.init(boundIt->rankObj);
if (elCode.getCode() != recvCodes[i].getCode()) {
TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
-
- bool b = fitElementToMeshCode(refineManager,
- recvCodes[i],
+ superCounter++;
+ clock_t fff = clock();
+
+ int refCounter = 0;
+ bool b = fitElementToMeshCode(recvCodes[i],
boundIt->rankObj.el,
boundIt->rankObj.ithObj,
- boundIt->rankObj.elType);
+ boundIt->rankObj.elType, refCounter);
+
+ MSG("SIZE OF ELINFO3d = %d\n", sizeof(ElInfo3d));
+ MSG("Code length %d with %d refs needed %.5f seconds\n",
+ recvCodes[i].getNumElements(),
+ refCounter,
+ TIME_USED(fff, clock()));
+
if (b)
meshFitTogether = false;
}
@@ -810,9 +872,84 @@ namespace AMDiS {
}
}
+ MSG("time for %d needed %.5f seconds\n", superCounter, TIME_USED(first, clock()));
+ MSG("CHECK_AND_ADAPT 4!\n");
+
return meshFitTogether;
}
-
+
+
+ bool ParallelDomainBase::fitElementToMeshCode(MeshStructure &code,
+ Element *el,
+ int ithSide,
+ int elType,
+ int &refCounter)
+ {
+ FUNCNAME("fitElementToMeshCode()");
+
+ if (code.empty())
+ return false;
+
+ int s1 = el->getSideOfChild(0, ithSide, elType);
+ int s2 = el->getSideOfChild(1, ithSide, elType);
+
+ TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");
+
+ if (s1 != -1 && s2 != -1)
+ return fitElementToMeshCode2(code, el, ithSide, elType, refCounter);
+
+ if (el->isLeaf()) {
+ if (code.getNumElements() == 1 && code.isLeafElement())
+ return false;
+
+ el->setMark(1);
+ refineManager->refineElement(el->getMesh(), el);
+ refCounter++;
+ }
+
+ if (s1 != -1)
+ return fitElementToMeshCode2(code,
+ el->getFirstChild(), s1, el->getChildType(elType), refCounter);
+ else
+ return fitElementToMeshCode2(code,
+ el->getSecondChild(), s2, el->getChildType(elType), refCounter);
+ }
+
+ bool ParallelDomainBase::fitElementToMeshCode2(MeshStructure &code,
+ Element *el,
+ int ithSide,
+ int elType, int &refCounter)
+ {
+ if (code.isLeafElement())
+ return false;
+
+ bool value = false;
+
+ if (el->isLeaf()) {
+ el->setMark(1);
+ refineManager->refineElement(el->getMesh(), el);
+ value = true;
+ refCounter++;
+ }
+
+ int s1 = el->getSideOfChild(0, ithSide, elType);
+ int s2 = el->getSideOfChild(1, ithSide, elType);
+
+ if (s1 != -1) {
+ code.nextElement();
+ value |= fitElementToMeshCode2(code, el->getFirstChild(),
+ s1, el->getChildType(elType), refCounter);
+ }
+
+ if (s2 != -1) {
+ code.nextElement();
+ value |= fitElementToMeshCode2(code, el->getSecondChild(),
+ s2, el->getChildType(elType), refCounter);
+ }
+
+ return value;
+ }
+
void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
{
@@ -872,32 +1009,55 @@ namespace AMDiS {
double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo)
{
- double localWeightSum = 0.0;
- int elNum = -1;
+ FUNCNAME("ParallelDomainBase::setElemWeights()");
+ double localWeightSum = 0.0;
elemWeights.clear();
+
+ std::string filename = "";
+ GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
+ if (filename != "") {
+ MSG("Read macro weights from %s\n", filename.c_str());
+
+ std::ifstream infile;
+ infile.open(filename.c_str(), std::ifstream::in);
+ while (!infile.eof()) {
+ int elNum, elWeight;
+ infile >> elNum;
+ if (infile.eof())
+ break;
+ infile >> elWeight;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- Element *element = elInfo->getElement();
-
- // get partition data
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
+ elemWeights[elNum] = elWeight;
+ localWeightSum += elWeight;
+ }
- if (partitionData && partitionData->getPartitionStatus() == IN) {
- if (partitionData->getLevel() == 0)
- elNum = element->getIndex();
+ infile.close();
+ } else {
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ Element *element = elInfo->getElement();
+
+ // get partition data
+ PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
+ (element->getElementData(PARTITION_ED));
- TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
- if (element->isLeaf()) {
- elemWeights[elNum] += 1.0;
- localWeightSum += 1.0;
+ if (partitionData && partitionData->getPartitionStatus() == IN) {
+ int elNum = -1;
+ if (partitionData->getLevel() == 0)
+ elNum = element->getIndex();
+
+ TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
+ if (element->isLeaf()) {
+ elemWeights[elNum] += 1.0;
+ localWeightSum += 1.0;
+ }
}
+
+ elInfo = stack.traverseNext(elInfo);
}
-
- elInfo = stack.traverseNext(elInfo);
}
return localWeightSum;
@@ -1229,7 +1389,7 @@ namespace AMDiS {
// Otherwise, it is part of the interior boundary and we add it to the set
// rankBoundaryEdges.
if (edgeOwner[edge] > mpiRank)
- boundIt->rankObj.notIncludedSubStructures.push_back(std::pair<GeoIndex, int>(EDGE, edgeNo));
+ boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
else
rankBoundaryEdges.insert(edge);
}
@@ -1240,7 +1400,7 @@ namespace AMDiS {
DegreeOfFreedom dof = localIndices[dofNo];
if (dofOwner[dof] > mpiRank)
- boundIt->rankObj.notIncludedSubStructures.push_back(std::pair<GeoIndex, int>(VERTEX, dofNo));
+ boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
else
rankBoundaryDofs.insert(dof);
}
@@ -1262,7 +1422,7 @@ namespace AMDiS {
GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));
if (edgeOwner[edge] > rankIt->first)
- boundIt->rankObj.notIncludedSubStructures.push_back(std::pair<GeoIndex, int>(EDGE, edgeNo));
+ boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
else
rankBoundaryEdges.insert(edge);
@@ -1273,7 +1433,7 @@ namespace AMDiS {
DegreeOfFreedom dof = localIndices[dofNo];
if (dofOwner[dof] > rankIt->first)
- boundIt->rankObj.notIncludedSubStructures.push_back(std::pair<GeoIndex, int>(VERTEX, dofNo));
+ boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
else
rankBoundaryDofs.insert(dof);
}
@@ -1816,8 +1976,7 @@ namespace AMDiS {
// === Send new DOF indices. ===
-#if 0
- StdMpi<std::vector<DegreeOfFreedom> > stdMpi(mpiComm);
+ StdMpi<std::vector<DegreeOfFreedom> > stdMpi(mpiComm, false);
for (RankToDofContainer::iterator sendIt = sendDofs.begin();
sendIt != sendDofs.end(); ++sendIt, i++) {
stdMpi.getSendData(sendIt->first).resize(0);
@@ -1826,56 +1985,18 @@ namespace AMDiS {
dofIt != sendIt->second.end(); ++dofIt)
stdMpi.getSendData(sendIt->first).push_back(rankDofsNewGlobalIndex[*dofIt]);
}
-#endif
-
- std::vector<int*> sendBuffers(sendDofs.size());
- std::vector<int*> recvBuffers(recvDofs.size());
-
- MPI::Request request[sendDofs.size() + recvDofs.size()];
- int requestCounter = 0;
-
- i = 0;
- for (RankToDofContainer::iterator sendIt = sendDofs.begin();
- sendIt != sendDofs.end(); ++sendIt, i++) {
- int nSendDofs = sendIt->second.size();
- sendBuffers[i] = new int[nSendDofs];
- int c = 0;
- for (DofContainer::iterator dofIt = sendIt->second.begin();
- dofIt != sendIt->second.end(); ++dofIt)
- sendBuffers[i][c++] = rankDofsNewGlobalIndex[*dofIt];
-
- request[requestCounter++] =
- mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
- }
-
- i = 0;
- for (RankToDofContainer::iterator recvIt = recvDofs.begin();
- recvIt != recvDofs.end(); ++recvIt, i++) {
- int nRecvDofs = recvIt->second.size();
- recvBuffers[i] = new int[nRecvDofs];
-
- request[requestCounter++] =
- mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
- }
-
- MPI::Request::Waitall(requestCounter, request);
-
- for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
- delete [] sendBuffers[j];
+ stdMpi.updateSendDataSize();
+ stdMpi.recv(recvDofs);
+ stdMpi.startCommunication<int>(MPI_INT);
- i = 0;
for (RankToDofContainer::iterator recvIt = recvDofs.begin();
recvIt != recvDofs.end(); ++recvIt) {
int j = 0;
for (DofContainer::iterator dofIt = recvIt->second.begin();
dofIt != recvIt->second.end(); ++dofIt) {
-
- rankDofsNewGlobalIndex[*dofIt] = recvBuffers[i][j];
+ rankDofsNewGlobalIndex[*dofIt] = stdMpi.getRecvData(recvIt->first)[j++];
isRankDof[rankDofsNewLocalIndex[*dofIt]] = false;
- j++;
}
-
- delete [] recvBuffers[i++];
}
// === Create now the local to global index and local to dof index mappings. ===
@@ -2017,23 +2138,18 @@ namespace AMDiS {
// Clear all periodic dof mappings calculated before. We do it from scratch.
periodicDof.clear();
- MPI::Request request[min(static_cast<int>(periodicBoundary.boundary.size() * 2), 4)];
- int requestCounter = 0;
- std::vector<int*> sendBuffers, recvBuffers;
+ StdMpi<std::vector<int> > stdMpi(mpiComm, false);
// === Each rank traverse its periodic boundaries and sends the dof indices ===
// === to the rank "on the other side" of the periodic boundary. ===
for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
it != periodicBoundary.boundary.end(); ++it) {
-
TEST_EXIT_DBG(it->first != mpiRank)
("Periodic interior boundary within the rank itself is not possible!\n");
- DofContainer dofs;
-
// Create dof indices on the boundary.
-
+ DofContainer dofs;
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs);
@@ -2041,27 +2157,16 @@ namespace AMDiS {
}
// Send the global indices to the rank on the other side.
-
- int *sendbuf = new int[dofs.size()];
+ stdMpi.getSendData(it->first).reserve(dofs.size());
for (int i = 0; i < static_cast<int>(dofs.size()); i++)
- sendbuf[i] = mapLocalGlobalDofs[*(dofs[i])];
-
- request[requestCounter++] =
- mpiComm.Isend(sendbuf, dofs.size(), MPI_INT, it->first, 0);
-
- sendBuffers.push_back(sendbuf);
+ stdMpi.getSendData(it->first).push_back(mapLocalGlobalDofs[*(dofs[i])]);
// Receive from this rank the same number of dofs.
-
- int *recvbuf = new int[dofs.size()];
- request[requestCounter++] =
- mpiComm.Irecv(recvbuf, dofs.size(), MPI_INT, it->first, 0);
-
- recvBuffers.push_back(recvbuf);
+ stdMpi.recv(it->first, dofs.size());
}
-
- MPI::Request::Waitall(requestCounter, request);
+ stdMpi.updateSendDataSize();
+ stdMpi.startCommunication<int>(MPI_INT);
// === The rank has received the dofs from the rank on the other side of ===
// === the boundary. Now it can use them to create the mapping between ===
@@ -2070,8 +2175,6 @@ namespace AMDiS {
std::map<DegreeOfFreedom, std::set<int> > dofFromRank;
-
- int i = 0;
for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
it != periodicBoundary.boundary.end(); ++it) {
DofContainer dofs;
@@ -2090,23 +2193,19 @@ namespace AMDiS {
// Added the received dofs to the mapping.
for (int j = 0; j < static_cast<int>(dofs.size()); j++) {
int globalDofIndex = mapLocalGlobalDofs[*(dofs[j])];
- periodicDof[globalDofIndex].insert(recvBuffers[i][j]);
+ periodicDof[globalDofIndex].insert(stdMpi.getRecvData(it->first)[j]);
dofFromRank[globalDofIndex].insert(it->first);
}
-
- delete [] sendBuffers[i];
- delete [] recvBuffers[i];
- i++;
}
- sendBuffers.clear();
- recvBuffers.clear();
-
if (dofFromRank.size() > 0)
TEST_EXIT_DBG(mesh->getDim() == 2)
("Periodic boundary corner problem must be generalized to 3d!\n");
- requestCounter = 0;
+ MPI::Request request[min(static_cast<int>(periodicBoundary.boundary.size() * 2), 4)];
+ int requestCounter = 0;
+ std::vector<int*> sendBuffers, recvBuffers;
+
for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin();
it != dofFromRank.end(); ++it) {
if (it->second.size() == 2) {
@@ -2138,7 +2237,7 @@ namespace AMDiS {
MPI::Request::Waitall(requestCounter, request);
- i = 0;
+ int i = 0;
for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin();
it != dofFromRank.end(); ++it) {
if (it->second.size() == 2) {
diff --git a/AMDiS/src/ParallelDomainBase.h b/AMDiS/src/ParallelDomainBase.h
index d04c90a1..32684301 100644
--- a/AMDiS/src/ParallelDomainBase.h
+++ b/AMDiS/src/ParallelDomainBase.h
@@ -101,6 +101,15 @@ namespace AMDiS {
virtual void exitParallelization(AdaptInfo *adaptInfo);
+ /** \brief
+ * This function checks if the mesh has changed on at least on rank. In this case,
+ * the interior boundaries are adapted on all ranks such that they fit together on
+ * all ranks. Furthermore the function \ref updateLocalGlobalNumbering() is called
+ * to update the dof numberings and mappings on all rank due to the new mesh
+ * structure.
+ */
+ void checkMeshChange();
+
void updateDofAdmins();
/** \brief
@@ -184,7 +193,6 @@ namespace AMDiS {
{
return NULL;
}
-
// Writes all data of this object to an output stream.
virtual void serialize(std::ostream &out);
@@ -265,6 +273,9 @@ namespace AMDiS {
DofToRank& boundaryDofs,
DofToBool& vertexDof);
+ /// Creates a new non zero pattern structure for the Petsc matrix.
+ void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
+
/** \brief
* Create a PETSc matrix and PETSc vectors. The given DOF matrices are used to
* create the nnz structure of the PETSc matrix and the values are transfered to it.
@@ -283,15 +294,6 @@ namespace AMDiS {
void setDofVector(Vec& petscVec, DOFVector<double>* vec,
int disMult = 1, int dispAdd = 0);
- /** \brief
- * This function checks if the mesh has changed on at least on rank. In this case,
- * the interior boundaries are adapted on all ranks such that they fit together on
- * all ranks. Furthermore the function \ref updateLocalGlobalNumbering() is called
- * to update the dof numberings and mappings on all rank due to the new mesh
- * structure.
- */
- void checkMeshChange();
-
/** \brief
* Checks for all given interior boundaries if the elements fit together on both
* sides of the boundaries. If this is not the case, the mesh is adapted. Because
@@ -376,6 +378,18 @@ namespace AMDiS {
*/
void synchVector(SystemVector &vec);
+ bool fitElementToMeshCode(MeshStructure &code,
+ Element *el,
+ int ithSide,
+ int elType,
+ int &refCounter);
+
+ bool fitElementToMeshCode2(MeshStructure &code,
+ Element *el,
+ int ithSide,
+ int elType,
+ int &refCounter);
+
/// Writes a vector of dof pointers to an output stream.
void serialize(std::ostream &out, DofContainer &data);
@@ -536,6 +550,9 @@ namespace AMDiS {
* temporary vector for calculating the final residuum.
*/
Vec petscRhsVec, petscSolVec, petscTmpVec;
+
+ /// Arrays definig the non zero pattern of Petsc's matrix.
+ int *d_nnz, *o_nnz;
/// Number of DOFs in the rank mesh.
int nRankDofs;
diff --git a/AMDiS/src/StdMpi.h b/AMDiS/src/StdMpi.h
index 41d72e13..56bd928a 100644
--- a/AMDiS/src/StdMpi.h
+++ b/AMDiS/src/StdMpi.h
@@ -297,6 +297,13 @@ namespace AMDiS {
return sendData[rank];
}
+ void updateSendDataSize()
+ {
+ for (typename std::map<int, SendT>::iterator it = sendData.begin();
+ it != sendData.end(); ++it)
+ sendDataSize[it->first] = intSizeOf(it->second);
+ }
+
void recv(int fromRank, int size = -1)
{
FUNCNAME("StdMpi::recv()");
diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc
index 8a34ec25..8899c1f2 100644
--- a/AMDiS/src/Tetrahedron.cc
+++ b/AMDiS/src/Tetrahedron.cc
@@ -215,7 +215,7 @@ namespace AMDiS {
{
if (parentVertices) {
- ERROR_EXIT("Einbau notIncludedSubStructures!\n");
+ ERROR_EXIT("Einbau excludedSubstructures!\n");
switch (bound.ithObj) {
case 0:
@@ -264,8 +264,8 @@ namespace AMDiS {
nextBound1.reverseMode = false;
bool addDof = true;
- for (std::vector<std::pair<GeoIndex, int> >::iterator it = bound.notIncludedSubStructures.begin();
- it != bound.notIncludedSubStructures.end(); ++it)
+ for (ExcludeList::iterator it = bound.excludedSubstructures.begin();
+ it != bound.excludedSubstructures.end(); ++it)
if (it->first == EDGE && it->second == 0)
addDof = false;
@@ -300,7 +300,7 @@ namespace AMDiS {
FUNCNAME("Tetrahedron::getVertexDofsAtEdge()");
if (parentVertices) {
- ERROR_EXIT("Einbau notIncludedSubStructures!\n");
+ ERROR_EXIT("Einbau excludedSubstructures!\n");
}
if (!child[0])
@@ -330,8 +330,8 @@ namespace AMDiS {
switch (bound.subObj) {
case FACE:
- for (std::vector<std::pair<GeoIndex, int> >::iterator it = bound.notIncludedSubStructures.begin();
- it != bound.notIncludedSubStructures.end(); ++it) {
+ for (ExcludeList::iterator it = bound.excludedSubstructures.begin();
+ it != bound.excludedSubstructures.end(); ++it) {
if (it->first == EDGE)
it->second = edgeOfChild[bound.elType][ithChild][it->second];
}
--
GitLab