diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc index 5f00577c1df538aeadeb674b39aaead29f591971..29356cc5e00a27455bfe8e05a667b84c5a197d08 100644 --- a/AMDiS/src/FirstOrderAssembler.cc +++ b/AMDiS/src/FirstOrderAssembler.cc @@ -118,25 +118,21 @@ namespace AMDiS { Lb.resize(nPoints); std::vector<double> phival(nCol); - 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].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 < nCol; i++) { - phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq)); - } + for (int i = 0; i < nCol; i++) + phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq)); for (int i = 0; i < nRow; i++) { (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi); - for (int j = 0; j < nCol; j++) { - mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j]; - } + for (int j = 0; j < nCol; j++) + mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j]; } } } @@ -150,12 +146,10 @@ namespace AMDiS { 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].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(); @@ -298,9 +292,8 @@ namespace AMDiS { k = q10->getKVec(i, j); values = q10->getValVec(i, j); double val = 0.0; - for (int m = 0; m < nEntries[i][j]; m++) { - val += values[m] * Lb[0][k[m]]; - } + for (int m = 0; m < nEntries[i][j]; m++) + val += values[m] * Lb[0][k[m]]; mat[i][j] += val; } } @@ -308,33 +301,33 @@ namespace AMDiS { Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad) : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI) - {} + { + VectorOfFixVecs<DimVec<double> > + grdPhi(assembler->getRowFESpace()->getMesh()->getDim(), nCol, NO_INIT); + tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi); + } void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat) { - VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT); - 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]; + VectorOfFixVecs<DimVec<double> >& Lb = tmpLb[myRank]; + VectorOfFixVecs<DimVec<double> >& grdPhi = tmpGrdPhi[myRank]; Lb.resize(nPoints); - for (int iq = 0; iq < nPoints; 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++) { + 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 < nCol; i++) { + for (int i = 0; i < nCol; i++) (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]); - } for (int i = 0; i < nRow; i++) { double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq)); @@ -370,12 +363,10 @@ namespace AMDiS { VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank]; Lb.resize(nPoints); - for (int iq = 0; iq < nPoints; 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++) { + 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(); @@ -383,10 +374,9 @@ namespace AMDiS { const double *psi = psiFast->getPhi(iq); grdPhi = phiFast->getGradient(iq); - for (int i = 0; i < nRow; i++) { + for (int i = 0; i < nRow; i++) for (int j = 0; j < nCol; j++) mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i]; - } } } diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h index 5f32e6c11eb055e601c2845d1cbd4f0c95d3fbc5..723ee31e6816af5478d8065ee67296209567499e 100644 --- a/AMDiS/src/FirstOrderAssembler.h +++ b/AMDiS/src/FirstOrderAssembler.h @@ -97,9 +97,13 @@ namespace AMDiS { void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat); /// Implements SubAssembler::calculateElementVector(). - void calculateElementVector(const ElInfo*, ElementVector&) { + void calculateElementVector(const ElInfo*, ElementVector&) + { ERROR_EXIT("should not be called\n"); } + + protected: + std::vector<VectorOfFixVecs<DimVec<double> > > tmpGrdPhi; };