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