4#ifndef DUNE_PDELAB_CONSTRAINTS_CONFORMING_HH
5#define DUNE_PDELAB_CONSTRAINTS_CONFORMING_HH
12#include <dune/geometry/referenceelements.hh>
17#include <dune/localfunctions/common/interfaceswitch.hh>
19#include <dune/typetree/typetree.hh>
21#include <dune/pdelab/common/geometrywrapper.hh>
22#include <dune/pdelab/constraints/common/constraintsparameters.hh>
23#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
40 enum { doBoundary =
true };
41 enum { doProcessor =
false };
42 enum { doSkeleton =
false };
43 enum { doVolume =
false };
52 template<
typename P,
typename IG,
typename LFS,
typename T>
53 void boundary (
const P& param,
const IG& ig,
const LFS& lfs, T& trafo)
const
56 typename LFS::Traits::FiniteElementType
60 const int face = ig.indexInInside();
67 typename T::RowType
empty;
69 const FaceCoord testpoint = face_refelem.position(0,0);
72 if (!param.isDirichlet(ig,testpoint))
76 i<std::size_t(FESwitch::coefficients(lfs.finiteElement()).size());
81 FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
83 if (codim==0)
continue;
85 for (
int j=0; j<refelem.size(face,1,codim); j++){
87 if (
static_cast<int>(FESwitch::coefficients(lfs.finiteElement()).
88 localKey(i).subEntity())
89 == refelem.subEntity(face,1,j,codim))
90 trafo[lfs.dofIndex(i)] =
empty;
100 enum { doProcessor =
true };
108 template<
typename IG,
typename LFS,
typename T>
109 void processor (
const IG& ig,
const LFS& lfs, T& trafo)
const
112 typename LFS::Traits::FiniteElementType
116 const int face = ig.indexInInside();
121 typename T::RowType
empty;
124 for (
size_t i=0; i<FESwitch::coefficients(lfs.finiteElement()).size();
129 FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
131 if (codim==0)
continue;
133 for (
int j=0; j<refelem.size(face,1,codim); j++)
134 if (FESwitch::coefficients(lfs.finiteElement()).localKey(i).
135 subEntity() == std::size_t(refelem.subEntity(face,1,j,codim)))
136 trafo[lfs.dofIndex(i)] =
empty;
142 template<
typename GV>
146 enum { doVolume =
true };
154 template<
typename P,
typename EG,
typename LFS,
typename T>
155 void volume (
const P& param,
const EG& eg,
const LFS& lfs, T& trafo)
const
158 typename LFS::Traits::FiniteElementType
161 auto& entity = eg.entity();
167 typedef typename FESwitch::Coefficients Coefficients;
168 const Coefficients& coeffs = FESwitch::coefficients(lfs.finiteElement());
171 typename T::RowType
empty;
176 for (
size_t i = 0; i < coeffs.size(); ++i)
178 size_t codim = coeffs.localKey(i).codim();
179 size_t sub_entity = coeffs.localKey(i).subEntity();
181 size_t entity_index = _gv.indexSet().subIndex(entity,sub_entity,codim);
184 size_t index = _gt_offsets[gt_index] + entity_index;
188 trafo[lfs.dofIndex(i)] =
empty;
194 void compute_ghosts (
const GFS& gfs)
196 std::fill(_gt_offsets.begin(),_gt_offsets.end(),0);
198 for (
size_t codim = 0; codim <= GV::dimension; ++codim)
200 if (gfs.ordering().contains(codim))
202 for (
auto gt : _gv.indexSet().types(codim))
207 std::partial_sum(_gt_offsets.begin(),_gt_offsets.end(),_gt_offsets.begin());
209 _ghosts.assign(_gt_offsets.back(),
true);
211 typedef typename GV::template Codim<0>::
212 template Partition<Interior_Partition>::Iterator Iterator;
214 for(Iterator it = _gv.template begin<0, Interior_Partition>(),
215 end = _gv.template end<0, Interior_Partition>();
223 for (
size_t codim = 0; codim <= GV::dimension; ++codim)
224 if (gfs.ordering().contains(codim))
226 for (
int i = 0; i < ref_el.size(codim); ++i)
228 size_t entity_index = _gv.indexSet().subIndex(entity,i,codim);
230 size_t index = _gt_offsets[gt_index] + entity_index;
232 _ghosts[index] =
false;
239 NonoverlappingConformingDirichletConstraints(
const GV& gv)
241 , _rank(gv.comm().rank())
242 , _gt_offsets(GlobalGeometryTypeIndex::
size(GV::dimension) + 1)
249 std::vector<bool> _ghosts;
250 std::vector<size_t> _gt_offsets;
vector space out of a tensor product of fields.
Definition: fvector.hh:91
static constexpr std::size_t index(const GeometryType >)
Compute the index for the given geometry type over all dimensions.
Definition: typeindex.hh:138
Base class for leaf nodes in a dune-typetree.
Definition: leafnode.hh:28
Different resources needed by all grid implementations.
A few common exception classes.
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr std::bool_constant<(sizeof...(II)==0)> empty(std::integer_sequence< T, II... >)
Checks whether the sequence is empty.
Definition: integersequence.hh:80
Switch for uniform treatment of finite element with either the local or the global interface.
Definition: interfaceswitch.hh:30
A unique label for each type of element that can occur in a grid.