Skip to content
Snippets Groups Projects
Commit fbe7cc69 authored by Peter Bastian's avatar Peter Bastian
Browse files

- use capabilities to avoid hanging node computations if there are not any 

- base class for local assembler

[[Imported from SVN: r4142]]
parent d70a5d09
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
disc/groundwater/p1groundwater.hh
+ 41
− 130
View file @ fbe7cc69
......@@ -66,7 +66,7 @@ namespace Dune
- RT type used for return values
*/
template<class G, class RT>
class GroundwaterEquationLocalStiffness
class GroundwaterEquationLocalStiffness : public P1LocalStiffness<typename G::ctype,RT,G::dimension,1>
{
// grid types
typedef typename G::ctype DT;
......@@ -74,25 +74,18 @@ namespace Dune
typedef typename G::template Codim<0>::EntityPointer EEntityPointer;
typedef typename G::Traits::IntersectionIterator IntersectionIterator;
// some other sizes
enum {n=G::dimension};
enum {SIZE=Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize, SIZEF=SIZE};
public:
// define the number of components of your system, this is used outside
// to allocate the correct size of (dense) blocks with a FieldMatrix
enum {n=G::dimension};
enum {m=1};
// types for matrics, vectors and boundary conditions
typedef FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
typedef FieldVector<RT,m> VBlockType; // one entry in the global vectors
typedef FixedArray<BoundaryConditions::Flags,m> BCBlockType; // componentwise boundary conditions
//! Constructor
GroundwaterEquationLocalStiffness (const GroundwaterEquationParameters<G,RT>& params, bool procBoundaryAsDirichlet_=true)
GroundwaterEquationLocalStiffness (const GroundwaterEquationParameters<G,RT>& params,
bool procBoundaryAsDirichlet_=true)
: problem(params)
{
currentsize = 0;
this->currentsize = 0;
procBoundaryAsDirichlet = procBoundaryAsDirichlet_;
}
......@@ -111,20 +104,17 @@ namespace Dune
{
// extract some important parameters
Dune::NewGeometryType gt = e.geometry().type();
const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type& sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,k);
currentsize = sfs.size();
// std::cout << "Entity:" << std::endl;
// for (int i=0; i<e.template count<n>(); i++)
// std::cout << "corner i=" << i << " pos=" << e.geometry()[i] << std::endl;
const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type&
sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,k);
this->currentsize = sfs.size();
// clear assemble data
for (int i=0; i<sfs.size(); i++)
for (int i=0; i<this->currentsize; i++)
{
b[i] = 0;
bctype[i][0] = BoundaryConditions::neumann;
for (int j=0; j<sfs.size(); j++)
A[i][j] = 0;
this->b[i] = 0;
this->bctype[i][0] = BoundaryConditions::neumann;
for (int j=0; j<this->currentsize; j++)
this->A[i][j] = 0;
}
// Loop over all quadrature points and assemble matrix and right hand side
......@@ -133,18 +123,20 @@ namespace Dune
if (k>1) p=2*(k-1);
for (size_t g=0; g<Dune::QuadratureRules<DT,n>::rule(gt,p).size(); ++g) // run through all quadrature points
{
const Dune::FieldVector<DT,n>& local = Dune::QuadratureRules<DT,n>::rule(gt,p)[g].position(); // pos of integration point
Dune::FieldVector<DT,n> global = e.geometry().global(local); // ip in global coordinates
const Dune::FieldMatrix<DT,n,n> jac = e.geometry().jacobianInverseTransposed(local); // eval jacobian inverse
const Dune::FieldMatrix<DT,n,n> K = problem.K(global,e,local); // eval diffusion tensor
double weight = Dune::QuadratureRules<DT,n>::rule(gt,p)[g].weight(); // weight of quadrature point
DT detjac = e.geometry().integrationElement(local); // determinant of jacobian
RT q = problem.q(global,e,local); // source term
const Dune::FieldVector<DT,n>&
local = Dune::QuadratureRules<DT,n>::rule(gt,p)[g].position(); // pos of integration point
Dune::FieldVector<DT,n> global = e.geometry().global(local); // ip in global coordinates
const Dune::FieldMatrix<DT,n,n>
jac = e.geometry().jacobianInverseTransposed(local); // eval jacobian inverse
const Dune::FieldMatrix<DT,n,n> K = problem.K(global,e,local); // eval diffusion tensor
double weight = Dune::QuadratureRules<DT,n>::rule(gt,p)[g].weight(); // weight of quadrature point
DT detjac = e.geometry().integrationElement(local); // determinant of jacobian
RT q = problem.q(global,e,local); // source term
RT factor = weight*detjac;
// evaluate gradients at Gauss points
Dune::FieldVector<DT,n> grad[SIZE], temp, gv;
for (int i=0; i<sfs.size(); i++)
Dune::FieldVector<DT,n> grad[this->SIZE], temp, gv;
for (int i=0; i<this->currentsize; i++)
{
for (int l=0; l<n; l++)
temp[l] = sfs[i].evaluateDerivative(0,l,local);
......@@ -152,19 +144,19 @@ namespace Dune
jac.umv(temp,grad[i]); // transform gradient to global ooordinates
}
for (int i=0; i<sfs.size(); i++) // loop over test function number
for (int i=0; i<this->currentsize; i++) // loop over test function number
{
// rhs
b[i] += q*sfs[i].evaluateFunction(0,local)*factor;
this->b[i] += q*sfs[i].evaluateFunction(0,local)*factor;
// matrix
gv = 0; K.umv(grad[i],gv); // multiply with diffusion tensor
A[i][i] += (grad[i]*gv)*factor;
this->A[i][i] += (grad[i]*gv)*factor;
for (int j=0; j<i; j++)
{
RT t = (grad[j]*gv)*factor;
A[i][j] += t;
A[j][i] += t;
this->A[i][j] += t;
this->A[j][i] += t;
}
}
}
......@@ -191,29 +183,16 @@ namespace Dune
FieldVector<DT,n> global = it.intersectionGlobal().global(facelocal);
bctypeface = problem.bctype(global,e,local); // eval bctype
// std::cout << "=== Boundary found"
// << " facenumber=" << it.numberInSelf()
// << " quadpoint=" << g
// << " facelocal=" << facelocal
// << " local=" << local
// << " global=" << global
// << std::endl;
if (bctypeface!=BoundaryConditions::neumann) break;
RT J = problem.J(global,e,local);
double weightface = Dune::QuadratureRules<DT,n-1>::rule(gtface,p)[g].weight();
DT detjacface = it.intersectionGlobal().integrationElement(facelocal);
for (int i=0; i<sfs.size(); i++) // loop over test function number
if (bctype[i][0]==BoundaryConditions::neumann)
for (int i=0; i<this->currentsize; i++) // loop over test function number
if (this->bctype[i][0]==BoundaryConditions::neumann)
{
b[i] -= J*sfs[i].evaluateFunction(0,local)*weightface*detjacface;
// std::cout << "Neumann BC found, accumulating"
// << -J*sfs[i].evaluateFunction(0,local)*weightface*detjacface
// << std::endl;
// std::cout << "J=" << J << " shapef=" << sfs[i].evaluateFunction(0,local)
// << " weight=" << weightface
// << " detjac=" << detjacface << std::endl;
this->b[i] -= J*sfs[i].evaluateFunction(0,local)*weightface*detjacface;
}
}
if (bctypeface==BoundaryConditions::neumann) continue; // was a neumann face, go to next face
......@@ -226,22 +205,22 @@ namespace Dune
continue; // then it acts like homogeneous Neumann
// now handle exterior or interior Dirichlet boundary
for (int i=0; i<sfs.size(); i++) // loop over test function number
for (int i=0; i<this->currentsize; i++) // loop over test function number
{
if (sfs[i].codim()==0) continue; // skip interior dof
if (sfs[i].codim()==1) // handle face dofs
{
if (sfs[i].entity()==it.numberInSelf())
{
if (bctype[i][0]<bctypeface)
if (this->bctype[i][0]<bctypeface)
{
bctype[i][0] = bctypeface;
this->bctype[i][0] = bctypeface;
if (bctypeface==BoundaryConditions::process)
b[i] = 0;
this->b[i] = 0;
if (bctypeface==BoundaryConditions::dirichlet)
{
Dune::FieldVector<DT,n> global = e.geometry().global(sfs[i].position());
b[i] = problem.g(global,e,sfs[i].position());
this->b[i] = problem.g(global,e,sfs[i].position());
}
}
}
......@@ -251,15 +230,15 @@ namespace Dune
for (int j=0; j<ReferenceElements<DT,n>::general(gt).size(it.numberInSelf(),1,sfs[i].codim()); j++)
if (sfs[i].entity()==ReferenceElements<DT,n>::general(gt).subEntity(it.numberInSelf(),1,j,sfs[i].codim()))
{
if (bctype[i][0]<bctypeface)
if (this->bctype[i][0]<bctypeface)
{
bctype[i][0] = bctypeface;
this->bctype[i][0] = bctypeface;
if (bctypeface==BoundaryConditions::process)
b[i] = 0;
this->b[i] = 0;
if (bctypeface==BoundaryConditions::dirichlet)
{
Dune::FieldVector<DT,n> global = e.geometry().global(sfs[i].position());
b[i] = problem.g(global,e,sfs[i].position());
this->b[i] = problem.g(global,e,sfs[i].position());
}
}
}
......@@ -267,78 +246,10 @@ namespace Dune
}
}
// print contents of local stiffness matrix
void print (std::ostream& s, int width, int precision)
{
// set the output format
s.setf(std::ios_base::scientific, std::ios_base::floatfield);
int oldprec = s.precision();
s.precision(precision);
for (int i=0; i<currentsize; i++)
{
s << "FEM"; // start a new row
s << " "; // space in front of each entry
s.width(4); // set width for counter
s << i; // number of first entry in a line
for (int j=0; j<currentsize; j++)
{
s << " "; // space in front of each entry
s.width(width); // set width for each entry anew
s << A[i][j]; // yeah, the number !
}
s << " "; // space in front of each entry
s.width(width); // set width for each entry anew
s << b[i];
s << " "; // space in front of each entry
s.width(width); // set width for each entry anew
s << bctype[i][0];
s << std::endl; // start a new line
}
// reset the output format
s.precision(oldprec);
s.setf(std::ios_base::fixed, std::ios_base::floatfield);
}
//! access local stiffness matrix
/*! Access elements of the local stiffness matrix. Elements are
undefined without prior call to the assemble method.
*/
const MBlockType& mat (int i, int j)
{
return A[i][j];
}
//! access right hand side
/*! Access elements of the right hand side vector. Elements are
undefined without prior call to the assemble method.
*/
const VBlockType& rhs (int i)
{
return b[i];
}
//! access boundary condition for each dof
/*! Access boundary condition type for each degree of freedom. Elements are
undefined without prior call to the assemble method.
*/
const BCBlockType bc (int i) const
{
return bctype[i];
}
private:
// parameters given in constructor
const GroundwaterEquationParameters<G,RT>& problem;
bool procBoundaryAsDirichlet;
// assembled data
int currentsize;
MBlockType A[SIZE][SIZE];
VBlockType b[SIZE];
BCBlockType bctype[SIZE];
};
/** @} */
......
disc/operators/p1operator.hh
+ 185
− 77
View file @ fbe7cc69
......@@ -12,6 +12,7 @@
#include <stdio.h>
#include <stdlib.h>
#include "dune/common/timer.hh"
#include "dune/common/fvector.hh"
#include "dune/common/exceptions.hh"
#include "dune/common/exceptions.hh"
......@@ -426,6 +427,7 @@ namespace Dune
extraDOFs = 0;
// resize the S vector needed for detecting hanging nodes
watch.reset();
hanging.resize(vertexmapper.size());
std::vector<unsigned char> S(vertexmapper.size());
for (int i=0; i<vertexmapper.size(); i++)
......@@ -434,95 +436,106 @@ namespace Dune
hanging[i] = false;
}
// LOOP 1 : Prepare hanging node detection
// collect links of border vertices
Iterator eendit = indexset.template end<0,All_Partition>();
for (Iterator it = indexset.template begin<0,All_Partition>(); it!=eendit; ++it)
// do the following only if grid may have hanging nodes
if (Dune::Capabilities::template hasHangingNodes<G>::v)
{
Dune::NewGeometryType gt = it->geometry().type();
const typename Dune::ReferenceElementContainer<DT,n>::value_type&
refelem = ReferenceElements<DT,n>::general(gt);
// compute S value in vertex
for (int i=0; i<refelem.size(n); i++)
// LOOP 1 : Prepare hanging node detection
Iterator eendit = indexset.template end<0,All_Partition>();
for (Iterator it = indexset.template begin<0,All_Partition>(); it!=eendit; ++it)
{
int alpha = vertexmapper.template map<n>(*it,i);
if (S[alpha]>it->level()) S[alpha] = it->level(); // compute minimum
}
}
// LOOP 2 : second stage of detecting hanging nodes
for (Iterator it = indexset.template begin<0,All_Partition>(); it!=eendit; ++it)
{
Dune::NewGeometryType gt = it->geometry().type();
const typename Dune::ReferenceElementContainer<DT,n>::value_type&
refelem = ReferenceElements<DT,n>::general(gt);
Dune::NewGeometryType gt = it->geometry().type();
const typename Dune::ReferenceElementContainer<DT,n>::value_type&
refelem = ReferenceElements<DT,n>::general(gt);
// detect hanging nodes
IntersectionIterator endiit = it->iend();
for (IntersectionIterator iit = it->ibegin(); iit!=endiit; ++iit)
if (iit.neighbor())
// compute S value in vertex
for (int i=0; i<refelem.size(n); i++)
{
// check if neighbor is on lower level
const EEntityPointer outside = iit.outside();
if (it->level()<=outside->level()) continue;
int alpha = vertexmapper.template map<n>(*it,i);
if (S[alpha]>it->level()) S[alpha] = it->level(); // compute minimum
}
}
// loop over all vertices of this face
for (int j=0; j<refelem.size(iit.numberInSelf(),1,n); j++)
// LOOP 2 : second stage of detecting hanging nodes
for (Iterator it = indexset.template begin<0,All_Partition>(); it!=eendit; ++it)
{
Dune::NewGeometryType gt = it->geometry().type();
const typename Dune::ReferenceElementContainer<DT,n>::value_type&
refelem = ReferenceElements<DT,n>::general(gt);
// detect hanging nodes
IntersectionIterator endiit = it->iend();
for (IntersectionIterator iit = it->ibegin(); iit!=endiit; ++iit)
if (iit.neighbor())
{
int alpha = vertexmapper.template map<n>(*it,refelem.subEntity(iit.numberInSelf(),1,j,n));
if (S[alpha]==it->level())
hanging[alpha] = true;
}
}
}
// check if neighbor is on lower level
const EEntityPointer outside = iit.outside();
if (it->level()<=outside->level()) continue;
// local to global maps
int l2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
int fl2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
// loop over all vertices of this face
for (int j=0; j<refelem.size(iit.numberInSelf(),1,n); j++)
{
int alpha = vertexmapper.template map<n>(*it,refelem.subEntity(iit.numberInSelf(),1,j,n));
if (S[alpha]==it->level())
hanging[alpha] = true;
}
}
}
// LOOP 3 : determine additional links due to hanging nodes
for (Iterator it = indexset.template begin<0,All_Partition>(); it!=eendit; ++it)
{
Dune::NewGeometryType gt = it->geometry().type();
const typename Dune::ReferenceElementContainer<DT,n>::value_type&
refelem = ReferenceElements<DT,n>::general(gt);
// local to global maps
int l2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
int fl2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
// build local to global map
bool hasHangingNodes = false; // flag set to true if this element has hanging nodes
for (int i=0; i<refelem.size(n); i++)
// LOOP 3 : determine additional links due to hanging nodes
for (Iterator it = indexset.template begin<0,All_Partition>(); it!=eendit; ++it)
{
l2g[i] = vertexmapper.template map<n>(*it,i);
if (hanging[l2g[i]]) hasHangingNodes=true;
}
if (!hasHangingNodes) continue;
Dune::NewGeometryType gt = it->geometry().type();
const typename Dune::ReferenceElementContainer<DT,n>::value_type&
refelem = ReferenceElements<DT,n>::general(gt);
// handle father element if hanging nodes were detected
// get father element
const EEntityPointer father = it->father();
// build local to global map
bool hasHangingNodes = false; // flag set to true if this element has hanging nodes
for (int i=0; i<refelem.size(n); i++)
{
l2g[i] = vertexmapper.template map<n>(*it,i);
if (hanging[l2g[i]]) hasHangingNodes=true;
}
if (!hasHangingNodes) continue;
// build local to global map for father
for (int i=0; i<refelem.size(n); i++)
fl2g[i] = vertexmapper.template map<n>(*father,i);
// handle father element if hanging nodes were detected
// get father element
const EEntityPointer father = it->father();
// a map that inverts l2g
std::map<int,int> g2l;
for (int i=0; i<refelem.size(n); i++)
g2l[l2g[i]] = i;
// build local to global map for father
for (int i=0; i<refelem.size(n); i++)
fl2g[i] = vertexmapper.template map<n>(*father,i);
// connect all fine nodes to all coarse nodes
for (int i=0; i<refelem.size(n); i++) // nodes in *it
for (int j=0; j<refelem.size(n); j++) // nodes in *father
if (g2l.find(fl2g[j])==g2l.end())
{
links.insert(P1OperatorLink(l2g[i],fl2g[j]));
links.insert(P1OperatorLink(fl2g[j],l2g[i]));
// std::cout << "link" << " gi=" << l2g[i] << " gj=" << fl2g[j] << std::endl;
// std::cout << "link" << " gi=" << fl2g[j] << " gj=" << l2g[i] << std::endl;
}
// a map that inverts l2g
std::map<int,int> g2l;
for (int i=0; i<refelem.size(n); i++)
g2l[l2g[i]] = i;
// connect all fine nodes to all coarse nodes
for (int i=0; i<refelem.size(n); i++) // nodes in *it
for (int j=0; j<refelem.size(n); j++) // nodes in *father
if (g2l.find(fl2g[j])==g2l.end())
{
links.insert(P1OperatorLink(l2g[i],fl2g[j]));
links.insert(P1OperatorLink(fl2g[j],l2g[i]));
// std::cout << "link" << " gi=" << l2g[i] << " gj=" << fl2g[j] << std::endl;
// std::cout << "link" << " gi=" << fl2g[j] << " gj=" << l2g[i] << std::endl;
}
}
}
// count hanging nodes
hangingnodes = 0;
for (int i=0; i<vertexmapper.size(); i++)
if (hanging[i]) hangingnodes++;
std::cout << "=== P1OperatorBase hanging node detection + add links " << watch.elapsed() << std::endl;
// compute additional links due to extended overlap
watch.reset();
if (extendOverlap)
{
// set of neighbors in global ids for border vertices
......@@ -545,11 +558,7 @@ namespace Dune
// Note: links contains now also connections that are standard.
// So below we have throw out these connections again!
// count hanging nodes, can be used whether grid has hanging nodes at all
hangingnodes = 0;
for (int i=0; i<vertexmapper.size(); i++)
if (hanging[i]) hangingnodes++;
std::cout << "=== P1OperatorBase parallel extend overlap " << watch.elapsed() << std::endl;
return true;
}
......@@ -675,6 +684,7 @@ namespace Dune
for (int i=0; i<allmapper.size(); i++) visited[i] = false;
// LOOP 4 : Compute row sizes
watch.reset();
Iterator eendit = is.template end<0,All_Partition>();
for (Iterator it = is.template begin<0,All_Partition>(); it!=eendit; ++it)
{
......@@ -729,11 +739,13 @@ namespace Dune
// now the row sizes have been set
A.endrowsizes();
std::cout << "=== P1OperatorBase compute row sizes " << watch.elapsed() << std::endl;
// clear the flags for the next round, actually that is not necessary because addindex takes care of this
for (int i=0; i<allmapper.size(); i++) visited[i] = false;
// LOOP 5 : insert the nonzeros
watch.reset();
for (Iterator it = is.template begin<0,All_Partition>(); it!=eendit; ++it)
{
Dune::NewGeometryType gt = it->geometry().type();
......@@ -783,6 +795,7 @@ namespace Dune
// now the matrix is ready for use
A.endindices();
std::cout << "=== P1OperatorBase index insertion " << watch.elapsed() << std::endl;
// delete additional links
links.clear();
......@@ -820,6 +833,7 @@ namespace Dune
}
protected:
Timer watch;
const G& grid;
const IS& is;
LC lc;
......@@ -944,8 +958,10 @@ namespace Dune
DUNE_THROW(MathError,"P1OperatorAssembler::assemble(): size mismatch");
// clear global stiffness matrix and right hand side
this->watch.reset();
this->A = 0;
*f = 0;
std::cout << "=== P1OperatorBase clear matrix " << this->watch.elapsed() << std::endl;
// allocate flag vector to hold flags for essential boundary conditions
std::vector<BCBlockType> essential(this->vertexmapper.size());
......@@ -1423,6 +1439,98 @@ namespace Dune
};
/*! @brief Base class for local P1 assemblers
*/
template<class DT, class RT, int n, int m>
class P1LocalStiffness
{
protected:
// some other sizes
enum {SIZE=Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize, SIZEF=SIZE};
public:
// types for matrics, vectors and boundary conditions
typedef FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
typedef FieldVector<RT,m> VBlockType; // one entry in the global vectors
typedef FixedArray<BoundaryConditions::Flags,m> BCBlockType; // componentwise boundary conditions
//! Constructor
P1LocalStiffness ()
{
currentsize = 0;
}
// print contents of local stiffness matrix
void print (std::ostream& s, int width, int precision)
{
// set the output format
s.setf(std::ios_base::scientific, std::ios_base::floatfield);
int oldprec = s.precision();
s.precision(precision);
for (int i=0; i<currentsize; i++)
{
s << "FEM"; // start a new row
s << " "; // space in front of each entry
s.width(4); // set width for counter
s << i; // number of first entry in a line
for (int j=0; j<currentsize; j++)
{
s << " "; // space in front of each entry
s.width(width); // set width for each entry anew
s << A[i][j]; // yeah, the number !
}
s << " "; // space in front of each entry
s.width(width); // set width for each entry anew
s << b[i];
s << " "; // space in front of each entry
s.width(width); // set width for each entry anew
s << bctype[i][0];
s << std::endl; // start a new line
}
// reset the output format
s.precision(oldprec);
s.setf(std::ios_base::fixed, std::ios_base::floatfield);
}
//! access local stiffness matrix
/*! Access elements of the local stiffness matrix. Elements are
undefined without prior call to the assemble method.
*/
const MBlockType& mat (int i, int j)
{
return A[i][j];
}
//! access right hand side
/*! Access elements of the right hand side vector. Elements are
undefined without prior call to the assemble method.
*/
const VBlockType& rhs (int i)
{
return b[i];
}
//! access boundary condition for each dof
/*! Access boundary condition type for each degree of freedom. Elements are
undefined without prior call to the assemble method.
*/
const BCBlockType bc (int i) const
{
return bctype[i];
}
protected:
// assembled data
int currentsize;
MBlockType A[SIZE][SIZE];
VBlockType b[SIZE];
BCBlockType bctype[SIZE];
};
/** @} */
}
......
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