Loading...
Searching...
No Matches
Poisson equation (C++)

The following explains how to solve the Poisson equation using dune-functions. The full example poisson-pq2.cc can be found in the examples/ subdirectory.

### Local assemblers

The program first includes the required headers and then defines free functions for assembling the Poisson problem. The function getLocalMatrix() implements the assembler of the local stiffness matrix for the bilinear form .

// Compute the stiffness matrix for a single element
template <class LocalView, class MatrixType>
void getLocalMatrix(const LocalView& localView, MatrixType& elementMatrix)
{
// Get the grid element from the local FE basis view
typedef typename LocalView::Element Element;
const Element& element = localView.element();
const int dim = Element::dimension;
auto geometry = element.geometry();
// Get set of shape functions for this element
const auto& localFiniteElement = localView.tree().finiteElement();
// Set all matrix entries to zero
elementMatrix.setSize(localFiniteElement.localBasis().size(),localFiniteElement.localBasis().size());
elementMatrix = 0; // fills the entire matrix with zeroes
// Get a quadrature rule
int order = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
// Loop over all quadrature points
for (size_t pt=0; pt < quad.size(); pt++) {
// Position of the current quadrature point in the reference element
const FieldVector<double,dim>& quadPos = quad[pt].position();
// The inverse Jacobian of the map from the reference element to the element
const auto& jacobianInverse = geometry.jacobianInverse(quadPos);
// The multiplicative factor in the integral transformation formula
const double integrationElement = geometry.integrationElement(quadPos);
// The gradients of the shape functions on the reference element
std::vector<FieldMatrix<double,1,dim> > referenceJacobians;
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceJacobians);
// Compute the shape function gradients on the real element
std::vector<FieldMatrix<double,1,dim> > jacobians(referenceJacobians.size());
for (size_t i=0; i<jacobians.size(); i++)
jacobians[i] = referenceJacobians[i] * jacobianInverse;
// 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] += (jacobians[i] * transpose(jacobians[j])) * quad[pt].weight() * integrationElement;
}
}

The getVolumeTerm() functions implements the local assembler for the volume right hand side term .

// Compute the source term for a single element
template <class LocalView, class LocalVolumeTerm>
void getVolumeTerm( const LocalView& localView,
BlockVector<double>& localRhs,
LocalVolumeTerm&& localVolumeTerm)
{
// Get the grid element from the local FE basis view
typedef typename LocalView::Element Element;
const Element& element = localView.element();
const int dim = Element::dimension;
// Get set of shape functions for this element
const auto& localFiniteElement = localView.tree().finiteElement();
// Set all entries to zero
localRhs.resize(localFiniteElement.localBasis().size());
localRhs = 0;
// A quadrature rule
int order = dim*localFiniteElement.localBasis().order();
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), order);
// Loop over all quadrature points
for ( size_t pt=0; pt < quad.size(); pt++ ) {
// Position of the current quadrature point in the reference element
const FieldVector<double,dim>& quadPos = quad[pt].position();
// The multiplicative factor in the integral transformation formula
const double integrationElement = element.geometry().integrationElement(quadPos);
double functionValue = localVolumeTerm(quadPos);
// Evaluate all shape function values at this point
std::vector<FieldVector<double,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
// Actually compute the vector entries
for (size_t i=0; i<localRhs.size(); i++)
localRhs[i] += shapeFunctionValues[i] * functionValue * quad[pt].weight() * integrationElement;
}
}

### Global assembler

Assemble the global matrix pattern.

