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;
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;
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;
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();
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::dimension-1> & 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::dimension-1> & 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;
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();
608 auto& is = es.indexSet();
611 for (
const auto& element : elements(es))
614 auto id = is.uniqueIndex(element);
619 using CL =
typename CG::LocalTransformation;
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())
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())
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 void set_nonconstrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:962
PowerConstraintsParameters()
Definition: constraints.hh:373
const IG & ig
Definition: constraints.hh:148
periodic boundary intersection (neighbor() == true && boundary() == true)
static const int dim
Definition: adaptivity.hh:83
CL & cl
Definition: constraints.hh:149
PowerConstraintsParameters(T &c0, T &c1)
Definition: constraints.hh:381
Wrap intersection.
Definition: geometrywrapper.hh:56
void copy_nonconstrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:989
domain boundary intersection (neighbor() == false && boundary() == true)
Definition: constraints.hh:333
XG & xg
Definition: interpolate.hh:71
const std::string s
Definition: function.hh:830
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1016
CompositeConstraintsOperator(Children &... children)
Definition: constraints.hh:338
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:798
skeleton intersection (neighbor() == true && boundary() == false)
Wrap element.
Definition: geometrywrapper.hh:15
Definition: function.hh:287
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3)
Definition: constraints.hh:392
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
void update()
Definition: lfsindexcache.hh:300
PowerConstraintsParameters(T &c0, T &c1, T &c2)
Definition: constraints.hh:386
TypeTree::CompositeNode< Children... > BaseT
Definition: constraints.hh:356
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5)
Definition: constraints.hh:407
CompositeConstraintsParameters(std::shared_ptr< Children >... children)
Definition: constraints.hh:362
composite functions
Definition: function.hh:528
Dune::TypeTree::GenericLeafNodeTransformation< LeafNode, GridFunctionToLocalViewTransformation, Imp::LocalGridViewFunctionAdapter< LeafNode > > registerNodeTransformation(LeafNode *l, GridFunctionToLocalViewTransformation *t, GridFunctionTag *tag)
PowerConstraintsParameters(T &c)
Definition: constraints.hh:377
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8)
Definition: constraints.hh:437
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
const P & p
Definition: constraints.hh:147
Definition: function.hh:356
Definition: constraintsparameters.hh:293
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:841
void constrain_residual(const CG &cg, XG &xg)
transform residual into transformed basis: r -> r~
Definition: constraints.hh:906
TypeTree::CompositeNode< Children... > BaseT
Definition: constraints.hh:336
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7)
Definition: constraints.hh:426
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4)
Definition: constraints.hh:399
PowerConstraintsParameters(const std::array< std::shared_ptr< T >, k > &children)
Definition: constraints.hh:462
Definition: function.hh:516
Definition: constraints.hh:368
void copy_constrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:938
TypeTree::PowerNode< T, k > BaseT
Definition: constraints.hh:371
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6)
Definition: constraints.hh:416
CompositeConstraintsOperator(std::shared_ptr< Children >... children)
Definition: constraints.hh:342
product of identical functions
Definition: function.hh:366
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:751
CompositeConstraintsParameters(Children &... children)
Definition: constraints.hh:358
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
processor boundary intersection (neighbor() == false && boundary() == false) or outside entity not in...
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8, T &c9)
Definition: constraints.hh:449
Definition: constraints.hh:353