4 #ifndef DUNE_PDELAB_CONFORMINGCONSTRAINTS_HH
5 #define DUNE_PDELAB_CONFORMINGCONSTRAINTS_HH
10 #include <dune/common/exceptions.hh>
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/type.hh>
15 #include <dune/grid/common/grid.hh>
17 #include <dune/localfunctions/common/interfaceswitch.hh>
19 #include <dune/typetree/typetree.hh>
51 template<
typename P,
typename IG,
typename LFS,
typename T>
52 void boundary (
const P& param,
const IG&
ig,
const LFS& lfs, T& trafo)
const
54 typedef FiniteElementInterfaceSwitch<
55 typename LFS::Traits::FiniteElementType
57 typedef FieldVector<
typename IG::ctype, IG::dimension-1> FaceCoord;
59 const int face = ig.indexInInside();
62 Dune::GeometryType gt = ig.inside()->type();
63 typedef typename IG::ctype DT;
64 const int dim = IG::Entity::Geometry::dimension;
65 const Dune::ReferenceElement<DT,dim>& refelem = Dune::ReferenceElements<DT,dim>::general(gt);
67 const Dune::ReferenceElement<DT,dim-1> &
68 face_refelem = Dune::ReferenceElements<DT,dim-1>::general(ig.geometry().type());
71 typename T::RowType empty;
73 const FaceCoord testpoint = face_refelem.position(0,0);
76 if (!param.isDirichlet(ig,testpoint))
80 i<std::size_t(FESwitch::coefficients(lfs.finiteElement()).size());
85 FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
87 if (codim==0)
continue;
89 for (
int j=0; j<refelem.size(face,1,codim); j++){
91 if (static_cast<int>(FESwitch::coefficients(lfs.finiteElement()).
92 localKey(i).subEntity())
93 == refelem.subEntity(face,1,j,codim))
94 trafo[lfs.dofIndex(i)] = empty;
112 template<
typename IG,
typename LFS,
typename T>
115 typedef FiniteElementInterfaceSwitch<
116 typename LFS::Traits::FiniteElementType
120 const int face = ig.indexInInside();
123 Dune::GeometryType gt = ig.inside()->type();
124 typedef typename IG::ctype DT;
125 const int dim = IG::Entity::Geometry::dimension;
127 const Dune::ReferenceElement<DT,dim>& refelem = Dune::ReferenceElements<DT,dim>::general(gt);
130 typename T::RowType empty;
133 for (
size_t i=0; i<FESwitch::coefficients(lfs.finiteElement()).size();
138 FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
140 if (codim==0)
continue;
142 for (
int j=0; j<refelem.size(face,1,codim); j++)
143 if (FESwitch::coefficients(lfs.finiteElement()).localKey(i).
144 subEntity() == std::size_t(refelem.subEntity(face,1,j,codim)))
145 trafo[lfs.dofIndex(i)] = empty;
151 template<
typename GV>
163 template<
typename EG,
typename LFS,
typename T>
164 void volume (
const EG&
eg,
const LFS& lfs, T& trafo)
const
166 typedef FiniteElementInterfaceSwitch<
167 typename LFS::Traits::FiniteElementType
171 if (eg.entity().partitionType()==Dune::InteriorEntity)
174 typedef typename FESwitch::Coefficients Coefficients;
175 const Coefficients& coeffs = FESwitch::coefficients(lfs.finiteElement());
178 typename T::RowType empty;
180 const ReferenceElement<typename GV::ctype,GV::dimension>& ref_el =
181 ReferenceElements<typename GV::ctype,GV::dimension>::general(eg.entity().type());
184 for (
size_t i = 0; i < coeffs.size(); ++i)
186 size_t codim = coeffs.localKey(i).codim();
187 size_t sub_entity = coeffs.localKey(i).subEntity();
189 size_t entity_index = _gv.indexSet().subIndex(eg.entity(),sub_entity,codim);
190 size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(sub_entity,codim));
192 size_t index = _gt_offsets[gt_index] + entity_index;
196 trafo[lfs.dofIndex(i)] = empty;
204 std::fill(_gt_offsets.begin(),_gt_offsets.end(),0);
206 typedef std::vector<GeometryType> GTVector;
208 for (
size_t codim = 0; codim <= GV::dimension; ++codim)
210 if (gfs.ordering().contains(codim))
212 const GTVector& geom_types = _gv.indexSet().geomTypes(codim);
213 for (GTVector::const_iterator it = geom_types.begin(),
214 end = geom_types.end();
217 _gt_offsets[GlobalGeometryTypeIndex::index(*it) + 1] = _gv.indexSet().size(*it);
221 std::partial_sum(_gt_offsets.begin(),_gt_offsets.end(),_gt_offsets.begin());
223 _ghosts.assign(_gt_offsets.back(),
true);
225 typedef typename GV::template Codim<0>::
226 template Partition<Interior_Partition>::Iterator Iterator;
228 for(Iterator it = _gv.template begin<0, Interior_Partition>(),
229 end = _gv.template end<0, Interior_Partition>();
233 const ReferenceElement<typename GV::ctype,GV::dimension>& ref_el =
234 ReferenceElements<typename GV::ctype,GV::dimension>::general(it->type());
236 for (
size_t codim = 0; codim <= GV::dimension; ++codim)
237 if (gfs.ordering().contains(codim))
239 for (
int i = 0; i < ref_el.size(codim); ++i)
241 size_t entity_index = _gv.indexSet().subIndex(*it,i,codim);
242 size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(i,codim));
243 size_t index = _gt_offsets[gt_index] + entity_index;
245 _ghosts[index] =
false;
254 , _rank(gv.comm().rank())
255 , _gt_offsets(GlobalGeometryTypeIndex::size(GV::dimension) + 1)
262 std::vector<bool> _ghosts;
263 std::vector<size_t> _gt_offsets;
Definition: conforming.hh:42
Definition: conforming.hh:41
void volume(const EG &eg, const LFS &lfs, T &trafo) const
volume constraints
Definition: conforming.hh:164
Definition: conforming.hh:39
Definition: conforming.hh:40
static const int dim
Definition: adaptivity.hh:82
void processor(const IG &ig, const LFS &lfs, T &trafo) const
processor constraints
Definition: conforming.hh:113
Definition: conforming.hh:155
const IG & ig
Definition: common/constraints.hh:146
NonoverlappingConformingDirichletConstraints(const GV &gv)
Definition: conforming.hh:252
Definition: conforming.hh:104
void compute_ghosts(const GFS &gfs)
Definition: conforming.hh:202
extend conforming constraints class by processor boundary
Definition: conforming.hh:152
const EG & eg
Definition: common/constraints.hh:277
void boundary(const P ¶m, const IG &ig, const LFS &lfs, T &trafo) const
boundary constraints
Definition: conforming.hh:52
extend conforming constraints class by processor boundary
Definition: conforming.hh:101
Dirichlet Constraints construction.
Definition: conforming.hh:36