Commit 33516d51 authored by Praetorius, Simon's avatar Praetorius, Simon

some openMP cleanup

parent 9284b5b6
......@@ -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.
......
......@@ -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];
}
......
......@@ -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");
}
......
......@@ -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();
......
......@@ -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");
......
......@@ -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)
......
......@@ -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,
......
This diff is collapsed.
......@@ -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;
/** \} */
};
......
......@@ -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++)
......
......@@ -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;
......
......@@ -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;
......
......@@ -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,