[Dune] LeafMultipleCodimMultipleGeomTypeMapper: question concerning DOF mapping (Dune 1.2.1 stable)

[Jan Gerno]

Hello,

I am trying to write some simple Poisson solver with P2 elements in Dune. I
thought, this could be easily done using
a LeafMultipleCodimMultipleGeomTypeMapper. I wrote some small helper class
to map from local element indices to global indices in my linear system.

The problem is: the mapping of the DOF's is all messed up. I think the codim
2 entities (nodes) are not mapped correctly. I am pretty sure, that the bug
hides somewhere in my get(Face)DOFMapping routines.

I have no idea what I could have done wrong here. I am trying to fix this
for about 3 days now and did not make significant progress. Does anyone here
see the problem? Any help is greatly appreciated! I am using Dune 1.2.1 and
a ALUSimplex2D mesh.

Kind regards,

J. Gerno




==== CODE SNIPPET THAT PRODUCES THE WRONG MAPPING
====================================================================

template<typename G>
class P22DSpace
{
typedef G GridType;
typedef typename GridType::ctype RealType;
typedef typename GridType::LeafGridView::template Codim<0>::Iterator::Entity
EntityType;
typedef Dune::P22DLocalFiniteElement<RealType, RealType> FiniteElementType;
typedef Dune::P21DLocalFiniteElement<RealType, RealType>
FaceFiniteElementType;
 Dune::ReferenceSimplex<double, 2> ref;
 template<int D>
struct Layout
{
bool contains(Dune::GeometryType gt)
{
return gt.dim() == D-1 || gt.dim() == D-2;
}
};
 GridType const& grid;
unsigned dofoff;
FiniteElementType localfe;
FaceFiniteElementType facefe;
Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType, Layout> mapper;
 public:
P22DSpace(GridType const& grid, unsigned dofoff = 0) :
grid(grid), dofoff(dofoff), mapper(grid)
{
}
 FiniteElementType const& localFiniteElement() const
{
return localfe;
}
 FaceFiniteElementType const& localFaceFiniteElement() const
{
return facefe;
}
 size_t getDOFMapping(EntityType const& e, IndexVector& dofno) const
{
dofno.resize(6);
for(int i=0; i<3; i++)
dofno[i] = dofoff + mapper.template map<2>(e, ref.subEntity(0,0,i,2));
for(int i=0; i<3; i++)
dofno[i+3] = dofoff + mapper.template map<1>(e, ref.subEntity(0,0,i,1));
return dofno.size();
}
 size_t getFaceDOFMapping(EntityType const& e, unsigned index, IndexVector&
dofno) const
{
dofno.resize(3);
dofno[0] = dofoff + mapper.template map<2>(e, ref.subEntity(index,1,0,2));
// node 1
dofno[1] = dofoff + mapper.template map<2>(e, ref.subEntity(index,1,1,2));
// node 2
dofno[2] = dofoff + mapper.template map<1>(e, ref.subEntity(0,0,index,1));
// edge dof
return dofno.size();
}
 unsigned size() const
{
return mapper.size();
}
 static unsigned const order()
{
return 2;
}
 static unsigned const dim()
{
return 1;
}
};
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20091026/e58ec2c2/attachment.htm>