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>
26 template<
typename C,
bool doIt>
27 struct ConstraintsCallBoundary
29 template<
typename F,
typename IG,
typename LFS,
typename T>
30 static void boundary (
const C& c,
const F&
f,
const IG&
ig,
const LFS& lfs, T& trafo)
34 template<
typename C,
bool doIt>
35 struct ConstraintsCallProcessor
37 template<
typename IG,
typename LFS,
typename T>
38 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
42 template<
typename C,
bool doIt>
43 struct ConstraintsCallSkeleton
45 template<
typename IG,
typename LFS,
typename T>
46 static void skeleton (
const C& c,
const IG&
ig,
47 const LFS& lfs_e,
const LFS& lfs_f,
48 T& trafo_e, T& trafo_f)
52 template<
typename C,
bool doIt>
53 struct ConstraintsCallVolume
55 template<
typename EG,
typename LFS,
typename T>
56 static void volume (
const C& c,
const EG&
eg,
const LFS& lfs, T& trafo)
63 struct ConstraintsCallBoundary<C,true>
65 template<
typename F,
typename IG,
typename LFS,
typename T>
66 static void boundary (
const C& c,
const F&
f,
const IG&
ig,
const LFS& lfs, T& trafo)
69 c.boundary(f,ig,lfs,trafo);
73 struct ConstraintsCallProcessor<C,true>
75 template<
typename IG,
typename LFS,
typename T>
76 static void processor (
const C& c,
const IG&
ig,
const LFS& lfs, T& trafo)
79 c.processor(ig,lfs,trafo);
83 struct ConstraintsCallSkeleton<C,true>
85 template<
typename IG,
typename LFS,
typename T>
86 static void skeleton (
const C& c,
const IG&
ig,
87 const LFS& lfs_e,
const LFS& lfs_f,
88 T& trafo_e, T& trafo_f)
90 if (lfs_e.size() || lfs_f.size())
91 c.skeleton(ig, lfs_e, lfs_f, trafo_e, trafo_f);
95 struct ConstraintsCallVolume<C,true>
97 template<
typename EG,
typename LFS,
typename T>
98 static void volume (
const C& c,
const EG&
eg,
const LFS& lfs, T& trafo)
101 c.volume(eg,lfs,trafo);
106 struct BoundaryConstraintsBase
107 :
public TypeTree::TreePairVisitor
113 template<
typename F,
typename LFS,
typename TreePath>
114 void leaf(
const F&
f,
const LFS& lfs, TreePath treePath)
const
116 dune_static_assert((AlwaysFalse<F>::Value),
117 "unsupported combination of function and LocalFunctionSpace");
122 template<
typename F,
typename IG,
typename CL>
123 struct BoundaryConstraintsForParametersLeaf
124 :
public TypeTree::TreeVisitor
125 ,
public TypeTree::DynamicTraversal
128 template<
typename LFS,
typename TreePath>
129 void leaf(
const LFS& lfs, TreePath treePath)
const
133 typedef typename LFS::Traits::ConstraintsType C;
136 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),
f,
ig,lfs,
cl);
139 BoundaryConstraintsForParametersLeaf(
const F& f_,
const IG& ig_, CL& cl_)
152 template<
typename IG,
typename CL>
153 struct BoundaryConstraints
154 :
public BoundaryConstraintsBase
155 ,
public TypeTree::DynamicTraversal
159 template<
typename F,
typename LFS,
typename TreePath>
160 typename enable_if<F::isLeaf && LFS::isLeaf>::type
161 leaf(
const F&
f,
const LFS& lfs, TreePath treePath)
const
164 typedef typename LFS::Traits::ConstraintsType C;
167 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),
f,
ig,lfs,
cl);
171 template<
typename F,
typename LFS,
typename TreePath>
172 typename enable_if<F::isLeaf && (!LFS::isLeaf)>::type
173 leaf(
const F& f,
const LFS& lfs, TreePath treePath)
const
176 TypeTree::applyToTree(lfs,BoundaryConstraintsForParametersLeaf<F,IG,CL>(f,ig,
cl));
179 BoundaryConstraints(
const IG& ig_, CL& cl_)
191 template<
typename IG,
typename CL>
192 struct ProcessorConstraints
193 :
public TypeTree::TreeVisitor
194 ,
public TypeTree::DynamicTraversal
197 template<
typename LFS,
typename TreePath>
198 void leaf(
const LFS& lfs, TreePath treePath)
const
201 typedef typename LFS::Traits::ConstraintsType C;
204 ConstraintsCallProcessor<C,C::doProcessor>::processor(lfs.constraints(),
ig,lfs,
cl);
207 ProcessorConstraints(
const IG& ig_, CL& cl_)
219 template<
typename IG,
typename CL>
220 struct SkeletonConstraints
221 :
public TypeTree::TreePairVisitor
222 ,
public TypeTree::DynamicTraversal
225 template<
typename LFS,
typename TreePath>
226 void leaf(
const LFS& lfs_e,
const LFS& lfs_f, TreePath treePath)
const
229 typedef typename LFS::Traits::ConstraintsType C;
234 const C & c = lfs_e.constraints();
237 ConstraintsCallSkeleton<C,C::doSkeleton>::skeleton(c,ig,lfs_e,lfs_f,
cl_e,
cl_f);
240 SkeletonConstraints(
const IG& ig_, CL& cl_e_, CL& cl_f_)
254 template<
typename EG,
typename CL>
255 struct VolumeConstraints
256 :
public TypeTree::TreeVisitor
257 ,
public TypeTree::DynamicTraversal
260 template<
typename LFS,
typename TreePath>
261 void leaf(
const LFS& lfs, TreePath treePath)
const
264 typedef typename LFS::Traits::ConstraintsType C;
265 const C & c = lfs.constraints();
268 ConstraintsCallVolume<C,C::doVolume>::volume(c,
eg,lfs,
cl);
271 VolumeConstraints(
const EG& eg_, CL& cl_)
285 template<DUNE_TYPETREE_COMPOSITENODE_TEMPLATE_CHILDREN_FOR_SPECIALIZATION>
287 :
public DUNE_TYPETREE_COMPOSITENODE_BASETYPE
289 typedef DUNE_TYPETREE_COMPOSITENODE_BASETYPE
BaseT;
292 :
BaseT(DUNE_TYPETREE_COMPOSITENODE_CHILDVARIABLES)
296 :
BaseT(DUNE_TYPETREE_COMPOSITENODE_CHILDVARIABLES)
305 template<DUNE_TYPETREE_COMPOSITENODE_TEMPLATE_CHILDREN_FOR_SPECIALIZATION>
307 :
public DUNE_TYPETREE_COMPOSITENODE_BASETYPE
309 typedef DUNE_TYPETREE_COMPOSITENODE_BASETYPE
BaseT;
312 :
BaseT(DUNE_TYPETREE_COMPOSITENODE_CHILDVARIABLES)
316 :
BaseT(DUNE_TYPETREE_COMPOSITENODE_CHILDVARIABLES)
320 template<
typename T, std::
size_t k>
322 :
public TypeTree::PowerNode<T,k>
324 typedef TypeTree::PowerNode<T,k>
BaseT;
357 :
BaseT(c0,c1,c2,c3,c4)
366 :
BaseT(c0,c1,c2,c3,c4,c5)
376 :
BaseT(c0,c1,c2,c3,c4,c5,c6)
387 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
399 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
412 :
BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
424 class OldStyleConstraintsWrapper
425 :
public TypeTree::LeafNode
427 shared_ptr<const F> _f;
431 template<
typename Transformation>
432 OldStyleConstraintsWrapper(shared_ptr<const F> f,
const Transformation& t,
unsigned int i=0)
437 template<
typename Transformation>
438 OldStyleConstraintsWrapper(
const F & f,
const Transformation& t,
unsigned int i=0)
439 : _f(stackobject_to_shared_ptr(f))
444 bool isDirichlet(
const I & intersection,
const FieldVector<typename I::ctype, I::dimension-1> & coord)
const
446 typename F::Traits::RangeType bctype;
447 _f->evaluate(intersection,coord,bctype);
448 return bctype[_i] > 0;
452 bool isNeumann(
const I & intersection,
const FieldVector<typename I::ctype, I::dimension-1> & coord)
const
454 typename F::Traits::RangeType bctype;
455 _f->evaluate(intersection,coord,bctype);
456 return bctype[_i] == 0;
461 class NoConstraintsParameters :
public TypeTree::LeafNode {};
464 struct gf_to_constraints {};
467 template<
typename F,
typename Transformation>
468 struct MultiComponentOldStyleConstraintsWrapperDescription
471 static const bool recursive =
false;
473 enum {
dim = F::Traits::dimRange };
474 typedef OldStyleConstraintsWrapper<F> node_type;
475 typedef PowerConstraintsParameters<node_type, dim> transformed_type;
476 typedef shared_ptr<transformed_type> transformed_storage_type;
478 static transformed_type transform(
const F&
s,
const Transformation& t)
480 shared_ptr<const F> sp = stackobject_to_shared_ptr(s);
481 array<shared_ptr<node_type>,
dim> childs;
482 for (
int i=0; i<
dim; i++)
483 childs[i] = make_shared<node_type>(sp,t,i);
484 return transformed_type(childs);
487 static transformed_storage_type transform_storage(shared_ptr<const F> s,
const Transformation& t)
489 array<shared_ptr<node_type>, dim> childs;
490 for (
int i=0; i<
dim; i++)
491 childs[i] = make_shared<node_type>(s,t,i);
492 return make_shared<transformed_type>(childs);
497 template<
typename Gr
idFunction>
498 typename conditional<
499 (GridFunction::Traits::dimRange == 1),
501 Dune::TypeTree::GenericLeafNodeTransformation<GridFunction,gf_to_constraints,OldStyleConstraintsWrapper<GridFunction> >,
503 MultiComponentOldStyleConstraintsWrapperDescription<GridFunction,gf_to_constraints>
508 template<
typename PowerGr
idFunction>
509 Dune::TypeTree::SimplePowerNodeTransformation<PowerGridFunction,gf_to_constraints,PowerConstraintsParameters>
513 #if HAVE_VARIADIC_TEMPLATES
514 template<
typename CompositeGr
idFunction>
515 Dune::TypeTree::SimpleVariadicCompositeNodeTransformation<CompositeGridFunction,gf_to_constraints,CompositeConstraintsParameters>
518 template<
typename CompositeGr
idFunction>
519 Dune::TypeTree::SimpleCompositeNodeTransformation<CompositeGridFunction,gf_to_constraints,CompositeConstraintsParameters>
534 template<
typename P,
typename GFS,
typename GV,
typename CG,
bool isFunction>
535 struct ConstraintsAssemblerHelper
553 assemble(
const P&
p,
const GFS& gfs,
const GV& gv, CG& cg,
const bool verbose)
556 typedef typename GV::Traits::template Codim<0>::Entity Element;
557 typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
558 typedef typename GV::IntersectionIterator IntersectionIterator;
559 typedef typename IntersectionIterator::Intersection Intersection;
562 typedef LocalFunctionSpace<GFS> LFS;
564 LFSIndexCache<LFS> lfs_cache_e(lfs_e);
566 LFSIndexCache<LFS> lfs_cache_f(lfs_f);
569 const typename GV::IndexSet& is=gv.indexSet();
572 const int chunk=1<<28;
574 std::map<Dune::GeometryType,int> gtoffset;
577 for (ElementIterator it = gv.template begin<0>();
578 it!=gv.template end<0>(); ++it)
581 if (gtoffset.find(it->type())==gtoffset.end())
583 gtoffset[it->type()] =
offset;
587 const typename GV::IndexSet::IndexType
id = is.index(*it)+gtoffset[it->type()];
592 typedef typename CG::LocalTransformation CL;
597 typedef ElementGeometry<Element> ElementWrapper;
598 TypeTree::applyToTree(lfs_e,VolumeConstraints<ElementWrapper,CL>(ElementWrapper(*it),cl_self));
601 unsigned int intersection_index = 0;
602 IntersectionIterator endit = gv.iend(*it);
603 for (IntersectionIterator iit = gv.ibegin(*it); iit!=endit; ++iit, ++intersection_index)
607 typedef IntersectionGeometry<Intersection> IntersectionWrapper;
608 TypeTree::applyToTreePair(p,lfs_e,BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(*iit,intersection_index),cl_self));
612 if ((!iit->boundary()) && (!iit->neighbor()))
614 typedef IntersectionGeometry<Intersection> IntersectionWrapper;
615 TypeTree::applyToTree(lfs_e,ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(*iit,intersection_index),cl_self));
619 if (iit->neighbor()){
621 Dune::GeometryType gtn = iit->outside()->type();
622 const typename GV::IndexSet::IndexType idn = is.index(*(iit->outside()))+gtoffset[gtn];
626 lfs_f.bind( *(iit->outside()) );
630 typedef IntersectionGeometry<Intersection> IntersectionWrapper;
631 TypeTree::applyToTreePair(lfs_e,lfs_f,SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(*iit,intersection_index),cl_self,cl_neighbor));
633 if (!cl_neighbor.empty())
635 lfs_cache_f.update();
636 cg.import_local_transformation(cl_neighbor,lfs_cache_f);
643 if (!cl_self.empty())
645 lfs_cache_e.update();
646 cg.import_local_transformation(cl_self,lfs_cache_e);
653 std::cout <<
"constraints:" << std::endl;
654 typedef typename CG::iterator global_col_iterator;
655 typedef typename CG::value_type::second_type global_row_type;
656 typedef typename global_row_type::iterator global_row_iterator;
658 std::cout << cg.size() <<
" constrained degrees of freedom" << std::endl;
660 for (global_col_iterator cit=cg.begin(); cit!=cg.end(); ++cit)
662 std::cout << cit->first <<
": ";
663 for (global_row_iterator rit=(cit->second).begin(); rit!=(cit->second).end(); ++rit)
664 std::cout <<
"(" << rit->first <<
"," << rit->second <<
") ";
665 std::cout << std::endl;
674 template<
typename F,
typename GFS,
typename GV>
675 struct ConstraintsAssemblerHelper<F, GFS, GV, EmptyTransformation, true>
677 static void assemble(
const F& f,
const GFS& gfs,
const GV& gv, EmptyTransformation& cg,
const bool verbose)
682 template<
typename F,
typename GFS,
typename GV>
683 struct ConstraintsAssemblerHelper<F, GFS, GV, EmptyTransformation, false>
685 static void assemble(
const F& f,
const GFS& gfs,
const GV& gv, EmptyTransformation& cg,
const bool verbose)
692 template<
typename F,
typename GFS,
typename GV,
typename CG>
693 struct ConstraintsAssemblerHelper<F, GFS, GV, CG, true>
696 assemble(
const F& f,
const GFS& gfs,
const GV& gv, CG& cg,
const bool verbose)
699 typedef typename Dune::TypeTree::TransformTree<F,gf_to_constraints> Transformation;
700 typedef typename Transformation::Type P;
702 P p = Transformation::transform(f);
722 template<
typename GFS,
typename CG>
724 const bool verbose =
false)
726 typedef typename GFS::Traits::GridViewType GV;
727 NoConstraintsParameters
p;
728 ConstraintsAssemblerHelper<NoConstraintsParameters, GFS, GV, CG, false>::assemble(p,gfs,gfs.gridView(),cg,verbose);
749 template<
typename P,
typename GFS,
typename CG>
751 const bool verbose =
false)
753 typedef typename GFS::Traits::GridViewType GV;
771 template<
typename CG,
typename XG>
773 typename XG::ElementType x,
776 typedef typename CG::const_iterator global_col_iterator;
777 for (global_col_iterator cit=cg.begin(); cit!=cg.end(); ++cit)
785 template<
typename XG>
787 typename XG::ElementType x,
815 template<
typename CG,
typename XG,
typename Cmp>
817 XG&
xg,
const Cmp& cmp = Cmp())
819 typedef typename CG::const_iterator global_col_iterator;
820 for (global_col_iterator cit=cg.begin(); cit!=cg.end(); ++cit)
821 if(cmp.ne(xg[cit->first], x))
845 template<
typename CG,
typename XG>
850 FloatCmpOps<typename XG::ElementType>());
857 template<
typename XG,
typename Cmp>
859 XG& xg,
const Cmp& cmp = Cmp())
865 template<
typename XG>
881 template<
typename CG,
typename XG>
884 typedef typename CG::const_iterator global_col_iterator;
885 typedef typename CG::value_type::second_type::const_iterator global_row_iterator;
887 for (global_col_iterator cit=cg.begin(); cit!=cg.end(); ++cit)
888 for(global_row_iterator rit = cit->second.begin(); rit!=cit->second.end(); ++rit)
889 xg[rit->first] += rit->second * xg[cit->first];
893 for (global_col_iterator cit=cg.begin(); cit!=cg.end(); ++cit)
894 xg[cit->first] =
typename XG::ElementType(0);
901 template<
typename XG>
916 template<
typename CG,
typename XG>
919 typedef typename CG::const_iterator global_col_iterator;
920 for (global_col_iterator cit=cg.begin(); cit!=cg.end(); ++cit)
922 xgout[cit->first] = xgin[cit->first];
930 template<
typename XG>
942 template<
typename CG,
typename XG>
961 template<
typename XG>
975 template<
typename CG,
typename XG>
994 template<
typename XG>
1008 template<
typename CG,
typename XG>
1016 typedef typename CG::const_iterator global_col_iterator;
1017 for (global_col_iterator cit=cg.begin(); cit!=cg.end(); ++cit)
1018 if (cit->second.size() == 0)
1020 tmp[cit->first] = xg[cit->first];
1040 template<
typename XG>
1041 void set_shifted_dofs (
const EmptyTransformation& cg,
typename XG::ElementType x, XG& xg)
1052 #endif // DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4)
Definition: common/constraints.hh:352
PowerConstraintsParameters(T &c)
Definition: common/constraints.hh:330
Definition: common/constraints.hh:286
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8, T &c9)
Definition: common/constraints.hh:402
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7)
Definition: common/constraints.hh:379
CL & cl
Definition: common/constraints.hh:147
CompositeConstraintsOperator(DUNE_TYPETREE_COMPOSITENODE_CONSTRUCTOR_SIGNATURE)
Definition: common/constraints.hh:291
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: common/constraints.hh:1009
CL & cl_f
Definition: common/constraints.hh:249
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6)
Definition: common/constraints.hh:369
PowerConstraintsParameters(T &c0, T &c1)
Definition: common/constraints.hh:334
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: common/constraints.hh:772
void constrain_residual(const CG &cg, XG &xg)
transform residual into transformed basis: r -> r~
Definition: common/constraints.hh:882
CompositeConstraintsOperator(DUNE_TYPETREE_COMPOSITENODE_STORAGE_CONSTRUCTOR_SIGNATURE)
Definition: common/constraints.hh:295
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: common/constraints.hh:816
void copy_nonconstrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: common/constraints.hh:976
static const int dim
Definition: adaptivity.hh:82
const std::size_t offset
Definition: localfunctionspace.hh:74
TypeTree::PowerNode< T, k > BaseT
Definition: common/constraints.hh:324
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3)
Definition: common/constraints.hh:345
XG & xg
Definition: interpolate.hh:67
PowerConstraintsParameters(T &c0, T &c1, T &c2)
Definition: common/constraints.hh:339
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: common/constraints.hh:723
CompositeConstraintsParameters(DUNE_TYPETREE_COMPOSITENODE_STORAGE_CONSTRUCTOR_SIGNATURE)
Definition: common/constraints.hh:315
const IG & ig
Definition: common/constraints.hh:146
DUNE_TYPETREE_COMPOSITENODE_BASETYPE BaseT
Definition: common/constraints.hh:289
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
Dune::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)
const F & f
Definition: common/constraints.hh:145
CompositeConstraintsParameters(DUNE_TYPETREE_COMPOSITENODE_CONSTRUCTOR_SIGNATURE)
Definition: common/constraints.hh:311
CL & cl_e
Definition: common/constraints.hh:248
Definition: common/constraints.hh:321
void copy_constrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: common/constraints.hh:917
const EG & eg
Definition: common/constraints.hh:277
DUNE_TYPETREE_COMPOSITENODE_BASETYPE BaseT
Definition: common/constraints.hh:309
Definition: common/constraints.hh:306
PowerConstraintsParameters()
Definition: common/constraints.hh:326
R p(int i, D x)
Definition: qkdg.hh:62
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5, T &c6, T &c7, T &c8)
Definition: common/constraints.hh:390
PowerConstraintsParameters(const array< shared_ptr< T >, k > &children)
Definition: common/constraints.hh:415
PowerConstraintsParameters(T &c0, T &c1, T &c2, T &c3, T &c4, T &c5)
Definition: common/constraints.hh:360
void set_nonconstrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: common/constraints.hh:943
const std::string s
Definition: function.hh:1103