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>
27 template<
typename C,
bool doIt>
28 struct ConstraintsCallBoundary
30 template<
typename P,
typename IG,
typename LFS,
typename T>
31 static void boundary (
const C& c,
const P&
p,
const IG&
ig,
const LFS& lfs, T& trafo)
35 template<
typename C,
bool doIt>
36 struct ConstraintsCallProcessor
38 template<
typename IG,
typename LFS,
typename T>
39 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
43 template<
typename C,
bool doIt>
44 struct ConstraintsCallSkeleton
46 template<
typename IG,
typename LFS,
typename T>
47 static void skeleton (
const C& c,
const IG&
ig,
48 const LFS& lfs_e,
const LFS& lfs_f,
49 T& trafo_e, T& trafo_f)
53 template<
typename C,
bool doIt>
54 struct ConstraintsCallVolume
56 template<
typename P,
typename EG,
typename LFS,
typename T>
57 static void volume (
const C& c,
const P&,
const EG&
eg,
const LFS& lfs, T& trafo)
64 struct ConstraintsCallBoundary<C,true>
66 template<
typename P,
typename IG,
typename LFS,
typename T>
67 static void boundary (
const C& c,
const P&
p,
const IG&
ig,
const LFS& lfs, T& trafo)
70 c.boundary(p,ig,lfs,trafo);
74 struct ConstraintsCallProcessor<C,true>
76 template<
typename IG,
typename LFS,
typename T>
77 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
80 c.processor(ig,lfs,trafo);
84 struct ConstraintsCallSkeleton<C,true>
86 template<
typename IG,
typename LFS,
typename T>
87 static void skeleton (
const C& c,
const IG&
ig,
88 const LFS& lfs_e,
const LFS& lfs_f,
89 T& trafo_e, T& trafo_f)
91 if (lfs_e.size() || lfs_f.size())
92 c.skeleton(ig, lfs_e, lfs_f, trafo_e, trafo_f);
96 struct ConstraintsCallVolume<C,true>
98 template<
typename P,
typename EG,
typename LFS,
typename T>
99 static void volume (
const C& c,
const P&
p,
const EG&
eg,
const LFS& lfs, T& trafo)
102 c.volume(p,eg,lfs,trafo);
107 struct ParameterizedConstraintsBase
108 :
public TypeTree::TreePairVisitor
114 template<
typename P,
typename LFS,
typename TreePath>
115 void leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
117 static_assert((AlwaysFalse<P>::Value),
118 "unsupported combination of function and LocalFunctionSpace");
123 template<
typename P,
typename IG,
typename CL>
124 struct BoundaryConstraintsForParametersLeaf
125 :
public TypeTree::TreeVisitor
126 ,
public TypeTree::DynamicTraversal
129 template<
typename LFS,
typename TreePath>
130 void leaf(
const LFS& lfs, TreePath treePath)
const
134 typedef typename LFS::Traits::ConstraintsType C;
137 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),
p,
ig,lfs,
cl);
140 BoundaryConstraintsForParametersLeaf(
const P& p_,
const IG& ig_, CL& cl_)
153 template<
typename IG,
typename CL>
154 struct BoundaryConstraints
155 :
public ParameterizedConstraintsBase
156 ,
public TypeTree::DynamicTraversal
160 template<
typename P,
typename LFS,
typename TreePath>
161 typename enable_if<P::isLeaf && LFS::isLeaf>::type
162 leaf(
const P&
p,
const LFS& lfs, TreePath treePath)
const
165 typedef typename LFS::Traits::ConstraintsType C;
168 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),
p,
ig,lfs,
cl);
172 template<
typename P,
typename LFS,
typename TreePath>
173 typename enable_if<P::isLeaf && (!LFS::isLeaf)>::type
174 leaf(
const P& p,
const LFS& lfs, TreePath treePath)
const
177 TypeTree::applyToTree(lfs,BoundaryConstraintsForParametersLeaf<P,IG,CL>(p,ig,
cl));
180 BoundaryConstraints(
const IG& ig_, CL& cl_)
192 template<
typename IG,
typename CL>
193 struct ProcessorConstraints
194 :
public TypeTree::TreeVisitor
195 ,
public TypeTree::DynamicTraversal
198 template<
typename LFS,
typename TreePath>
199 void leaf(
const LFS& lfs, TreePath treePath)
const
202 typedef typename LFS::Traits::ConstraintsType C;
205 ConstraintsCallProcessor<C,C::doProcessor>::processor(lfs.constraints(),
ig,lfs,
cl);
208 ProcessorConstraints(
const IG& ig_, CL& cl_)
220 template<
typename IG,
typename CL>
221 struct SkeletonConstraints
222 :
public TypeTree::TreePairVisitor
223 ,
public TypeTree::DynamicTraversal
226 template<
typename LFS,
typename TreePath>
227 void leaf(
const LFS& lfs_e,
const LFS& lfs_f, TreePath treePath)
const
230 typedef typename LFS::Traits::ConstraintsType C;
235 const C & c = lfs_e.constraints();
238 ConstraintsCallSkeleton<C,C::doSkeleton>::skeleton(c,ig,lfs_e,lfs_f,
cl_e,
cl_f);
241 SkeletonConstraints(
const IG& ig_, CL& cl_e_, CL& cl_f_)
255 template<
typename P,
typename EG,
typename CL>
256 struct VolumeConstraintsForParametersLeaf
257 :
public TypeTree::TreeVisitor
258 ,
public TypeTree::DynamicTraversal
261 template<
typename LFS,
typename TreePath>
262 void leaf(
const LFS& lfs, TreePath treePath)
const
265 typedef typename LFS::Traits::ConstraintsType C;
266 const C & c = lfs.constraints();
269 ConstraintsCallVolume<C,C::doVolume>::volume(c,p,
eg,lfs,
cl);
272 VolumeConstraintsForParametersLeaf(
const P& p_,
const EG& eg_, CL& cl_)
286 template<
typename EG,
typename CL>
287 struct VolumeConstraints
288 :
public ParameterizedConstraintsBase
289 ,
public TypeTree::DynamicTraversal
293 template<
typename P,
typename LFS,
typename TreePath>
294 typename enable_if<P::isLeaf && LFS::isLeaf>::type
295 leaf(
const P& p,
const LFS& lfs, TreePath treePath)
const
301 typedef typename LFS::Traits::ConstraintsType C;
302 const C & c = lfs.constraints();
305 ConstraintsCallVolume<C,C::doVolume>::volume(c,p,
eg,lfs,cl);
309 template<
typename P,
typename LFS,
typename TreePath>
310 typename enable_if<P::isLeaf && (!LFS::isLeaf)>::type
311 leaf(
const P& p,
const LFS& lfs, TreePath treePath)
const
314 TypeTree::applyToTree(lfs,VolumeConstraintsForParametersLeaf<P,EG,CL>(p,
eg,
cl));
317 VolumeConstraints(
const EG& eg_, CL& cl_)
331 template<
typename... Children>
333 :
public TypeTree::CompositeNode<Children...>
335 typedef TypeTree::CompositeNode<Children...>
BaseT;
351 template<
typename... Children>
353 :
public TypeTree::CompositeNode<Children...>
355 typedef TypeTree::CompositeNode<Children...>
BaseT;
366 template<
typename T, std::
size_t k>
368 :
public TypeTree::PowerNode<T,k>
370 typedef TypeTree::PowerNode<T,k>
BaseT;
403 : BaseT(c0,c1,c2,c3,c4)
412 : BaseT(c0,c1,c2,c3,c4,c5)
422 : BaseT(c0,c1,c2,c3,c4,c5,c6)
433 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
445 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
458 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
470 class OldStyleConstraintsWrapper
471 :
public TypeTree::LeafNode
473 std::shared_ptr<const F> _f;
477 template<
typename Transformation>
478 OldStyleConstraintsWrapper(std::shared_ptr<const F> f,
const Transformation& t,
unsigned int i=0)
483 template<
typename Transformation>
484 OldStyleConstraintsWrapper(
const F & f,
const Transformation& t,
unsigned int i=0)
485 : _f(stackobject_to_shared_ptr(f))
490 bool isDirichlet(
const I & intersection,
const FieldVector<typename I::ctype, I::dimension-1> & coord)
const
492 typename F::Traits::RangeType bctype;
493 _f->evaluate(intersection,coord,bctype);
494 return bctype[_i] > 0;
498 bool isNeumann(
const I & intersection,
const FieldVector<typename I::ctype, I::dimension-1> & coord)
const
500 typename F::Traits::RangeType bctype;
501 _f->evaluate(intersection,coord,bctype);
502 return bctype[_i] == 0;
507 struct gf_to_constraints {};
510 template<
typename F,
typename Transformation>
511 struct MultiComponentOldStyleConstraintsWrapperDescription
514 static const bool recursive =
false;
516 enum {
dim = F::Traits::dimRange };
517 typedef OldStyleConstraintsWrapper<F> node_type;
518 typedef PowerConstraintsParameters<node_type, dim> transformed_type;
519 typedef std::shared_ptr<transformed_type> transformed_storage_type;
521 static transformed_type transform(
const F&
s,
const Transformation& t)
523 std::shared_ptr<const F> sp = stackobject_to_shared_ptr(s);
524 std::array<std::shared_ptr<node_type>,
dim> childs;
525 for (
int i=0; i<
dim; i++)
526 childs[i] = std::make_shared<node_type>(sp,t,i);
527 return transformed_type(childs);
530 static transformed_storage_type transform_storage(std::shared_ptr<const F> s,
const Transformation& t)
532 std::array<std::shared_ptr<node_type>, dim> childs;
533 for (
int i=0; i<
dim; i++)
534 childs[i] = std::make_shared<node_type>(s,t,i);
535 return std::make_shared<transformed_type>(childs);
540 template<
typename Gr
idFunction>
541 typename conditional<
542 (GridFunction::Traits::dimRange == 1),
544 TypeTree::GenericLeafNodeTransformation<GridFunction,gf_to_constraints,OldStyleConstraintsWrapper<GridFunction> >,
546 MultiComponentOldStyleConstraintsWrapperDescription<GridFunction,gf_to_constraints>
551 template<
typename PowerGr
idFunction>
552 TypeTree::SimplePowerNodeTransformation<PowerGridFunction,gf_to_constraints,PowerConstraintsParameters>
556 template<
typename CompositeGr
idFunction>
557 TypeTree::SimpleCompositeNodeTransformation<CompositeGridFunction,gf_to_constraints,CompositeConstraintsParameters>
571 template<
typename P,
typename GFS,
typename GV,
typename CG,
bool isFunction>
572 struct ConstraintsAssemblerHelper
590 assemble(
const P& p,
const GFS& gfs,
const GV& gv, CG& cg,
const bool verbose)
593 typedef typename GV::Traits::template Codim<0>::Entity Element;
594 typedef typename GV::Intersection Intersection;
597 typedef LocalFunctionSpace<GFS> LFS;
599 LFSIndexCache<LFS> lfs_cache_e(lfs_e);
601 LFSIndexCache<LFS> lfs_cache_f(lfs_f);
604 const typename GV::IndexSet& is=gv.indexSet();
607 const int chunk=1<<28;
609 std::map<Dune::GeometryType,int> gtoffset;
612 for (
const auto& cell : elements(gv))
615 if (gtoffset.find(cell.type())==gtoffset.end())
617 gtoffset[cell.type()] =
offset;
621 const typename GV::IndexSet::IndexType
id = is.index(cell)+gtoffset[cell.type()];
626 typedef typename CG::LocalTransformation CL;
630 typedef ElementGeometry<Element> ElementWrapper;
631 TypeTree::applyToTreePair(p,lfs_e,VolumeConstraints<ElementWrapper,CL>(ElementWrapper(cell),cl_self));
634 unsigned int intersection_index = 0;
635 for (
const auto& intersection : intersections(gv,cell))
637 if (intersection.boundary())
639 typedef IntersectionGeometry<Intersection> IntersectionWrapper;
640 TypeTree::applyToTreePair(p,lfs_e,BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
644 if ((!intersection.boundary()) && (!intersection.neighbor()))
646 typedef IntersectionGeometry<Intersection> IntersectionWrapper;
647 TypeTree::applyToTree(lfs_e,ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
651 if (intersection.neighbor()){
653 auto outside_cell = intersection.outside();
654 Dune::GeometryType gtn = outside_cell.type();
655 const typename GV::IndexSet::IndexType idn = is.index(outside_cell)+gtoffset[gtn];
659 lfs_f.bind(outside_cell);
663 typedef IntersectionGeometry<Intersection> IntersectionWrapper;
664 TypeTree::applyToTreePair(lfs_e,lfs_f,SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
666 if (!cl_neighbor.empty())
668 lfs_cache_f.update();
669 cg.import_local_transformation(cl_neighbor,lfs_cache_f);
674 ++intersection_index;
677 if (!cl_self.empty())
679 lfs_cache_e.update();
680 cg.import_local_transformation(cl_self,lfs_cache_e);
687 std::cout <<
"constraints:" << std::endl;
689 std::cout << cg.size() <<
" constrained degrees of freedom" << std::endl;
691 for (
const auto& col : cg)
693 std::cout << col.first <<
": ";
694 for (
const auto& row : col.second)
695 std::cout <<
"(" << row.first <<
"," << row.second <<
") ";
696 std::cout << std::endl;
705 template<
typename F,
typename GFS,
typename GV>
706 struct ConstraintsAssemblerHelper<F, GFS, GV, EmptyTransformation, true>
708 static void assemble(
const F& f,
const GFS& gfs,
const GV& gv, EmptyTransformation& cg,
const bool verbose)
713 template<
typename F,
typename GFS,
typename GV>
714 struct ConstraintsAssemblerHelper<F, GFS, GV, EmptyTransformation, false>
716 static void assemble(
const F& f,
const GFS& gfs,
const GV& gv, EmptyTransformation& cg,
const bool verbose)
723 template<
typename F,
typename GFS,
typename GV,
typename CG>
724 struct ConstraintsAssemblerHelper<F, GFS, GV, CG, true>
727 assemble(
const F& f,
const GFS& gfs,
const GV& gv, CG& cg,
const bool verbose)
730 typedef typename TypeTree::TransformTree<F,gf_to_constraints> Transformation;
731 typedef typename Transformation::Type P;
733 P p = Transformation::transform(f);
753 template<
typename GFS,
typename CG>
755 const bool verbose =
false)
757 typedef typename GFS::Traits::GridViewType GV;
759 ConstraintsAssemblerHelper<NoConstraintsParameters, GFS, GV, CG, false>::assemble(p,gfs,gfs.gridView(),cg,verbose);
780 template<
typename P,
typename GFS,
typename CG>
782 const bool verbose =
false)
784 typedef typename GFS::Traits::GridViewType GV;
802 template<
typename CG,
typename XG>
804 typename XG::ElementType x,
807 for (
const auto& col : cg)
815 template<
typename XG>
817 typename XG::ElementType x,
845 template<
typename CG,
typename XG,
typename Cmp>
847 XG&
xg,
const Cmp& cmp = Cmp())
849 for (
const auto& col : cg)
850 if(cmp.ne(xg[col.first], x))
874 template<
typename CG,
typename XG>
879 FloatCmpOps<typename XG::ElementType>());
886 template<
typename XG,
typename Cmp>
888 XG& xg,
const Cmp& cmp = Cmp())
894 template<
typename XG>
910 template<
typename CG,
typename XG>
913 for (
const auto& col : cg)
914 for (
const auto& row : col.second)
915 xg[row.first] += row.second * xg[col.first];
919 for (
const auto& col : cg)
920 xg[col.first] =
typename XG::ElementType(0);
927 template<
typename XG>
942 template<
typename CG,
typename XG>
945 for (
const auto& col : cg)
946 xgout[col.first] = xgin[col.first];
953 template<
typename XG>
966 template<
typename CG,
typename XG>
978 template<
typename XG>
993 template<
typename CG,
typename XG>
1005 template<
typename XG>
1020 template<
typename CG,
typename XG>
1026 for (
const auto& col : cg)
1029 if (col.second.size() == 0)
1031 tmp[col.first] = xg[col.first];
1042 template<
typename XG>
1043 void set_shifted_dofs (
const EmptyTransformation& cg,
typename XG::ElementType x, XG& xg)
1054 #endif // DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
const std::size_t offset
Definition: localfunctionspace.hh:74
CL & cl_f
Definition: constraints.hh:250
TypeTree::PowerNode< T, k > BaseT
Definition: constraints.hh:370
Definition: constraints.hh:352
void copy_nonconstrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:994
void copy_constrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:943
XG & xg
Definition: interpolate.hh:67
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7)
Definition: constraints.hh:425
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:803
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4)
Definition: constraints.hh:398
PowerConstraintsParameters(T &c)
Definition: constraints.hh:376
Definition: constraints.hh:367
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6)
Definition: constraints.hh:415
void set_nonconstrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:967
CompositeConstraintsOperator(Children &...children)
Definition: constraints.hh:337
PowerConstraintsParameters(const std::array< std::shared_ptr< T >, k > &children)
Definition: constraints.hh:461
CL & cl
Definition: constraints.hh:148
const IG & ig
Definition: constraints.hh:147
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
CL & cl_e
Definition: constraints.hh:249
PowerConstraintsParameters(T &c0, T &c1)
Definition: constraints.hh:380
TypeTree::CompositeNode< Children...> BaseT
Definition: constraints.hh:355
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3)
Definition: constraints.hh:391
CompositeConstraintsParameters(Children &...children)
Definition: constraints.hh:357
PowerConstraintsParameters(T &c0, T &c1, T &c2)
Definition: constraints.hh:385
CompositeConstraintsParameters(std::shared_ptr< Children >...children)
Definition: constraints.hh:361
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8)
Definition: constraints.hh:436
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:846
TypeTree::CompositeNode< Children...> BaseT
Definition: constraints.hh:335
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:754
Definition: adaptivity.hh:27
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1021
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)
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5)
Definition: constraints.hh:406
static const int dim
Definition: adaptivity.hh:83
PowerConstraintsParameters()
Definition: constraints.hh:372
const std::string s
Definition: function.hh:1102
Definition: constraintsparameters.hh:293
const P & p
Definition: constraints.hh:146
void constrain_residual(const CG &cg, XG &xg)
transform residual into transformed basis: r -> r~
Definition: constraints.hh:911
Definition: constraints.hh:332
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8, T &c9)
Definition: constraints.hh:448
const EG & eg
Definition: constraints.hh:280
CompositeConstraintsOperator(std::shared_ptr< Children >...children)
Definition: constraints.hh:341