Skip to content
Snippets Groups Projects
Commit 02403d46 authored by graeser's avatar graeser
Browse files

Extract reusable functionality to header

parent cd9137e6
No related branches found
No related tags found
No related merge requests found
......@@ -41,200 +41,9 @@
#include <dune/fu-tutorial/referenceelementutility.hh>
template<class Grid>
auto createCubeGrid()
{
static const int dim = Grid::dimension;
auto lowerLeft = Dune::FieldVector<double,dim>(0);
auto upperRight = Dune::FieldVector<double,dim>(1);
auto elementsPerDirection = std::array<unsigned int,dim>();
for(auto& e : elementsPerDirection)
e = 2;
return Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, elementsPerDirection);
}
template<class Grid>
auto createSimplexGrid()
{
static const int dim = Grid::dimension;
auto lowerLeft = Dune::FieldVector<double,dim>(0);
auto upperRight = Dune::FieldVector<double,dim>(1);
auto elementsPerDirection = std::array<unsigned int,dim>();
for(auto& e : elementsPerDirection)
e = 2;
return Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, elementsPerDirection);
}
template<class GridView>
void writeGridView(const GridView& gridView, std::string postFix)
{
Dune::VTKWriter<GridView> vtkWriter(gridView);
vtkWriter.write(std::string("04-gridviews-")+postFix);
}
// Compute local per element matrix
template <class Element, class LocalFiniteElement, class MatrixType>
void assembleLocalStiffnessMatrix(const Element& element, const LocalFiniteElement& localFiniteElement, MatrixType& elementMatrix)
{
auto geometryType = element.type();
const auto& geometry = element.geometry();
const int dim = Element::dimension;
std::size_t localSize = localFiniteElement.localBasis().size();
// Resize and initialize matrix
elementMatrix.setSize(localSize, localSize);
elementMatrix = 0; // fills the entire matrix with zeroes
// Get a quadrature rule
int order = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
// Loop over all quadrature points
for (size_t k=0; k < quad.size(); ++k)
{
const auto& quadPos = quad[k].position();
const auto& weight = quad[k].weight();
// The transposed inverse Jacobian of the map from the reference element to the element
const auto& jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
const auto& integrationElement = geometry.integrationElement(quadPos);
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
// Compute the shape function gradients on the real basis functions
std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
// Compute the actual matrix entries
for (size_t i=0; i<elementMatrix.N(); i++)
for (size_t j=0; j<elementMatrix.M(); j++ )
elementMatrix[i][j] += ( gradients[i] * gradients[j] ) * weight * integrationElement;
}
}
// Compute local per element right hand side
template <class Element, class LocalFiniteElement, class LocalRHSVector, class RHSFunction>
void assembleLocalRHS(
const Element& element,
const LocalFiniteElement& localFiniteElement,
LocalRHSVector& localRhs,
const RHSFunction& rhsFunction)
{
auto geometryType = element.type();
const auto& geometry = element.geometry();
const int dim = Element::dimension;
std::size_t localSize = localFiniteElement.localBasis().size();
// Resize and initialize vector
localRhs.resize(localSize);
localRhs = 0;
// Get a quadrature rule
int order = 2*(dim*localFiniteElement.localBasis().order());
const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
// Loop over all quadrature points
for (size_t k=0; k < quad.size(); ++k)
{
const auto& quadPos = quad[k].position();
const auto& weight = quad[k].weight();
// Quadrature position in global world coordinates
const auto& globalQuadPos = geometry.global(quadPos);
// The multiplicative factor in the integral transformation formula
const auto& integrationElement = geometry.integrationElement(quadPos);
double functionValue = rhsFunction(globalQuadPos);
// Evaluate all shape function values at this point
std::vector<Dune::FieldVector<double,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
// Actually compute the vector entries
for (size_t i=0; i<localSize; i++)
localRhs[i] += shapeFunctionValues[i] * functionValue * weight * integrationElement;
}
}
template<class GridView, class RHSFunction>
void assemblePoissonProblem(
const GridView& gridView,
Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>& matrix,
Dune::BlockVector<Dune::FieldVector<double,1>>& rhs,
const RHSFunction& rhsFunction)
{
static const int dim = GridView::dimension;
const auto& indexSet = gridView.indexSet();
std::size_t size = indexSet.size(dim);
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
using ElementMatrix = Dune::Matrix<Dune::FieldMatrix<double,1,1>>;
using ElementRhs = Dune::BlockVector<Dune::FieldVector<double,1>>;
rhs.resize(size);
rhs = 0;
Dune::ImplicitMatrixBuilder<Matrix> matrixBuilder(matrix, size, size, 10, 0.1);
Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache;
ElementMatrix elementMatrix;
ElementRhs elementRhs;
for(const auto& e: Dune::elements(gridView))
{
const auto& localFiniteElement = feCache.get(e.type());
std::size_t localSize = localFiniteElement.localBasis().size();
assembleLocalStiffnessMatrix(e, localFiniteElement, elementMatrix);
assembleLocalRHS(e, localFiniteElement, elementRhs, rhsFunction);
for(std::size_t i=0; i<localSize; ++i)
{
std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim);
for(std::size_t j=0; j<localSize; ++j)
{
std::size_t globalJ = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(j).subEntity(), dim);
matrixBuilder[globalI][globalJ] += elementMatrix[i][j];
}
rhs[globalI] = elementRhs[i];
}
}
matrix.compress();
}
template<class GridView, class BitVector>
void computeBoundaryVertices(const GridView& gridView, BitVector& isBoundary)
{
using namespace Dune;
using namespace Dune::FuTutorial;
const auto& indexSet = gridView.indexSet();
for(const auto& element: elements(gridView))
for(const auto& intersection: intersections(gridView, element))
if (intersection.boundary())
for(const auto& vertex: subVertices(referenceElement(element), insideFacet(intersection)))
isBoundary[subIndex(indexSet, element, vertex)] = true;
}
// included dune-fu-tutorial headers
#include "04-gridviews.hh"
#include "05-poisson-problem.hh"
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUTUTORIAL_SRC_05POISSON_PROBLEM_HH
#define DUNE_FUTUTORIAL_SRC_05POISSON_PROBLEM_HH
// included standard library headers
#include <iostream>
#include <array>
// included dune-common headers
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/stringutility.hh>
// included dune-geometry headers
#include <dune/geometry/quadraturerules.hh>
// included dune-grid headers
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/uggrid.hh>
// included dune-istl headers
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
// included dune-localfunctions headers
#include <dune/localfunctions/lagrange/pqkfactory.hh>
// included dune-fu-tutorial headers
#include <dune/fu-tutorial/referenceelementutility.hh>
// Compute local per element matrix
template <class Element, class LocalFiniteElement, class MatrixType>
void assembleLocalStiffnessMatrix(const Element& element, const LocalFiniteElement& localFiniteElement, MatrixType& elementMatrix)
{
auto geometryType = element.type();
const auto& geometry = element.geometry();
const int dim = Element::dimension;
std::size_t localSize = localFiniteElement.localBasis().size();
// Resize and initialize matrix
elementMatrix.setSize(localSize, localSize);
elementMatrix = 0; // fills the entire matrix with zeroes
// Get a quadrature rule
int order = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
// Loop over all quadrature points
for (size_t k=0; k < quad.size(); ++k)
{
const auto& quadPos = quad[k].position();
const auto& weight = quad[k].weight();
// The transposed inverse Jacobian of the map from the reference element to the element
const auto& jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
const auto& integrationElement = geometry.integrationElement(quadPos);
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
// Compute the shape function gradients on the real basis functions
std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
// Compute the actual matrix entries
for (size_t i=0; i<elementMatrix.N(); i++)
for (size_t j=0; j<elementMatrix.M(); j++ )
elementMatrix[i][j] += ( gradients[i] * gradients[j] ) * weight * integrationElement;
}
}
// Compute local per element right hand side
template <class Element, class LocalFiniteElement, class LocalRHSVector, class RHSFunction>
void assembleLocalRHS(
const Element& element,
const LocalFiniteElement& localFiniteElement,
LocalRHSVector& localRhs,
const RHSFunction& rhsFunction)
{
auto geometryType = element.type();
const auto& geometry = element.geometry();
const int dim = Element::dimension;
std::size_t localSize = localFiniteElement.localBasis().size();
// Resize and initialize vector
localRhs.resize(localSize);
localRhs = 0;
// Get a quadrature rule
int order = 2*(dim*localFiniteElement.localBasis().order());
const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
// Loop over all quadrature points
for (size_t k=0; k < quad.size(); ++k)
{
const auto& quadPos = quad[k].position();
const auto& weight = quad[k].weight();
// Quadrature position in global world coordinates
const auto& globalQuadPos = geometry.global(quadPos);
// The multiplicative factor in the integral transformation formula
const auto& integrationElement = geometry.integrationElement(quadPos);
double functionValue = rhsFunction(globalQuadPos);
// Evaluate all shape function values at this point
std::vector<Dune::FieldVector<double,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
// Actually compute the vector entries
for (size_t i=0; i<localSize; i++)
localRhs[i] += shapeFunctionValues[i] * functionValue * weight * integrationElement;
}
}
template<class GridView, class RHSFunction>
void assemblePoissonProblem(
const GridView& gridView,
Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>& matrix,
Dune::BlockVector<Dune::FieldVector<double,1>>& rhs,
const RHSFunction& rhsFunction)
{
static const int dim = GridView::dimension;
const auto& indexSet = gridView.indexSet();
std::size_t size = indexSet.size(dim);
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
using ElementMatrix = Dune::Matrix<Dune::FieldMatrix<double,1,1>>;
using ElementRhs = Dune::BlockVector<Dune::FieldVector<double,1>>;
rhs.resize(size);
rhs = 0;
Dune::ImplicitMatrixBuilder<Matrix> matrixBuilder(matrix, size, size, 10, 0.1);
Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache;
ElementMatrix elementMatrix;
ElementRhs elementRhs;
for(const auto& e: Dune::elements(gridView))
{
const auto& localFiniteElement = feCache.get(e.type());
std::size_t localSize = localFiniteElement.localBasis().size();
assembleLocalStiffnessMatrix(e, localFiniteElement, elementMatrix);
assembleLocalRHS(e, localFiniteElement, elementRhs, rhsFunction);
for(std::size_t i=0; i<localSize; ++i)
{
std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim);
for(std::size_t j=0; j<localSize; ++j)
{
std::size_t globalJ = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(j).subEntity(), dim);
matrixBuilder[globalI][globalJ] += elementMatrix[i][j];
}
rhs[globalI] = elementRhs[i];
}
}
matrix.compress();
}
template<class GridView, class BitVector>
void computeBoundaryVertices(const GridView& gridView, BitVector& isBoundary)
{
using namespace Dune;
using namespace Dune::FuTutorial;
const auto& indexSet = gridView.indexSet();
for(const auto& element: elements(gridView))
for(const auto& intersection: intersections(gridView, element))
if (intersection.boundary())
for(const auto& vertex: subVertices(referenceElement(element), insideFacet(intersection)))
isBoundary[subIndex(indexSet, element, vertex)] = true;
}
#endif // DUNE_FUTUTORIAL_SRC_05POISSON_PROBLEM_HH
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment