Commit 1f3eed71 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Bugfixing around SolutionDataStorage

parent 535ca95d
......@@ -53,13 +53,14 @@ namespace AMDiS {
usedCount = src.usedCount;
holeCount = src.holeCount;
sizeUsed = src.sizeUsed;
for (int i = 0; i < mesh->getDim(); i++) {
for (int i = 0; i <= mesh->getDim(); i++) {
nrDOF[i] = src.nrDOF[i];
nr0DOF[i] = src.nr0DOF[i];
}
dofIndexedList = src.dofIndexedList;
dofContainerList = src.dofContainerList;
}
return *this;
}
/****************************************************************************/
......
......@@ -76,32 +76,33 @@ namespace AMDiS {
elementVector(NULL),
boundaryManager(NULL),
nBasFcts(0)
{};
{}
DOFVectorBase(const FiniteElemSpace *f, std::string n);
virtual ~DOFVectorBase();
virtual const T *getLocalVector(const Element *el, T* localVec) const;
virtual const T *getLocalVector(const Element *el,
T* localVec) const;
const T *getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const T *getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
T *vecAtQPs) const;
T *vecAtQPs) const;
const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
WorldVector<T> *grdAtQPs) const;
WorldVector<T> *grdAtQPs) const;
const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
WorldMatrix<T> *d2AtQPs) const;
WorldMatrix<T> *d2AtQPs) const;
virtual const FiniteElemSpace* getFESpace() const {
inline const FiniteElemSpace* getFESpace() const {
return feSpace;
};
}
ElementVector *assemble(T factor, ElInfo *elInfo,
const BoundaryType *bound,
......@@ -119,32 +120,32 @@ namespace AMDiS {
operators.push_back(op);
operatorFactor.push_back(factor);
operatorEstFactor.push_back(estFactor);
};
}
inline std::vector<double*>::iterator getOperatorFactorBegin() {
return operatorFactor.begin();
};
}
inline std::vector<double*>::iterator getOperatorFactorEnd() {
return operatorFactor.end();
};
}
inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
return operatorEstFactor.begin();
};
}
inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
return operatorEstFactor.end();
};
}
inline std::vector<Operator*>::iterator getOperatorsBegin() {
return operators.begin();
};
}
inline std::vector<Operator*>::iterator getOperatorsEnd() {
return operators.end();
};
}
Flag getAssembleFlag();
......@@ -158,30 +159,30 @@ namespace AMDiS {
inline std::vector<Operator*>& getOperators() {
return operators;
};
}
inline std::vector<double*>& getOperatorFactor() {
return operatorFactor;
};
}
inline std::vector<double*>& getOperatorEstFactor() {
return operatorEstFactor;
};
}
/** \brief
* Returns \ref name
*/
inline const std::string& getName() const {
return name;
};
}
inline BoundaryManager* getBoundaryManager() const {
return boundaryManager;
};
}
inline void setBoundaryManager(BoundaryManager *bm) {
boundaryManager = bm;
};
}
protected:
/** \brief
......@@ -285,26 +286,28 @@ namespace AMDiS {
public:
Iterator(DOFIndexed<T> *c, DOFIteratorType type)
: DOFIterator<T>(c, type)
{};
{}
Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
: DOFIterator<T>(admin, c, type)
{};
{}
};
class Creator : public CreatorInterface<DOFVector<T> > {
public:
MEMORY_MANAGED(Creator);
Creator(FiniteElemSpace *feSpace_) : feSpace(feSpace_) {};
Creator(FiniteElemSpace *feSpace_)
: feSpace(feSpace_)
{}
DOFVector<T> *create() {
return NEW DOFVector<T>(feSpace, "");
};
}
void free(DOFVector<T> *vec) {
DELETE vec;
};
}
private:
FiniteElemSpace *feSpace;
......@@ -351,14 +354,14 @@ namespace AMDiS {
*/
typename std::vector<T>::iterator begin() {
return vec.begin();
};
}
/** \brief
* Returns iterator to the end of \ref vec
*/
typename std::vector<T>::iterator end() {
return vec.end();
};
}
/** \brief
* Used by DOFAdmin to compress this DOFVector. Implementation of
......@@ -469,7 +472,7 @@ namespace AMDiS {
*/
inline double L2Norm(Quadrature* q = NULL) const {
return sqrt(L2NormSquare());
};
}
/** \brief
* Calculates square of L2 norm of this DOFVector
......@@ -481,7 +484,7 @@ namespace AMDiS {
*/
inline double H1Norm(Quadrature* q = NULL) const {
return sqrt(H1NormSquare());
};
};
/** \brief
* Calculates square of H1 norm of this DOFVector
......@@ -575,6 +578,9 @@ namespace AMDiS {
*/
void print() const;
/** \brief
*
*/
int calcMemoryUsage() const;
/** \brief
......@@ -692,11 +698,11 @@ namespace AMDiS {
inline double min(const DOFVector<double>& v) {
return v.min();
};
}
inline double max(const DOFVector<double>& v) {
return v.max();
};
}
// ===========================================================================
// ===== class DOFVectorDOF ==================================================
......@@ -720,14 +726,14 @@ namespace AMDiS {
: DOFVector<DegreeOfFreedom>(feSpace_, name_)
{
feSpace->getAdmin()->addDOFContainer(this);
};
}
/** \brief
* Deregisters itself at DOFAdmin.
*/
~DOFVectorDOF() {
feSpace->getAdmin()->removeDOFContainer(this);
};
}
/** \brief
* Implements DOFContainer::operator[]() by calling
......@@ -735,14 +741,14 @@ namespace AMDiS {
*/
DegreeOfFreedom& operator[](DegreeOfFreedom i) {
return DOFVector<DegreeOfFreedom>::operator[](i);
};
}
/** \brief
* Implements DOFIndexedBase::getSize()
*/
int getSize() const {
return DOFVector<DegreeOfFreedom>::getSize();
};
}
/** \brief
* Implements DOFIndexedBase::resize()
......@@ -765,22 +771,22 @@ namespace AMDiS {
template<typename T>
double norm(DOFVector<T> *vec) {
return vec->nrm2();
};
}
template<typename T>
double L2Norm(DOFVector<T> *vec) {
return vec->L2Norm();
};
}
template<typename T>
double H1Norm(DOFVector<T> *vec) {
return vec->H1Norm();
};
}
template<typename T>
void print(DOFVector<T> *vec) {
vec->print();
};
}
// point wise multiplication
template<typename T>
......@@ -834,7 +840,9 @@ namespace AMDiS {
void xpay(double a,const DOFVector<T>& x,DOFVector<T>& y);
template<typename T>
inline void scal(T a, DOFVector<T>& y) {y*=a;};
inline void scal(T a, DOFVector<T>& y) {
y *= a;
}
template<typename T>
inline const DOFVector<T>& mult(double scal,
......@@ -861,21 +869,19 @@ namespace AMDiS {
template<typename T>
inline void set(DOFVector<T>& vec, T d)
{
inline void set(DOFVector<T>& vec, T d) {
vec.set(d);
};
}
template<typename T>
inline void setValue(DOFVector<T>& vec, T d)
{
inline void setValue(DOFVector<T>& vec, T d) {
vec.set(d);
};
}
template<typename T>
inline int size(DOFVector<T> *vec) {
return vec->getUsedSize();
};
}
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *result);
......
......@@ -23,7 +23,6 @@ namespace AMDiS {
// === get admin ===
localAdmin_ = const_cast<DOFAdmin*>(feSpace->getAdmin());
// === create vertex info vector ===
vertexInfos_ = NEW DOFVector< std::list<VertexInfo> >(feSpace, "vertex infos");
interpPointInd_ = NEW DOFVector<int>(feSpace, "interpolation point indices");
......@@ -100,8 +99,9 @@ namespace AMDiS {
// Traverse elements to create element information
elInfo = stack.traverseFirst(mesh_, level_, flag);
while (elInfo) {
if (!writeElem_ || writeElem_(elInfo))
if (!writeElem_ || writeElem_(elInfo)) {
addElementData(elInfo);
}
elInfo = stack.traverseNext(elInfo);
}
......@@ -124,7 +124,7 @@ namespace AMDiS {
basisFcts_ = const_cast<BasisFunction*>(feSpace_->getBasisFcts());
nBasisFcts_ = basisFcts_->getNumber();
localDOFs_ = GET_MEMORY(DegreeOfFreedom, nBasisFcts_);
TraverseStack stack;
// Traverse elements to add value information and to mark
......@@ -133,8 +133,9 @@ namespace AMDiS {
level_,
traverseFlag_ | Mesh::FILL_COORDS);
while (elInfo) {
if (!writeElem_ || writeElem_(elInfo))
if (!writeElem_ || writeElem_(elInfo)) {
addValueData(elInfo);
}
elInfo = stack.traverseNext(elInfo);
}
......@@ -158,7 +159,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
}
FREE_MEMORY(localDOFs_, DegreeOfFreedom, nBasisFcts_);
valueDataCollected_ = true;
......@@ -231,12 +232,10 @@ namespace AMDiS {
} else {
elementInfo.elementRegion = -1;
}
// read surface regions to element info
ed = elInfo->getElement()->getElementData(SURFACE_REGION);
elementInfo.surfaceRegions.set(-1);
while (ed) {
SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed);
elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion();
......@@ -250,7 +249,7 @@ namespace AMDiS {
// get dof index of this vertex
vertexDOF = dof[i][nPreDofs_];
// search for coords at this dof
std::list<VertexInfo>::iterator it =
find((*vertexInfos_)[vertexDOF].begin(),
......@@ -280,6 +279,7 @@ namespace AMDiS {
outputIndices_[elInfo->getNeighbour(i)->getIndex()] :
-1;
}
if (dim_ == 3) {
elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());
......@@ -294,12 +294,11 @@ namespace AMDiS {
int DataCollector::addValueData(ElInfo *elInfo)
{
FUNCNAME("DataCollector::addValueData()");
basisFcts_->getLocalIndices(elInfo->getElement(), localAdmin_, localDOFs_);
// WorldVector<double> vertexCoords;
// First, traverse all DOFs at the vertices of the element, determine
// their coordinates and add them to the corresponding entry in dofCoords_.
// their coordinates and add them to the corresponding entry in dofCoords_.
for (int i = 0; i < mesh_->getGeo(VERTEX); i++) {
DegreeOfFreedom dofi = localDOFs_[i];
......@@ -307,7 +306,7 @@ namespace AMDiS {
// get coords of this vertex
*vertexCoords = elInfo->getCoord(i);
// search for coords at this dof
std::list<WorldVector<double> >::iterator it =
find((*dofCoords_)[dofi].begin(),
......@@ -321,8 +320,6 @@ namespace AMDiS {
}
}
// Then, traverse all interpolation DOFs of the element, determine
// their coordinates and add them to the corresponding entry in
// interpPointCoords_.
......@@ -348,7 +345,7 @@ namespace AMDiS {
}
}
}
return(0);
}
......
......@@ -121,10 +121,10 @@ namespace AMDiS {
/* =========== And clone the children ============= */
if (child[0]) {
el->child[0] = child[0]->clone();
el->child[0] = child[0]->cloneWithDOFs();
}
if (child[1]) {
el->child[1] = child[1]->clone();
el->child[1] = child[1]->cloneWithDOFs();
}
return el;
......
......@@ -58,6 +58,7 @@ namespace AMDiS {
mesh = NEW Mesh(feSpace.mesh->getName(), feSpace.mesh->getDim());
*mesh = *(feSpace.mesh);
admin = &(const_cast<DOFAdmin&>(mesh->getDOFAdmin(0)));
basFcts = feSpace.basFcts;
TEST_EXIT(feSpace.admin == &(feSpace.mesh->getDOFAdmin(0)))
("Gut, dass muss ich mir noch mal ueberlegen!\n");
......
......@@ -708,7 +708,7 @@ namespace AMDiS {
}
int vertex[3];
int** dof = const_cast<int**>( el->getDOF());
int** dof = const_cast<int**>(el->getDOF());
int verticesOfPosition = dimOfPosition + 1;
for (int i = 0; i < verticesOfPosition; i++) {
......@@ -931,12 +931,13 @@ namespace AMDiS {
for (int i = 0; i < num; node0++, i++) {
indi = orderOfPositionIndices(el, posIndex, i);
for (int k = 0; k < nrDOFs; k++)
for (int k = 0; k < nrDOFs; k++) {
result[j++] = dof[node0][n0 + indi[k]];
}
}
}
}
return result;
}
......
......@@ -170,6 +170,8 @@ namespace AMDiS {
int insertCounter = 0;
macroElements.clear();
// Go through all MacroElements of mesh m, and create for every a new
// MacroElement in this mesh.
for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
......
......@@ -65,6 +65,15 @@ namespace AMDiS {
bool pop(T **solution,
typename SolutionHelper<T>::type feSpace,
double *timestep);
T* getSolution() {
return solutions[0];
}
/** \brief
* Deletes all pointers and empties all internal vectors.
*/
void clear();
/** \brief
......@@ -109,11 +118,6 @@ namespace AMDiS {
}
}
/** \brief
* Deletes all pointers and empties all internal vectors.
*/
void cleanup();
/** \brief
* Number of MBytes of memory that can be used for solution storage.
*/
......
......@@ -23,7 +23,7 @@ namespace AMDiS {
template<typename T>
SolutionDataStorage<T>::~SolutionDataStorage()
{
cleanup();
clear();
}
template<typename T>
......@@ -45,7 +45,7 @@ namespace AMDiS {
{
// If pop was the last operation, cleanup and reset the data storage.
if (poped) {
cleanup();
clear();
poped = false;
}
......@@ -110,21 +110,25 @@ namespace AMDiS {
template<typename T>
void SolutionDataStorage<T>::cleanup()
void SolutionDataStorage<T>::clear()
{
for (int i = 0; i < static_cast<int>(solutions.size()); i++) {
DELETE solutions[i];
}
std::cout << "MARK 1" << std::endl;
if (!fixedFESpace) {
for (int i = 0; i < static_cast<int>(feSpaces.size()); i++) {
deleteFeSpace(feSpaces[i]);
}
}
solutions.empty();
feSpaces.empty();
timestamps.empty();
std::cout << "MARK 2" << std::endl;
solutions.clear();
feSpaces.clear();
timestamps.clear();
lastPos = -1;
}
......
......@@ -198,12 +198,15 @@ namespace AMDiS {
}
/** \brief
* Returns \ref feSpace.
* Returns the fe space for a given component.
*/
inline FiniteElemSpace *getFESpace(int i) const {
return feSpace[i];
}
/** \brief
* Returns the fe spaces for all components.
*/
inline std::vector<FiniteElemSpace*> getFESpaces() const {
return feSpace;
}
......
......@@ -79,6 +79,7 @@ namespace AMDiS {
const char *filename)
{
DataCollector *dc = NEW DataCollector(values->getFESpace(), values);
std::vector<DataCollector*> dcList(0);
dcList.push_back(dc);
......
......@@ -21,7 +21,6 @@ namespace AMDiS {
nVertices += (*dataCollector)[0]->getNumberInterpPoints();
nElements *= 16;
}
file << "<?xml version=\"1.0\"?>\n";
file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
file << " <UnstructuredGrid>\n";
......@@ -62,7 +61,7 @@ namespace AMDiS {
file << " </DataArray>\n";
file << " <DataArray type=\"Int32\" Name=\"connectivity\">\n";
writeConnectivity(file);
file << " </DataArray>\n";
......@@ -72,7 +71,7 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(dataCollector->size()); i++) {
file << " <DataArray type=\"Float32\" Name=\"value" << i
<< "\" format=\"ascii\">\n";
writeVertexValues(file, i);
file << " </DataArray>\n";
......