Commit cd220992 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

ElementMatrix and ElementVector are now MTL4 data structures.

parent 386a66c2
#include <boost/numeric/mtl/mtl.hpp>
#include "CompositeFEMOperator.h" #include "CompositeFEMOperator.h"
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "OpenMP.h" #include "OpenMP.h"
#include "SubElementAssembler.h" #include "SubElementAssembler.h"
#include "SubElInfo.h" #include "SubElInfo.h"
#include "SubPolytope.h" #include "SubPolytope.h"
void void
CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo, CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
ElementMatrix *userMat, ElementMatrix& userMat,
double factor) double factor)
{ {
FUNCNAME("CompositeFEMOperator::getElementMatrix"); FUNCNAME("CompositeFEMOperator::getElementMatrix");
...@@ -18,9 +15,6 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo, ...@@ -18,9 +15,6 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
VectorOfFixVecs<DimVec<double> > *intersecPoints = NULL; VectorOfFixVecs<DimVec<double> > *intersecPoints = NULL;
SubPolytope *subPolytope = NULL; SubPolytope *subPolytope = NULL;
double levelSetSubPolytope; double levelSetSubPolytope;
ElementMatrix *elMat;
ElementMatrix *subPolMat1;
ElementMatrix *subPolMat2;
DimVec<double> subElVertexBarCoords(elInfo->getMesh()->getDim()); DimVec<double> subElVertexBarCoords(elInfo->getMesh()->getDim());
/** /**
...@@ -99,9 +93,9 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo, ...@@ -99,9 +93,9 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
ERROR_EXIT("cannot get position of subpolytope\n"); ERROR_EXIT("cannot get position of subpolytope\n");
} }
subPolMat1 = NEW ElementMatrix(subElementAssembler->getNRow(), ElementMatrix subPolMat1(subElementAssembler->getNRow(),
subElementAssembler->getNCol()); subElementAssembler->getNCol());
subPolMat1->set(0.0); set_to_zero(subPolMat1);
subElementAssembler->getSubPolytopeMatrix(subPolytope, subElementAssembler->getSubPolytopeMatrix(subPolytope,
subElementAssembler, subElementAssembler,
elInfo, elInfo,
...@@ -110,12 +104,12 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo, ...@@ -110,12 +104,12 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
/** /**
* Integration on second subpolytope produced by the intersection. * Integration on second subpolytope produced by the intersection.
*/ */
elMat = NEW ElementMatrix(subElementAssembler->getNRow(), ElementMatrix elMat(subElementAssembler->getNRow(),
subElementAssembler->getNCol()); subElementAssembler->getNCol());
elMat->set(0.0); set_to_zero(elMat);
subPolMat2 = NEW ElementMatrix(subElementAssembler->getNRow(), ElementMatrix subPolMat2(subElementAssembler->getNRow(),
subElementAssembler->getNCol()); subElementAssembler->getNCol());
subPolMat2->set(0.0); set_to_zero(subPolMat2);
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
...@@ -138,28 +132,21 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo, ...@@ -138,28 +132,21 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
elInfo, elInfo,
subPolMat2); subPolMat2);
axpy(-1.0, *subPolMat2, *elMat); elMat -= subPolMat2;
// Get integral on element as sum of the two integrals on subpolytopes. // Get integral on element as sum of the two integrals on subpolytopes.
axpy(1.0, *subPolMat1, *elMat); elMat += subPolMat1;
// Add integral to userMat. // Add integral to userMat.
axpy(factor, *elMat, *userMat); userMat += factor * elMat;
/**
* Free data.
*/
DELETE subPolytope;
DELETE elMat;
DELETE subPolMat1;
DELETE subPolMat2;
return; // Free data
delete subPolytope;
} }
void void
CompositeFEMOperator::getElementVector(const ElInfo *elInfo, CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
ElementVector *userVec, ElementVector& userVec,
double factor) double factor)
{ {
FUNCNAME("CompositeFEMOperator::getElementVector"); FUNCNAME("CompositeFEMOperator::getElementVector");
...@@ -167,9 +154,6 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo, ...@@ -167,9 +154,6 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
VectorOfFixVecs<DimVec<double> >*intersecPoints = NULL; VectorOfFixVecs<DimVec<double> >*intersecPoints = NULL;
SubPolytope *subPolytope = NULL; SubPolytope *subPolytope = NULL;
double levelSetSubPolytope; double levelSetSubPolytope;
ElementVector *elVec;
ElementVector *subPolVec1;
ElementVector *subPolVec2;
DimVec<double> subElVertexBarCoords(elInfo->getMesh()->getDim()); DimVec<double> subElVertexBarCoords(elInfo->getMesh()->getDim());
/** /**
...@@ -247,8 +231,8 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo, ...@@ -247,8 +231,8 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
ERROR_EXIT("cannot get position of subpolytope\n"); ERROR_EXIT("cannot get position of subpolytope\n");
} }
subPolVec1 = NEW ElementVector(subElementAssembler->getNRow()); ElementVector subPolVec1(subElementAssembler->getNRow());
subPolVec1->set(0.0); set_to_zero(subPolVec1);
subElementAssembler->getSubPolytopeVector(subPolytope, subElementAssembler->getSubPolytopeVector(subPolytope,
subElementAssembler, subElementAssembler,
elInfo, elInfo,
...@@ -257,10 +241,10 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo, ...@@ -257,10 +241,10 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
/** /**
* Integration on second subpolytope produced by the intersection. * Integration on second subpolytope produced by the intersection.
*/ */
elVec = NEW ElementVector(subElementAssembler->getNRow()); ElementVector elVec(subElementAssembler->getNRow());
elVec->set(0.0); set_to_zero(elVec);
subPolVec2 = NEW ElementVector(subElementAssembler->getNRow()); ElementVector subPolVec2(subElementAssembler->getNRow());
subPolVec2->set(0.0); set_to_zero(subPolVec2);
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
...@@ -282,21 +266,14 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo, ...@@ -282,21 +266,14 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
elInfo, elInfo,
subPolVec2); subPolVec2);
axpy(-1.0, *subPolVec2, *elVec); elVec -= subPolVec2;
// Get integral on element as sum of the two integrals on subpolytopes. // Get integral on element as sum of the two integrals on subpolytopes.
axpy(1.0, *subPolVec1, *elVec); elVec += subPolVec1;
// Add integral to userVec. // Add integral to userVec.
axpy(factor, *elVec, *userVec); userVec += factor * elVec;
/** // Free data
* Free data. delete subPolytope;
*/
DELETE subPolytope;
DELETE elVec;
DELETE subPolVec1;
DELETE subPolVec2;
return;
} }
...@@ -6,7 +6,6 @@ ...@@ -6,7 +6,6 @@
#include "FixVec.h" #include "FixVec.h"
#include "Flag.h" #include "Flag.h"
#include "ElementLevelSet.h" #include "ElementLevelSet.h"
#include "MemoryManager.h"
#include "Operator.h" #include "Operator.h"
#include "ElementLevelSet.h" #include "ElementLevelSet.h"
...@@ -35,10 +34,7 @@ using namespace std; ...@@ -35,10 +34,7 @@ using namespace std;
class CompositeFEMOperator : public Operator class CompositeFEMOperator : public Operator
{ {
public: public:
MEMORY_MANAGED(CompositeFEMOperator); /// Constructor.
/**
* Constructor.
*/
CompositeFEMOperator(Flag operatorType_, CompositeFEMOperator(Flag operatorType_,
ElementLevelSet *elLS_, ElementLevelSet *elLS_,
const FiniteElemSpace *rowFESpace_, const FiniteElemSpace *rowFESpace_,
...@@ -67,7 +63,7 @@ public: ...@@ -67,7 +63,7 @@ public:
* integration domain. * integration domain.
*/ */
void getElementMatrix(const ElInfo *elInfo, void getElementMatrix(const ElInfo *elInfo,
ElementMatrix *userMat, ElementMatrix& userMat,
double factor = 1.0); double factor = 1.0);
/** /**
...@@ -78,7 +74,7 @@ public: ...@@ -78,7 +74,7 @@ public:
* integration domain. * integration domain.
*/ */
void getElementVector(const ElInfo *elInfo, void getElementVector(const ElInfo *elInfo,
ElementVector *userVec, ElementVector& userVec,
double factor = 1.0); double factor = 1.0);
protected: protected:
......
#include "PenaltyOperator.h" #include "PenaltyOperator.h"
#include "SurfaceOperator.h" #include "SurfaceOperator.h"
double double
...@@ -17,7 +16,7 @@ PenaltyOperator::getPenaltyCoeff(const ElInfo *elInfo) ...@@ -17,7 +16,7 @@ PenaltyOperator::getPenaltyCoeff(const ElInfo *elInfo)
void void
PenaltyOperator::getElementMatrix(const ElInfo *elInfo, PenaltyOperator::getElementMatrix(const ElInfo *elInfo,
ElementMatrix *userMat, ElementMatrix& userMat,
double factor) double factor)
{ {
FUNCNAME("PenaltyOperator::getElementMatrix"); FUNCNAME("PenaltyOperator::getElementMatrix");
...@@ -135,7 +134,7 @@ PenaltyOperator::getElementMatrix(const ElInfo *elInfo, ...@@ -135,7 +134,7 @@ PenaltyOperator::getElementMatrix(const ElInfo *elInfo,
void void
PenaltyOperator::getElementVector(const ElInfo *elInfo, PenaltyOperator::getElementVector(const ElInfo *elInfo,
ElementVector *userVec, ElementVector& userVec,
double factor) double factor)
{ {
FUNCNAME("PenaltyOperator::getElementVector"); FUNCNAME("PenaltyOperator::getElementVector");
...@@ -189,12 +188,12 @@ PenaltyOperator::getElementVector(const ElInfo *elInfo, ...@@ -189,12 +188,12 @@ PenaltyOperator::getElementVector(const ElInfo *elInfo,
else { else {
surfaceOp->adaptSurfaceOperator((*tempCoords)); surfaceOp->adaptSurfaceOperator((*tempCoords));
} }
surfaceOp->getElementVector(elInfo, userVec, factor*penaltyCoeff); surfaceOp->getElementVector(elInfo, userVec, factor * penaltyCoeff);
// Treat second simplex. // Treat second simplex.
(*tempCoords)[0] = (*intersecPoints)[3]; (*tempCoords)[0] = (*intersecPoints)[3];
surfaceOp->adaptSurfaceOperator((*tempCoords)); surfaceOp->adaptSurfaceOperator((*tempCoords));
surfaceOp->getElementVector(elInfo, userVec, factor*penaltyCoeff); surfaceOp->getElementVector(elInfo, userVec, factor * penaltyCoeff);
} }
else { else {
......
...@@ -5,10 +5,8 @@ ...@@ -5,10 +5,8 @@
#include "ElementLevelSet.h" #include "ElementLevelSet.h"
#include "Flag.h" #include "Flag.h"
#include "MemoryManager.h"
#include "Operator.h" #include "Operator.h"
#include "SurfaceOperator.h" #include "SurfaceOperator.h"
#include "ElementLevelSet.h" #include "ElementLevelSet.h"
namespace AMDiS { namespace AMDiS {
...@@ -29,7 +27,6 @@ using namespace std; ...@@ -29,7 +27,6 @@ using namespace std;
class PenaltyOperator : public Operator class PenaltyOperator : public Operator
{ {
public: public:
MEMORY_MANAGED(PenaltyOperator);
/** /**
* Constructor. * Constructor.
*/ */
...@@ -76,7 +73,7 @@ public: ...@@ -76,7 +73,7 @@ public:
* domain. * domain.
*/ */
void getElementMatrix(const ElInfo *elInfo, void getElementMatrix(const ElInfo *elInfo,
ElementMatrix *userMat, ElementMatrix& userMat,
double factor = 1.0); double factor = 1.0);
/** /**
...@@ -87,7 +84,7 @@ public: ...@@ -87,7 +84,7 @@ public:
* the integration domain. * the integration domain.
*/ */
void getElementVector(const ElInfo *elInfo, void getElementVector(const ElInfo *elInfo,
ElementVector *userVec, ElementVector& userVec,
double factor = 1.0); double factor = 1.0);
protected: protected:
......
#include <vector> #include <vector>
#include "ElementMatrix.h" #include <boost/numeric/mtl/mtl.hpp>
#include "ElementVector.h"
#include "SubElementAssembler.h" #include "SubElementAssembler.h"
#include "ScalableQuadrature.h" #include "ScalableQuadrature.h"
#include "SubPolytope.h" #include "SubPolytope.h"
...@@ -72,7 +71,7 @@ namespace AMDiS { ...@@ -72,7 +71,7 @@ namespace AMDiS {
void SubElementAssembler::getSubElementVector(SubElInfo *subElInfo, void SubElementAssembler::getSubElementVector(SubElInfo *subElInfo,
const ElInfo *elInfo, const ElInfo *elInfo,
ElementVector *userVec) ElementVector& userVec)
{ {
/** /**
* Manipulate the quadratures of the SubAssemblers for subelement. * Manipulate the quadratures of the SubAssemblers for subelement.
...@@ -87,14 +86,13 @@ namespace AMDiS { ...@@ -87,14 +86,13 @@ namespace AMDiS {
* the result must be corrected with respect to subelement. * the result must be corrected with respect to subelement.
*/ */
double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet()); double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++)
(*userVec)[i] *= corrFactor; userVec[i] *= corrFactor;
}
} }
void SubElementAssembler::getSubElementMatrix(SubElInfo *subElInfo, void SubElementAssembler::getSubElementMatrix(SubElInfo *subElInfo,
const ElInfo *elInfo, const ElInfo *elInfo,
ElementMatrix *userMat) ElementMatrix& userMat)
{ {
/** /**
* Manipulate the quadratures of the SubAssemblers for subelement. * Manipulate the quadratures of the SubAssemblers for subelement.
...@@ -114,7 +112,7 @@ namespace AMDiS { ...@@ -114,7 +112,7 @@ namespace AMDiS {
double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet()); double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) { for (int j = 0; j < nCol; j++) {
(*userMat)[i][j] *= corrFactor; userMat[i][j] *= corrFactor;
} }
} }
} }
...@@ -122,44 +120,34 @@ namespace AMDiS { ...@@ -122,44 +120,34 @@ namespace AMDiS {
void SubElementAssembler::getSubPolytopeVector(SubPolytope *subPolytope, void SubElementAssembler::getSubPolytopeVector(SubPolytope *subPolytope,
SubElementAssembler *subElementAssembler, SubElementAssembler *subElementAssembler,
const ElInfo *elInfo, const ElInfo *elInfo,
ElementVector *subPolVec) ElementVector& subPolVec)
{ {
/** /// Note: There is no reset of subPolVec.
* Note: There is no reset of subPolVec.
*/
std::vector<SubElInfo *>::iterator it; std::vector<SubElInfo *>::iterator it;
ElementVector *subElVec = NEW ElementVector(nRow); ElementVector subElVec(nRow);
/** /// Assemble for each subelement of subpolytope.
* Assemble for each subelement of subpolytope.
*/
for (it = subPolytope->getSubElementsBegin(); for (it = subPolytope->getSubElementsBegin();
it != subPolytope->getSubElementsEnd(); it != subPolytope->getSubElementsEnd();
it++) { it++) {
subElVec->set(0.0); set_to_zero(subElVec);
subElementAssembler->getSubElementVector(*it, elInfo, subElVec); subElementAssembler->getSubElementVector(*it, elInfo, subElVec);
/** /// Add results for subelement to total result for subpolytope.
* Add results for subelement to total result for subpolytope. subPolVec += subElVec;
*/
for (int i = 0; i < nRow; i++) {
(*subPolVec)[i] += (*subElVec)[i];
} }
} }
DELETE subElVec;
}
void SubElementAssembler::getSubPolytopeMatrix(SubPolytope *subPolytope, void SubElementAssembler::getSubPolytopeMatrix(SubPolytope *subPolytope,
SubElementAssembler *subElementAssembler, SubElementAssembler *subElementAssembler,
const ElInfo *elInfo, const ElInfo *elInfo,
ElementMatrix *subPolMat) ElementMatrix& subPolMat)
{ {
/** /**
* Note: There is no reset of subPolMat. * Note: There is no reset of subPolMat.
*/ */
std::vector<SubElInfo *>::iterator it; std::vector<SubElInfo *>::iterator it;
ElementMatrix *subElMat = NEW ElementMatrix(nRow, nCol); ElementMatrix subElMat(nRow, nCol);
/** /**
* Assemble for each subelement of subpolytope. * Assemble for each subelement of subpolytope.
...@@ -167,20 +155,16 @@ namespace AMDiS { ...@@ -167,20 +155,16 @@ namespace AMDiS {
for (it = subPolytope->getSubElementsBegin(); for (it = subPolytope->getSubElementsBegin();
it != subPolytope->getSubElementsEnd(); it != subPolytope->getSubElementsEnd();
it++) { it++) {
subElMat->set(0.0); set_to_zero(subElMat);
subElementAssembler->getSubElementMatrix(*it, elInfo, subElMat); subElementAssembler->getSubElementMatrix(*it, elInfo, subElMat);
/** /**
* Add results for subelement to total result for subpolytope. * Add results for subelement to total result for subpolytope.
*/ */
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++)
for (int j = 0; j < nCol; j++) { for (int j = 0; j < nCol; j++)
(*subPolMat)[i][j] += (*subElMat)[i][j]; subPolMat[i][j] += subElMat[i][j];
}
} }
} }
DELETE subElMat;
}
} }
...@@ -3,7 +3,6 @@ ...@@ -3,7 +3,6 @@
#ifndef AMDIS_SUBELEMENTASSEMBLER_H #ifndef AMDIS_SUBELEMENTASSEMBLER_H
#define AMDIS_SUBELEMENTASSEMBLER_H #define AMDIS_SUBELEMENTASSEMBLER_H
#include "MemoryManager.h"
#include "Assembler.h" #include "Assembler.h"
#include "SubElInfo.h" #include "SubElInfo.h"
#include "ScalableQuadrature.h" #include "ScalableQuadrature.h"
...@@ -54,8 +53,6 @@ namespace AMDiS { ...@@ -54,8 +53,6 @@ namespace AMDiS {
class SubElementAssembler : public StandardAssembler class SubElementAssembler : public StandardAssembler
{ {
public: public:
MEMORY_MANAGED(SubElementAssembler);
SubElementAssembler(Operator *op, SubElementAssembler(Operator *op,
const FiniteElemSpace *rowFESpace, const FiniteElemSpace *rowFESpace,
const FiniteElemSpace *colFESpace = NULL); const FiniteElemSpace *colFESpace = NULL);
...@@ -63,34 +60,34 @@ namespace AMDiS { ...@@ -63,34 +60,34 @@ namespace AMDiS {
virtual ~SubElementAssembler() virtual ~SubElementAssembler()
{ {
if (zeroOrderScalableQuadrature) if (zeroOrderScalableQuadrature)
DELETE zeroOrderScalableQuadrature; delete zeroOrderScalableQuadrature;
if (firstOrderGrdPsiScalableQuadrature) if (firstOrderGrdPsiScalableQuadrature)
DELETE firstOrderGrdPsiScalableQuadrature; delete firstOrderGrdPsiScalableQuadrature;
if (firstOrderGrdPhiScalableQuadrature) if (firstOrderGrdPhiScalableQuadrature)
DELETE firstOrderGrdPhiScalableQuadrature; delete firstOrderGrdPhiScalableQuadrature;
if (secondOrderScalableQuadrature) if (secondOrderScalableQuadrature)
DELETE secondOrderScalableQuadrature; delete secondOrderScalableQuadrature;
}; }
void scaleQuadratures(const SubElInfo& subElInfo); void scaleQuadratures(const SubElInfo& subElInfo);
void getSubElementVector(SubElInfo *subElInfo, void getSubElementVector(SubElInfo *subElInfo,
const ElInfo *elInfo,