Skip to content
GitLab
Commit 450e3553 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

next step for boundary operators

parent ed9107af
Hide whitespace changes
Inline Side-by-side
dune/amdis/Operator.hpp
View file @ 450e3553
......@@ -123,8 +123,8 @@ namespace AMDiS
/// Calculates the needed quadrature degree for the given term-order \p order.
/// This takes into account the polynomial degree of trial and test functions
/// and the polynomial degree of the coefficient functions
template <class Element>
int getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType = GRD_PHI);
template <class Geometry>
int getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, FirstOrderType firstOrderType = GRD_PHI);
template <class RowView, class ColView>
......@@ -133,12 +133,24 @@ namespace AMDiS
ElementMatrix& elementMatrix,
double* factor = NULL);
template <class Intersection, class RowView, class ColView>
bool getBoundaryElementMatrix(Intersection const& intersection,
RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor = NULL);
template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor = NULL);
template <class Intersection, class RowView>
bool getBoundaryElementVector(Intersection const& intersection,
RowView const& rowView,
ElementVector& elementVector,
double* factor = NULL);
protected:
template <class RowView, class ColView>
void assembleZeroOrderTerms(RowView const& rowView,
......
dune/amdis/Operator.inc.hpp
View file @ 450e3553
......@@ -27,14 +27,11 @@ namespace AMDiS
template <class MeshView>
template <class Element>
int Operator<MeshView>::getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType)
template <class Geometry>
int Operator<MeshView>::getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, FirstOrderType firstOrderType)
{
std::list<OperatorTermType*>* terms = NULL;
auto const& t = element.type();
auto geometry = element.geometry();
switch(order)
{
case 0:
......@@ -86,8 +83,6 @@ namespace AMDiS
ElementMatrix& elementMatrix,
double* factor)
{
AMDIS_FUNCNAME("Operator::getElementMatrix()");
double fac = factor ? *factor : 1.0;
if (fac == 0.0)
......@@ -128,6 +123,27 @@ namespace AMDiS
return true;
}
template <class MeshView>
template <class Intersection, class RowView, class ColView>
bool Operator<MeshView>::getBoundaryElementMatrix(Intersection const& intersection,
RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor)
{
const int dim = Intersection::mydimension;
auto const& element = intersection.inside();
auto geometryInInside = intersection.geometryInInside();
int order = getQuadratureDegree(intersection.type(), intersection.geometry(), 0);
auto const& quad = Dune::QuadratureRules<double, dim-1>::rule(intersection.type(), order);
int facet = intersection.indexInInside();
// TODO:
// - we need an abstraction for quadrature, quad-order, intersection-restriction
}
template <class MeshView>
template <class RowView>
......@@ -172,6 +188,15 @@ namespace AMDiS
}
template <class MeshView>
template <class Intersection, class RowView>
bool Operator<MeshView>::getBoundaryElementVector(Intersection const& intersection,
RowView const& rowView,
ElementVector& elementVector,
double* factor)
{
}
template <class MeshView>
template <class RowView, class ColView>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
......@@ -188,7 +213,7 @@ namespace AMDiS
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(element, 0);
int order = getQuadratureDegree(element.type(), geometry, 0);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
......@@ -237,7 +262,7 @@ namespace AMDiS
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(element, 0);
int order = getQuadratureDegree(element.type(), geometry, 0);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
......@@ -278,7 +303,7 @@ namespace AMDiS
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(element, 1, GRD_PHI);
int order = getQuadratureDegree(element.type(), geometry, 1, GRD_PHI);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPhi)
......@@ -336,7 +361,7 @@ namespace AMDiS
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(element, 1, GRD_PSI);
int order = getQuadratureDegree(element.type(), geometry, 1, GRD_PSI);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPsi)
......@@ -393,7 +418,7 @@ namespace AMDiS
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(element, 1, GRD_PSI);
int order = getQuadratureDegree(element.type(), geometry, 1, GRD_PSI);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPsi)
......@@ -444,7 +469,7 @@ namespace AMDiS
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
// auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(element, 2);
int order = getQuadratureDegree(element.type(), geometry, 2);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : secondOrder)
......
dune/amdis/ProblemStat.hpp
View file @ 450e3553
......@@ -34,6 +34,7 @@
namespace AMDiS
{
struct BoundaryType { int b; };
template <class Traits>
class ProblemStatSeq
......@@ -76,6 +77,32 @@ namespace AMDiS
using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
typename SystemVectorType::MultiVector>;
struct ScaledOperator
{
ScaledOperator(std::shared_ptr<OperatorType> const& op, double* factor = nullptr, double* estFactor = nullptr)
: op(op)
, factor(factor)
, estFactor(estFactor)
{}
std::shared_ptr<OperatorType> op;
double* factor;
double* estFactor;
};
struct ScaledBoundaryOperator
{
ScaledBoundaryOperator(std::shared_ptr<OperatorType> const& op, BoundaryType b, double* factor = nullptr, double* estFactor = nullptr)
: op(op)
, b(b)
, factor(factor)
, estFactor(estFactor)
{}
std::shared_ptr<OperatorType> op;
BoundaryType b;
double* factor;
double* estFactor;
};
public:
/**
* \brief Constructor. Takes the name of the problem that is used to
......@@ -106,21 +133,17 @@ namespace AMDiS
* MESH[0]->global refinements: nr of initial global refinements
**/
void initialize(Flag initFlag,
Self* adoptProblem = NULL,
Self* adoptProblem = nullptr,
Flag adoptFlag = INIT_NOTHING);
/// Adds an operator to \ref A.
/** @{ */
void addMatrixOperator(OperatorType& op,
int i, int j,
double* factor = NULL,
double* estFactor = NULL);
void addMatrixOperator(OperatorType const& op, int i, int j,
double* factor = nullptr, double* estFactor = nullptr);
void addMatrixOperator(std::shared_ptr<OperatorType> op,
int i, int j,
double* factor = NULL,
double* estFactor = NULL);
void addBoundaryMatrixOperator(BoundaryType b, OperatorType const& op, int i, int j,
double* factor = nullptr, double* estFactor = nullptr);
void addMatrixOperator(std::map< std::pair<int,int>,
std::shared_ptr<OperatorType> > ops);
......@@ -131,15 +154,11 @@ namespace AMDiS
/// Adds an operator to \ref rhs.
/** @{ */
void addVectorOperator(OperatorType& op,
int i,
double* factor = NULL,
double* estFactor = NULL);
void addVectorOperator(OperatorType const& op, int i,
double* factor = nullptr, double* estFactor = nullptr);
void addVectorOperator(std::shared_ptr<OperatorType> op,
int i,
double* factor = NULL,
double* estFactor = NULL);
void addBoundaryVectorOperator(BoundaryType b, OperatorType const& op, int i,
double* factor = nullptr, double* estFactor = nullptr);
void addVectorOperator(std::map< int, std::shared_ptr<OperatorType> > ops);
void addVectorOperator(std::map< int, std::pair<std::shared_ptr<OperatorType>, bool> > ops);
......@@ -314,14 +333,14 @@ namespace AMDiS
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
std::list<std::shared_ptr<OperatorType>>& operators,
std::list<double*> const& factors);
std::list<ScaledOperator>& operators,
std::list<ScaledBoundaryOperator>& boundary_operators);
template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
std::list<std::shared_ptr<OperatorType>>& operators,
std::list<double*> const& factors);
std::list<ScaledOperator>& operators,
std::list<ScaledBoundaryOperator>& boundary_operators);
template <class Matrix, class RowIndexSet, class ColIndexSet>
......@@ -374,29 +393,30 @@ namespace AMDiS
template <class T>
using MatrixEntries = std::map< std::pair<int,int>, std::list<T> >;
using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
template <class T>
using VectorEntries = std::map< int, std::list<T> >;
using VectorEntries = Dune::FieldVector<T, nComponents>;
using BoundaryCondition = std::list<std::shared_ptr<DirichletBC<WorldVector>>>;
/// A map (i,j) -> list<DirichleBC> string a boundary conditions for
/// each matrix block
MatrixEntries<std::shared_ptr<DirichletBC<WorldVector>>> dirichletBc;
MatrixEntries<BoundaryCondition> dirichletBc;
/// A map (i,j) -> list<OperatorType> string the differential operators for
/// each matrix block
MatrixEntries<std::shared_ptr<OperatorType>> matrixOperators;
MatrixEntries<double*> matrixFactors;
std::map< std::pair<int,int>, bool > matrixAssembled; // if false, do reassemble
std::map< std::pair<int,int>, bool > matrixChanging; // if true, or vectorAssembled false, do reassemble
MatrixEntries<std::list<ScaledOperator>> matrixOperators;
MatrixEntries<std::list<ScaledBoundaryOperator>> matrixBoundaryOperators;
MatrixEntries<bool> matrixAssembled; // if false, do reassemble
MatrixEntries<bool> matrixChanging; // if true, or vectorAssembled false, do reassemble
/// A map (i) -> list<OperatorType> string the differential operators for
/// each rhs block
VectorEntries<std::shared_ptr<OperatorType>> vectorOperators;
VectorEntries<double*> vectorFactors;
std::map< int, bool > vectorAssembled; // if false, do reassemble
std::map< int, bool > vectorChanging; // if true, or vectorAssembled false, do reassemble
VectorEntries<std::list<ScaledOperator>> vectorOperators;
VectorEntries<std::list<ScaledBoundaryOperator>> vectorBoundaryOperators;
VectorEntries<bool> vectorAssembled; // if false, do reassemble
VectorEntries<bool> vectorChanging; // if true, or vectorAssembled false, do reassemble
template <int R, int C>
struct MatrixData
......@@ -404,12 +424,13 @@ namespace AMDiS
using DOFMatrixType =
std::tuple_element_t<C, std::tuple_element_t<R, typename SystemMatrixType::DOFMatrices>>;
DOFMatrixType& matrix;
std::list<std::shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors;
bool assemble = true;
DOFMatrixType& matrix;
std::list<ScaledOperator>& operators;
std::list<ScaledBoundaryOperator>& boundary_operators;
bool assemble = true;
std::pair<int,int> row_col = {R, C};
int const row = R;
int const col = C;
};
template <int I>
......@@ -417,12 +438,12 @@ namespace AMDiS
{
using DOFVectorType = std::tuple_element_t<I, typename SystemVectorType::DOFVectors>;
DOFVectorType& vector;
std::list<std::shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors;
bool assemble = true;
DOFVectorType& vector;
std::list<ScaledOperator>& operators;
std::list<ScaledBoundaryOperator>& boundary_operators;
bool assemble = true;
int row = I;
int const index = I;
};
};
......
dune/amdis/ProblemStat.inc.hpp
View file @ 450e3553
......@@ -107,83 +107,75 @@ namespace AMDiS
template <class Traits>
void ProblemStatSeq<Traits>::addMatrixOperator(std::shared_ptr<OperatorType> op,
int i, int j,
double* factor,
double* estFactor)
void ProblemStatSeq<Traits>::
addMatrixOperator(OperatorType const& op, int i, int j, double* factor, double* estFactor)
{
// TODO: currently the factors are ignored
assert( estFactor == NULL );
matrixOperators[{i,j}].push_back(op);
matrixFactors[{i,j}].push_back(factor);
matrixChanging[{i,j}] = true;
matrixOperators[i][j].emplace_back(std::make_shared<OperatorType>(op), factor, estFactor);
matrixChanging[i][j] = true;
}
template <class Traits>
void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op,
int i, int j,
double* factor,
double* estFactor)
void ProblemStatSeq<Traits>::
addBoundaryMatrixOperator(BoundaryType b, OperatorType const& op, int i, int j, double* factor, double* estFactor)
{
addMatrixOperator(std::make_shared<OperatorType>(op), i, j, factor, estFactor);
matrixBoundaryOperators[i][j].emplace_back(std::make_shared<OperatorType>(op), b, factor, estFactor);
matrixChanging[i][j] = true;
}
template <class Traits>
void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::shared_ptr<OperatorType> > ops)
{
for (auto op : ops)
addMatrixOperator(op.second, op.first.first, op.first.second);
for (auto op : ops){
int r = op.first.first;
int c = op.first.second;
matrixOperators[r][c].emplace_back(op.second);
matrixChanging[r][c] = true;
}
}
template <class Traits>
void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::pair<std::shared_ptr<OperatorType>,bool> > ops)
{
for (auto op : ops) {
auto row_col = op.first;
matrixOperators[row_col].push_back(op.second.first);
matrixFactors[row_col].push_back(NULL);
matrixChanging[row_col] = matrixChanging[row_col] || op.second.second;
int r = op.first.first;
int c = op.first.second;
matrixOperators[r][c].emplace_back(op.second.first);
matrixChanging[r][c] = matrixChanging[r][c] || op.second.second;
}
}
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(std::shared_ptr<OperatorType> op,
int i,
double* factor,
double* estFactor)
void ProblemStatSeq<Traits>::
addVectorOperator(OperatorType const& op, int i, double* factor, double* estFactor)
{
// TODO: currently the est-factors are ignored
assert( estFactor == NULL );
vectorOperators[i].push_back(op);
vectorFactors[i].push_back(factor);
vectorOperators[i].emplace_back(std::make_shared<OperatorType>(op), factor, estFactor);
vectorChanging[i] = true;
}
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op,
int i,
double* factor,
double* estFactor)
void ProblemStatSeq<Traits>::
addBoundaryVectorOperator(BoundaryType b, OperatorType const& op, int i, double* factor, double* estFactor)
{
addVectorOperator(std::make_shared<OperatorType>(op), i, factor, estFactor);
vectorBoundaryOperators[i].emplace_back(std::make_shared<OperatorType>(op), b, factor, estFactor);
vectorChanging[i] = true;
}
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::shared_ptr<OperatorType> > ops)
{
for (auto op : ops)
addVectorOperator(op.second, op.first);
for (auto op : ops) {
vectorOperators[op.first].emplace_back(op.second);
vectorChanging[op.first] = true;
}
}
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::pair<std::shared_ptr<OperatorType>, bool> > ops)
{
for (auto op : ops) {
vectorOperators[op.first].push_back(op.second.first);
vectorFactors[op.first].push_back(NULL);
vectorOperators[op.first].emplace_back(op.second.first);
vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
}
}
......@@ -191,9 +183,8 @@ namespace AMDiS
/// Adds a Dirichlet boundary condition
template <class Traits>
template <class Predicate, class Values>
void ProblemStatSeq<Traits>::addDirichletBC(Predicate const& predicate,
int row, int col,
Values const& values)
void ProblemStatSeq<Traits>::
addDirichletBC(Predicate const& predicate, int row, int col, Values const& values)
{
static_assert( Concepts::Callable<Predicate, WorldVector>,
"Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
......@@ -207,14 +198,13 @@ namespace AMDiS
// boundaryConditionSet = true;
using BcType = DirichletBC<WorldVector>;
dirichletBc[{row, col}].emplace_back( new BcType{predicate, values} );
dirichletBc[row][col].emplace_back( new BcType{predicate, values} );
}
template <class Traits>
void ProblemStatSeq<Traits>::solve(AdaptInfo& adaptInfo,
bool createMatrixData,
bool storeMatrixData)
void ProblemStatSeq<Traits>::
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
Timer t;
......@@ -253,8 +243,8 @@ namespace AMDiS
template <class Traits>
void ProblemStatSeq<Traits>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
bool asmMatrix_, bool asmVector_)
void ProblemStatSeq<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix_, bool asmVector_)
{
Timer t;
......@@ -289,45 +279,35 @@ namespace AMDiS
auto& rhs_vec = (*rhs)[_r];
auto row_col = std::make_pair(R, C);
bool assembleMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
bool assembleMatrix = asmMatrix_ && (!matrixAssembled[R][C] || matrixChanging[R][C]);
bool assembleVector = asmVector && R == C;
int r_ = 0, c_ = 0;
if (assembleMatrix) {
// init boundary condition
for (auto bc_list : dirichletBc) {
std::tie(r_, c_) = bc_list.first;
if (r_ == R) {
for (auto bc : bc_list.second)
bc->init(c_ == C, dofmatrix, solution_vec, rhs_vec);
}
}
for (int c = 0; c < nComponents; ++c)
for (auto bc : dirichletBc[R][c])
bc->init(c == C, dofmatrix, solution_vec, rhs_vec);
}
auto mat = MatrixData{dofmatrix, matrixOperators[row_col], matrixFactors[row_col], assembleMatrix};
auto vec = VectorData{rhs_vec, vectorOperators[R], vectorFactors[R], assembleVector};
auto mat = MatrixData{dofmatrix, matrixOperators[R][C], matrixBoundaryOperators[R][C], assembleMatrix};
auto vec = VectorData{rhs_vec, vectorOperators[R], vectorBoundaryOperators[R], assembleVector};
// assemble the DOFMatrix block and the corresponding rhs vector, of r==c
dofmatrix.init(assembleMatrix);
this->assemble(mat, vec, this->leafGridView());
// TODO: assemble boundary conditions
dofmatrix.finish();
if (assembleMatrix)
matrixAssembled[row_col] = true;
matrixAssembled[R][C] = true;
if (assembleVector)
vectorAssembled[R] = true;
if (assembleMatrix) {
// finish boundary condition
for (auto bc_list : dirichletBc) {
std::tie(r_, c_) = bc_list.first;
if (r_ == R) {
for (auto bc : bc_list.second)
bc->finish(c_ == C, dofmatrix, solution_vec, rhs_vec);
}
}
for (int c = 0; c < nComponents; ++c)
for (auto bc : dirichletBc[R][c])
bc->finish(c == C, dofmatrix, solution_vec, rhs_vec);
nnz += dofmatrix.getNnz();
}
......@@ -351,8 +331,7 @@ namespace AMDiS
template <class Traits>
template <class LhsData, class RhsData, class GV>
void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs,
GV const& gridView)
void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs, GV const& gridView)
{
auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
auto const& colFeSpace = lhs.matrix.getColFeSpace();
......@@ -362,10 +341,10 @@ namespace AMDiS
return; // nothing to do
for (auto op : lhs.operators)
op->init(rowFeSpace, colFeSpace);
for (auto op : rhs.operators)
op->init(rowFeSpace, colFeSpace);
for (auto scaledOp : lhs.operators)
scaledOp.op->init(rowFeSpace, colFeSpace);
for (auto scaledOp : rhs.operators)
scaledOp.op->init(rowFeSpace, colFeSpace);
auto rowLocalView = rowFeSpace.localView();
auto rowIndexSet = rowFeSpace.localIndexSet();
......@@ -383,7 +362,7 @@ namespace AMDiS
if (lhs.assemble) {
ElementMatrix elementMatrix;
bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
lhs.operators, lhs.factors);
lhs.operators, lhs.boundary_operators);
if (add)
addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
}
......@@ -391,7 +370,7 @@ namespace AMDiS
if (rhs.assemble) {
ElementVector elementVector;
bool add = getElementVector(rowLocalView, elementVector,
rhs.operators, rhs.factors);
rhs.operators, rhs.boundary_operators);
if (add)
addElementVector(rhs.vector, rowIndexSet, elementVector);
}
......@@ -404,45 +383,44 @@ namespace AMDiS
bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
ColView const& colLocalView,