Commit d1ff5073 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Argh, some error, go back to revision 1027.

parent bfa16bad
...@@ -252,11 +252,12 @@ namespace AMDiS { ...@@ -252,11 +252,12 @@ namespace AMDiS {
{ {
FUNCNAME("Assembler::matVecAssemble()"); FUNCNAME("Assembler::matVecAssemble()");
Element *el = elInfo->getElement();
double *uhOldLoc = new double[nRow]; 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); set_to_zero(elementMatrix);
calculateElementMatrix(elInfo, elementMatrix); calculateElementMatrix(elInfo, elementMatrix);
} }
......
...@@ -275,7 +275,7 @@ namespace AMDiS { ...@@ -275,7 +275,7 @@ namespace AMDiS {
virtual void coarseRestr(DOFVector<WorldVector<double> >*, RCNeighbourList*, int) 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, virtual const DegreeOfFreedom *getLocalIndices(const Element *el,
const DOFAdmin *admin, const DOFAdmin *admin,
DegreeOfFreedom *dofPtr) const DegreeOfFreedom *dofPtr) const
...@@ -283,11 +283,17 @@ namespace AMDiS { ...@@ -283,11 +283,17 @@ namespace AMDiS {
return NULL; return NULL;
} }
/// Calculates local dof indices of the element for the given DOF admin. inline void getLocalIndices(const Element *el,
virtual void getLocalIndices(const Element *el,
const DOFAdmin *admin, const DOFAdmin *admin,
std::vector<DegreeOfFreedom> &indices) const 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, virtual void getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin, const DOFAdmin *admin,
......
...@@ -15,6 +15,7 @@ namespace AMDiS { ...@@ -15,6 +15,7 @@ namespace AMDiS {
BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace) BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
{ {
localBounds.resize(omp_get_overall_max_threads()); localBounds.resize(omp_get_overall_max_threads());
dofIndices.resize(omp_get_overall_max_threads());
allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber(); allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber();
for (int i = 0; i < static_cast<int>(localBounds.size()); i++) for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds]; localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
...@@ -25,6 +26,7 @@ namespace AMDiS { ...@@ -25,6 +26,7 @@ namespace AMDiS {
localBCs = bm.localBCs; localBCs = bm.localBCs;
allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds; allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds;
localBounds.resize(bm.localBounds.size()); localBounds.resize(bm.localBounds.size());
dofIndices.resize(bm.localBounds.size());
for (int i = 0; i < static_cast<int>(localBounds.size()); i++) for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds]; localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
} }
...@@ -47,27 +49,22 @@ namespace AMDiS { ...@@ -47,27 +49,22 @@ namespace AMDiS {
return result; return result;
} }
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
DOFVectorBase<double> *vec) DOFVectorBase<double> *vec)
{ {
FUNCNAME("BoundaryManager::fillBoundaryConditions()"); if (localBCs.size() > 0) {
if (localBCs.size() <= 0)
return;
const FiniteElemSpace *feSpace = vec->getFESpace(); const FiniteElemSpace *feSpace = vec->getFESpace();
std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts(); const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber(); int nBasFcts = basisFcts->getNumber();
dofVec.resize(nBasFcts);
// get boundaries of all DOFs // get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()]; BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound); basisFcts->getBound(elInfo, localBound);
// get dof indices // get dof indices
std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace); basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
("Local index vector is too small!\n");
// apply non dirichlet boundary conditions // apply non dirichlet boundary conditions
for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it) for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
...@@ -81,27 +78,25 @@ namespace AMDiS { ...@@ -81,27 +78,25 @@ namespace AMDiS {
(*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0], (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
localBound, nBasFcts); localBound, nBasFcts);
} }
}
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, DOFMatrix *mat) void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, DOFMatrix *mat)
{ {
FUNCNAME("BoundaryManager::fillBoundaryConditions()");
if (localBCs.size() <= 0) if (localBCs.size() <= 0)
return; return;
const FiniteElemSpace *feSpace = mat->getRowFESpace(); const FiniteElemSpace *feSpace = mat->getRowFESpace();
std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts(); const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber(); int nBasFcts = basisFcts->getNumber();
dofVec.resize(nBasFcts);
// get boundaries of all DOFs // get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()]; BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound); basisFcts->getBound(elInfo, localBound);
// get dof indices // get dof indices
std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace); basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
("Local index vector is too small!\n");
// apply non dirichlet boundary conditions // apply non dirichlet boundary conditions
for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it) for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
......
...@@ -125,6 +125,9 @@ namespace AMDiS { ...@@ -125,6 +125,9 @@ namespace AMDiS {
/// Temporary thread-safe variable for functions fillBoundaryconditions. /// Temporary thread-safe variable for functions fillBoundaryconditions.
std::vector<BoundaryType*> localBounds; std::vector<BoundaryType*> localBounds;
/// Temporary thread-safe variable for functions fillBoundaryconditions.
std::vector<std::vector<DegreeOfFreedom> > dofIndices;
/** \brief /** \brief
* Stores the number of byte that were allocated in the constructor for * Stores the number of byte that were allocated in the constructor for
* each localBounds value. Is used to free the memory in the destructor. * each localBounds value. Is used to free the memory in the destructor.
......
...@@ -168,12 +168,9 @@ namespace AMDiS { ...@@ -168,12 +168,9 @@ namespace AMDiS {
// === Get indices mapping from local to global matrix indices. === // === Get indices mapping from local to global matrix indices. ===
std::vector<DegreeOfFreedom> &rowIndices2 = rowElInfo->getLocalIndices(rowFESpace); rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
std::vector<DegreeOfFreedom> &colIndices2 = rowIndices2; rowFESpace->getAdmin(),
rowIndices);
// rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
// rowFESpace->getAdmin(),
// rowIndices);
if (rowFESpace == colFESpace) { if (rowFESpace == colFESpace) {
colIndices = rowIndices; colIndices = rowIndices;
} else { } else {
...@@ -191,8 +188,21 @@ namespace AMDiS { ...@@ -191,8 +188,21 @@ namespace AMDiS {
using namespace mtl; 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++) { for (int i = 0; i < nRow; i++) {
DegreeOfFreedom row = rowIndices2[i]; DegreeOfFreedom row = rowIndices[i];
BoundaryCondition *condition = BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL; bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
...@@ -200,15 +210,20 @@ namespace AMDiS { ...@@ -200,15 +210,20 @@ namespace AMDiS {
if (condition && condition->isDirichlet()) { if (condition && condition->isDirichlet()) {
if (condition->applyBoundaryCondition()) { if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
if ((*rankDofs)[row]) if ((*rankDofs)[rowIndices[i]])
applyDBCs.insert(static_cast<int>(row)); applyDBCs.insert(static_cast<int>(row));
#else #else
applyDBCs.insert(static_cast<int>(row)); applyDBCs.insert(static_cast<int>(row));
#endif #endif
} }
} else { } else {
for (int j = 0; j < nCol; j++) for (int j = 0; j < nCol; j++) {
ins[row][colIndices2[j]] += elMat[i][j]; DegreeOfFreedom col = colIndices[j];
double entry = elMat[i][j];
// std::cout << "ADD at " << row << " " << col << std::endl;
ins[row][col] += entry;
}
} }
} }
} }
......
...@@ -108,16 +108,13 @@ namespace AMDiS { ...@@ -108,16 +108,13 @@ namespace AMDiS {
// traverse mesh // traverse mesh
std::vector<bool> visited(getUsedSize(), false); std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack; TraverseStack stack;
stack.addFeSpace(feSpace); Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
Flag fillFlag =
Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_COORDS | Mesh::FILL_LOCAL_INDICES;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag); ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) { while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda(); const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
getLocalVector(elInfo, localUh); getLocalVector(elInfo->getElement(), localUh);
int localDOFNr = 0; int localDOFNr = 0;
for (int i = 0; i < nNodes; i++) { // for all nodes for (int i = 0; i < nNodes; i++) { // for all nodes
...@@ -180,11 +177,9 @@ namespace AMDiS { ...@@ -180,11 +177,9 @@ namespace AMDiS {
// traverse mesh // traverse mesh
Mesh *mesh = feSpace->getMesh(); Mesh *mesh = feSpace->getMesh();
TraverseStack stack; TraverseStack stack;
stack.addFeSpace(feSpace);
ElInfo *elInfo = stack.traverseFirst(mesh, -1, ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
Mesh::FILL_LOCAL_INDICES);
double *localUh = new double[basFcts->getNumber()]; double *localUh = new double[basFcts->getNumber()];
...@@ -192,7 +187,7 @@ namespace AMDiS { ...@@ -192,7 +187,7 @@ namespace AMDiS {
double det = elInfo->getDet(); double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda(); const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
getLocalVector(elInfo, localUh); getLocalVector(elInfo->getElement(), localUh);
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd); basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
for (int i = 0; i < dim + 1; i++) { for (int i = 0; i < dim + 1; i++) {
...@@ -254,7 +249,7 @@ namespace AMDiS { ...@@ -254,7 +249,7 @@ namespace AMDiS {
} }
double *localVec = localVectors[myRank]; double *localVec = localVectors[myRank];
getLocalVector(elInfo, localVec); getLocalVector(elInfo->getElement(), localVec);
DimVec<double> &grd1 = *grdTmp[myRank]; DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim); int parts = Global::getGeo(PARTS, dim);
...@@ -339,7 +334,7 @@ namespace AMDiS { ...@@ -339,7 +334,7 @@ namespace AMDiS {
} }
double *localVec = localVectors[myRank]; double *localVec = localVectors[myRank];
getLocalVector(largeElInfo, localVec); getLocalVector(largeElInfo->getElement(), localVec);
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree()); mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
...@@ -396,6 +391,8 @@ namespace AMDiS { ...@@ -396,6 +391,8 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts()) TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions"); ("invalid basis functions");
Element *el = elInfo->getElement();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
int dow = Global::getGeo(WORLD); int dow = Global::getGeo(WORLD);
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
...@@ -415,7 +412,7 @@ namespace AMDiS { ...@@ -415,7 +412,7 @@ namespace AMDiS {
} }
double *localVec = localVectors[myRank]; double *localVec = localVectors[myRank];
getLocalVector(elInfo, localVec); getLocalVector(el, localVec);
DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0); DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
int parts = Global::getGeo(PARTS, dim); int parts = Global::getGeo(PARTS, dim);
...@@ -779,7 +776,7 @@ namespace AMDiS { ...@@ -779,7 +776,7 @@ namespace AMDiS {
} }
double *localVec = localVectors[omp_get_thread_num()]; double *localVec = localVectors[omp_get_thread_num()];
getLocalVector(largeElInfo, localVec); getLocalVector(largeElInfo->getElement(), localVec);
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree()); mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
......
...@@ -67,9 +67,7 @@ namespace AMDiS { ...@@ -67,9 +67,7 @@ namespace AMDiS {
* For the given element, this function returns an array of all DOFs of this * For the given element, this function returns an array of all DOFs of this
* DOFVector that are defined on this element. * DOFVector that are defined on this element.
*/ */
const T *getLocalVector(const Element *el, T* localVec) const; virtual const T *getLocalVector(const Element *el, T* localVec) const;
const T *getLocalVector(const ElInfo *elInfo, T* localVec) const;
/** \brief /** \brief
* Evaluates the DOF vector at a set of quadrature points defined on the * Evaluates the DOF vector at a set of quadrature points defined on the
......
...@@ -147,9 +147,9 @@ namespace AMDiS { ...@@ -147,9 +147,9 @@ namespace AMDiS {
{ {
FUNCNAME("DOFVector::addElementVector()"); FUNCNAME("DOFVector::addElementVector()");
std::vector<DegreeOfFreedom> &indices = elInfo->getLocalIndices(feSpace); std::vector<DegreeOfFreedom> indices(nBasFcts);
TEST_EXIT_DBG(static_cast<int>(indices.size()) == nBasFcts) feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
("Local index vector is too small!\n"); indices);
for (DegreeOfFreedom i = 0; i < nBasFcts; i++) { for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
BoundaryCondition *condition = BoundaryCondition *condition =
...@@ -454,11 +454,9 @@ namespace AMDiS { ...@@ -454,11 +454,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints(); int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints); std::vector<T> uh_vec(nPoints);
TraverseStack stack; TraverseStack stack;
stack.addFeSpace(this->feSpace);
ElInfo *elInfo = ElInfo *elInfo =
stack.traverseFirst(mesh, -1, stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) { while (elInfo) {
double det = elInfo->getDet(); double det = elInfo->getDet();
double normT = 0.0; double normT = 0.0;
...@@ -492,11 +490,9 @@ namespace AMDiS { ...@@ -492,11 +490,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints(); int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints); std::vector<T> uh_vec(nPoints);
TraverseStack stack; TraverseStack stack;
stack.addFeSpace(this->feSpace);
ElInfo *elInfo = ElInfo *elInfo =
stack.traverseFirst(mesh, -1, stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) { while (elInfo) {
double det = elInfo->getDet(); double det = elInfo->getDet();
double normT = 0.0; double normT = 0.0;
...@@ -530,11 +526,9 @@ namespace AMDiS { ...@@ -530,11 +526,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints(); int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints); std::vector<T> uh_vec(nPoints);
TraverseStack stack; TraverseStack stack;
stack.addFeSpace(this->feSpace);
ElInfo *elInfo = ElInfo *elInfo =
stack.traverseFirst(mesh, -1, stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) { while (elInfo) {
double det = elInfo->getDet(); double det = elInfo->getDet();
double normT = 0.0; double normT = 0.0;
...@@ -992,49 +986,6 @@ namespace AMDiS { ...@@ -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> template<typename T>
const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad, const Quadrature *quad,
...@@ -1053,6 +1004,7 @@ namespace AMDiS { ...@@ -1053,6 +1004,7 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts()) TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions"); ("invalid basis functions");
Element *el = elInfo->getElement();
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad; const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
...@@ -1076,7 +1028,7 @@ namespace AMDiS { ...@@ -1076,7 +1028,7 @@ namespace AMDiS {
} }
T *localVec = localVectors[omp_get_thread_num()]; T *localVec = localVectors[omp_get_thread_num()];
getLocalVector(elInfo, localVec); getLocalVector(el, localVec);
for (int i = 0; i < nPoints; i++) { for (int i = 0; i < nPoints; i++) {
result[i] = 0.0; result[i] = 0.0;
...@@ -1141,11 +1093,9 @@ namespace AMDiS { ...@@ -1141,11 +1093,9 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints(); int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints); std::vector<T> uh_vec(nPoints);
TraverseStack stack; TraverseStack stack;
stack.addFeSpace(this->feSpace);
ElInfo *elInfo = ElInfo *elInfo =
stack.traverseFirst(mesh, -1, stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) { while (elInfo) {
double det = elInfo->getDet(); double det = elInfo->getDet();
double normT = 0.0; double normT = 0.0;
......