Commit e10fa30c authored by Thomas Witkowski's avatar Thomas Witkowski

And some more bug fixes for parallel code.

parent 3f1c87ca
......@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
# Libtool was configured on host deimos102:
# Libtool was configured on host p1q024:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -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 deimos102:
# Libtool was configured on host p1q024:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
# Libtool was configured on host deimos102:
# Libtool was configured on host p1q024:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......
......@@ -190,6 +190,36 @@ namespace AMDiS {
}
int MeshStructure::lookAhead(unsigned int n)
{
FUNCNAME("MeshStructure::lookAhead()");
int returnValue = 0;
int tmp_pos = pos;
int tmp_currentElement = currentElement;
int tmp_currentIndex = currentIndex;
uint64_t tmp_currentCode = currentCode;
for (unsigned int i = 0; i < n; i++) {
if (nextElement() == false) {
returnValue = -1;
break;
}
}
if (returnValue != -1)
returnValue = static_cast<int>(!isLeafElement());
pos = tmp_pos;
currentElement = tmp_currentElement;
currentIndex = tmp_currentIndex;
currentCode = tmp_currentCode;
return returnValue;
}
bool MeshStructure::skipBranch(MeshStructure *insert)
{
FUNCNAME("MeshStructure::skipBranch()");
......
......@@ -91,6 +91,8 @@ namespace AMDiS {
bool nextElement(MeshStructure *insert = NULL);
int lookAhead(unsigned int n = 1);
inline bool isLeafElement()
{
return (currentCode & 1) == 0;
......
......@@ -675,6 +675,8 @@ namespace AMDiS {
bool foundEdge = false;
while (elInfo2) {
MSG("TRY TO FIND ON EL %d\n", elInfo2->getElement()->getIndex());
for (int i = 0; i < 6; i++) {
DofEdge edge2 = elInfo2->getElement()->getEdge(i);
if (edge2.first == *(edge[0]) &&
......@@ -688,13 +690,14 @@ namespace AMDiS {
// must be refined at least once to get a refinement edge.
if (i == 0) {
MSG("AND FOUND!\n");
// Edge is refinement edge, so add it to refine list.
refineList.setElType(n_neigh, elInfo2->getType());
refineList.setElement(n_neigh, elInfo2->getElement());
n_neigh++;
foundEdge = true;
TraverseStack *tmpStack = stack;
stack = &stack2;
......@@ -704,11 +707,12 @@ namespace AMDiS {
}
stack = tmpStack;
foundEdge = true;
break;
} else {
// Edge i not refinement edge, so refine the element further.
Element *el2 = elInfo->getElement();
Element *el2 = elInfo2->getElement();
el2->setMark(std::max(el2->getMark(), 1));
TraverseStack *tmpStack = stack;
......
......@@ -10,6 +10,11 @@
// See also license.opensource.txt in the distribution.
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/convenience.hpp>
#include "ElementFileWriter.h"
#include "BasisFunction.h"
#include "Parameters.h"
......@@ -18,9 +23,11 @@
namespace AMDiS {
ElementFileWriter::ElementFileWriter(std::string name_,
using namespace std;
ElementFileWriter::ElementFileWriter(string name_,
Mesh *mesh_,
std::map<int, double> &mapvec)
map<int, double> &mapvec)
: name(name_),
amdisMeshDatExt(".elem.mesh"),
vtkExt(".vtu"),
......@@ -59,7 +66,7 @@ namespace AMDiS {
if (timestepNumber != 0 && !force)
return;
std::string fn = filename;
string fn = filename;
if (appendIndex) {
TEST_EXIT(indexLength <= 99)("index lenght > 99\n");
......@@ -92,9 +99,9 @@ namespace AMDiS {
}
void ElementFileWriter::writeFile(std::map<int, double> &vec,
void ElementFileWriter::writeFile(map<int, double> &vec,
Mesh *mesh,
std::string filename,
string filename,
int level)
{
ElementFileWriter efw("", mesh, vec);
......@@ -102,27 +109,34 @@ namespace AMDiS {
}
void ElementFileWriter::writeMeshDatValues(std::string filename, double time)
void ElementFileWriter::writeMeshDatValues(string filename, double time)
{
FUNCNAME("ElementFileWriter::writeMeshDatValues()");
std::ofstream fout(filename.c_str());
TEST_EXIT(fout)("Could not open file %s !\n", filename.c_str());
boost::iostreams::filtering_ostream file;
{
//boost::iostreams seems not to truncate the file
ofstream swapfile(filename.c_str(), ios::out | ios::trunc);
TEST_EXIT(swapfile.is_open())
("Cannot open file %s for writing!\n", name.c_str());
swapfile.close();
}
file.push(boost::iostreams::file_descriptor_sink(filename, ios::trunc));
int dim = mesh->getDim();
double val;
// === Write header. ===
fout << "mesh name: " << mesh->getName().c_str() << "\n\n";
fout << "time: " << time << "\n\n";
fout << "DIM: " << dim << "\n";
fout << "DIM_OF_WORLD: " << Global::getGeo(WORLD) << "\n\n";
fout << "number of vertices: " << (dim+1)*mesh->getNumberOfLeaves() << "\n";
fout << "number of elements: " << mesh->getNumberOfLeaves() << "\n\n";
file << "mesh name: " << mesh->getName().c_str() << "\n\n";
file << "time: " << time << "\n\n";
file << "DIM: " << dim << "\n";
file << "DIM_OF_WORLD: " << Global::getGeo(WORLD) << "\n\n";
file << "number of vertices: " << (dim+1)*mesh->getNumberOfLeaves() << "\n";
file << "number of elements: " << mesh->getNumberOfLeaves() << "\n\n";
// === Write vertex coordinates (every vertex for every element). ===
fout << "vertex coordinates:\n";
file << "vertex coordinates:\n";
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh,
......@@ -135,9 +149,9 @@ namespace AMDiS {
// Write coordinates of all element vertices.
for (int i = 0; i <= dim; ++i) {
for (int j = 0; j < dim; ++j) {
fout << elInfo->getCoord(i)[j] << " ";
file << elInfo->getCoord(i)[j] << " ";
}
fout << "\n";
file << "\n";
}
elInfo = stack.traverseNext(elInfo);
......@@ -147,33 +161,33 @@ namespace AMDiS {
// === Write elements. ===
int numLeaves = mesh->getNumberOfLeaves();
int vertCntr = 0;
fout << "\n";
fout << "element vertices:\n";
file << "\n";
file << "element vertices:\n";
for (int i = 0; i < numLeaves; ++i) {
for (int j = 0; j <= dim; ++j) {
fout << vertCntr << " ";
file << vertCntr << " ";
++vertCntr;
}
fout << "\n";
file << "\n";
}
// === Write values. ===
// Write values header.
fout << "\n";
fout << "number of values: 1\n\n";
fout << "value description: " << name.c_str() << "\n";
fout << "number of interpolation points: 0" << "\n";
fout << "type: scalar" << "\n";
fout << "interpolation type: lagrange" << "\n";
fout << "interpolation degree: 1" << "\n";
fout << "end of description: " << name.c_str() << "\n\n";
file << "\n";
file << "number of values: 1\n\n";
file << "value description: " << name.c_str() << "\n";
file << "number of interpolation points: 0" << "\n";
file << "type: scalar" << "\n";
file << "interpolation type: lagrange" << "\n";
file << "interpolation degree: 1" << "\n";
file << "end of description: " << name.c_str() << "\n\n";
// Write values.
fout << "vertex values: " << name.c_str() << "\n";
file << "vertex values: " << name.c_str() << "\n";
fout.setf(std::ios::scientific,std::ios::floatfield);
file.setf(ios::scientific,ios::floatfield);
elInfo = stack.traverseFirst(mesh,
-1,
......@@ -185,7 +199,7 @@ namespace AMDiS {
// Write value for each vertex of each element.
for (int i = 0; i <= dim; ++i) {
fout << val << "\n";
file << val << "\n";
}
elInfo = stack.traverseNext(elInfo);
......@@ -193,21 +207,25 @@ namespace AMDiS {
// Write values trailor.
fout << "\n";
fout << "interpolation values: " << name.c_str() << "\n\n\n";
fout << "element interpolation points: " << name.c_str() << "\n";
fout.close();
file << "\n";
file << "interpolation values: " << name.c_str() << "\n\n\n";
file << "element interpolation points: " << name.c_str() << "\n";
}
void ElementFileWriter::writeVtkValues(std::string filename, int level)
void ElementFileWriter::writeVtkValues(string filename, int level)
{
FUNCNAME("ElementFileWriter::writeVtkValues()");
std::ofstream fout(filename.c_str());
TEST_EXIT(fout)("Could not open file %s !\n", filename.c_str());
boost::iostreams::filtering_ostream file;
{
//boost::iostreams seems not to truncate the file
ofstream swapfile(filename.c_str(), ios::out | ios::trunc);
TEST_EXIT(swapfile.is_open())
("Cannot open file %s for writing!\n", name.c_str());
swapfile.close();
}
file.push(boost::iostreams::file_descriptor_sink(filename, ios::trunc));
int dim = mesh->getDim();
int dimOfWorld = Global::getGeo(WORLD);
......@@ -226,13 +244,13 @@ namespace AMDiS {
}
}
fout << "<?xml version=\"1.0\"?>" << std::endl;
fout << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
fout << " <UnstructuredGrid>" << std::endl;
fout << " <Piece NumberOfPoints=\"" << (dim + 1) * nElements
<< "\" NumberOfCells=\"" << nElements << "\">" << std::endl;
fout << " <Points>" << std::endl;
fout << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
file << "<?xml version=\"1.0\"?>\n";
file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
file << " <UnstructuredGrid>\n";
file << " <Piece NumberOfPoints=\"" << (dim + 1) * nElements
<< "\" NumberOfCells=\"" << nElements << "\">\n";
file << " <Points>\n";
file << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
// === Write vertex coordinates (every vertex for every element). ===
......@@ -244,61 +262,61 @@ namespace AMDiS {
// Write coordinates of all element vertices.
for (int i = 0; i <= dim; i++) {
for (int j = 0; j < dimOfWorld; j++)
fout << elInfo->getCoord(i)[j] << " ";
file << elInfo->getCoord(i)[j] << " ";
for (int j = dimOfWorld; j < 3; j++)
fout << "0.0";
file << "0.0";
fout << "\n";
file << "\n";
}
elInfo = stack.traverseNext(elInfo);
}
fout << " </DataArray>" << std::endl;
fout << " </Points>" << std::endl;
fout << " <Cells>" << std::endl;
file << " </DataArray>\n";
file << " </Points>\n";
file << " <Cells>\n";
fout << " <DataArray type=\"Int32\" Name=\"offsets\">" << std::endl;
file << " <DataArray type=\"Int32\" Name=\"offsets\">\n";
for (int i = 0; i < nElements; i++)
fout << " " << (i + 1) * vertices << std::endl;
fout << " </DataArray>" << std::endl;
file << " " << (i + 1) * vertices << "\n";
file << " </DataArray>\n";
fout << " <DataArray type=\"UInt8\" Name=\"types\">" << std::endl;
file << " <DataArray type=\"UInt8\" Name=\"types\">\n";
for (int i = 0; i < nElements; i++) {
switch (vertices) {
case 2:
fout << " 3" << std::endl;
file << " 3\n";
break;
case 3:
fout << " 5" << std::endl;
file << " 5\n";
break;
case 4:
fout << " 10" << std::endl;
file << " 10\n";
break;
default:
break;
}
}
fout << " </DataArray>" << std::endl;
file << " </DataArray>\n";
fout << " <DataArray type=\"Int32\" Name=\"connectivity\">" << std::endl;
file << " <DataArray type=\"Int32\" Name=\"connectivity\">\n";
int vertCntr = 0;
for (int i = 0; i < nElements; ++i) {
for (int j = 0; j <= dim; ++j) {
fout << vertCntr << " ";
file << vertCntr << " ";
++vertCntr;
}
fout << std::endl;
file << "\n";
}
fout << " </DataArray>" << std::endl;
file << " </DataArray>\n";
fout << " </Cells>" << std::endl;
fout << " <PointData>" << std::endl;
fout << " <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">" << std::endl;
file << " </Cells>\n";
file << " <PointData>\n";
file << " <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">\n";
fout.setf(std::ios::scientific,std::ios::floatfield);
file.setf(ios::scientific,ios::floatfield);
elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
int vc = 0;
......@@ -308,20 +326,18 @@ namespace AMDiS {
// Write value for each vertex of each element.
for (int i = 0; i <= dim; i++)
fout << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
file << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
vc++;
elInfo = stack.traverseNext(elInfo);
}
fout << " </DataArray>" << std::endl;
file << " </DataArray>\n";
fout << " </PointData>" << std::endl;
fout << " </Piece>" << std::endl;
fout << " </UnstructuredGrid>" << std::endl;
fout << "</VTKFile>" << std::endl;
fout.close();
file << " </PointData>\n";
file << " </Piece>\n";
file << " </UnstructuredGrid>\n";
file << "</VTKFile>\n";
}
}
......@@ -53,19 +53,9 @@ namespace AMDiS {
std::ofstream swapfile(name.c_str(), std::ios::out | std::ios::trunc);
swapfile.close();
}
file.push(boost::iostreams::file_descriptor_sink(name, std::ios::trunc));
file.push(boost::iostreams::file_descriptor_sink(name, std::ios::trunc));
writeFileToStream(file);
#if 0
std::ofstream file;
file.open(name.c_str());
TEST_EXIT(file.is_open())("Cannot open file %s for writing\n", name.c_str());
writeFileToStream(file);
file.close();
#endif
return 0;
}
......
......@@ -818,6 +818,10 @@ namespace AMDiS {
nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
repartitionMesh();
nMeshChangesAfterLastRepartitioning = 0;
} else {
MSG_DBG("Repartitioning not tried because tryRepartitioning = %d reparttitioningAllowed = %d nMeshChange =%d repartitionIthChange = %d\n",
tryRepartition, repartitioningAllowed,
nMeshChangesAfterLastRepartitioning, repartitionIthChange);
}
// === Print imbalance factor. ===
......@@ -881,6 +885,10 @@ namespace AMDiS {
elCode.init(boundIt->rankObj);
if (elCode.getCode() != recvCodes[i].getCode()) {
// MSG("MACRO EL %d %d %d DOES NOT FIT WITH MACRO EL %d %d %d ON RANK %d\n",
// boundIt->rankObj.elIndex, boundIt->rankObj.subObj, boundIt->rankObj.ithObj,
// boundIt->neighObj.elIndex, boundIt->neighObj.subObj, boundIt->neighObj.ithObj,
// it->first);
TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
MeshManipulation mm(feSpace);
......@@ -1051,7 +1059,6 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
// === Run mesh partitioner to calculate a new mesh partitioning. ===
partitioner->setLocalGlobalDofMap(&mapLocalGlobalDofs);
......
......@@ -374,82 +374,111 @@ namespace AMDiS {
bool meshChanged = false;
Element *el = elInfo->getElement();
int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
// === If element is leaf (and the code is not), refine the element. ===
bool elementRefined = true;
if (el->isLeaf()) {
TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");
el->setMark(1);
refineManager->setMesh(el->getMesh());
refineManager->setStack(&stack);
refineManager->refineFunction(elInfo);
meshChanged = true;
}
// === Continue fitting the mesh structure code to the children of the ===
// === current element. ===
// In some situations refinement of an element can be omitted, altough
// it is defined in the mesh structure code. One case is when the
// following holds:
// - the code is used to adapt a face of a tetrahedron
// - only one of the children of the current element contains this face
// - the face to be adapted of the current element is either the
// local face 0 or 1
// - the next structure code is 0, i.e., we would be finished after
// the refinement of this element
// In this scenario, the element must be refined due to the structure
// code, but the refinement does not introduce new DOFs on the face,
// that should be adapted. Thus, we can ommit the refinement.
if (subObj == FACE) {
if (s0 != -1 && s1 == -1 || s0 == -1 && s1 != -1) {
if (ithObj <= 1 && code.lookAhead() == 0) {
elementRefined = false;
code.nextElement();
stack.traverseNext(elInfo);
}
}
}
int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
Element *child0 = el->getFirstChild();
Element *child1 = el->getSecondChild();
if (reverseMode) {
swap(s0, s1);
swap(child0, child1);
if (elementRefined) {
MeshStructure t = code;
el->setMark(1);
refineManager->setMesh(el->getMesh());
refineManager->setStack(&stack);
refineManager->refineFunction(elInfo);
meshChanged = true;
}
}
// === Traverse left child. ===
if (s0 != -1) {
// The edge/face is contained in the left child, therefore fit this
// child to the mesh structure code.
stack.traverseNext(elInfo);
code.nextElement();
meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
elInfo = stack.getElInfo();
} else {
// The edge/face is not contained in the left child. Hence we need
// to traverse through all subelements of the left child until we get
// the second child of the current element.
do {
elInfo = stack.traverseNext(elInfo);
} while (elInfo && elInfo->getElement() != child1);
TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
}
TEST_EXIT_DBG(elInfo->getElement() == child1)
("Got wrong child with idx = %d! Searched for child idx = %d\n",
elInfo->getElement()->getIndex(), child1->getIndex());
// === Traverse right child. ===
if (s1 != -1) {
// The edge/face is contained in the right child, therefore fit this
// child to the mesh structure code.
code.nextElement();
meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
} else {
// The edge/face is not contained in the right child. Hence we need
// to traverse through all subelements of the right child until we have
// finished traversing the current element with all its subelements.
int level = elInfo->getLevel();
if (elementRefined) {
// === Continue fitting the mesh structure code to the children of the ===
// === current element. ===
Element *child0 = el->getFirstChild();
Element *child1 = el->getSecondChild();
if (reverseMode) {
swap(s0, s1);
swap(child0, child1);
}
do {
elInfo = stack.traverseNext(elInfo);
} while (elInfo && elInfo->getLevel() > level);
// === Traverse left child. ===
if (s0 != -1) {
// The edge/face is contained in the left child, therefore fit this
// child to the mesh structure code.
stack.traverseNext(elInfo);
code.nextElement();
meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);