diff --git a/AMDiS/src/AbstractFunction.h b/AMDiS/src/AbstractFunction.h
index bb5925fd5f0bb660baa56896bcdd5c768f9a80e0..1ce4a63b3cd0493d35a7499380891c5e5e4db383 100644
--- a/AMDiS/src/AbstractFunction.h
+++ b/AMDiS/src/AbstractFunction.h
@@ -49,19 +49,23 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
- AbstractFunction(int degree = 0) : degree_(degree) {};
+ AbstractFunction(int degree = 0) :
+ degree_(degree)
+ {};
virtual ~AbstractFunction() {};
/** \brief
* Returns \ref degree_.
*/
- inline int getDegree() const { return degree_; };
+ inline int getDegree() const {
+ return degree_;
+ };
/** \brief
* Deligates the evaluation to overriden method f.
*/
- virtual const ReturnType& operator()(const ArgumentType& x) const = 0;
+ virtual ReturnType operator()(const ArgumentType& x) const = 0;
protected:
int degree_;
@@ -86,7 +90,9 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
- BinaryAbstractFunction(int degree = 0) : degree_(degree) {};
+ BinaryAbstractFunction(int degree = 0) :
+ degree_(degree)
+ {};
virtual ~BinaryAbstractFunction() {};
@@ -100,8 +106,8 @@ namespace AMDiS {
/** \brief
* Deligates the evaluation to overriden method f.
*/
- virtual const ReturnType& operator()(const ArgumentType1& x,
- const ArgumentType2& y) const = 0;
+ virtual ReturnType operator()(const ArgumentType1& x,
+ const ArgumentType2& y) const = 0;
protected:
int degree_;
@@ -127,7 +133,9 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
- TertiaryAbstractFunction(int degree = 0) : degree_(degree) {};
+ TertiaryAbstractFunction(int degree = 0) :
+ degree_(degree)
+ {};
virtual ~TertiaryAbstractFunction() {};
@@ -141,9 +149,9 @@ namespace AMDiS {
/** \brief
* function evaluation.
*/
- virtual const ReturnType& operator()(const ArgumentType1& x,
- const ArgumentType2& y,
- const ArgumentType3& z) const = 0;
+ virtual ReturnType operator()(const ArgumentType1& x,
+ const ArgumentType2& y,
+ const ArgumentType3& z) const = 0;
protected:
int degree_;
@@ -171,7 +179,9 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
- QuartAbstractFunction(int degree = 0) : degree_(degree) {};
+ QuartAbstractFunction(int degree = 0) :
+ degree_(degree)
+ {};
virtual ~QuartAbstractFunction() {};
@@ -185,10 +195,10 @@ namespace AMDiS {
/** \brief
* function evaluation.
*/
- virtual const ReturnType& operator()(const ArgumentType1& x,
- const ArgumentType2& y,
- const ArgumentType3& z,
- const ArgumentType4& u) const = 0;
+ virtual ReturnType operator()(const ArgumentType1& x,
+ const ArgumentType2& y,
+ const ArgumentType3& z,
+ const ArgumentType4& u) const = 0;
protected:
int degree_;
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index 2ea2f69a7c31bea27d2680240d84bbfa3424eaf3..6b9170d0e48c5fa3d6e0f23ba57888a4bf6c168d 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -37,6 +37,7 @@
#include "FirstOrderAssembler.h"
#include "SecondOrderAssembler.h"
#include "ElInfo.h"
+#include "OpenMP.h"
namespace AMDiS {
@@ -218,9 +219,9 @@ namespace AMDiS {
* Checks whether this is a new travese.
*/
inline void checkForNewTraverse() {
- if (lastTraverseId != ElInfo::traverseId) {
+ if (lastTraverseId != ElInfo::traverseId[omp_get_thread_num()]) {
lastVecEl = lastMatEl = NULL;
- lastTraverseId = ElInfo::traverseId;
+ lastTraverseId = ElInfo::traverseId[omp_get_thread_num()];
}
};
diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc
index fe9f303d3aa4d735f234f34816b232c0471a9cd2..6575877e367b27bfea14feeb86b5c54121d1af68 100755
--- a/AMDiS/src/DOFAdmin.cc
+++ b/AMDiS/src/DOFAdmin.cc
@@ -15,13 +15,6 @@ namespace AMDiS {
const int DOFAdmin::sizeIncrement = 10;
- void DOFAdmin::init()
- {
- firstHole = size = usedCount = holeCount = sizeUsed = 0;
- dofFree.clear();
- }
-
-
DOFAdmin::DOFAdmin(Mesh* m)
: mesh(m),
nrDOF(mesh->getDim(), NO_INIT),
@@ -39,6 +32,12 @@ namespace AMDiS {
init();
}
+ void DOFAdmin::init()
+ {
+ firstHole = size = usedCount = holeCount = sizeUsed = 0;
+ dofFree.clear();
+ }
+
DOFAdmin& DOFAdmin::operator=(const DOFAdmin& src)
{
if (this != &src) {
@@ -75,7 +74,8 @@ namespace AMDiS {
DOFAdmin::DOFAdmin(const DOFAdmin&)
- {}
+ {
+ }
void DOFAdmin::freeDOFIndex(int dof) {
FUNCNAME("DOFAdmin::freeDOFIndex()");
@@ -184,29 +184,42 @@ namespace AMDiS {
void DOFAdmin::addDOFIndexed(DOFIndexedBase* dofIndexed) {
FUNCNAME("DOFAdmin::addDOFIndexed()");
- TEST_EXIT_DBG(dofIndexed)("no dofIndexed\n");
- if (dofIndexed->getSize() < size) {
- dofIndexed->resize(size);
- }
+ TEST_EXIT(dofIndexed)("no dofIndexed\n");
- dofIndexedList.push_back(dofIndexed);
+#ifdef _OPENMP
+#pragma omp critical (dofIndexAccess)
+#endif
+ {
+ if (dofIndexed->getSize() < size) {
+ dofIndexed->resize(size);
+ }
+
+ dofIndexedList.push_back(dofIndexed);
+ }
}
void DOFAdmin::removeDOFIndexed(DOFIndexedBase* dofIndexed)
{
FUNCNAME("DOFAdmin::removeDOFIndexed()");
- std::list<DOFIndexedBase*>::iterator it;
- std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end();
- for (it = dofIndexedList.begin(); it != end; ++it) {
- if (*it == dofIndexed) {
- dofIndexedList.erase(it);
- return;
+ bool removed = false;
+#ifdef _OPENMP
+#pragma omp critical (dofIndexAccess)
+#endif
+ {
+ std::list<DOFIndexedBase*>::iterator it;
+ std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end();
+ for (it = dofIndexedList.begin(); it != end; ++it) {
+ if (*it == dofIndexed) {
+ dofIndexedList.erase(it);
+ removed = true;
+ break;
+ }
}
}
- ERROR("DOFIndexed not in list\n");
+ TEST_EXIT(removed)("DOFIndexed not in list\n");
}
void DOFAdmin::addDOFContainer(DOFContainer* cont)
@@ -315,7 +328,8 @@ namespace AMDiS {
}
DOFAdmin::~DOFAdmin()
- {}
+ {
+ }
void DOFAdmin::serialize(std::ostream &out)
{
diff --git a/AMDiS/src/DOFAdmin.h b/AMDiS/src/DOFAdmin.h
index bb215a856f4bf55132b4724f22fcaa1183e4efcd..097fe7dc82a77f015f168c0c9240261299118591 100644
--- a/AMDiS/src/DOFAdmin.h
+++ b/AMDiS/src/DOFAdmin.h
@@ -36,6 +36,7 @@
#include "FixVec.h"
#include "MemoryManager.h"
#include "Serializable.h"
+#include "OpenMP.h"
#include <vector>
#include <memory>
#include <list>
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 2398e57f14e3f70534f1028c2959234a5045519e..ec34e98bff01e3817f4141c0812846b71fdc0733 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -441,10 +441,10 @@ namespace AMDiS {
if (!op && operators.size() == 0) {
return;
}
-
+
Operator *operat = op ? op : operators[0];
operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
-
+
if (op) {
op->getElementMatrix(elInfo, elementMatrix);
} else {
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index 62a422c963cf1a5f2c879421d5e4ad5e413d221a..4e064fe7e5427829a100c80e00f10321197a263c 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -14,7 +14,8 @@
namespace AMDiS {
- int ElInfo::traverseId = -1;
+ int ElInfo::traverseId[16] = {-1, -1, -1, -1, -1, -1, -1, -1,
+ -1, -1, -1, -1, -1, -1, -1, -1};
int ElInfo::traverseIdCounter = 0;
ElInfo::ElInfo(Mesh *aMesh)
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index aeaea4af7e524408c9f1129060d407a3f4443216..f72bf4c5ad3e7f1b2171f9c124c5b136915614e0 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -586,9 +586,10 @@ namespace AMDiS {
static const int childEdge[3][2][6];
/** \brief
- * Used to determine to which traverse an ElInfo belongs
+ * Used to determine to which traverse an ElInfo belongs.
+ * Stores the value for each thread, maximum value 16!.
*/
- static int traverseId;
+ static int traverseId[16];
/** \brief
* Used to determine to which traverse an ElInfo belongs
diff --git a/AMDiS/src/Error.h b/AMDiS/src/Error.h
index d7988b1c7e68ae9203a2bdc101ed9b5324517f8e..65583b8868809625244950221e897d023e967506 100644
--- a/AMDiS/src/Error.h
+++ b/AMDiS/src/Error.h
@@ -83,9 +83,9 @@ namespace AMDiS {
static int relFct(ElInfo* elinfo);
public:
- static const T& errUFct(const DimVec<double>& lambda);
+ static T errUFct(const DimVec<double>& lambda);
- static const WorldVector<T>& grdErrUFct(const DimVec<double>& lambda);
+ static WorldVector<T> grdErrUFct(const DimVec<double>& lambda);
/** \brief
*
@@ -95,12 +95,12 @@ namespace AMDiS {
public:
AbstrFctErrU() : AbstractFunction<T, DimVec<double> >(0) {};
- inline const T& operator()(const DimVec<double> & x) const {
- return Error<T>::errUFct(x);
+ inline T operator()(const DimVec<double> & x) const {
+ return Error<T>::errUFct(x);
};
};
- static AbstrFctErrU errU;
+ static AbstrFctErrU errU;
/** \brief
*
@@ -110,8 +110,8 @@ namespace AMDiS {
public:
AbstrFctGrdErrU() : AbstractFunction<WorldVector<T>, DimVec<double> >(0) {};
- inline const WorldVector<T>& operator()(const DimVec<double> & x) const {
- return Error<T> :: grdErrUFct (x);
+ inline WorldVector<T> operator()(const DimVec<double> & x) const {
+ return Error<T>::grdErrUFct(x);
};
};
diff --git a/AMDiS/src/Error.hh b/AMDiS/src/Error.hh
index 91cc803f9a97567d7652185c8bc62cfda85aa2c0..4179336f01bde4c78b378657dfae2b71b05d0b39 100644
--- a/AMDiS/src/Error.hh
+++ b/AMDiS/src/Error.hh
@@ -5,26 +5,18 @@
namespace AMDiS {
template<typename T>
- const T& Error<T>::errUFct(const DimVec<double>& lambda)
+ T Error<T>::errUFct(const DimVec<double>& lambda)
{
WorldVector<double> x;
- // if (el_parametric) {
- // ERROR_EXIT("not yet\n");
- // } else {
elinfo->coordToWorld(lambda, &x);
- // }
return((*pU)(x));
}
template<typename T>
- const WorldVector<T>& Error<T>::grdErrUFct(const DimVec<double>& lambda)
+ WorldVector<T> Error<T>::grdErrUFct(const DimVec<double>& lambda)
{
WorldVector<double> x;
- // if (el_parametric) {
- // ERROR_EXIT("not yet\n");
- // } else {
elinfo->coordToWorld(lambda, &x);
- // }
return((*pGrdU)(x));
}
@@ -32,23 +24,11 @@ namespace AMDiS {
template<typename T>
int Error<T>::maxErrAtQpFct(ElInfo* el_info)
{
- int i;
double err;
const double *u_vec, *uh_vec;
- // Parametric *parametric = el_info->getMesh()->getParametric();
elinfo = el_info;
- // if (parametric && parametric->initElement(el_info)) {
- // el_parametric = parametric;
- // } else {
- // el_parametric = NULL;
- // }
-
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
-
- //uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
- //uh_vec = quadFast->uhAtQp(uh_el, NULL);
-
uh_vec = errUh->getVecAtQPs(el_info,
NULL,
quadFast,
@@ -56,11 +36,10 @@ namespace AMDiS {
int numPoints = quadFast->getNumPoints();
- for (i = 0; i < numPoints; i++)
- {
- err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i];
- maxErr = max(maxErr, err);
- }
+ for (int i = 0; i < numPoints; i++) {
+ err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i];
+ maxErr = max(maxErr, err);
+ }
return;
}
@@ -70,14 +49,14 @@ namespace AMDiS {
const DOFVector<T>& uh,
const Quadrature* q)
{
- FUNCNAME("Error<T>::maxErrAtQp");
+ FUNCNAME("Error<T>::maxErrAtQp()");
const FiniteElemSpace *fe_space;
- if (!(pU = &u))
- {
- ERROR("no function u specified; doing nothing\n");
- return(-1.0);
- }
+ if (!(pU = &u)) {
+ ERROR("no function u specified; doing nothing\n");
+ return(-1.0);
+ }
+
if (!(errUh = &uh) || !(fe_space = uh->getFESpace()))
{
ERROR("no discrete function or no fe_space for it; doing nothing\n");
@@ -113,32 +92,11 @@ namespace AMDiS {
{
int i, j;
double err, err_2, h1_err_el, norm_el, norm2, det, exact;
- //int dim = el_info->getMesh()->getDim();
- //const DimVec<WorldVector<double> > &Lambda = el_info->getGrdLambda();
- //const T *uh_el;
const WorldVector<double> *grdu_vec, *grduh_vec;
- // const Parametric *parametric = el_info->getMesh()->getParametric();
elinfo = el_info;
- // if (parametric && parametric->initElement(el_info)) {
- // el_parametric = parametric;
- // }
- // else {
- // el_parametric = NULL;
- // }
-
grdu_vec = quadFast->getQuadrature()->grdFAtQp(grdErrU, NULL);
-
- //uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
-
- // if (el_parametric) {
- // el_parametric->grd_lambda(el_info, NULL, 1, NULL, &Lambda, &det);
- // }
- // else {
det = el_info->getDet();
- // }
-
- //grduh_vec = quadFast->grdUhAtQp(Lambda, uh_el, NULL);
grduh_vec = errUh->getGrdAtQPs(elinfo, NULL, quadFast, NULL);
int numPoints = quadFast->getNumPoints();
@@ -177,7 +135,6 @@ namespace AMDiS {
double Error<T>::H1Err(
const AbstractFunction<WorldVector<T>, WorldVector<double> >& grdU,
const DOFVector<T>& uh,
- //const Quadrature* quad,
int relErr,
double* max,
bool writeLeafData,
@@ -212,7 +169,6 @@ namespace AMDiS {
int deg = grdU.getDegree();
int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
- //if (!quad)
q = Quadrature::provideQuadrature(dim, degree);
quadFast = FastQuadrature::provideFastQuadrature(basFct,
*q,
@@ -250,34 +206,13 @@ namespace AMDiS {
int i;
double err, det, l2_err_el, norm_el, exact;
const double *u_vec, *uh_vec;
- // Parametric *parametric = const_cast<Parametric*>(el_info->getMesh()->getParametric());
-
elinfo = el_info;
- // if (parametric && parametric->initElement(el_info)) {
- // el_parametric = parametric;
- // }
- // else {
- // el_parametric = NULL;
- // }
-
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
-
- //uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
- //uh_vec = quadFast->uhAtQp(uh_el, NULL);
-
uh_vec = errUh->getVecAtQPs(el_info,
NULL,
quadFast,
NULL);
-
- // if (el_parametric) {
- // el_parametric->det(el_info, NULL, 1, NULL, &det);
- // }
- // else {
- // det = el_info->calcDet();
- // }
-
det = el_info->getDet();
int numPoints = quadFast->getNumPoints();
diff --git a/AMDiS/src/FileWriter.cc b/AMDiS/src/FileWriter.cc
index 1f715b00c71d4b5563f0eb91451b35228e0ff7d4..aa53fe1dc621018556c294fe2768b6a5d5b94753 100644
--- a/AMDiS/src/FileWriter.cc
+++ b/AMDiS/src/FileWriter.cc
@@ -28,8 +28,6 @@ namespace AMDiS {
solutionVecs_.resize(1);
solutionVecs_[0] = vec;
-
- dataCollectors_.resize(1);
}
@@ -50,8 +48,6 @@ namespace AMDiS {
feSpace = vecs[0]->getFESpace();
solutionVecs_ = vecs;
-
- dataCollectors_.resize(vecs.size());
}
@@ -83,8 +79,6 @@ namespace AMDiS {
}
feSpace = vec->getFESpace();
-
- dataCollectors_.resize(nTmpSolutions_);
}
@@ -97,12 +91,6 @@ namespace AMDiS {
DELETE solutionVecs_[i];
}
}
-
- for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
- if (dataCollectors_[i]) {
- DELETE dataCollectors_[i];
- }
- }
}
@@ -180,15 +168,18 @@ namespace AMDiS {
ERROR_EXIT("This should not happen!\n");
}
+ // Containers, which store the data to be written;
+ std::vector< DataCollector* > dataCollectors(solutionVecs_.size());
+
if (writeElem) {
- for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
- dataCollectors_[i] = NEW DataCollector(feSpace, solutionVecs_[i],
- level, flag, writeElem);
+ for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
+ dataCollectors[i] = NEW DataCollector(feSpace, solutionVecs_[i],
+ level, flag, writeElem);
}
} else {
- for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
- dataCollectors_[i] = NEW DataCollector(feSpace, solutionVecs_[i],
- traverseLevel, flag | traverseFlag, writeElement);
+ for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
+ dataCollectors[i] = NEW DataCollector(feSpace, solutionVecs_[i],
+ traverseLevel, flag | traverseFlag, writeElement);
}
}
@@ -214,8 +205,8 @@ namespace AMDiS {
ERROR_EXIT("Delay writing only supported for ParaView file format!\n");
}
- for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
- dataCollectors_[i]->fillAllData();
+ for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
+ dataCollectors[i]->fillAllData();
}
writingIsDelayed_ = true;
@@ -236,26 +227,26 @@ namespace AMDiS {
if (writeAMDiSFormat) {
TEST_EXIT(mesh)("no mesh\n");
- MacroWriter::writeMacro(dataCollectors_[0],
+ MacroWriter::writeMacro(dataCollectors[0],
const_cast<char*>((fn + amdisMeshExt).c_str()),
adaptInfo ? adaptInfo->getTime() : 0.0);
MSG("macro file written to %s\n", (fn + amdisMeshExt).c_str());
- ValueWriter::writeValues(dataCollectors_[0],
+ ValueWriter::writeValues(dataCollectors[0],
(fn + amdisDataExt).c_str(),
adaptInfo ? adaptInfo->getTime() : 0.0);
MSG("value file written to %s\n", (fn + amdisDataExt).c_str());
}
if (writePeriodicFormat) {
- MacroWriter::writePeriodicFile(dataCollectors_[0],
+ MacroWriter::writePeriodicFile(dataCollectors[0],
(fn + periodicFileExt).c_str());
MSG("periodic file written to %s\n", (fn + periodicFileExt).c_str());
}
if (writeParaViewFormat) {
- VtkWriter vtkWriter(&dataCollectors_);
+ VtkWriter vtkWriter(&dataCollectors);
vtkWriter.setCompression(compression);
vtkWriter.writeFile(const_cast<char*>((fn + paraViewFileExt).c_str()));
@@ -263,41 +254,20 @@ namespace AMDiS {
}
if (writeParaViewAnimation) {
- VtkWriter vtkWriter(&dataCollectors_);
+ VtkWriter vtkWriter(&dataCollectors);
vtkWriter.updateAnimationFile(fn + paraViewFileExt,
¶ViewAnimationFrames_,
const_cast<char*>((filename + ".pvd").c_str()));
}
- for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
- DELETE dataCollectors_[i];
+ for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
+ DELETE dataCollectors[i];
}
}
void FileWriter::writeDelayedFiles()
{
- if (!writingIsDelayed_) {
- return;
- }
-
- if (writeParaViewFormat) {
- VtkWriter vtkWriter(&dataCollectors_);
- vtkWriter.writeFile(const_cast<char*>((delayedFilename_ + paraViewFileExt).c_str()));
- }
-
- if (writeParaViewAnimation) {
- VtkWriter vtkWriter(&dataCollectors_);
- vtkWriter.updateAnimationFile(delayedFilename_ + paraViewFileExt,
- ¶ViewAnimationFrames_,
- const_cast<char*>((filename + ".pvd").c_str()));
- }
-
- for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
- DELETE dataCollectors_[i];
- }
-
- writingIsDelayed_ = false;
- delayedFilename_ = "";
+ ERROR_EXIT("no more!\n");
}
}
diff --git a/AMDiS/src/FileWriter.h b/AMDiS/src/FileWriter.h
index 7d748f6d7df6fce755505976a823969dc0d3b1cb..40469ecc759460fe7169e7dd2b47f48652a65716 100644
--- a/AMDiS/src/FileWriter.h
+++ b/AMDiS/src/FileWriter.h
@@ -297,11 +297,6 @@ namespace AMDiS {
*/
int nTmpSolutions_;
- /** \brief
- * Containers, which store the data to be written;
- */
- std::vector< DataCollector* > dataCollectors_;
-
/** \brief
* If set to 1, the FileWriter will delay the file writing to the future, where
* it can be executed in parallel with some other independent calculations.
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 05181fc9f7dfdcd287dedcfdb58ade0af00fe8c6..868318f6a19fe5551d0ef70271db3a876e951b83 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -159,11 +159,16 @@ namespace AMDiS {
const double *phi;
if (firstCall) {
- const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
- basFcts = owner->getColFESpace()->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
+ basFcts = owner->getColFESpace()->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
int nPoints = quadrature->getNumPoints();
@@ -205,12 +210,17 @@ namespace AMDiS {
const double *values;
if (firstCall) {
- q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(),
- owner->getColFESpace()->getBasisFcts(),
- quadrature);
- q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
- quadrature);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(),
+ owner->getColFESpace()->getBasisFcts(),
+ quadrature);
+ q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
+ quadrature);
+ firstCall = false;
+ }
}
const int **nEntries = q10->getNumberEntries();
@@ -310,11 +320,16 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > *grdPhi;
if (firstCall) {
- const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
- basFcts = owner->getColFESpace()->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
+ basFcts = owner->getColFESpace()->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
+ firstCall = false;
+ }
}
int nPoints = quadrature->getNumPoints();
@@ -346,11 +361,16 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > *grdPsi;
if (firstCall) {
- const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
- basFcts = owner->getColFESpace()->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
+ basFcts = owner->getColFESpace()->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
int nPoints = quadrature->getNumPoints();
@@ -388,12 +408,17 @@ namespace AMDiS {
const double *values;
if (firstCall) {
- q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFESpace()->getBasisFcts(),
- owner->getColFESpace()->getBasisFcts(),
- quadrature);
- q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
- quadrature);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFESpace()->getBasisFcts(),
+ owner->getColFESpace()->getBasisFcts(),
+ quadrature);
+ q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
+ quadrature);
+ firstCall = false;
+ }
}
const int **nEntries = q01->getNumberEntries();
@@ -427,12 +452,17 @@ namespace AMDiS {
const double *values;
if (firstCall) {
- q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(),
- owner->getColFESpace()->getBasisFcts(),
- quadrature);
- q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
- quadrature);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(),
+ owner->getColFESpace()->getBasisFcts(),
+ quadrature);
+ q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
+ quadrature);
+ firstCall = false;
+ }
}
const int *nEntries = q1->getNumberEntries();
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index a10e8c12c4d4dc9edb1e59c4089aa86c83557289..5141e7eb3bc4d74c942e722fda2d53d6ddfda020 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -74,7 +74,7 @@ namespace AMDiS {
double factor)
{
int i, j, k, l;
- static const int dimOfWorld = Global::getGeo(WORLD);
+ const int dimOfWorld = Global::getGeo(WORLD);
int dim = LALt.getNumRows() - 1;
double val = 0.0;
@@ -128,7 +128,7 @@ namespace AMDiS {
DimMat<double>& LALt,
double factor)
{
- static const int dimOfWorld = Global::getGeo(WORLD);
+ const int dimOfWorld = Global::getGeo(WORLD);
int dim = LALt.getNumRows() - 1;
double val = 0.0;
@@ -157,7 +157,7 @@ namespace AMDiS {
{
int i, j;
int dim = Lb.getSize() - 1;
- static const int dimOfWorld = Global::getGeo(WORLD);
+ const int dimOfWorld = Global::getGeo(WORLD);
double val = 0.0;
@@ -253,16 +253,21 @@ namespace AMDiS {
Quadrature *quad0)
{
if (!assembler[rank])
- 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);
+#ifdef _OPENMP
+#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);
+ }
}
}
diff --git a/AMDiS/src/QPsiPhi.cc b/AMDiS/src/QPsiPhi.cc
index cb54ab5db381242d0b4b9e1d55357b20408ede40..2ccb4d74e6413e3e351e8e0b7eb3a69bd3144977 100644
--- a/AMDiS/src/QPsiPhi.cc
+++ b/AMDiS/src/QPsiPhi.cc
@@ -544,10 +544,13 @@ namespace AMDiS {
Q00PsiPhi::Q00PsiPhi(const BasisFunction *ps,
const BasisFunction *ph,
- const Quadrature *quadrat)
- : psi(ps), phi(ph), quadrature(quadrat)
+ const Quadrature *quadrat)
+ : psi(ps),
+ phi(ph),
+ quadrature(quadrat)
{
- FUNCNAME("Q00PsiPhi::Q00PsiPhi");
+ FUNCNAME("Q00PsiPhi::Q00PsiPhi()");
+
const FastQuadrature *q_phi, *q_psi;
int i, j,iq, n_psi, n_phi;
double val;
@@ -617,7 +620,7 @@ namespace AMDiS {
{
FUNCNAME("Q00PsiPhi::provideQ00PsiPhi");
std::list<Q00PsiPhi*>::iterator list;
-
+
if (!ps && !ph) return NULL;
if (!ps) ps = ph;
diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc
index 25288401e2c08a86f9a99ac1cb715731ff140b30..0c1c5502fae8cd4c61731faebc96434ac48053db 100644
--- a/AMDiS/src/Quadrature.cc
+++ b/AMDiS/src/Quadrature.cc
@@ -1462,30 +1462,37 @@ namespace AMDiS {
{
FUNCNAME("FastQuadrature::FastQuadrature()");
- std::list<FastQuadrature*>::iterator fast = fastQuadList.begin();
-
FastQuadrature *quad_fast;
- for (fast = fastQuadList.begin(); fast != fastQuadList.end(); fast++)
- if ((*fast)->basisFunctions == bas_fcts && (*fast)->quadrature == &quad)
- break;
-
- if ((fast != fastQuadList.end()) && (((*fast)->init_flag & init_flag) == init_flag))
- return(*fast);
-
-
- if (fast == fastQuadList.end()) {
- quad_fast = NEW FastQuadrature(const_cast<BasisFunction*>( bas_fcts), const_cast<Quadrature*>( &quad), 0);
- fastQuadList.push_front(quad_fast);
+#ifdef _OPENMP
+#pragma omp critical (provideFastQuad)
+#endif
+ {
+ std::list<FastQuadrature*>::iterator fast = fastQuadList.begin();
+
+
+ for (fast = fastQuadList.begin(); fast != fastQuadList.end(); fast++)
+ if ((*fast)->basisFunctions == bas_fcts && (*fast)->quadrature == &quad)
+ break;
+
+ if ((fast != fastQuadList.end()) && (((*fast)->init_flag & init_flag) == init_flag)) {
+ quad_fast = *fast;
+ } else {
+ if (fast == fastQuadList.end()) {
+ quad_fast = NEW FastQuadrature(const_cast<BasisFunction*>(bas_fcts), const_cast<Quadrature*>(&quad), 0);
+
+ fastQuadList.push_front(quad_fast);
+
+ max_points = std::max(max_points, quad.getNumPoints());
+ } else {
+ quad_fast = (*fast);
+ }
+ }
- max_points = std::max(max_points, quad.getNumPoints());
- } else {
- quad_fast = (*fast);
+ quad_fast->init(init_flag);
}
- quad_fast->init(init_flag);
-
return quad_fast;
}
diff --git a/AMDiS/src/Recovery.h b/AMDiS/src/Recovery.h
index c34df0a3501078077bc62104e2619a332897b4e7..d5b155e76d83ffc23393e013f934409a37724444 100644
--- a/AMDiS/src/Recovery.h
+++ b/AMDiS/src/Recovery.h
@@ -52,13 +52,12 @@ public:
virtual ~Monomial() {};
- const double& operator()(const WorldVector<double>& y,
- const WorldVector<double>& z) const
+ double operator()(const WorldVector<double>& y,
+ const WorldVector<double>& z) const
{
- static double result;
- result = pow(y[0]-z[0], exponent[0]);
- for (int i=1; i<exponent.size(); i++)
- result *= pow(y[i]-z[i], exponent[i]);
+ double result = pow(y[0] - z[0], exponent[0]);
+ for (int i = 1; i < exponent.size(); i++)
+ result *= pow(y[i] - z[i], exponent[i]);
return result;
};
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 6d1c61701049a79cf0929660846287e21add1bea..ae40b69a0341c9617276aa7c215e84f81f7e6d14 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -49,7 +49,7 @@ namespace AMDiS {
FUNCNAME("TraverseStack::traverseNext()");
ElInfo *elinfo = NULL;
- ElInfo::traverseId = id;
+ ElInfo::traverseId[omp_get_thread_num()] = id;
Parametric *parametric = traverse_mesh->getParametric();
if (stack_used) {
@@ -121,7 +121,7 @@ namespace AMDiS {
Element *el = elinfo->getElement();
int mg_level, sum = 0;
Parametric *parametric = mesh->getParametric();
- ElInfo::traverseId = id;
+ ElInfo::traverseId[omp_get_thread_num()] = id;
if (flag.isSet(Mesh::CALL_LEAF_EL)) {
if (el->getFirstChild()) {
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 01090252986b33667e50d1da43cc0dc8e8da46fc..00b764cf52a46d09e172cef1fa6e447f227a9fb5 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -34,6 +34,7 @@
#include "ElInfo.h"
#include "ElInfoStack.h"
#include "MemoryManager.h"
+#include "OpenMP.h"
#include <vector>
#include <deque>
@@ -72,6 +73,9 @@ namespace AMDiS {
myThreadId_(0),
maxThreads_(1)
{
+#ifdef _OPENMP
+#pragma omp critical
+#endif
id = ElInfo::traverseIdCounter++;
};
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 6ff161b385143d76a6b28bf80503daab209412a0..d9df629ebc8152e707cd049a534611956e3b0008 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -36,9 +36,7 @@ namespace AMDiS {
ZeroOrderAssembler *newAssembler;
std::vector<SubAssembler*> *subAssemblers =
- optimized ?
- &optimizedSubAssemblers :
- &standardSubAssemblers;
+ optimized ? &optimizedSubAssemblers : &standardSubAssemblers;
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
@@ -245,12 +243,17 @@ namespace AMDiS {
int myRank = omp_get_thread_num();
if (firstCall) {
- cPtrs[myRank] = GET_MEMORY(double, nPoints);
- const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
- basFcts = owner->getColFESpace()->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ cPtrs[myRank] = GET_MEMORY(double, nPoints);
+ const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
+ basFcts = owner->getColFESpace()->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
double *c = cPtrs[myRank];
@@ -296,11 +299,16 @@ namespace AMDiS {
void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
if (firstCall) {
- const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
- basFcts = owner->getColFESpace()->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
+ basFcts = owner->getColFESpace()->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
int nPoints = quadrature->getNumPoints();
@@ -335,12 +343,17 @@ namespace AMDiS {
void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
if (firstCall) {
- q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
- owner->getColFESpace()->getBasisFcts(),
- quadrature);
- q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
- quadrature);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
+ owner->getColFESpace()->getBasisFcts(),
+ quadrature);
+ q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
+ quadrature);
+ firstCall = false;
+ }
}
double c = 0.0;
@@ -372,12 +385,17 @@ namespace AMDiS {
void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
if (firstCall) {
- q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
- owner->getColFESpace()->getBasisFcts(),
- quadrature);
- q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
- quadrature);
- firstCall = false;
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+ {
+ q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
+ owner->getColFESpace()->getBasisFcts(),
+ quadrature);
+ q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
+ quadrature);
+ firstCall = false;
+ }
}
std::vector<OperatorTerm*>::iterator termIt;