Skip to content
GitLab
Commit 6caf6e9c authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Simplification of assembling on different meshes.

parent 7cca5f50
Hide whitespace changes
Inline Side-by-side
AMDiS/src/AMDiS_fwd.h
View file @ 6caf6e9c
......@@ -69,6 +69,11 @@ namespace AMDiS {
class Projection;
class PreconditionerScal;
class Quadrature;
class Q00PsiPhi;
class Q0Psi;
class Q10PsiPhi;
class Q01PsiPhi;
class Q1Psi;
class RCNeighbourList;
class RefinementManager;
class RobinBC;
......
AMDiS/src/Assembler.cc
View file @ 6caf6e9c
......@@ -109,19 +109,23 @@ namespace AMDiS {
ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
TEST_EXIT(zeroOrderAssembler)("Not yet implemented for not zero order assembler!\n");
if (secondOrderAssembler)
secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo, mat);
secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
if (firstOrderAssemblerGrdPsi)
firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo, mat);
firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
if (firstOrderAssemblerGrdPhi)
firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo, mat);
if (zeroOrderAssembler)
zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo, mat);
firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
if (zeroOrderAssembler)
zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
ElementMatrix tmpMat(nRow, nRow);
tmpMat = m * mat;
mat = tmpMat;
if (rememberElMat)
userMat += factor * elementMatrix;
}
......
AMDiS/src/FirstOrderAssembler.cc
View file @ 6caf6e9c
......@@ -147,59 +147,6 @@ namespace AMDiS {
}
}
void Stand10::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
FUNCNAME("Stand10::calculateElementMatrix()");
TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
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(smallElInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++)
Lb[iq] *= smallElInfo->getDet();
DimVec<double> val1(dim, DEFAULT_VALUE, 0.0);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
for (int iq = 0; iq < nPoints; iq++) {
val1.set(0.0);
for (int k = 0; k < nCol; k++) {
(*(psi->getGrdPhi(k)))(quadrature->getLambda(iq), grdPsi);
grdPsi *= m[i][k];
val1 += grdPsi;
}
mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * val1) *
(*(phi->getPhi(j)))(quadrature->getLambda(iq));
}
}
}
}
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
......@@ -365,7 +312,6 @@ namespace AMDiS {
}
}
Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
{}
......@@ -404,55 +350,6 @@ namespace AMDiS {
}
}
void Stand01::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
FUNCNAME("Stand01::calculateElementMatrix()");
TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");
DimVec<double> grdPhi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
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(smallElInfo, nPoints, Lb);
for (int iq = 0; iq < nPoints; iq++)
Lb[iq] *= smallElInfo->getDet();
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
for (int iq = 0; iq < nPoints; iq++) {
double val0 = 0.0;
for (int k = 0; k < nCol; k++)
val0 = m[i][k] * (*(psi->getPhi(k)))(quadrature->getLambda(iq));
(*(phi->getGrdPhi(j)))(quadrature->getLambda(iq), grdPhi);
mat[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * val0) * grdPhi);
}
}
}
}
Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{}
......
AMDiS/src/FirstOrderAssembler.h
View file @ 6caf6e9c
#ifndef AMDIS_FIRST_ORDER_ASSEMBLER_H
#define AMDIS_FIRST_ORDER_ASSEMBLER_H
#include "AMDiS_fwd.h"
#include "SubAssembler.h"
namespace AMDiS {
class Q10PsiPhi;
class Q01PsiPhi;
class Q1Psi;
/**
* \ingroup Assembler
*
......@@ -79,13 +76,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat);
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *, ElementVector&);
};
......@@ -106,13 +96,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat);
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo*, ElementVector&) {
ERROR_EXIT("should not be called\n");
......@@ -135,16 +118,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
ERROR_EXIT("CEM not yet\n");
}
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *, ElementVector&);
};
......@@ -164,16 +137,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
ERROR_EXIT("CEM not yet\n");
}
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo*, ElementVector&) {
ERROR_EXIT("should not be called\n");
......@@ -196,16 +159,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
ERROR_EXIT("CEM not yet\n");
}
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo*, ElementVector&);
......@@ -235,16 +188,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
ERROR_EXIT("CEM not yet\n");
}
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo*, ElementVector&) {
ERROR_EXIT("should not be called\n");
......
AMDiS/src/SecondOrderAssembler.cc
View file @ 6caf6e9c
......@@ -70,11 +70,11 @@ namespace AMDiS {
if (!optimized) {
newAssembler = NEW Stand2(op, assembler, quad);
} else {
if (pwConst) {
newAssembler = NEW Pre2(op, assembler, quad);
} else {
newAssembler = NEW Quad2(op, assembler, quad);
}
if (pwConst) {
newAssembler = NEW Pre2(op, assembler, quad);
} else {
newAssembler = NEW Quad2(op, assembler, quad);
}
}
subAssemblers->push_back(newAssembler);
......@@ -321,114 +321,4 @@ namespace AMDiS {
delete [] LALt;
}
void Stand2::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
FUNCNAME("Stand2::calculateElementMatrix()");
TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");
bool p = false;
/*
if (smallElInfo->getLevel() == largeElInfo->getLevel() + 1) {
p = true;
std::cout << "--------------------------------------" << std::endl;
std::cout << "i-th child = " << smallElInfo->getIChild() << std::endl;
std::cout << "Mesh row = " << owner->getRowFESpace()->getMesh() << std::endl;
std::cout << "Mesh col = " << owner->getColFESpace()->getMesh() << std::endl;
std::cout << "Mesh main = " << rowElInfo->getMesh() << std::endl;
std::cout << "Mesh aux = " << colElInfo->getMesh() << std::endl;
}
*/
DimVec<double> grdPsi(dim, NO_INIT);
VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
int nPoints = quadrature->getNumPoints();
DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
for (int iq = 0; iq < nPoints; iq++) {
LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
LALt[iq]->set(0.0);
}
int myRank = omp_get_thread_num();
if (p) {
std::cout << "Terms = " << terms[myRank].size() << std::endl;
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<SecondOrderTerm*>(terms[myRank][i]))
->getLALt(smallElInfo, nPoints, LALt);
}
if (p) {
for (int i = 0; i < nPoints; i++) {
std::cout << "LALt for iq = " << i << std::endl;
LALt[i]->print();
}
}
DimVec<double> val1(dim, DEFAULT_VALUE, 0.0);
if (p)
std::cout << "Multiply LALt with det = " << smallElInfo->getDet() << std::endl;
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= smallElInfo->getDet();
}
/*
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
for (int iq = 0; iq < nPoints; iq++) {
val1.set(0.0);
for (int k = 0; k < nCol; k++) {
(*(phi->getGrdPhi(k)))(quadrature->getLambda(iq), grdPsi);
grdPsi *= m[j][k];
val1 += grdPsi;
}
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
mat[i][j] += quadrature->getWeight(iq) *
(grdPsi * ((*LALt[iq]) * val1));
}
}
}*/
for (int iq = 0; iq < nPoints; iq++) {
for (int i = 0; i < nCol; i++) {
(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
if (p) {
std::cout << "grdPhi[" << i << "]" << std::endl;
grdPhi[i].print();
}
}
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) *
(grdPsi * ((*LALt[iq]) * grdPhi[j]));
}
}
}
for (int iq = 0; iq < nPoints; iq++)
delete LALt[iq];
delete [] LALt;
}
}
AMDiS/src/SecondOrderAssembler.h
View file @ 6caf6e9c
......@@ -62,13 +62,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat);
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *, ElementVector&) {
ERROR_EXIT("should not be called\n");
......@@ -92,16 +85,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
ERROR_EXIT("CEM not yet\n");
}
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *, ElementVector&) {
ERROR_EXIT("should not be called\n");
......@@ -133,16 +116,6 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
///
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
ERROR_EXIT("Not yet implemented!\n");
}
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *, ElementVector&) {
ERROR_EXIT("should not be called\n");
......
AMDiS/src/SubAssembler.h
View file @ 6caf6e9c
......@@ -49,12 +49,6 @@ namespace AMDiS {
virtual void calculateElementMatrix(const ElInfo *elInfo,
ElementMatrix& mat) = 0;
virtual void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat) = 0;
/** \brief
* Calculates the element vector for elInfo and adds it to vec. Memory
* for vec must be provided by the caller.
......
AMDiS/src/ZeroOrderAssembler.cc
View file @ 6caf6e9c
......@@ -137,72 +137,6 @@ namespace AMDiS {
}
}
void StandardZOA::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
FUNCNAME("StandardZOA::calculateElementMatrix()");
TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");
calculateElementMatrix(smallElInfo, mat);
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
ElementMatrix tmpMat(nRow, nRow);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nRow; j++) {
tmpMat[i][j] = 0.0;
for (int k = 0; k < nRow; k++) {
tmpMat[i][j] += m[i][j] * mat[j][i];
}
}
}
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nRow; j++) {
mat[i][j] = tmpMat[i][j];
}
}
/*
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
int myRank = omp_get_thread_num();
int nPoints = quadrature->getNumPoints();
std::vector<double> c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(smallElInfo, nPoints, c);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
for (int iq = 0; iq < nPoints; iq++) {
double val0 = c[iq] * smallElInfo->getDet() * quadrature->getWeight(iq) *
(*(phi->getPhi(i)))(quadrature->getLambda(iq));
double val1 = 0.0;
for (int k = 0; k < nCol; k++)
val1 += m[j][k] * (*(psi->getPhi(k)))(quadrature->getLambda(iq));
val0 *= val1;
mat[i][j] += val0;
}
}
}
*/
}
void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
int nPoints = quadrature->getNumPoints();
......@@ -285,62 +219,6 @@ namespace AMDiS {
}
}
void FastQuadZOA::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
ElementMatrix& mat)
{
FUNCNAME("FastQuadZOA::calculateElementMatrix()");
ERROR_EXIT("CHECK FIRST, IF IT WILL REALLY WORK!\n");
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
if (firstCall) {
#ifdef _OPENMP
#pragma omp critical
#endif
{
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
firstCall = false;
}
}
std::vector<double> c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(rowElInfo, nPoints, c);
}
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= smallElInfo->getDet();
const double *psi = psiFast->getPhi(iq);
const double *phi = phiFast->getPhi(iq);
for (int i = 0; i < nCol; i++) {
for (int j = 0; j < nRow; j++) {
double val = quadrature->getWeight(iq) * c[iq] * psi[i];
double tmpval = 0.0;
for (int k = 0; k < nCol; k++) {
tmpval += m[j][k] * phi[k];
}
val *= tmpval;
mat[j][i] += val;
}
}
}
}