diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 5c7e0796aadf5710277ebd4d2d6c02caa8d1410e..d14a984699e66f5e8cc700eb070f2f57a0f19614 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -24,7 +24,11 @@ namespace AMDiS {
bool optimized,
FirstOrderType type)
: SubAssembler(op, assembler, quad, 1, optimized, type)
- {}
+ {
+ VectorOfFixVecs<DimVec<double> >
+ Lb(assembler->getRowFESpace()->getMesh()->getDim(), 1, NO_INIT);
+ tmpLb.resize(omp_get_overall_max_threads(), Lb);
+ }
FirstOrderAssembler*
FirstOrderAssembler::getSubAssembler(Operator* op,
@@ -95,8 +99,8 @@ namespace AMDiS {
} else {
newAssembler =
(type == GRD_PSI) ?
- dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) :
- dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad));
+ dynamic_cast<FirstOrderAssembler*>(NEW Quad10(op, assembler, quad)) :
+ dynamic_cast<FirstOrderAssembler*>(NEW Quad01(op, assembler, quad));
}
}
@@ -108,23 +112,22 @@ namespace AMDiS {
: FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
{}
-
void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
- double *phival = GET_MEMORY(double, nCol);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
+ std::vector<double> phival(nCol);
- VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
-
- int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
@@ -143,7 +146,6 @@ namespace AMDiS {
}
}
}
- FREE_MEMORY(phival, double, nCol);
}
void Stand10::calculateElementMatrix(const ElInfo *rowElInfo,
@@ -165,13 +167,14 @@ namespace AMDiS {
TEST_EXIT(m)("No subElemCoordsMat!\n");
int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
- VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
- int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->
getLb(smallElInfo, nPoints, Lb);
@@ -202,17 +205,41 @@ namespace AMDiS {
}
}
+ void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+ {
+ DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
+ const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
+ int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
+
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].set(0.0);
+ }
+ for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
+ }
+
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq] *= elInfo->getDet();
+
+ for (int i = 0; i < nRow; i++) {
+ (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
+ (*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
+ }
+ }
+ }
+
Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
- {
- }
+ {}
void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
- const double *phi;
if (firstCall) {
#ifdef _OPENMP
@@ -228,13 +255,14 @@ namespace AMDiS {
}
int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
- VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
- int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
@@ -242,7 +270,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
- phi = phiFast->getPhi(iq);
+ const double *phi = phiFast->getPhi(iq);
grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) {
@@ -252,6 +280,45 @@ namespace AMDiS {
}
}
+ void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+ {
+ VectorOfFixVecs<DimVec<double> > *grdPsi;
+
+ if (firstCall) {
+#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();
+ int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
+
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].set(0.0);
+ }
+ for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
+ }
+
+ for (int iq = 0; iq < nPoints; iq++) {
+
+ Lb[iq] *= elInfo->getDet();
+ grdPsi = psiFast->getGradient(iq);
+
+ for (int i = 0; i < nRow; i++) {
+ (*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
+ }
+ }
+ }
Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
@@ -261,7 +328,6 @@ namespace AMDiS {
void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
- VectorOfFixVecs<DimVec<double> > Lb(dim, 1, NO_INIT);
const int *k;
const double *values;
@@ -281,6 +347,8 @@ namespace AMDiS {
const int **nEntries = q10->getNumberEntries();
int myRank = omp_get_thread_num();
+ // Do not need do resize Lb, because it's size is always at least one.
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
@@ -307,7 +375,6 @@ namespace AMDiS {
: FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
{}
-
void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
@@ -316,8 +383,9 @@ namespace AMDiS {
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
@@ -347,7 +415,7 @@ namespace AMDiS {
const ElInfo *largeElInfo,
ElementMatrix *mat)
{
- FUNCNAME("Stand10::calculateElementMatrix()");
+ FUNCNAME("Stand01::calculateElementMatrix()");
TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");
@@ -360,13 +428,14 @@ namespace AMDiS {
TEST_EXIT(m)("No subElemCoordsMat!\n");
int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
- VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
- int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->
getLb(smallElInfo, nPoints, Lb);
@@ -396,35 +465,9 @@ namespace AMDiS {
}
}
- void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
- {
- DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
- const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
- int nPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
- int myRank = omp_get_thread_num();
-
- for (int iq = 0; iq < nPoints; iq++) {
- Lb[iq].set(0.0);
- }
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
- }
-
- for (int iq = 0; iq < nPoints; iq++) {
- Lb[iq] *= elInfo->getDet();
-
- for (int i = 0; i < nRow; i++) {
- (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
- (*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
- }
- }
- }
-
Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
- {
- }
+ {}
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
@@ -444,8 +487,9 @@ namespace AMDiS {
}
int nPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
+ Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
@@ -467,45 +511,6 @@ namespace AMDiS {
}
}
- void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
- {
- VectorOfFixVecs<DimVec<double> > *grdPsi;
-
- if (firstCall) {
-#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();
- VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
- int myRank = omp_get_thread_num();
-
- for (int iq = 0; iq < nPoints; iq++) {
- Lb[iq].set(0.0);
- }
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
- }
-
- for (int iq = 0; iq < nPoints; iq++) {
-
- Lb[iq] *= elInfo->getDet();
- grdPsi = psiFast->getGradient(iq);
-
- for (int i = 0; i < nRow; i++) {
- (*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
- }
- }
- }
-
Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{
@@ -513,8 +518,6 @@ namespace AMDiS {
void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
- VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
-
const int *l;
const double *values;
@@ -534,6 +537,8 @@ namespace AMDiS {
const int **nEntries = q01->getNumberEntries();
int myRank = omp_get_thread_num();
+ // Do not need to resize Lb, because it's size is always at least one!
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
@@ -557,8 +562,6 @@ namespace AMDiS {
void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
- VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
-
const int *k;
const double *values;
@@ -578,6 +581,8 @@ namespace AMDiS {
const int *nEntries = q1->getNumberEntries();
int myRank = omp_get_thread_num();
+ // Do not need to resize Lb, because it's size is always at least one!
+ VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h
index 16236732272bb31b861e62be725784cefcdacef8..0d8de4b9815dc37122ed45e77c8443a93a37e5a7 100644
--- a/AMDiS/src/FirstOrderAssembler.h
+++ b/AMDiS/src/FirstOrderAssembler.h
@@ -46,6 +46,12 @@ namespace AMDiS {
protected:
+ /** \brief
+ * Thread safe temporary vector of DimMats for calculation in
+ * function calculateElementMatrix().
+ */
+ std::vector<VectorOfFixVecs<DimVec<double> > > tmpLb;
+
/// List of all yet created optimized zero order assemblers for grdPsi.
static std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
@@ -71,7 +77,7 @@ namespace AMDiS {
public:
MEMORY_MANAGED(Stand10);
- /// Constructor.
+ /// Constructor
Stand10(Operator *op, Assembler *assembler, Quadrature *quad);
/// Implements SubAssembler::calculateElementMatrix().
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index 6fcb65d4a1fc2a60b54dbf4548ab61f2d12dfc21..0d0ed8cabf78e4ca4944713fe4f662becf59c85b 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -108,33 +108,23 @@ namespace AMDiS {
this->set(ini);
}
- /** \brief
- * Initialisation for dim.
- */
- inline void init(int dim)
- {
+ /// Initialisation for dim.
+ inline void init(int dim) {
this->resize(calcSize(dim));
}
- /** \brief
- * Initialisation for size
- */
- inline void initSize(int size)
- {
+ /// Initialisation for size
+ inline void initSize(int size) {
this->resize(size);
}
- /** \brief
- * Returns the \ref size_ of the FixVec.
- */
+ /// Returns the \ref size_ of the FixVec.
inline int size() const {
return this->getSize();
}
protected:
- /** \brief
- * Determines needed vector size.
- */
+ /// Determines needed vector size.
static int calcSize(int dim) {
if (dim < 0) {
return Global::getGeo(WORLD);
@@ -172,15 +162,15 @@ namespace AMDiS {
* FixVec's constructors. size_ is the number of contained FixVecs. initType
* must be NO_INIT.
*/
- VectorOfFixVecs(int dim, int size, InitType initType)
- : size_(size)
+ VectorOfFixVecs(int d, int s, InitType initType)
+ : size(s),
+ dim(d)
{
TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
- vec = GET_MEMORY(FixVecType*, size_);
- for (FixVecType** i = &vec[0]; i < &vec[size_]; i++) {
- *i = NEW FixVecType(dim, NO_INIT);
- }
+ vec.resize(size);
+ for (int i = 0; i < size; i++)
+ vec[i] = NEW FixVecType(dim, NO_INIT);
}
/** \brief
@@ -188,15 +178,15 @@ namespace AMDiS {
* FixVec's constructors. size_ is the number of contained FixVecs. initType
* must be VALUE_LIST. ini contains the initialisation values.
*/
- VectorOfFixVecs(int dim, int s, InitType initType, const FixVecType* ini)
- : size_(s)
+ VectorOfFixVecs(int d, int s, InitType initType, const FixVecType* ini)
+ : size(s),
+ dim(d)
{
TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
- vec = GET_MEMORY(FixVecType*, size_);
- for (FixVecType** i = &vec[0]; i < &vec[size_]; i++) {
- *i = NEW FixVecType(ini[i]);
- }
+ vec.resize(size);
+ for (int i = 0; i < size; i++)
+ vec[i] = NEW FixVecType(ini[i]);
}
/** \brief
@@ -204,100 +194,96 @@ namespace AMDiS {
* FixVec's constructors. size_ is the number of contained FixVecs. initType
* must be DEFAULT_VALUE. All entries are set to ini.
*/
- VectorOfFixVecs(int /*dim*/, int s, InitType initType, const FixVecType& ini)
- : size_(s)
+ VectorOfFixVecs(int d, int s, InitType initType, const FixVecType& ini)
+ : size(s),
+ dim(d)
{
- TEST_EXIT_DBG(initType == DEFAULT_VALUE)
- ("wrong initType or wrong initializer\n");
- vec = GET_MEMORY(FixVecType*, size_);
- for (FixVecType** i = &vec[0]; i < &vec[size_]; i++) {
- *i = NEW FixVecType(ini);
- }
+ TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
+
+ vec.resize(size);
+ for (int i = 0; i < size; i++)
+ vec[i] = NEW FixVecType(ini);
}
- /** \brief
- * copy constructor
- */
+ /// Copy constructor
VectorOfFixVecs(const VectorOfFixVecs<FixVecType>& rhs)
{
- size_ = rhs.getSize();
- vec = GET_MEMORY(FixVecType*, size_);
+ size = rhs.getSize();
+ dim = rhs.getDim();
- for (int i = 0; i < size_; i++) {
+ vec.resize(size);
+ for (int i = 0; i < size; i++)
vec[i] = NEW FixVecType(*(rhs.vec[i]));
- }
- };
+ }
- /** \brief
- * destructor
- */
+ /// Destructor
virtual ~VectorOfFixVecs()
{
- for (FixVecType** i = &vec[0]; i<&vec[size_]; i++) {
- DELETE *i;
- }
- FREE_MEMORY(vec, FixVecType*, size_);
- //delete [] vec;
- };
+ for (int i = 0; i < size; i++)
+ DELETE vec[i];
- /** \brief
- * Allows the access like in a normal array via []. Used for const objects.
- */
+ vec.clear();
+ }
+
+ /// Allows the access like in a normal array via []. Used for const objects.
inline const FixVecType& operator[](int index) const
{
- TEST_EXIT_DBG(index >= 0 && index < size_)("invalid index\n");
+ TEST_EXIT_DBG(index >= 0 && index < size)("invalid index\n");
return *(vec[index]);
- };
+ }
- /** \brief
- * Allows the access like in a normal array via [].
- */
+ /// Allows the access like in a normal array via [].
inline FixVecType& operator[](int index)
{
- TEST_EXIT_DBG(index >= 0 && index < size_)("invalid index\n");
+ TEST_EXIT_DBG(index >= 0 && index < size)("invalid index\n");
return *(vec[index]);
- };
+ }
- /** \brief
- * assignment operator
- */
+ /// Assignment operator
VectorOfFixVecs<FixVecType>&
operator=(const VectorOfFixVecs<FixVecType>& rhs)
{
- TEST_EXIT_DBG(size_ == rhs.size_)("vectors of different size\n");
- if (this!=&rhs) {
- FixVecType **i, **j;
- for(i=&vec[0], j=&(rhs.vec[0]); i < &vec[size_]; i++, j++) {
- **i = **j;
- }
+ TEST_EXIT_DBG(size == rhs.size)("vectors of different size\n");
+ if (this != &rhs) {
+ for (int i = 0; i < size; i++)
+ *(vec[i]) = *(rhs.vec[i]);
}
return *this;
- };
+ }
- /** \brief
- * Returns the \ref size of this VectorOfFixVecs
- */
+ /// Resize the vector
+ inline void resize(int newsize)
+ {
+ vec.resize(newsize);
+ for (int i = size; i < newsize; i++)
+ vec[i] = NEW FixVecType(dim, NO_INIT);
+ size = newsize;
+ }
+
+ /// Returns the \ref size of this VectorOfFixVecs
inline int getSize() const {
- return size_;
- };
+ return size;
+ }
- /** \brief
- * Returns the size of the contained FixVecs
- */
+ /// Returns \ref dim
+ inline int getDim() const {
+ return dim;
+ }
+
+ /// Returns the size of the contained FixVecs
inline int getSizeOfFixVec() const {
return vec[0]->getSize();
}
protected:
- /** \brief
- * number of contained FixVecs
- */
- int size_;
+ /// number of contained FixVecs
+ int size;
- /** \brief
- * array of pointers to FixVecs
- */
- FixVecType **vec;
+ /// dimension of world
+ int dim;
+
+ /// array of pointers to FixVecs
+ std::vector<FixVecType*> vec;
};
// ===========================================================================
@@ -325,70 +311,52 @@ namespace AMDiS {
{
TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
vec = GET_MEMORY(VectorOfFixVecs<FixVecType>*, rows);
- for(VectorOfFixVecs<FixVecType>** i=&vec[0]; i<&vec[rows]; i++) {
+ for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++) {
*i = NEW VectorOfFixVecs<FixVecType>(dim, columns, NO_INIT);
}
- };
-
+ }
- /** \brief
- * destructor
- */
+ /// destructor
virtual ~MatrixOfFixVecs()
{
- for(VectorOfFixVecs<FixVecType>** i=&vec[0]; i<&vec[rows]; i++) {
+ for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++) {
DELETE *i;
}
FREE_MEMORY(vec, VectorOfFixVecs<FixVecType>*, rows);
- };
+ }
- /** \brief
- * assignment operator
- */
+ /// assignment operator
inline VectorOfFixVecs<FixVecType >& operator[](int index)
{
TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
return *(vec[index]);
- };
+ }
- /** \brief
- * assignment operator
- */
+ /// assignment operator
inline const VectorOfFixVecs<FixVecType >& operator[](int index) const
{
TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
return *(vec[index]);
- };
+ }
- /** \brief
- * Returns \ref rows
- */
+ /// Returns \ref rows
inline int getNumberOfRows() const {
return rows;
- };
+ }
- /** \brief
- * Returns \ref columns
- */
+ /// Returns \ref columns
inline int getNumberOfColumns() const {
return columns;
- };
+ }
private:
- /** \brief
- * number of rows of the matrix
- */
-
+ /// number of rows of the matrix
int rows;
- /** \brief
- * number of columns of the matrix
- */
+ /// number of columns of the matrix
int columns;
- /** \brief
- * array of pointers to VectorOfFixVecs
- */
+ /// array of pointers to VectorOfFixVecs
VectorOfFixVecs<FixVecType> **vec;
};
@@ -408,32 +376,24 @@ namespace AMDiS {
public:
MEMORY_MANAGED(DimVec<T>);
- DimVec() {};
+ DimVec() {}
- /** \brief
- * Calls the corresponding constructor of FixVec
- */
- DimVec(int dim, InitType initType=NO_INIT)
+ /// Calls the corresponding constructor of FixVec
+ DimVec(int dim, InitType initType = NO_INIT)
: FixVec<T,PARTS>(dim, initType)
{}
- /** \brief
- * Calls the corresponding constructor of FixVec
- */
+ /// Calls the corresponding constructor of FixVec
DimVec(int dim, InitType initType, T* ini)
: FixVec<T,PARTS>(dim, initType, ini)
{}
- /** \brief
- * Calls the corresponding constructor of FixVec
- */
+ /// Calls the corresponding constructor of FixVec
DimVec(int dim, InitType initType, const T& ini)
: FixVec<T,PARTS>(dim, initType, ini)
{}
- /** \brief
- * Multplication of a matrix with a vector.
- */
+ /// Multplication of a matrix with a vector.
void multMatrixVec(DimMat<T> &m, DimVec<T> &v);
};
@@ -564,22 +524,18 @@ namespace AMDiS {
public:
MEMORY_MANAGED(WorldMatrix<T>);
- /** \brief
- * Calls the corresponding constructor of FixVec
- */
+ /// Calls the corresponding constructor of FixVec
WorldMatrix()
: Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
- {};
+ {}
- /** \brief
- * Calls the corresponding constructor of FixVec
- */
+ /// Calls the corresponding constructor of FixVec
WorldMatrix(InitType initType, T* ini)
: Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
{
TEST_EXIT_DBG(initType == VALUE_LIST)("???\n");
setValues(ini);
- };
+ }
/** \brief
* Calls the corresponding constructor of FixVec and sets all matrix entries
@@ -590,21 +546,15 @@ namespace AMDiS {
{
TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
this->set(ini);
- };
+ }
- /** \brief
- * Returns true if the matrix is a diagonal matrix, returns false otherwise.
- */
+ /// Returns true if the matrix is a diagonal matrix, returns false otherwise.
bool isDiagMatrix() const;
- /** \brief
- * Returns true if the matrix is symmetric, returns false otherwise.
- */
+ /// Returns true if the matrix is symmetric, returns false otherwise.
bool isSymmetric() const;
- /** \brief
- * Sets the diagonal entries to the given value.
- */
+ /// Sets the diagonal entries to the given value.
void setDiag(T value);
/** \brief
@@ -625,40 +575,32 @@ namespace AMDiS {
// ===========================================================================
- /** \brief
- * returns the euclidian distance of a and b
- */
+ /// returns the euclidian distance of a and b
template<typename T, GeoIndex d>
double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b);
- /** \brief
- * FixVec operator for stream output
- */
+ /// FixVec operator for stream output
template<typename T, GeoIndex d>
std::ostream& operator <<(std::ostream& out, const FixVec<T,d>& fixvec)
{
- for(int i=0; i < fixvec.getSize()-1; i++) {
+ for (int i = 0; i < fixvec.getSize()-1; i++) {
out << fixvec[i] << " ";
}
out << fixvec[fixvec.getSize()-1];
return out;
}
- /** \brief
- * creates and inits a VectorOfFixVecs<DimVec<double> >
- */
+ /// creates and inits a VectorOfFixVecs<DimVec<double> >
VectorOfFixVecs<DimVec<double> > *createAndInit(int dim, int size, ...);
- /** \brief
- * creates and inits and double array
- */
+ /// creates and inits and double array
double *createAndInitArray(int size, ...);
inline WorldVector<double> operator*(const WorldVector<double>& v, double d) {
WorldVector<double> result = v;
result *= d;
return result;
- };
+ }
inline WorldVector<double> operator+(const WorldVector<double>& v1,
const WorldVector<double>& v2)
@@ -666,7 +608,7 @@ namespace AMDiS {
WorldVector<double> result = v1;
result += v2;
return result;
- };
+ }
inline WorldVector<int> operator+(const WorldVector<int>& v1,
const WorldVector<int>& v2)
@@ -674,7 +616,7 @@ namespace AMDiS {
WorldVector<int> result = v1;
result += v2;
return result;
- };
+ }
inline WorldVector<double> operator-(const WorldVector<double>& v1,
const WorldVector<double>& v2)
@@ -682,28 +624,29 @@ namespace AMDiS {
WorldVector<double> result = v1;
result -= v2;
return result;
- };
+ }
inline bool operator<(const WorldVector<double>& v1,
const WorldVector<double>& v2)
{
- int i, dow = Global::getGeo(WORLD);
- for(i = 0; i < dow; i++) {
- if(abs(v1[i] - v2[i]) < DBL_TOL) continue;
+ int dow = Global::getGeo(WORLD);
+ for (int i = 0; i < dow; i++) {
+ if (abs(v1[i] - v2[i]) < DBL_TOL)
+ continue;
return v1[i] < v2[i];
}
return false;
- };
+ }
inline bool operator==(const WorldVector<double>& v1,
const WorldVector<double>& v2)
{
- int i, dow = Global::getGeo(WORLD);
- for(i = 0; i < dow; i++) {
- if(abs(v1[i] - v2[i]) > DBL_TOL) return false;
+ int dow = Global::getGeo(WORLD);
+ for (int i = 0; i < dow; i++) {
+ if (abs(v1[i] - v2[i]) > DBL_TOL) return false;
}
return true;
- };
+ }
template<typename T>
const WorldMatrix<T>& operator*=(WorldMatrix<T>& m, T scal);
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index c4e6d759e37834b306f9968ccf983e60f75a8023..35fb49f653915aea763d5372bad63ee58ca944a7 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -1328,7 +1328,7 @@ namespace AMDiS {
}
}
- void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
if (f) {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq]);
@@ -1359,7 +1359,7 @@ namespace AMDiS {
}
- void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
}
@@ -1378,7 +1378,7 @@ namespace AMDiS {
}
}
- void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
}
@@ -1397,7 +1397,7 @@ namespace AMDiS {
}
}
- void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
}
@@ -1419,7 +1419,8 @@ namespace AMDiS {
}
- void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
}
@@ -1439,7 +1440,8 @@ namespace AMDiS {
}
- void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
}
@@ -1458,7 +1460,8 @@ namespace AMDiS {
}
}
- void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
}
@@ -1494,14 +1497,16 @@ namespace AMDiS {
}
}
- void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const
+ void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
{
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
}
}
- void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
}
@@ -1519,7 +1524,8 @@ namespace AMDiS {
fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
}
- void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const
+ void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
@@ -1536,7 +1542,8 @@ namespace AMDiS {
result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
}
- void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
}
@@ -1553,7 +1560,7 @@ namespace AMDiS {
fac * (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]) * uhAtQP[iq];
}
- void FctGradient_ZOT::getC(const ElInfo *, int nPoints, double *C) const {
+ void FctGradient_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(gradAtQPs[iq]);
}
@@ -1569,7 +1576,8 @@ namespace AMDiS {
result[iq] += fac * (*f)(gradAtQPs[iq]) * uhAtQP[iq];
}
- void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, double *C) const {
+ void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints,
+ std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*g)(coordsAtQPs[iq]);
}
@@ -2257,7 +2265,8 @@ namespace AMDiS {
}
}
- void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const
+ void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
{
int size = static_cast<int>(vecs.size());
std::vector<double> arg(size);
@@ -2310,7 +2319,8 @@ namespace AMDiS {
}
- void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints, double *C) const
+ void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
{
int size = static_cast<int>(vecs.size());
std::vector<WorldVector<double>*> arg(size);
@@ -2341,7 +2351,8 @@ namespace AMDiS {
}
}
- void VecDivergence_ZOT::getC(const ElInfo *, int nPoints, double *C) const
+ void VecDivergence_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
{
int size = static_cast<int>(vecs.size());
@@ -2618,7 +2629,8 @@ namespace AMDiS {
}
}
- void General_ZOT::getC(const ElInfo *, int nPoints, double *C) const
+ void General_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
@@ -2865,7 +2877,8 @@ namespace AMDiS {
}
}
- void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, double *C) const
+ void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
@@ -2920,7 +2933,7 @@ namespace AMDiS {
}
void VecAndVecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
- double *C) const
+ std::vector<double> &C) const
{
int size = static_cast<int>(vecs.size());
std::vector<WorldVector<double>*> arg(size);
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 43d31440c2b857bdaa72dffbb0cf8c6a3f186283..236872019c820362a275270cce1b754bae145090 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -196,7 +196,7 @@ namespace AMDiS {
double factor)
{
int dim = Lb.getSize() - 1;
- static const int dimOfWorld = Global::getGeo(WORLD);
+ const int dimOfWorld = Global::getGeo(WORLD);
for (int i = 0; i <= dim; i++) {
double val = 0.0;
@@ -2041,7 +2041,7 @@ namespace AMDiS {
*/
virtual void getC(const ElInfo *elInfo,
int nPoints,
- double *result) const = 0;
+ std::vector<double> &C) const = 0;
};
@@ -2057,10 +2057,8 @@ namespace AMDiS {
/// Constructor.
Simple_ZOT(double f = 1.0) : ZeroOrderTerm(0), factor(f) {}
- /** \brief
- * Implements ZeroOrderTerm::getC().
- */
- inline void getC(const ElInfo *, int nPoints, double *C) const
+ /// Implements ZeroOrderTerm::getC().
+ inline void getC(const ElInfo *, int nPoints, std::vector<double> &C) const
{
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += factor;
@@ -2113,10 +2111,8 @@ namespace AMDiS {
SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements ZeroOrderTerm::getC().
- */
- void getC(const ElInfo *, int nPoints, double *C) const;
+ /// Implements ZeroOrderTerm::getC().
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2169,7 +2165,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2226,7 +2222,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
@@ -2276,7 +2272,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2332,7 +2328,7 @@ namespace AMDiS {
/** \brief
* Implements SecondOrderTerm::getC().
*/
- void getC(const ElInfo *elInfo, int nPoints, double *C) const;
+ void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2383,7 +2379,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2440,7 +2436,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2496,7 +2492,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2555,7 +2551,7 @@ namespace AMDiS {
/** \brief
* Implements SecondOrderTerm::getC().
*/
- void getC(const ElInfo *elInfo, int nPoints, double *C) const;
+ void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2604,7 +2600,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2725,7 +2721,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2898,7 +2894,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2955,7 +2951,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -3021,7 +3017,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -3066,7 +3062,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -3114,7 +3110,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -3157,7 +3153,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -3213,7 +3209,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -3263,7 +3259,7 @@ namespace AMDiS {
/** \brief
* Implements ZeroOrderTerm::getC().
*/
- void getC(const ElInfo *, int nPoints, double *C) const;
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -3627,7 +3623,7 @@ namespace AMDiS {
* Calls getC() for each term in \ref zeroOrder
* and adds the results to c.
*/
- void getC(const ElInfo *elInfo, int nPoints, double *c) const
+ void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c) const
{
int myRank = omp_get_thread_num();
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index b65fd0575d181173e31b093c8b86e3cd72d13470..a510dddf43f9f65df2fbf22a1c87bb95cc944ac8 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -725,16 +725,11 @@ namespace AMDiS {
result, -factor);
}
} else {
- double *fx = GET_MEMORY(double, nPoints);
- for (int iq = 0; iq < nPoints; iq++) {
- fx[iq] = 0.0;
- }
+ std::vector<double> fx(nPoints, 0.0);
(*it)->getC(elInfo, nPoints, fx);
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
result[iq] -= factor * fx[iq];
- }
- FREE_MEMORY(fx, double, nPoints);
}
}
}
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index 436d82b0c72320765de550f02583d1e4dca9bad3..cb6dd5f38b36b6103130111295c5d0c8eac233f4 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -292,17 +292,14 @@ namespace AMDiS {
NULL,
NULL);
- double *f = GET_MEMORY(double, numPoints);
- for (iq = 0; iq < numPoints; iq++) {
- f[iq] = 0.0;
- }
+ std::vector<double> f(numPoints, 0.0);
if(robinOperators) {
(*robinOperators)[face]->evalZeroOrder(numPoints,
uhAtQp,
NULL,
NULL,
- f,
+ &f[0],
1.0);
}
@@ -322,7 +319,7 @@ namespace AMDiS {
}
- if(neumannOperators)
+ if (neumannOperators)
(*neumannOperators)[face]->getC(elInfo, numPoints, f);
for (val = iq = 0; iq < numPoints; iq++) {
@@ -332,8 +329,6 @@ namespace AMDiS {
DELETE [] grdUh;
DELETE [] A_grdUh;
-
- FREE_MEMORY(f, double, numPoints);
}
}
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 8a023e9a6b5b221d53542e5b858fbf7104f14485..727edcce5d136f1d7633404ba5233c31e660cc2e 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -114,7 +114,7 @@ namespace AMDiS {
DimMat<double> &tmpMat = *LALt[0];
tmpMat.set(0.0);
- for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
}
@@ -122,6 +122,8 @@ namespace AMDiS {
nEntries = q11->getNumberEntries();
if (symmetric) {
+ TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
+
for (int i = 0; i < nRow; i++) {
k = q11->getKVec(i, i);
l = q11->getLVec(i, i);
@@ -210,6 +212,8 @@ namespace AMDiS {
VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
if (symmetric) {
+ TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
+
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
@@ -272,6 +276,8 @@ namespace AMDiS {
}
if (symmetric) {
+ TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
+
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
diff --git a/AMDiS/src/SecondOrderAssembler.h b/AMDiS/src/SecondOrderAssembler.h
index e82cc09fa1d5ac657657db9e5598988cc4f70e10..0e2f113bb21af189b55bde3012f3f574916e3bb6 100644
--- a/AMDiS/src/SecondOrderAssembler.h
+++ b/AMDiS/src/SecondOrderAssembler.h
@@ -157,10 +157,7 @@ namespace AMDiS {
}
protected:
- /** \brief
- * Integral of the product of the derivative of psi and the derivative
- * of phi.
- */
+ /// Integral of the product of the derivative of psi and the derivative of phi.
const Q11PsiPhi *q11;
/** \brief
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index f8dadcc0790f703203611a90cce38e6a478c13e1..f1c33085295b9b7f5f78717f3c6549f310690a82 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -52,14 +52,21 @@ namespace AMDiS {
break;
}
- // check if all terms are symmetric
- symmetric = true;
- for (int i = 0; i < static_cast<int>(terms[0].size()); i++) {
- if (!(terms[0][i])->isSymmetric()) {
- symmetric = false;
- break;
- }
- }
+ // If the number of basis functions in row and col fe space are equal,
+ // the element matrix may be symmetric if all operator terms are symmetric.
+ // If the row and col fe space are different, the element matrix is neither
+ // quadratic, and therefore cannot be symmetric.
+ if (nRow == nCol) {
+ symmetric = true;
+ for (int i = 0; i < static_cast<int>(terms[0].size()); i++) {
+ if (!(terms[0][i])->isSymmetric()) {
+ symmetric = false;
+ break;
+ }
+ }
+ } else {
+ symmetric = false;
+ }
dim = assembler->rowFESpace->getMesh()->getDim();
}
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 8fdb4605e68ea79872beab3abaa11dc147c5bdfd..cdd5f726f5d6e44d26f0ab624ba4c96ab3d5d842 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -89,13 +89,9 @@ namespace AMDiS {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
- double *phival = GET_MEMORY(double, nCol);
int nPoints = quadrature->getNumPoints();
- double *c = GET_MEMORY(double, nPoints);
-
- for (int iq = 0; iq < nPoints; iq++) {
- c[iq] = 0.0;
- }
+ std::vector<double> c(nPoints, 0.0);
+ std::vector<double> phival(nCol);
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
@@ -104,6 +100,8 @@ namespace AMDiS {
}
if (symmetric) {
+ TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
+
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
@@ -139,9 +137,6 @@ namespace AMDiS {
}
}
}
-
- FREE_MEMORY(phival, double, nCol);
- FREE_MEMORY(c, double, nPoints);
}
void StandardZOA::calculateElementMatrix(const ElInfo *rowElInfo,
@@ -162,10 +157,7 @@ namespace AMDiS {
int myRank = omp_get_thread_num();
int nPoints = quadrature->getNumPoints();
- double *c = GET_MEMORY(double, nPoints);
- for (int iq = 0; iq < nPoints; iq++) {
- c[iq] = 0.0;
- }
+ std::vector<double> c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
@@ -191,20 +183,12 @@ namespace AMDiS {
}
}
-
-
- FREE_MEMORY(c, double, nPoints);
}
void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
int nPoints = quadrature->getNumPoints();
-
- double *c = GET_MEMORY(double, nPoints);
-
- for (int iq = 0; iq < nPoints; iq++) {
- c[iq] = 0.0;
- }
+ std::vector<double> c(nPoints, 0.0);
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
@@ -221,22 +205,11 @@ namespace AMDiS {
(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi;
}
}
-
- FREE_MEMORY(c, double, nPoints);
}
FastQuadZOA::FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
- {
- cPtrs.resize(omp_get_overall_max_threads());
- }
-
- FastQuadZOA::~FastQuadZOA()
- {
- for (int i = 0; i < omp_get_overall_max_threads(); i++) {
- FREE_MEMORY(cPtrs[i], double, quadrature->getNumPoints());
- }
- }
+ {}
void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
@@ -248,7 +221,6 @@ namespace AMDiS {
#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();
@@ -257,17 +229,15 @@ namespace AMDiS {
}
}
- double *c = cPtrs[myRank];
- for (int iq = 0; iq < nPoints; iq++) {
- c[iq] = 0.0;
- }
-
+ std::vector<double> c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
if (symmetric) {
+ TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
+
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
@@ -316,7 +286,6 @@ namespace AMDiS {
#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();
@@ -325,11 +294,7 @@ namespace AMDiS {
}
}
- double *c = cPtrs[myRank];
- for (int iq = 0; iq < nPoints; iq++) {
- c[iq] = 0.0;
- }
-
+ std::vector<double> c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(rowElInfo, nPoints, c);
@@ -368,7 +333,6 @@ namespace AMDiS {
#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();
@@ -377,11 +341,7 @@ namespace AMDiS {
}
}
- double *c = cPtrs[myRank];
- for (int iq = 0; iq < nPoints; iq++) {
- c[iq] = 0.0;
- }
-
+ std::vector<double> c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
@@ -418,30 +378,31 @@ namespace AMDiS {
}
}
- double c = 0.0;
+ std::vector<double> c(1, 0.0);
int myRank = omp_get_thread_num();
int size = static_cast<int>(terms[myRank].size());
for (int i = 0; i < size; i++) {
- (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, &c);
+ (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, c);
}
- c *= elInfo->getDet();
+ c[0] *= elInfo->getDet();
if (symmetric) {
for (int i = 0; i < nRow; i++) {
- (*mat)[i][i] += c * q00->getValue(i,i);
+ (*mat)[i][i] += c[0] * q00->getValue(i,i);
for (int j = i + 1; j < nCol; j++) {
- double val = c * q00->getValue(i, j);
+ double val = c[0] * q00->getValue(i, j);
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
}
} else {
- for (int i = 0; i < nRow; i++)
+ for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
- (*mat)[i][j] += c * q00->getValue(i, j);
+ (*mat)[i][j] += c[0] * q00->getValue(i, j);
}
+ }
}
}
@@ -464,15 +425,15 @@ namespace AMDiS {
std::vector<OperatorTerm*>::iterator termIt;
int myRank = omp_get_thread_num();
- double c = 0.0;
+ std::vector<double> c(1, 0.0);
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
- (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, &c);
+ (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, c);
}
- c *= elInfo->getDet();
+ c[0] *= elInfo->getDet();
for (int i = 0; i < nRow; i++)
- (*vec)[i] += c * q0->getValue(i);
+ (*vec)[i] += c[0] * q0->getValue(i);
}
}
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index 4c9f0350b4e478a15f77b75844aa4cfa48943691..2d89690ef37d3c64d625b15cc6260cc6c5670c3b 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -94,8 +94,6 @@ namespace AMDiS {
/// Constructor.
FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad);
- ~FastQuadZOA();
-
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
@@ -107,9 +105,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *elInfo, ElementVector *vec);
-
- protected:
- std::vector<double*> cPtrs;
};