dune-pdelab  2.7-git
constraints.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_COMMON_CONSTRAINTS_HH
5 #define DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/shared_ptr.hh>
9 #include<dune/common/float_cmp.hh>
10 
18 
19 namespace Dune {
20  namespace PDELab {
21 
25 
26  namespace { // hide internals
27  // do method invocation only if class has the method
28 
29  template<typename C, bool doIt>
30  struct ConstraintsCallBoundary
31  {
32  template<typename P, typename IG, typename LFS, typename T>
33  static void boundary (const C& c, const P& p, const IG& ig, const LFS& lfs, T& trafo)
34  {
35  }
36  };
37  template<typename C, bool doIt>
38  struct ConstraintsCallProcessor
39  {
40  template<typename IG, typename LFS, typename T>
41  static void processor (const C& c, const IG& ig, const LFS& lfs, T& trafo)
42  {
43  }
44  };
45  template<typename C, bool doIt>
46  struct ConstraintsCallSkeleton
47  {
48  template<typename IG, typename LFS, typename T>
49  static void skeleton (const C& c, const IG& ig,
50  const LFS& lfs_e, const LFS& lfs_f,
51  T& trafo_e, T& trafo_f)
52  {
53  }
54  };
55  template<typename C, bool doIt>
56  struct ConstraintsCallVolume
57  {
58  template<typename P, typename EG, typename LFS, typename T>
59  static void volume (const C& c, const P&, const EG& eg, const LFS& lfs, T& trafo)
60  {
61  }
62  };
63 
64 
65  template<typename C>
66  struct ConstraintsCallBoundary<C,true>
67  {
68  template<typename P, typename IG, typename LFS, typename T>
69  static void boundary (const C& c, const P& p, const IG& ig, const LFS& lfs, T& trafo)
70  {
71  if (lfs.size())
72  c.boundary(p,ig,lfs,trafo);
73  }
74  };
75  template<typename C>
76  struct ConstraintsCallProcessor<C,true>
77  {
78  template<typename IG, typename LFS, typename T>
79  static void processor (const C& c, const IG& ig, const LFS& lfs, T& trafo)
80  {
81  if (lfs.size())
82  c.processor(ig,lfs,trafo);
83  }
84  };
85  template<typename C>
86  struct ConstraintsCallSkeleton<C,true>
87  {
88  template<typename IG, typename LFS, typename T>
89  static void skeleton (const C& c, const IG& ig,
90  const LFS& lfs_e, const LFS& lfs_f,
91  T& trafo_e, T& trafo_f)
92  {
93  if (lfs_e.size() || lfs_f.size())
94  c.skeleton(ig, lfs_e, lfs_f, trafo_e, trafo_f);
95  }
96  };
97  template<typename C>
98  struct ConstraintsCallVolume<C,true>
99  {
100  template<typename P, typename EG, typename LFS, typename T>
101  static void volume (const C& c, const P& p, const EG& eg, const LFS& lfs, T& trafo)
102  {
103  if (lfs.size())
104  c.volume(p,eg,lfs,trafo);
105  }
106  };
107 
108 
109  struct ParameterizedConstraintsBase
110  : public TypeTree::TreePairVisitor
111  {
112  // This acts as a catch-all for unsupported leaf- / non-leaf combinations in the two
113  // trees. It is necessary because otherwise, the visitor would fall back to the default
114  // implementation in TreeVisitor, which simply does nothing. The resulting bugs would
115  // probably be hell to find...
116  template<typename P, typename LFS, typename TreePath>
117  void leaf(const P& p, const LFS& lfs, TreePath treePath) const
118  {
119  static_assert((AlwaysFalse<P>::Value),
120  "unsupported combination of function and LocalFunctionSpace");
121  }
122  };
123 
124 
125  template<typename P, typename IG, typename CL>
126  struct BoundaryConstraintsForParametersLeaf
127  : public TypeTree::TreeVisitor
128  , public TypeTree::DynamicTraversal
129  {
130 
131  template<typename LFS, typename TreePath>
132  void leaf(const LFS& lfs, TreePath treePath) const
133  {
134 
135  // extract constraints type
136  typedef typename LFS::Traits::ConstraintsType C;
137 
138  // iterate over boundary, need intersection iterator
139  ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),p,ig,lfs,cl);
140  }
141 
142  BoundaryConstraintsForParametersLeaf(const P& p_, const IG& ig_, CL& cl_)
143  : p(p_)
144  , ig(ig_)
145  , cl(cl_)
146  {}
147 
148  const P& p;
149  const IG& ig;
150  CL& cl;
151 
152  };
153 
154 
155  template<typename IG, typename CL>
156  struct BoundaryConstraints
157  : public ParameterizedConstraintsBase
158  , public TypeTree::DynamicTraversal
159  {
160 
161  // standard case - leaf in both trees
162  template<typename P, typename LFS, typename TreePath>
163  typename std::enable_if<P::isLeaf && LFS::isLeaf>::type
164  leaf(const P& p, const LFS& lfs, TreePath treePath) const
165  {
166  // extract constraints type
167  typedef typename LFS::Traits::ConstraintsType C;
168 
169  // iterate over boundary, need intersection iterator
170  ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),p,ig,lfs,cl);
171  }
172 
173  // reuse constraints parameter information from p for all LFS children
174  template<typename P, typename LFS, typename TreePath>
175  typename std::enable_if<P::isLeaf && (!LFS::isLeaf)>::type
176  leaf(const P& p, const LFS& lfs, TreePath treePath) const
177  {
178  // traverse LFS tree and reuse parameter information
179  TypeTree::applyToTree(lfs,BoundaryConstraintsForParametersLeaf<P,IG,CL>(p,ig,cl));
180  }
181 
182  BoundaryConstraints(const IG& ig_, CL& cl_)
183  : ig(ig_)
184  , cl(cl_)
185  {}
186 
187  private:
188  const IG& ig;
189  CL& cl;
190 
191  };
192 
193 
194  template<typename IG, typename CL>
195  struct ProcessorConstraints
196  : public TypeTree::TreeVisitor
197  , public TypeTree::DynamicTraversal
198  {
199 
200  template<typename LFS, typename TreePath>
201  void leaf(const LFS& lfs, TreePath treePath) const
202  {
203  // extract constraints type
204  typedef typename LFS::Traits::ConstraintsType C;
205 
206  // iterate over boundary, need intersection iterator
207  ConstraintsCallProcessor<C,C::doProcessor>::processor(lfs.constraints(),ig,lfs,cl);
208  }
209 
210  ProcessorConstraints(const IG& ig_, CL& cl_)
211  : ig(ig_)
212  , cl(cl_)
213  {}
214 
215  private:
216  const IG& ig;
217  CL& cl;
218 
219  };
220 
221 
222  template<typename IG, typename CL>
223  struct SkeletonConstraints
224  : public TypeTree::TreePairVisitor
225  , public TypeTree::DynamicTraversal
226  {
227 
228  template<typename LFS, typename TreePath>
229  void leaf(const LFS& lfs_e, const LFS& lfs_f, TreePath treePath) const
230  {
231  // extract constraints type
232  typedef typename LFS::Traits::ConstraintsType C;
233 
234  // as LFS::constraints() just returns the constraints of the
235  // GridFunctionSpace, lfs_e.constraints() is equivalent to
236  // lfs_f.constraints()
237  const C & c = lfs_e.constraints();
238 
239  // iterate over boundary, need intersection iterator
240  ConstraintsCallSkeleton<C,C::doSkeleton>::skeleton(c,ig,lfs_e,lfs_f,cl_e,cl_f);
241  }
242 
243  SkeletonConstraints(const IG& ig_, CL& cl_e_, CL& cl_f_)
244  : ig(ig_)
245  , cl_e(cl_e_)
246  , cl_f(cl_f_)
247  {}
248 
249  private:
250  const IG& ig;
251  CL& cl_e;
252  CL& cl_f;
253 
254  };
255 
256 
257  template<typename P, typename EG, typename CL>
258  struct VolumeConstraintsForParametersLeaf
259  : public TypeTree::TreeVisitor
260  , public TypeTree::DynamicTraversal
261  {
262 
263  template<typename LFS, typename TreePath>
264  void leaf(const LFS& lfs, TreePath treePath) const
265  {
266  // extract constraints type
267  typedef typename LFS::Traits::ConstraintsType C;
268  const C & c = lfs.constraints();
269 
270  // iterate over boundary, need intersection iterator
271  ConstraintsCallVolume<C,C::doVolume>::volume(c,p,eg,lfs,cl);
272  }
273 
274  VolumeConstraintsForParametersLeaf(const P& p_, const EG& eg_, CL& cl_)
275  : p(p_)
276  , eg(eg_)
277  , cl(cl_)
278  {}
279 
280  private:
281  const P& p;
282  const EG& eg;
283  CL& cl;
284 
285  };
286 
287 
288  template<typename EG, typename CL>
289  struct VolumeConstraints
290  : public ParameterizedConstraintsBase
291  , public TypeTree::DynamicTraversal
292  {
293 
294  // standard case - leaf in both trees
295  template<typename P, typename LFS, typename TreePath>
296  typename std::enable_if<P::isLeaf && LFS::isLeaf>::type
297  leaf(const P& p, const LFS& lfs, TreePath treePath) const
298  {
299  // extract constraints type
300  typedef typename LFS::Traits::ConstraintsType C;
301  const C & c = lfs.constraints();
302 
303  // iterate over boundary, need intersection iterator
304  ConstraintsCallVolume<C,C::doVolume>::volume(c,p,eg,lfs,cl);
305  }
306 
307  // reuse constraints parameter information from p for all LFS children
308  template<typename P, typename LFS, typename TreePath>
309  typename std::enable_if<P::isLeaf && (!LFS::isLeaf)>::type
310  leaf(const P& p, const LFS& lfs, TreePath treePath) const
311  {
312  // traverse LFS tree and reuse parameter information
313  TypeTree::applyToTree(lfs,VolumeConstraintsForParametersLeaf<P,EG,CL>(p,eg,cl));
314  }
315 
316  VolumeConstraints(const EG& eg_, CL& cl_)
317  : eg(eg_)
318  , cl(cl_)
319  {}
320 
321  private:
322  const EG& eg;
323  CL& cl;
324 
325  };
326 
327 
328  } // anonymous namespace
329 
330  template<typename... Children>
332  : public TypeTree::CompositeNode<Children...>
333  {
334  typedef TypeTree::CompositeNode<Children...> BaseT;
335 
336  CompositeConstraintsOperator(Children&... children)
337  : BaseT(children...)
338  {}
339 
340  CompositeConstraintsOperator(std::shared_ptr<Children>... children)
341  : BaseT(children...)
342  {}
343 
344  // aggregate flags
345 
346  // forward methods to childs
347 
348  };
349 
350  template<typename... Children>
352  : public TypeTree::CompositeNode<Children...>
353  {
354  typedef TypeTree::CompositeNode<Children...> BaseT;
355 
356  CompositeConstraintsParameters(Children&... children)
357  : BaseT(children...)
358  {}
359 
360  CompositeConstraintsParameters(std::shared_ptr<Children>... children)
361  : BaseT(children...)
362  {}
363  };
364 
365  template<typename T, std::size_t k>
367  : public TypeTree::PowerNode<T,k>
368  {
369  typedef TypeTree::PowerNode<T,k> BaseT;
370 
372  : BaseT()
373  {}
374 
376  : BaseT(c)
377  {}
378 
380  T& c1)
381  : BaseT(c0,c1)
382  {}
383 
385  T& c1,
386  T& c2)
387  : BaseT(c0,c1,c2)
388  {}
389 
391  T& c1,
392  T& c2,
393  T& c3)
394  : BaseT(c0,c1,c2,c3)
395  {}
396 
398  T& c1,
399  T& c2,
400  T& c3,
401  T& c4)
402  : BaseT(c0,c1,c2,c3,c4)
403  {}
404 
406  T& c1,
407  T& c2,
408  T& c3,
409  T& c4,
410  T& c5)
411  : BaseT(c0,c1,c2,c3,c4,c5)
412  {}
413 
415  T& c1,
416  T& c2,
417  T& c3,
418  T& c4,
419  T& c5,
420  T& c6)
421  : BaseT(c0,c1,c2,c3,c4,c5,c6)
422  {}
423 
425  T& c1,
426  T& c2,
427  T& c3,
428  T& c4,
429  T& c5,
430  T& c6,
431  T& c7)
432  : BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
433  {}
434 
436  T& c1,
437  T& c2,
438  T& c3,
439  T& c4,
440  T& c5,
441  T& c6,
442  T& c7,
443  T& c8)
444  : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
445  {}
446 
448  T& c1,
449  T& c2,
450  T& c3,
451  T& c4,
452  T& c5,
453  T& c6,
454  T& c7,
455  T& c8,
456  T& c9)
457  : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
458  {}
459 
460  PowerConstraintsParameters (const std::array<std::shared_ptr<T>,k>& children)
461  : BaseT(children)
462  {}
463  };
464 
465 #ifndef DOXYGEN
466 
468  template<typename F>
469  class OldStyleConstraintsWrapper
470  : public TypeTree::LeafNode
471  {
472  std::shared_ptr<const F> _f;
473  unsigned int _i;
474  public:
475 
476  template<typename Transformation>
477  OldStyleConstraintsWrapper(std::shared_ptr<const F> f, const Transformation& t, unsigned int i=0)
478  : _f(f)
479  , _i(i)
480  {}
481 
482  template<typename Transformation>
483  OldStyleConstraintsWrapper(const F & f, const Transformation& t, unsigned int i=0)
484  : _f(stackobject_to_shared_ptr(f))
485  , _i(i)
486  {}
487 
488  template<typename I>
489  bool isDirichlet(const I & intersection, const FieldVector<typename I::ctype, I::mydimension> & coord) const
490  {
491  typename F::Traits::RangeType bctype;
492  _f->evaluate(intersection,coord,bctype);
493  return bctype[_i] > 0;
494  }
495 
496  template<typename I>
497  bool isNeumann(const I & intersection, const FieldVector<typename I::ctype, I::mydimension> & coord) const
498  {
499  typename F::Traits::RangeType bctype;
500  _f->evaluate(intersection,coord,bctype);
501  return bctype[_i] == 0;
502  }
503  };
504 
505  // Tag to name trafo GridFunction -> OldStyleConstraintsWrapper
506  struct gf_to_constraints {};
507 
508  // register trafos GridFunction -> OldStyleConstraintsWrapper
509  template<typename F, typename Transformation>
510  struct MultiComponentOldStyleConstraintsWrapperDescription
511  {
512 
513  static const bool recursive = false;
514 
515  enum { dim = F::Traits::dimRange };
516  typedef OldStyleConstraintsWrapper<F> node_type;
517  typedef PowerConstraintsParameters<node_type, dim> transformed_type;
518  typedef std::shared_ptr<transformed_type> transformed_storage_type;
519 
520  static transformed_type transform(const F& s, const Transformation& t)
521  {
522  std::shared_ptr<const F> sp = stackobject_to_shared_ptr(s);
523  std::array<std::shared_ptr<node_type>, dim> childs;
524  for (int i=0; i<dim; i++)
525  childs[i] = std::make_shared<node_type>(sp,t,i);
526  return transformed_type(childs);
527  }
528 
529  static transformed_storage_type transform_storage(std::shared_ptr<const F> s, const Transformation& t)
530  {
531  std::array<std::shared_ptr<node_type>, dim> childs;
532  for (int i=0; i<dim; i++)
533  childs[i] = std::make_shared<node_type>(s,t,i);
534  return std::make_shared<transformed_type>(childs);
535  }
536 
537  };
538  // trafos for leaf nodes
539  template<typename GridFunction>
540  typename std::conditional<
541  (GridFunction::Traits::dimRange == 1),
542  // trafo for scalar leaf nodes
543  TypeTree::GenericLeafNodeTransformation<GridFunction,gf_to_constraints,OldStyleConstraintsWrapper<GridFunction> >,
544  // trafo for multi component leaf nodes
545  MultiComponentOldStyleConstraintsWrapperDescription<GridFunction,gf_to_constraints>
546  >::type
547  registerNodeTransformation(GridFunction*, gf_to_constraints*, GridFunctionTag*);
548 
549  // trafo for power nodes
550  template<typename PowerGridFunction>
551  TypeTree::SimplePowerNodeTransformation<PowerGridFunction,gf_to_constraints,PowerConstraintsParameters>
552  registerNodeTransformation(PowerGridFunction*, gf_to_constraints*, PowerGridFunctionTag*);
553 
554  // trafos for composite nodes
555  template<typename CompositeGridFunction>
556  TypeTree::SimpleCompositeNodeTransformation<CompositeGridFunction,gf_to_constraints,CompositeConstraintsParameters>
557  registerNodeTransformation(CompositeGridFunction*, gf_to_constraints*, CompositeGridFunctionTag*);
558 
560 
570  template<typename P, typename GFS, typename CG, bool isFunction>
571  struct ConstraintsAssemblerHelper
572  {
574 
588  static void
589  assemble(const P& p, const GFS& gfs, CG& cg, const bool verbose)
590  {
591  // get some types
592  using ES = typename GFS::Traits::EntitySet;
593  using Element = typename ES::Traits::Element;
594  using Intersection = typename ES::Traits::Intersection;
595 
596  ES es = gfs.entitySet();
597 
598  // make local function space
599  using LFS = LocalFunctionSpace<GFS>;
600  LFS lfs_e(gfs);
601  LFSIndexCache<LFS> lfs_cache_e(lfs_e);
602  LFS lfs_f(gfs);
603  LFSIndexCache<LFS> lfs_cache_f(lfs_f);
604 
605  // get index set
606  auto& is = es.indexSet();
607 
608  // loop once over the grid
609  for (const auto& element : elements(es))
610  {
611 
612  auto id = is.uniqueIndex(element);
613 
614  // bind local function space to element
615  lfs_e.bind(element);
616 
617  using CL = typename CG::LocalTransformation;
618 
619  CL cl_self;
620 
621  using ElementWrapper = ElementGeometry<Element>;
622  using IntersectionWrapper = IntersectionGeometry<Intersection>;
623 
624  TypeTree::applyToTreePair(p,lfs_e,VolumeConstraints<ElementWrapper,CL>(ElementWrapper(element),cl_self));
625 
626  // iterate over intersections and call metaprogram
627  unsigned int intersection_index = 0;
628  for (const auto& intersection : intersections(es,element))
629  {
630 
631  auto intersection_data = classifyIntersection(es,intersection);
632  auto intersection_type = std::get<0>(intersection_data);
633  auto& outside_element = std::get<1>(intersection_data);
634 
635  switch (intersection_type) {
636 
639  {
640  auto idn = is.uniqueIndex(outside_element);
641 
642  if(id > idn){
643  // bind local function space to element in neighbor
644  lfs_f.bind(outside_element);
645 
646  CL cl_neighbor;
647 
648  TypeTree::applyToTreePair(lfs_e,lfs_f,SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
649 
650  if (!cl_neighbor.empty())
651  {
652  lfs_cache_f.update();
653  cg.import_local_transformation(cl_neighbor,lfs_cache_f);
654  }
655 
656  }
657  break;
658  }
659 
661  TypeTree::applyToTreePair(p,lfs_e,BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
662  break;
663 
665  TypeTree::applyToTree(lfs_e,ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
666  break;
667 
668  }
669  ++intersection_index;
670  }
671 
672  if (!cl_self.empty())
673  {
674  lfs_cache_e.update();
675  cg.import_local_transformation(cl_self,lfs_cache_e);
676  }
677 
678  }
679 
680  // print result
681  if(verbose){
682  std::cout << "constraints:" << std::endl;
683 
684  std::cout << cg.size() << " constrained degrees of freedom" << std::endl;
685 
686  for (const auto& col : cg)
687  {
688  std::cout << col.first << ": "; // col index
689  for (const auto& row : col.second)
690  std::cout << "(" << row.first << "," << row.second << ") "; // row index , value
691  std::cout << std::endl;
692  }
693  }
694  }
695  }; // end ConstraintsAssemblerHelper
696 
697 
698 
699  // Disable constraints assembly for empty transformation
700  template<typename F, typename GFS>
701  struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, true>
702  {
703  static void assemble(const F& f, const GFS& gfs, EmptyTransformation& cg, const bool verbose)
704  {}
705  };
706 
707  // Disable constraints assembly for empty transformation
708  template<typename F, typename GFS>
709  struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, false>
710  {
711  static void assemble(const F& f, const GFS& gfs, EmptyTransformation& cg, const bool verbose)
712  {}
713  };
714 
715 
716 
717  // Backwards compatibility shim
718  template<typename F, typename GFS, typename CG>
719  struct ConstraintsAssemblerHelper<F, GFS, CG, true>
720  {
721  static void
722  assemble(const F& f, const GFS& gfs, CG& cg, const bool verbose)
723  {
724  // type of transformed tree
725  typedef typename TypeTree::TransformTree<F,gf_to_constraints> Transformation;
726  typedef typename Transformation::Type P;
727  // transform tree
728  P p = Transformation::transform(f);
729  // call parameter based implementation
731  }
732  };
733 #endif
734 
736 
748  template<typename GFS, typename CG>
749  void constraints(const GFS& gfs, CG& cg,
750  const bool verbose = false)
751  {
753  ConstraintsAssemblerHelper<NoConstraintsParameters, GFS, CG, false>::assemble(p,gfs,cg,verbose);
754  }
755 
757 
774  template<typename P, typename GFS, typename CG>
775  void constraints(const P& p, const GFS& gfs, CG& cg,
776  const bool verbose = false)
777  {
778  // clear global constraints
779  cg.clear();
781  }
782 
784 
795  template<typename CG, typename XG>
796  void set_constrained_dofs(const CG& cg,
797  typename XG::ElementType x,
798  XG& xg)
799  {
800  for (const auto& col : cg)
801  xg[col.first] = x;
802  }
803 
804 
805 #ifndef DOXYGEN
806 
807  // Specialized version for unconstrained spaces
808  template<typename XG>
809  void set_constrained_dofs(const EmptyTransformation& cg,
810  typename XG::ElementType x,
811  XG& xg)
812  {}
813 
814 #endif // DOXYGEN
815 
816 
818 
838  template<typename CG, typename XG, typename Cmp>
839  bool check_constrained_dofs(const CG& cg, typename XG::ElementType x,
840  XG& xg, const Cmp& cmp = Cmp())
841  {
842  for (const auto& col : cg)
843  if(cmp.ne(xg[col.first], x))
844  return false;
845  return true;
846  }
847 
849 
867  template<typename CG, typename XG>
868  bool check_constrained_dofs(const CG& cg, typename XG::ElementType x,
869  XG& xg)
870  {
871  return check_constrained_dofs(cg, x, xg,
872  FloatCmpOps<typename XG::ElementType>());
873  }
874 
875 
876 #ifndef DOXYGEN
877 
878  // Specialized version for unconstrained spaces
879  template<typename XG, typename Cmp>
880  bool check_constrained_dofs(const EmptyTransformation& cg, typename XG::ElementType x,
881  XG& xg, const Cmp& cmp = Cmp())
882  {
883  return true;
884  }
885 
886  // Specialized version for unconstrained spaces
887  template<typename XG>
888  bool check_constrained_dofs(const EmptyTransformation& cg, typename XG::ElementType x,
889  XG& xg)
890  {
891  return true;
892  }
893 
894 #endif // DOXYGEN
895 
896 
898 
903  template<typename CG, typename XG>
904  void constrain_residual (const CG& cg, XG& xg)
905  {
906  for (const auto& col : cg)
907  for (const auto& row : col.second)
908  xg[row.first] += row.second * xg[col.first];
909 
910  // extra loop because constrained dofs might have contributions
911  // to constrained dofs
912  for (const auto& col : cg)
913  xg[col.first] = typename XG::ElementType(0);
914  }
915 
916 
917 #ifndef DOXYGEN
918 
919  // Specialized version for unconstrained spaces
920  template<typename XG>
921  void constrain_residual (const EmptyTransformation& cg, XG& xg)
922  {}
923 
924 #endif // DOXYGEN
925 
929 
935  template<typename CG, typename XG>
936  void copy_constrained_dofs (const CG& cg, const XG& xgin, XG& xgout)
937  {
938  for (const auto& col : cg)
939  xgout[col.first] = xgin[col.first];
940  }
941 
942 
943 #ifndef DOXYGEN
944 
945  // Specialized version for unconstrained spaces
946  template<typename XG>
947  void copy_constrained_dofs (const EmptyTransformation& cg, const XG& xgin, XG& xgout)
948  {}
949 
950 #endif // DOXYGEN
951 
952 
959  template<typename CG, typename XG>
960  void set_nonconstrained_dofs (const CG& cg, typename XG::ElementType x, XG& xg)
961  {
962  XG tmp(xg);
963  xg = x;
964  copy_constrained_dofs(cg,tmp,xg);
965  }
966 
967 
968 #ifndef DOXYGEN
969 
970  // Specialized version for unconstrained spaces
971  template<typename XG>
972  void set_nonconstrained_dofs (const EmptyTransformation& cg, typename XG::ElementType x, XG& xg)
973  {
974  xg = x;
975  }
976 
977 #endif // DOXYGEN
978 
979 
986  template<typename CG, typename XG>
987  void copy_nonconstrained_dofs (const CG& cg, const XG& xgin, XG& xgout)
988  {
989  XG tmp(xgin);
990  copy_constrained_dofs(cg,xgout,tmp);
991  xgout = tmp;
992  }
993 
994 
995 #ifndef DOXYGEN
996 
997  // Specialized version for unconstrained spaces
998  template<typename XG>
999  void copy_nonconstrained_dofs (const EmptyTransformation& cg, const XG& xgin, XG& xgout)
1000  {
1001  xgout = xgin;
1002  }
1003 
1004 #endif // DOXYGEN
1005 
1006 
1013  template<typename CG, typename XG>
1014  void set_shifted_dofs (const CG& cg, typename XG::ElementType x, XG& xg)
1015  {
1016  XG tmp(xg);
1017  tmp = x;
1018 
1019  for (const auto& col : cg)
1020  {
1021  // this is our marker for Dirichlet constraints
1022  if (col.second.size() == 0)
1023  {
1024  tmp[col.first] = xg[col.first];
1025  }
1026  }
1027 
1028  xg = tmp;
1029  }
1030 
1031 
1032 #ifndef DOXYGEN
1033 
1034  // Specialized version for unconstrained spaces
1035  template<typename XG>
1036  void set_shifted_dofs (const EmptyTransformation& cg, typename XG::ElementType x, XG& xg)
1037  {}
1038 
1039 #endif // DOXYGEN
1040 
1042 
1044  } // namespace PDELab
1045 } // namespace Dune
1046 
1047 #endif // DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
static const int dim
Definition: adaptivity.hh:84
const std::string s
Definition: function.hh:843
const P & p
Definition: constraints.hh:148
const IG & ig
Definition: constraints.hh:149
CL & cl
Definition: constraints.hh:150
XG & xg
Definition: interpolate.hh:55
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:749
void set_nonconstrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:960
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
void constrain_residual(const CG &cg, XG &xg)
transform residual into transformed basis: r -> r~
Definition: constraints.hh:904
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1014
bool check_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg, const Cmp &cmp=Cmp())
check that constrained dofs match a certain value
Definition: constraints.hh:839
void copy_constrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:936
void copy_nonconstrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:987
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
@ processor
processor boundary intersection (neighbor() == false && boundary() == false) or outside entity not in...
@ boundary
domain boundary intersection (neighbor() == false && boundary() == true)
@ skeleton
skeleton intersection (neighbor() == true && boundary() == false)
@ periodic
periodic boundary intersection (neighbor() == true && boundary() == true)
std::tuple< IntersectionType, typename EntitySet::Element > classifyIntersection(const EntitySet &entity_set, const Intersection &is)
Classifies the type of an intersection wrt to the passed EntitySet.
Definition: intersectiontype.hh:37
Dune::TypeTree::GenericLeafNodeTransformation< LeafNode, GridFunctionToLocalViewTransformation, Imp::LocalGridViewFunctionAdapter< LeafNode > > registerNodeTransformation(LeafNode *l, GridFunctionToLocalViewTransformation *t, GridFunctionTag *tag)
Definition: constraints.hh:333
TypeTree::CompositeNode< Children... > BaseT
Definition: constraints.hh:334
CompositeConstraintsOperator(Children &... children)
Definition: constraints.hh:336
CompositeConstraintsOperator(std::shared_ptr< Children >... children)
Definition: constraints.hh:340
Definition: constraints.hh:353
TypeTree::CompositeNode< Children... > BaseT
Definition: constraints.hh:354
CompositeConstraintsParameters(std::shared_ptr< Children >... children)
Definition: constraints.hh:360
CompositeConstraintsParameters(Children &... children)
Definition: constraints.hh:356
Definition: constraints.hh:368
PowerConstraintsParameters(T &c)
Definition: constraints.hh:375
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8, T &c9)
Definition: constraints.hh:447
PowerConstraintsParameters(T &c0, T &c1)
Definition: constraints.hh:379
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3)
Definition: constraints.hh:390
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6)
Definition: constraints.hh:414
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5)
Definition: constraints.hh:405
PowerConstraintsParameters()
Definition: constraints.hh:371
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8)
Definition: constraints.hh:435
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7)
Definition: constraints.hh:424
TypeTree::PowerNode< T, k > BaseT
Definition: constraints.hh:369
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4)
Definition: constraints.hh:397
PowerConstraintsParameters(const std::array< std::shared_ptr< T >, k > &children)
Definition: constraints.hh:460
PowerConstraintsParameters(T &c0, T &c1, T &c2)
Definition: constraints.hh:384
Definition: constraintsparameters.hh:293
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139