Commit 3222f37e authored by Nestler, Michael's avatar Nestler, Michael
Browse files

Anpassung in Tetrahedron.cc in getHigherOrderDofs

Quadratische Ansatzfunktionen erste Tests bestanden

In ParallelDebug.cc einen SVN Writer für DOFVektoren hinzugefügt
parent 66434cfd
......@@ -227,8 +227,9 @@ namespace AMDiS {
{
FUNCNAME("Tetrahedron::getNodeDofsAtFace()");
if (!child[0])
if (!child[0]){
return;
}
switch (bound.ithObj) {
case 0:
......@@ -292,8 +293,10 @@ namespace AMDiS {
{
FUNCNAME("Tetrahedron::getNodeDofsAtEdge()");
if (!child[0])
return;
if (!child[0]){
return;
}
int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX));
BoundaryObject nextBound0 = bound, nextBound1 = bound;
......@@ -306,13 +309,15 @@ namespace AMDiS {
nextBound1.reverseMode = false;
if (bound.reverseMode) {
child[1]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[0]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
child[1]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[0]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
} else {
child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
}
break;
......@@ -333,6 +338,7 @@ namespace AMDiS {
if (nextBound1.ithObj != -1)
child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
}
}
......@@ -365,21 +371,30 @@ namespace AMDiS {
// element is also refined. Then we have go down further in refinement
// hierarchie.
if (bound.reverseMode) {
if (nextBound1.ithObj >= 0)
child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
baseDofPtr, dofGeoIndex);
if (nextBound0.ithObj >= 0)
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
baseDofPtr, dofGeoIndex);
// Additonal Check if Boundary Edge is NOT Refinement Edge of BoundaryObject
// AND the complete Edge is part of both childs
// therefore only recursion on ONE child is necessary
if ( (bound.ithObj > 0) && (nextBound0.ithObj >= 0 && nextBound1.ithObj >= 0) ){
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
baseDofPtr, dofGeoIndex);
} else {
if (nextBound0.ithObj >= 0)
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
baseDofPtr, dofGeoIndex);
if (nextBound1.ithObj >= 0)
child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
baseDofPtr, dofGeoIndex);
if (bound.reverseMode) {
if (nextBound1.ithObj >= 0)
child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
baseDofPtr, dofGeoIndex);
if (nextBound0.ithObj >= 0)
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
baseDofPtr, dofGeoIndex);
} else {
if (nextBound0.ithObj >= 0)
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
baseDofPtr, dofGeoIndex);
if (nextBound1.ithObj >= 0)
child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
baseDofPtr, dofGeoIndex);
}
}
} else {
// Either the edge is not contained in further refined children, or
// the element is not refined further on this edge. Then we can get
......@@ -401,9 +416,10 @@ namespace AMDiS {
do {
if (elDofIter.getCurrentPos() == 1 &&
elDofIter.getCurrentElementPos() == bound.ithObj) {
dofs.push_back(elDofIter.getDofPtr());
if (dofGeoIndex != NULL)
dofGeoIndex->push_back(EDGE);
dofs.push_back(elDofIter.getDofPtr());
if (dofGeoIndex != NULL)
dofGeoIndex->push_back(EDGE);
}
} while (elDofIter.next());
}
......@@ -474,8 +490,11 @@ namespace AMDiS {
addDof = true;
else
addDof = false;
break;
case 3:
// On CENTER, dof can not be part on internal boundary
addDof = false;
break;
default:
ERROR_EXIT("Should not happen!\n");
}
......@@ -483,12 +502,13 @@ namespace AMDiS {
if (addDof) {
if (baseDofPtr)
dofs.push_back(elDofIter.getBaseDof());
dofs.push_back(elDofIter.getBaseDof());
else
dofs.push_back(elDofIter.getDofPtr());
dofs.push_back(elDofIter.getDofPtr());
if (dofGeoIndex != NULL)
dofGeoIndex->push_back(elDofIter.getPosIndex());
dofGeoIndex->push_back(elDofIter.getPosIndex());
}
next = (baseDofPtr ? elDofIter.nextStrict() : elDofIter.next());
......@@ -499,6 +519,7 @@ namespace AMDiS {
default:
ERROR_EXIT("Should not happen!\n");
}
}
......
......@@ -63,31 +63,32 @@ namespace AMDiS {
for (unsigned int i = 0; i < feSpaces.size(); i++)
for (int level = 0; level < nLevel; level++)
for (InteriorBoundary::iterator it(boundary, level); !it.end(); ++it)
it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj,
data[level][it.getRank()][feSpaces[i]]);
for (InteriorBoundary::iterator it(boundary, level); !it.end(); ++it){
it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj,
data[level][it.getRank()][feSpaces[i]]);
}
// === Remove empty data containers. ===
for (unsigned int i = 0; i < data.size(); i++) {
DataIter dit = data[i].begin();
while (dit != data[i].end()) {
FeMapIter it = dit->second.begin();
while (it != dit->second.end()) {
if (it->second.size() == 0) {
const FiniteElemSpace *fe = it->first;
++it;
dit->second.erase(fe);
} else
++it;
}
DataIter dit = data[i].begin();
while (dit != data[i].end()) {
FeMapIter it = dit->second.begin();
while (it != dit->second.end()) {
if (it->second.size() == 0) {
const FiniteElemSpace *fe = it->first;
++it;
dit->second.erase(fe);
} else
++it;
}
if (dit->second.size() == 0)
data[i].erase(dit++);
else
++dit;
}
if (dit->second.size() == 0)
data[i].erase(dit++);
else
++dit;
}
}
}
......
......@@ -228,7 +228,7 @@ namespace AMDiS {
Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
if (writePartMesh > 0) {
debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");
debug::writeElementIndexMesh(mesh , debugOutputDir + "elementIndex");
ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
}
}
......@@ -1803,13 +1803,23 @@ namespace AMDiS {
}
// === Update DOF admins due to new number of DOFs. ===
lastMeshChangeIndex = mesh->getChangeIndex();
#if (DEBUG != 0)
static int fileNumber(0); //improvised counter for adapt Iteration
stringstream ss;
ss << debugOutputDir << "elementMaps." << fileNumber ;
//write local Element Maps to csv file
ParallelDebug::writeCsvElementMap(feSpaces[0],
*(dofMaps[0]),
ss.str(), "csv");
fileNumber++;
ParallelDebug::testDofContainerCommunication(*this);
MSG("------------- Debug information -------------\n");
......@@ -1838,11 +1848,20 @@ namespace AMDiS {
int test = 0;
Parameters::get("parallel->remove periodic boundary", test);
if (test == 0) {
ParallelDebug::testCommonDofs(*this, true);
ParallelDebug::testGlobalIndexByCoords(*this);
// if (test == 0) {
//121106 replaced if statement
if(true){
//write latest mesh properties to vtu file
debug::writeElementIndexMesh(mesh , debugOutputDir + "elementIndex");
ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
ParallelDebug::testCommonDofs(*this, true);
ParallelDebug::testGlobalIndexByCoords(*this);
}
#else
for (int i = 0; i < static_cast<int>(dofMaps.size()); i++) {
vector<const FiniteElemSpace*>& dofMapSpaces = dofMaps[i]->getFeSpaces();
......
......@@ -19,6 +19,7 @@
#include "StdMpi.h"
#include "Debug.h"
#include "io/VtkWriter.h"
#include "ElementDofIterator.h"
namespace AMDiS {
......@@ -914,4 +915,87 @@ namespace AMDiS {
bound.neighObj.elIndex, mesh->getDim(),
code.toStr().c_str());
}
void ParallelDebug::writeCsvElementMap(const FiniteElemSpace *feSpace,
ParallelDofMapping &dofMap,
string prefix,
string postfix)
{
FUNCNAME("ParallelDebug::writeCsvElementMap()");
MSG("writing local Element map to CSV File \n");
Mesh *mesh = feSpace->getMesh();
stringstream filename;
filename << prefix << "-" << MPI::COMM_WORLD.Get_rank() << "." << postfix;
ofstream file;
file.open(filename.str().c_str());
DOFVector<WorldVector<double> > dofCoords(feSpace, "DOF coords");
mesh->getDofIndexCoords(dofCoords);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh,
-1,
Mesh::CALL_EVERY_EL_POSTORDER);
while (elInfo) {
/*
MSG("Start traverse element %d in level %d\n",
elInfo->getElement()->getIndex(),
elInfo->getLevel());
*/
// check if Element is Leaflevel
// because higherOrderDofs(=NON VERTICES DOFS) are not saved in nonleaf Elements
if (elInfo->getElement()->isLeaf()) {
//traverse ALL DOFS
ElementDofIterator elDofIter(feSpace, true);
elDofIter.reset(elInfo->getElement());
int locDofNumber = 0;
do {
WorldVector<double> c = dofCoords[elDofIter.getDof()];
file << elInfo->getElement()->getIndex() << "\t"; //element number
file << elInfo->getLevel() << "\t"; //element Level
file << elDofIter.getDof() << "\t"; //dof number
file << elDofIter.getCurrentPos()+1 << "\t"; //dof type (+1 to pseudo convert to geoIndex, does not work for CENTER DOFS( geoIndex=0))
file << c[0] << "\t"; //dof coords
file << c[1] << "\t";
file << c[2] << "\t";
file << locDofNumber << "\t"; //local Dof number
file << elInfo->getType() << "\n"; //element Type
locDofNumber++;
} while (elDofIter.next());
} else {
// traverse only VERTEX DOFS
for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
DegreeOfFreedom dof = elInfo->getElement()->getDof(i, 0);
WorldVector<double> c = dofCoords[dof];
file << elInfo->getElement()->getIndex() << "\t"; //element number
file << elInfo->getLevel() << "\t"; //element Level
file << dof << "\t"; //dof number
file << VERTEX << "\t"; //dof type
file << c[0] << "\t"; //dof coords
file << c[1] << "\t";
file << c[2] << "\t";
file << i << "\t"; //local Dof number
file << elInfo->getType() << "\n"; //element Type
}
}
elInfo = stack.traverseNext(elInfo);
}
}
}
......@@ -189,6 +189,16 @@ namespace AMDiS {
static void followBoundary(Mesh *mesh,
AtomicBoundary &bound,
MeshStructure &code);
/** \brief
* Writes Element Map of local Rank
* Map containes for each DOF in each Element (resulting in massive doubling of DOFs):
* localElementNumber elementLevel localDOFNumber dofType dofCoords(0-2) elementLocalDOFNumber typeOfElement (0-2)
*/
static void writeCsvElementMap(const FiniteElemSpace *feSpace,
ParallelDofMapping &dofMap,
string prefix,
string postfix);
};
} // namespace AMDiS
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment