dune-pdelab  2.5-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_CONSTRAINTS_CONFORMING_HH
5 #define DUNE_PDELAB_CONSTRAINTS_CONFORMING_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::mydimension> FaceCoord;
58 
59  const int face = ig.indexInInside();
60 
61  // find all local indices of this face
62  auto refelem = referenceElement(ig.inside().geometry());
63  auto face_refelem = referenceElement(ig.geometry());
64 
65  // empty map means Dirichlet constraint
66  typename T::RowType empty;
67 
68  const FaceCoord testpoint = face_refelem.position(0,0);
69 
70  // Abort if this isn't a Dirichlet boundary
71  if (!param.isDirichlet(ig,testpoint))
72  return;
73 
74  for (std::size_t i=0;
75  i<std::size_t(FESwitch::coefficients(lfs.finiteElement()).size());
76  i++)
77  {
78  // The codim to which this dof is attached to
79  unsigned int codim =
80  FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
81 
82  if (codim==0) continue;
83 
84  for (int j=0; j<refelem.size(face,1,codim); j++){
85 
86  if (static_cast<int>(FESwitch::coefficients(lfs.finiteElement()).
87  localKey(i).subEntity())
88  == refelem.subEntity(face,1,j,codim))
89  trafo[lfs.dofIndex(i)] = empty;
90  }
91  }
92  }
93  };
94 
97  {
98  public:
99  enum { doProcessor = true };
100 
102 
107  template<typename IG, typename LFS, typename T>
108  void processor (const IG& ig, const LFS& lfs, T& trafo) const
109  {
110  typedef FiniteElementInterfaceSwitch<
111  typename LFS::Traits::FiniteElementType
112  > FESwitch;
113 
114  // determine face
115  const int face = ig.indexInInside();
116 
117  auto refelem = referenceElement(ig.inside().geometry());
118 
119  // empty map means Dirichlet constraint
120  typename T::RowType empty;
121 
122  // loop over all degrees of freedom and check if it is on given face
123  for (size_t i=0; i<FESwitch::coefficients(lfs.finiteElement()).size();
124  i++)
125  {
126  // The codim to which this dof is attached to
127  unsigned int codim =
128  FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
129 
130  if (codim==0) continue;
131 
132  for (int j=0; j<refelem.size(face,1,codim); j++)
133  if (FESwitch::coefficients(lfs.finiteElement()).localKey(i).
134  subEntity() == std::size_t(refelem.subEntity(face,1,j,codim)))
135  trafo[lfs.dofIndex(i)] = empty;
136  }
137  }
138  };
139 
141  template<typename GV>
143  {
144  public:
145  enum { doVolume = true };
146 
148 
153  template<typename P, typename EG, typename LFS, typename T>
154  void volume (const P& param, const EG& eg, const LFS& lfs, T& trafo) const
155  {
156  typedef FiniteElementInterfaceSwitch<
157  typename LFS::Traits::FiniteElementType
158  > FESwitch;
159 
160  auto& entity = eg.entity();
161 
162  // nothing to do for interior entities
163  if (entity.partitionType()==Dune::InteriorEntity)
164  return;
165 
166  typedef typename FESwitch::Coefficients Coefficients;
167  const Coefficients& coeffs = FESwitch::coefficients(lfs.finiteElement());
168 
169  // empty map means Dirichlet constraint
170  typename T::RowType empty;
171 
172  auto ref_el = referenceElement(entity.geometry());
173 
174  // loop over all degrees of freedom and check if it is not owned by this processor
175  for (size_t i = 0; i < coeffs.size(); ++i)
176  {
177  size_t codim = coeffs.localKey(i).codim();
178  size_t sub_entity = coeffs.localKey(i).subEntity();
179 
180  size_t entity_index = _gv.indexSet().subIndex(entity,sub_entity,codim);
181  size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(sub_entity,codim));
182 
183  size_t index = _gt_offsets[gt_index] + entity_index;
184 
185  if (_ghosts[index])
186  {
187  trafo[lfs.dofIndex(i)] = empty;
188  }
189  }
190  }
191 
192  template<class GFS>
193  void compute_ghosts (const GFS& gfs)
194  {
195  std::fill(_gt_offsets.begin(),_gt_offsets.end(),0);
196 
197  for (size_t codim = 0; codim <= GV::dimension; ++codim)
198  {
199  if (gfs.ordering().contains(codim))
200  {
201  for (auto gt : _gv.indexSet().types(codim))
202  _gt_offsets[GlobalGeometryTypeIndex::index(gt) + 1] = _gv.indexSet().size(gt);
203  }
204  }
205 
206  std::partial_sum(_gt_offsets.begin(),_gt_offsets.end(),_gt_offsets.begin());
207 
208  _ghosts.assign(_gt_offsets.back(),true);
209 
210  typedef typename GV::template Codim<0>::
211  template Partition<Interior_Partition>::Iterator Iterator;
212 
213  for(Iterator it = _gv.template begin<0, Interior_Partition>(),
214  end = _gv.template end<0, Interior_Partition>();
215  it != end;
216  ++it)
217  {
218  auto entity = *it;
219 
220  auto ref_el = referenceElement(entity.geometry());
221 
222  for (size_t codim = 0; codim <= GV::dimension; ++codim)
223  if (gfs.ordering().contains(codim))
224  {
225  for (int i = 0; i < ref_el.size(codim); ++i)
226  {
227  size_t entity_index = _gv.indexSet().subIndex(entity,i,codim);
228  size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(i,codim));
229  size_t index = _gt_offsets[gt_index] + entity_index;
230 
231  _ghosts[index] = false;
232  }
233  }
234  }
235 
236  }
237 
239  : _gv(gv)
240  , _rank(gv.comm().rank())
241  , _gt_offsets(GlobalGeometryTypeIndex::size(GV::dimension) + 1)
242  {}
243 
244  private:
245 
246  GV _gv;
247  int _rank;
248  std::vector<bool> _ghosts;
249  std::vector<size_t> _gt_offsets;
250  };
252 
253  }
254 }
255 
256 #endif // DUNE_PDELAB_CONSTRAINTS_CONFORMING_HH
Dune::PDELab::OverlappingConformingDirichletConstraints::doProcessor
@ doProcessor
Definition: conforming.hh:99
geometrywrapper.hh
gridfunctionspace.hh
index
std::size_t index
Definition: interpolate.hh:118
Dune::PDELab::OverlappingConformingDirichletConstraints::processor
void processor(const IG &ig, const LFS &lfs, T &trafo) const
processor constraints
Definition: conforming.hh:108
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::NonoverlappingConformingDirichletConstraints
extend conforming constraints class by processor boundary
Definition: conforming.hh:142
Dune::PDELab::ConformingDirichletConstraints::doSkeleton
@ doSkeleton
Definition: conforming.hh:41
localfunctionspacetags.hh
constraintsparameters.hh
Dune::PDELab::OverlappingConformingDirichletConstraints
extend conforming constraints class by processor boundary
Definition: conforming.hh:96
Dune::PDELab::ConformingDirichletConstraints::doVolume
@ doVolume
Definition: conforming.hh:42
Dune::PDELab::NonoverlappingConformingDirichletConstraints::NonoverlappingConformingDirichletConstraints
NonoverlappingConformingDirichletConstraints(const GV &gv)
Definition: conforming.hh:238
Dune::PDELab::ConformingDirichletConstraints::doProcessor
@ doProcessor
Definition: conforming.hh:40
localvector.hh
ig
const IG & ig
Definition: constraints.hh:148
Dune::PDELab::NonoverlappingConformingDirichletConstraints::doVolume
@ doVolume
Definition: conforming.hh:145
Dune::PDELab::ConformingDirichletConstraints::boundary
void boundary(const P &param, const IG &ig, const LFS &lfs, T &trafo) const
boundary constraints
Definition: conforming.hh:52
Dune::PDELab::NonoverlappingConformingDirichletConstraints::compute_ghosts
void compute_ghosts(const GFS &gfs)
Definition: conforming.hh:193
Dune::PDELab::ConformingDirichletConstraints::doBoundary
@ doBoundary
Definition: conforming.hh:39
Dune::PDELab::NonoverlappingConformingDirichletConstraints::volume
void volume(const P &param, const EG &eg, const LFS &lfs, T &trafo) const
volume constraints
Definition: conforming.hh:154
Dune::PDELab::ConformingDirichletConstraints
Dirichlet Constraints construction.
Definition: conforming.hh:36