From 33516d515c089c3ef5a28966c3f04c1d0f331bbf Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Tue, 6 Jun 2017 15:02:14 +0200
Subject: [PATCH] some openMP cleanup
---
AMDiS/src/AdaptInfo.h | 7 +-
AMDiS/src/CouplingProblemStat.h | 13 +-
AMDiS/src/DOFAdmin.cc | 49 +++---
AMDiS/src/FirstOrderAssembler.cc | 90 +++++++----
AMDiS/src/MeshStructure.cc | 16 +-
AMDiS/src/ProblemIterationInterface.h | 6 +
AMDiS/src/ProblemStat.h | 16 +-
AMDiS/src/Quadrature.cc | 216 ++++++++------------------
AMDiS/src/Quadrature.h | 200 ++++++++++++++----------
AMDiS/src/SecondOrderAssembler.cc | 27 ++--
AMDiS/src/StandardProblemIteration.h | 6 +
AMDiS/src/ZeroOrderAssembler.cc | 60 ++++---
AMDiS/src/expressions/valueOf.hpp | 8 +-
13 files changed, 383 insertions(+), 331 deletions(-)
diff --git a/AMDiS/src/AdaptInfo.h b/AMDiS/src/AdaptInfo.h
index b13206d9..b7434de7 100644
--- a/AMDiS/src/AdaptInfo.h
+++ b/AMDiS/src/AdaptInfo.h
@@ -156,11 +156,8 @@ namespace AMDiS {
rosenbrockMode(false)
{
init();
- char number[5];
- for (int i = 0; i < size; i++) {
- sprintf(number, "[%d]", i);
- scalContents[i] = new ScalContent(name + std::string(number));
- }
+ for (int i = 0; i < size; i++)
+ scalContents[i] = new ScalContent(name + "[" + std::to_string(i) + "]");
}
/// Destructor.
diff --git a/AMDiS/src/CouplingProblemStat.h b/AMDiS/src/CouplingProblemStat.h
index a32073a9..d52fceef 100644
--- a/AMDiS/src/CouplingProblemStat.h
+++ b/AMDiS/src/CouplingProblemStat.h
@@ -347,9 +347,20 @@ namespace AMDiS {
{
return problems[number];
}
+
+ virtual ProblemStatType const* getProblem(int number = 0) const
+ {
+ return problems[number];
+ }
/// Returns \ref meshes[i]
- inline Mesh* getMesh(int number = 0)
+ Mesh* getMesh(int number = 0)
+ {
+ return meshes[number];
+ }
+
+ /// Returns \ref meshes[i]
+ Mesh const* getMesh(int number = 0) const
{
return meshes[number];
}
diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc
index f87dabf3..ec7a8ad0 100644
--- a/AMDiS/src/DOFAdmin.cc
+++ b/AMDiS/src/DOFAdmin.cc
@@ -215,10 +215,15 @@ namespace AMDiS {
TEST_EXIT(dofIndexed)("no dofIndexed\n");
- if (dofIndexed->getSize() < size)
- dofIndexed->resize(size);
-
- dofIndexedList.push_back(dofIndexed);
+ #ifdef _OPENMP
+ #pragma omp critical (addDOFIndexed)
+ #endif
+ {
+ if (dofIndexed->getSize() < size)
+ dofIndexed->resize(size);
+
+ dofIndexedList.push_back(dofIndexed);
+ }
}
@@ -226,18 +231,17 @@ namespace AMDiS {
{
FUNCNAME("DOFAdmin::removeDOFIndexed()");
- bool removed = false;
- 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;
+ #ifdef _OPENMP
+ #pragma omp critical (removeDOFIndexed)
+ #endif
+ {
+ auto it = std::find(dofIndexedList.begin(), dofIndexedList.end(), dofIndexed);
+ if (it != dofIndexedList.end())
+ dofIndexedList.erase(it);
+ else {
+ ERROR_EXIT("DOFIndexed not in list!\n");
}
}
-
- TEST_EXIT(removed)("DOFIndexed not in list!\n");
}
@@ -255,16 +259,17 @@ namespace AMDiS {
{
FUNCNAME("DOFAdmin::removeDOFContainer()");
- std::list<DOFContainer*>::iterator it;
- std::list<DOFContainer*>::iterator end = dofContainerList.end();
- for (it = dofContainerList.begin(); it != end; ++it) {
- if (*it == cont) {
- dofContainerList.erase(it);
- return;
+ #ifdef _OPENMP
+ #pragma omp critical (removeDOFContainer)
+ #endif
+ {
+ auto it = std::find(dofContainerList.begin(), dofContainerList.end(), cont);
+ if (it != dofContainerList.end())
+ dofContainerList.erase(it);
+ else {
+ ERROR_EXIT("Container not in list!\n");
}
}
-
- ERROR("Container not in list!\n");
}
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index e332728e..99686387 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -204,11 +204,16 @@ namespace AMDiS {
void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
if (firstCall) {
- const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
- basFcts = colFeSpace->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
+ basFcts = colFeSpace->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
int nPoints = quadrature->getNumPoints();
@@ -240,11 +245,16 @@ namespace AMDiS {
void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
if (firstCall) {
- const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
- basFcts = colFeSpace->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
+ basFcts = colFeSpace->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
int nPoints = quadrature->getNumPoints();
@@ -280,11 +290,16 @@ namespace AMDiS {
const double *values;
if (firstCall) {
- q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
- colFeSpace->getBasisFcts(),
- quadrature);
- q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
+ colFeSpace->getBasisFcts(),
+ quadrature);
+ q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
+ firstCall = false;
+ }
}
const int **nEntries = q10->getNumberEntries();
@@ -362,11 +377,16 @@ namespace AMDiS {
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
if (firstCall) {
- const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
- basFcts = colFeSpace->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
+ basFcts = colFeSpace->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
+ firstCall = false;
+ }
}
int nPoints = quadrature->getNumPoints();
@@ -410,11 +430,16 @@ namespace AMDiS {
const double *values;
if (firstCall) {
- q01 = Q01PsiPhi::provideQ01PsiPhi(rowFeSpace->getBasisFcts(),
- colFeSpace->getBasisFcts(),
- quadrature);
- q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ q01 = Q01PsiPhi::provideQ01PsiPhi(rowFeSpace->getBasisFcts(),
+ colFeSpace->getBasisFcts(),
+ quadrature);
+ q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
+ firstCall = false;
+ }
}
const int **nEntries = q01->getNumberEntries();
@@ -446,11 +471,16 @@ namespace AMDiS {
const double *values;
if (firstCall) {
- q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
- colFeSpace->getBasisFcts(),
- quadrature);
- q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
+ colFeSpace->getBasisFcts(),
+ quadrature);
+ q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
+ firstCall = false;
+ }
}
const int *nEntries = q1->getNumberEntries();
diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc
index b23b4911..4cadb1d4 100644
--- a/AMDiS/src/MeshStructure.cc
+++ b/AMDiS/src/MeshStructure.cc
@@ -280,25 +280,25 @@ namespace AMDiS {
bool cont = true;
while (cont) {
bool cont1;
-#ifndef NDEBUG
+// #ifndef NDEBUG
bool cont2;
-#endif
+// #endif
if (structure1->isLeafElement() == structure2->isLeafElement()) {
cont1 = structure1->nextElement(result);
-#ifndef NDEBUG
+// #ifndef NDEBUG
cont2 = structure2->nextElement();
-#endif
+// #endif
} else {
if (structure1->isLeafElement()) {
cont1 = structure1->nextElement();
-#ifndef NDEBUG
+// #ifndef NDEBUG
cont2 = structure2->skipBranch(result);
-#endif
+// #endif
} else {
cont1 = structure1->skipBranch(result);
-#ifndef NDEBUG
+// #ifndef NDEBUG
cont2 = structure2->nextElement();
-#endif
+// #endif
}
}
TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
diff --git a/AMDiS/src/ProblemIterationInterface.h b/AMDiS/src/ProblemIterationInterface.h
index b3240ff6..7e0c07a5 100644
--- a/AMDiS/src/ProblemIterationInterface.h
+++ b/AMDiS/src/ProblemIterationInterface.h
@@ -75,6 +75,12 @@ namespace AMDiS {
* is managed by this master problem, the number hasn't to be given.
*/
virtual ProblemStatBase *getProblem(int number = 0) = 0;
+
+ // dirty hack
+ virtual ProblemStatBase const* getProblem(int number = 0) const
+ {
+ return const_cast<ProblemStatBase const*>(const_cast<ProblemIterationInterface*>(this)->getProblem(number));
+ }
/// Returns the problem with the given name.
virtual ProblemStatBase *getProblem(std::string name)
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 449a1de3..f6062518 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -317,6 +317,14 @@ namespace AMDiS {
("invalid component number\n");
return componentMeshes[comp];
}
+
+ Mesh const* getMesh(int comp = 0) const
+ {
+ FUNCNAME("ProblemStatSeq::getMesh()");
+ TEST_EXIT(comp < static_cast<int>(componentMeshes.size()) && comp >= 0)
+ ("invalid component number\n");
+ return componentMeshes[comp];
+ }
/// Returns \ref meshes
inline std::vector<Mesh*> getMeshes()
@@ -325,7 +333,7 @@ namespace AMDiS {
}
/// Returns \ref feSpace_.
- inline const FiniteElemSpace* getFeSpace(int comp = 0)
+ inline const FiniteElemSpace* getFeSpace(int comp = 0) const
{
FUNCNAME("ProblemStatSeq::getFeSpace()");
TEST_EXIT(comp < static_cast<int>(componentSpaces.size()) && comp >= 0)
@@ -726,6 +734,12 @@ namespace AMDiS {
{
return this;
}
+
+ /// Returns the problem with the given number.
+ virtual ProblemStatBase const* getProblem(int number = 0) const override
+ {
+ return this;
+ }
/// Determines the execution order of the single adaption steps. If adapt is
/// true, mesh adaption will be performed. This allows to avoid mesh adaption,
diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc
index dce37e56..bbc08352 100644
--- a/AMDiS/src/Quadrature.cc
+++ b/AMDiS/src/Quadrature.cc
@@ -34,21 +34,13 @@ namespace AMDiS {
int FastQuadrature::max_points = 0;
Quadrature::Quadrature(const Quadrature& q)
- {
- name = q.name;
- degree = q.degree;
- dim = q.dim;
- n_points = q.n_points;
-
- // copy barycentric coordinates
- lambda = new VectorOfFixVecs<DimVec<double> >(*(q.lambda));
-
- // copy weights
- w = new double[n_points];
-
- for (int i = 0; i < n_points; i++)
- w[i] = q.w[i];
- }
+ : name(q.name)
+ , degree(q.degree)
+ , dim(q.dim)
+ , n_points(q.n_points)
+ , lambda(new VectorOfFixVecs<DimVec<double> >(*(q.lambda)))
+ , w(q.w)
+ {}
/****************************************************************************/
/* initialize gradient values of a function f in local coordinates at the */
@@ -115,76 +107,7 @@ namespace AMDiS {
}
- Quadrature **Quadrature::quad_nd[4];
- Quadrature *Quadrature::quad_0d[1];
- Quadrature *Quadrature::quad_1d[20];
- Quadrature *Quadrature::quad_2d[18];
- Quadrature *Quadrature::quad_3d[8];
-
- VectorOfFixVecs<DimVec<double> > *Quadrature::x_0d;
- double *Quadrature::w_0d;
-
- VectorOfFixVecs<DimVec<double> > *Quadrature::x0_1d = NULL;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x1_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x2_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x3_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x4_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x5_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x6_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x7_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x8_1d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x9_1d;
- double *Quadrature::w0_1d;
- double *Quadrature::w1_1d;
- double *Quadrature::w2_1d;
- double *Quadrature::w3_1d;
- double *Quadrature::w4_1d;
- double *Quadrature::w5_1d;
- double *Quadrature::w6_1d;
- double *Quadrature::w7_1d;
- double *Quadrature::w8_1d;
- double *Quadrature::w9_1d;
-
- VectorOfFixVecs<DimVec<double> > *Quadrature::x1_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x2_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x3_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x4_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x5_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x7_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x8_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x9_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x10_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x11_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x12_2d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x17_2d;
- double *Quadrature::w1_2d;
- double *Quadrature::w2_2d;
- double *Quadrature::w3_2d;
- double *Quadrature::w4_2d;
- double *Quadrature::w5_2d;
- double *Quadrature::w7_2d;
- double *Quadrature::w8_2d;
- double *Quadrature::w9_2d;
- double *Quadrature::w10_2d;
- double *Quadrature::w11_2d;
- double *Quadrature::w12_2d;
- double *Quadrature::w17_2d;
-
- VectorOfFixVecs<DimVec<double> > *Quadrature::x1_3d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x2_3d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x3_3d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x4_3d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x5_3d;
- VectorOfFixVecs<DimVec<double> > *Quadrature::x7_3d;
- double *Quadrature::w1_3d;
- double *Quadrature::w2_3d;
- double *Quadrature::w3_3d;
- double *Quadrature::w4_3d;
- double *Quadrature::w5_3d;
- double *Quadrature::w7_3d;
-
-
- void Quadrature::initStaticQuadratures()
+ QuadratureFactory::QuadratureFactory()
{
TEST_EXIT(x0_1d == NULL)("static quadratures already initialized\n");
@@ -201,7 +124,7 @@ namespace AMDiS {
w_0d = createAndInitArray(1, StdVol * 1.0);
- Quadrature::quad_0d[0] = new Quadrature("0d", 0, 0, 1, Quadrature::x_0d, Quadrature::w_0d);
+ quad_0d[0].reset(new Quadrature("0d", 0, 0, 1, x_0d, w_0d));
/****************************************************************************/
/* 1d quadrature formulas using 2 barycentric coordinates */
@@ -386,26 +309,26 @@ namespace AMDiS {
StdVol * 0.074725674575291,
StdVol * 0.033335672154344);
- Quadrature::quad_1d[0]= new Quadrature("1d-Gauss: P_1", 1, 1, 1, Quadrature::x0_1d, Quadrature::w0_1d); /* P_0 */
- Quadrature::quad_1d[1]= new Quadrature("1d-Gauss: P_1", 1, 1, 1, Quadrature::x0_1d, Quadrature::w0_1d); /* P_1 */
- Quadrature::quad_1d[2]= new Quadrature("1d-Gauss: P_3", 3, 1, 2, Quadrature::x1_1d, Quadrature::w1_1d); /* P_2 */
- Quadrature::quad_1d[3]= new Quadrature("1d-Gauss: P_3", 3, 1, 2, Quadrature::x1_1d, Quadrature::w1_1d); /* P_3 */
- Quadrature::quad_1d[4]= new Quadrature("1d-Gauss: P_5", 5, 1, 3, Quadrature::x2_1d, Quadrature::w2_1d); /* P_4 */
- Quadrature::quad_1d[5]= new Quadrature("1d-Gauss: P_5", 5, 1, 3, Quadrature::x2_1d, Quadrature::w2_1d); /* P_5 */
- Quadrature::quad_1d[6]= new Quadrature("1d-Gauss: P_7", 7, 1, 4, Quadrature::x3_1d, Quadrature::w3_1d); /* P_6 */
- Quadrature::quad_1d[7]= new Quadrature("1d-Gauss: P_7", 7, 1, 4, Quadrature::x3_1d, Quadrature::w3_1d); /* P_7 */
- Quadrature::quad_1d[8]= new Quadrature("1d-Gauss: P_9", 9, 1, 5, Quadrature::x4_1d, Quadrature::w4_1d); /* P_8 */
- Quadrature::quad_1d[9]= new Quadrature("1d-Gauss: P_9", 9, 1, 5, Quadrature::x4_1d, Quadrature::w4_1d); /* P_9 */
- Quadrature::quad_1d[10]= new Quadrature("1d-Gauss: P_11", 11, 1, 6, Quadrature::x5_1d, Quadrature::w5_1d); /* P_10 */
- Quadrature::quad_1d[11]= new Quadrature("1d-Gauss: P_11", 11, 1, 6, Quadrature::x5_1d, Quadrature::w5_1d); /* P_11 */
- Quadrature::quad_1d[12]= new Quadrature("1d-Gauss: P_13", 13, 1, 7, Quadrature::x6_1d, Quadrature::w6_1d); /* P_12 */
- Quadrature::quad_1d[13]= new Quadrature("1d-Gauss: P_13", 13, 1, 7, Quadrature::x6_1d, Quadrature::w6_1d); /* P_13 */
- Quadrature::quad_1d[14]= new Quadrature("1d-Gauss: P_15", 15, 1, 8, Quadrature::x7_1d, Quadrature::w7_1d); /* P_14 */
- Quadrature::quad_1d[15]= new Quadrature("1d-Gauss: P_15", 15, 1, 8, Quadrature::x7_1d, Quadrature::w7_1d); /* P_15 */
- Quadrature::quad_1d[16]= new Quadrature("1d-Gauss: P_17", 17, 1, 9, Quadrature::x8_1d, Quadrature::w8_1d); /* P_16 */
- Quadrature::quad_1d[17]= new Quadrature("1d-Gauss: P_17", 17, 1, 9, Quadrature::x8_1d, Quadrature::w8_1d); /* P_17 */
- Quadrature::quad_1d[18]= new Quadrature("1d-Gauss: P_19", 19, 1, 10, Quadrature::x9_1d, Quadrature::w9_1d); /* P_18 */
- Quadrature::quad_1d[19]= new Quadrature("1d-Gauss: P_19", 19, 1, 10, Quadrature::x9_1d, Quadrature::w9_1d);/* P_19 */
+ quad_1d[0].reset(new Quadrature("1d-Gauss: P_1", 1, 1, 1, x0_1d, w0_1d)); /* P_0 */
+ quad_1d[1].reset(new Quadrature("1d-Gauss: P_1", 1, 1, 1, x0_1d, w0_1d)); /* P_1 */
+ quad_1d[2].reset(new Quadrature("1d-Gauss: P_3", 3, 1, 2, x1_1d, w1_1d)); /* P_2 */
+ quad_1d[3].reset(new Quadrature("1d-Gauss: P_3", 3, 1, 2, x1_1d, w1_1d)); /* P_3 */
+ quad_1d[4].reset(new Quadrature("1d-Gauss: P_5", 5, 1, 3, x2_1d, w2_1d)); /* P_4 */
+ quad_1d[5].reset(new Quadrature("1d-Gauss: P_5", 5, 1, 3, x2_1d, w2_1d)); /* P_5 */
+ quad_1d[6].reset(new Quadrature("1d-Gauss: P_7", 7, 1, 4, x3_1d, w3_1d)); /* P_6 */
+ quad_1d[7].reset(new Quadrature("1d-Gauss: P_7", 7, 1, 4, x3_1d, w3_1d)); /* P_7 */
+ quad_1d[8].reset(new Quadrature("1d-Gauss: P_9", 9, 1, 5, x4_1d, w4_1d)); /* P_8 */
+ quad_1d[9].reset(new Quadrature("1d-Gauss: P_9", 9, 1, 5, x4_1d, w4_1d)); /* P_9 */
+ quad_1d[10].reset(new Quadrature("1d-Gauss: P_11", 11, 1, 6, x5_1d, w5_1d)); /* P_10 */
+ quad_1d[11].reset(new Quadrature("1d-Gauss: P_11", 11, 1, 6, x5_1d, w5_1d)); /* P_11 */
+ quad_1d[12].reset(new Quadrature("1d-Gauss: P_13", 13, 1, 7, x6_1d, w6_1d)); /* P_12 */
+ quad_1d[13].reset(new Quadrature("1d-Gauss: P_13", 13, 1, 7, x6_1d, w6_1d)); /* P_13 */
+ quad_1d[14].reset(new Quadrature("1d-Gauss: P_15", 15, 1, 8, x7_1d, w7_1d)); /* P_14 */
+ quad_1d[15].reset(new Quadrature("1d-Gauss: P_15", 15, 1, 8, x7_1d, w7_1d)); /* P_15 */
+ quad_1d[16].reset(new Quadrature("1d-Gauss: P_17", 17, 1, 9, x8_1d, w8_1d)); /* P_16 */
+ quad_1d[17].reset(new Quadrature("1d-Gauss: P_17", 17, 1, 9, x8_1d, w8_1d)); /* P_17 */
+ quad_1d[18].reset(new Quadrature("1d-Gauss: P_19", 19, 1, 10, x9_1d, w9_1d)); /* P_18 */
+ quad_1d[19].reset(new Quadrature("1d-Gauss: P_19", 19, 1, 10, x9_1d, w9_1d));/* P_19 */
#undef StdVol
/****************************************************************************/
@@ -1131,24 +1054,24 @@ namespace AMDiS {
#undef w14
#undef w15
- Quadrature::quad_2d[0] = new Quadrature("2d-P_1", 1, 2, N1, Quadrature::x1_2d, Quadrature::w1_2d); /* P 0 */
- Quadrature::quad_2d[1] = new Quadrature("2d-P_1", 1, 2, N1, Quadrature::x1_2d, Quadrature::w1_2d); /* P 1 */
- Quadrature::quad_2d[2] = new Quadrature("2d Stroud: P_2", 2, 2, N2, Quadrature::x2_2d, Quadrature::w2_2d); /* P 2 */
- Quadrature::quad_2d[3] = new Quadrature("2d Stroud: P_3", 3, 2, N3, Quadrature::x3_2d, Quadrature::w3_2d); /* P 3 */
- Quadrature::quad_2d[4] = new Quadrature("2d Dunavant: P_4", 4, 2, N4, Quadrature::x4_2d, Quadrature::w4_2d); /* P 4 */
- Quadrature::quad_2d[5] = new Quadrature("2d Dunavant: P_5", 5, 2, N5, Quadrature::x5_2d, Quadrature::w5_2d); /* P 5 */
- Quadrature::quad_2d[6] = new Quadrature("2d Gattermann: P_7", 7, 2, N7, Quadrature::x7_2d, Quadrature::w7_2d); /* P 6 */
- Quadrature::quad_2d[7] = new Quadrature("2d Gattermann: P_7", 7, 2, N7, Quadrature::x7_2d, Quadrature::w7_2d); /* P 7 */
- Quadrature::quad_2d[8] = new Quadrature("2d Dunavant: P_8", 8, 2, N8, Quadrature::x8_2d, Quadrature::w8_2d); /* P 8 */
- Quadrature::quad_2d[9] = new Quadrature("2d Dunavant: P_9", 9, 2, N9, Quadrature::x9_2d, Quadrature::w9_2d); /* P 9 */
- Quadrature::quad_2d[10] = new Quadrature("2d Dunavant: P_10", 10, 2, N10, Quadrature::x10_2d, Quadrature::w10_2d);/* P 10 */
- Quadrature::quad_2d[11] = new Quadrature("2d Dunavant: P_11", 11, 2, N11, Quadrature::x11_2d, Quadrature::w11_2d);/* P 11 */
- Quadrature::quad_2d[12] = new Quadrature("2d Dunavant: P_12", 12, 2, N12, Quadrature::x12_2d, Quadrature::w12_2d);/* P 12 */
- Quadrature::quad_2d[13] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 13 */
- Quadrature::quad_2d[14] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 14 */
- Quadrature::quad_2d[15] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 15 */
- Quadrature::quad_2d[16] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 16 */
- Quadrature::quad_2d[17] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d); /* P 17 */
+ quad_2d[0].reset(new Quadrature("2d-P_1", 1, 2, N1, x1_2d, w1_2d)); /* P 0 */
+ quad_2d[1].reset(new Quadrature("2d-P_1", 1, 2, N1, x1_2d, w1_2d)); /* P 1 */
+ quad_2d[2].reset(new Quadrature("2d Stroud: P_2", 2, 2, N2, x2_2d, w2_2d)); /* P 2 */
+ quad_2d[3].reset(new Quadrature("2d Stroud: P_3", 3, 2, N3, x3_2d, w3_2d)); /* P 3 */
+ quad_2d[4].reset(new Quadrature("2d Dunavant: P_4", 4, 2, N4, x4_2d, w4_2d)); /* P 4 */
+ quad_2d[5].reset(new Quadrature("2d Dunavant: P_5", 5, 2, N5, x5_2d, w5_2d)); /* P 5 */
+ quad_2d[6].reset(new Quadrature("2d Gattermann: P_7", 7, 2, N7, x7_2d, w7_2d)); /* P 6 */
+ quad_2d[7].reset(new Quadrature("2d Gattermann: P_7", 7, 2, N7, x7_2d, w7_2d)); /* P 7 */
+ quad_2d[8].reset(new Quadrature("2d Dunavant: P_8", 8, 2, N8, x8_2d, w8_2d)); /* P 8 */
+ quad_2d[9].reset(new Quadrature("2d Dunavant: P_9", 9, 2, N9, x9_2d, w9_2d)); /* P 9 */
+ quad_2d[10].reset(new Quadrature("2d Dunavant: P_10", 10, 2, N10, x10_2d, w10_2d));/* P 10 */
+ quad_2d[11].reset(new Quadrature("2d Dunavant: P_11", 11, 2, N11, x11_2d, w11_2d));/* P 11 */
+ quad_2d[12].reset(new Quadrature("2d Dunavant: P_12", 12, 2, N12, x12_2d, w12_2d));/* P 12 */
+ quad_2d[13].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 13 */
+ quad_2d[14].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 14 */
+ quad_2d[15].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 15 */
+ quad_2d[16].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 16 */
+ quad_2d[17].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d)); /* P 17 */
#undef StdVol
@@ -1401,14 +1324,14 @@ namespace AMDiS {
/* use that of degree (only on function evaluation also) */
/****************************************************************************/
- Quadrature::quad_3d[0] = new Quadrature("3d Stroud: P_1", 1, 3, 1, Quadrature::x1_3d, Quadrature::w1_3d); /* P_0 */
- Quadrature::quad_3d[1] = new Quadrature("3d Stroud: P_1", 1, 3, 1, Quadrature::x1_3d, Quadrature::w1_3d); /* P_1 */
- Quadrature::quad_3d[2] = new Quadrature("3d Stroud: P_2", 2, 3, 4, Quadrature::x2_3d, Quadrature::w2_3d); /* P_2 */
- Quadrature::quad_3d[3] = new Quadrature("3d Stroud: P_3", 3, 3, 8, Quadrature::x3_3d, Quadrature::w3_3d); /* P_3 */
- Quadrature::quad_3d[4] = new Quadrature("3d ???: P_5", 5, 3, 15, Quadrature::x5_3d, Quadrature::w5_3d); /* P_4 */
- Quadrature::quad_3d[5] = new Quadrature("3d ???: P_5", 5, 3, 15, Quadrature::x5_3d, Quadrature::w5_3d); /* P_5 */
- Quadrature::quad_3d[6] = new Quadrature("3d ???: P_7", 7, 3, 64, Quadrature::x7_3d, Quadrature::w7_3d); /* P_6 */
- Quadrature::quad_3d[7] = new Quadrature("3d ???: P_7", 7, 3, 64, Quadrature::x7_3d, Quadrature::w7_3d); /* P_7 */
+ quad_3d[0].reset(new Quadrature("3d Stroud: P_1", 1, 3, 1, x1_3d, w1_3d)); /* P_0 */
+ quad_3d[1].reset(new Quadrature("3d Stroud: P_1", 1, 3, 1, x1_3d, w1_3d)); /* P_1 */
+ quad_3d[2].reset(new Quadrature("3d Stroud: P_2", 2, 3, 4, x2_3d, w2_3d)); /* P_2 */
+ quad_3d[3].reset(new Quadrature("3d Stroud: P_3", 3, 3, 8, x3_3d, w3_3d)); /* P_3 */
+ quad_3d[4].reset(new Quadrature("3d ???: P_5", 5, 3, 15, x5_3d, w5_3d)); /* P_4 */
+ quad_3d[5].reset(new Quadrature("3d ???: P_5", 5, 3, 15, x5_3d, w5_3d)); /* P_5 */
+ quad_3d[6].reset(new Quadrature("3d ???: P_7", 7, 3, 64, x7_3d, w7_3d)); /* P_6 */
+ quad_3d[7].reset(new Quadrature("3d ???: P_7", 7, 3, 64, x7_3d, w7_3d)); /* P_7 */
#undef StdVol
@@ -1417,16 +1340,17 @@ namespace AMDiS {
/* integration in different dimensions */
/****************************************************************************/
- Quadrature::quad_nd[0] = Quadrature::quad_0d;
- Quadrature::quad_nd[1] = Quadrature::quad_1d;
- Quadrature::quad_nd[2] = Quadrature::quad_2d;
- Quadrature::quad_nd[3] = Quadrature::quad_3d;
+// quad_nd[0] = quad_0d;
+// quad_nd[1] = quad_1d;
+// quad_nd[2] = quad_2d;
+// quad_nd[3] = quad_3d;
}
-
+
Quadrature* Quadrature::provideQuadrature(int dim_, int degree_)
{
FUNCNAME("Quadrature::provideQuadrature()");
+ static QuadratureFactory factory{};
switch (dim_) {
case 0:
@@ -1445,17 +1369,7 @@ namespace AMDiS {
ERROR_EXIT("invalid dim\n");
}
- if (x0_1d == NULL)
- initStaticQuadratures();
-
- return (quad_nd[dim_][degree_]);
- }
-
-
- Quadrature::~Quadrature()
- {
- delete lambda;
- delete [] w;
+ return factory.get(dim_, degree_);
}
@@ -1482,8 +1396,10 @@ namespace AMDiS {
Flag init_flag)
{
FastQuadrature *quad_fast = NULL;
-
-// #pragma omp critical
+
+ #ifdef _OPENMP
+ #pragma omp critical (provideFastQuad)
+ #endif
{
list<FastQuadrature*>::iterator fast = fastQuadList.begin();
for (; fast != fastQuadList.end(); fast++)
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index d6a7419e..085764ce 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -25,8 +25,10 @@
#ifndef AMDIS_QUADRATURE_H
#define AMDIS_QUADRATURE_H
+#include <array>
#include <list>
#include <vector>
+#include <memory>
#include <boost/numeric/mtl/mtl.hpp>
#include "AMDiS_fwd.h"
#include "BasisFunction.h"
@@ -34,7 +36,9 @@
#include "FixVec.h"
namespace AMDiS {
-
+
+ class QuadratureFactory;
+
/**
* \ingroup Assembler
*
@@ -50,6 +54,8 @@ namespace AMDiS {
*/
class Quadrature
{
+ friend class QuadratureFactory;
+
protected:
/// Avoids call of default constructor
Quadrature();
@@ -71,7 +77,7 @@ namespace AMDiS {
dim(dim_),
n_points(n_points_),
lambda(lambda_),
- w(w_)
+ w(w_, w_+n_points_)
{}
public:
@@ -79,7 +85,7 @@ namespace AMDiS {
Quadrature(const Quadrature&);
/// Destructor
- ~Quadrature();
+// ~Quadrature();
/// Returns a Quadrature for dimension dim exact for degree degree.
static Quadrature *provideQuadrature(int dim, int degree);
@@ -115,9 +121,15 @@ namespace AMDiS {
}
/// Returns \ref w.
- inline double* getWeight() const
+ inline double const* getWeight() const
+ {
+ return w.data();
+ }
+
+ /// Returns \ref w.
+ inline double* getWeight()
{
- return w;
+ return w.data();
}
/// Returns \ref dim
@@ -199,84 +211,114 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > *lambda;
/// Vector of quadrature weights
- double *w;
-
+ std::vector<double> w;
+ };
+
+
+ class QuadratureFactory
+ {
+ friend class Quadrature;
+
+ private:
+
+ QuadratureFactory();
+
+ public:
+
+ Quadrature* get(const int dim, const int degree) const
+ {
+ switch (dim) {
+ case 0:
+ return quad_0d[degree].get(); break;
+ case 1:
+ return quad_1d[degree].get(); break;
+ case 2:
+ return quad_2d[degree].get(); break;
+ case 3:
+ return quad_3d[degree].get(); break;
+ default:
+ ERROR_EXIT("Unknown dimension!");
+ return NULL;
+ }
+ }
+
protected:
- /// Initialisation of all static Quadrature objects which will be returned
- /// by \ref provideQuadrature()
- static void initStaticQuadratures();
- /** \name static quadratures, used weights, and barycentric coords
+ /** \name (static) quadratures, used weights, and barycentric coords
* \{
*/
- static Quadrature **quad_nd[4];
- static Quadrature *quad_0d[1];
- static Quadrature *quad_1d[20];
- static Quadrature *quad_2d[18];
- static Quadrature *quad_3d[8];
-
- static VectorOfFixVecs<DimVec<double> > *x_0d;
- static double *w_0d;
-
- static VectorOfFixVecs<DimVec<double> > *x0_1d;
- static VectorOfFixVecs<DimVec<double> > *x1_1d;
- static VectorOfFixVecs<DimVec<double> > *x2_1d;
- static VectorOfFixVecs<DimVec<double> > *x3_1d;
- static VectorOfFixVecs<DimVec<double> > *x4_1d;
- static VectorOfFixVecs<DimVec<double> > *x5_1d;
- static VectorOfFixVecs<DimVec<double> > *x6_1d;
- static VectorOfFixVecs<DimVec<double> > *x7_1d;
- static VectorOfFixVecs<DimVec<double> > *x8_1d;
- static VectorOfFixVecs<DimVec<double> > *x9_1d;
- static double *w0_1d;
- static double *w1_1d;
- static double *w2_1d;
- static double *w3_1d;
- static double *w4_1d;
- static double *w5_1d;
- static double *w6_1d;
- static double *w7_1d;
- static double *w8_1d;
- static double *w9_1d;
-
- static VectorOfFixVecs<DimVec<double> > *x1_2d;
- static VectorOfFixVecs<DimVec<double> > *x2_2d;
- static VectorOfFixVecs<DimVec<double> > *x3_2d;
- static VectorOfFixVecs<DimVec<double> > *x4_2d;
- static VectorOfFixVecs<DimVec<double> > *x5_2d;
- static VectorOfFixVecs<DimVec<double> > *x7_2d;
- static VectorOfFixVecs<DimVec<double> > *x8_2d;
- static VectorOfFixVecs<DimVec<double> > *x9_2d;
- static VectorOfFixVecs<DimVec<double> > *x10_2d;
- static VectorOfFixVecs<DimVec<double> > *x11_2d;
- static VectorOfFixVecs<DimVec<double> > *x12_2d;
- static VectorOfFixVecs<DimVec<double> > *x17_2d;
- static double *w1_2d;
- static double *w2_2d;
- static double *w3_2d;
- static double *w4_2d;
- static double *w5_2d;
- static double *w7_2d;
- static double *w8_2d;
- static double *w9_2d;
- static double *w10_2d;
- static double *w11_2d;
- static double *w12_2d;
- static double *w17_2d;
-
- static VectorOfFixVecs<DimVec<double> > *x1_3d;
- static VectorOfFixVecs<DimVec<double> > *x2_3d;
- static VectorOfFixVecs<DimVec<double> > *x3_3d;
- static VectorOfFixVecs<DimVec<double> > *x4_3d;
- static VectorOfFixVecs<DimVec<double> > *x5_3d;
- static VectorOfFixVecs<DimVec<double> > *x7_3d;
- static double *w1_3d;
- static double *w2_3d;
- static double *w3_3d;
- static double *w4_3d;
- static double *w5_3d;
- static double *w7_3d;
-
+// Quadrature **quad_nd[4];
+ std::array<std::unique_ptr<Quadrature>,1> quad_0d;
+ std::array<std::unique_ptr<Quadrature>,20> quad_1d;
+ std::array<std::unique_ptr<Quadrature>,18> quad_2d;
+ std::array<std::unique_ptr<Quadrature>,8> quad_3d;
+// Quadrature *quad_0d[1];
+// Quadrature *quad_1d[20];
+// Quadrature *quad_2d[18];
+// Quadrature *quad_3d[8];
+
+ VectorOfFixVecs<DimVec<double> > *x_0d;
+ double *w_0d;
+
+ VectorOfFixVecs<DimVec<double> > *x0_1d;
+ VectorOfFixVecs<DimVec<double> > *x1_1d;
+ VectorOfFixVecs<DimVec<double> > *x2_1d;
+ VectorOfFixVecs<DimVec<double> > *x3_1d;
+ VectorOfFixVecs<DimVec<double> > *x4_1d;
+ VectorOfFixVecs<DimVec<double> > *x5_1d;
+ VectorOfFixVecs<DimVec<double> > *x6_1d;
+ VectorOfFixVecs<DimVec<double> > *x7_1d;
+ VectorOfFixVecs<DimVec<double> > *x8_1d;
+ VectorOfFixVecs<DimVec<double> > *x9_1d;
+ double *w0_1d;
+ double *w1_1d;
+ double *w2_1d;
+ double *w3_1d;
+ double *w4_1d;
+ double *w5_1d;
+ double *w6_1d;
+ double *w7_1d;
+ double *w8_1d;
+ double *w9_1d;
+
+ VectorOfFixVecs<DimVec<double> > *x1_2d;
+ VectorOfFixVecs<DimVec<double> > *x2_2d;
+ VectorOfFixVecs<DimVec<double> > *x3_2d;
+ VectorOfFixVecs<DimVec<double> > *x4_2d;
+ VectorOfFixVecs<DimVec<double> > *x5_2d;
+ VectorOfFixVecs<DimVec<double> > *x7_2d;
+ VectorOfFixVecs<DimVec<double> > *x8_2d;
+ VectorOfFixVecs<DimVec<double> > *x9_2d;
+ VectorOfFixVecs<DimVec<double> > *x10_2d;
+ VectorOfFixVecs<DimVec<double> > *x11_2d;
+ VectorOfFixVecs<DimVec<double> > *x12_2d;
+ VectorOfFixVecs<DimVec<double> > *x17_2d;
+ double *w1_2d;
+ double *w2_2d;
+ double *w3_2d;
+ double *w4_2d;
+ double *w5_2d;
+ double *w7_2d;
+ double *w8_2d;
+ double *w9_2d;
+ double *w10_2d;
+ double *w11_2d;
+ double *w12_2d;
+ double *w17_2d;
+
+ VectorOfFixVecs<DimVec<double> > *x1_3d;
+ VectorOfFixVecs<DimVec<double> > *x2_3d;
+ VectorOfFixVecs<DimVec<double> > *x3_3d;
+ VectorOfFixVecs<DimVec<double> > *x4_3d;
+ VectorOfFixVecs<DimVec<double> > *x5_3d;
+ VectorOfFixVecs<DimVec<double> > *x7_3d;
+ double *w1_3d;
+ double *w2_3d;
+ double *w3_3d;
+ double *w4_3d;
+ double *w5_3d;
+ double *w7_3d;
+
/** \} */
};
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index db720fcf..0c400e89 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -187,17 +187,22 @@ namespace AMDiS {
const int nPoints = quadrature->getNumPoints();
if (firstCall) {
- dimVec.change_dim(dim + 1);
- LALt.resize(nPoints);
- for (int j = 0; j < nPoints; j++)
- LALt[j].change_dim(dim + 1, dim + 1);
-
- psiFast =
- updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
- phiFast =
- updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
-
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ dimVec.change_dim(dim + 1);
+ LALt.resize(nPoints);
+ for (int j = 0; j < nPoints; j++)
+ LALt[j].change_dim(dim + 1, dim + 1);
+
+ psiFast =
+ updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
+ phiFast =
+ updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
+
+ firstCall = false;
+ }
}
for (int i = 0; i < nPoints; i++)
diff --git a/AMDiS/src/StandardProblemIteration.h b/AMDiS/src/StandardProblemIteration.h
index 95967c70..024134e8 100644
--- a/AMDiS/src/StandardProblemIteration.h
+++ b/AMDiS/src/StandardProblemIteration.h
@@ -61,6 +61,12 @@ namespace AMDiS {
{
return problem;
}
+
+ /// Implementation of \ref ProblemIterationInterface::getProblem(int)
+ ProblemStatBase const* getProblem(int number = 0) const override
+ {
+ return problem;
+ }
/// Returns the name of the problem.
std::string getName() override;
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 12fa11ed..db100f8d 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -195,11 +195,16 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
if (firstCall) {
- const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
- basFcts = colFeSpace->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
+ basFcts = colFeSpace->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
if (num_rows(c) < static_cast<unsigned int>(nPoints))
@@ -248,11 +253,16 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
if (firstCall) {
- const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
- basFcts = colFeSpace->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
+ basFcts = colFeSpace->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
}
ElementVector c(nPoints, 0.0);
@@ -282,11 +292,16 @@ namespace AMDiS {
ElementMatrix& mat)
{
if (firstCall) {
- q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(),
- colFeSpace->getBasisFcts(),
- quadrature);
- q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(),
+ colFeSpace->getBasisFcts(),
+ quadrature);
+ q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);
+ firstCall = false;
+ }
}
ElementVector c(1, 0.0);
@@ -318,11 +333,16 @@ namespace AMDiS {
ElementVector& vec)
{
if (firstCall) {
- q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(),
- colFeSpace->getBasisFcts(),
- quadrature);
- q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);
- firstCall = false;
+ #ifdef _OPENMP
+ #pragma omp critical
+ #endif
+ {
+ q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(),
+ colFeSpace->getBasisFcts(),
+ quadrature);
+ q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);
+ firstCall = false;
+ }
}
vector<OperatorTerm*>::iterator termIt;
diff --git a/AMDiS/src/expressions/valueOf.hpp b/AMDiS/src/expressions/valueOf.hpp
index 0306d536..b7d35738 100644
--- a/AMDiS/src/expressions/valueOf.hpp
+++ b/AMDiS/src/expressions/valueOf.hpp
@@ -395,13 +395,13 @@ namespace AMDiS
template<typename Name, template<class> class Matrix, typename T>
typename boost::enable_if<typename traits::is_matrix<Matrix<T> >::type,
expressions::ValueOf<Matrix<DOFVector<T>*>, Name > >::type
- valueOf(Matrix<DOFVector<T>*> &mat)
+ valueOf(Matrix<DOFVector<T>*> const& mat)
{ return expressions::ValueOf<Matrix<DOFVector<T>*>, Name >(mat); }
template<typename Name, template<class> class Vector, typename T>
typename boost::enable_if<typename traits::is_vector<Vector<T> >::type,
expressions::ValueOf<Vector<DOFVector<T>*>, Name > >::type
- valueOf(Vector<DOFVector<T>*> &vector)
+ valueOf(Vector<DOFVector<T>*> const& vector)
{ return expressions::ValueOf<Vector<DOFVector<T>*>, Name >(vector); }
@@ -417,13 +417,13 @@ namespace AMDiS
template<template<class> class Matrix, typename T>
typename boost::enable_if<typename traits::is_matrix<Matrix<T> >::type,
expressions::ValueOf<Matrix<DOFVector<T>*>, _unknown > >::type
- valueOf(Matrix<DOFVector<T>*> &mat)
+ valueOf(Matrix<DOFVector<T>*> const& mat)
{ return expressions::ValueOf<Matrix<DOFVector<T>*>, _unknown >(mat); }
template<template<class...> class Vector, typename T>
typename boost::enable_if<typename traits::is_vector<Vector<T> >::type,
expressions::ValueOf<Vector<DOFVector<T>*>, _unknown > >::type
- valueOf(Vector<DOFVector<T>*> &vector)
+ valueOf(Vector<DOFVector<T>*> const& vector)
{ return expressions::ValueOf<Vector<DOFVector<T>*>, _unknown >(vector); }
--
GitLab