diff --git a/AMDiS/src/BasisFunction.cc b/AMDiS/src/BasisFunction.cc
index 06461d8fb2b8b94fb3d046b0875a229dab04479f..878615f391f9ffe0862cbe6746f541abd9c72cee 100644
--- a/AMDiS/src/BasisFunction.cc
+++ b/AMDiS/src/BasisFunction.cc
@@ -4,6 +4,7 @@
#include "DOFVector.h"
#include "BasisFunction.h"
#include "Lagrange.h"
+#include "OpenMP.h"
namespace AMDiS {
@@ -14,19 +15,34 @@ namespace AMDiS {
/* are those corresponding to these barycentric coordinates. */
/****************************************************************************/
- BasisFunction::~BasisFunction()
- {
- DELETE nDOF;
- }
-
-
BasisFunction::BasisFunction(const ::std::string& name_, int dim_, int degree_)
- : name(name_), degree(degree_), dim(dim_)
+ : name(name_),
+ degree(degree_),
+ dim(dim_)
{
FUNCNAME("BasisFunction::BasisFunction()");
+
nDOF = NEW DimVec<int>(dim, DEFAULT_VALUE, -1);
+ dow = Global::getGeo(WORLD);
+
+ grdTmpVec1.resize(omp_get_max_threads());
+ grdTmpVec2.resize(omp_get_max_threads());
+
+ for (int i = 0; i < omp_get_max_threads(); i++) {
+ grdTmpVec1[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
+ grdTmpVec2[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
+ }
};
+ BasisFunction::~BasisFunction()
+ {
+ DELETE nDOF;
+
+ for (int i = 0; i < grdTmpVec1.size(); i++) {
+ DELETE grdTmpVec1[i];
+ DELETE grdTmpVec2[i];
+ }
+ }
/****************************************************************************/
/* some routines for evaluation of a finite element function, its gradient */
@@ -50,9 +66,8 @@ namespace AMDiS {
const WorldVector<double> *uh_loc,
WorldVector<double>* values) const
{
- static WorldVector<double> Values(DEFAULT_VALUE, 0.);
+ static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
WorldVector<double> *val = (NULL != values) ? values : &Values;
- int dow = Global::getGeo(WORLD);
for (int n = 0; n < dow; n++)
(*val)[n] = 0;
@@ -70,32 +85,27 @@ namespace AMDiS {
const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& grd_lambda,
const double *uh_loc,
- WorldVector<double>* grd_uh) const
+ WorldVector<double>* val) const
{
- static WorldVector<double> grd;
+ TEST_EXIT_DBG(val)("return value is NULL\n");
+
+ DimVec<double> *grdTmp1 = grdTmpVec1[omp_get_thread_num()];
+ DimVec<double> *grdTmp2 = grdTmpVec2[omp_get_thread_num()];
- DimVec<double> grd_b(dim, DEFAULT_VALUE, 0.0);
- WorldVector<double> *val;
- DimVec<double> grd1(dim, DEFAULT_VALUE, 0.);
-
- val = grd_uh ? grd_uh : &grd;
-
for (int j = 0; j < dim + 1; j++)
- grd1[j] = 0.0;
+ (*grdTmp2)[j] = 0.0;
for (int i = 0; i < nBasFcts; i++) {
- (*(*grdPhi)[i])(lambda, grd_b);
+ (*(*grdPhi)[i])(lambda, *grdTmp1);
for (int j = 0; j < dim + 1; j++)
- grd1[j] += uh_loc[i] * grd_b[j];
+ (*grdTmp2)[j] += uh_loc[i] * (*grdTmp1)[j];
}
- int dow = Global::getGeo(WORLD);
-
for (int i = 0; i < dow; i++) {
- (*val)[i] = 0;
+ (*val)[i] = 0.0;
for (int j = 0; j < dim + 1; j++)
- (*val)[i] += grd_lambda[j][i] * grd1[j];
+ (*val)[i] += grd_lambda[j][i] * (*grdTmp2)[j];
}
return ((*val));
@@ -105,7 +115,7 @@ namespace AMDiS {
const DimVec<WorldVector<double> >& grd_lambda,
const double *uh_loc, WorldMatrix<double>* D2_uh) const
{
- static WorldMatrix<double> D2(DEFAULT_VALUE, 0.);
+ static WorldMatrix<double> D2(DEFAULT_VALUE, 0.0);
DimMat<double> D2_b(dim, DEFAULT_VALUE, 0.0);
DimMat<double> D2_tmp(dim, DEFAULT_VALUE, 0.0);
WorldMatrix<double> *val = D2_uh ? D2_uh : &D2;
@@ -117,8 +127,6 @@ namespace AMDiS {
D2_tmp[k][l] += uh_loc[i] * D2_b[k][l];
}
- int dow = Global::getGeo(WORLD);
-
for (int i = 0; i < dow; i++)
for (int j = 0; j < dow; j++) {
(*val)[i][j] = 0.0;
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index e5c2485ac84d8a0597d154613319045ebbb2aeed..87fea8fafc346e916cab917de5554ba8d4f8426d 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -432,6 +432,11 @@ namespace AMDiS {
*/
int dim;
+ /** \brief
+ * Dimension of the world.
+ */
+ int dow;
+
/** \brief
* Number of DOFs at the different positions
*/
@@ -450,7 +455,20 @@ namespace AMDiS {
/** \brief
* Vector of second derivatives
*/
- ::std::vector<D2BasFctType*> *d2Phi;
+ ::std::vector<D2BasFctType*> *d2Phi;
+
+
+ /** \brief
+ * Is used by function evalGrdUh. To make it thread safe, we need a
+ * temporary DimVec for every thread.
+ */
+ ::std::vector<DimVec<double>* > grdTmpVec1;
+
+ /** \brief
+ * Is used by function evalGrdUh. To make it thread safe, we need a
+ * temporary DimVec for every thread.
+ */
+ ::std::vector<DimVec<double>* > grdTmpVec2;
};
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index d9ba341974217ec652e1b4b3c04208137eeeba17..aacea37647cb9d61f56fd1f46ff03596edc8b0a7 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -53,7 +53,7 @@ namespace AMDiS {
// define result vector
static DOFVector<WorldVector<double> > *result = NULL;
- if(grad) {
+ if (grad) {
result = grad;
} else {
if(result && result->getFESpace() != feSpace) {
@@ -62,8 +62,6 @@ namespace AMDiS {
}
}
- int i, j;
-
Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim();
@@ -80,11 +78,11 @@ namespace AMDiS {
int numNodes = 0;
int numDOFs = 0;
- for(i = 0; i < dim + 1; i++) {
+ for (int i = 0; i < dim + 1; i++) {
GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
int numPositionNodes = mesh->getGeo(geoIndex);
int numPreDOFs = admin->getNumberOfPreDOFs(i);
- for(j = 0; j < numPositionNodes; j++) {
+ for (int j = 0; j < numPositionNodes; j++) {
int dofs = basFcts->getNumberOfDOFs(geoIndex);
numNodeDOFs.push_back(dofs);
numDOFs += dofs;
@@ -99,7 +97,7 @@ namespace AMDiS {
TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
("number of dofs != number of basis functions\n");
- for(i = 0; i < numDOFs; i++) {
+ for (int i = 0; i < numDOFs; i++) {
bary.push_back(basFcts->getCoords(i));
}
@@ -112,22 +110,19 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
- while(elInfo) {
+ while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
int localDOFNr = 0;
- for(i = 0; i < numNodes; i++) { // for all nodes
- for(j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
+ for (int i = 0; i < numNodes; i++) { // for all nodes
+ for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
- if(!visited[dofIndex]) {
-
- result[dofIndex] = basFcts->evalGrdUh(*(bary[localDOFNr]),
- grdLambda,
- localUh,
- NULL);
+ if (!visited[dofIndex]) {
+ basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda,
+ localUh, &((*result)[dofIndex]));
visited[dofIndex] = true;
}
@@ -168,8 +163,6 @@ namespace AMDiS {
DOFVector<double> volume(feSpace, "volume");
volume.set(0.0);
- int i;
-
Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim();
@@ -190,17 +183,15 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
- while(elInfo) {
+ 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();
- const WorldVector<double> &grd = basFcts->evalGrdUh(bary,
- grdLambda,
- localUh,
- NULL);
+ WorldVector<double> grd;
+ basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
- for(i = 0; i < dim + 1; i++) {
+ for (int i = 0; i < dim + 1; i++) {
DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
(*result)[dofIndex] += grd * det;
volume[dofIndex] += det;
@@ -212,8 +203,8 @@ namespace AMDiS {
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);
}
}
@@ -248,7 +239,6 @@ namespace AMDiS {
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
- int nBasFcts = basFcts->getNumber();
static WorldVector<double> *grd = NULL;
@@ -349,7 +339,6 @@ namespace AMDiS {
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
- int nBasFcts = basFcts->getNumber();
int i, j, k, l, iq;
static WorldMatrix<double> *vec = NULL;
@@ -403,7 +392,6 @@ namespace AMDiS {
for (i = 0; i < nBasFcts; i++) {
WARNING("not tested after index correction\n");
- //(*(basFcts->getD2Phi(j)))(quad->getLambda(i), D2Phi);
(*(basFcts->getD2Phi(i)))(quad->getLambda(iq), D2Phi);
for (k = 0; k < parts; k++)
@@ -441,7 +429,7 @@ namespace AMDiS {
this->set(0.0);
- DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, nBasisFcts);
+ DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
double *sourceLocalCoeffs = GET_MEMORY(double, nSourceBasisFcts);
if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
@@ -454,13 +442,13 @@ namespace AMDiS {
while (elInfo) {
Element *el = elInfo->getElement();
- basisFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
+ basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
source->getLocalVector(el, sourceLocalCoeffs);
for (int i = 0; i < nBasisFcts; i++) {
- if (vec[localIndices[i]] == 0.0) {
+ if (vec[myLocalIndices[i]] == 0.0) {
coords = basisFcts->getCoords(i);
- vec[localIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
+ vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
}
}
elInfo = stack.traverseNext(elInfo);
@@ -483,12 +471,12 @@ namespace AMDiS {
while (nextTraverse) {
basisFcts->getLocalIndices(elInfo1->getElement(),
feSpace->getAdmin(),
- localIndices);
+ myLocalIndices);
source->getLocalVector(elInfo2->getElement(),
sourceLocalCoeffs);
for (int i = 0; i < nBasisFcts; i++) {
- if (vec[localIndices[i]] == 0.0) {
+ if (vec[myLocalIndices[i]] == 0.0) {
coords1 = basisFcts->getCoords(i);
elInfo1->coordToWorld(*coords1, &worldVec);
elInfo2->worldToCoord(worldVec, &coords2);
@@ -502,7 +490,7 @@ namespace AMDiS {
}
if (isPositive) {
- vec[localIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);
+ vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);
}
}
}
@@ -512,7 +500,6 @@ namespace AMDiS {
}
}
- FREE_MEMORY(localIndices, DegreeOfFreedom, nBasisFcts);
FREE_MEMORY(sourceLocalCoeffs, double, nSourceBasisFcts);
}
@@ -539,7 +526,7 @@ namespace AMDiS {
int vNumBasFcts = vBasFcts->getNumber();
if (feSpace->getMesh() == vFESpace->getMesh()) {
- DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, numBasFcts);
+ DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
WorldVector<double> *vLocalCoeffs = NEW WorldVector<double>[vNumBasFcts];
Mesh *mesh = feSpace->getMesh();
TraverseStack stack;
@@ -550,20 +537,19 @@ namespace AMDiS {
while (elInfo) {
Element *el = elInfo->getElement();
- basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
+ basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
v->getLocalVector(el, vLocalCoeffs);
for (int i = 0; i < numBasFcts; i++) {
- if (vec[localIndices[i]] == nul) {
+ if (vec[myLocalIndices[i]] == nul) {
coords = basFcts->getCoords(i);
- vec[localIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs,NULL) * factor;
+ vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs,NULL) * factor;
}
}
elInfo = stack.traverseNext(elInfo);
}
- FREE_MEMORY(localIndices, DegreeOfFreedom, numBasFcts);
DELETE [] vLocalCoeffs;
} else {
ERROR_EXIT("not yet for dual traverse\n");
@@ -651,10 +637,8 @@ namespace AMDiS {
DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
if (!visited[dofIndex]) {
- grd = basFcts->evalGrdUh(*(bary[localDOFNr]),
- grdLambda,
- localUh,
- NULL);
+ basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda,
+ localUh, &grd);
for (int k = 0; k < dow; k++) {
(*result)[k][dofIndex] = grd[k];
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index d13173b1b1f01edc373559440439532ccd62da41..5ae335c960dc61a83eaca0bce0c8827562cb8525 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -38,6 +38,7 @@
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h"
+#include "BasisFunction.h"
#include <vector>
#include <memory>
#include <list>
@@ -69,18 +70,17 @@ namespace AMDiS {
class DOFVectorBase : public DOFIndexed<T>
{
public:
+
DOFVectorBase()
: feSpace(NULL),
elementVector(NULL),
- boundaryManager(NULL)
+ boundaryManager(NULL),
+ nBasFcts(0)
{};
+
+ DOFVectorBase(const FiniteElemSpace *f, ::std::string n);
- DOFVectorBase(const FiniteElemSpace *f, ::std::string n)
- : feSpace(f),
- name(n),
- elementVector(NULL),
- boundaryManager(NULL)
- {};
+ virtual ~DOFVectorBase();
virtual const T *getLocalVector(const Element *el, T* localVec) const;
@@ -156,9 +156,14 @@ namespace AMDiS {
*/
T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);
- inline ::std::vector<Operator*>& getOperators() { return operators; };
+ inline ::std::vector<Operator*>& getOperators() {
+ return operators;
+ };
+
+ inline ::std::vector<double*>& getOperatorFactor() {
+ return operatorFactor;
+ };
- inline ::std::vector<double*>& getOperatorFactor() { return operatorFactor; };
inline ::std::vector<double*>& getOperatorEstFactor() {
return operatorEstFactor;
};
@@ -166,45 +171,66 @@ namespace AMDiS {
/** \brief
* Returns \ref name
*/
- inline const ::std::string& getName() const { return name; };
+ inline const ::std::string& getName() const {
+ return name;
+ };
- inline BoundaryManager* getBoundaryManager() const { return boundaryManager; };
+ inline BoundaryManager* getBoundaryManager() const {
+ return boundaryManager;
+ };
inline void setBoundaryManager(BoundaryManager *bm) {
boundaryManager = bm;
};
protected:
+ /** \brief
+ *
+ */
const FiniteElemSpace *feSpace;
+ /** \brief
+ *
+ */
::std::string name;
+ /** \brief
+ *
+ */
ElementVector *elementVector;
+ /** \brief
+ *
+ */
::std::vector<Operator*> operators;
+ /** \brief
+ *
+ */
::std::vector<double*> operatorFactor;
+ /** \brief
+ *
+ */
::std::vector<double*> operatorEstFactor;
+ /** \brief
+ *
+ */
BoundaryManager *boundaryManager;
- };
+ /** \brief
+ * Number of basis functions of the used finite element space.
+ */
+ int nBasFcts;
- // template<typename T>
- // class QPEvaluatable : public DOFVectorBase<T>
- // {
- // public:
- // QPEvaluatable()
- // : DOFVectorBase<T>()
- // {};
-
- // QPEvaluatable(const FiniteElemSpace *f, ::std::string n)
- // : DOFVectorBase<T>(f, n)
- // {};
+ /** \brief
+ * Are used to store temporary local dofs of an element.
+ * Thread safe.
+ */
+ ::std::vector<DegreeOfFreedom* > localIndices;
+ };
- // protected:
- // };
// ===========================================================================
// ===== defs ================================================================
@@ -280,7 +306,6 @@ namespace AMDiS {
*/
DOFVector()
: DOFVectorBase<T>(),
- //elementVector(NULL),
refineInter(false),
feSpace(NULL),
coarsenOperation(NO_OPERATION)
@@ -289,7 +314,7 @@ namespace AMDiS {
/** \brief
* Constructs a DOFVector with name n belonging to FiniteElemSpace f
*/
- DOFVector(const FiniteElemSpace* f,::std::string n);
+ DOFVector(const FiniteElemSpace* f, ::std::string n);
/** \brief
* Initialization.
@@ -300,9 +325,9 @@ namespace AMDiS {
* Copy Constructor
*/
DOFVector(const DOFVector& rhs) {
- *this=rhs;
- name=rhs.name+"copy";
- if(feSpace && feSpace->getAdmin()) {
+ *this = rhs;
+ name = rhs.name + "copy";
+ if (feSpace && feSpace->getAdmin()) {
(dynamic_cast<DOFAdmin*>(feSpace->getAdmin()))->addDOFIndexed(this);
}
};
@@ -374,16 +399,10 @@ namespace AMDiS {
/** \brief
* Returns \ref vec
*/
- // ::std::vector<T>& getVector() const {return (::std::vector<T>&) vec;} ;
::std::vector<T>& getVector() {
- return vec;
+ return vec;
};
- // /** \brief
- // * Returns \ref feSpace
- // */
- // inline const FiniteElemSpace* getFESpace() const { return feSpace; };
-
/** \brief
* Returns size of \ref vec
*/
@@ -439,9 +458,7 @@ namespace AMDiS {
("Illegal vector index %d.\n",i);
return vec[i];
};
-
- // Andreas Ergaenzungen:
-
+
/** \brief
* Calculates Integral of this DOFVector
*/
@@ -451,10 +468,7 @@ namespace AMDiS {
* Calculates L1 norm of this DOFVector
*/
double L1Norm(Quadrature* q = NULL) const;
-
- // Ende Andreas Ergaenzungen (hier) unten geht es weiter
-
-
+
/** \brief
* Calculates L2 norm of this DOFVector
*/
@@ -507,16 +521,12 @@ namespace AMDiS {
inline double l1norm() const {
return asum();
};
-
- // hier Andreas Ergaenzung
-
+
/** \brief
* Calculates the sum of this DOFVector
*/
T sum() const;
-
- // bis hier Andreas E.
-
+
/** \brief
* Sets \ref vec[i] = val, i=0 , ... , size
*/
@@ -535,58 +545,11 @@ namespace AMDiS {
*/
DOFVector<T>& operator=(const DOFVector<T>& );
-
- //virtual DOFVector<T>& operator-=(const DOFVector<T>& rhs) {
- // axpy(-1.0, rhs); return *this;
- //};
-
- //virtual DOFVector<T>& operator+=(const DOFVector<T>& rhs) {
- // axpy(1.0, rhs); return *this;
- //};
-
- /** \brief
- * Multiplies all entries of \ref vec with s
- */
- //virtual void scal(T s);
-
- /** \brief
- * Multiplication with a scalar
- */
- //virtual const DOFVector<T>& operator*(T d) const {
- // static DOFVector<T> vec(this->feSpace, "operator*-vector");
- // vec = *this;
- // vec.scal(d);
- // return vec;
- //};
-
- /** \brief
- * \ref vec[i] = d * \ref vec[i]
- */
- //virtual const DOFVector<T>& operator*=(T d) {
- // scal(d);
- // return (*this);
- //};
-
/** \brief
* vec[i] = v.vec[i]
*/
void copy(const DOFVector<T>& v);
- /** \brief
- * vec[i] = a * x.vec[i] + vec[i]
- */
- //void axpy(T a,const DOFVector<T>& x);
-
- /** \brief
- * vec[i] = a * x.vec[i] + y.vec[i]
- */
- //void axpy(T a,const DOFVector<T>& x, const DOFVector<T>& y);
-
- /** \brief
- * vec[i] = x.vec[i] + a * vec[i]
- */
- //void xpay(T a,const DOFVector<T>& x);
-
/** \brief
* Returns minimum of DOFVector
*/
@@ -638,10 +601,7 @@ namespace AMDiS {
getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
- protected:
-
- //Hier auch Ergaenzungen (Andreas)
-
+ protected:
/** \brief
* Used by Int while mesh traversal
*/
@@ -651,11 +611,7 @@ namespace AMDiS {
* Used by L1Norm while mesh traversal
*/
static int L1Norm_fct(ElInfo* elinfo);
-
- // Letzte Erganzung in diesem File (Andreas)
-
-
-
+
/** \brief
* Used by L2Norm while mesh traversal
*/
@@ -670,7 +626,7 @@ namespace AMDiS {
/** \brief
* Name of this DOFVector
*/
- ::std::string name;
+ ::std::string name;
/** \brief
* FiniteElemSpace of the vector
@@ -682,8 +638,6 @@ namespace AMDiS {
*/
::std::vector<T> vec;
- // ElementVector *elementVector;
-
/** \brief
* Specifies whether interpolation should be performed after refinement
*/
@@ -694,12 +648,6 @@ namespace AMDiS {
*/
CoarsenOperation coarsenOperation;
- // ::std::vector<Operator*> operators;
-
- // ::std::vector<double*> operatorFactor;
-
- // ::std::vector<double*> operatorEstFactor;
-
/** \brief
* Used by \ref interpol
*/
@@ -711,11 +659,6 @@ namespace AMDiS {
static DOFVector<T> *traverseVector;
protected:
- /** \brief
- * Used for mesh traversal
- */
- // static DOFVector<T> *traverseVector;
-
/** \brief
* Used while calculating vector norms
*/
@@ -732,14 +675,21 @@ namespace AMDiS {
static int dim;
};
+ // ===============================================================================
+
template<>
void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
template<>
void DOFVector<double>::coarseRestrict(RCNeighbourList&, int);
- inline double min(const DOFVector<double>& v) {return v.min();};
- inline double max(const DOFVector<double>& v) {return v.max();};
+ inline double min(const DOFVector<double>& v) {
+ return v.min();
+ };
+
+ inline double max(const DOFVector<double>& v) {
+ return v.max();
+ };
// ===========================================================================
// ===== class DOFVectorDOF ==================================================
@@ -800,6 +750,11 @@ namespace AMDiS {
DOFVectorDOF();
};
+
+
+ // ===============================================================================
+
+
template<typename T>
double norm(DOFVector<T> *vec) {
return vec->nrm2();
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 9e1b315c412b2c2216e2e19ca6df35f04c331334..fb1e841098bf61ee9de10640e0d86b4aaaecf30a 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -18,29 +18,31 @@
namespace AMDiS {
template<typename T>
- void DOFVectorBase<T>::addElementVector(T factor,
- const ElementVector &elVec,
- const BoundaryType *bound,
- bool add)
+ DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, ::std::string n)
+ : feSpace(f),
+ name(n),
+ elementVector(NULL),
+ boundaryManager(NULL)
{
- FUNCNAME("DOFVector::addElementVector()");
+ nBasFcts = feSpace->getBasisFcts()->getNumber();
- int n_row = elVec.getSize();
+ localIndices.resize(omp_get_max_threads());
- for (DegreeOfFreedom i = 0; i < n_row; i++) {
- BoundaryCondition *condition =
- bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;
-
- if(!(condition && condition->isDirichlet())) {
- DegreeOfFreedom irow = elVec.dofIndices[i];
- (*this)[irow] = (add ? (*this)[irow] : 0.0);
- (*this)[irow] += factor * elVec[i];
- }
+ for (int i = 0; i < omp_get_max_threads(); i++) {
+ localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
}
}
-
+
template<typename T>
- DOFVector<T>::DOFVector(const FiniteElemSpace* f,::std::string n)
+ DOFVectorBase<T>::~DOFVectorBase()
+ {
+ for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
+ FREE_MEMORY(localIndices[i], DegreeOfFreedom, this->nBasFcts);
+ }
+ }
+
+ template<typename T>
+ DOFVector<T>::DOFVector(const FiniteElemSpace* f, ::std::string n)
: DOFVectorBase<T>(f, n),
refineInter(false),
coarsenOperation(NO_OPERATION)
@@ -49,13 +51,14 @@ namespace AMDiS {
}
template<typename T>
- void DOFVector<T>::init(const FiniteElemSpace* f,::std::string n)
+ void DOFVector<T>::init(const FiniteElemSpace* f, ::std::string n)
{
name = n;
feSpace = f;
if (feSpace && feSpace->getAdmin()) {
(feSpace->getAdmin())->addDOFIndexed(this);
- }
+ }
+
this->boundaryManager = NEW BoundaryManager;
}
@@ -83,6 +86,28 @@ namespace AMDiS {
template<typename T>
int DOFVector<T>::dim = 0;
+ template<typename T>
+ void DOFVectorBase<T>::addElementVector(T factor,
+ const ElementVector &elVec,
+ const BoundaryType *bound,
+ bool add)
+ {
+ FUNCNAME("DOFVector::addElementVector()");
+
+ int n_row = elVec.getSize();
+
+ for (DegreeOfFreedom i = 0; i < n_row; i++) {
+ BoundaryCondition *condition =
+ bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;
+
+ if(!(condition && condition->isDirichlet())) {
+ DegreeOfFreedom irow = elVec.dofIndices[i];
+ (*this)[irow] = (add ? (*this)[irow] : 0.0);
+ (*this)[irow] += factor * elVec[i];
+ }
+ }
+ }
+
template<typename T>
double DOFVector<T>::nrm2() const
{
@@ -1129,11 +1154,7 @@ namespace AMDiS {
const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
{
static T* localVec = NULL;
- const DOFAdmin* admin = feSpace->getAdmin();
- int nBasFcts = feSpace->getBasisFcts()->getNumber();
-
- T *result;
-
+ T *result;
if (d) {
result = d;
} else {
@@ -1143,16 +1164,14 @@ namespace AMDiS {
result = localVec;
}
- DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
-
- feSpace->getBasisFcts()->getLocalIndices(el, admin, localIndices);
+ DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
+ feSpace->getBasisFcts()->getLocalIndices(el, feSpace->getAdmin(),
+ myLocalIndices);
for (int i = 0; i < nBasFcts; i++) {
- result[i] = (*this)[localIndices[i]];
+ result[i] = (*this)[myLocalIndices[i]];
}
-
- delete [] localIndices;
-
+
return result;
}
@@ -1180,7 +1199,6 @@ namespace AMDiS {
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
- int nBasFcts = basFcts->getNumber();
static T *localvec = NULL;
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 2e284a2ef55555a0f83656cf14fd63c5cb23fe63..ad8f8c94c32d51deeb0185366301bd1692476a93 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -404,11 +404,6 @@ namespace AMDiS {
}
}
- // if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- // det = calcGrdLambda(*Lambda);
- // }
-
-
if (fillFlag__local.isSet(Mesh::FILL_NEIGH) ||
fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
diff --git a/AMDiS/src/FileWriter.cc b/AMDiS/src/FileWriter.cc
index 2831fa89770fa4e56784bc0febcd2e1dcec01482..1c9b9147ce770e41c3b6677efbd1088475f2f112 100644
--- a/AMDiS/src/FileWriter.cc
+++ b/AMDiS/src/FileWriter.cc
@@ -207,6 +207,9 @@ namespace AMDiS {
writingIsDelayed_ = true;
delayedFilename_ = fn;
+
+ MSG("Delayed writing of ParaView file %s\n", (fn + paraViewFileExt).c_str());
+
return;
}
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index f91407c8a9ffe4ae3c1fd6eb6926406a79173525..a1a8d835ef11c627522f207624086d361e202e9e 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -891,7 +891,6 @@ namespace AMDiS {
el_info->testFlag(Mesh::FILL_COORDS);
- int dow = Global::getGeo(WORLD);
int vertices = Global::getGeo(VERTEX, dim);
if (b_no) {
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 26b15aab6307607bac36286f558ff1a1cb9f582d..a96f3e2f1e060fab532a1dc15b565ef9a58b5661 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -616,8 +616,6 @@ namespace AMDiS {
INFO(info_, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first,clock()));
- // ::std::cout << "Speicherverbrauch: " << memSizeStr(systemMatrix_->memsize()) << ::std::endl;
-
return;
}
diff --git a/AMDiS/src/Recovery.cc b/AMDiS/src/Recovery.cc
index e65cc77f1257c5078618c34d8c8ea3dcbd1940aa..5b2a5428249506bcbe31d4cc811379054c4368b5 100644
--- a/AMDiS/src/Recovery.cc
+++ b/AMDiS/src/Recovery.cc
@@ -1021,43 +1021,37 @@ Recovery::recovery(DOFVector<double> *uh,
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
- while (elInfo)
- {
- double det = elInfo->getDet();
- const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
- const double *localUh = uh->getLocalVector(elInfo->getElement(), NULL);
- const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- const WorldVector<double> &grd = basFcts->evalGrdUh(bary,
- grdLambda,
- localUh,
- NULL);
-
- double fAtBary = 1.0;
-
- if (f_vec)
- {
- elInfo->coordToWorld(bary, &barycenter);
- fAtBary = (*f_vec)(barycenter);
- }
-
- if (f_scal)
- {
- if (aux_vec)
- localUh = aux_vec->getLocalVector(elInfo->getElement(), NULL);
-
- fAtBary = basFcts->evalUh(bary, localUh);
- fAtBary = (*f_scal)(fAtBary);
- }
-
- for (i=0; i<dim+1; i++)
- {
- DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
- (*result)[dofIndex] += grd * fAtBary * det;
- volume[dofIndex] += det;
- }
-
- elInfo = stack.traverseNext(elInfo);
+ while (elInfo) {
+ double det = elInfo->getDet();
+ const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
+ const double *localUh = uh->getLocalVector(elInfo->getElement(), NULL);
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ WorldVector<double> grd;
+ basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
+
+ double fAtBary = 1.0;
+
+ if (f_vec) {
+ elInfo->coordToWorld(bary, &barycenter);
+ fAtBary = (*f_vec)(barycenter);
+ }
+
+ if (f_scal) {
+ if (aux_vec)
+ localUh = aux_vec->getLocalVector(elInfo->getElement(), NULL);
+
+ fAtBary = basFcts->evalUh(bary, localUh);
+ fAtBary = (*f_scal)(fAtBary);
}
+
+ for (i = 0; i < dim + 1; i++) {
+ DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
+ (*result)[dofIndex] += grd * fAtBary * det;
+ volume[dofIndex] += det;
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ }
DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
diff --git a/AMDiS/src/RobinBC.hh b/AMDiS/src/RobinBC.hh
index 883a1331bfff3ff912eacec69b6dbb8b95a0ad46..f9a49ba86932be557a2edaa78a3b0f3a4ba95177 100644
--- a/AMDiS/src/RobinBC.hh
+++ b/AMDiS/src/RobinBC.hh
@@ -94,7 +94,8 @@ double RobinBC<T>::boundResidual(ElInfo *elInfo, Estimator<T> *estimator)
(*neumannOperators)[face]->getAssembler()->getZeroOrderAssembler()->
initElement(elInfo);
- for (val = iq = 0; iq < numPoints; iq++) {
+ val = 0.0;
+ for (iq = 0; iq < numPoints; iq++) {
lambda = quadrature->getLambda(iq);
estimator->bas_fcts->evalGrdUh(lambda, Lambda, estimator->uhEl, &grdUh);