// Get the occupation pattern of the stiffness matrix
template <class FEBasis>
void getOccupationPattern(const FEBasis& feBasis, MatrixIndexSet& nb)
{
// Total number of grid vertices
auto n = feBasis.size();
nb.resize(n, n);
// A view on the FE basis on a single element
auto localView = feBasis.localView();
// Loop over all leaf elements
for(const auto& e : elements(feBasis.gridView()))
{
// Bind the local FE basis view to the current element
localView.bind(e);
// There is a matrix entry a_ij if the i-th and j-th vertex are connected in the grid
for (size_t i=0; i<localView.tree().size(); i++) {
for (size_t j=0; j<localView.tree().size(); j++) {
auto iIdx = localView.index(i);
auto jIdx = localView.index(j);
// Add a nonzero entry to the matrix
nb.add(iIdx, jIdx);
}
}
}
}

Assembly of matrix and right-hand-side.

template <class FEBasis, class VolumeTerm>
void assembleLaplaceMatrix(const FEBasis& feBasis,
BCRSMatrix<double>& matrix,
BlockVector<double>& rhs,
VolumeTerm&& volumeTerm)
{
// Get the grid view from the finite element basis
typedef typename FEBasis::GridView GridView;
GridView gridView = feBasis.gridView();
auto localVolumeTerm = localFunction(Functions::makeGridViewFunction(volumeTerm, gridView));
// MatrixIndexSets store the occupation pattern of a sparse matrix.
// They are not particularly efficient, but simple to use.
MatrixIndexSet occupationPattern;
getOccupationPattern(feBasis, occupationPattern);
// ... and give it the occupation pattern we want.
occupationPattern.exportIdx(matrix);
// set rhs to correct length -- the total number of basis vectors in the basis
rhs.resize(feBasis.size());
// Set all entries to zero
matrix = 0;
rhs = 0;
// A view on the FE basis on a single element
auto localView = feBasis.localView();
// A loop over all elements of the grid
for(const auto& e : elements(gridView))
{
// Bind the local FE basis view to the current element
localView.bind(e);
// Now let's get the element stiffness matrix
// A dense matrix is used for the element stiffness matrix
Matrix<double> elementMatrix;
getLocalMatrix(localView, elementMatrix);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t i=0; i<elementMatrix.N(); i++) {
// The global index of the i-th local degree of freedom of the element 'e'
auto row = localView.index(i);
for (size_t j=0; j<elementMatrix.M(); j++ ) {
// The global index of the j-th local degree of freedom of the element 'e'
auto col = localView.index(j);
matrix[row][col] += elementMatrix[i][j];
}
}
// Now get the local contribution to the right-hand side vector
BlockVector<double> localRhs;
localVolumeTerm.bind(e);
getVolumeTerm(localView, localRhs, localVolumeTerm);
for (size_t i=0; i<localRhs.size(); i++) {
// The global index of the i-th vertex of the element 'e'
auto row = localView.index(i);
rhs[row] += localRhs[i];
}
}
}

### Helper functions

Treatment of boundary condition.

template <class FEBasis>
void boundaryTreatment (const FEBasis& feBasis, std::vector<char>& dirichletNodes )
{
dirichletNodes.clear();
dirichletNodes.resize(feBasis.size(), false);
Functions::forEachBoundaryDOF(feBasis, [&] (auto&& index) {
dirichletNodes[index] = true;
});
}

Create a mixed grid containing triangles and quadrilaterals.

auto createMixedGrid()
{
using Grid = Dune::UGGrid<2>;
Dune::GridFactory<Grid> factory;
for(unsigned int k : Dune::range(9))
factory.insertVertex({0.5*(k%3), 0.5*(k/3)});
factory.insertElement(Dune::GeometryTypes::cube(2), {0, 1, 3, 4});
factory.insertElement(Dune::GeometryTypes::cube(2), {1, 2, 4, 5});
factory.insertElement(Dune::GeometryTypes::simplex(2), {3, 4, 6});
factory.insertElement(Dune::GeometryTypes::simplex(2), {4, 7, 6});
factory.insertElement(Dune::GeometryTypes::simplex(2), {4, 5, 7});
factory.insertElement(Dune::GeometryTypes::simplex(2), {5, 8, 7});
return factory.createGrid();
}

