4 #ifndef DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
5 #define DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/float_cmp.hh>
28 template<
typename C,
bool doIt>
29 struct ConstraintsCallBoundary
31 template<
typename P,
typename IG,
typename LFS,
typename T>
32 static void boundary (
const C& c,
const P&
p,
const IG&
ig,
const LFS& lfs, T& trafo)
36 template<
typename C,
bool doIt>
37 struct ConstraintsCallProcessor
39 template<
typename IG,
typename LFS,
typename T>
40 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
44 template<
typename C,
bool doIt>
45 struct ConstraintsCallSkeleton
47 template<
typename IG,
typename LFS,
typename T>
48 static void skeleton (
const C& c,
const IG&
ig,
49 const LFS& lfs_e,
const LFS& lfs_f,
50 T& trafo_e, T& trafo_f)
54 template<
typename C,
bool doIt>
55 struct ConstraintsCallVolume
57 template<
typename P,
typename EG,
typename LFS,
typename T>
58 static void volume (
const C& c,
const P&,
const EG& eg,
const LFS& lfs, T& trafo)
65 struct ConstraintsCallBoundary<C,true>
67 template<
typename P,
typename IG,
typename LFS,
typename T>
68 static void boundary (
const C& c,
const P&
p,
const IG&
ig,
const LFS& lfs, T& trafo)
71 c.boundary(
p,
ig,lfs,trafo);
75 struct ConstraintsCallProcessor<C,true>
77 template<
typename IG,
typename LFS,
typename T>
78 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
81 c.processor(
ig,lfs,trafo);
85 struct ConstraintsCallSkeleton<C,true>
87 template<
typename IG,
typename LFS,
typename T>
88 static void skeleton (
const C& c,
const IG&
ig,
89 const LFS& lfs_e,
const LFS& lfs_f,
90 T& trafo_e, T& trafo_f)
92 if (lfs_e.size() || lfs_f.size())
93 c.skeleton(
ig, lfs_e, lfs_f, trafo_e, trafo_f);
97 struct ConstraintsCallVolume<C,true>
99 template<
typename P,
typename EG,
typename LFS,
typename T>
100 static void volume (
const C& c,
const P&
p,
const EG& eg,
const LFS& lfs, T& trafo)
103 c.volume(
p,eg,lfs,trafo);
108 struct ParameterizedConstraintsBase
109 :
public TypeTree::TreePairVisitor
115 template<
typename P,
typename LFS,
typename TreePath>
116 void leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
118 static_assert((AlwaysFalse<P>::Value),
119 "unsupported combination of function and LocalFunctionSpace");
124 template<
typename P,
typename IG,
typename CL>
125 struct BoundaryConstraintsForParametersLeaf
126 :
public TypeTree::TreeVisitor
127 ,
public TypeTree::DynamicTraversal
130 template<
typename LFS,
typename TreePath>
131 void leaf(
const LFS& lfs, TreePath treePath)
const
135 typedef typename LFS::Traits::ConstraintsType C;
138 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),
p,
ig,lfs,
cl);
141 BoundaryConstraintsForParametersLeaf(
const P& p_,
const IG& ig_, CL& cl_)
154 template<
typename IG,
typename CL>
155 struct BoundaryConstraints
156 :
public ParameterizedConstraintsBase
157 ,
public TypeTree::DynamicTraversal
161 template<
typename P,
typename LFS,
typename TreePath>
162 typename std::enable_if<P::isLeaf && LFS::isLeaf>::type
163 leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
166 typedef typename LFS::Traits::ConstraintsType C;
169 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),
p,
ig,lfs,
cl);
173 template<
typename P,
typename LFS,
typename TreePath>
174 typename std::enable_if<P::isLeaf && (!LFS::isLeaf)>::type
175 leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
178 TypeTree::applyToTree(lfs,BoundaryConstraintsForParametersLeaf<P,IG,CL>(
p,
ig,
cl));
181 BoundaryConstraints(
const IG& ig_, CL& cl_)
193 template<
typename IG,
typename CL>
194 struct ProcessorConstraints
195 :
public TypeTree::TreeVisitor
196 ,
public TypeTree::DynamicTraversal
199 template<
typename LFS,
typename TreePath>
200 void leaf(
const LFS& lfs, TreePath treePath)
const
203 typedef typename LFS::Traits::ConstraintsType C;
206 ConstraintsCallProcessor<C,C::doProcessor>::processor(lfs.constraints(),
ig,lfs,
cl);
209 ProcessorConstraints(
const IG& ig_, CL& cl_)
221 template<
typename IG,
typename CL>
222 struct SkeletonConstraints
223 :
public TypeTree::TreePairVisitor
224 ,
public TypeTree::DynamicTraversal
227 template<
typename LFS,
typename TreePath>
228 void leaf(
const LFS& lfs_e,
const LFS& lfs_f, TreePath treePath)
const
231 typedef typename LFS::Traits::ConstraintsType C;
236 const C & c = lfs_e.constraints();
239 ConstraintsCallSkeleton<C,C::doSkeleton>::skeleton(c,
ig,lfs_e,lfs_f,cl_e,cl_f);
242 SkeletonConstraints(
const IG& ig_, CL& cl_e_, CL& cl_f_)
256 template<
typename P,
typename EG,
typename CL>
257 struct VolumeConstraintsForParametersLeaf
258 :
public TypeTree::TreeVisitor
259 ,
public TypeTree::DynamicTraversal
262 template<
typename LFS,
typename TreePath>
263 void leaf(
const LFS& lfs, TreePath treePath)
const
266 typedef typename LFS::Traits::ConstraintsType C;
267 const C & c = lfs.constraints();
270 ConstraintsCallVolume<C,C::doVolume>::volume(c,
p,eg,lfs,
cl);
273 VolumeConstraintsForParametersLeaf(
const P& p_,
const EG& eg_, CL& cl_)
287 template<
typename EG,
typename CL>
288 struct VolumeConstraints
289 :
public ParameterizedConstraintsBase
290 ,
public TypeTree::DynamicTraversal
294 template<
typename P,
typename LFS,
typename TreePath>
295 typename std::enable_if<P::isLeaf && LFS::isLeaf>::type
296 leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
302 typedef typename LFS::Traits::ConstraintsType C;
303 const C & c = lfs.constraints();
306 ConstraintsCallVolume<C,C::doVolume>::volume(c,
p,eg,lfs,
cl);
310 template<
typename P,
typename LFS,
typename TreePath>
311 typename std::enable_if<P::isLeaf && (!LFS::isLeaf)>::type
312 leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
315 TypeTree::applyToTree(lfs,VolumeConstraintsForParametersLeaf<P,EG,CL>(
p,eg,
cl));
318 VolumeConstraints(
const EG& eg_, CL& cl_)
332 template<
typename... Children>
334 :
public TypeTree::CompositeNode<Children...>
336 typedef TypeTree::CompositeNode<Children...>
BaseT;
352 template<
typename... Children>
354 :
public TypeTree::CompositeNode<Children...>
356 typedef TypeTree::CompositeNode<Children...>
BaseT;
367 template<
typename T, std::
size_t k>
369 :
public TypeTree::PowerNode<T,k>
371 typedef TypeTree::PowerNode<T,k>
BaseT;
404 :
BaseT(c0,c1,c2,c3,c4)
413 :
BaseT(c0,c1,c2,c3,c4,c5)
423 :
BaseT(c0,c1,c2,c3,c4,c5,c6)
434 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
446 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
459 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
471 class OldStyleConstraintsWrapper
472 :
public TypeTree::LeafNode
474 std::shared_ptr<const F> _f;
478 template<
typename Transformation>
479 OldStyleConstraintsWrapper(std::shared_ptr<const F> f,
const Transformation& t,
unsigned int i=0)
484 template<
typename Transformation>
485 OldStyleConstraintsWrapper(
const F & f,
const Transformation& t,
unsigned int i=0)
486 : _f(stackobject_to_shared_ptr(f))
491 bool isDirichlet(
const I & intersection,
const FieldVector<typename I::ctype, I::mydimension> & coord)
const
493 typename F::Traits::RangeType bctype;
494 _f->evaluate(intersection,coord,bctype);
495 return bctype[_i] > 0;
499 bool isNeumann(
const I & intersection,
const FieldVector<typename I::ctype, I::mydimension> & coord)
const
501 typename F::Traits::RangeType bctype;
502 _f->evaluate(intersection,coord,bctype);
503 return bctype[_i] == 0;
508 struct gf_to_constraints {};
511 template<
typename F,
typename Transformation>
512 struct MultiComponentOldStyleConstraintsWrapperDescription
515 static const bool recursive =
false;
517 enum {
dim = F::Traits::dimRange };
518 typedef OldStyleConstraintsWrapper<F> node_type;
519 typedef PowerConstraintsParameters<node_type, dim> transformed_type;
520 typedef std::shared_ptr<transformed_type> transformed_storage_type;
522 static transformed_type transform(
const F&
s,
const Transformation& t)
524 std::shared_ptr<const F> sp = stackobject_to_shared_ptr(
s);
525 std::array<std::shared_ptr<node_type>,
dim> childs;
526 for (
int i=0; i<
dim; i++)
527 childs[i] = std::make_shared<node_type>(sp,t,i);
528 return transformed_type(childs);
531 static transformed_storage_type transform_storage(std::shared_ptr<const F>
s,
const Transformation& t)
533 std::array<std::shared_ptr<node_type>,
dim> childs;
534 for (
int i=0; i<
dim; i++)
535 childs[i] = std::make_shared<node_type>(
s,t,i);
536 return std::make_shared<transformed_type>(childs);
541 template<
typename Gr
idFunction>
542 typename std::conditional<
543 (GridFunction::Traits::dimRange == 1),
545 TypeTree::GenericLeafNodeTransformation<GridFunction,gf_to_constraints,OldStyleConstraintsWrapper<GridFunction> >,
547 MultiComponentOldStyleConstraintsWrapperDescription<GridFunction,gf_to_constraints>
552 template<
typename PowerGr
idFunction>
553 TypeTree::SimplePowerNodeTransformation<PowerGridFunction,gf_to_constraints,PowerConstraintsParameters>
557 template<
typename CompositeGr
idFunction>
558 TypeTree::SimpleCompositeNodeTransformation<CompositeGridFunction,gf_to_constraints,CompositeConstraintsParameters>
572 template<
typename P,
typename GFS,
typename CG,
bool isFunction>
573 struct ConstraintsAssemblerHelper
591 assemble(
const P&
p,
const GFS& gfs, CG& cg,
const bool verbose)
594 using ES =
typename GFS::Traits::EntitySet;
595 using Element =
typename ES::Traits::Element;
596 using Intersection =
typename ES::Traits::Intersection;
598 ES es = gfs.entitySet();
601 using LFS = LocalFunctionSpace<GFS>;
603 LFSIndexCache<LFS> lfs_cache_e(lfs_e);
605 LFSIndexCache<LFS> lfs_cache_f(lfs_f);
608 auto& is = es.indexSet();
611 for (
const auto& element : elements(es))
614 auto id = is.uniqueIndex(element);
619 using CL =
typename CG::LocalTransformation;
623 using ElementWrapper = ElementGeometry<Element>;
624 using IntersectionWrapper = IntersectionGeometry<Intersection>;
626 TypeTree::applyToTreePair(
p,lfs_e,VolumeConstraints<ElementWrapper,CL>(ElementWrapper(element),cl_self));
629 unsigned int intersection_index = 0;
630 for (
const auto& intersection : intersections(es,element))
634 auto intersection_type = std::get<0>(intersection_data);
635 auto& outside_element = std::get<1>(intersection_data);
637 switch (intersection_type) {
642 auto idn = is.uniqueIndex(outside_element);
646 lfs_f.bind(outside_element);
650 TypeTree::applyToTreePair(lfs_e,lfs_f,SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
652 if (!cl_neighbor.empty())
654 lfs_cache_f.update();
655 cg.import_local_transformation(cl_neighbor,lfs_cache_f);
663 TypeTree::applyToTreePair(
p,lfs_e,BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
667 TypeTree::applyToTree(lfs_e,ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
671 ++intersection_index;
674 if (!cl_self.empty())
676 lfs_cache_e.update();
677 cg.import_local_transformation(cl_self,lfs_cache_e);
684 std::cout <<
"constraints:" << std::endl;
686 std::cout << cg.size() <<
" constrained degrees of freedom" << std::endl;
688 for (
const auto& col : cg)
690 std::cout << col.first <<
": ";
691 for (
const auto& row : col.second)
692 std::cout <<
"(" << row.first <<
"," << row.second <<
") ";
693 std::cout << std::endl;
702 template<
typename F,
typename GFS>
703 struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, true>
705 static void assemble(
const F& f,
const GFS& gfs, EmptyTransformation& cg,
const bool verbose)
710 template<
typename F,
typename GFS>
711 struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, false>
713 static void assemble(
const F& f,
const GFS& gfs, EmptyTransformation& cg,
const bool verbose)
720 template<
typename F,
typename GFS,
typename CG>
721 struct ConstraintsAssemblerHelper<F, GFS, CG, true>
724 assemble(
const F& f,
const GFS& gfs, CG& cg,
const bool verbose)
727 typedef typename TypeTree::TransformTree<F,gf_to_constraints> Transformation;
728 typedef typename Transformation::Type P;
730 P
p = Transformation::transform(f);
750 template<
typename GFS,
typename CG>
752 const bool verbose =
false)
755 ConstraintsAssemblerHelper<NoConstraintsParameters, GFS, CG, false>::assemble(
p,gfs,cg,verbose);
776 template<
typename P,
typename GFS,
typename CG>
778 const bool verbose =
false)
797 template<
typename CG,
typename XG>
799 typename XG::ElementType x,
802 for (
const auto& col : cg)
810 template<
typename XG>
812 typename XG::ElementType x,
840 template<
typename CG,
typename XG,
typename Cmp>
842 XG&
xg,
const Cmp& cmp = Cmp())
844 for (
const auto& col : cg)
845 if(cmp.ne(
xg[col.first], x))
869 template<
typename CG,
typename XG>
874 FloatCmpOps<typename XG::ElementType>());
881 template<
typename XG,
typename Cmp>
883 XG&
xg,
const Cmp& cmp = Cmp())
889 template<
typename XG>
905 template<
typename CG,
typename XG>
908 for (
const auto& col : cg)
909 for (
const auto& row : col.second)
910 xg[row.first] += row.second *
xg[col.first];
914 for (
const auto& col : cg)
915 xg[col.first] =
typename XG::ElementType(0);
922 template<
typename XG>
937 template<
typename CG,
typename XG>
940 for (
const auto& col : cg)
941 xgout[col.first] = xgin[col.first];
948 template<
typename XG>
961 template<
typename CG,
typename XG>
973 template<
typename XG>
988 template<
typename CG,
typename XG>
1000 template<
typename XG>
1015 template<
typename CG,
typename XG>
1021 for (
const auto& col : cg)
1024 if (col.second.size() == 0)
1026 tmp[col.first] =
xg[col.first];
1037 template<
typename XG>
1038 void set_shifted_dofs (
const EmptyTransformation& cg,
typename XG::ElementType x, XG&
xg)
1049 #endif // DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH