From d1ff5073612710851f71fa099bee4578796fff6a Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 30 Mar 2010 08:39:30 +0000
Subject: [PATCH] Argh, some error, go back to revision 1027.
---
AMDiS/src/Assembler.cc | 5 +-
AMDiS/src/BasisFunction.h | 18 +++--
AMDiS/src/BoundaryManager.cc | 67 +++++++++---------
AMDiS/src/BoundaryManager.h | 3 +
AMDiS/src/DOFMatrix.cc | 35 +++++++---
AMDiS/src/DOFVector.cc | 23 +++----
AMDiS/src/DOFVector.h | 4 +-
AMDiS/src/DOFVector.hh | 68 +++----------------
AMDiS/src/DualTraverse.h | 13 ----
AMDiS/src/ElInfo.h | 20 +-----
AMDiS/src/Error.hh | 55 +++++++--------
AMDiS/src/Estimator.cc | 15 -----
AMDiS/src/Lagrange.cc | 41 ++---------
AMDiS/src/Lagrange.h | 5 --
AMDiS/src/Mesh.cc | 32 ++++-----
AMDiS/src/Mesh.h | 3 -
AMDiS/src/ProblemVec.cc | 32 +++------
AMDiS/src/QPInfo.h | 120 +++++++++++++++++++++++++--------
AMDiS/src/Recovery.cc | 47 ++++++-------
AMDiS/src/RecoveryEstimator.cc | 4 +-
AMDiS/src/ResidualEstimator.cc | 19 +++---
AMDiS/src/RobinBC.cc | 15 +++--
AMDiS/src/Traverse.cc | 80 +++++++++++-----------
AMDiS/src/Traverse.h | 44 ------------
AMDiS/src/TraverseParallel.h | 6 --
25 files changed, 323 insertions(+), 451 deletions(-)
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 13e5893c..611fd67a 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -252,11 +252,12 @@ namespace AMDiS {
{
FUNCNAME("Assembler::matVecAssemble()");
+ Element *el = elInfo->getElement();
double *uhOldLoc = new double[nRow];
- operat->uhOld->getLocalVector(elInfo, uhOldLoc);
+ operat->uhOld->getLocalVector(el, uhOldLoc);
- if (elInfo->getElement() != lastMatEl) {
+ if (el != lastMatEl) {
set_to_zero(elementMatrix);
calculateElementMatrix(elInfo, elementMatrix);
}
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index ce045c73..6203bb6c 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -275,7 +275,7 @@ namespace AMDiS {
virtual void coarseRestr(DOFVector<WorldVector<double> >*, RCNeighbourList*, int)
{}
- /// Returns local dof indices of the element for the given DOF admin.
+ /// Returns local dof indices of the element for the given fe space.
virtual const DegreeOfFreedom *getLocalIndices(const Element *el,
const DOFAdmin *admin,
DegreeOfFreedom *dofPtr) const
@@ -283,11 +283,17 @@ namespace AMDiS {
return NULL;
}
- /// Calculates local dof indices of the element for the given DOF admin.
- virtual void getLocalIndices(const Element *el,
- const DOFAdmin *admin,
- std::vector<DegreeOfFreedom> &indices) const
- {}
+ inline void getLocalIndices(const Element *el,
+ const DOFAdmin *admin,
+ std::vector<DegreeOfFreedom> &indices) const
+ {
+ FUNCNAME("BasisFunction::getLocalIndices()");
+
+ TEST_EXIT_DBG(static_cast<int>(indices.size()) >= nBasFcts)
+ ("Index vector is too small!\n");
+
+ getLocalIndices(el, admin, &(indices[0]));
+ }
virtual void getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
diff --git a/AMDiS/src/BoundaryManager.cc b/AMDiS/src/BoundaryManager.cc
index 57a464a8..6fd1a426 100644
--- a/AMDiS/src/BoundaryManager.cc
+++ b/AMDiS/src/BoundaryManager.cc
@@ -15,6 +15,7 @@ namespace AMDiS {
BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
{
localBounds.resize(omp_get_overall_max_threads());
+ dofIndices.resize(omp_get_overall_max_threads());
allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber();
for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
@@ -25,6 +26,7 @@ namespace AMDiS {
localBCs = bm.localBCs;
allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds;
localBounds.resize(bm.localBounds.size());
+ dofIndices.resize(bm.localBounds.size());
for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
}
@@ -47,61 +49,54 @@ namespace AMDiS {
return result;
}
-
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
DOFVectorBase<double> *vec)
{
- FUNCNAME("BoundaryManager::fillBoundaryConditions()");
-
- if (localBCs.size() <= 0)
- return;
-
- const FiniteElemSpace *feSpace = vec->getFESpace();
- const BasisFunction *basisFcts = feSpace->getBasisFcts();
- int nBasFcts = basisFcts->getNumber();
-
- // get boundaries of all DOFs
- BoundaryType *localBound = localBounds[omp_get_thread_num()];
- basisFcts->getBound(elInfo, localBound);
-
- // get dof indices
- std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace);
- TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
- ("Local index vector is too small!\n");
-
- // apply non dirichlet boundary conditions
- for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second && !(*it).second->isDirichlet())
- (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
- localBound, nBasFcts);
-
- // apply dirichlet boundary conditions
- for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second && (*it).second->isDirichlet())
- (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
- localBound, nBasFcts);
+ if (localBCs.size() > 0) {
+ const FiniteElemSpace *feSpace = vec->getFESpace();
+ std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
+ const BasisFunction *basisFcts = feSpace->getBasisFcts();
+ int nBasFcts = basisFcts->getNumber();
+ dofVec.resize(nBasFcts);
+
+ // get boundaries of all DOFs
+ BoundaryType *localBound = localBounds[omp_get_thread_num()];
+ basisFcts->getBound(elInfo, localBound);
+
+ // get dof indices
+ basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
+
+ // apply non dirichlet boundary conditions
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && !(*it).second->isDirichlet())
+ (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
+ localBound, nBasFcts);
+
+ // apply dirichlet boundary conditions
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && (*it).second->isDirichlet())
+ (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
+ localBound, nBasFcts);
+ }
}
-
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, DOFMatrix *mat)
{
- FUNCNAME("BoundaryManager::fillBoundaryConditions()");
-
if (localBCs.size() <= 0)
return;
const FiniteElemSpace *feSpace = mat->getRowFESpace();
+ std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
+ dofVec.resize(nBasFcts);
// get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound);
// get dof indices
- std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace);
- TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
- ("Local index vector is too small!\n");
+ basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
// apply non dirichlet boundary conditions
for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h
index 98615d23..c9f391a0 100644
--- a/AMDiS/src/BoundaryManager.h
+++ b/AMDiS/src/BoundaryManager.h
@@ -125,6 +125,9 @@ namespace AMDiS {
/// Temporary thread-safe variable for functions fillBoundaryconditions.
std::vector<BoundaryType*> localBounds;
+ /// Temporary thread-safe variable for functions fillBoundaryconditions.
+ std::vector<std::vector<DegreeOfFreedom> > dofIndices;
+
/** \brief
* Stores the number of byte that were allocated in the constructor for
* each localBounds value. Is used to free the memory in the destructor.
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 163f7aa6..fcc63e76 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -168,12 +168,9 @@ namespace AMDiS {
// === Get indices mapping from local to global matrix indices. ===
- std::vector<DegreeOfFreedom> &rowIndices2 = rowElInfo->getLocalIndices(rowFESpace);
- std::vector<DegreeOfFreedom> &colIndices2 = rowIndices2;
-
-// rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
-// rowFESpace->getAdmin(),
-// rowIndices);
+ rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
+ rowFESpace->getAdmin(),
+ rowIndices);
if (rowFESpace == colFESpace) {
colIndices = rowIndices;
} else {
@@ -191,8 +188,21 @@ namespace AMDiS {
using namespace mtl;
+#if 0
+ std::cout << "----- PRINT MAT --------" << std::endl;
+ std::cout << elMat << std::endl;
+ std::cout << "rows: ";
+ for (int i = 0; i < rowIndices.size(); i++)
+ std::cout << rowIndices[i] << " ";
+ std::cout << std::endl;
+ std::cout << "cols: ";
+ for (int i = 0; i < colIndices.size(); i++)
+ std::cout << colIndices[i] << " ";
+ std::cout << std::endl;
+#endif
+
for (int i = 0; i < nRow; i++) {
- DegreeOfFreedom row = rowIndices2[i];
+ DegreeOfFreedom row = rowIndices[i];
BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
@@ -200,15 +210,20 @@ namespace AMDiS {
if (condition && condition->isDirichlet()) {
if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- if ((*rankDofs)[row])
+ if ((*rankDofs)[rowIndices[i]])
applyDBCs.insert(static_cast<int>(row));
#else
applyDBCs.insert(static_cast<int>(row));
#endif
}
} else {
- for (int j = 0; j < nCol; j++)
- ins[row][colIndices2[j]] += elMat[i][j];
+ for (int j = 0; j < nCol; j++) {
+ DegreeOfFreedom col = colIndices[j];
+ double entry = elMat[i][j];
+
+ // std::cout << "ADD at " << row << " " << col << std::endl;
+ ins[row][col] += entry;
+ }
}
}
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 61389406..cb978648 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -108,16 +108,13 @@ namespace AMDiS {
// traverse mesh
std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
- stack.addFeSpace(feSpace);
- Flag fillFlag =
- Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_COORDS | Mesh::FILL_LOCAL_INDICES;
+ Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- getLocalVector(elInfo, localUh);
+ getLocalVector(elInfo->getElement(), localUh);
int localDOFNr = 0;
for (int i = 0; i < nNodes; i++) { // for all nodes
@@ -180,11 +177,9 @@ namespace AMDiS {
// traverse mesh
Mesh *mesh = feSpace->getMesh();
TraverseStack stack;
- stack.addFeSpace(feSpace);
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS |
- Mesh::FILL_LOCAL_INDICES);
+ Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
double *localUh = new double[basFcts->getNumber()];
@@ -192,7 +187,7 @@ namespace AMDiS {
double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- getLocalVector(elInfo, localUh);
+ getLocalVector(elInfo->getElement(), localUh);
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
for (int i = 0; i < dim + 1; i++) {
@@ -254,7 +249,7 @@ namespace AMDiS {
}
double *localVec = localVectors[myRank];
- getLocalVector(elInfo, localVec);
+ getLocalVector(elInfo->getElement(), localVec);
DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim);
@@ -339,7 +334,7 @@ namespace AMDiS {
}
double *localVec = localVectors[myRank];
- getLocalVector(largeElInfo, localVec);
+ getLocalVector(largeElInfo->getElement(), localVec);
const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
@@ -396,6 +391,8 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
+ Element *el = elInfo->getElement();
+
int myRank = omp_get_thread_num();
int dow = Global::getGeo(WORLD);
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
@@ -415,7 +412,7 @@ namespace AMDiS {
}
double *localVec = localVectors[myRank];
- getLocalVector(elInfo, localVec);
+ getLocalVector(el, localVec);
DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
int parts = Global::getGeo(PARTS, dim);
@@ -779,7 +776,7 @@ namespace AMDiS {
}
double *localVec = localVectors[omp_get_thread_num()];
- getLocalVector(largeElInfo, localVec);
+ getLocalVector(largeElInfo->getElement(), localVec);
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 7e88f4b9..2d42e128 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -67,9 +67,7 @@ namespace AMDiS {
* For the given element, this function returns an array of all DOFs of this
* DOFVector that are defined on this element.
*/
- const T *getLocalVector(const Element *el, T* localVec) const;
-
- const T *getLocalVector(const ElInfo *elInfo, T* localVec) const;
+ virtual const T *getLocalVector(const Element *el, T* localVec) const;
/** \brief
* Evaluates the DOF vector at a set of quadrature points defined on the
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index b8435cca..d97dede0 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -147,9 +147,9 @@ namespace AMDiS {
{
FUNCNAME("DOFVector::addElementVector()");
- std::vector<DegreeOfFreedom> &indices = elInfo->getLocalIndices(feSpace);
- TEST_EXIT_DBG(static_cast<int>(indices.size()) == nBasFcts)
- ("Local index vector is too small!\n");
+ std::vector<DegreeOfFreedom> indices(nBasFcts);
+ feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
+ indices);
for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
BoundaryCondition *condition =
@@ -454,11 +454,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
- stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
- Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
@@ -492,11 +490,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
- stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
- Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
@@ -530,11 +526,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
- stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
- Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
@@ -992,49 +986,6 @@ namespace AMDiS {
}
- template<typename T>
- const T *DOFVectorBase<T>::getLocalVector(const ElInfo *elInfo, T *d) const
- {
- FUNCNAME("DOFVectorBase<T>::getLocalVector()");
-
- TEST_EXIT_DBG(feSpace->getMesh() == elInfo->getElement()->getMesh())
- ("Element is defined on a different mesh than the DOF vector!\n");
-
- static T* localVec = NULL;
- static int localVecSize = 0;
- T *result;
-
- if (d) {
- result = d;
- } else {
-#ifdef _OPENMP
- ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
-#endif
- if (localVec && nBasFcts > localVecSize) {
- delete [] localVec;
- localVec = new T[nBasFcts];
- }
-
- if (!localVec)
- localVec = new T[nBasFcts];
-
- localVecSize = nBasFcts;
- result = localVec;
- }
-
- std::vector<DegreeOfFreedom> &localIndices =
- const_cast<ElInfo*>(elInfo)->getLocalIndices(feSpace);
-
- TEST_EXIT_DBG(static_cast<int>(localIndices.size()) == nBasFcts)
- ("Local indices vector is too small!\n");
-
- for (int i = 0; i < nBasFcts; i++)
- result[i] = (*this)[localIndices[i]];
-
- return result;
- }
-
-
template<typename T>
const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
@@ -1053,6 +1004,7 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
+ Element *el = elInfo->getElement();
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nPoints = quadrature->getNumPoints();
@@ -1076,7 +1028,7 @@ namespace AMDiS {
}
T *localVec = localVectors[omp_get_thread_num()];
- getLocalVector(elInfo, localVec);
+ getLocalVector(el, localVec);
for (int i = 0; i < nPoints; i++) {
result[i] = 0.0;
@@ -1141,11 +1093,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
- stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
- Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h
index d5edd66b..82264140 100644
--- a/AMDiS/src/DualTraverse.h
+++ b/AMDiS/src/DualTraverse.h
@@ -125,19 +125,6 @@ namespace AMDiS {
* \param[in] largeFace A local edge/face number on the large element.
*/
static int getFace(DualElInfo *dualElInfo, int largeFace);
-
- void addFeSpace(int stack, const FiniteElemSpace *feSpace)
- {
- FUNCNAME("DualTraverse::addFeSpace()");
-
- TEST_EXIT_DBG(stack == 0 || stack == 1)("Wrong stack number!\n");
-
- if (stack == 0)
- stack1.addFeSpace(feSpace);
-
- if (stack == 1)
- stack2.addFeSpace(feSpace);
- }
protected:
/** \brief
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 6e0cc6f6..bd0f3955 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -227,17 +227,6 @@ namespace AMDiS {
return parametric;
}
- /// Returns the local indices, \ref localIndices, for given FE space.
- inline std::vector<DegreeOfFreedom>& getLocalIndices(const FiniteElemSpace* feSpace)
- {
- FUNCNAME("ElInfo::getLocalIndices()");
-
- TEST_EXIT_DBG(element->isLeaf() == true)
- ("Local indices are computed only for leaf elements!\n");
-
- return localIndices[feSpace];
- }
-
/// Returns \ref refinementPath
inline unsigned long getRefinementPath() const
{
@@ -404,7 +393,7 @@ namespace AMDiS {
* parent data parentInfo.
* pure virtual => must be overriden in sub-class.
*/
- virtual void fillElInfo(int iChild, const ElInfo *parentInfo) = 0;
+ virtual void fillElInfo(int ichild, const ElInfo *parentInfo) = 0;
/** \brief
* calculates the Jacobian of the barycentric coordinates on \element and stores
@@ -536,13 +525,6 @@ namespace AMDiS {
/// Gradient of lambda.
DimVec<WorldVector<double> > grdLambda;
- /** \brief
- * If the local indices of elements should be set during mesh traverse, they are
- * stored in this variable for all differnt FE spaces that are defined on the
- * mesh. Node that local indices are computed only for leaf elements of the mesh.
- */
- std::map<const FiniteElemSpace*, std::vector<DegreeOfFreedom> > localIndices;
-
/// True, if this elInfo stores parametrized information. False, otherwise.
bool parametric;
diff --git a/AMDiS/src/Error.hh b/AMDiS/src/Error.hh
index 11df84a7..9dd255b1 100644
--- a/AMDiS/src/Error.hh
+++ b/AMDiS/src/Error.hh
@@ -28,33 +28,29 @@ namespace AMDiS {
{
FUNCNAME("Error<T>::maxErrAtQp()");
- const FiniteElemSpace *feSpace = uh->getFESpace();
-
+ const FiniteElemSpace *fe_space;
if (!(pU = &u)) {
ERROR("no function u specified; doing nothing\n");
return(-1.0);
}
- if (!(errUh = &uh) || !feSpace) {
- ERROR("no discrete function or no feSpace for it; doing nothing\n");
+ if (!(errUh = &uh) || !(fe_space = uh->getFESpace())) {
+ ERROR("no discrete function or no fe_space for it; doing nothing\n");
return(-1.0);
}
- if (!(basFct = feSpace->getBasisFcts())) {
+ if (!(basFct = fe_space->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(-1.0);
}
if (!q)
- q = Quadrature::provideQuadrature(feSpace->getMesh()->getDim(),
- 2 * feSpace->getBasisFcts()->getDegree() - 2);
+ q = Quadrature::provideQuadrature(fe_space->getMesh()->getDim(),
+ 2 * fe_space->getBasisFcts()->getDegree() - 2);
quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
double maxErr = 0.0;
TraverseStack stack;
- stack.addFeSpace(feSpace);
- ElInfo *elInfo =
- stack.traverseFirst(feSpace->getMesh(), -1,
- Mesh::FILL_COORDS | Mesh::FILL_LOCAL_INDICES |
- Mesh::CALL_LEAF_EL);
+ ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
+ Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL);
while (elInfo) {
elinfo = elInfo;
double err = 0.0;
@@ -85,25 +81,25 @@ namespace AMDiS {
{
FUNCNAME("Error<T>::H1Err()");
- const FiniteElemSpace *feSpace;
+ const FiniteElemSpace *fe_space;
writeInLeafData = writeLeafData;
component = comp;
Quadrature *q = NULL;
pGrdU = &grdU;
errUh = &uh;
- if (!(feSpace = uh.getFESpace())) {
- ERROR("no feSpace for uh; doing nothing\n");
+ if (!(fe_space = uh.getFESpace())) {
+ ERROR("no fe_space for uh; doing nothing\n");
return(0.0);
}
- if (!(basFct = feSpace->getBasisFcts())) {
+ if (!(basFct = fe_space->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(0.0);
}
- int dim = feSpace->getMesh()->getDim();
+ int dim = fe_space->getMesh()->getDim();
int deg = grdU.getDegree();
- int degree = deg ? deg : 2 * feSpace->getBasisFcts()->getDegree() - 2;
+ int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
q = Quadrature::provideQuadrature(dim, degree);
quadFast = FastQuadrature::provideFastQuadrature(basFct,
@@ -115,7 +111,7 @@ namespace AMDiS {
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
+ ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL |
Mesh::FILL_DET |
@@ -163,7 +159,7 @@ namespace AMDiS {
if (relative) {
double relNorm2 = h1Norm2 + 1.e-15;
- elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
double exact = elInfo->getElement()->getEstimation(component) / relNorm2;
if (writeInLeafData)
@@ -191,7 +187,7 @@ namespace AMDiS {
{
FUNCNAME("Error<T>::L2Err()");
- const FiniteElemSpace *feSpace;
+ const FiniteElemSpace *fe_space;
Quadrature *q = NULL;
writeInLeafData = writeLeafData;
component = comp;
@@ -201,19 +197,19 @@ namespace AMDiS {
return(0.0);
}
- if (!(errUh = &uh) || !(feSpace = uh.getFESpace())) {
- ERROR("no discrete function or no feSpace for it; doing nothing\n");
+ if (!(errUh = &uh) || !(fe_space = uh.getFESpace())) {
+ ERROR("no discrete function or no fe_space for it; doing nothing\n");
return(0.0);
}
- if (!(basFct = feSpace->getBasisFcts())) {
+ if (!(basFct = fe_space->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(0.0);
}
- int dim = feSpace->getMesh()->getDim();
+ int dim = fe_space->getMesh()->getDim();
int deg = u.getDegree();
- int degree = deg ? deg : 2 * feSpace->getBasisFcts()->getDegree() - 2;
+ int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
q = Quadrature::provideQuadrature(dim, degree);
quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
@@ -223,12 +219,11 @@ namespace AMDiS {
double *u_vec = new double[quadFast->getQuadrature()->getNumPoints()];
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
+ ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL |
Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_LOCAL_INDICES);
+ Mesh::FILL_GRD_LAMBDA);
while (elInfo) {
elinfo = elInfo;
const double *uh_vec;
@@ -266,7 +261,7 @@ namespace AMDiS {
if (relative) {
double relNorm2 = l2Norm2 + 1.e-15;
- elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
double exact = elInfo->getElement()->getEstimation(component) / relNorm2;
if (writeInLeafData)
diff --git a/AMDiS/src/Estimator.cc b/AMDiS/src/Estimator.cc
index c35f2f44..ba4ca4bc 100644
--- a/AMDiS/src/Estimator.cc
+++ b/AMDiS/src/Estimator.cc
@@ -76,11 +76,6 @@ namespace AMDiS {
FUNCNAME("Estimator::singleMeshTraverse()");
TraverseStack stack;
-
- if (traverseFlag.isSet(Mesh::FILL_LOCAL_INDICES))
- for (unsigned int i = 0; i < matrix.size(); i++)
- stack.addFeSpace(matrix[i]->getFESpace());
-
ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
while (elInfo) {
estimateElement(elInfo);
@@ -96,16 +91,6 @@ namespace AMDiS {
DualTraverse dualTraverse;
DualElInfo dualElInfo;
- if (traverseFlag.isSet(Mesh::FILL_LOCAL_INDICES)) {
- for (unsigned int i = 0; i < matrix.size(); i++) {
- if (matrix[i]->getFESpace()->getMesh() == mesh)
- dualTraverse.addFeSpace(0, matrix[i]->getFESpace());
-
- if (matrix[i]->getFESpace()->getMesh() == auxMesh)
- dualTraverse.addFeSpace(1, matrix[i]->getFESpace());
- }
- }
-
bool cont = dualTraverse.traverseFirst(mesh, auxMesh, -1, -1,
traverseFlag, traverseFlag,
dualElInfo);
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 8c50657a..9d10f8cc 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -692,7 +692,6 @@ namespace AMDiS {
return getLocalIndices(el, &admin, idof);
}
-
int* Lagrange::orderOfPositionIndices(const Element* el,
GeoIndex position,
int positionIndex) const
@@ -718,10 +717,11 @@ namespace AMDiS {
int vertex[3];
int** dof = const_cast<int**>(el->getDOF());
- const Element *refElement = Global::getReferenceElement(dim);
+ int verticesOfPosition = dimOfPosition + 1;
- for (int i = 0; i <= dimOfPosition; i++)
- vertex[i] = refElement->getVertexOfPosition(position, positionIndex, i);
+ for (int i = 0; i < verticesOfPosition; i++)
+ vertex[i] = Global::getReferenceElement(dim)->
+ getVertexOfPosition(position, positionIndex, i);
if (dimOfPosition == 1) {
if (degree == 3) {
@@ -769,7 +769,6 @@ namespace AMDiS {
return NULL;
}
-
void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
{
el_info->testFlag(Mesh::FILL_BOUND);
@@ -954,38 +953,6 @@ namespace AMDiS {
return result;
}
-
- void Lagrange::getLocalIndices(const Element* el,
- const DOFAdmin *admin,
- std::vector<DegreeOfFreedom> &indices) const
- {
- indices.resize(nBasFcts);
-
- const DegreeOfFreedom **dof = el->getDOF();
-
- int nrDOFs, n0, node0, num = 0;
- GeoIndex posIndex;
-
- for (int pos = 0, j = 0; pos <= dim; pos++) {
- posIndex = INDEX_OF_DIM(pos, dim);
- nrDOFs = admin->getNumberOfDOFs(posIndex);
-
- if (nrDOFs) {
- n0 = admin->getNumberOfPreDOFs(posIndex);
- node0 = admin->getMesh()->getNode(posIndex);
- num = Global::getGeo(posIndex, dim);
-
- for (int i = 0; i < num; node0++, i++) {
- const int *indi = orderOfPositionIndices(el, posIndex, i);
-
- for (int k = 0; k < nrDOFs; k++)
- indices[j++] = dof[node0][n0 + indi[k]];
- }
- }
- }
- }
-
-
void Lagrange::getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
std::vector<const DegreeOfFreedom*>& vec) const
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index 4cb390c4..f4a7bddb 100644
--- a/AMDiS/src/Lagrange.h
+++ b/AMDiS/src/Lagrange.h
@@ -123,11 +123,6 @@ namespace AMDiS {
const DOFAdmin *admin,
DegreeOfFreedom *dofs) const;
- /// Implements BasisFunction::getLocalIndices().
- void getLocalIndices(const Element *el,
- const DOFAdmin *admin,
- std::vector<DegreeOfFreedom> &indices) const;
-
void getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
std::vector<const DegreeOfFreedom*>& vec) const;
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 0bae77b6..5b03933e 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -33,21 +33,20 @@ namespace AMDiS {
// flags, which information should be present in the elInfo structure
//**************************************************************************
- const Flag Mesh::FILL_NOTHING = 0X000L;
- const Flag Mesh::FILL_COORDS = 0X001L;
- const Flag Mesh::FILL_BOUND = 0X002L;
- const Flag Mesh::FILL_NEIGH = 0X004L;
- const Flag Mesh::FILL_OPP_COORDS = 0X008L;
- const Flag Mesh::FILL_ORIENTATION = 0X010L;
- const Flag Mesh::FILL_DET = 0X020L;
- const Flag Mesh::FILL_GRD_LAMBDA = 0X040L;
- const Flag Mesh::FILL_LOCAL_INDICES = 0X080L;
- const Flag Mesh::FILL_ADD_ALL = 0X100L;
-
-
- const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X100L);
- const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X100L);
- const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X100L);
+ const Flag Mesh::FILL_NOTHING = 0X00L;
+ const Flag Mesh::FILL_COORDS = 0X01L;
+ const Flag Mesh::FILL_BOUND = 0X02L;
+ const Flag Mesh::FILL_NEIGH = 0X04L;
+ const Flag Mesh::FILL_OPP_COORDS = 0X08L;
+ const Flag Mesh::FILL_ORIENTATION= 0X10L;
+ const Flag Mesh::FILL_DET = 0X20L;
+ const Flag Mesh::FILL_GRD_LAMBDA = 0X40L;
+ const Flag Mesh::FILL_ADD_ALL = 0X80L;
+
+
+ const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
+ const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
+ const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L);
//**************************************************************************
// flags for Mesh traversal
@@ -61,6 +60,9 @@ namespace AMDiS {
const Flag Mesh::CALL_EL_LEVEL = 0X2000L;
const Flag Mesh::CALL_MG_LEVEL = 0X4000L ; // used in mg methods
+
+ // const Flag Mesh::USE_PARAMETRIC = 0X8000L ; // used in mg methods
+
std::vector<DegreeOfFreedom> Mesh::dof_used;
const int Mesh::MAX_DOF = 100;
std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index 87b62a7c..57369844 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -616,9 +616,6 @@ namespace AMDiS {
///
static const Flag FILL_GRD_LAMBDA;
- ///
- static const Flag FILL_LOCAL_INDICES;
-
//**************************************************************************
// flags for Mesh traversal
//**************************************************************************
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 68e7af28..827d6c5c 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -564,7 +564,7 @@ namespace AMDiS {
// here is reached already because of time adaption
allowFirstRefinement();
- TEST_EXIT_DBG(static_cast<unsigned int>(nComponents) == marker.size())
+ TEST_EXIT_DBG(static_cast<unsigned int>(nComponents == marker.size()))
("Wrong number of markers!\n");
Flag markFlag = 0;
@@ -655,11 +655,10 @@ namespace AMDiS {
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_LOCAL_INDICES |
- Mesh::FILL_GRD_LAMBDA |
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH;
if (useGetBound)
@@ -808,11 +807,10 @@ namespace AMDiS {
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_LOCAL_INDICES |
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH;
if (useGetBound)
@@ -896,8 +894,6 @@ namespace AMDiS {
useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
DualTraverse dualTraverse;
- dualTraverse.addFeSpace(0, feSpaces[0]);
- dualTraverse.addFeSpace(1, feSpaces[1]);
DualElInfo dualElInfo;
int oldElIndex0 = -1;
int oldElIndex1 = -1;
@@ -1268,12 +1264,9 @@ namespace AMDiS {
#ifdef _OPENMP
TraverseParallelStack stack(0, 1);
#else
- TraverseStack stack;
+ TraverseStack stack;
#endif
- for (unsigned int i = 0; i < feSpaces.size(); i++)
- stack.addFeSpace(feSpaces[i]);
-
// == Initialize matrix and vector. If we have to assemble in parallel, ===
// == temporary thread owned matrix and vector must be created. ===
@@ -1430,8 +1423,6 @@ namespace AMDiS {
// === Parallel boundary assemblage ===
TraverseParallelStack stack;
- for (unsigned int i = 0; i < feSpaces.size(); i++)
- stack.addFeSpace(feSpaces[i]);
#pragma omp parallel
{
@@ -1478,9 +1469,6 @@ namespace AMDiS {
// === Sequentiell boundary assemblage ===
TraverseStack stack;
- for (unsigned int i = 0; i < feSpaces.size(); i++)
- stack.addFeSpace(feSpaces[i]);
-
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) {
if (rhs->getBoundaryManager())
diff --git a/AMDiS/src/QPInfo.h b/AMDiS/src/QPInfo.h
index a7182945..dac5de36 100644
--- a/AMDiS/src/QPInfo.h
+++ b/AMDiS/src/QPInfo.h
@@ -34,10 +34,14 @@ namespace AMDiS {
class QPInfo
{
public:
- /// Sets \ref currentElInfo_ to elInfo and all valid flags to false.
+ /** \brief
+ * Sets \ref currentElInfo_ to elInfo and all valid flags to false.
+ */
void initElement(const ElInfo *elInfo);
- /// Returns coordinates at quadrature points.
+ /** \brief
+ * Returns coordinates at quadrature points.
+ */
WorldVector<double> *getCoordsAtQPs(int numPoints);
/** \brief
@@ -64,33 +68,53 @@ namespace AMDiS {
int numPoints,
const FastQuadrature *quadFast = NULL);
- /// Returns element normals at quadrature points.
+
+ /** \brief
+ * Returns element normals at quadrature points.
+ */
WorldVector<double> **getElementNormalAtQPs(int numPoints);
- ///
+ /** \brief
+ *
+ */
DimVec<WorldVector<double> > **getGrdLambdaAtQPs(int numPoints);
- /// Returns a QPInfo instance for the given quadrature.
+ /** \brief
+ * Returns a QPInfo instance for the given quadrature.
+ */
static QPInfo *provideQPInfo(const Quadrature*, const FastQuadrature*);
- /// Deletes the QPInfo instance for the given quadrature.
+ /** \brief
+ * Deletes the QPInfo instance for the given quadrature.
+ */
static void clearQPInfo(const Quadrature*, const FastQuadrature*);
- /// Deletes all QPInfo instances.
+ /** \brief
+ * Deletes all QPInfo instances.
+ */
static void clearAllQPInfos();
protected:
- /// Constructor. Called by \ref provideQPInfo().
+ /** \brief
+ * Constructor. Called by \ref provideQPInfo().
+ */
QPInfo(const Quadrature*);
- /// Destructor. Called by \ref clearQPInfo() and \ref clearAllQPInfos().
+ /** \brief
+ * Destructor. Called by \ref clearQPInfo() and \ref clearAllQPInfos().
+ */
~QPInfo();
protected:
- /// Structure which stores infos about one DOFVector.
+ /** \brief
+ * Structure which stores infos about one DOFVector.
+ */
class VecQPInfo
{
public:
+ /** \brief
+ * Constructor.
+ */
VecQPInfo()
: valAtQPs_(NULL),
grdAtQPs_(NULL),
@@ -100,62 +124,100 @@ namespace AMDiS {
D2NumPointsValid_(0)
{}
- /// Values at quadrature points
+ /** \brief
+ * Values at quadrature points
+ */
double *valAtQPs_;
- /// Gradients at quadrature points
+ /** \brief
+ * Gradients at quadrature points
+ */
WorldVector<double> *grdAtQPs_;
- /// D2 at quadrature points
+ /** \brief
+ * D2 at quadrature points
+ */
WorldMatrix<double> *D2AtQPs_;
- /// Valid flag for values
+ /** \brief
+ * valid flag for values
+ */
int valNumPointsValid_;
- /// Valid flag for gradients
+ /** \brief
+ * valid flag for gradients
+ */
int grdNumPointsValid_;
- /// Valid flag for D2
+ /** \brief
+ * valid flag for D2
+ */
bool D2NumPointsValid_;
};
- /// Quadrature of this QPInfo
+ /** \brief
+ * Quadrature of this QPInfo
+ */
const Quadrature *quadrature_;
- /// Set to \ref quadrature_->getNumPoints().
+ /** \brief
+ * Set to \ref quadrature_->getNumPoints().
+ */
int numPoints_;
- /// ElInfo of the current element
+ /** \brief
+ * ElInfo of the current element
+ */
const ElInfo *currentElInfo_;
- /// Coords at quadrature points
+ /** \brief
+ * Coords at quadrature points
+ */
WorldVector<double> *coordsAtQPs_;
- /// Valid flag for coords
+ /** \brief
+ * Valid flag for coords
+ */
int coordsNumPointsValid_;
- /// Map of all vector infos
+ /** \brief
+ * Map of all vector infos
+ */
std::map<const DOFVector<double>*, VecQPInfo*> vecQPInfos_;
- /// element normal at quadrature points (array of pointers)
+ /** \brief
+ * element normal at quadrature points (array of pointers)
+ */
WorldVector<double> **elementNormalAtQPs_;
- /// for constant values at all QPs (all entries point to same memory)
+ /** \brief
+ * for constant values at all QPs (all entries point to same memory)
+ */
WorldVector<double> **elementNormalConst_;
- /// valid flag for element normals
+ /** \brief
+ * valid flag for element normals
+ */
int elementNormalNumPointsValid_;
- /// gradient of barycentric coordinates at QPs (array of pointers)
+ /** \brief
+ * gradient of barycentric coordinates at QPs (array of pointers)
+ */
DimVec<WorldVector<double> > **grdLambdaAtQPs_;
- /// for constant values at all QPs (all entries point to same memory)
+ /** \brief
+ * for constant values at all QPs (all entries point to same memory)
+ */
DimVec<WorldVector<double> > **grdLambdaConst_;
- /// number of valid points of grdLambdaAtQPs_
+ /** \brief
+ * number of valid points of grdLambdaAtQPs_
+ */
int grdLambdaNumPointsValid_;
- /// Static map of all QPInfos. Used by \ref provideQPInfo().
+ /** \brief
+ * Static map of all QPInfos. Used by \ref provideQPInfo().
+ */
static std::map<const Quadrature*, QPInfo*> qpInfos_;
};
diff --git a/AMDiS/src/Recovery.cc b/AMDiS/src/Recovery.cc
index 8e05184c..f9a66903 100644
--- a/AMDiS/src/Recovery.cc
+++ b/AMDiS/src/Recovery.cc
@@ -501,23 +501,21 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
Mesh::FILL_COORDS |
Mesh::FILL_BOUND |
Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_LOCAL_INDICES;
+ Mesh::FILL_GRD_LAMBDA;
if (degree == 2 && dim > 1)
fill_flag |= Mesh::FILL_NEIGH | Mesh::FILL_OPP_COORDS;
TraverseStack stack;
- stack.addFeSpace(feSpace);
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, fill_flag);
+ ElInfo *el_info = stack.traverseFirst(mesh, -1, fill_flag);
- while (elInfo) { // traversing the mesh.
- const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
+ while (el_info) { // traversing the mesh.
+ const DegreeOfFreedom **dof = el_info->getElement()->getDOF();
int n_neighbors = 0; // counting interior vertices of element
for (int i = 0; i < n_vertices; i++) {
DegreeOfFreedom k = dof[i][preDOFs[VERTEX]];
- if (elInfo->getBoundary(VERTEX, i) == INTERIOR)
+ if (el_info->getBoundary(VERTEX, i) == INTERIOR)
interior_vertices[n_neighbors++] = k;
}
@@ -530,10 +528,10 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Setting world coordinates of node.
if (!(*struct_vec)[k].coords) {
(*struct_vec)[k].coords = new WorldVector<double>;
- *(*struct_vec)[k].coords = elInfo->getCoord(i);
+ *(*struct_vec)[k].coords = el_info->getCoord(i);
}
- if (elInfo->getBoundary(VERTEX, i) == INTERIOR) {
+ if (el_info->getBoundary(VERTEX, i) == INTERIOR) {
// Allocating memory for matrix and right hand side.
if (!(*struct_vec)[k].A) {
(*struct_vec)[k].A = new Matrix<double>(n_monomials, n_monomials);
@@ -551,20 +549,20 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Computing the integrals.
if (method)
- compute_integrals(uh, elInfo, &(*struct_vec)[k],
+ compute_integrals(uh, el_info, &(*struct_vec)[k],
f_vec, f_scal, aux_vec);
else if (gradient)
- compute_interior_sums(uh, elInfo, &(*struct_vec)[k], quad,
+ compute_interior_sums(uh, el_info, &(*struct_vec)[k], quad,
f_vec, f_scal, aux_vec);
else {
if (!(*struct_vec)[k].neighbors)
(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
if (degree == 2 && dim > 1)
- compute_sums_linear(uh, elInfo, &(*struct_vec)[k],
+ compute_sums_linear(uh, el_info, &(*struct_vec)[k],
i, pre_dofs, n_vertices);
else
- compute_node_sums(uh, elInfo, &(*struct_vec)[k],
+ compute_node_sums(uh, el_info, &(*struct_vec)[k],
pre_dofs, n_vertices, n_edges, n_faces);
}
} else { // Setting list of adjacent interior vertices.
@@ -585,7 +583,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
if (!(*struct_vec)[k].coords) {
// Setting world coordinates of node.
- elInfo->coordToWorld(*basis_fcts->getCoords(n_count),
+ el_info->coordToWorld(*basis_fcts->getCoords(n_count),
coordinates);
(*struct_vec)[k].coords = new WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -593,10 +591,10 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Setting list of adjacent interior vertices.
(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
- if (elInfo->getBoundary(EDGE, i) == INTERIOR)
+ if (el_info->getBoundary(EDGE, i) == INTERIOR)
for (int m = 0; m < 2; m++) {
int l = Global::getReferenceElement(dim)->getVertexOfEdge(i, m);
- if (elInfo->getBoundary(VERTEX, l) == INTERIOR)
+ if (el_info->getBoundary(VERTEX, l) == INTERIOR)
(*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
} else
for (int m = 0; m < n_neighbors; m++)
@@ -613,7 +611,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
if (!(*struct_vec)[k].coords) {
// Setting world coordinates of node.
- elInfo->coordToWorld(*basis_fcts->getCoords(n_count),
+ el_info->coordToWorld(*basis_fcts->getCoords(n_count),
coordinates);
(*struct_vec)[k].coords = new WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -621,10 +619,10 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Setting list of adjacent interior vertices.
(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
- if (elInfo->getBoundary(FACE, i) == INTERIOR)
+ if (el_info->getBoundary(FACE, i) == INTERIOR)
for (int m = 0; m < 3; m++) {
int l = Global::getReferenceElement(dim)->getVertexOfPosition(FACE, i, m);
- if (elInfo->getBoundary(VERTEX, l) == INTERIOR)
+ if (el_info->getBoundary(VERTEX, l) == INTERIOR)
(*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
}
else
@@ -640,7 +638,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
DegreeOfFreedom k = dof[n_vertices+n_edges+n_faces][preDOFs[CENTER] + j];
// Setting world coordinates of node.
- elInfo->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
+ el_info->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
(*struct_vec)[k].coords = new WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -653,7 +651,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
n_count++;
}
- elInfo = stack.traverseNext(elInfo);
+ el_info = stack.traverseNext(el_info);
}
}
@@ -886,11 +884,8 @@ Recovery::recovery(DOFVector<double> *uh,
// traverse mesh
TraverseStack stack;
- stack.addFeSpace(feSpace);
- Flag fillFlag =
- Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS |
- Mesh::FILL_LOCAL_INDICES;
+ Flag fillFlag =
+ Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) {
diff --git a/AMDiS/src/RecoveryEstimator.cc b/AMDiS/src/RecoveryEstimator.cc
index 4c28a587..7eac009f 100644
--- a/AMDiS/src/RecoveryEstimator.cc
+++ b/AMDiS/src/RecoveryEstimator.cc
@@ -94,7 +94,6 @@ namespace AMDiS {
// clear error indicators
TraverseStack stack;
- stack.addFeSpace(uh_->getFESpace());
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
elInfo->getElement()->setEstimation(0.0, row);
@@ -106,8 +105,7 @@ namespace AMDiS {
elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
- Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_LOCAL_INDICES);
+ Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
while (elInfo) {
Element *el = elInfo->getElement();
double det = elInfo->getDet();
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index 387ee37b..271a6ae9 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -102,15 +102,14 @@ namespace AMDiS {
est_max = 0.0;
est_t_sum = 0.0;
est_t_max = 0.0;
-
+
traverseFlag =
- Mesh::FILL_NEIGH |
- Mesh::FILL_COORDS |
- Mesh::FILL_OPP_COORDS |
- Mesh::FILL_BOUND |
- Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_DET |
- Mesh::FILL_LOCAL_INDICES |
+ Mesh::FILL_NEIGH |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_OPP_COORDS |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_DET |
Mesh::CALL_LEAF_EL;
neighInfo = mesh->createNewElInfo();
@@ -287,7 +286,7 @@ namespace AMDiS {
if (timestep && uhOld[system]) {
TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
- uhOld[system]->getLocalVector(elInfo, uhOldEl[system]);
+ uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
// === Compute time error. ===
@@ -446,7 +445,7 @@ namespace AMDiS {
if (matrix[system] == NULL || secondOrderTerms[system] == false)
continue;
- uh[system]->getLocalVector(elInfo, uhEl[system]);
+ uh[system]->getLocalVector(el, uhEl[system]);
uh[system]->getLocalVector(neigh, uhNeigh[system]);
for (int iq = 0; iq < nPointsSurface; iq++) {
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index 6a4ec41e..9df83ac1 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -22,10 +22,10 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim);
- coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
+ coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
// for all element sides
- for (int i = 0; i < dim + 1; i++) {
+ for (int i = 0; i < dim+1; i++) {
coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
DimVec<double>(dim,
DEFAULT_VALUE,
@@ -74,7 +74,7 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim);
- coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
+ coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
// for all element sides
for (int i = 0; i < dim + 1; i++) {
@@ -95,7 +95,7 @@ namespace AMDiS {
jOp->addZeroOrderTerm(new VecAtQP_ZOT(j, NULL));
neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
- for (int i=0; i < dim + 1; i++)
+ for (int i=0; i < dim+1; i++)
(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
delete jOp;
@@ -126,7 +126,7 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim);
- coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
+ coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
// for all element sides
for (int i = 0; i < dim + 1; i++) {
@@ -248,7 +248,10 @@ namespace AMDiS {
(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->getZeroOrderAssembler()->
initElement(elInfo);
- const double *uhAtQp = dv->getVecAtQPs(elInfo, quadrature, NULL, NULL);
+ const double *uhAtQp = dv->getVecAtQPs(elInfo,
+ quadrature,
+ NULL,
+ NULL);
std::vector<double> f(numPoints, 0.0);
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 9b848227..8c9cb62f 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -187,7 +187,7 @@ namespace AMDiS {
int i = info_stack[stack_used];
el = const_cast<Element*>(((i == 0) ? el->getFirstChild() : el->getSecondChild()));
info_stack[stack_used]++;
- fillElInfo(i);
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
@@ -279,7 +279,7 @@ namespace AMDiS {
int i = info_stack[stack_used];
info_stack[stack_used]++;
- fillElInfo(i);
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
@@ -347,7 +347,7 @@ namespace AMDiS {
int i = info_stack[stack_used];
info_stack[stack_used]++;
- fillElInfo(i);
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
@@ -412,12 +412,12 @@ namespace AMDiS {
while (elinfo_stack[stack_used]->getElement()->getFirstChild() &&
info_stack[stack_used] < 2) {
- if (stack_used >= stack_size - 1)
+ if (stack_used >= stack_size-1)
enlargeTraverseStack();
int i = info_stack[stack_used];
info_stack[stack_used]++;
- fillElInfo(i);
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
info_stack[stack_used] = 0;
}
@@ -482,7 +482,7 @@ namespace AMDiS {
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
i = 1 - neighbour;
- fillElInfo(i);
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
info_stack[stack_used] = i + 1;
stack_used++;
info_stack[stack_used] = 0;
@@ -540,11 +540,11 @@ namespace AMDiS {
elinfo2 = save_elinfo_stack[stack2_used];
el2 = elinfo2->getElement();
- if (stack_used >= stack_size - 1)
+ if (stack_used >= stack_size-1)
enlargeTraverseStack();
i = 2 - info_stack[stack_used];
- info_stack[stack_used] = i + 1;
- fillElInfo(i);
+ info_stack[stack_used] = i+1;
+ elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
info_stack[stack_used] = 0;
nb = 0;
@@ -558,8 +558,8 @@ namespace AMDiS {
if (nb < 2) { /* go down one level in hierarchy */
if (stack_used >= stack_size-1)
enlargeTraverseStack();
- info_stack[stack_used] = 2 - nb;
- fillElInfo(1 - nb);
+ info_stack[stack_used] = 2-nb;
+ elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
stack_used++;
info_stack[stack_used] = 0;
elinfo = elinfo_stack[stack_used];
@@ -594,10 +594,10 @@ namespace AMDiS {
nb = 2;
}
- info_stack[stack_used] = i + 1;
- if (stack_used >= stack_size - 1)
+ info_stack[stack_used] = i+1;
+ if (stack_used >= stack_size-1)
enlargeTraverseStack();
- fillElInfo(i);
+ elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
info_stack[stack_used] = 0;
@@ -630,11 +630,11 @@ namespace AMDiS {
if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
MSG(" looking for neighbour %d of element %d at %p\n",
- neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>(old_elinfo->getElement()));
+ neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>( old_elinfo->getElement()));
MSG(" originally: neighbour %d of element %d at %p\n",
- sav_neighbour, sav_index, reinterpret_cast<void*>(old_elinfo->getElement()));
+ sav_neighbour, sav_index, reinterpret_cast<void*>( old_elinfo->getElement()));
MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
- elinfo->getElement()->getIndex(), reinterpret_cast<void*>(elinfo->getElement()),
+ elinfo->getElement()->getIndex(), reinterpret_cast<void*>( elinfo->getElement()),
opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
("didn't succeed !?!?!?\n");
@@ -679,16 +679,15 @@ namespace AMDiS {
TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");
elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
- el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo_stack[stack_used]->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo_stack[stack_used]->getElement()));
sav_index = el->getIndex();
sav_el = el;
/* first, goto to leaf level, if necessary... */
if (!(el->isLeaf()) && (neighbour < 2)) {
- if (stack_used >= static_cast<int>(elinfo_stack.size()) - 1)
- enlargeTraverseStack();
+ if (stack_used >= static_cast<int>( elinfo_stack.size())-1) enlargeTraverseStack();
i = 1 - neighbour;
- fillElInfo(i);
+ elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
info_stack[stack_used] = i + 1;
stack_used++;
neighbour = 2;
@@ -737,7 +736,7 @@ namespace AMDiS {
stack2_used = 1;
}
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
i = traverse_mel->getOppVertex(nb);
traverse_mel = traverse_mel->getNeighbour(nb);
@@ -746,7 +745,7 @@ namespace AMDiS {
nb = i;
stack_used = 1;
- elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>(traverse_mel));
+ elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
info_stack[stack_used] = 0;
} else { /* goto other child */
@@ -756,17 +755,17 @@ namespace AMDiS {
stack2_used++; /* go down one level in OLD hierarchy */
}
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
if (stack_used >= stack_size-1) {
enlargeTraverseStack();
}
i = 2 - info_stack[stack_used];
- info_stack[stack_used] = i + 1;
- fillElInfo(i);
+ info_stack[stack_used] = i+1;
+ elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
- nb = 1 - i;
+ nb = 1-i;
}
/****************************************************************************/
@@ -776,15 +775,15 @@ namespace AMDiS {
/****************************************************************************/
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
while (el->getFirstChild()) {
if (nb < 2) { /* go down one level in hierarchy */
- if (stack_used >= stack_size - 1)
+ if (stack_used >= stack_size-1)
enlargeTraverseStack();
- fillElInfo(1 - nb);
- info_stack[stack_used] = 2 - nb;
+ elinfo_stack[stack_used+1]->fillElInfo(1-nb, elinfo_stack[stack_used]);
+ info_stack[stack_used] = 2-nb;
stack_used++;
nb = 2;
}
@@ -795,23 +794,25 @@ namespace AMDiS {
/* use child i, neighbour of el2->child[nb-1] */
i = 2 - save_info_stack[stack2_used];
TEST_EXIT_DBG(i < 2)("invalid OLD refinement?");
- info_stack[stack_used] = i + 1;
- fillElInfo(i);
+ info_stack[stack_used] = i+1;
+ elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
nb = i;
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
stack2_used++;
- if (save_stack_used > stack2_used)
- stack2_used++; /* go down one level in OLD hierarchy */
+ if (save_stack_used > stack2_used) {
+ stack2_used++; /* go down one level in OLD hierarchy */
+ }
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
- } else { /* now we're done... */
+ }
+ else { /* now we're done... */
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
}
}
@@ -873,4 +874,5 @@ namespace AMDiS {
elInfo.setRefinementPath(rPath);
elInfo.setRefinementPathLength(levelDif);
}
+
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 19479de1..b4ded5bd 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -38,8 +38,6 @@
#include "ElInfo.h"
#include "ElInfoStack.h"
#include "OpenMP.h"
-#include "Mesh.h"
-#include "BasisFunction.h"
#include "AMDiS_fwd.h"
namespace AMDiS {
@@ -138,42 +136,6 @@ namespace AMDiS {
return stack_used;
}
- /** \brief
- * Adds new FE space to \ref feSpaces. Is used, when local indices of leaf
- * elements should be computed during mesh traverse.
- */
- void addFeSpace(const FiniteElemSpace *feSpace)
- {
- feSpaces.insert(feSpace);
- }
-
- private:
- /** \brief
- * Fills the element info data of a new ElInfo object on the top of the stack.
- *
- * \param[in] iChild Defines if the new element on the stack is the left or the
- * right children of its parent element.
- */
- inline void fillElInfo(int iChild)
- {
- FUNCNAME("Traverse::fillElInfo()");
-
- ElInfo *elInfo = elinfo_stack[stack_used + 1];
- elInfo->fillElInfo(iChild, elinfo_stack[stack_used]);
-
- if (traverse_fill_flag.isSet(Mesh::FILL_LOCAL_INDICES) &&
- elInfo->getElement()->isLeaf()) {
-
- TEST_EXIT_DBG(feSpaces.size() > 0)("There are no FE space defined!\n");
-
- for (std::set<const FiniteElemSpace*>::iterator it = feSpaces.begin();
- it != feSpaces.end(); ++it)
- (*it)->getBasisFcts()->getLocalIndices(elInfo->getElement(),
- (*it)->getAdmin(),
- elInfo->getLocalIndices(*it));
- }
- }
-
private:
/// Enlargement of the stack
void enlargeTraverseStack();
@@ -262,12 +224,6 @@ namespace AMDiS {
///
int id;
- /** \brief
- * If local indices should be field for leaf elements, here all used FE spaces
- * are stored, for which the local indices should be computed.
- */
- std::set<const FiniteElemSpace*> feSpaces;
-
/** \brief
* Is used for parallel mesh traverse. The thread with the id
* myThreadId is only allowed to access coarse elements, which id
diff --git a/AMDiS/src/TraverseParallel.h b/AMDiS/src/TraverseParallel.h
index aa994797..44046eaa 100644
--- a/AMDiS/src/TraverseParallel.h
+++ b/AMDiS/src/TraverseParallel.h
@@ -44,12 +44,6 @@ namespace AMDiS {
ElInfo* traverseNext(ElInfo* elInfoOld);
- void addFeSpace(const FiniteElemSpace* feSpace)
- {
- for (unsigned int i = 0; i < stacks.size(); i++)
- stacks[i].addFeSpace(feSpace);
- }
-
private:
/// Number of threads using the stack in parallel.
int nThreads;
--
GitLab