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

* Splitted Assember.h and Assembler.cc in several files

* First support for independent meshes for different components
parent ec239f4e
......@@ -88,6 +88,10 @@ $(SOURCE_DIR)/Serializable.h $(SOURCE_DIR)/BallProject.h $(SOURCE_DIR)/CylinderP
$(SOURCE_DIR)/MacroReader.h $(SOURCE_DIR)/MacroReader.cc \
$(SOURCE_DIR)/ValueReader.h $(SOURCE_DIR)/ValueReader.cc \
$(SOURCE_DIR)/Projection.h $(SOURCE_DIR)/Projection.cc \
$(SOURCE_DIR)/SubAssembler.h $(SOURCE_DIR)/SubAssembler.cc \
$(SOURCE_DIR)/ZeroOrderAssembler.h $(SOURCE_DIR)/ZeroOrderAssembler.cc \
$(SOURCE_DIR)/FirstOrderAssembler.h $(SOURCE_DIR)/FirstOrderAssembler.cc \
$(SOURCE_DIR)/SecondOrderAssembler.h $(SOURCE_DIR)/SecondOrderAssembler.cc \
$(SOURCE_DIR)/Assembler.h $(SOURCE_DIR)/Assembler.cc \
$(SOURCE_DIR)/AdaptInfo.h $(SOURCE_DIR)/AdaptInfo.cc \
$(SOURCE_DIR)/Marker.h $(SOURCE_DIR)/Marker.cc \
......@@ -221,6 +225,8 @@ $(SOURCE_DIR)/ElementInfo.h \
$(SOURCE_DIR)/VertexInfo.h \
$(SOURCE_DIR)/PeriodicInfo.h \
$(SOURCE_DIR)/OpenMP.h \
$(SOURCE_DIR)/ScalableQuadrature.h $(SOURCE_DIR)/ScalableQuadrature.cc \
$(SOURCE_DIR)/SubElInfo.h $(SOURCE_DIR)/SubElInfo.cc \
$(SOURCE_DIR)/parareal/ProblemBase.h \
$(SOURCE_DIR)/parareal/AdaptParaReal.h $(SOURCE_DIR)/parareal/AdaptParaReal.cc
......@@ -245,9 +251,6 @@ $(COMPOSITE_SOURCE_DIR)/ElementLevelSet.cc \
$(COMPOSITE_SOURCE_DIR)/CompositeFEMOperator.h \
$(COMPOSITE_SOURCE_DIR)/CompositeFEMOperator.cc \
$(COMPOSITE_SOURCE_DIR)/SubPolytope.h $(COMPOSITE_SOURCE_DIR)/SubPolytope.cc \
$(COMPOSITE_SOURCE_DIR)/ScalableQuadrature.h \
$(COMPOSITE_SOURCE_DIR)/ScalableQuadrature.cc \
$(COMPOSITE_SOURCE_DIR)/SubElInfo.h $(COMPOSITE_SOURCE_DIR)/SubElInfo.cc \
$(COMPOSITE_SOURCE_DIR)/SubElementAssembler.h \
$(COMPOSITE_SOURCE_DIR)/SubElementAssembler.cc \
$(COMPOSITE_SOURCE_DIR)/TranslateLsFct.h
......
This diff is collapsed.
......@@ -16,11 +16,9 @@ namespace AMDiS {
// false - no
////////////////////////////////////////////////////////////////////////////
int zeroCounter;
for (int i = 0; i < numIntPoints; i++) {
zeroCounter = 0;
for (int j = 0; j < dim+1; j++) {
int zeroCounter = 0;
for (int j = 0; j < dim + 1; j++) {
if (fabs((*intPoints)[i][j]) <= 1.e-15 ) {
zeroCounter++;
}
......@@ -119,7 +117,7 @@ namespace AMDiS {
TEST_EXIT(dim == 1 && numIntPoints == 1)("invalid call of this routine\n");
VectorOfFixVecs<DimVec<double> > *subElVertices =
NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
NEW VectorOfFixVecs<DimVec<double> >(dim, dim + 1, NO_INIT);
DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
/**
......@@ -132,8 +130,7 @@ namespace AMDiS {
* subelement.
*/
vertex[1] = 0.0;
}
else {
} else {
/**
* The vertex in element with barycentric coordinates (0,1) is in
* subelement.
......@@ -189,7 +186,7 @@ namespace AMDiS {
("invalid call of this routine\n");
VectorOfFixVecs<DimVec<double> >*subElVertices =
NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
NEW VectorOfFixVecs<DimVec<double> >(dim, dim + 1, NO_INIT);
DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
/**
......@@ -197,7 +194,7 @@ namespace AMDiS {
* a subelement of element.
*/
for (int i = 0; i < numIntPoints; i++) {
for (int j = 0; j < dim+1; j++) {
for (int j = 0; j < dim + 1; j++) {
if ( fabs((*intPoints)[i][j]) <= 1.e-15 ) {
vertex[j] = 0.0;
};
......@@ -232,7 +229,7 @@ namespace AMDiS {
int SubPolytope::getIndexSecondFaceIntPoint0(int indexFirstFace, int dim)
{
for (int i = 0; i < dim+1; i++) {
for (int i = 0; i < dim + 1; i++) {
if ( fabs((*intPoints)[0][i]) <= 1.e-15 && i != indexFirstFace ) {
return i;
}
......@@ -276,7 +273,7 @@ namespace AMDiS {
("invalid index for vertex of a tetrahedron");
VectorOfFixVecs<DimVec<double> > *subElVertices =
NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
NEW VectorOfFixVecs<DimVec<double> >(dim, dim + 1, NO_INIT);
DimVec<double> vertexA(dim, DEFAULT_VALUE, 0.0);
DimVec<double> vertexB(dim, DEFAULT_VALUE, 0.0);
......@@ -319,7 +316,7 @@ namespace AMDiS {
// Get the edges including the intersection points.
for (int i = 0; i < numIntPoints; i++) {
int k = 0;
for (int j = 0; j < dim+1; j++) {
for (int j = 0; j < dim + 1; j++) {
if (fabs((*intPoints)[i][j]) > 1.e-15 ) {
indexEdge[k] = j;
k++;
......@@ -332,7 +329,7 @@ namespace AMDiS {
// Get the vertex of element adjacent with indexElVertInPol1 whose
// common edge with indexElVertInPol1 doesn't contain an
// intersection point, and store it in indexElVertInPol2.
for (int i = 0; i < dim+1; i++) {
for (int i = 0; i < dim + 1; i++) {
if (intPointOnEdge[indexElVertInPol1][i] == false &&
i != indexElVertInPol1 ) {
indexElVertInPol2 = i;
......
This diff is collapsed.
This diff is collapsed.
......@@ -463,6 +463,36 @@ namespace AMDiS {
addElementMatrix(factor, *elementMatrix, bound);
}
void DOFMatrix::assemble(double factor, ElInfo *rowElInfo, ElInfo *colElInfo,
const BoundaryType *bound, Operator *op)
{
FUNCNAME("DOFMatrix::assemble()");
if (!op && operators.size() == 0) {
return;
}
Operator *operat = op ? op : operators[0];
operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, rowElInfo);
if (op) {
op->getElementMatrix(rowElInfo, colElInfo, elementMatrix);
} else {
std::vector<Operator*>::iterator it;
std::vector<double*>::iterator factorIt;
for (it = operators.begin(), factorIt = operatorFactor.begin();
it != operators.end();
++it, ++factorIt) {
(*it)->getElementMatrix(rowElInfo, colElInfo,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
}
addElementMatrix(factor, *elementMatrix, bound);
}
Flag DOFMatrix::getAssembleFlag()
{
Flag fillFlag(0);
......
......@@ -364,6 +364,10 @@ namespace AMDiS {
const BoundaryType *bound, Operator *op = NULL);
void assemble(double factor, ElInfo *rowElInfo, ElInfo *colElInfo,
const BoundaryType *bound, Operator *op = NULL);
/** \brief
* Adds an element matrix to \ref matrix
*/
......
......@@ -15,6 +15,7 @@
#include "ElementVector.h"
#include "Assembler.h"
#include "OpenMP.h"
#include "Operator.h"
namespace AMDiS {
......
#include <vector>
#include "Assembler.h"
#include "FirstOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"
#include "ElementMatrix.h"
#include "OpenMP.h"
namespace AMDiS {
std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPsi;
std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPhi;
std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPsi;
FirstOrderAssembler::FirstOrderAssembler(Operator *op,
Assembler *assembler,
Quadrature *quad,
bool optimized,
FirstOrderType type)
: SubAssembler(op, assembler, quad, 1, optimized, type)
{}
FirstOrderAssembler*
FirstOrderAssembler::getSubAssembler(Operator* op,
Assembler *assembler,
Quadrature *quad,
FirstOrderType type,
bool optimized)
{
std::vector<SubAssembler*> *subAssemblers =
optimized ?
(type == GRD_PSI ?
&optimizedSubAssemblersGrdPsi :
&optimizedSubAssemblersGrdPhi) :
(type == GRD_PSI ?
&standardSubAssemblersGrdPsi :
&standardSubAssemblersGrdPhi);
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*> opTerms
= (type == GRD_PSI) ?
op->firstOrderGrdPsi[myRank] :
op->firstOrderGrdPhi[myRank];
// check if a assembler is needed at all
if (opTerms.size() == 0) {
return NULL;
}
sort(opTerms.begin(), opTerms.end());
FirstOrderAssembler *newAssembler;
// check if a new assembler is needed
for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
sort(assTerms.begin(), assTerms.end());
if ((opTerms == assTerms) &&
((*subAssemblers)[i]->getQuadrature() == quad)) {
return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
}
}
// check if all terms are pw_const
bool pwConst = true;
for (int i = 0; i < static_cast<int>( opTerms.size()); i++) {
if (!(opTerms[i])->isPWConst()) {
pwConst = false;
break;
}
}
// create new assembler
if (!optimized) {
newAssembler =
(type == GRD_PSI) ?
dynamic_cast<FirstOrderAssembler*>(NEW Stand10(op, assembler, quad)) :
dynamic_cast<FirstOrderAssembler*>(NEW Stand01(op, assembler, quad));
} else {
if (pwConst) {
newAssembler =
(type == GRD_PSI) ?
dynamic_cast<FirstOrderAssembler*>(NEW Pre10(op, assembler, quad)) :
dynamic_cast<FirstOrderAssembler*>(NEW Pre01(op, assembler, quad));
} else {
newAssembler =
(type == GRD_PSI) ?
dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) :
dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad));
}
}
subAssemblers->push_back(newAssembler);
return newAssembler;
};
Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
{}
void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
double *phival = GET_MEMORY(double, nCol);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
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 < 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];
}
}
}
FREE_MEMORY(phival, double, nCol);
}
Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
{
}
void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
const double *phi;
if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
firstCall = false;
}
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
int myRank = omp_get_thread_num();
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
phi = phiFast->getPhi(iq);
grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++)
(*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
}
}
}
Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
{
}
void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > Lb(dim, 1, NO_INIT);
const int *k;
const double *values;
if (firstCall) {
q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
quadrature);
firstCall = false;
}
const int **nEntries = q10->getNumberEntries();
int myRank = omp_get_thread_num();
Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
}
Lb[0] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
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]];
}
(*mat)[i][j] += val;
}
}
}
Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
{}
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();
VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < 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));
for (int j = 0; j < nCol; j++)
(*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
}
}
}
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
}
}
}
Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{
}
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > *grdPhi;
if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
firstCall = false;
}
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
const double *psi = psiFast->getPhi(iq);
grdPhi = phiFast->getGradient(iq);
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];
}
}
}
void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
firstCall = false;
}
int nPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
int myRank = omp_get_thread_num();
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) {
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
}
}
}
Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{
}
void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
const int *l;
const double *values;
if (firstCall) {
q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
quadrature);
firstCall = false;
}
const int **nEntries = q01->getNumberEntries();
int myRank = omp_get_thread_num();
Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);