Commit 8eeb64a1 authored by Carsten Gräser's avatar Carsten Gräser

Merge branch 'feature/fix-pq0nodalbasis-2.4' into 'releases/2.4-compatible'

Feature/fix pq0nodalbasis 2.4

Cherry-picked Elias patches from !6 for #1 .

See merge request !28
parents 0fd59758 4bc23f03
......@@ -56,14 +56,22 @@ public:
// Precompute the number of dofs per entity type
const static int dofsPerEdge = k-1;
const static int dofsPerTriangle = (k-1)*(k-2)/2;
const static int dofsPerQuad = (k-1)*(k-1);
const static int dofsPerTetrahedron = ((k-3)*(k-2)*(k-1)/6 > 0) ? (k-3)*(k-2)*(k-1)/6 : 0;
const static int dofsPerPrism = (k-1)*(k-1)*(k-2)/2;
const static int dofsPerHexahedron = (k-1)*(k-1)*(k-1);
const static int dofsPerPyramid = ((k-2)*(k-1)*(2*k-3))/6;
const static int dofsPerVertex =
k == 0 ? (dim == 0 ? 1 : 0) : 1;
const static int dofsPerEdge =
k == 0 ? (dim == 1 ? 1 : 0) : k-1;
const static int dofsPerTriangle =
k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-2)/2;
const static int dofsPerQuad =
k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-1);
const static int dofsPerTetrahedron =
k == 0 ? (dim == 3 ? 1 : 0) : (k-3)*(k-2)*(k-1)/6;
const static int dofsPerPrism =
k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-2)/2;
const static int dofsPerHexahedron =
k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-1);
const static int dofsPerPyramid =
k == 0 ? (dim == 3 ? 1 : 0) : (k-2)*(k-1)*(2*k-3)/6;
template<class TP>
using Node = PQkNode<GV, k, size_type, TP>;
......@@ -85,7 +93,7 @@ public:
void initializeIndices()
{
vertexOffset_ = 0;
edgeOffset_ = vertexOffset_ + gridView_.size(dim);
edgeOffset_ = vertexOffset_ + dofsPerVertex * gridView_.size(dim);
triangleOffset_ = edgeOffset_ + dofsPerEdge * gridView_.size(dim-1);
GeometryType triangle;
......@@ -135,14 +143,18 @@ public:
switch (dim)
{
case 1:
return gridView_.size(1) + dofsPerEdge*gridView_.size(0);
return dofsPerVertex * gridView_.size(dim)
+ dofsPerEdge*gridView_.size(dim-1);
case 2:
{
GeometryType triangle, quad;
triangle.makeTriangle();
quad.makeQuadrilateral();
return gridView_.size(dim) + dofsPerEdge*gridView_.size(1)
+ dofsPerTriangle*gridView_.size(triangle) + dofsPerQuad*gridView_.size(quad);
return dofsPerVertex * gridView_.size(dim)
+ dofsPerEdge * gridView_.size(dim-1)
+ dofsPerTriangle * gridView_.size(triangle)
+ dofsPerQuad * gridView_.size(quad);
}
case 3:
{
......@@ -153,10 +165,14 @@ public:
pyramid.makePyramid();
prism.makePrism();
hexahedron.makeCube(3);
return gridView_.size(dim) + dofsPerEdge*gridView_.size(2)
+ dofsPerTriangle*gridView_.size(triangle) + dofsPerQuad*gridView_.size(quad)
+ dofsPerTetrahedron*gridView_.size(tetrahedron) + dofsPerPyramid*gridView_.size(pyramid)
+ dofsPerPrism*gridView_.size(prism) + dofsPerHexahedron*gridView_.size(hexahedron);
return dofsPerVertex * gridView_.size(dim)
+ dofsPerEdge * gridView_.size(dim-1)
+ dofsPerTriangle * gridView_.size(triangle)
+ dofsPerQuad * gridView_.size(quad)
+ dofsPerTetrahedron * gridView_.size(tetrahedron)
+ dofsPerPyramid * gridView_.size(pyramid)
+ dofsPerPrism * gridView_.size(prism)
+ dofsPerHexahedron * gridView_.size(hexahedron);
}
}
DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
......@@ -306,13 +322,15 @@ public:
size_t dofDim = dim - localKey.codim();
if (dofDim==0) { // vertex dof
return { gridIndexSet.subIndex(element,localKey.subEntity(),dim) };
return {{ gridIndexSet.subIndex(element,localKey.subEntity(),dim) }};
}
if (dofDim==1)
{ // edge dof
if (dim==1) // element dof -- any local numbering is fine
return { nodeFactory_->edgeOffset_ + (k-1)*gridIndexSet.subIndex(element,0,0) + localKey.index() };
if (dim==1) // element dof -- any local numbering is fine
return {{ nodeFactory_->edgeOffset_
+ nodeFactory_->dofsPerEdge * gridIndexSet.subIndex(element,0,0)
+ localKey.index() }};
else
{
const Dune::ReferenceElement<double,dim>& refElement
......@@ -323,9 +341,13 @@ public:
size_t v0 = gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
size_t v1 = gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
bool flip = (v0 > v1);
return { (flip)
? nodeFactory_->edgeOffset_ + (k-1)*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) + (k-2)-localKey.index()
: nodeFactory_->edgeOffset_ + (k-1)*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) + localKey.index() };
return {{ (flip)
? nodeFactory_->edgeOffset_
+ nodeFactory_->dofsPerEdge*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())
+ (nodeFactory_->dofsPerEdge-1)-localKey.index()
: nodeFactory_->edgeOffset_
+ nodeFactory_->dofsPerEdge*gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())
+ localKey.index() }};
}
}
......@@ -336,12 +358,12 @@ public:
if (element.type().isTriangle())
{
const int interiorLagrangeNodesPerTriangle = (k-1)*(k-2)/2;
return { nodeFactory_->triangleOffset_ + interiorLagrangeNodesPerTriangle*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->triangleOffset_ + interiorLagrangeNodesPerTriangle*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
}
else if (element.type().isQuadrilateral())
{
const int interiorLagrangeNodesPerQuadrilateral = (k-1)*(k-1);
return { nodeFactory_->quadrilateralOffset_ + interiorLagrangeNodesPerQuadrilateral*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->quadrilateralOffset_ + interiorLagrangeNodesPerQuadrilateral*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
}
else
DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
......@@ -356,7 +378,7 @@ public:
if (k==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
DUNE_THROW(Dune::NotImplemented, "PQkNodalBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
return { nodeFactory_->triangleOffset_ + gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) };
return {{ nodeFactory_->triangleOffset_ + gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()) }};
}
}
......@@ -365,13 +387,13 @@ public:
if (dim==3) // element dof -- any local numbering is fine
{
if (element.type().isTetrahedron())
return { nodeFactory_->tetrahedronOffset_ + NodeFactory::dofsPerTetrahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->tetrahedronOffset_ + NodeFactory::dofsPerTetrahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else if (element.type().isHexahedron())
return { nodeFactory_->hexahedronOffset_ + NodeFactory::dofsPerHexahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->hexahedronOffset_ + NodeFactory::dofsPerHexahedron*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else if (element.type().isPrism())
return { nodeFactory_->prismOffset_ + NodeFactory::dofsPerPrism*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->prismOffset_ + NodeFactory::dofsPerPrism*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else if (element.type().isPyramid())
return { nodeFactory_->pyramidOffset_ + NodeFactory::dofsPerPyramid*gridIndexSet.subIndex(element,0,0) + localKey.index() };
return {{ nodeFactory_->pyramidOffset_ + NodeFactory::dofsPerPyramid*gridIndexSet.subIndex(element,0,0) + localKey.index() }};
else
DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids");
} else
......
......@@ -270,6 +270,10 @@ void testOnStructuredGrid()
if (dim<3) // Currently not implemented for dim >= 3
testScalarBasis(pq4Basis, true);
// Test PQkNodalBasis for the corner case k == 0
PQkNodalBasis<GridView, 0> pq0Basis(gridView);
testScalarBasis(pq0Basis, true);
// Test LagrangeDGBasis for k==1
LagrangeDGBasis<GridView, 1> lagrangeDG1Basis(gridView);
testScalarBasis(lagrangeDG1Basis, true);
......
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