"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "finiteelements.cc" between
dune-grid-howto-2.6.0.tar.gz and dune-grid-howto-2.7.0.tar.gz

About: dune-grid-howto - DUNE (Distributed and Unified Numerics Environment) is a modular toolbox for solving partial differential equations (PDEs) with grid-based methods: DUNE tutorial on the grid interface.

finiteelements.cc  (dune-grid-howto-2.6.0):finiteelements.cc  (dune-grid-howto-2.7.0)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2: // vi: set et ts=4 sw=2 sts=2:
#include <config.h> #include <config.h>
#include <iostream> #include <iostream>
#include <vector> #include <vector>
#include <set> #include <set>
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh> #include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/albertagrid.hh> #include <dune/grid/albertagrid.hh>
#if HAVE_DUNE_ISTL #if HAVE_DUNE_ISTL
#include <dune/istl/bvector.hh> #include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh> #include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/ilu.hh> #include <dune/istl/ilu.hh>
#include <dune/istl/operators.hh> #include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh> #include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh> #include <dune/istl/preconditioners.hh>
skipping to change at line 84 skipping to change at line 85
void P1Elements<GV, F>::determineAdjacencyPattern() void P1Elements<GV, F>::determineAdjacencyPattern()
{ {
const int N = gv.size(dim); const int N = gv.size(dim);
adjacencyPattern.resize(N); adjacencyPattern.resize(N);
const LeafIndexSet& set = gv.indexSet(); const LeafIndexSet& set = gv.indexSet();
const LeafIterator itend = gv.template end<0>(); const LeafIterator itend = gv.template end<0>();
for (LeafIterator it = gv.template begin<0>(); it != itend; ++it) for (LeafIterator it = gv.template begin<0>(); it != itend; ++it)
{ {
Dune::GeometryType gt = it->type(); auto geo = it->geometry();
const Dune::template ReferenceElement<ctype,dim> &ref = auto ref = referenceElement(geo);
Dune::ReferenceElements<ctype,dim>::general(gt);
// traverse all codim-1-entities of the current element and store all // traverse all codim-1-entities of the current element and store all
// pairs of vertices in adjacencyPattern // pairs of vertices in adjacencyPattern
const IntersectionIterator isend = gv.iend(*it); const IntersectionIterator isend = gv.iend(*it);
for (IntersectionIterator is = gv.ibegin(*it) ; is != isend ; ++is) for (IntersectionIterator is = gv.ibegin(*it) ; is != isend ; ++is)
{ {
int vertexsize = ref.size(is->indexInInside(),1,dim); int vertexsize = ref.size(is->indexInInside(),1,dim);
for (int i=0; i < vertexsize; i++) for (int i=0; i < vertexsize; i++)
{ {
int indexi = set.subIndex(*it,ref.subEntity(is->indexInInside(),1,i,dim) ,dim); int indexi = set.subIndex(*it,ref.subEntity(is->indexInInside(),1,i,dim) ,dim);
skipping to change at line 119 skipping to change at line 119
{ {
const int N = gv.size(dim); const int N = gv.size(dim);
const LeafIndexSet& set = gv.indexSet(); const LeafIndexSet& set = gv.indexSet();
const LeafIterator itend = gv.template end<0>(); const LeafIterator itend = gv.template end<0>();
// set sizes of A and b // set sizes of A and b
#if HAVE_DUNE_ISTL #if HAVE_DUNE_ISTL
A.setSize(N, N, N + 2*gv.size(dim-1)); A.setSize(N, N, N + 2*gv.size(dim-1));
A.setBuildMode(Matrix::random); A.setBuildMode(Matrix::random);
b.resize(N, false); b.resize(N);
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++)
A.setrowsize(i,adjacencyPattern[i].size()); A.setrowsize(i,adjacencyPattern[i].size());
A.endrowsizes(); A.endrowsizes();
// set sparsity pattern of A with the information gained in determineAdjacency Pattern // set sparsity pattern of A with the information gained in determineAdjacency Pattern
for (int i = 0; i < N; i++) /*@\label{fem:setpattern}@*/ for (int i = 0; i < N; i++) /*@\label{fem:setpattern}@*/
{ {
std::template set<int>::iterator setend = adjacencyPattern[i].end(); std::template set<int>::iterator setend = adjacencyPattern[i].end();
for (std::template set<int>::iterator setit = adjacencyPattern[i].begin(); for (std::template set<int>::iterator setit = adjacencyPattern[i].begin();
skipping to change at line 149 skipping to change at line 149
// initialize A and b // initialize A and b
A = 0.0; A = 0.0;
b = 0.0; b = 0.0;
// get a set of P1 shape functions // get a set of P1 shape functions
const P1ShapeFunctionSet<ctype,ctype,dim>& basis = P1ShapeFunctionSet<ctype,ct ype,dim>::instance(); const P1ShapeFunctionSet<ctype,ctype,dim>& basis = P1ShapeFunctionSet<ctype,ct ype,dim>::instance();
for (LeafIterator it = gv.template begin<0>(); it != itend; ++it) /*@\label{ fem:loop1}@*/ for (LeafIterator it = gv.template begin<0>(); it != itend; ++it) /*@\label{ fem:loop1}@*/
{ {
// determine geometry type of the current element and get the matching refer // determine geometry of the current element and get the matching reference
ence element element
Dune::GeometryType gt = it->type(); auto geo = it->geometry();
const Dune::template ReferenceElement<ctype,dim> &ref = auto ref = referenceElement(geo);
Dune::ReferenceElements<ctype,dim>::general(gt);
int vertexsize = ref.size(dim); int vertexsize = ref.size(dim);
// get a quadrature rule of order one for the given geometry type // get a quadrature rule of order one for the given geometry type
Dune::GeometryType gt = it->type();
const Dune::QuadratureRule<ctype,dim>& rule = Dune::QuadratureRules<ctype,di m>::rule(gt,1); const Dune::QuadratureRule<ctype,dim>& rule = Dune::QuadratureRules<ctype,di m>::rule(gt,1);
for (typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule.begin (); for (typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule.begin ();
r != rule.end() ; ++r) r != rule.end() ; ++r)
{ {
// compute the jacobian inverse transposed to transform the gradients // compute the jacobian inverse transposed to transform the gradients
JacobianInverseTransposed jacInvTra = JacobianInverseTransposed jacInvTra =
it->geometry().jacobianInverseTransposed(r->position()); it->geometry().jacobianInverseTransposed(r->position());
// get the weight at the current quadrature point // get the weight at the current quadrature point
ctype weight = r->weight(); ctype weight = r->weight();
skipping to change at line 210 skipping to change at line 210
} }
} /*@\label{fem:loop2}@*/ } /*@\label{fem:loop2}@*/
// Dirichlet boundary conditions: // Dirichlet boundary conditions:
// replace lines in A related to Dirichlet vertices by trivial lines // replace lines in A related to Dirichlet vertices by trivial lines
for ( LeafIterator it = gv.template begin<0>() ; it != itend ; ++it) /*@\lab el{fem:boundary1}@*/ for ( LeafIterator it = gv.template begin<0>() ; it != itend ; ++it) /*@\lab el{fem:boundary1}@*/
{ {
const IntersectionIterator isend = gv.iend(*it); const IntersectionIterator isend = gv.iend(*it);
for (IntersectionIterator is = gv.ibegin(*it) ; is != isend ; ++is) for (IntersectionIterator is = gv.ibegin(*it) ; is != isend ; ++is)
{ {
// determine geometry type of the current element and get the matching ref // determine geometry of the current element and get the matching referenc
erence element e element
Dune::GeometryType gt = it->type(); auto geo = it->geometry();
const Dune::template ReferenceElement<ctype,dim> &ref = auto ref = referenceElement(geo);
Dune::ReferenceElements<ctype,dim>::general(gt);
// check whether current intersection is on the boundary // check whether current intersection is on the boundary
if ( is->boundary() ) if ( is->boundary() )
{ {
// traverse all vertices the intersection consists of // traverse all vertices the intersection consists of
for (int i=0; i < ref.size(is->indexInInside(),1,dim); i++) 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 // and replace the associated line of A and b with a trivial one
int indexi = set.subIndex(*it,ref.subEntity(is->indexInInside(),1,i,di m),dim); int indexi = set.subIndex(*it,ref.subEntity(is->indexInInside(),1,i,di m),dim);
skipping to change at line 241 skipping to change at line 240
} }
#if HAVE_DUNE_ISTL #if HAVE_DUNE_ISTL
template<class GV, class E> template<class GV, class E>
void P1Elements<GV, E>::solve() void P1Elements<GV, E>::solve()
{ {
// make linear operator from A // make linear operator from A
Dune::MatrixAdapter<Matrix,ScalarField,ScalarField> op(A); Dune::MatrixAdapter<Matrix,ScalarField,ScalarField> op(A);
// initialize preconditioner // initialize preconditioner
Dune::SeqILUn<Matrix,ScalarField,ScalarField> ilu1(A, 1, 0.92); Dune::SeqILU<Matrix,ScalarField,ScalarField> ilu1(A, 1, 0.92);
// the inverse operator // the inverse operator
Dune::BiCGSTABSolver<ScalarField> bcgs(op, ilu1, 1e-15, 5000, 0); Dune::BiCGSTABSolver<ScalarField> bcgs(op, ilu1, 1e-15, 5000, 0);
Dune::InverseOperatorResult r; Dune::InverseOperatorResult r;
// initialize u to some arbitrary value to avoid u being the exact // initialize u to some arbitrary value to avoid u being the exact
// solution // solution
u.resize(b.N(), false); u.resize(b.N());
u = 2.0; u = 2.0;
// finally solve the system // finally solve the system
bcgs.apply(u, b, r); bcgs.apply(u, b, r);
} }
#endif // HAVE_DUNE_ISTL #endif // HAVE_DUNE_ISTL
// an example right hand side function // an example right hand side function
template<class ctype, int dim> template<class ctype, int dim>
class Bump { class Bump {
skipping to change at line 274 skipping to change at line 273
for (int i=0 ; i < dim ; i++) for (int i=0 ; i < dim ; i++)
result += 2.0 * x[i]* (1-x[i]); result += 2.0 * x[i]* (1-x[i]);
return result; return result;
} }
}; };
int main(int argc, char** argv) int main(int argc, char** argv)
{ {
#if HAVE_ALBERTA && ALBERTA_DIM==2 #if HAVE_ALBERTA && ALBERTA_DIM==2
static const int dim = 2; /*@\label{fem:dim}@*/ static const int dim = 2; /*@\label{fem:dim}@*/
const char* gridfile = "grids/2dgrid.al"; /*@\label{fem:file}@*/ std::stringstream gridfile;
gridfile << DUNE_GRID_HOWTO_EXAMPLE_GRIDS_PATH
<< "2dgrid.al"; /*@\label{fem:file}@*/
typedef Dune::AlbertaGrid<dim,dim> GridType; typedef Dune::AlbertaGrid<dim,dim> GridType;
typedef GridType::LeafGridView GV; typedef GridType::LeafGridView GV;
typedef GridType::ctype ctype; typedef GridType::ctype ctype;
typedef Bump<ctype,dim> Func; typedef Bump<ctype,dim> Func;
GridType grid(gridfile); GridType grid(gridfile.str());
const GV& gv = grid.leafGridView(); const GV& gv = grid.leafGridView();
Func f; Func f;
P1Elements<GV,Func> p1(gv, f); P1Elements<GV,Func> p1(gv, f);
#if HAVE_DUNE_ISTL #if HAVE_DUNE_ISTL
grid.globalRefine(16); grid.globalRefine(16);
#else #else
grid.globalRefine(10); grid.globalRefine(10);
#endif // HAVE_DUNE_ISTL #endif // HAVE_DUNE_ISTL
 End of changes. 10 change blocks. 
18 lines changed or deleted 19 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)