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