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/shared_ptr.hh>
9#include<dune/common/float_cmp.hh>
29 template<
typename C,
bool doIt>
30 struct ConstraintsCallBoundary
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)
37 template<
typename C,
bool doIt>
38 struct ConstraintsCallProcessor
40 template<
typename IG,
typename LFS,
typename T>
41 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
45 template<
typename C,
bool doIt>
46 struct ConstraintsCallSkeleton
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)
55 template<
typename C,
bool doIt>
56 struct ConstraintsCallVolume
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)
66 struct ConstraintsCallBoundary<C,true>
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)
72 c.boundary(
p,
ig,lfs,trafo);
76 struct ConstraintsCallProcessor<C,true>
78 template<
typename IG,
typename LFS,
typename T>
79 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
82 c.processor(
ig,lfs,trafo);
86 struct ConstraintsCallSkeleton<C,true>
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)
93 if (lfs_e.size() || lfs_f.size())
94 c.skeleton(
ig, lfs_e, lfs_f, trafo_e, trafo_f);
98 struct ConstraintsCallVolume<C,true>
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)
104 c.volume(
p,eg,lfs,trafo);
109 struct ParameterizedConstraintsBase
110 :
public TypeTree::TreePairVisitor
116 template<
typename P,
typename LFS,
typename TreePath>
117 void leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
119 static_assert((AlwaysFalse<P>::Value),
120 "unsupported combination of function and LocalFunctionSpace");
125 template<
typename P,
typename IG,
typename CL>
126 struct BoundaryConstraintsForParametersLeaf
127 :
public TypeTree::TreeVisitor
128 ,
public TypeTree::DynamicTraversal
131 template<
typename LFS,
typename TreePath>
132 void leaf(
const LFS& lfs, TreePath treePath)
const
136 typedef typename LFS::Traits::ConstraintsType C;
139 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),p,ig,lfs,cl);
142 BoundaryConstraintsForParametersLeaf(
const P& p_,
const IG& ig_, CL& cl_)
155 template<
typename IG,
typename CL>
156 struct BoundaryConstraints
157 :
public ParameterizedConstraintsBase
158 ,
public TypeTree::DynamicTraversal
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
167 typedef typename LFS::Traits::ConstraintsType C;
170 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),
p,
ig,lfs,
cl);
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
179 TypeTree::applyToTree(lfs,BoundaryConstraintsForParametersLeaf<P,IG,CL>(
p,ig,cl));
182 BoundaryConstraints(
const IG& ig_, CL& cl_)
194 template<
typename IG,
typename CL>
195 struct ProcessorConstraints
196 :
public TypeTree::TreeVisitor
197 ,
public TypeTree::DynamicTraversal
200 template<
typename LFS,
typename TreePath>
201 void leaf(
const LFS& lfs, TreePath treePath)
const
204 typedef typename LFS::Traits::ConstraintsType C;
207 ConstraintsCallProcessor<C,C::doProcessor>::processor(lfs.constraints(),ig,lfs,cl);
210 ProcessorConstraints(
const IG& ig_, CL& cl_)
222 template<
typename IG,
typename CL>
223 struct SkeletonConstraints
224 :
public TypeTree::TreePairVisitor
225 ,
public TypeTree::DynamicTraversal
228 template<
typename LFS,
typename TreePath>
229 void leaf(
const LFS& lfs_e,
const LFS& lfs_f, TreePath treePath)
const
232 typedef typename LFS::Traits::ConstraintsType C;
237 const C & c = lfs_e.constraints();
240 ConstraintsCallSkeleton<C,C::doSkeleton>::skeleton(c,ig,lfs_e,lfs_f,cl_e,cl_f);
243 SkeletonConstraints(
const IG& ig_, CL& cl_e_, CL& cl_f_)
257 template<
typename P,
typename EG,
typename CL>
258 struct VolumeConstraintsForParametersLeaf
259 :
public TypeTree::TreeVisitor
260 ,
public TypeTree::DynamicTraversal
263 template<
typename LFS,
typename TreePath>
264 void leaf(
const LFS& lfs, TreePath treePath)
const
267 typedef typename LFS::Traits::ConstraintsType C;
268 const C & c = lfs.constraints();
271 ConstraintsCallVolume<C,C::doVolume>::volume(c,p,eg,lfs,cl);
274 VolumeConstraintsForParametersLeaf(
const P& p_,
const EG& eg_, CL& cl_)
288 template<
typename EG,
typename CL>
289 struct VolumeConstraints
290 :
public ParameterizedConstraintsBase
291 ,
public TypeTree::DynamicTraversal
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
300 typedef typename LFS::Traits::ConstraintsType C;
301 const C & c = lfs.constraints();
304 ConstraintsCallVolume<C,C::doVolume>::volume(c,
p,eg,lfs,cl);
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
313 TypeTree::applyToTree(lfs,VolumeConstraintsForParametersLeaf<P,EG,CL>(
p,eg,cl));
316 VolumeConstraints(
const EG& eg_, CL& cl_)
330 template<
typename... Children>
332 :
public TypeTree::CompositeNode<Children...>
334 typedef TypeTree::CompositeNode<Children...>
BaseT;
350 template<
typename... Children>
352 :
public TypeTree::CompositeNode<Children...>
354 typedef TypeTree::CompositeNode<Children...>
BaseT;
365 template<
typename T, std::
size_t k>
367 :
public TypeTree::PowerNode<T,k>
369 typedef TypeTree::PowerNode<T,k>
BaseT;
402 :
BaseT(c0,c1,c2,c3,c4)
411 :
BaseT(c0,c1,c2,c3,c4,c5)
421 :
BaseT(c0,c1,c2,c3,c4,c5,c6)
432 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
444 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
457 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
469 class OldStyleConstraintsWrapper
470 :
public TypeTree::LeafNode
472 std::shared_ptr<const F> _f;
476 template<
typename Transformation>
477 OldStyleConstraintsWrapper(std::shared_ptr<const F> f,
const Transformation& t,
unsigned int i=0)
482 template<
typename Transformation>
483 OldStyleConstraintsWrapper(
const F & f,
const Transformation& t,
unsigned int i=0)
484 : _f(stackobject_to_shared_ptr(f))
489 bool isDirichlet(
const I & intersection,
const FieldVector<typename I::ctype, I::mydimension> & coord)
const
491 typename F::Traits::RangeType bctype;
492 _f->evaluate(intersection,coord,bctype);
493 return bctype[_i] > 0;
497 bool isNeumann(
const I & intersection,
const FieldVector<typename I::ctype, I::mydimension> & coord)
const
499 typename F::Traits::RangeType bctype;
500 _f->evaluate(intersection,coord,bctype);
501 return bctype[_i] == 0;
506 struct gf_to_constraints {};
509 template<
typename F,
typename Transformation>
510 struct MultiComponentOldStyleConstraintsWrapperDescription
513 static const bool recursive =
false;
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;
520 static transformed_type transform(
const F&
s,
const Transformation& t)
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);
529 static transformed_storage_type transform_storage(std::shared_ptr<const F>
s,
const Transformation& t)
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);
539 template<
typename Gr
idFunction>
540 typename std::conditional<
541 (GridFunction::Traits::dimRange == 1),
543 TypeTree::GenericLeafNodeTransformation<GridFunction,gf_to_constraints,OldStyleConstraintsWrapper<GridFunction> >,
545 MultiComponentOldStyleConstraintsWrapperDescription<GridFunction,gf_to_constraints>
550 template<
typename PowerGr
idFunction>
551 TypeTree::SimplePowerNodeTransformation<PowerGridFunction,gf_to_constraints,PowerConstraintsParameters>
555 template<
typename CompositeGr
idFunction>
556 TypeTree::SimpleCompositeNodeTransformation<CompositeGridFunction,gf_to_constraints,CompositeConstraintsParameters>
570 template<
typename P,
typename GFS,
typename CG,
bool isFunction>
571 struct ConstraintsAssemblerHelper
589 assemble(
const P&
p,
const GFS& gfs, CG& cg,
const bool verbose)
592 using ES =
typename GFS::Traits::EntitySet;
593 using Element =
typename ES::Traits::Element;
594 using Intersection =
typename ES::Traits::Intersection;
596 ES es = gfs.entitySet();
599 using LFS = LocalFunctionSpace<GFS>;
601 LFSIndexCache<LFS> lfs_cache_e(lfs_e);
603 LFSIndexCache<LFS> lfs_cache_f(lfs_f);
606 auto& is = es.indexSet();
609 for (
const auto& element : elements(es))
612 auto id = is.uniqueIndex(element);
617 using CL =
typename CG::LocalTransformation;
621 using ElementWrapper = ElementGeometry<Element>;
622 using IntersectionWrapper = IntersectionGeometry<Intersection>;
624 TypeTree::applyToTreePair(
p,lfs_e,VolumeConstraints<ElementWrapper,CL>(ElementWrapper(element),cl_self));
627 unsigned int intersection_index = 0;
628 for (
const auto& intersection : intersections(es,element))
632 auto intersection_type = std::get<0>(intersection_data);
633 auto& outside_element = std::get<1>(intersection_data);
635 switch (intersection_type) {
637 case IntersectionType::skeleton:
638 case IntersectionType::periodic:
640 auto idn = is.uniqueIndex(outside_element);
644 lfs_f.bind(outside_element);
648 TypeTree::applyToTreePair(lfs_e,lfs_f,SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
650 if (!cl_neighbor.empty())
652 lfs_cache_f.update();
653 cg.import_local_transformation(cl_neighbor,lfs_cache_f);
660 case IntersectionType::boundary:
661 TypeTree::applyToTreePair(
p,lfs_e,BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
664 case IntersectionType::processor:
665 TypeTree::applyToTree(lfs_e,ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
669 ++intersection_index;
672 if (!cl_self.empty())
674 lfs_cache_e.update();
675 cg.import_local_transformation(cl_self,lfs_cache_e);
682 std::cout <<
"constraints:" << std::endl;
684 std::cout << cg.size() <<
" constrained degrees of freedom" << std::endl;
686 for (
const auto& col : cg)
688 std::cout << col.first <<
": ";
689 for (
const auto& row : col.second)
690 std::cout <<
"(" << row.first <<
"," << row.second <<
") ";
691 std::cout << std::endl;
700 template<
typename F,
typename GFS>
701 struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, true>
703 static void assemble(
const F& f,
const GFS& gfs, EmptyTransformation& cg,
const bool verbose)
708 template<
typename F,
typename GFS>
709 struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, false>
711 static void assemble(
const F& f,
const GFS& gfs, EmptyTransformation& cg,
const bool verbose)
718 template<
typename F,
typename GFS,
typename CG>
719 struct ConstraintsAssemblerHelper<F, GFS, CG, true>
722 assemble(
const F& f,
const GFS& gfs, CG& cg,
const bool verbose)
725 typedef typename TypeTree::TransformTree<F,gf_to_constraints> Transformation;
726 typedef typename Transformation::Type P;
728 P
p = Transformation::transform(f);
730 ConstraintsAssemblerHelper<P, GFS, CG, IsGridFunction<P>::value>::assemble(
p,gfs,cg,verbose);
748 template<
typename GFS,
typename CG>
750 const bool verbose =
false)
753 ConstraintsAssemblerHelper<NoConstraintsParameters, GFS, CG, false>::assemble(
p,gfs,cg,verbose);
774 template<
typename P,
typename GFS,
typename CG>
776 const bool verbose =
false)
780 ConstraintsAssemblerHelper<P, GFS, CG, IsGridFunction<P>::value>::assemble(
p,gfs,cg,verbose);
795 template<
typename CG,
typename XG>
797 typename XG::ElementType x,
800 for (
const auto& col : cg)
808 template<
typename XG>
810 typename XG::ElementType x,
838 template<
typename CG,
typename XG,
typename Cmp>
840 XG&
xg,
const Cmp& cmp = Cmp())
842 for (
const auto& col : cg)
843 if(cmp.ne(
xg[col.first], x))
867 template<
typename CG,
typename XG>
872 FloatCmpOps<typename XG::ElementType>());
879 template<
typename XG,
typename Cmp>
881 XG&
xg,
const Cmp& cmp = Cmp())
887 template<
typename XG>
903 template<
typename CG,
typename XG>
906 for (
const auto& col : cg)
907 for (
const auto& row : col.second)
908 xg[row.first] += row.second *
xg[col.first];
912 for (
const auto& col : cg)
913 xg[col.first] =
typename XG::ElementType(0);
920 template<
typename XG>
935 template<
typename CG,
typename XG>
938 for (
const auto& col : cg)
939 xgout[col.first] = xgin[col.first];
946 template<
typename XG>
959 template<
typename CG,
typename XG>
971 template<
typename XG>
986 template<
typename CG,
typename XG>
998 template<
typename XG>
1013 template<
typename CG,
typename XG>
1019 for (
const auto& col : cg)
1022 if (col.second.size() == 0)
1024 tmp[col.first] =
xg[col.first];
1035 template<
typename XG>
1036 void set_shifted_dofs (
const EmptyTransformation& cg,
typename XG::ElementType x, XG&
xg)
const std::string s
Definition: function.hh:843
static const int dim
Definition: adaptivity.hh:84
XG & xg
Definition: interpolate.hh:55
const P & p
Definition: constraints.hh:148
const IG & ig
Definition: constraints.hh:149
CL & cl
Definition: constraints.hh:150
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)
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