Commit 5b2a3ca8 authored by Sven Marnach's avatar Sven Marnach

Added P1 finite element example (code by Dominic Kempf and Sven Marnach)

[[Imported from SVN: r271]]
parent 298fab97
......@@ -17,6 +17,7 @@ stamp-h1
adaptivefinitevolume
adaptiveintegration
finitevolume
finiteelements
gettingstarted
groundwater
integration
......
......@@ -12,15 +12,16 @@ examples_HEADERS = elementdata.hh parfvdatahandle.hh unitcube_sgrid.hh \
finitevolumeadapt.hh transportproblem.hh unitcube_yaspgrid.hh \
functors.hh unitcube_albertagrid.hh vertexdata.hh \
initialize.hh unitcube_alugrid.hh vtkout.hh integrateentity.hh \
unitcube.hh parevolve.hh unitcube_onedgrid.hh basicunitcube.hh
unitcube.hh parevolve.hh unitcube_onedgrid.hh basicunitcube.hh \
shapefunctions.hh
examples_PROGRAMS = gettingstarted traversal integration othergrids\
adaptiveintegration finitevolume adaptivefinitevolume parfinitevolume\
visualization
visualization finiteelements
examples_DATA = gettingstarted.cc traversal.cc integration.cc othergrids.cc \
adaptiveintegration.cc finitevolume.cc adaptivefinitevolume.cc \
parfinitevolume.cc visualization.cc
parfinitevolume.cc visualization.cc finiteelements.cc
gettingstarted_SOURCES = gettingstarted.cc
gettingstarted_CXXFLAGS = $(MPI_CPPFLAGS) $(ALL_PKG_CPPFLAGS)
......@@ -58,6 +59,10 @@ visualization_SOURCES = visualization.cc
visualization_CXXFLAGS = $(MPI_CPPFLAGS) $(ALL_PKG_CPPFLAGS)
visualization_LDADD = $(DUNE_LDFLAGS) $(DUNE_LIBS) $(ALL_PKG_LDFLAGS) $(ALL_PKG_LIBS) $(MPI_LDFLAGS)
finiteelements_SOURCES = finiteelements.cc
finiteelements_CXXFLAGS = $(MPI_CPPFLAGS) $(ALL_PKG_CPPFLAGS)
finiteelements_LDADD = $(DUNE_LDFLAGS) $(DUNE_LIBS) $(ALL_PKG_LDFLAGS) $(ALL_PKG_LIBS) $(MPI_LDFLAGS)
# don't follow the full GNU-standard
# we need automake 1.5
AUTOMAKE_OPTIONS = foreign 1.5
......
......@@ -2,4 +2,4 @@
Module: dune-grid-howto
Maintainer: dune@dune-project.org
Version: 1.3svn
Depends: dune-common (>= 1.3) dune-grid (>= 1.3)
Depends: dune-common (>= 1.3) dune-grid (>= 1.3) dune-istl (>= 1.3)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
#include <iostream>
#include <vector>
#include <set>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/grid/common/quadraturerules.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/ilu.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/io.hh>
#include <dune/grid/albertagrid.hh>
#include "shapefunctions.hh"
// P1Elements:
// a P1 finite element discretization for elliptic problems Dirichlet
// boundary conditions on simplicial conforming grids
template<class GV, class F>
class P1Elements
{
public:
static const int dim = GV::dimension;
typedef typename GV::ctype ctype;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<ctype,1,1> > Matrix;
typedef Dune::BlockVector<Dune::FieldVector<ctype,1> > ScalarField;
private:
typedef typename GV::template Codim<0>::Iterator LeafIterator;
typedef typename GV::IntersectionIterator IntersectionIterator;
typedef typename GV::IndexSet LeafIndexSet;
const GV& gv;
const F& f;
public:
Matrix A;
ScalarField b;
ScalarField u;
std::vector< std::set<int> > adjacencyPattern;
P1Elements(const GV& gv_, const F& f_) : gv(gv_), f(f_) {}
// store adjacency information in a vector of sets
void determineAdjacencyPattern();
// assemble stiffness matrix A and right side b
void assemble();
// finally solve Au = b for u
void solve();
};
template<class GV, class F>
void P1Elements<GV, F>::determineAdjacencyPattern()
{
const int N = gv.size(dim);
adjacencyPattern.resize(N);
const LeafIndexSet& set = gv.indexSet();
const LeafIterator itend = gv.template end<0>();
for (LeafIterator it = gv.template begin<0>(); it != itend; ++it)
{
Dune::GeometryType gt = it->type();
const Dune::template GenericReferenceElement<ctype,dim> &ref =
Dune::GenericReferenceElements<ctype,dim>::general(gt);
// traverse all codim-1-entities of the current element and store all
// pairs of vertices in adjacencyPattern
const IntersectionIterator isend = gv.iend(*it);
for (IntersectionIterator is = gv.ibegin(*it) ; is != isend ; ++is)
{
int vertexsize = ref.size(is->indexInInside(),1,dim);
for (int i=0; i < vertexsize; i++)
{
int indexi = set.subIndex(*it,ref.subEntity(is->indexInInside(),1,i,dim),dim);
for (int j=0; j < vertexsize; j++)
{
int indexj = set.subIndex(*it,ref.subEntity(is->indexInInside(),1,j,dim),dim);
adjacencyPattern[indexi].insert(indexj);
}
}
}
}
}
template<class GV, class F>
void P1Elements<GV, F>::assemble()
{
const int N = gv.size(dim);
const LeafIndexSet& set = gv.indexSet();
const LeafIterator itend = gv.template end<0>();
// set sizes of A and b
A.setSize(N, N, N + 2*gv.size(dim-1));
A.setBuildMode(Matrix::random);
b.resize(N, false);
for (int i = 0; i < N; i++)
A.setrowsize(i,adjacencyPattern[i].size());
A.endrowsizes();
// set sparsity pattern of A with the information gained in determineAdjacencyPattern
for (int i = 0; i < N; i++)
{
std::template set<int>::iterator setend = adjacencyPattern[i].end();
for (std::template set<int>::iterator setit = adjacencyPattern[i].begin();
setit != setend; ++setit)
A.addindex(i,*setit);
}
A.endindices();
// initialize A and b
A = 0.0;
b = 0.0;
// get a set of P1 shape functions
P1ShapeFunctionSet<ctype,ctype,dim> basis = P1ShapeFunctionSet<ctype,ctype,dim>::instance();
for (LeafIterator it = gv.template begin<0>(); it != itend; ++it)
{
// determine geometry type of the current element and get the matching reference element
Dune::GeometryType gt = it->type();
const Dune::template GenericReferenceElement<ctype,dim> &ref =
Dune::GenericReferenceElements<ctype,dim>::general(gt);
int vertexsize = ref.size(dim);
// get a quadrature rule of order one for the given geometry type
const Dune::QuadratureRule<ctype,dim>& rule = Dune::QuadratureRules<ctype,dim>::rule(gt,1);
for (typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule.begin();
r != rule.end() ; ++r)
{
// compute the jacobian inverse transposed to transform the gradients
Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
it->geometry().jacobianInverseTransposed(r->position());
// get the weight at the current quadrature point
ctype weight = r->weight();
// compute Jacobian determinant for the transformation formula
ctype detjac = it->geometry().integrationElement(r->position());
for (int i = 0; i < vertexsize; i++)
{
// compute transformed gradients
Dune::FieldVector<ctype,dim> grad1;
jacInvTra.mv(basis[i].evaluateGradient(r->position()),grad1);
for (int j = 0; j < vertexsize; j++)
{
Dune::FieldVector<ctype,dim> grad2;
jacInvTra.mv(basis[j].evaluateGradient(r->position()),grad2);
// gain global inidices of vertices i and j and update associated matrix entry
A[set.subIndex(*it,i,dim)][set.subIndex(*it,j,dim)]
+= (grad1*grad2) * weight * detjac;
}
}
}
// get a quadrature rule of order two for the given geometry type
const Dune::QuadratureRule<ctype,dim>& rule2 = Dune::QuadratureRules<ctype,dim>::rule(gt,2);
for (typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule2.begin();
r != rule2.end() ; ++r)
{
ctype weight = r->weight();
ctype detjac = it->geometry().integrationElement(r->position());
for (int i = 0 ; i<vertexsize; i++)
{
// evaluate the integrand of the right side
ctype fval = basis[i].evaluateFunction(it->geometry().global(r->position()))
* f(it->geometry().global(r->position())) ;
b[set.subIndex(*it,i,dim)] += fval * weight * detjac;
}
}
}
// Dirichlet boundary conditions:
// replace lines in A related to Dirichlet vertices by trivial lines
for ( LeafIterator it = gv.template begin<0>() ; it != itend ; ++it)
{
const IntersectionIterator isend = gv.iend(*it);
for (IntersectionIterator is = gv.ibegin(*it) ; is != isend ; ++is)
{
// determine geometry type of the current element and get the matching reference element
Dune::GeometryType gt = it->type();
const Dune::template GenericReferenceElement<ctype,dim> &ref =
Dune::GenericReferenceElements<ctype,dim>::general(gt);
// check whether current intersection is on the boundary
if ( is->boundary() )
{
// traverse all vertices the intersection consists of
for (int i=0; i < ref.size(is->indexInInside(),1,dim); i++)
{
// and replace the associated line of A and b with a trivial one
int indexi = set.subIndex(*it,ref.subEntity(is->indexInInside(),1,i,dim),dim);
A[indexi] = 0.0;
A[indexi][indexi] = 1.0;
b[indexi] = 0.0;
}
}
}
}
}
template<class GV, class E>
void P1Elements<GV, E>::solve()
{
// make linear operator from A
Dune::MatrixAdapter<Matrix,ScalarField,ScalarField> op(A);
// initialize preconditioner
Dune::SeqILUn<Matrix,ScalarField,ScalarField> ilu1(A, 1, 0.92);
// the inverse operator
Dune::BiCGSTABSolver<ScalarField> bcgs(op, ilu1, 1e-15, 5000, 0);
Dune::InverseOperatorResult r;
// initialue u to some arbitrary value to avoid u being the exact
// solution
u.resize(b.N(), false);
u = 2.0;
// finally solve the system
bcgs.apply(u, b, r);
}
// an example right hand side function
template<class ctype, int dim>
class Bump {
public:
ctype operator() (Dune::FieldVector<ctype,dim> x) const
{
return 2.0 * (x[0]*(1-x[0]) + x[1]*(1-x[1]));
}
};
int main(int argc, char** argv)
{
static const int dim = 2;
const char* gridfile = "grids/2dgrid.al";
typedef Dune::AlbertaGrid<dim,dim> GridType;
typedef GridType::LeafGridView GV;
typedef GridType::ctype ctype;
typedef Bump<ctype,dim> Func;
GridType grid(gridfile);
const GV& gv=grid.leafView();
Func f;
P1Elements<GV,Func> p1(gv, f);
grid.globalRefine(1);
std::cout << "-----------------------------------" << "\n";
std::cout << "number of unknowns: " << grid.size(dim) << "\n";
std::cout << "determine adjacency pattern..." << "\n";
p1.determineAdjacencyPattern();
std::cout << "assembling..." << "\n";
p1.assemble();
std::cout << "solving..." << "\n";
p1.solve();
std::cout << "visualizing..." << "\n";
Dune::VTKWriter<GridType::LeafGridView> vtkwriter(grid.leafView());
vtkwriter.addVertexData(p1.u, "u");
vtkwriter.write("test", Dune::VTKOptions::binaryappended);
}
DIM: 3
DIM_OF_WORLD: 3
number of vertices: 8
number of elements: 6
vertex coordinates:
0.0 0.0 0.0
1.0 0.0 0.0
0.0 0.0 1.0
1.0 0.0 1.0
1.0 1.0 0.0
1.0 1.0 1.0
0.0 1.0 0.0
0.0 1.0 1.0
element vertices:
0 5 4 1
0 5 3 1
0 5 3 2
0 5 4 6
0 5 7 6
0 5 7 2
element boundaries:
1 1 0 0
1 1 0 0
1 1 0 0
1 1 0 0
1 1 0 0
1 1 0 0
element neighbours:
-1 -1 1 3
-1 -1 0 2
-1 -1 5 1
-1 -1 4 0
-1 -1 3 5
-1 -1 2 4
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef SHAPEFUNCTIONS_HH
#define SHAPEFUNCTIONS_HH
#include <dune/common/fvector.hh>
// LinearShapeFunction:
// represents a shape function and provides methods to evaluate the function
// and its gradient
template<class ctype, class rtype, int dim>
class LinearShapeFunction
{
public:
enum { dimension = dim };
LinearShapeFunction() : coeff0(0.0), coeff1(0.0) {}
LinearShapeFunction(rtype coeff0_, const Dune::FieldVector<rtype,dim>& coeff1_)
: coeff0(coeff0_), coeff1(coeff1_) {}
void setCoeff(rtype coeff0_, const Dune::FieldVector<rtype,dim>& coeff1_)
{
coeff0 = coeff0_;
coeff1 = coeff1_;
}
rtype evaluateFunction(const Dune::FieldVector<ctype,dim>& local) const
{
ctype result = coeff0;
for (int i = 0; i < dim; ++i)
result += coeff1[i] * local[i];
return result;
}
Dune::FieldVector<rtype,dim>
evaluateGradient(const Dune::FieldVector<ctype,dim>& local) const
{
return coeff1;
}
private:
rtype coeff0;
Dune::FieldVector<rtype,dim> coeff1;
};
// P1ShapeFunctionSet
// initializes one and only one set of LinearShapeFunction
template<class ctype, class rtype, int dim>
class P1ShapeFunctionSet
{
public:
enum { n = dim + 1 };
typedef LinearShapeFunction<ctype,rtype,dim> ShapeFunction;
typedef rtype resulttype;
// get the only instance of this class
static const P1ShapeFunctionSet& instance()
{
static const P1ShapeFunctionSet sfs;
return sfs;
}
const ShapeFunction& operator[](int i) const
{
if (!i)
return f0;
else
return f1[i - 1];
}
private:
// private constructor prevents additional instances
P1ShapeFunctionSet()
{
Dune::FieldVector<rtype,dim> e(-1.0);
f0.setCoeff(1.0, e);
for (int i = 0; i < dim; ++i)
{
Dune::FieldVector<rtype,dim> e(0.0);
e[i] = 1.0;
f1[i].setCoeff(0.0, e);
}
}
ShapeFunction f0;
ShapeFunction f1[dim];
};
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment