Skip to content
Snippets Groups Projects
Commit 1428d2d3 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Implement GlobalP2Mapper for cube grids without dune-parmg

Unfortunately, this code is used even in sequential situations,
but previously the code only supported 1d- and 2d grids.
And as we still support installations without dune-parmg,
this commit adds an implementation for hypercube grids
of all dimensions.
parent 78cc500c
No related branches found
No related tags found
1 merge request!146Implement GlobalP2Mapper for 3d grids even without dune-parmg
#ifndef DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH #ifndef DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
#define DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH #define DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
#include <memory>
#include <dune/functions/functionspacebases/lagrangebasis.hh> #include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/gfe/parallel/globalmapper.hh> #include <dune/gfe/parallel/globalmapper.hh>
...@@ -22,96 +24,44 @@ namespace Dune { ...@@ -22,96 +24,44 @@ namespace Dune {
p2Mapper_(gridView) p2Mapper_(gridView)
{ {
#if !HAVE_DUNE_PARMG #if !HAVE_DUNE_PARMG
static_assert(GridView::dimension<=2, "Only implemented for one- and two-dimensional grids"); if (gridView.size(GeometryTypes::simplex(GridView::dimension))>=1)
DUNE_THROW(NotImplemented, "GlobalP2Mapper only works for grids without simplex elements!");
// Indexed by *co*dimension!
std::array<std::unique_ptr<GlobalIndexSet<GridView> >, GridView::dimension+1> globalIndexSets;
if (gridView.size(GeometryTypes::triangle)>1) for (int i=0; i<=GridView::dimension; ++i)
DUNE_THROW(NotImplemented, "GlobalP2Mapper only works for quad grids!"); globalIndexSets[i] = std::make_unique<GlobalIndexSet<GridView> >(gridView, i);
GlobalIndexSet<GridView> globalVertexIndex(gridView,GridView::dimension); // total number of degrees of freedom
GlobalIndexSet<GridView> globalElementIndex(gridView,0); this->size_ = 0;
for (int i=0; i<=GridView::dimension; ++i)
this->size_ += globalIndexSets[i]->size(i);
auto localView = p2Mapper_.localView(); auto localView = p2Mapper_.localView();
if (GridView::dimension==1) // Loop over all elements, as a means to loop over all degrees of freedom
for (const auto& element : elements(gridView))
{ {
// total number of degrees of freedom localView.bind(element);
this->size_ = globalVertexIndex.size(GridView::dimension) + globalElementIndex.size(0);
// Determine // Loop over all local degrees of freedom
for (const auto& element : elements(gridView)) for (size_t i=0; i<localView.size(); i++)
{ {
localView.bind(element); auto codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
auto entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
// Loop over all local degrees of freedom
for (size_t i=0; i<localView.size(); i++)
{
int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
int entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
auto localIndex = localView.index(i);
int globalIndex;
switch (codim)
{
case 1 : // vertex dofs
globalIndex = globalVertexIndex.index(element.template subEntity<1>(entity));
break;
case 0 : // element dofs
globalIndex = globalElementIndex.index(element.template subEntity<0>(entity))
+ globalVertexIndex.size(1);
break;
default :
DUNE_THROW(Dune::Exception, "Impossible codimension!");
}
localGlobalMap_[localIndex] = globalIndex;
}
}
} auto localIndex = localView.index(i);
else
{
GlobalIndexSet<GridView> globalEdgeIndex(gridView,1);
// total number of degrees of freedom // Compute a global index for this degree of freedom
this->size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0); int globalIndex = 0;
// Determine for (int j=codim+1; j<=GridView::dimension; ++j)
for (const auto& element : elements(gridView)) globalIndex += globalIndexSets[j]->size(j);
{
localView.bind(element); globalIndex += globalIndexSets[codim]->subIndex(element, entity, codim);
// Loop over all local degrees of freedom localGlobalMap_[localIndex] = globalIndex;
for (size_t i=0; i<localView.size(); i++)
{
int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
int entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
auto localIndex = localView.index(i);
int globalIndex;
switch (codim)
{
case 2 : // vertex dofs
globalIndex = globalVertexIndex.index(element.template subEntity<GridView::dimension>(entity));
break;
case 1 : // edge dofs
globalIndex = globalEdgeIndex.index(element.template subEntity<1>(entity)) + globalVertexIndex.size(2);
break;
case 0 : // element dofs
globalIndex = globalElementIndex.index(element.template subEntity<0>(entity))
+ globalEdgeIndex.size(1)
+ globalVertexIndex.size(2);
break;
default :
DUNE_THROW(Dune::Exception, "Impossible codimension!");
}
localGlobalMap_[localIndex] = globalIndex;
}
} }
} }
#endif #endif
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment