DUNE PDELab (git)

conforming.hh
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
11
12#include <dune/geometry/referenceelements.hh>
13#include <dune/geometry/type.hh>
14
16
17#include <dune/localfunctions/common/interfaceswitch.hh>
18
19#include <dune/typetree/typetree.hh>
20
21#include <dune/pdelab/common/geometrywrapper.hh>
22#include <dune/pdelab/constraints/common/constraintsparameters.hh>
23#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
26
27namespace Dune {
28 namespace PDELab {
29
33
35 // works in any dimension and on all element types
37 : public TypeTree::LeafNode
38 {
39 public:
40 enum { doBoundary = true };
41 enum { doProcessor = false };
42 enum { doSkeleton = false };
43 enum { doVolume = false };
44
46
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
54 {
56 typename LFS::Traits::FiniteElementType
57 > FESwitch;
59
60 const int face = ig.indexInInside();
61
62 // find all local indices of this face
63 auto refelem = referenceElement(ig.inside().geometry());
64 auto face_refelem = referenceElement(ig.geometry());
65
66 // empty map means Dirichlet constraint
67 typename T::RowType empty;
68
69 const FaceCoord testpoint = face_refelem.position(0,0);
70
71 // Abort if this isn't a Dirichlet boundary
72 if (!param.isDirichlet(ig,testpoint))
73 return;
74
75 for (std::size_t i=0;
76 i<std::size_t(FESwitch::coefficients(lfs.finiteElement()).size());
77 i++)
78 {
79 // The codim to which this dof is attached to
80 unsigned int codim =
81 FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
82
83 if (codim==0) continue;
84
85 for (int j=0; j<refelem.size(face,1,codim); j++){
86
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;
91 }
92 }
93 }
94 };
95
98 {
99 public:
100 enum { doProcessor = true };
101
103
108 template<typename IG, typename LFS, typename T>
109 void processor (const IG& ig, const LFS& lfs, T& trafo) const
110 {
112 typename LFS::Traits::FiniteElementType
113 > FESwitch;
114
115 // determine face
116 const int face = ig.indexInInside();
117
118 auto refelem = referenceElement(ig.inside().geometry());
119
120 // empty map means Dirichlet constraint
121 typename T::RowType empty;
122
123 // loop over all degrees of freedom and check if it is on given face
124 for (size_t i=0; i<FESwitch::coefficients(lfs.finiteElement()).size();
125 i++)
126 {
127 // The codim to which this dof is attached to
128 unsigned int codim =
129 FESwitch::coefficients(lfs.finiteElement()).localKey(i).codim();
130
131 if (codim==0) continue;
132
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;
137 }
138 }
139 };
140
142 template<typename GV>
144 {
145 public:
146 enum { doVolume = true };
147
149
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
156 {
158 typename LFS::Traits::FiniteElementType
159 > FESwitch;
160
161 auto& entity = eg.entity();
162
163 // nothing to do for interior entities
164 if (entity.partitionType()==Dune::InteriorEntity)
165 return;
166
167 typedef typename FESwitch::Coefficients Coefficients;
168 const Coefficients& coeffs = FESwitch::coefficients(lfs.finiteElement());
169
170 // empty map means Dirichlet constraint
171 typename T::RowType empty;
172
173 auto ref_el = referenceElement(entity.geometry());
174
175 // loop over all degrees of freedom and check if it is not owned by this processor
176 for (size_t i = 0; i < coeffs.size(); ++i)
177 {
178 size_t codim = coeffs.localKey(i).codim();
179 size_t sub_entity = coeffs.localKey(i).subEntity();
180
181 size_t entity_index = _gv.indexSet().subIndex(entity,sub_entity,codim);
182 size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(sub_entity,codim));
183
184 size_t index = _gt_offsets[gt_index] + entity_index;
185
186 if (_ghosts[index])
187 {
188 trafo[lfs.dofIndex(i)] = empty;
189 }
190 }
191 }
192
193 template<class GFS>
194 void compute_ghosts (const GFS& gfs)
195 {
196 std::fill(_gt_offsets.begin(),_gt_offsets.end(),0);
197
198 for (size_t codim = 0; codim <= GV::dimension; ++codim)
199 {
200 if (gfs.ordering().contains(codim))
201 {
202 for (auto gt : _gv.indexSet().types(codim))
203 _gt_offsets[GlobalGeometryTypeIndex::index(gt) + 1] = _gv.indexSet().size(gt);
204 }
205 }
206
207 std::partial_sum(_gt_offsets.begin(),_gt_offsets.end(),_gt_offsets.begin());
208
209 _ghosts.assign(_gt_offsets.back(),true);
210
211 typedef typename GV::template Codim<0>::
212 template Partition<Interior_Partition>::Iterator Iterator;
213
214 for(Iterator it = _gv.template begin<0, Interior_Partition>(),
215 end = _gv.template end<0, Interior_Partition>();
216 it != end;
217 ++it)
218 {
219 auto entity = *it;
220
221 auto ref_el = referenceElement(entity.geometry());
222
223 for (size_t codim = 0; codim <= GV::dimension; ++codim)
224 if (gfs.ordering().contains(codim))
225 {
226 for (int i = 0; i < ref_el.size(codim); ++i)
227 {
228 size_t entity_index = _gv.indexSet().subIndex(entity,i,codim);
229 size_t gt_index = GlobalGeometryTypeIndex::index(ref_el.type(i,codim));
230 size_t index = _gt_offsets[gt_index] + entity_index;
231
232 _ghosts[index] = false;
233 }
234 }
235 }
236
237 }
238
239 NonoverlappingConformingDirichletConstraints(const GV& gv)
240 : _gv(gv)
241 , _rank(gv.comm().rank())
242 , _gt_offsets(GlobalGeometryTypeIndex::size(GV::dimension) + 1)
243 {}
244
245 private:
246
247 GV _gv;
248 int _rank;
249 std::vector<bool> _ghosts;
250 std::vector<size_t> _gt_offsets;
251 };
253
254 }
255}
256
257#endif // DUNE_PDELAB_CONSTRAINTS_CONFORMING_HH
vector space out of a tensor product of fields.
Definition: fvector.hh:91
static constexpr std::size_t index(const GeometryType &gt)
Compute the index for the given geometry type over all dimensions.
Definition: typeindex.hh:138
Dirichlet Constraints construction.
Definition: conforming.hh:38
void boundary(const P &param, const IG &ig, const LFS &lfs, T &trafo) const
boundary constraints
Definition: conforming.hh:53
extend conforming constraints class by processor boundary
Definition: conforming.hh:144
void volume(const P &param, const EG &eg, const LFS &lfs, T &trafo) const
volume constraints
Definition: conforming.hh:155
extend conforming constraints class by processor boundary
Definition: conforming.hh:98
void processor(const IG &ig, const LFS &lfs, T &trafo) const
processor constraints
Definition: conforming.hh:109
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
@ InteriorEntity
all interior entities
Definition: gridenums.hh:31
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.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)