### Setup and initialization

Initialize MPI

int main (int argc, char *argv[]) try
{
// Set up MPI, if available
MPIHelper::instance(argc, argv);

### Create grid and finite element basis

Create a mixed grid and obtain a leaf grid view.

Note
  A GridView is a view to a subset of the grid's elements, vertices, ... that should be stored by value and can be copied cheaply. Grids in dune are in general hierarchical and composed by elements on several levels. The discretization usually lives on the set of most refined elements that is denote the leaf GridView in dune.

auto grid = createMixedGrid();
grid->globalRefine(4);
auto gridView = grid->leafGridView();
using GridView = decltype(gridView);

As a next step the program creates a global finite element basis on the GridView.

typedef Functions::LagrangeBasis<GridView,2> FEBasis;
FEBasis feBasis(gridView);

And here is the rest of the file:

using VectorType = BlockVector<double>;
using MatrixType = BCRSMatrix<double>;
VectorType rhs;
MatrixType stiffnessMatrix;
// Assemble the system
std::cout << "Number of DOFs is " << feBasis.dimension() << std::endl;
auto rightHandSide = [] (const auto& x) { return 10;};
Dune::Timer timer;
assembleLaplaceMatrix(feBasis, stiffnessMatrix, rhs, rightHandSide);
std::cout << "Assembling the problem took " << timer.elapsed() << "s" << std::endl;
// Choose an initial iterate
VectorType x(feBasis.size());
x = 0;
// Determine Dirichlet dofs
std::vector<char> dirichletNodes;
boundaryTreatment(feBasis, dirichletNodes);
// Don't trust on non-standard M_PI.
auto pi = std::acos(-1.0);
auto dirichletValueFunction = [pi](const auto& x){ return std::sin(2*pi*x[0]); };
// Interpolate dirichlet values at the boundary nodes
interpolate(feBasis, x, dirichletValueFunction, dirichletNodes);
// Incorporate Dirichlet values in a symmetric way
// Compute residual for non-homogeneous Dirichlet values
// stored in x. Since x is zero for non-Dirichlet DOFs,
// we can simply multiply by the matrix.
stiffnessMatrix.mmv(x, rhs);
// Change Dirichlet matrix rows and columns to the identity.
for (size_t i=0; i<stiffnessMatrix.N(); i++) {
if (dirichletNodes[i]) {
rhs[i] = x[i];
// loop over nonzero matrix entries in current row
for (auto&& [entry, idx] : sparseRange(stiffnessMatrix[i]))
entry = (i==idx) ? 1.0 : 0.0;
}
else {
// loop over nonzero matrix entries in current row
for (auto&& [entry, idx] : sparseRange(stiffnessMatrix[i]))
if (dirichletNodes[idx])
entry = 0.0;
}
}
// Compute solution
// Technicality: turn the matrix into a linear operator
MatrixAdapter<MatrixType,VectorType,VectorType> op(stiffnessMatrix);
// Sequential incomplete LU decomposition as the preconditioner
SeqILDL<MatrixType,VectorType,VectorType> ildl(stiffnessMatrix,1.0);
// Preconditioned conjugate-gradient solver
CGSolver<VectorType> cg(op,
ildl, // preconditioner
1e-4, // desired residual reduction factor
100, // maximum number of iterations
2); // verbosity of the solver
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
// Solve!
cg.apply(x, rhs, statistics);
// Make a discrete function from the FE basis and the coefficient vector
auto xFunction = Functions::makeDiscreteGlobalBasisFunction<double>(feBasis, x);
// Write result to VTK file
// We need to subsample, because VTK cannot natively display real second-order functions
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));
vtkWriter.addVertexData(xFunction, VTK::FieldInfo("x", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("poisson-pq2");
}
// Error handling
catch (Exception& e) {
std::cout << e.what() << std::endl;
}

Previous | Next
--- | ---
Examples | Poisson equation (Python)