From d3b88dd1a94c756a40cdf6a4cb2540afb249bf55 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski <thomas.witkowski@gmx.de> Date: Tue, 2 Sep 2008 11:27:14 +0000 Subject: [PATCH] * This and that --- AMDiS/src/Assembler.cc | 18 +-- AMDiS/src/Assembler.h | 22 ++- AMDiS/src/ElInfo2d.cc | 16 +- AMDiS/src/Operator.cc | 341 ++++++++++++++++------------------------- AMDiS/src/Operator.h | 91 ++++++----- AMDiS/src/Quadrature.h | 4 +- 6 files changed, 209 insertions(+), 283 deletions(-) diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc index cca6cd46..287cc65b 100644 --- a/AMDiS/src/Assembler.cc +++ b/AMDiS/src/Assembler.cc @@ -336,7 +336,7 @@ namespace AMDiS { if (!optimized) { newAssembler = NEW Stand0(op, assembler, quad); } else { - if(pwConst) { + if (pwConst) { newAssembler = NEW Pre0(op, assembler, quad); } else { newAssembler = NEW Quad0(op, assembler, quad); @@ -493,12 +493,9 @@ namespace AMDiS { void Stand0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) { - double val; - const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts(); const BasisFunction *phi = owner->getColFESpace()->getBasisFcts(); - double psival; double *phival = GET_MEMORY(double, nCol); int nPoints = quadrature->getNumPoints(); double *c = GET_MEMORY(double, nPoints); @@ -523,10 +520,10 @@ namespace AMDiS { } for (int i = 0; i < nRow; i++) { - psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq)); + double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq)); (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i]; for (int j = i + 1; j < nCol; j++) { - val = quadrature->getWeight(iq) * c[iq] * psival * phival[j]; + double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j]; (*mat)[i][j] += val; (*mat)[j][i] += val; } @@ -542,7 +539,7 @@ namespace AMDiS { } for (int i = 0; i < nRow; i++) { - psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq)); + double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq)); for (int j = 0; j < nCol; j++) { (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j]; } @@ -1141,7 +1138,6 @@ namespace AMDiS { } tmpMat *= elInfo->getDet(); - nEntries = q11->getNumberEntries(); if (symmetric) { @@ -1274,7 +1270,6 @@ namespace AMDiS { void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat) { - double val; DimVec<double> grdPsi(dim, NO_INIT); VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT); @@ -1309,7 +1304,7 @@ namespace AMDiS { (grdPsi * ((*LALt[iq]) * grdPhi[i])); for (int j = i + 1; j < nCol; j++) { - val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j])); + double val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j])); (*mat)[i][j] += val; (*mat)[j][i] += val; } @@ -1517,13 +1512,10 @@ namespace AMDiS { // create sub assemblers secondOrderAssembler = SecondOrderAssembler::getSubAssembler(op, this, quad2, opt); - firstOrderAssemblerGrdPsi = FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt); - firstOrderAssemblerGrdPhi = FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt); - zeroOrderAssembler = ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt); } diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h index e782eb04..d1b07fae 100644 --- a/AMDiS/src/Assembler.h +++ b/AMDiS/src/Assembler.h @@ -53,13 +53,9 @@ namespace AMDiS { class Q0Psi; class Q1Psi; class Q2Psi; - // class Operator; - // class OperatorTerm; template<typename T> class DOFVectorBase; - // enum FirstOrderType; - // ============================================================================ // ===== class SubAssembler =================================================== // ============================================================================ @@ -1121,17 +1117,17 @@ namespace AMDiS { /** \brief * SubAssembler for the first order terms (grdPsi) */ - FirstOrderAssembler *firstOrderAssemblerGrdPsi; + FirstOrderAssembler *firstOrderAssemblerGrdPsi; /** \brief * SubAssembler for the first order terms (grdPhi) */ - FirstOrderAssembler *firstOrderAssemblerGrdPhi; + FirstOrderAssembler *firstOrderAssemblerGrdPhi; /** \brief * SubAssembler for the zero order terms */ - ZeroOrderAssembler *zeroOrderAssembler; + ZeroOrderAssembler *zeroOrderAssembler; bool remember; @@ -1196,13 +1192,13 @@ namespace AMDiS { /** \brief * Constructor. */ - StandardAssembler(Operator *op, + StandardAssembler(Operator *op, Quadrature *quad2, Quadrature *quad1GrdPsi, Quadrature *quad1GrdPhi, Quadrature *quad0, - const FiniteElemSpace *rowFESpace_, - const FiniteElemSpace *colFESpace_=NULL); + const FiniteElemSpace *rowFESpace, + const FiniteElemSpace *colFESpace = NULL); }; // ============================================================================ @@ -1223,13 +1219,13 @@ namespace AMDiS { /** \brief * Constructor. */ - OptimizedAssembler(Operator *op, + OptimizedAssembler(Operator *op, Quadrature *quad2, Quadrature *quad1GrdPsi, Quadrature *quad1GrdPhi, Quadrature *quad0, - const FiniteElemSpace *rowFESpace_, - const FiniteElemSpace *colFESpace_=NULL); + const FiniteElemSpace *rowFESpace, + const FiniteElemSpace *colFESpace = NULL); }; } diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc index f8c90468..062bf6f2 100644 --- a/AMDiS/src/ElInfo2d.cc +++ b/AMDiS/src/ElInfo2d.cc @@ -59,8 +59,7 @@ namespace AMDiS { MacroElement *macroNeighbour = mel->getNeighbour(i); if (macroNeighbour) { - neighbour_[i] = macroNeighbour->getElement(); - + neighbour_[i] = macroNeighbour->getElement(); Element *nb = const_cast<Element*>(neighbour_[i]); int edgeNo = oppVertex_[i] = mel->getOppVertex(i); @@ -110,6 +109,19 @@ namespace AMDiS { break; default: + std::cout << "------------- Error --------------" << std::endl; + std::cout << " Element index = " << element_->getIndex() << "\n\n"; + for (int j = 0; j < neighbours; j++) { + if (mel->getNeighbour(j)) { + std::cout << " Neighbour " << j << ": " + << mel->getNeighbour(j)->getElement()->getIndex() + << std::endl; + } else { + std::cout << " Neighbour " << j << ": not existing" << std::endl; + } + std::cout << " OppVertex " << j << ": " << static_cast<int>(mel->getOppVertex(j)) << std::endl; + std::cout << std::endl; + } ERROR_EXIT("should not happen!\n"); break; } diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc index b21dd755..06460c68 100644 --- a/AMDiS/src/Operator.cc +++ b/AMDiS/src/Operator.cc @@ -261,7 +261,7 @@ namespace AMDiS { void VecAtQP_SOT::initElement(const ElInfo* elInfo, SubAssembler* subAssembler, - Quadrature *quad) + Quadrature* quad) { vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad); } @@ -450,8 +450,7 @@ namespace AMDiS { void MatrixFct_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, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0); } } @@ -502,8 +501,7 @@ namespace AMDiS { void MatrixGradientAndCoords_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], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0); } } @@ -511,17 +509,15 @@ namespace AMDiS { void VecGradCoordsAtQP_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], coordsAtQPs[iq])); } } void VecAtQP_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); - int iq; - for(iq = 0; iq < numPoints; iq++) { - if(b) + for (int iq = 0; iq < numPoints; iq++) { + if (b) lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq])); else l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq])); @@ -530,8 +526,7 @@ namespace AMDiS { void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { l1(Lambda, Lb[iq], (*g)(coordsAtQPs[iq])); } } @@ -540,22 +535,19 @@ namespace AMDiS { VectorOfFixVecs<DimVec<double> >& Lb) const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); - int iq; - for(iq = 0; iq < numPoints; iq++) - { - lb(Lambda, b, Lb[iq], (*g)(coordsAtQPs[iq])); - } + for (int iq = 0; iq < numPoints; iq++) { + lb(Lambda, b, Lb[iq], (*g)(coordsAtQPs[iq])); + } } void VectorGradient_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); - int iq; - if(f) { - for(iq = 0; iq < numPoints; iq++) { + if (f) { + for (int iq = 0; iq < numPoints; iq++) { lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0); } } else { - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0); } } @@ -563,16 +555,14 @@ namespace AMDiS { void VectorFct_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0); } } void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda(); - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0); } } @@ -596,8 +586,7 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq]; } } @@ -615,13 +604,9 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { - result[iq] += - fac * - (*f1)(vecAtQPs1[iq]) * - (*f2)(vecAtQPs2[iq]) * - uhAtQP[iq]; + for (int iq = 0; iq < numPoints; iq++) { + result[iq] += fac * (*f1)(vecAtQPs1[iq]) * + (*f2)(vecAtQPs2[iq]) * uhAtQP[iq]; } } @@ -638,11 +623,8 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { - result[iq] += - fac * - (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * + for (int iq = 0; iq < numPoints; iq++) { + result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * uhAtQP[iq]; } } @@ -660,11 +642,8 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { - result[iq] += - fac * - (*f)(vecAtQPs[iq], coordsAtQPs[iq]) * + for (int iq = 0; iq < numPoints; iq++) { + result[iq] += fac * (*f)(vecAtQPs[iq], coordsAtQPs[iq]) * uhAtQP[iq]; } } @@ -683,11 +662,8 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { - result[iq] += - fac * - (*f)(gradAtQPs[iq], coordsAtQPs[iq]) * + for (int iq = 0; iq < numPoints; iq++) { + result[iq] += fac * (*f)(gradAtQPs[iq], coordsAtQPs[iq]) * uhAtQP[iq]; } } @@ -705,8 +681,7 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]) * @@ -714,8 +689,6 @@ namespace AMDiS { } } - // bis hierhin - void VecAndGradAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const { for (int iq = 0; iq < numPoints; iq++) { C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]); @@ -729,8 +702,7 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * @@ -752,8 +724,7 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * @@ -774,8 +745,7 @@ namespace AMDiS { double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { result[iq] += fac * (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]) * @@ -826,23 +796,22 @@ namespace AMDiS { double *result, double factor) const { - int i, j, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { double resultQP = 0.0; WorldMatrix<double> A = (*matrixFct)(vecAtQPs[iq]); - if(D2UhAtQP) { - for(i=0; i < dow; i++) { - for(j=0; j < dow; j++) { + if (D2UhAtQP) { + for (int i = 0; i < dow; i++) { + for (int j = 0; j < dow; j++) { resultQP += A[i][j] * D2UhAtQP[iq][j][i]; } } } - if(grdUhAtQP) { + if (grdUhAtQP) { resultQP += (*divFct)(A) * grdUhAtQP[iq]; } @@ -854,10 +823,9 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - int iq; - if(grdUhAtQP) { + if (grdUhAtQP) { WorldMatrix<double> A; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { A = (*matrixFct)(vecAtQPs[iq]); result[iq] += A * grdUhAtQP[iq]; } @@ -872,14 +840,13 @@ namespace AMDiS { double *result, double factor) const { - int i, j, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double resultQP = 0.0; - for(i=0; i < dow; i++) { - for(j=0; j < dow; j++) { + for (int i = 0; i < dow; i++) { + for (int j = 0; j < dow; j++) { resultQP += matrix[i][j] * D2UhAtQP[iq][j][i]; } } @@ -892,9 +859,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { result[iq] += matrix * grdUhAtQP[iq]; } } @@ -909,23 +875,22 @@ namespace AMDiS { double *result, double factor) const { - int i, j, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { double resultQP = 0.0; WorldMatrix<double> A = (*f)(gradAtQPs[iq]); - if(D2UhAtQP) { - for(i=0; i < dow; i++) { - for(j=0; j < dow; j++) { + if (D2UhAtQP) { + for (int i = 0; i < dow; i++) { + for (int j = 0; j < dow; j++) { resultQP += A[i][j] * D2UhAtQP[iq][j][i]; } } } - if(grdUhAtQP) { + if (grdUhAtQP) { resultQP += (*divFct)(A) * grdUhAtQP[iq]; } @@ -937,10 +902,9 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; + if (grdUhAtQP) { WorldMatrix<double> A; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { A = (*f)(gradAtQPs[iq]); result[iq] += A * grdUhAtQP[iq]; } @@ -954,14 +918,13 @@ namespace AMDiS { double *result, double fac) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += fac * factor * resultQP; @@ -973,9 +936,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq]); axpy(factor, grdUhAtQP[iq], result[iq]); } @@ -989,14 +951,13 @@ namespace AMDiS { double *result, double f) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += factor * f * resultQP; @@ -1008,9 +969,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); axpy(factor, grdUhAtQP[iq], result[iq]); } @@ -1024,15 +984,13 @@ namespace AMDiS { double *result, double fac) const { - int i, dow = Global::getGeo(WORLD); - - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(gradAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += fac * factor * resultQP; @@ -1044,9 +1002,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(gradAtQPs[iq]); axpy(factor, grdUhAtQP[iq], result[iq]); } @@ -1060,15 +1017,13 @@ namespace AMDiS { double *result, double fac) const { - int i, dow = Global::getGeo(WORLD); - - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += fac * factor * resultQP; @@ -1080,9 +1035,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]); axpy(factor, grdUhAtQP[iq], result[iq]); } @@ -1098,14 +1052,13 @@ namespace AMDiS { double *result, double fac) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += fac * factor * resultQP; @@ -1117,9 +1070,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]); axpy(factor, grdUhAtQP[iq], result[iq]); } @@ -1133,17 +1085,16 @@ namespace AMDiS { double *result, double factor) const { - int i, j, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { double resultQP = 0.0; WorldMatrix<double> A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]); - if(D2UhAtQP) { - for(i=0; i < dow; i++) { - for(j=0; j < dow; j++) { + if (D2UhAtQP) { + for (int i = 0; i < dow; i++) { + for (int j = 0; j < dow; j++) { resultQP += A[i][j] * D2UhAtQP[iq][j][i]; } } @@ -1161,10 +1112,9 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; + if (grdUhAtQP) { WorldMatrix<double> A; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]); result[iq] += A * grdUhAtQP[iq]; } @@ -1178,15 +1128,13 @@ namespace AMDiS { double *result, double fac) const { - int i, dow = Global::getGeo(WORLD); - - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += fac * factor * resultQP; @@ -1198,9 +1146,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]); axpy(factor, grdUhAtQP[iq], result[iq]); } @@ -1214,14 +1161,13 @@ namespace AMDiS { double *result, double fac) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(grdUhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += grdUhAtQP[iq][i]; } result[iq] += fac * factor * resultQP; @@ -1236,14 +1182,13 @@ namespace AMDiS { double *result, double f) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(grdUhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += grdUhAtQP[iq][i]; } result[iq] += f * factor * resultQP; @@ -1258,14 +1203,13 @@ namespace AMDiS { double *result, double f) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(grdUhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += grdUhAtQP[iq][i]; } result[iq] += f * factor * resultQP; @@ -1279,17 +1223,15 @@ namespace AMDiS { const WorldMatrix<double> *D2UhAtQP, double *result, double factor) const - { - int iq; - - if(grdUhAtQP) { - if(f) { - for(iq = 0; iq < numPoints; iq++) { + { + if (grdUhAtQP) { + if (f) { + for (int iq = 0; iq < numPoints; iq++) { WorldVector<double> b = (*f)(gradAtQPs[iq]); result[iq] += b * grdUhAtQP[iq] * factor; } } else { - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < numPoints; iq++) { result[iq] += gradAtQPs[iq] * grdUhAtQP[iq] * factor; } } @@ -1303,10 +1245,8 @@ namespace AMDiS { double *result, double factor) const { - int iq; - - if(grdUhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { WorldVector<double> b = (*vecFct)(vecAtQPs[iq]); result[iq] += b * grdUhAtQP[iq] * factor; } @@ -1320,14 +1260,13 @@ namespace AMDiS { double *result, double f) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(grdUhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double resultQP = 0.0; const WorldVector<double> &b = (*g)(coordsAtQPs[iq]); - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += b[i] * grdUhAtQP[iq][i]; } result[iq] += f * resultQP; @@ -1361,8 +1300,7 @@ namespace AMDiS { 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_kl(Lambda, xi, xj, *(LALt[iq]), (*g)(coordsAtQPs[iq])); } } @@ -1374,10 +1312,8 @@ namespace AMDiS { double *result, double f) const { - int iq; - - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); result[iq] = D2UhAtQP[iq][xi][xj] * factor * f; } @@ -1388,9 +1324,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); result[iq][xi] += grdUhAtQP[iq][xj] * factor; } @@ -1409,8 +1344,7 @@ namespace AMDiS { 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_kl(Lambda, xi, xj, *(LALt[iq]), (*f)(vecAtQPs[iq])); } } @@ -1422,10 +1356,8 @@ namespace AMDiS { double *result, double fac) const { - int iq; - - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq]); result[iq] = D2UhAtQP[iq][xi][xj] * factor * fac; } @@ -1436,9 +1368,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - if(grdUhAtQP) { - int iq; - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double factor = (*f)(vecAtQPs[iq]); result[iq][xi] += grdUhAtQP[iq][xj] * factor; } @@ -1450,8 +1381,8 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - int i ,size = static_cast<int>(vecs.size()); - for(i = 0; i < size; i++) { + int size = static_cast<int>(vecs.size()); + for (int i = 0; i < size; i++) { vecsAtQPs[i] = subAssembler->getVectorAtQPs(vecs[i], elInfo, quad); } } @@ -1476,13 +1407,11 @@ namespace AMDiS { double *result, double fac) const { - int i ,size = static_cast<int>(vecs.size()); - + int size = static_cast<int>(vecs.size()); std::vector<double> arg(size); - int iq; - for(iq = 0; iq < numPoints; iq++) { - for(i = 0; i < size; i++) { + for (int iq = 0; iq < numPoints; iq++) { + for (int i = 0; i < size; i++) { arg[i] = vecsAtQPs[i][iq]; } result[iq] += fac * (*f)(arg) * uhAtQP[iq]; @@ -1504,8 +1433,8 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - int i ,size = static_cast<int>(vecs.size()); - for(i = 0; i < size; i++) { + int size = static_cast<int>(vecs.size()); + for (int i = 0; i < size; i++) { gradsAtQPs[i] = subAssembler->getGradientsAtQPs(vecs[i], elInfo, quad); } } diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h index 55aefe30..9942094d 100644 --- a/AMDiS/src/Operator.h +++ b/AMDiS/src/Operator.h @@ -92,8 +92,8 @@ namespace AMDiS { * and coordinates at quadrature points can be calculated. */ virtual void initElement(const ElInfo*, - SubAssembler* , - Quadrature *quad= NULL) + SubAssembler*, + Quadrature *quad = NULL) {}; /** \brief @@ -124,11 +124,11 @@ namespace AMDiS { /** \brief * Evaluation of the OperatorTerm at all quadrature points. */ - virtual void eval(int numPoints, - const double *uhAtQP, + virtual void eval(int nPoints, + const double *uhAtQP, const WorldVector<double> *grdUhAtQP, const WorldMatrix<double> *D2UhAtQP, - double *result, + double *result, double factor) const = 0; protected: @@ -136,10 +136,10 @@ namespace AMDiS { * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$. */ static void 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); /** \brief * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for \f$ A \f$ @@ -148,35 +148,35 @@ namespace AMDiS { */ static void lalt_kl(const DimVec<WorldVector<double> >& Lambda, int k, int l, - DimMat<double>& LALt, - double factor); + DimMat<double>& LALt, + double factor); /** \brief * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the * identity. */ static void l1lt(const DimVec<WorldVector<double> >& Lambda, - DimMat<double>& LALt, - double factor); + DimMat<double>& LALt, + double factor); /** \brief * Evaluation of \f$ \Lambda \cdot b\f$. */ static void lb(const DimVec<WorldVector<double> >& Lambda, - const WorldVector<double>& b, - DimVec<double>& Lb, - double factor); + const WorldVector<double>& b, + DimVec<double>& Lb, + double factor); /** \brief * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in * each component. */ static void l1(const DimVec<WorldVector<double> >& Lambda, - DimVec<double>& Lb, - double factor) + DimVec<double>& Lb, + double factor) { - int dim = Lb.getSize() - 1; - static const int dimOfWorld = Global::getGeo(WORLD); + int dim = Lb.getSize() - 1; + static const int dimOfWorld = Global::getGeo(WORLD); for (int i = 0; i <= dim; i++) { double val = 0.0; @@ -250,13 +250,13 @@ namespace AMDiS { * Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points. */ virtual void getLALt(const ElInfo *elInfo, - int numPoints, + int nPoints, DimMat<double> **result) const = 0; /** \brief * Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points. */ - virtual void weakEval(int numPoints, + virtual void weakEval(int nPoints, const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const = 0; @@ -297,17 +297,17 @@ namespace AMDiS { /** \brief * Implementation of SecondOrderTerm::eval(). */ - inline void eval(int numPoints, - const double * // uhAtQP - , const WorldVector<double>*// grdUhAtQP - , const WorldMatrix<double> *D2UhAtQP, - double *result, - double factor) const + inline void eval(int nPoints, + const double * , // uhAtQP + const WorldVector<double> * , // grdUhAtQP + const WorldMatrix<double> *D2UhAtQP, + double *result, + double factor) const { int dow = Global::getGeo(WORLD); if (D2UhAtQP) { - for (int iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < nPoints; iq++) { double resultQP = 0.0; for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; @@ -320,12 +320,12 @@ namespace AMDiS { /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ - void weakEval(int numPoints, + void weakEval(int nPoints, const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { if (grdUhAtQP) { - for (int iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < nPoints; iq++) { result[iq] += grdUhAtQP[iq]; } } @@ -387,13 +387,12 @@ namespace AMDiS { double *result, double f) const { - int i, dow = Global::getGeo(WORLD); - int iq; + int dow = Global::getGeo(WORLD); - if(D2UhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (D2UhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { double resultQP = 0.0; - for(i=0; i < dow; i++) { + for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += resultQP * f * (*factor); @@ -408,9 +407,8 @@ namespace AMDiS { const WorldVector<double> *grdUhAtQP, WorldVector<double> *result) const { - int iq; - if(grdUhAtQP) { - for(iq = 0; iq < numPoints; iq++) { + if (grdUhAtQP) { + for (int iq = 0; iq < numPoints; iq++) { axpy(*factor, grdUhAtQP[iq], result[iq]); } } @@ -2150,9 +2148,9 @@ namespace AMDiS { /** \brief * Implements ZeroOrderTerm::getC(). */ - inline void getC(const ElInfo *, int numPoints, double *C) const + inline void getC(const ElInfo *, int nPoints, double *C) const { - for (int iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < nPoints; iq++) { C[iq] += factor; } }; @@ -2160,15 +2158,14 @@ namespace AMDiS { /** \brief * Implements ZeroOrderTerm::eval(). */ - inline void eval(int numPoints, - const double *uhAtQP, + inline void eval(int nPoints, + const double *uhAtQP, const WorldVector<double> *, const WorldMatrix<double> *, double *result, double fac) const { - int iq; - for(iq = 0; iq < numPoints; iq++) { + for (int iq = 0; iq < nPoints; iq++) { result[iq] += fac * factor * uhAtQP[iq]; } }; @@ -3471,8 +3468,8 @@ namespace AMDiS { * or on both sides of the system. */ Operator(Flag operatorType, - const FiniteElemSpace *rowFESpace_, - const FiniteElemSpace *colFESpace_ = NULL); + const FiniteElemSpace *rowFESpace, + const FiniteElemSpace *colFESpace = NULL); /** \brief * Destructor. diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h index 2c5d7c2f..b2403a53 100644 --- a/AMDiS/src/Quadrature.h +++ b/AMDiS/src/Quadrature.h @@ -95,9 +95,9 @@ namespace AMDiS { Quadrature(const Quadrature&); /** \brief - * Returns a Quadrature for dimension dim_ exact for degree degree_. + * Returns a Quadrature for dimension dim exact for degree degree. */ - static Quadrature *provideQuadrature(int dim_, int degree_); + static Quadrature *provideQuadrature(int dim, int degree); /** \brief * Approximates an integral by the numerical quadrature described by quad; -- GitLab