From a6ae4ebcf8f5ec01d04a868e10d14616c4a88c41 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 29 Oct 2009 07:47:31 +0000
Subject: [PATCH] Fixed a bug when using openmp parallelization for assemlbing
matrices.
---
AMDiS/src/BasisFunction.cc | 1 +
AMDiS/src/DOFMatrix.cc | 2 +-
AMDiS/src/DOFVector.cc | 71 +++++++++++++++----------------
AMDiS/src/DOFVector.hh | 26 +++++++----
AMDiS/src/ElementFunction.h | 15 ++++---
AMDiS/src/Lagrange.cc | 9 ++--
AMDiS/src/Operator.cc | 41 +++++++++---------
AMDiS/src/Operator.h | 45 +++++++-------------
AMDiS/src/ProblemVec.cc | 2 +-
AMDiS/src/Quadrature.cc | 23 ++++++----
AMDiS/src/SecondOrderAssembler.cc | 22 ++++++----
AMDiS/src/SubAssembler.cc | 7 ++-
12 files changed, 133 insertions(+), 131 deletions(-)
diff --git a/AMDiS/src/BasisFunction.cc b/AMDiS/src/BasisFunction.cc
index bf1cfa7e..3ab8ed8f 100644
--- a/AMDiS/src/BasisFunction.cc
+++ b/AMDiS/src/BasisFunction.cc
@@ -140,6 +140,7 @@ namespace AMDiS {
return ((*val));
}
+
int BasisFunction::getNumberOfDOFs(int i) const
{
return (*nDOF)[i];
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 4f60974c..bd60a0e2 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -18,6 +18,7 @@
namespace AMDiS {
using namespace mtl;
+ using boost::lexical_cast;
DOFMatrix *DOFMatrix::traversePtr = NULL;
@@ -186,7 +187,6 @@ namespace AMDiS {
}
}
-
for (int i = 0; i < nRow; i++) {
DegreeOfFreedom row = rowIndices[i];
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index fb0740e0..78f4c2b5 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -71,11 +71,8 @@ namespace AMDiS {
}
Mesh *mesh = feSpace->getMesh();
-
int dim = mesh->getDim();
-
const BasisFunction *basFcts = feSpace->getBasisFcts();
-
DOFAdmin *admin = feSpace->getAdmin();
// count number of nodes and dofs per node
@@ -99,30 +96,24 @@ namespace AMDiS {
numNodes += numPositionNodes;
}
- // TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
- // ("invalid number of nodes\n");
-
TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
("number of dofs != number of basis functions\n");
- for (int i = 0; i < numDOFs; i++) {
- bary.push_back(basFcts->getCoords(i));
- }
+ for (int i = 0; i < numDOFs; i++)
+ bary.push_back(basFcts->getCoords(i));
+
+ double *localUh = new double[basFcts->getNumber()];
// traverse mesh
std::vector<bool> visited(getUsedSize(), false);
-
TraverseStack stack;
-
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 double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ getLocalVector(elInfo->getElement(), localUh);
int localDOFNr = 0;
for (int i = 0; i < numNodes; i++) { // for all nodes
@@ -141,6 +132,8 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
+ delete [] localUh;
+
return result;
}
@@ -150,6 +143,10 @@ namespace AMDiS {
{
FUNCNAME("DOFVector<double>::getRecoveryGradient()");
+#ifdef _OPENMP
+ ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
+#endif
+
// define result vector
static DOFVector<WorldVector<double> > *vec = NULL;
@@ -160,9 +157,9 @@ namespace AMDiS {
delete vec;
vec = NULL;
}
- if (!vec) {
+ if (!vec)
vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
- }
+
result = vec;
}
@@ -183,11 +180,13 @@ namespace AMDiS {
Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
+ double *localUh = new double[basFcts->getNumber()];
+
while (elInfo) {
double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
- const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ getLocalVector(elInfo->getElement(), localUh);
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
for (int i = 0; i < dim + 1; i++) {
@@ -199,14 +198,14 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
+ delete [] localUh;
+
DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
- for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt) {
- if (*volIt != 0.0) {
+ for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
+ if (*volIt != 0.0)
*grdIt *= 1.0 / (*volIt);
- }
- }
return result;
}
@@ -242,9 +241,9 @@ namespace AMDiS {
delete [] grd;
grd = new WorldVector<double>[nPoints];
- for (int i = 0; i < nPoints; i++) {
+ for (int i = 0; i < nPoints; i++)
grd[i].set(0.0);
- }
+
result = grd;
}
@@ -612,6 +611,10 @@ namespace AMDiS {
{
FUNCNAME("DOFVector<double>::getGradient()");
+#ifdef _OPENMP
+ ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
+#endif
+
Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim();
int dow = Global::getGeo(WORLD);
@@ -660,15 +663,11 @@ namespace AMDiS {
numNodes += numPositionNodes;
}
-// TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
-// ("invalid number of nodes\n");
-
TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
("number of dofs != number of basis functions\n");
- for (int i = 0; i < numDOFs; i++) {
- bary.push_back(basFcts->getCoords(i));
- }
+ for (int i = 0; i < numDOFs; i++)
+ bary.push_back(basFcts->getCoords(i));
// traverse mesh
std::vector<bool> visited(getUsedSize(), false);
@@ -691,9 +690,8 @@ namespace AMDiS {
if (!visited[dofIndex]) {
basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
- for (int k = 0; k < dow; k++) {
- (*(*result)[k])[dofIndex] = grd[k];
- }
+ for (int k = 0; k < dow; k++)
+ (*(*result)[k])[dofIndex] = grd[k];
visited[dofIndex] = true;
}
@@ -779,9 +777,8 @@ namespace AMDiS {
if (localvec)
delete [] localvec;
localvec = new double[nPoints];
- for (int i = 0; i < nPoints; i++) {
- localvec[i] = 0.0;
- }
+ for (int i = 0; i < nPoints; i++)
+ localvec[i] = 0.0;
result = localvec;
}
@@ -795,9 +792,9 @@ namespace AMDiS {
result[i] = 0.0;
for (int j = 0; j < nBasFcts; j++) {
double val = 0.0;
- for (int k = 0; k < nBasFcts; k++) {
+ for (int k = 0; k < nBasFcts; k++)
val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
- }
+
result[i] += localVec[j] * val;
}
}
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 198715e6..237c3503 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -962,6 +962,9 @@ namespace AMDiS {
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];
@@ -1002,35 +1005,40 @@ namespace AMDiS {
Element *el = elInfo->getElement();
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
- int numPoints = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
int nBasFcts = basFcts->getNumber();
- int j;
static T *localvec = NULL;
T *result;
if (vecAtQPs) {
result = vecAtQPs;
} else {
+#ifdef _OPENMP
+ ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
+#endif
if (localvec)
delete [] localvec;
- localvec = new T[numPoints];
- for (int i = 0; i < numPoints; i++)
+ localvec = new T[nPoints];
+ for (int i = 0; i < nPoints; i++)
localvec[i] = 0.0;
result = localvec;
}
- const T *localVec = getLocalVector(el, NULL);
+ T *localVec = new T[nBasFcts];
+ getLocalVector(el, localVec);
- for (int i = 0; i < numPoints; i++) {
- for (result[i] = j = 0; j < nBasFcts; j++) {
+ for (int i = 0; i < nPoints; i++) {
+ result[i] = 0.0;
+ for (int j = 0; j < nBasFcts; j++)
result[i] += localVec[j] *
(quadFast ?
(quadFast->getPhi(i, j)) :
- ((*(basFcts->getPhi(j)))(quad->getLambda(i))));
- }
+ ((*(basFcts->getPhi(j)))(quad->getLambda(i))));
}
+ delete [] localVec;
+
return const_cast<const T*>(result);
}
diff --git a/AMDiS/src/ElementFunction.h b/AMDiS/src/ElementFunction.h
index 5a6232cb..53fdf2e0 100644
--- a/AMDiS/src/ElementFunction.h
+++ b/AMDiS/src/ElementFunction.h
@@ -83,24 +83,25 @@ namespace AMDiS {
{
public:
/// constructor.
- ElementFunctionDOFVec(const DOFVector<T> *dofVector)
+ ElementFunctionDOFVec(const DOFVector<T> *vec)
: ElementFunction<T>(),
- dofVector_(dofVector)
+ dofVector(vec)
{}
/// evaluation at given coordinates.
T operator()(const DimVec<double>& bary) const
{
- static T t;
- const T* localVec =
- dofVector_->getLocalVector(this->elInfo_->getElement(), NULL);
- t = dofVector_->getFESpace()->getBasisFcts()->evalUh(bary, localVec);
+ T* localVec = new T[dofVector->getFESpace()->getBasisFcts()->getNumber()];
+ dofVector->getLocalVector(this->elInfo_->getElement(), localVec);
+ T t = dofVector->getFESpace()->getBasisFcts()->evalUh(bary, localVec);
+
+ delete [] localVec;
return t;
}
protected:
/// DOFVector to be interpolated.
- const DOFVector<T> *dofVector_;
+ const DOFVector<T> *dofVector;
};
}
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index fa06efd6..4a2757bb 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -2774,11 +2774,11 @@ namespace AMDiS {
DegreeOfFreedom pdof[4], dof9;
const DOFAdmin *admin;
- if (n < 1) return;
- el = list->getElement(0);
+ if (n < 1)
+ return;
+ el = list->getElement(0);
admin = drv->getFESpace()->getAdmin();
-
basFct->getLocalIndices(el, admin, pdof);
/****************************************************************************/
@@ -2824,7 +2824,8 @@ namespace AMDiS {
(*drv)[cdof[5]] + 0.9375*(*drv)[cdof[6]] + 0.1875*(*drv)[cdof[9]];
(*drv)[pdof[9]] += 0.75*(*drv)[cdof[9]];
- if (n <= 1) return;
+ if (n <= 1)
+ return;
/****************************************************************************/
/* adjust the value on the neihgbour */
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index af4974cc..97053258 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -6,8 +6,13 @@
#include "Quadrature.h"
#include "OpenMP.h"
+#include "boost/lexical_cast.hpp"
+#include <fstream>
+
namespace AMDiS {
+ using boost::lexical_cast;
+
const Flag OperatorTerm::PW_CONST = 1;
const Flag OperatorTerm::SYMMETRIC = 2;
@@ -71,7 +76,7 @@ namespace AMDiS {
{
FUNCNAME("OperatorTerm::getVectorAtQPs()");
- TEST_EXIT(elInfo->getMesh() == vec->getFESpace()->getMesh())
+ TEST_EXIT_DBG(elInfo->getMesh() == vec->getFESpace()->getMesh())
("There is something wrong!\n");
return subAssembler->getVectorAtQPs(vec, elInfo, quad);
@@ -94,23 +99,20 @@ namespace AMDiS {
// Both elements have the same size, so we can use the simple procedure
// to determine the vecAtQPs.
- if (vec->getFESpace()->getMesh() == smallElInfo->getMesh()) {
+ if (vec->getFESpace()->getMesh() == smallElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
- } else {
- return subAssembler->getVectorAtQPs(vec, largeElInfo, quad);
- }
+ else
+ return subAssembler->getVectorAtQPs(vec, largeElInfo, quad);
} else {
// The two elements are different. If the vector is defined on the mesh of the
// small element, we can still use the simple procedure to determine the vecAtQPs.
- if (vec->getFESpace()->getMesh() == largeElInfo->getMesh()) {
+ if (vec->getFESpace()->getMesh() == largeElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, largeElInfo, quad);
- } else {
+ else
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
- }
-
}
}
@@ -165,7 +167,7 @@ namespace AMDiS {
const WorldMatrix<double>& matrix,
DimMat<double>& LALt,
bool symm,
- double factor)
+ double factor) const
{
int j, k, l;
const int dimOfWorld = Global::getGeo(WORLD);
@@ -219,7 +221,7 @@ namespace AMDiS {
void OperatorTerm::l1lt(const DimVec<WorldVector<double> >& Lambda,
DimMat<double>& LALt,
- double factor)
+ double factor) const
{
const int dimOfWorld = Global::getGeo(WORLD);
int dim = LALt.getNumRows() - 1;
@@ -342,15 +344,14 @@ namespace AMDiS {
#pragma omp critical (initAssembler)
#endif
{
- if (optimized) {
- assembler[rank] = new OptimizedAssembler(this,
- quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFESpace, colFESpace);
- } else {
- assembler[rank] = new StandardAssembler(this,
- quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFESpace, colFESpace);
- }
+ if (optimized)
+ assembler[rank] =
+ new OptimizedAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
+ rowFESpace, colFESpace);
+ else
+ assembler[rank] =
+ new StandardAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
+ rowFESpace, colFESpace);
}
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index b456bd31..1094f876 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -144,11 +144,11 @@ namespace AMDiS {
protected:
/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$.
- static void lalt(const DimVec<WorldVector<double> >& Lambda,
- const WorldMatrix<double>& matrix,
- DimMat<double>& LALt,
- bool symm,
- double factor);
+ void lalt(const DimVec<WorldVector<double> >& Lambda,
+ const WorldMatrix<double>& matrix,
+ DimMat<double>& LALt,
+ bool symm,
+ double factor) const;
/** \brief
* Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for \f$ A \f$
@@ -161,9 +161,9 @@ namespace AMDiS {
double factor);
/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
- static void l1lt(const DimVec<WorldVector<double> >& Lambda,
- DimMat<double>& LALt,
- double factor);
+ void l1lt(const DimVec<WorldVector<double> >& Lambda,
+ DimMat<double>& LALt,
+ double factor) const;
/// Evaluation of \f$ \Lambda \cdot b\f$.
static inline void lb(const DimVec<WorldVector<double> >& Lambda,
@@ -628,20 +628,14 @@ namespace AMDiS {
Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -649,30 +643,21 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
protected:
- /** \brief
- * DOFVector to be evaluated at quadrature points.
- */
+ /// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
- /** \brief
- * Pointer to an array containing the DOFVector evaluated at quadrature
- * points.
- */
+ /// Pointer to an array containing the DOFVector evaluated at quadrature points.
double* vecAtQPs1;
double* vecAtQPs2;
- /** \brief
- * Function evaluated at \ref vecAtQPs.
- */
+ /// Function evaluated at \ref vecAtQPs.
BinaryAbstractFunction<double, double, double> *f;
};
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 6681e8fd..952458d7 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -770,7 +770,7 @@ namespace AMDiS {
matrix->getBoundaryManager()->exitMatrix(matrix);
if (matrix)
- nnz += matrix->getBaseMatrix().nnz();
+ nnz += matrix->getBaseMatrix().nnz();
}
// And now assemble boundary conditions on the vectors
diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc
index f9645b11..b33c77d0 100644
--- a/AMDiS/src/Quadrature.cc
+++ b/AMDiS/src/Quadrature.cc
@@ -39,7 +39,7 @@ namespace AMDiS {
{
FUNCNAME("Quadrature::grdFAtQp()");
- static WorldVector<double> *quad_vec_d = NULL;
+ static WorldVector<double> *quad_vec_d = NULL;
static size_t size = 0;
WorldVector<double> *val;
WorldVector<double> grd;
@@ -47,7 +47,11 @@ namespace AMDiS {
if (vec) {
val = vec;
} else {
- if (static_cast<int>( size) < n_points) {
+#ifdef _OPENMP
+ ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
+#endif
+
+ if (static_cast<int>(size) < n_points) {
size_t new_size = std::max(maxNQuadPoints[dim], n_points);
delete [] quad_vec_d;
quad_vec_d = new WorldVector<double>[new_size];
@@ -67,9 +71,8 @@ namespace AMDiS {
return val;
}
- const double
- *Quadrature::fAtQp(const AbstractFunction<double, DimVec<double> >& f,
- double *vec) const
+ const double *Quadrature::fAtQp(const AbstractFunction<double, DimVec<double> >& f,
+ double *vec) const
{
FUNCNAME("Quadrature::fAtQp()");
@@ -80,6 +83,10 @@ namespace AMDiS {
if (vec) {
val = vec;
} else {
+#ifdef _OPENMP
+ ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
+#endif
+
if (static_cast<int>( size) < n_points) {
size_t new_size = std::max(maxNQuadPoints[dim], n_points);
delete [] quad_vec;
@@ -1532,9 +1539,8 @@ namespace AMDiS {
// fill memory
for (int i = 0; i< nPoints; i++) {
lambda = quadrature->getLambda(i);
- for (int j = 0; j < nBasFcts; j++) {
+ for (int j = 0; j < nBasFcts; j++)
(*(basisFunctions->getGrdPhi(j)))(lambda, (*(grdPhi))[i][j]);
- }
}
// update flag
@@ -1551,9 +1557,8 @@ namespace AMDiS {
// fill memory
for (int i = 0; i < nPoints; i++) {
lambda = quadrature->getLambda(i);
- for (int j = 0; j < nBasFcts; j++) {
+ for (int j = 0; j < nBasFcts; j++)
(*(basisFunctions->getD2Phi(j)))(lambda, (*(D2Phi))[i][j]);
- }
}
// update flag
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 61e1d648..7ac17dc9 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -186,11 +186,16 @@ namespace AMDiS {
tmpLALt[myRank] = new DimMat<double>*[nPoints];
for (int j = 0; j < nPoints; j++)
tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);
-
- const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
- basFcts = owner->getColFESpace()->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
+
+#ifdef _OPENMP
+#pragma omp critical (dofIndexAccess)
+#endif
+ {
+ psiFast = updateFastQuadrature(psiFast, owner->getRowFESpace()->getBasisFcts(),
+ INIT_GRD_PHI);
+ phiFast = updateFastQuadrature(phiFast, owner->getRowFESpace()->getBasisFcts(),
+ INIT_GRD_PHI);
+ }
firstCall = false;
}
@@ -211,7 +216,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
-
+
grdPsi = psiFast->getGradient(iq);
grdPhi = phiFast->getGradient(iq);
@@ -267,9 +272,8 @@ namespace AMDiS {
int myRank = omp_get_thread_num();
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
- (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
- }
+ for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
+ (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
if (symmetric) {
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index e616e470..7fe786d9 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -156,8 +156,8 @@ namespace AMDiS {
const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();
TEST_EXIT_DBG(vec)("no dof vector!\n");
-
- TEST_EXIT(elInfo->getMesh() == vec->getFESpace()->getMesh())("Vector and fe space do not fit together!\n");
+ TEST_EXIT_DBG(elInfo->getMesh() == vec->getFESpace()->getMesh())
+ ("Vector and fe space do not fit together!\n");
Quadrature *localQuad = quad ? quad : quadrature;
@@ -170,7 +170,6 @@ namespace AMDiS {
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
double *values = valuesAtQPs[vec]->values.getValArray();
-
bool sameFESpaces =
(vec->getFESpace() == owner->rowFESpace) ||
(vec->getFESpace() == owner->colFESpace);
@@ -191,7 +190,7 @@ namespace AMDiS {
if (opt && !quad && sameFESpaces) {
if (psiFast->getBasisFunctions() == basFcts) {
vec->getVecAtQPs(elInfo, NULL, psiFast, values);
- } else if(phiFast->getBasisFunctions() == basFcts) {
+ } else if (phiFast->getBasisFunctions() == basFcts) {
vec->getVecAtQPs(elInfo, NULL, phiFast, values);
} else {
vec->getVecAtQPs(elInfo, localQuad, NULL, values);
--
GitLab