diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 044632a02c2052cab572c6846197f3b652e4eae9..7686e186e1237207efc3aba1987f0a38570c6948 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -122,7 +122,7 @@ namespace AMDiS {
{
Quadrature *localQuad = quad ? quad : quadrature;
- const int numPoints = localQuad->getNumPoints();
+ const int nPoints = localQuad->getNumPoints();
// already calculated for this element ?
if (coordsValid) {
@@ -130,19 +130,19 @@ namespace AMDiS {
}
if (coordsAtQPs) {
- if (coordsNumAllocated != numPoints) {
+ if (coordsNumAllocated != nPoints) {
DELETE [] coordsAtQPs;
- coordsAtQPs = NEW WorldVector<double>[numPoints];
- coordsNumAllocated = numPoints;
+ coordsAtQPs = NEW WorldVector<double>[nPoints];
+ coordsNumAllocated = nPoints;
}
} else {
- coordsAtQPs = NEW WorldVector<double>[numPoints];
- coordsNumAllocated = numPoints;
+ coordsAtQPs = NEW WorldVector<double>[nPoints];
+ coordsNumAllocated = nPoints;
}
// set new values
WorldVector<double>* k = &(coordsAtQPs[0]);
- for (int l = 0; k < &(coordsAtQPs[numPoints]); ++k, ++l) {
+ for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l) {
elInfo->coordToWorld(localQuad->getLambda(l), k);
}
@@ -223,7 +223,7 @@ namespace AMDiS {
Quadrature *localQuad = quad ? quad : quadrature;
- if(!gradientsAtQPs[vec]) {
+ if (!gradientsAtQPs[vec]) {
gradientsAtQPs[vec] = new GradientsAtQPs;
}
gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());
@@ -249,7 +249,7 @@ namespace AMDiS {
const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();
if (opt && !quad && sameFESpaces) {
- if(psiFast->getBasisFunctions() == basFcts) {
+ if (psiFast->getBasisFunctions() == basFcts) {
vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
} else {
vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
@@ -500,21 +500,21 @@ namespace AMDiS {
double psival;
double *phival = GET_MEMORY(double, nCol);
- int numPoints = quadrature->getNumPoints();
- double *c = GET_MEMORY(double, numPoints);
+ int nPoints = quadrature->getNumPoints();
+ double *c = GET_MEMORY(double, nPoints);
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
- (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
+ (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
if (symmetric) {
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
@@ -533,7 +533,7 @@ namespace AMDiS {
}
}
} else { // non symmetric assembling
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
@@ -551,26 +551,26 @@ namespace AMDiS {
}
FREE_MEMORY(phival, double, nCol);
- FREE_MEMORY(c, double, numPoints);
+ FREE_MEMORY(c, double, nPoints);
}
void Stand0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
- int numPoints = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
- double *c = GET_MEMORY(double, numPoints);
+ double *c = GET_MEMORY(double, nPoints);
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
- (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
+ (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
@@ -580,19 +580,29 @@ namespace AMDiS {
}
}
- FREE_MEMORY(c, double, numPoints);
+ FREE_MEMORY(c, double, nPoints);
}
Quad0::Quad0(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
{
+ cPtrs.resize(omp_get_max_threads());
+ }
+
+ Quad0::~Quad0()
+ {
+ for (int i = 0; i < omp_get_max_threads(); i++) {
+ FREE_MEMORY(cPtrs[i], double, quadrature->getNumPoints());
+ }
}
void Quad0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
- const double *psi, *phi;
+ int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
if (firstCall) {
+ cPtrs[myRank] = GET_MEMORY(double, nPoints);
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -600,25 +610,22 @@ namespace AMDiS {
firstCall = false;
}
- int numPoints = quadrature->getNumPoints();
- double *c = GET_MEMORY(double, numPoints);
-
- for (int iq = 0; iq < numPoints; iq++) {
+ double *c = cPtrs[myRank];
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
- int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
- (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
+ (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
if (symmetric) {
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
- psi = psiFast->getPhi(iq);
- phi = phiFast->getPhi(iq);
+ const double *psi = psiFast->getPhi(iq);
+ const double *phi = phiFast->getPhi(iq);
for (int i = 0; i < nRow; i++) {
(*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
for (int j = i + 1; j < nCol; j++) {
@@ -629,11 +636,11 @@ namespace AMDiS {
}
}
} else { /* non symmetric assembling */
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
- psi = psiFast->getPhi(iq);
- phi = phiFast->getPhi(iq);
+ const double *psi = psiFast->getPhi(iq);
+ const double *phi = phiFast->getPhi(iq);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
@@ -641,7 +648,6 @@ namespace AMDiS {
}
}
}
- FREE_MEMORY(c, double, numPoints);
}
void Quad0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
@@ -654,20 +660,20 @@ namespace AMDiS {
firstCall = false;
}
- int numPoints = quadrature->getNumPoints();
- double *c = GET_MEMORY(double, numPoints);
+ int nPoints = quadrature->getNumPoints();
+ double *c = GET_MEMORY(double, nPoints);
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
int myRank = omp_get_thread_num();
::std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
- (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
+ (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
const double *psi = psiFast->getPhi(iq);
@@ -675,7 +681,7 @@ namespace AMDiS {
(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
}
}
- FREE_MEMORY(c, double, numPoints);
+ FREE_MEMORY(c, double, nPoints);
}
Pre0::Pre0(Operator *op, Assembler *assembler, Quadrature *quad)
@@ -759,19 +765,19 @@ namespace AMDiS {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
- int numPoints = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
- for (int iq = 0; iq < numPoints; iq++) {
+ 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, numPoints, Lb);
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nCol; i++) {
@@ -808,19 +814,19 @@ namespace AMDiS {
firstCall = false;
}
- int numPoints = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
- for (int iq = 0; iq < numPoints; iq++) {
+ 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, numPoints, Lb);
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
phi = phiFast->getPhi(iq);
@@ -891,18 +897,18 @@ namespace AMDiS {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
- int numPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
+ int nPoints = quadrature->getNumPoints();
+ VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
int myRank = omp_get_thread_num();
- for (int iq = 0; iq < numPoints; iq++) {
+ 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, numPoints, Lb);
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nCol; i++) {
@@ -921,18 +927,18 @@ namespace AMDiS {
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
- int numPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
+ int nPoints = quadrature->getNumPoints();
+ VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
- for (int iq = 0; iq < numPoints; iq++) {
+ 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, numPoints, Lb);
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
@@ -959,18 +965,18 @@ namespace AMDiS {
firstCall = false;
}
- int numPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
+ int nPoints = quadrature->getNumPoints();
+ VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
- for (int iq = 0; iq < numPoints; iq++) {
+ 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, numPoints, Lb);
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
const double *psi = psiFast->getPhi(iq);
@@ -995,18 +1001,18 @@ namespace AMDiS {
firstCall = false;
}
- int numPoints = quadrature->getNumPoints();
- VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
+ int nPoints = quadrature->getNumPoints();
+ VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
- for (int iq = 0; iq < numPoints; iq++) {
+ 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, numPoints, Lb);
+ (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
@@ -1179,14 +1185,32 @@ namespace AMDiS {
Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, true)
- {}
+ {
+ tmpLALt.resize(omp_get_max_threads());
+ }
+
+ Quad2::~Quad2()
+ {
+ int nPoints = quadrature->getNumPoints();
+ for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
+ for (int j = 0; j < nPoints; j++) {
+ DELETE tmpLALt[i][j];
+ }
+ DELETE [] tmpLALt[i];
+ }
+ }
void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
- double val;
- VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;
+ int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
if (firstCall) {
+ tmpLALt[myRank] = NEW DimMat<double>*[nPoints];
+ for (int j = 0; j < nPoints; j++) {
+ tmpLALt[myRank][j] = NEW DimMat<double>(dim, NO_INIT);
+ }
+
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -1194,20 +1218,18 @@ namespace AMDiS {
firstCall = false;
}
- int nPoints = quadrature->getNumPoints();
- DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
- int myRank = omp_get_thread_num();
+ DimMat<double> **LALt = tmpLALt[myRank];
- for (int i = 0; i < nPoints;i++) {
- LALt[i] = NEW DimMat<double>(dim, NO_INIT);
- }
- for (int iq = 0; iq < nPoints; iq++) {
- LALt[iq]->set(0.0);
+ for (int i = 0; i < nPoints; i++) {
+ LALt[i]->set(0.0);
}
+
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
}
+ VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;
+
if (symmetric) {
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
@@ -1220,14 +1242,13 @@ namespace AMDiS {
((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));
for (int j = i + 1; j < nCol; j++) {
- val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
+ double val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
}
}
- }
- else { /* non symmetric assembling */
+ } else { /* non symmetric assembling */
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
@@ -1242,11 +1263,6 @@ namespace AMDiS {
}
}
}
-
- for (int i = 0;i < nPoints; i++)
- DELETE LALt[i];
-
- DELETE [] LALt;
}
Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad)
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index df347ea93d915840d4e146941aac179c31081c76..e782eb0435e1039f29fa5a348ce3a6e3aea7c7cd 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -117,7 +117,7 @@ namespace AMDiS {
/** \brief
* Returns \ref terms
*/
- inline ::std::vector<OperatorTerm*> *getTerms() {
+ inline std::vector<OperatorTerm*> *getTerms() {
return &terms[omp_get_thread_num()];
};
@@ -227,12 +227,12 @@ namespace AMDiS {
/** \brief
* Used for \ref getVectorAtQPs().
*/
- ::std::map<const DOFVectorBase<double>*, ValuesAtQPs*> valuesAtQPs;
+ std::map<const DOFVectorBase<double>*, ValuesAtQPs*> valuesAtQPs;
/** \brief
* Used for \ref getGradientsAtQPs().
*/
- ::std::map<const DOFVectorBase<double>*, GradientsAtQPs*> gradientsAtQPs;
+ std::map<const DOFVectorBase<double>*, GradientsAtQPs*> gradientsAtQPs;
/** \brief
* Set and updated by \ref initElement() for each ElInfo.
@@ -278,7 +278,7 @@ namespace AMDiS {
/** \brief
* List of all terms with a contribution to this SubAssembler
*/
- ::std::vector< ::std::vector<OperatorTerm*> > terms;
+ std::vector< std::vector<OperatorTerm*> > terms;
/** \brief
*
@@ -339,12 +339,12 @@ namespace AMDiS {
/** \brief
* List of all yet created optimized SubAssembler objects.
*/
- static ::std::vector<SubAssembler*> optimizedSubAssemblers;
+ static std::vector<SubAssembler*> optimizedSubAssemblers;
/** \brief
* List of all yet created standard SubAssembler objects.
*/
- static ::std::vector<SubAssembler*> standardSubAssemblers;
+ static std::vector<SubAssembler*> standardSubAssemblers;
};
// ============================================================================
@@ -398,6 +398,8 @@ namespace AMDiS {
*/
Quad0(Operator *op, Assembler *assembler, Quadrature *quad);
+ ~Quad0();
+
/** \brief
* Implements SubAssembler::calculateElementMatrix().
*/
@@ -407,6 +409,9 @@ namespace AMDiS {
* Implements SubAssembler::calculateElementVector().
*/
void calculateElementVector(const ElInfo *elInfo, ElementVector *vec);
+
+ protected:
+ std::vector<double*> cPtrs;
};
@@ -502,22 +507,22 @@ namespace AMDiS {
/** \brief
* List of all yet created optimized zero order assemblers for grdPsi.
*/
- static ::std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
+ static std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
/** \brief
* List of all yet created standard zero order assemblers for grdPsi.
*/
- static ::std::vector<SubAssembler*> standardSubAssemblersGrdPsi;
+ static std::vector<SubAssembler*> standardSubAssemblersGrdPsi;
/** \brief
* List of all yet created optimized zero order assemblers for grdPhi.
*/
- static ::std::vector<SubAssembler*> optimizedSubAssemblersGrdPhi;
+ static std::vector<SubAssembler*> optimizedSubAssemblersGrdPhi;
/** \brief
* List of all yet created standard zero order assemblers for grdPhi.
*/
- static ::std::vector<SubAssembler*> standardSubAssemblersGrdPhi;
+ static std::vector<SubAssembler*> standardSubAssemblersGrdPhi;
};
// ============================================================================
@@ -771,12 +776,12 @@ namespace AMDiS {
/** \brief
* List of all yet created optimized second order assemblers.
*/
- static ::std::vector<SubAssembler*> optimizedSubAssemblers;
+ static std::vector<SubAssembler*> optimizedSubAssemblers;
/** \brief
* List of all yet created standard second order assemblers.
*/
- static ::std::vector<SubAssembler*> standardSubAssemblers;
+ static std::vector<SubAssembler*> standardSubAssemblers;
};
// ============================================================================
@@ -829,6 +834,8 @@ namespace AMDiS {
*/
Quad2(Operator *op, Assembler *assembler, Quadrature *quad);
+ ~Quad2();
+
/** \brief
* Implements SubAssembler::calculateElementMatrix().
*/
@@ -840,6 +847,12 @@ namespace AMDiS {
void calculateElementVector(const ElInfo *, ElementVector */*vec*/) {
ERROR_EXIT("should not be called\n");
};
+
+ protected:
+ /** \brief
+ * Thread safe temporary vector of DimMats for calculation in calculateElementMatrix().
+ */
+ std::vector< DimMat<double>** > tmpLALt;
};
// ============================================================================
@@ -886,7 +899,7 @@ namespace AMDiS {
/** \brief
* Thread safe temporary vector of DimMats for calculation in calculateElementMatrix().
*/
- ::std::vector< DimMat<double>** > tmpLALt;
+ std::vector< DimMat<double>** > tmpLALt;
friend class SecondOrderAssembler;
};
diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc
index 799d46c2a05731578eceed9cf74a996ad289981f..5a624099418143fbaf029deadc843359e9f85a84 100644
--- a/AMDiS/src/DualTraverse.cc
+++ b/AMDiS/src/DualTraverse.cc
@@ -4,19 +4,19 @@
namespace AMDiS {
- bool DualTraverse::traverseFirst(Mesh *mesh1,
- Mesh *mesh2,
- int level1,
- int level2,
- Flag flag1,
- Flag flag2,
+ bool DualTraverse::traverseFirst(Mesh *mesh1,
+ Mesh *mesh2,
+ int level1,
+ int level2,
+ Flag flag1,
+ Flag flag2,
ElInfo **elInfo1,
ElInfo **elInfo2,
ElInfo **elInfoSmall,
ElInfo **elInfoLarge)
{
// replace CALL_EL_LEVEL by CALL_MG_LEVEL (covers whole domain)
- if(flag1.isSet(Mesh::CALL_EL_LEVEL)) {
+ if (flag1.isSet(Mesh::CALL_EL_LEVEL)) {
flag1 &= ~Mesh::CALL_EL_LEVEL;
flag1 |= Mesh::CALL_MG_LEVEL;
level1_ = level1;
@@ -24,7 +24,7 @@ namespace AMDiS {
level1_ = -1;
}
- if(flag2.isSet(Mesh::CALL_EL_LEVEL)) {
+ if (flag2.isSet(Mesh::CALL_EL_LEVEL)) {
flag2 &= ~Mesh::CALL_EL_LEVEL;
flag2 |= Mesh::CALL_MG_LEVEL;
level2_ = level2;
@@ -33,7 +33,7 @@ namespace AMDiS {
}
// replace CALL_LEAF_EL_LEVEL by CALL_MG_LEVEL (covers whole domain)
- if(flag1.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
+ if (flag1.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
flag1 &= ~Mesh::CALL_LEAF_EL_LEVEL;
flag1 |= Mesh::CALL_MG_LEVEL;
level1_ = level1;
@@ -43,7 +43,7 @@ namespace AMDiS {
callLeafElLevel1_ = false;
}
- if(flag2.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
+ if (flag2.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
flag2 &= ~Mesh::CALL_LEAF_EL_LEVEL;
flag2 |= Mesh::CALL_MG_LEVEL;
level2_ = level2;
@@ -55,15 +55,17 @@ namespace AMDiS {
// call standard traverse
*elInfo1 = stack1.traverseFirst(mesh1, level1, flag1);
- while(*elInfo1 != NULL && skipEl1(*elInfo1))
+ while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
*elInfo1 = stack1.traverseNext(*elInfo1);
+ }
*elInfo2 = stack2.traverseFirst(mesh2, level2, flag2);
- while(*elInfo2 != NULL && skipEl2(*elInfo2))
+ while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
*elInfo2 = stack2.traverseNext(*elInfo2);
+ }
// finished ?
if (*elInfo1 == NULL || *elInfo2 == NULL) {
- TEST_EXIT_DBG(*elInfo1 == *elInfo2)("invalid dual traverse\n");
+ TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
return false;
}
@@ -73,14 +75,13 @@ namespace AMDiS {
elInfoSmall, elInfoLarge);
// check for non domain covering level traverse
- if(!accepted ||
- (level1_ != -1 && (*elInfo1)->getLevel() != level1_) ||
- (callLeafElLevel1_ && !(*elInfo1)->getElement()->isLeaf()) ||
- (level2_ != -1 && (*elInfo2)->getLevel() != level2_) ||
- (callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf()))
- {
- return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
- }
+ if (!accepted ||
+ (level1_ != -1 && (*elInfo1)->getLevel() != level1_) ||
+ (callLeafElLevel1_ && !(*elInfo1)->getElement()->isLeaf()) ||
+ (level2_ != -1 && (*elInfo2)->getLevel() != level2_) ||
+ (callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf())) {
+ return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
+ }
return true;
}
@@ -91,12 +92,12 @@ namespace AMDiS {
ElInfo **elInfoLarge)
{
// call standard traverse
- if(inc1_) {
+ if (inc1_) {
do {
*elInfo1 = stack1.traverseNext(*elInfo1);
} while(*elInfo1 != NULL && skipEl1(*elInfo1));
}
- if(inc2_) {
+ if (inc2_) {
do {
*elInfo2 = stack2.traverseNext(*elInfo2);
} while (*elInfo2 != NULL && skipEl2(*elInfo2));
@@ -104,7 +105,7 @@ namespace AMDiS {
// finished ?
if (*elInfo1 == NULL || *elInfo2 == NULL) {
- TEST_EXIT_DBG(*elInfo1 == *elInfo2)("invalid dual traverse\n");
+ TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
return false;
}
@@ -112,15 +113,14 @@ namespace AMDiS {
elInfoSmall, elInfoLarge);
// check for non domain covering level traverse
- if(!accepted ||
- (level1_ != -1 && (*elInfo1)->getLevel() != level1_) ||
- (callLeafElLevel1_ && !(*elInfo1)->getElement()->isLeaf()) ||
- (level2_ != -1 && (*elInfo2)->getLevel() != level2_) ||
- (callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf()))
- {
- return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
- }
-
+ if (!accepted ||
+ (level1_ != -1 && (*elInfo1)->getLevel() != level1_) ||
+ (callLeafElLevel1_ && !(*elInfo1)->getElement()->isLeaf()) ||
+ (level2_ != -1 && (*elInfo2)->getLevel() != level2_) ||
+ (callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf())) {
+ return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
+ }
+
return true;
}
@@ -145,7 +145,7 @@ namespace AMDiS {
// update rest
rest_ -= 1.0 / (1 << (smallLevel - largeLevel));
- if(rest_ < 1e-32) {
+ if (rest_ < 1e-32) {
// large element finished -> increment both elements
rest_ = 1.0;
inc1_ = true;
diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h
index 78b76b9a673de77baa809576a23ef303ed7266d7..ec321921fa8436012af26845ec57f6823633e52f 100644
--- a/AMDiS/src/DualTraverse.h
+++ b/AMDiS/src/DualTraverse.h
@@ -48,12 +48,12 @@ namespace AMDiS {
/** \brief
* Start dual traversal
*/
- bool traverseFirst(Mesh *mesh1,
- Mesh *mesh2,
- int level1,
- int level2,
- Flag flag1,
- Flag flag2,
+ bool traverseFirst(Mesh *mesh1,
+ Mesh *mesh2,
+ int level1,
+ int level2,
+ Flag flag1,
+ Flag flag2,
ElInfo **elInfo1,
ElInfo **elInfo2,
ElInfo **elInfoSmall,
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index a447efc6081c58619d1d330ba92c2ea2051e0f04..c6733b0efc50e664721caa95607062b2947711b8 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -39,7 +39,7 @@ namespace AMDiS {
DegreeOfFreedom **dof = macroInfo->dof;
// === read periodic data =================================
- if (periodicFile&&(strcmp(periodicFile, "") != 0)) {
+ if (periodicFile && (strcmp(periodicFile, "") != 0)) {
WARNING("periodic boundaries may lead to errors in small meshes if element neighbours not set\n");
FILE *file = fopen(periodicFile, "r");
@@ -209,8 +209,7 @@ namespace AMDiS {
for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
if (i != (*(assoc->second))[i])
MSG("association %d: vertex %d -> vertex %d\n",
- assoc->first,
- i, (*(assoc->second))[i]);
+ assoc->first, i, (*(assoc->second))[i]);
}
}
@@ -265,7 +264,7 @@ namespace AMDiS {
macroInfo->dirichletBoundary();
}
- if(mesh->getDim() > 1)
+ if (mesh->getDim() > 1)
boundaryDOFs(mesh);
// initial boundary projections
@@ -936,58 +935,54 @@ namespace AMDiS {
void MacroReader::computeNeighbours(Mesh *mesh)
{
- FUNCNAME("MacroReader::computeNeighbours");
- int i, j, k, l, m;
- int dim = mesh->getDim();
+ FUNCNAME("MacroReader::computeNeighbours()");
+ int dim = mesh->getDim();
FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
- for (i = 0; i < mesh->getNumberOfLeaves(); i++)
- {
- for (k = 0; k < mesh->getGeo(NEIGH); k++)
- {
- mesh->getMacroElement(i)->setOppVertex(k, AMDIS_UNDEFINED);
- mesh->getMacroElement(i)->setNeighbour(k, NULL);
- }
+ for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
+ for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
+ mesh->getMacroElement(i)->setOppVertex(k, AMDIS_UNDEFINED);
+ mesh->getMacroElement(i)->setNeighbour(k, NULL);
}
+ }
- for (i = 0; i < mesh->getNumberOfLeaves(); i++)
- {
- // MSG("search neighbour for element %d\n", i);
- for (k = 0; k < mesh->getGeo(NEIGH); k++)
- {
- if(mesh->getMacroElement(i)->getBoundary(k) != INTERIOR) {
- mesh->getMacroElement(i)->setNeighbour(k, NULL);
- mesh->getMacroElement(i)->setOppVertex(k, -1);
- continue;
- }
-
- if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED)
- {
- if(dim == 1)
- dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
- getElement()->getDOF(k));
- else
- for (l = 0; l < dim; l++)
- dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
- getElement()->
- getDOF((k+l+1)%(dim+1)));
-
- for (j = i + 1; j < mesh->getNumberOfLeaves(); j++) {
- if ((m = mesh->getMacroElement(j)->getElement()->oppVertex(dof)) != -1) {
- mesh->getMacroElement(i)->setNeighbour(k, mesh->getMacroElement(j));
- mesh->getMacroElement(j)->setNeighbour(m, mesh->getMacroElement(i));
- mesh->getMacroElement(i)->setOppVertex(k, m);
- mesh->getMacroElement(j)->setOppVertex(m, k);
- break;
- }
- }
+ for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
+ for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
+ if (mesh->getMacroElement(i)->getBoundary(k) != INTERIOR) {
+ mesh->getMacroElement(i)->setNeighbour(k, NULL);
+ mesh->getMacroElement(i)->setOppVertex(k, -1);
+ continue;
+ }
- TEST_EXIT(j < mesh->getNumberOfLeaves())
- ("could not find neighbour %d of element %d\n", k, i);
- }
+ if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED) {
+ if (dim == 1) {
+ dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
+ getElement()->getDOF(k));
+ } else {
+ for (int l = 0; l < dim; l++)
+ dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
+ getElement()->
+ getDOF((k+l+1)%(dim+1)));
}
+
+ int j = 0;
+ for (j = i + 1; j < mesh->getNumberOfLeaves(); j++) {
+ int m = mesh->getMacroElement(j)->getElement()->oppVertex(dof);
+ if (m != -1) {
+ mesh->getMacroElement(i)->setNeighbour(k, mesh->getMacroElement(j));
+ mesh->getMacroElement(j)->setNeighbour(m, mesh->getMacroElement(i));
+ mesh->getMacroElement(i)->setOppVertex(k, m);
+ mesh->getMacroElement(j)->setOppVertex(m, k);
+ break;
+ }
+ }
+
+ TEST_EXIT(j < mesh->getNumberOfLeaves())
+ ("could not find neighbour %d of element %d\n", k, i);
+ }
}
+ }
return;
@@ -1173,25 +1168,20 @@ namespace AMDiS {
void MacroReader::macroTest(Mesh *mesh, const char *nameneu)
{
- FUNCNAME("MacroReader::macroTest");
-
- int i;
+ FUNCNAME("MacroReader::macroTest()");
- i = macrotest(mesh);
+ int i = macrotest(mesh);
- if (i >= 0)
- {
- ERROR("There is a cycle beginning in macro element %d\n", i);
- ERROR("Entries in MacroElement structures get reordered\n");
- umb(NULL,mesh, umbVkantMacro);
+ if (i >= 0) {
+ ERROR("There is a cycle beginning in macro element %d\n", i);
+ ERROR("Entries in MacroElement structures get reordered\n");
+ umb(NULL, mesh, umbVkantMacro);
- if (nameneu)
- {
- ERROR_EXIT("mesh->feSpace\n");
- //MacroWriter::writeMacro(mesh, nameneu);
- MSG("Reordered mesh written to file %s.\n",nameneu);
- }
+ if (nameneu) {
+ ERROR_EXIT("mesh->feSpace\n");
}
+ }
+
return;
}
@@ -1207,89 +1197,69 @@ namespace AMDiS {
int MacroReader::macrotest(Mesh *mesh)
{
- FUNCNAME("MacroReader::macrotest");
+ FUNCNAME("MacroReader::macrotest()");
int *test;
- int i;
int *zykl;
::std::deque<MacroElement*>::const_iterator macro,mac;
int flg;
::std::deque<MacroElement*>::const_iterator macrolfd;
int zykstart;
-
int dim = mesh->getDim();
- test=GET_MEMORY(int, mesh->getNumberOfMacros());
- zykl=GET_MEMORY(int, mesh->getNumberOfMacros());
+ test = GET_MEMORY(int, mesh->getNumberOfMacros());
+ zykl = GET_MEMORY(int, mesh->getNumberOfMacros());
- for (i=0; i < mesh->getNumberOfMacros(); i++)
- test[i]=0;
-
- zykstart=-1;
+ for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
+ test[i] = 0;
+ }
- macrolfd=mesh->firstMacroElement();
+ zykstart = -1;
+ macrolfd = mesh->firstMacroElement();
- while (macrolfd != mesh->endOfMacroElements())
- {
- if (test[(*macrolfd)->getIndex()]==1)
- {
- macrolfd++;
- }
- else
- {
- for (i=0; i < mesh->getNumberOfMacros(); i++)
- zykl[i]=0;
+ while (macrolfd != mesh->endOfMacroElements()) {
+ if (test[(*macrolfd)->getIndex()] == 1) {
+ macrolfd++;
+ } else {
+ for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
+ zykl[i] = 0;
+ }
- macro=macrolfd;
-
- flg=2;
- do
- {
- if (zykl[(*macro)->getIndex()] == 1)
- {
- flg=0;
- zykstart=(*macro)->getIndex();
- }
- else
- {
- zykl[(*macro)->getIndex()]=1;
+ macro = macrolfd;
+ flg = 2;
+ do {
+ if (zykl[(*macro)->getIndex()] == 1) {
+ flg = 0;
+ zykstart = (*macro)->getIndex();
+ } else {
+ zykl[(*macro)->getIndex()] = 1;
- if (test[(*macro)->getIndex()]==1)
- {
- flg=1;
- }
- else if ((*macro)->getNeighbour(dim) == NULL)
- {
- flg=1;
- test[(*macro)->getIndex()]=1;
- }
- else if ((*macro) == (*macro)->getNeighbour(dim)->getNeighbour(dim))
- {
- flg=1;
- test[(*macro)->getIndex()]=1;
- test[(*macro)->getNeighbour(dim)->getIndex()]=1;
- }
- else
- {
- for(mac=mesh->firstMacroElement();
- (*mac)!=(*macro)->getNeighbour(dim);
- mac++);
- macro=mac;
- }
- }
-
- } while(flg == 2);
+ if (test[(*macro)->getIndex()] == 1) {
+ flg = 1;
+ } else if ((*macro)->getNeighbour(dim) == NULL) {
+ flg = 1;
+ test[(*macro)->getIndex()] = 1;
+ }
+ else if ((*macro) == (*macro)->getNeighbour(dim)->getNeighbour(dim)) {
+ flg = 1;
+ test[(*macro)->getIndex()] = 1;
+ test[(*macro)->getNeighbour(dim)->getIndex()] = 1;
+ } else {
+ for (mac = mesh->firstMacroElement();
+ (*mac)!=(*macro)->getNeighbour(dim);
+ mac++);
+ macro = mac;
+ }
+ }
+ } while(flg == 2);
- if (flg == 1)
- {
- macrolfd++;
- }
- else
- {
- macrolfd=mesh->endOfMacroElements();
- }
- }
+ if (flg == 1) {
+ macrolfd++;
+ } else {
+ macrolfd=mesh->endOfMacroElements();
+ }
}
+ }
FREE_MEMORY(zykl, int, mesh->getNumberOfMacros());
FREE_MEMORY(test, int, mesh->getNumberOfMacros());
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index af4800640c06a4dc70e0289130cecb33e9483406..eecdbf4e1cd2b3e8b576f4ebe1dbb5f874a80e90 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -55,7 +55,7 @@ namespace AMDiS {
void OperatorTerm::setSymmetric(bool symm)
{
- if(symm) {
+ if (symm) {
properties.setFlag(SYMMETRIC);
} else {
properties.unsetFlag(SYMMETRIC);
@@ -68,27 +68,29 @@ namespace AMDiS {
}
void OperatorTerm::lalt(const DimVec<WorldVector<double> >& Lambda,
- const WorldMatrix<double>& matrix,
- DimMat<double>& LALt,
- bool symm,
- double factor)
+ const WorldMatrix<double>& matrix,
+ DimMat<double>& LALt,
+ bool symm,
+ double factor)
{
- int i, j, k, l;
- static const int dimOfWorld = Global::getGeo(WORLD);
- int dim = LALt.getNumRows() - 1;
+ int i, j, k, l;
+ static const int dimOfWorld = Global::getGeo(WORLD);
+ int dim = LALt.getNumRows() - 1;
double val = 0.0;
- if(symm) {
+ if (symm) {
for (i = 0; i <= dim; i++) {
- for (val = k = 0; k < dimOfWorld; k++)
- for(l=0; l < dimOfWorld; l++)
+ val = 0.0;
+ for (k = 0; k < dimOfWorld; k++)
+ for (l = 0; l < dimOfWorld; l++)
val += Lambda[i][k] * matrix[k][l] * Lambda[i][l];
val *= factor;
LALt[i][i] += val;
for (j = i+1; j <= dim; j++) {
- for (val = k = 0; k < dimOfWorld; k++)
- for(l=0; l < dimOfWorld; l++)
+ val = 0.0;
+ for (k = 0; k < dimOfWorld; k++)
+ for (l = 0; l < dimOfWorld; l++)
val += Lambda[i][k] * matrix[k][l] * Lambda[j][l];
val *= factor;
LALt[i][j] += val;
@@ -98,8 +100,9 @@ namespace AMDiS {
} else {
for (i = 0; i <= dim; i++) {
for (j = 0; j <= dim; j++) {
- for (val = k = 0; k < dimOfWorld; k++)
- for(l=0; l < dimOfWorld; l++)
+ val = 0.0;
+ for (k = 0; k < dimOfWorld; k++)
+ for (l = 0; l < dimOfWorld; l++)
val += Lambda[i][k] * matrix[k][l] * Lambda[j][l];
val *= factor;
LALt[i][j] += val;
@@ -455,32 +458,28 @@ namespace AMDiS {
void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
}
}
void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
}
}
void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
lalt(Lambda, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
}
void FctGradient_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
}
}
@@ -488,16 +487,14 @@ namespace AMDiS {
void VecAndGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
}
}
void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt)
const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
}
}
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 8ad3d0e643c9eced7d85e658b4d861423dfb7e16..6e832de2aee450f3e60259a3d6b9dc98a6cb56c0 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -1,7 +1,6 @@
#include "ProblemVec.h"
#include "RecoveryEstimator.h"
#include "Serializer.h"
-#include "DualTraverse.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc
index 867fec4255dd62f6c01c322157b31f1244f1e6ce..de3e250a6724e8fadc5d3bb70766fa7a99411751 100644
--- a/AMDiS/src/Quadrature.cc
+++ b/AMDiS/src/Quadrature.cc
@@ -7,7 +7,7 @@ namespace AMDiS {
const int Quadrature::maxNQuadPoints[4] = {0, 10, 61, 64};
- ::std::list<FastQuadrature*> FastQuadrature::fastQuadList;
+ std::list<FastQuadrature*> FastQuadrature::fastQuadList;
int FastQuadrature::max_points = 0;
Quadrature::Quadrature(const Quadrature& q)
@@ -49,7 +49,7 @@ namespace AMDiS {
val = vec;
} else {
if (static_cast<int>( size) < n_points) {
- size_t new_size = ::std::max(maxNQuadPoints[dim], n_points);
+ size_t new_size = std::max(maxNQuadPoints[dim], n_points);
DELETE [] quad_vec_d;
quad_vec_d = NEW WorldVector<double>[new_size];
size = new_size;
@@ -74,7 +74,7 @@ namespace AMDiS {
{
FUNCNAME("Quadrature::fAtQp()");
- static double *quad_vec = NULL;
+ static double *quad_vec = NULL;
static size_t size = 0;
double *val;
@@ -82,7 +82,7 @@ namespace AMDiS {
val = vec;
} else {
if (static_cast<int>( size) < n_points) {
- size_t new_size = ::std::max(maxNQuadPoints[dim], n_points);
+ size_t new_size = std::max(maxNQuadPoints[dim], n_points);
FREE_MEMORY(quad_vec, double, size);
quad_vec = GET_MEMORY(double, new_size);
size = new_size;
@@ -1414,13 +1414,13 @@ namespace AMDiS {
degree_ = 0;
break;
case 1:
- degree_ = ::std::min(degree_, 19);
+ degree_ = std::min(degree_, 19);
break;
case 2:
- degree_ = ::std::min(degree_, 17);
+ degree_ = std::min(degree_, 17);
break;
case 3:
- degree_ = ::std::min(degree_, 7);
+ degree_ = std::min(degree_, 7);
break;
default:
ERROR_EXIT("invalid dim\n");
@@ -1447,7 +1447,7 @@ namespace AMDiS {
return 0;
}
- double result = 0;
+ double result = 0.0;
// calculate weighted sum over all quadrature-points
for (int i = 0; i < n_points; i++) {
result += w[i] * (*f)((*lambda)[i]);
@@ -1461,7 +1461,7 @@ namespace AMDiS {
{
FUNCNAME("FastQuadrature::FastQuadrature()");
- ::std::list<FastQuadrature*>::iterator fast = fastQuadList.begin();
+ std::list<FastQuadrature*>::iterator fast = fastQuadList.begin();
FastQuadrature *quad_fast;
@@ -1478,7 +1478,7 @@ namespace AMDiS {
fastQuadList.push_front(quad_fast);
- max_points = ::std::max(max_points, quad.getNumPoints());
+ max_points = std::max(max_points, quad.getNumPoints());
} else {
quad_fast = (*fast);
}
@@ -1491,7 +1491,7 @@ namespace AMDiS {
void FastQuadrature::init(Flag init_flag)
{
int dim = quadrature->getDim();
- int n_points = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
int nBasFcts = basisFunctions->getNumber();
DimVec<double> lambda(dim, NO_INIT);
@@ -1501,13 +1501,13 @@ namespace AMDiS {
if (!phi && init_flag.isSet(INIT_PHI)) { // check flag
// allocate memory
- phi = GET_MEMORY(double*, n_points);
- for (int i = 0; i < n_points; i++) {
+ phi = GET_MEMORY(double*, nPoints);
+ for (int i = 0; i < nPoints; i++) {
phi[i] = GET_MEMORY(double, nBasFcts);
}
// fill memory
- for (int i = 0; i< n_points; i++) {
+ for (int i = 0; i< nPoints; i++) {
lambda = quadrature->getLambda(i);
for (int j = 0; j < nBasFcts; j++) {
phi[i][j] = (*(basisFunctions->getPhi(j)))(lambda);
@@ -1523,10 +1523,10 @@ namespace AMDiS {
if (!grdPhi && init_flag.isSet(INIT_GRD_PHI)) { // check flag
// allocate memory
- grdPhi = NEW MatrixOfFixVecs<DimVec<double> >(dim, n_points, nBasFcts, NO_INIT);
+ grdPhi = NEW MatrixOfFixVecs<DimVec<double> >(dim, nPoints, nBasFcts, NO_INIT);
// fill memory
- for (int i = 0; i< n_points; i++) {
+ for (int i = 0; i< nPoints; i++) {
lambda = quadrature->getLambda(i);
for (int j = 0; j < nBasFcts; j++) {
(*(basisFunctions->getGrdPhi(j)))(lambda, (*(grdPhi))[i][j]);
@@ -1542,10 +1542,10 @@ namespace AMDiS {
if (!D2Phi && init_flag.isSet(INIT_D2_PHI)) { // check flag
// allocate memory
- D2Phi = NEW MatrixOfFixVecs<DimMat<double> >(dim, n_points, nBasFcts, NO_INIT);
+ D2Phi = NEW MatrixOfFixVecs<DimMat<double> >(dim, nPoints, nBasFcts, NO_INIT);
// fill memory
- for (int i = 0; i < n_points; i++) {
+ for (int i = 0; i < nPoints; i++) {
lambda = quadrature->getLambda(i);
for (int j = 0; j < nBasFcts; j++) {
(*(basisFunctions->getD2Phi(j)))(lambda, (*(D2Phi))[i][j]);
@@ -1572,12 +1572,12 @@ namespace AMDiS {
basisFunctions = fastQuad.basisFunctions;
quadrature = fastQuad.quadrature;
- int n_points = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
int nBasFcts = basisFunctions->getNumber();
if (fastQuad.phi) {
- phi = GET_MEMORY(double*, n_points);
- for (int i = 0; i < n_points; i++) {
+ phi = GET_MEMORY(double*, nPoints);
+ for (int i = 0; i < nPoints; i++) {
phi[i] = GET_MEMORY(double, nBasFcts);
for (int j = 0; j < nBasFcts; j++) {
phi[i][j] = fastQuad.phi[i][j];
@@ -1586,8 +1586,8 @@ namespace AMDiS {
}
if (fastQuad.grdPhi) {
- grdPhi = NEW MatrixOfFixVecs<DimVec<double> >(dim, n_points, nBasFcts, NO_INIT);
- for (int i = 0; i < n_points; i++) {
+ grdPhi = NEW MatrixOfFixVecs<DimVec<double> >(dim, nPoints, nBasFcts, NO_INIT);
+ for (int i = 0; i < nPoints; i++) {
for (int j = 0; j < nBasFcts; j++) {
(*grdPhi)[i][j] = (*(fastQuad.grdPhi))[i][j];
}
@@ -1595,8 +1595,8 @@ namespace AMDiS {
}
if (fastQuad.D2Phi) {
- D2Phi = NEW MatrixOfFixVecs<DimMat<double> >(dim, n_points, nBasFcts, NO_INIT);
- for (int i = 0; i < n_points; i++) {
+ D2Phi = NEW MatrixOfFixVecs<DimMat<double> >(dim, nPoints, nBasFcts, NO_INIT);
+ for (int i = 0; i < nPoints; i++) {
for (int j = 0; j < nBasFcts; j++) {
(*D2Phi)[i][j] = (*(fastQuad.D2Phi))[i][j];
}
@@ -1606,30 +1606,22 @@ namespace AMDiS {
FastQuadrature::~FastQuadrature()
{
- int n_points = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
- for (int i = 0; i < n_points; i++) {
+ for (int i = 0; i < nPoints; i++) {
FREE_MEMORY(phi[i], double, basisFunctions->getNumber());
}
- FREE_MEMORY(phi, double*, n_points);
+ FREE_MEMORY(phi, double*, nPoints);
DELETE grdPhi;
DELETE D2Phi;
}
const double FastQuadrature::getSecDer(int q,int i ,int j, int m) const {
- return (D2Phi)?(*D2Phi)[q][i][j][m]:0;
+ return (D2Phi) ? (*D2Phi)[q][i][j][m] : 0.0;
}
const VectorOfFixVecs<DimMat<double> > *FastQuadrature::getSecDer(int q) const {
return D2Phi ? (&((*D2Phi)[q])) : NULL;
}
-
- const double FastQuadrature::getGradient(int q,int i ,int j) const {
- return (grdPhi)?(*grdPhi)[q][i][j]:0;
- }
-
- VectorOfFixVecs<DimVec<double> >* FastQuadrature::getGradient(int q) const {
- return (grdPhi)?&((*grdPhi)[q]):NULL;
- }
}
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index 53793b4a739f5c389ac71ad82f81f3359b45f90a..2c5d7c2fa4ff38510e69ac9b76baad8f64b511d4 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -114,7 +114,7 @@ namespace AMDiS {
/** \brief
* Returns \ref name
*/
- inline const ::std::string& getName() {
+ inline const std::string& getName() {
return name;
};
@@ -212,7 +212,7 @@ namespace AMDiS {
/** \brief
* Name of this Quadrature
*/
- ::std::string name;
+ std::string name;
/** \brief
* Quadrature is exact of this degree
@@ -453,12 +453,16 @@ namespace AMDiS {
/** \brief
* Returns (*\ref grdPhi)[q][i][j]
*/
- const double getGradient(int q, int i ,int j) const;
+ inline const double getGradient(int q, int i ,int j) const {
+ return (grdPhi) ? (*grdPhi)[q][i][j] : 0.0;
+ };
/** \brief
* Returns (*\ref grdPhi)[q]
*/
- VectorOfFixVecs<DimVec<double> >* getGradient(int q) const;
+ inline VectorOfFixVecs<DimVec<double> >* getGradient(int q) const {
+ return (grdPhi) ? &((*grdPhi)[q]) : NULL;
+ };
/** \brief
* Returns \ref phi[q][i]
@@ -586,7 +590,7 @@ namespace AMDiS {
/** \brief
* List of all used FastQuadratures
*/
- static ::std::list<FastQuadrature*> fastQuadList;
+ static std::list<FastQuadrature*> fastQuadList;
/** \brief
* Maximal number of quadrature points for all yet initialised FastQuadrature