Commit 7908ff20 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

3d adaptivity for parallelization.

parent 82dffb9a
......@@ -154,7 +154,9 @@ namespace AMDiS {
}
void getVertexDofs(FiniteElemSpace* feSpace, int ithEdge, int elType,
DofContainer& dofs, bool parentVertices = 0) const
DofContainer& dofs,
bool reverseMode,
bool parentVertices = false) const
{
FUNCNAME("Line::getVertexDofs()");
......
......@@ -28,7 +28,7 @@ namespace AMDiS {
MacroElement::~MacroElement()
{
if (element)
delete element;
delete element;
}
MacroElement& MacroElement::operator=(const MacroElement &el)
......
......@@ -262,7 +262,7 @@ namespace AMDiS {
me->setIndex(macroElements.size());
}
void Mesh::removeMacroElements(std::vector<MacroElement*>& macros,
void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
const FiniteElemSpace *feSpace)
{
FUNCNAME("Mesh::removeMacroElement()");
......@@ -276,11 +276,13 @@ namespace AMDiS {
ElementDofIterator elDofIter(feSpace);
// Map that stores for each dof pointer (which may have a list of dofs)
// all macro element indices that own the dof.
// all macro element indices that own this dof.
DofElMap dofsOwner;
DofPosMap dofsPosIndex;
// Determine all dof owner macro elements.
// === Determine to all dofs the macro elements where the dof is part of. ===
for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin();
macroIt != macroElements.end(); ++macroIt) {
elDofIter.reset((*macroIt)->getElement());
......@@ -290,16 +292,31 @@ namespace AMDiS {
} while (elDofIter.nextStrict());
}
// Remove all the given macro elements.
for (std::vector<MacroElement*>::iterator macroIt = macros.begin();
macroIt != macros.end();
++macroIt) {
// === Remove macro elements from mesh macro element list. ===
std::deque<MacroElement*>::iterator mEl =
find(macroElements.begin(), macroElements.end(), *macroIt);
TEST_EXIT(mEl != macroElements.end())
("Cannot find MacroElement that should be removed!\n");
macroElements.erase(mEl);
// Removing arbitrary elements from an std::deque is very slow. Therefore, we
// create a new deque with all macro elements that should not be deleted. The
// macro element deque is than replaced by the new created one.
std::deque<MacroElement*> newMacroElements;
for (std::deque<MacroElement*>::iterator elIter = macroElements.begin();
elIter != macroElements.end(); ++elIter) {
// If the current mesh macro element should not be deleted, i.e., it is not a
// member of the list of macro elements to be deleted, is is inserted to the new
// macro element list.
if (macros.find(*elIter) == macros.end())
newMacroElements.push_back(*elIter);
}
// And replace the macro element list with the new one.
macroElements.clear();
macroElements = newMacroElements;
// === For all macro elements to be deleted, delete them also to be neighbours ===
// === of some other macro elements. ===
for (std::set<MacroElement*>::iterator macroIt = macros.begin();
macroIt != macros.end(); ++macroIt) {
// Go through all neighbours of the macro element and remove this macro element
// to be neighbour of some other macro element.
......@@ -315,33 +332,44 @@ namespace AMDiS {
}
}
// Decrease also the number of elements of the mesh.
nLeaves--;
nElements--;
}
// Remove this macro element from the dof owner list.
for (DofElMap::iterator dofsIt = dofsOwner.begin();
dofsIt != dofsOwner.end(); ++dofsIt) {
std::set<MacroElement*>::iterator mIt = dofsIt->second.find(*macroIt);
if (mIt != dofsIt->second.end())
dofsIt->second.erase(mIt);
}
// And remove the macro element from memory
delete *macroIt;
}
// === Check now all the dofs, that have no owner anymore and therefore have ===
// === to be removed. ===
int nRemainDofs = 0;
// Check now all the dofs, that have no owner anymore and therefore have to
// be removed.
for (DofElMap::iterator dofsIt = dofsOwner.begin();
dofsIt != dofsOwner.end(); ++dofsIt) {
if (dofsIt->second.size() == 0)
dofsIt != dofsOwner.end(); ++dofsIt) {
bool deleteDof = true;
for (std::set<MacroElement*>::iterator elIter = dofsIt->second.begin();
elIter != dofsIt->second.end(); ++elIter) {
std::set<MacroElement*>::iterator mIt = macros.find(*elIter);
if (mIt == macros.end()) {
deleteDof = false;
break;
}
}
if (deleteDof)
freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first),
dofsPosIndex[dofsIt->first]);
else
nRemainDofs++;
}
// === Finally, remove the macro elements from memory. ===
for (std::set<MacroElement*>::iterator macroIt = macros.begin();
macroIt != macros.end(); ++macroIt)
delete *macroIt;
nVertices = nRemainDofs;
}
......@@ -784,9 +812,9 @@ namespace AMDiS {
int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
ichild, c_outside, c_lambda[0],c_lambda[1],c_lambda[2],c_lambda[3]);
ichild, c_outside, c_lambda[0], c_lambda[1], c_lambda[2], c_lambda[3]);
MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
1-ichild, c_outside2, c_lambda2[0],c_lambda2[1],c_lambda2[2],
1 - ichild, c_outside2, c_lambda2[0], c_lambda2[1], c_lambda2[2],
c_lambda2[3]);
if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
......@@ -830,23 +858,24 @@ namespace AMDiS {
return inside;
}
bool Mesh::getDofIndexCoords(const DegreeOfFreedom* dof,
bool Mesh::getDofIndexCoords(DegreeOfFreedom dof,
const FiniteElemSpace* feSpace,
WorldVector<double>& coords)
{
DimVec<double>* baryCoords;
bool found = false;
TraverseStack stack;
Vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
std::vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
ElInfo *elInfo = stack.traverseFirst(this, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
feSpace->getAdmin(),
&dofVec);
feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
feSpace->getAdmin(),
dofVec);
for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
if (dofVec[i] == *dof) {
if (dofVec[i] == dof) {
baryCoords = feSpace->getBasisFcts()->getCoords(i);
elInfo->coordToWorld(*baryCoords, coords);
found = true;
......@@ -863,39 +892,30 @@ namespace AMDiS {
return found;
}
bool Mesh::getDofIndexCoords(DegreeOfFreedom dof,
const FiniteElemSpace* feSpace,
WorldVector<double>& coords)
void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace,
DOFVector<WorldVector<double> >& coords)
{
DimVec<double>* baryCoords;
bool found = false;
TraverseStack stack;
Vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
const BasisFunction* basFcts = feSpace->getBasisFcts();
int nBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> dofVec(nBasFcts);
ElInfo *elInfo = stack.traverseFirst(this, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
TraverseStack stack;
ElInfo *elInfo =
stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
feSpace->getAdmin(),
&dofVec);
for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
if (dofVec[i] == dof) {
baryCoords = feSpace->getBasisFcts()->getCoords(i);
elInfo->coordToWorld(*baryCoords, coords);
found = true;
break;
}
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
for (int i = 0; i < nBasFcts; i++) {
DimVec<double> *baryCoords = basFcts->getCoords(i);
elInfo->coordToWorld(*baryCoords, coords[dofVec[i]]);
}
if (found)
break;
elInfo = stack.traverseNext(elInfo);
}
return found;
}
void Mesh::setDiameter(const WorldVector<double>& w)
{
diam = w;
......
......@@ -414,7 +414,7 @@ namespace AMDiS {
* that there are no global or local refinements, i.e., all macro elements have
* no children.
*/
void removeMacroElements(std::vector<MacroElement*>& macros,
void removeMacroElements(std::set<MacroElement*>& macros,
const FiniteElemSpace* feSpace);
/// Frees the array of DOF pointers (see \ref createDOFPtrs)
......@@ -485,12 +485,31 @@ namespace AMDiS {
*/
bool getDofIndexCoords(const DegreeOfFreedom* dof,
const FiniteElemSpace* feSpace,
WorldVector<double>& coords);
WorldVector<double>& coords)
{
return getDofIndexCoords(*dof, feSpace, coords);
}
/** \brief
* This function is equal to \ref getDofIndexCoords as defined above, but takes
* a DOF index instead of a DOF pointer.
*/
bool getDofIndexCoords(DegreeOfFreedom dof,
const FiniteElemSpace* feSpace,
WorldVector<double>& coords);
/** \brief
* Traverse the whole mesh and stores to each DOF of the given finite element space
* the coordinates in a given DOFVector. Works in the same way as the function
* \ref getDofIndexCoords defined above.
*
* @param[in] feSpace The fe soace to be used for the search.
* @param[out] coords DOF vector that stores the coordinates to each dof.
*/
void getDofIndexCoords(const FiniteElemSpace* feSpace,
DOFVector<WorldVector<double> >& coords);
/// Returns FILL_ANY_?D
inline static const Flag& getFillAnyFlag(int dim)
{
......
......@@ -6,6 +6,7 @@
#include "Traverse.h"
#include "ElInfo.h"
#include "RefinementManager.h"
#include "InteriorBoundary.h"
namespace AMDiS {
......@@ -59,7 +60,46 @@ namespace AMDiS {
clear();
addAlongSide(el, ithSide, elType, reverseOrder);
int s1 = el->getSideOfChild(0, ithSide, elType);
int s2 = el->getSideOfChild(1, ithSide, elType);
TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
if (!el->isLeaf()) {
if (s1 == -1)
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder);
else if (s2 == -1)
addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder);
else
addAlongSide(el, ithSide, elType, reverseOrder);
}
commit();
}
void MeshStructure::init(BoundaryObject &bound)
{
FUNCNAME("MeshStructure::init()");
Element *el = bound.el;
clear();
int s1 = el->getSideOfChild(0, bound.ithObj, bound.elType);
int s2 = el->getSideOfChild(1, bound.ithObj, bound.elType);
TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
if (!el->isLeaf()) {
if (s1 == -1)
addAlongSide(el->getSecondChild(), s2,
el->getChildType(bound.elType), bound.reverseMode);
else if (s2 == -1)
addAlongSide(el->getFirstChild(), s1,
el->getChildType(bound.elType), bound.reverseMode);
else
addAlongSide(el, bound.ithObj, bound.elType, bound.reverseMode);
}
commit();
}
......@@ -67,6 +107,14 @@ namespace AMDiS {
void MeshStructure::addAlongSide(Element *el, int ithSide, int elType,
bool reverseOrder)
{
FUNCNAME("MeshStructure::addAlongSide()");
if (debugMode) {
MSG("addAlondSide(%d, %d, %d, %d)\n",
el->getIndex(), ithSide, elType, reverseOrder);
MSG("Element is leaf: %d\n", el->isLeaf());
}
insertElement(el->isLeaf());
if (!el->isLeaf()) {
......@@ -75,14 +123,14 @@ namespace AMDiS {
if (!reverseOrder) {
if (s1 != -1)
addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder);
addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), false);
if (s2 != -1)
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder);
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), false);
} else {
if (s2 != -1)
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder);
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), false);
if (s1 != -1)
addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder);
addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), false);
}
}
}
......@@ -90,9 +138,13 @@ namespace AMDiS {
void MeshStructure::reset()
{
currentIndex = 0;
currentCode = code[0];
pos = 0;
currentElement = 0;
if (code.size() > 0)
currentCode = code[0];
else
currentCode = 0;
}
bool MeshStructure::nextElement(MeshStructure *insert)
......@@ -238,4 +290,8 @@ namespace AMDiS {
} while (!finished);
}
bool MeshStructure::compare(MeshStructure &other)
{
return (other.getCode() == code);
}
}
......@@ -36,7 +36,8 @@ namespace AMDiS {
currentCode(0),
pos(0),
currentElement(0),
nElements(0)
nElements(0),
debugMode(false)
{}
void clear();
......@@ -46,6 +47,8 @@ namespace AMDiS {
void init(Element *el, int ithSide, int elType, bool reverseOrder);
void init(BoundaryObject &bound);
void init(const std::vector<unsigned long int>& initCode, int n)
{
code = initCode;
......@@ -59,6 +62,12 @@ namespace AMDiS {
*/
void reset();
/// Returns whether the code is empty or not.
inline bool empty()
{
return (nElements == 0);
}
inline void commit()
{
if (pos > 0)
......@@ -104,6 +113,11 @@ namespace AMDiS {
{
FUNCNAME("MeshStructure::print()");
if (empty()) {
std::cout << "-" << std::endl;
return;
}
reset();
bool cont = true;
while (cont) {
......@@ -133,6 +147,14 @@ namespace AMDiS {
return currentElement;
}
void setDebugMode(bool b)
{
debugMode = b;
}
/// Returns true, if the given mesh structure code is equal to this one.
bool compare(MeshStructure &other);
protected:
/// Insert a new element to the structure code. Is used by the init function.
void insertElement(bool isLeaf);
......@@ -157,6 +179,8 @@ namespace AMDiS {
int nElements;
bool debugMode;
static const int unsignedLongSize;
};
......
......@@ -24,8 +24,7 @@ namespace AMDiS {
int dow = Global::getGeo(WORLD);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_EVERY_EL_PREORDER);
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
while (elInfo) {
// get partition data
PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
......@@ -33,9 +32,8 @@ namespace AMDiS {
if (partitionData &&
partitionData->getPartitionStatus() == IN &&
partitionData->getLevel() == 0) {
elementCounter++;
}
partitionData->getLevel() == 0)
elementCounter++;
elInfo = stack.traverseNext(elInfo);
}
......@@ -75,8 +73,7 @@ namespace AMDiS {
elementCounter = 0;
elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_EVERY_EL_PREORDER |
Mesh::FILL_COORDS);
Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_COORDS);
while (elInfo) {
Element *element = elInfo->getElement();
int index = element->getIndex();
......@@ -169,8 +166,7 @@ namespace AMDiS {
void ParMetisPartitioner::deletePartitionData()
{
TraverseStack stack;
ElInfo *elInfo;
elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER);
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
while (elInfo) {
Element *element = elInfo->getElement();
element->deleteElementData(PARTITION_ED);
......@@ -178,16 +174,17 @@ namespace AMDiS {
}
}
void ParMetisPartitioner::createPartitionData() {
void ParMetisPartitioner::createPartitionData()
{
int mpiRank = mpiComm->Get_rank();
int mpiSize = mpiComm->Get_size();
int nLeaves = mesh->getNumberOfLeaves();
int elPerRank = nLeaves / mpiSize;
TraverseStack stack;
ElInfo *elInfo;
// === Create initial partitioning of the AMDiS mesh. ===
// === create initial partitioning on AMDiS mesh ===
int totalElementCounter = 0;
elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_LEAF_EL);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
Element *element = elInfo->getElement();
......@@ -198,13 +195,12 @@ namespace AMDiS {
new PartitionElementData(element->getElementData());
element->setElementData(elData);
if (totalElementCounter % mpiSize == mpiRank)
if (element->getIndex() >= mpiRank * elPerRank &&
element->getIndex() < (mpiRank + 1) * elPerRank)
elData->setPartitionStatus(IN);
else
elData->setPartitionStatus(UNDEFINED);
totalElementCounter++;
elInfo = stack.traverseNext(elInfo);
}
}
......@@ -213,15 +209,15 @@ namespace AMDiS {
PartitionMode mode,
float itr)
{
int mpiSize = mpiComm->Get_size();
FUNCNAME("ParMetisPartitioner::partition()");
TraverseStack stack;
ElInfo *elInfo;
int mpiSize = mpiComm->Get_size();
// === create parmetis mesh ===
if (parMetisMesh)
delete parMetisMesh;
parMetisMesh = new ParMetisMesh(mesh_, mpiComm);
parMetisMesh = new ParMetisMesh(mesh, mpiComm);
int nElements = parMetisMesh->getNumElements();
......@@ -231,7 +227,8 @@ namespace AMDiS {
float maxWgt = 0.0;
float *ptr_floatWgts = floatWgts;
elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
while (elInfo) {
Element *element = elInfo->getElement();
......@@ -242,6 +239,7 @@ namespace AMDiS {
if (partitionData &&
partitionData->getPartitionStatus() == IN &&
partitionData->getLevel() == 0) {
int index = element->getIndex();
// get weight
......@@ -267,6 +265,7 @@ namespace AMDiS {
int numflag = 0; // c numbering style!
int ncon = elemWeights ? 1 : 0; // one weight at each vertex!
int nparts = mpiSize; // number of partitions
float *tpwgts = elemWeights ? new float[mpiSize] : NULL;
float ubvec = 1.05;
int options[3] = {0, 0, 0}; // default options
......@@ -278,7 +277,7 @@ namespace AMDiS {
for (int i = 0; i < mpiSize; i++)
tpwgts[i] = 1.0 / nparts;
float scale = 10000 / maxWgt;
float scale = 10000.0 / maxWgt;
// scale wgts
for (int i = 0; i < nElements; i++)
wgts[i] = static_cast<int>(floatWgts[i] * scale);
......@@ -288,12 +287,14 @@ 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,
&wgtflag,
NULL,
NULL,
&numflag,
&ncon,
&nparts,
......@@ -370,8 +371,8 @@ namespace AMDiS {
partitionVec->clear();