dune-pdelab  2.4-dev
conforming.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_CONFORMINGCONSTRAINTS_HH
5 #define DUNE_PDELAB_CONFORMINGCONSTRAINTS_HH
6 
7 #include <cstddef>
8 #include <algorithm>
9 
10 #include <dune/common/exceptions.hh>
11 
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/type.hh>
14 
15 #include <dune/grid/common/grid.hh>
16 
17 #include <dune/localfunctions/common/interfaceswitch.hh>
18 
19 #include <dune/typetree/typetree.hh>
20 
26 
27 namespace Dune {
28  namespace PDELab {
29 
33 
35  // works in any dimension and on all element types
37  {
38  public:
39  enum { doBoundary = true };
40  enum { doProcessor = false };
41  enum { doSkeleton = false };
42  enum { doVolume = false };
43 
45 
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
53  {
54  typedef FiniteElementInterfaceSwitch<
55  typename LFS::Traits::FiniteElementType
56  > FESwitch;
57  typedef FieldVector<typename IG::ctype, IG::dimension-1> FaceCoord;
58 
59  const int face = ig.indexInInside();
60 
61  // find all local indices of this face
62  Dune::GeometryType gt = ig.inside().type();
63  typedef typename IG::ctype DT;
64  const int dim = IG::Entity::dimension;
65  const Dune::ReferenceElement<DT,dim>& refelem = Dune::ReferenceElements<DT,dim>::general(gt);
66 
67  const Dune::ReferenceElement<DT,dim-1> &
68  face_refelem = Dune::ReferenceElements<DT,dim-1>::general(ig.geometry().type());
69 
70  // empty map means Dirichlet constraint
71  typename T::RowType empty;
72 
73  const FaceCoord testpoint = face_refelem.position(0,0);
74 
75  // Abort if this isn't a Dirichlet boundary
76  if (!param.isDirichlet(ig,testpoint))
77  return;
78 
79  for (std::size_t i=0;
80  i<std::size_t(FESwitch::coefficients(lfs.finiteElement()).size());
81  i++)
82  {
83  // The codim to which this dof is attached to
84  unsigned int codim =
85  FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
86 
87  if (codim==0) continue;
88 
89  for (int j=0; j<refelem.size(face,1,codim); j++){
90 
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;
95  }
96  }
97  }
98  };
99 
102  {
103  public:
104  enum { doProcessor = true };
105 
107 
112  template<typename IG, typename LFS, typename T>
113  void processor (const IG& ig, const LFS& lfs, T& trafo) const
114  {
115  typedef FiniteElementInterfaceSwitch<
116  typename LFS::Traits::FiniteElementType
117  > FESwitch;
118 
119  // determine face
120  const int face = ig.indexInInside();
121 
122  // find all local indices of this face
123  Dune::GeometryType gt = ig.inside().type();
124  typedef typename IG::ctype DT;
125  const int dim = IG::Entity::dimension;
126 
127  const Dune::ReferenceElement<DT,dim>& refelem = Dune::ReferenceElements<DT,dim>::general(gt);
128 
129  // empty map means Dirichlet constraint
130  typename T::RowType empty;
131 
132  // loop over all degrees of freedom and check if it is on given face
133  for (size_t i=0; i<FESwitch::coefficients(lfs.finiteElement()).size();
134  i++)
135  {
136  // The codim to which this dof is attached to
137  unsigned int codim =
138  FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
139 
140  if (codim==0) continue;
141 
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;
146  }
147  }
148  };
149 
151  template<typename GV>
153  {
154  public:
155  enum { doVolume = true };
156 
158 
163  template<typename P, typename EG, typename LFS, typename T>
164  void volume (const P& param, const EG& eg, const LFS& lfs, T& trafo) const
165  {
166  typedef FiniteElementInterfaceSwitch<
167  typename LFS::Traits::FiniteElementType
168  > FESwitch;
169 
170  auto entity = eg.entity();
171 
172  // nothing to do for interior entities
173  if (entity.partitionType()==Dune::InteriorEntity)
174  return;
175 
176  typedef typename FESwitch::Coefficients Coefficients;
177  const Coefficients& coeffs = FESwitch::coefficients(lfs.finiteElement());
178 
179  // empty map means Dirichlet constraint
180  typename T::RowType empty;
181 
182  const ReferenceElement<typename GV::ctype,GV::dimension>& ref_el =
183  ReferenceElements<typename GV::ctype,GV::dimension>::general(entity.type());
184 
185  // loop over all degrees of freedom and check if it is not owned by this processor
186  for (size_t i = 0; i < coeffs.size(); ++i)
187  {
188  size_t codim = coeffs.localKey(i).codim();
189  size_t sub_entity = coeffs.localKey(i).subEntity();
190 
191  size_t entity_index = _gv.indexSet().subIndex(entity,sub_entity,codim);
192  size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(sub_entity,codim));
193 
194  size_t index = _gt_offsets[gt_index] + entity_index;
195 
196  if (_ghosts[index])
197  {
198  trafo[lfs.dofIndex(i)] = empty;
199  }
200  }
201  }
202 
203  template<class GFS>
204  void compute_ghosts (const GFS& gfs)
205  {
206  std::fill(_gt_offsets.begin(),_gt_offsets.end(),0);
207 
208  for (size_t codim = 0; codim <= GV::dimension; ++codim)
209  {
210  if (gfs.ordering().contains(codim))
211  {
212  for (auto gt : _gv.indexSet().types(codim))
213  _gt_offsets[GlobalGeometryTypeIndex::index(gt) + 1] = _gv.indexSet().size(gt);
214  }
215  }
216 
217  std::partial_sum(_gt_offsets.begin(),_gt_offsets.end(),_gt_offsets.begin());
218 
219  _ghosts.assign(_gt_offsets.back(),true);
220 
221  typedef typename GV::template Codim<0>::
222  template Partition<Interior_Partition>::Iterator Iterator;
223 
224  for(Iterator it = _gv.template begin<0, Interior_Partition>(),
225  end = _gv.template end<0, Interior_Partition>();
226  it != end;
227  ++it)
228  {
229  auto entity = *it;
230 
231  const ReferenceElement<typename GV::ctype,GV::dimension>& ref_el =
232  ReferenceElements<typename GV::ctype,GV::dimension>::general(entity.type());
233 
234  for (size_t codim = 0; codim <= GV::dimension; ++codim)
235  if (gfs.ordering().contains(codim))
236  {
237  for (int i = 0; i < ref_el.size(codim); ++i)
238  {
239  size_t entity_index = _gv.indexSet().subIndex(entity,i,codim);
240  size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(i,codim));
241  size_t index = _gt_offsets[gt_index] + entity_index;
242 
243  _ghosts[index] = false;
244  }
245  }
246  }
247 
248  }
249 
251  : _gv(gv)
252  , _rank(gv.comm().rank())
253  , _gt_offsets(GlobalGeometryTypeIndex::size(GV::dimension) + 1)
254  {}
255 
256  private:
257 
258  GV _gv;
259  int _rank;
260  std::vector<bool> _ghosts;
261  std::vector<size_t> _gt_offsets;
262  };
264 
265  }
266 }
267 
268 #endif
extend conforming constraints class by processor boundary
Definition: conforming.hh:152
Dirichlet Constraints construction.
Definition: conforming.hh:36
const IG & ig
Definition: constraints.hh:147
void volume(const P &param, const EG &eg, const LFS &lfs, T &trafo) const
volume constraints
Definition: conforming.hh:164
void compute_ghosts(const GFS &gfs)
Definition: conforming.hh:204
Definition: adaptivity.hh:27
static const int dim
Definition: adaptivity.hh:83
void processor(const IG &ig, const LFS &lfs, T &trafo) const
processor constraints
Definition: conforming.hh:113
void boundary(const P &param, const IG &ig, const LFS &lfs, T &trafo) const
boundary constraints
Definition: conforming.hh:52
NonoverlappingConformingDirichletConstraints(const GV &gv)
Definition: conforming.hh:250
extend conforming constraints class by processor boundary
Definition: conforming.hh:101
const EG & eg
Definition: constraints.hh:280