4#ifndef DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
5#define DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
11#include<dune/pdelab/common/function.hh>
12#include<dune/pdelab/common/geometrywrapper.hh>
13#include<dune/pdelab/common/typetraits.hh>
14#include<dune/pdelab/common/intersectiontype.hh>
15#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
16#include"constraintstransformation.hh"
17#include"constraintsparameters.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
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
316 VolumeConstraints(
const EG& eg_, CL& cl_)
330 template<
typename... Children>
331 struct CompositeConstraintsOperator
332 :
public TypeTree::CompositeNode<Children...>
334 typedef TypeTree::CompositeNode<Children...> BaseT;
336 CompositeConstraintsOperator(Children&... children)
340 CompositeConstraintsOperator(std::shared_ptr<Children>... children)
350 template<
typename... Children>
351 struct CompositeConstraintsParameters
352 :
public TypeTree::CompositeNode<Children...>
354 typedef TypeTree::CompositeNode<Children...> BaseT;
356 CompositeConstraintsParameters(Children&... children)
360 CompositeConstraintsParameters(std::shared_ptr<Children>... children)
365 template<
typename T, std::
size_t k>
366 struct PowerConstraintsParameters
367 :
public TypeTree::PowerNode<T,k>
369 typedef TypeTree::PowerNode<T,k> BaseT;
371 PowerConstraintsParameters()
375 PowerConstraintsParameters(T& c)
379 PowerConstraintsParameters (T& c0,
384 PowerConstraintsParameters (T& c0,
390 PowerConstraintsParameters (T& c0,
397 PowerConstraintsParameters (T& c0,
402 : BaseT(c0,c1,c2,c3,c4)
405 PowerConstraintsParameters (T& c0,
411 : BaseT(c0,c1,c2,c3,c4,c5)
414 PowerConstraintsParameters (T& c0,
421 : BaseT(c0,c1,c2,c3,c4,c5,c6)
424 PowerConstraintsParameters (T& c0,
432 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
435 PowerConstraintsParameters (T& c0,
444 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
447 PowerConstraintsParameters (T& c0,
457 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
460 PowerConstraintsParameters (
const std::array<std::shared_ptr<T>,k>& children)
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)
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)
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>;
627 unsigned int intersection_index = 0;
628 for (
const auto& intersection : intersections(es,element))
631 auto intersection_data = classifyIntersection(es,intersection);
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>
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)
Class encapsulating a default epsilon.
Definition: float_cmp.hh:288
Definition: constraintsparameters.hh:293
A few common exception classes.
Various ways to compare floating-point numbers.
bool check_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
check that constrained dofs match a certain value
Definition: constraints.hh:868
void set_nonconstrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:960
void constraints(const P &p, const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints from given constraints parameter tree
Definition: constraints.hh:775
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
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
constexpr auto treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:326
void applyToTreePair(Tree1 &&tree1, Tree2 &&tree2, Visitor &&visitor)
Apply visitor to a pair of TypeTrees.
Definition: pairtraversal.hh:128
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:239
Dune namespace.
Definition: alignedallocator.hh:13
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:72
This file implements several utilities related to std::shared_ptr.