DUNE PDELab (2.7)

constraints.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4#ifndef DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
5#define DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
6
7#include<dune/common/exceptions.hh>
8#include<dune/common/shared_ptr.hh>
9#include<dune/common/float_cmp.hh>
10
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"
18
19namespace Dune {
20 namespace PDELab {
21
25
26 namespace { // hide internals
27 // do method invocation only if class has the method
28
29 template<typename C, bool doIt>
30 struct ConstraintsCallBoundary
31 {
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)
34 {
35 }
36 };
37 template<typename C, bool doIt>
38 struct ConstraintsCallProcessor
39 {
40 template<typename IG, typename LFS, typename T>
41 static void processor (const C& c, const IG& ig, const LFS& lfs, T& trafo)
42 {
43 }
44 };
45 template<typename C, bool doIt>
46 struct ConstraintsCallSkeleton
47 {
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)
52 {
53 }
54 };
55 template<typename C, bool doIt>
56 struct ConstraintsCallVolume
57 {
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)
60 {
61 }
62 };
63
64
65 template<typename C>
66 struct ConstraintsCallBoundary<C,true>
67 {
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)
70 {
71 if (lfs.size())
72 c.boundary(p,ig,lfs,trafo);
73 }
74 };
75 template<typename C>
76 struct ConstraintsCallProcessor<C,true>
77 {
78 template<typename IG, typename LFS, typename T>
79 static void processor (const C& c, const IG& ig, const LFS& lfs, T& trafo)
80 {
81 if (lfs.size())
82 c.processor(ig,lfs,trafo);
83 }
84 };
85 template<typename C>
86 struct ConstraintsCallSkeleton<C,true>
87 {
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)
92 {
93 if (lfs_e.size() || lfs_f.size())
94 c.skeleton(ig, lfs_e, lfs_f, trafo_e, trafo_f);
95 }
96 };
97 template<typename C>
98 struct ConstraintsCallVolume<C,true>
99 {
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)
102 {
103 if (lfs.size())
104 c.volume(p,eg,lfs,trafo);
105 }
106 };
107
108
109 struct ParameterizedConstraintsBase
110 : public TypeTree::TreePairVisitor
111 {
112 // This acts as a catch-all for unsupported leaf- / non-leaf combinations in the two
113 // trees. It is necessary because otherwise, the visitor would fall back to the default
114 // implementation in TreeVisitor, which simply does nothing. The resulting bugs would
115 // probably be hell to find...
116 template<typename P, typename LFS, typename TreePath>
117 void leaf(const P& p, const LFS& lfs, TreePath treePath) const
118 {
119 static_assert((AlwaysFalse<P>::Value),
120 "unsupported combination of function and LocalFunctionSpace");
121 }
122 };
123
124
125 template<typename P, typename IG, typename CL>
126 struct BoundaryConstraintsForParametersLeaf
127 : public TypeTree::TreeVisitor
128 , public TypeTree::DynamicTraversal
129 {
130
131 template<typename LFS, typename TreePath>
132 void leaf(const LFS& lfs, TreePath treePath) const
133 {
134
135 // extract constraints type
136 typedef typename LFS::Traits::ConstraintsType C;
137
138 // iterate over boundary, need intersection iterator
139 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),p,ig,lfs,cl);
140 }
141
142 BoundaryConstraintsForParametersLeaf(const P& p_, const IG& ig_, CL& cl_)
143 : p(p_)
144 , ig(ig_)
145 , cl(cl_)
146 {}
147
148 const P& p;
149 const IG& ig;
150 CL& cl;
151
152 };
153
154
155 template<typename IG, typename CL>
156 struct BoundaryConstraints
157 : public ParameterizedConstraintsBase
158 , public TypeTree::DynamicTraversal
159 {
160
161 // standard case - leaf in both trees
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
165 {
166 // extract constraints type
167 typedef typename LFS::Traits::ConstraintsType C;
168
169 // iterate over boundary, need intersection iterator
170 ConstraintsCallBoundary<C,C::doBoundary>::boundary(lfs.constraints(),p,ig,lfs,cl);
171 }
172
173 // reuse constraints parameter information from p for all LFS children
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
177 {
178 // traverse LFS tree and reuse parameter information
179 TypeTree::applyToTree(lfs,BoundaryConstraintsForParametersLeaf<P,IG,CL>(p,ig,cl));
180 }
181
182 BoundaryConstraints(const IG& ig_, CL& cl_)
183 : ig(ig_)
184 , cl(cl_)
185 {}
186
187 private:
188 const IG& ig;
189 CL& cl;
190
191 };
192
193
194 template<typename IG, typename CL>
195 struct ProcessorConstraints
196 : public TypeTree::TreeVisitor
197 , public TypeTree::DynamicTraversal
198 {
199
200 template<typename LFS, typename TreePath>
201 void leaf(const LFS& lfs, TreePath treePath) const
202 {
203 // extract constraints type
204 typedef typename LFS::Traits::ConstraintsType C;
205
206 // iterate over boundary, need intersection iterator
207 ConstraintsCallProcessor<C,C::doProcessor>::processor(lfs.constraints(),ig,lfs,cl);
208 }
209
210 ProcessorConstraints(const IG& ig_, CL& cl_)
211 : ig(ig_)
212 , cl(cl_)
213 {}
214
215 private:
216 const IG& ig;
217 CL& cl;
218
219 };
220
221
222 template<typename IG, typename CL>
223 struct SkeletonConstraints
224 : public TypeTree::TreePairVisitor
225 , public TypeTree::DynamicTraversal
226 {
227
228 template<typename LFS, typename TreePath>
229 void leaf(const LFS& lfs_e, const LFS& lfs_f, TreePath treePath) const
230 {
231 // extract constraints type
232 typedef typename LFS::Traits::ConstraintsType C;
233
234 // as LFS::constraints() just returns the constraints of the
235 // GridFunctionSpace, lfs_e.constraints() is equivalent to
236 // lfs_f.constraints()
237 const C & c = lfs_e.constraints();
238
239 // iterate over boundary, need intersection iterator
240 ConstraintsCallSkeleton<C,C::doSkeleton>::skeleton(c,ig,lfs_e,lfs_f,cl_e,cl_f);
241 }
242
243 SkeletonConstraints(const IG& ig_, CL& cl_e_, CL& cl_f_)
244 : ig(ig_)
245 , cl_e(cl_e_)
246 , cl_f(cl_f_)
247 {}
248
249 private:
250 const IG& ig;
251 CL& cl_e;
252 CL& cl_f;
253
254 };
255
256
257 template<typename P, typename EG, typename CL>
258 struct VolumeConstraintsForParametersLeaf
259 : public TypeTree::TreeVisitor
260 , public TypeTree::DynamicTraversal
261 {
262
263 template<typename LFS, typename TreePath>
264 void leaf(const LFS& lfs, TreePath treePath) const
265 {
266 // extract constraints type
267 typedef typename LFS::Traits::ConstraintsType C;
268 const C & c = lfs.constraints();
269
270 // iterate over boundary, need intersection iterator
271 ConstraintsCallVolume<C,C::doVolume>::volume(c,p,eg,lfs,cl);
272 }
273
274 VolumeConstraintsForParametersLeaf(const P& p_, const EG& eg_, CL& cl_)
275 : p(p_)
276 , eg(eg_)
277 , cl(cl_)
278 {}
279
280 private:
281 const P& p;
282 const EG& eg;
283 CL& cl;
284
285 };
286
287
288 template<typename EG, typename CL>
289 struct VolumeConstraints
290 : public ParameterizedConstraintsBase
291 , public TypeTree::DynamicTraversal
292 {
293
294 // standard case - leaf in both trees
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
298 {
299 // extract constraints type
300 typedef typename LFS::Traits::ConstraintsType C;
301 const C & c = lfs.constraints();
302
303 // iterate over boundary, need intersection iterator
304 ConstraintsCallVolume<C,C::doVolume>::volume(c,p,eg,lfs,cl);
305 }
306
307 // reuse constraints parameter information from p for all LFS children
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
311 {
312 // traverse LFS tree and reuse parameter information
313 TypeTree::applyToTree(lfs,VolumeConstraintsForParametersLeaf<P,EG,CL>(p,eg,cl));
314 }
315
316 VolumeConstraints(const EG& eg_, CL& cl_)
317 : eg(eg_)
318 , cl(cl_)
319 {}
320
321 private:
322 const EG& eg;
323 CL& cl;
324
325 };
326
327
328 } // anonymous namespace
329
330 template<typename... Children>
331 struct CompositeConstraintsOperator
332 : public TypeTree::CompositeNode<Children...>
333 {
334 typedef TypeTree::CompositeNode<Children...> BaseT;
335
336 CompositeConstraintsOperator(Children&... children)
337 : BaseT(children...)
338 {}
339
340 CompositeConstraintsOperator(std::shared_ptr<Children>... children)
341 : BaseT(children...)
342 {}
343
344 // aggregate flags
345
346 // forward methods to childs
347
348 };
349
350 template<typename... Children>
351 struct CompositeConstraintsParameters
352 : public TypeTree::CompositeNode<Children...>
353 {
354 typedef TypeTree::CompositeNode<Children...> BaseT;
355
356 CompositeConstraintsParameters(Children&... children)
357 : BaseT(children...)
358 {}
359
360 CompositeConstraintsParameters(std::shared_ptr<Children>... children)
361 : BaseT(children...)
362 {}
363 };
364
365 template<typename T, std::size_t k>
366 struct PowerConstraintsParameters
367 : public TypeTree::PowerNode<T,k>
368 {
369 typedef TypeTree::PowerNode<T,k> BaseT;
370
371 PowerConstraintsParameters()
372 : BaseT()
373 {}
374
375 PowerConstraintsParameters(T& c)
376 : BaseT(c)
377 {}
378
379 PowerConstraintsParameters (T& c0,
380 T& c1)
381 : BaseT(c0,c1)
382 {}
383
384 PowerConstraintsParameters (T& c0,
385 T& c1,
386 T& c2)
387 : BaseT(c0,c1,c2)
388 {}
389
390 PowerConstraintsParameters (T& c0,
391 T& c1,
392 T& c2,
393 T& c3)
394 : BaseT(c0,c1,c2,c3)
395 {}
396
397 PowerConstraintsParameters (T& c0,
398 T& c1,
399 T& c2,
400 T& c3,
401 T& c4)
402 : BaseT(c0,c1,c2,c3,c4)
403 {}
404
405 PowerConstraintsParameters (T& c0,
406 T& c1,
407 T& c2,
408 T& c3,
409 T& c4,
410 T& c5)
411 : BaseT(c0,c1,c2,c3,c4,c5)
412 {}
413
414 PowerConstraintsParameters (T& c0,
415 T& c1,
416 T& c2,
417 T& c3,
418 T& c4,
419 T& c5,
420 T& c6)
421 : BaseT(c0,c1,c2,c3,c4,c5,c6)
422 {}
423
424 PowerConstraintsParameters (T& c0,
425 T& c1,
426 T& c2,
427 T& c3,
428 T& c4,
429 T& c5,
430 T& c6,
431 T& c7)
432 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
433 {}
434
435 PowerConstraintsParameters (T& c0,
436 T& c1,
437 T& c2,
438 T& c3,
439 T& c4,
440 T& c5,
441 T& c6,
442 T& c7,
443 T& c8)
444 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
445 {}
446
447 PowerConstraintsParameters (T& c0,
448 T& c1,
449 T& c2,
450 T& c3,
451 T& c4,
452 T& c5,
453 T& c6,
454 T& c7,
455 T& c8,
456 T& c9)
457 : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
458 {}
459
460 PowerConstraintsParameters (const std::array<std::shared_ptr<T>,k>& children)
461 : BaseT(children)
462 {}
463 };
464
465#ifndef DOXYGEN
466
468 template<typename F>
469 class OldStyleConstraintsWrapper
470 : public TypeTree::LeafNode
471 {
472 std::shared_ptr<const F> _f;
473 unsigned int _i;
474 public:
475
476 template<typename Transformation>
477 OldStyleConstraintsWrapper(std::shared_ptr<const F> f, const Transformation& t, unsigned int i=0)
478 : _f(f)
479 , _i(i)
480 {}
481
482 template<typename Transformation>
483 OldStyleConstraintsWrapper(const F & f, const Transformation& t, unsigned int i=0)
484 : _f(stackobject_to_shared_ptr(f))
485 , _i(i)
486 {}
487
488 template<typename I>
489 bool isDirichlet(const I & intersection, const FieldVector<typename I::ctype, I::mydimension> & coord) const
490 {
491 typename F::Traits::RangeType bctype;
492 _f->evaluate(intersection,coord,bctype);
493 return bctype[_i] > 0;
494 }
495
496 template<typename I>
497 bool isNeumann(const I & intersection, const FieldVector<typename I::ctype, I::mydimension> & coord) const
498 {
499 typename F::Traits::RangeType bctype;
500 _f->evaluate(intersection,coord,bctype);
501 return bctype[_i] == 0;
502 }
503 };
504
505 // Tag to name trafo GridFunction -> OldStyleConstraintsWrapper
506 struct gf_to_constraints {};
507
508 // register trafos GridFunction -> OldStyleConstraintsWrapper
509 template<typename F, typename Transformation>
510 struct MultiComponentOldStyleConstraintsWrapperDescription
511 {
512
513 static const bool recursive = false;
514
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;
519
520 static transformed_type transform(const F& s, const Transformation& t)
521 {
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);
527 }
528
529 static transformed_storage_type transform_storage(std::shared_ptr<const F> s, const Transformation& t)
530 {
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);
535 }
536
537 };
538 // trafos for leaf nodes
539 template<typename GridFunction>
540 typename std::conditional<
541 (GridFunction::Traits::dimRange == 1),
542 // trafo for scalar leaf nodes
543 TypeTree::GenericLeafNodeTransformation<GridFunction,gf_to_constraints,OldStyleConstraintsWrapper<GridFunction> >,
544 // trafo for multi component leaf nodes
545 MultiComponentOldStyleConstraintsWrapperDescription<GridFunction,gf_to_constraints>
546 >::type
547 registerNodeTransformation(GridFunction*, gf_to_constraints*, GridFunctionTag*);
548
549 // trafo for power nodes
550 template<typename PowerGridFunction>
551 TypeTree::SimplePowerNodeTransformation<PowerGridFunction,gf_to_constraints,PowerConstraintsParameters>
552 registerNodeTransformation(PowerGridFunction*, gf_to_constraints*, PowerGridFunctionTag*);
553
554 // trafos for composite nodes
555 template<typename CompositeGridFunction>
556 TypeTree::SimpleCompositeNodeTransformation<CompositeGridFunction,gf_to_constraints,CompositeConstraintsParameters>
557 registerNodeTransformation(CompositeGridFunction*, gf_to_constraints*, CompositeGridFunctionTag*);
558
560
570 template<typename P, typename GFS, typename CG, bool isFunction>
571 struct ConstraintsAssemblerHelper
572 {
574
588 static void
589 assemble(const P& p, const GFS& gfs, CG& cg, const bool verbose)
590 {
591 // get some types
592 using ES = typename GFS::Traits::EntitySet;
593 using Element = typename ES::Traits::Element;
594 using Intersection = typename ES::Traits::Intersection;
595
596 ES es = gfs.entitySet();
597
598 // make local function space
599 using LFS = LocalFunctionSpace<GFS>;
600 LFS lfs_e(gfs);
601 LFSIndexCache<LFS> lfs_cache_e(lfs_e);
602 LFS lfs_f(gfs);
603 LFSIndexCache<LFS> lfs_cache_f(lfs_f);
604
605 // get index set
606 auto& is = es.indexSet();
607
608 // loop once over the grid
609 for (const auto& element : elements(es))
610 {
611
612 auto id = is.uniqueIndex(element);
613
614 // bind local function space to element
615 lfs_e.bind(element);
616
617 using CL = typename CG::LocalTransformation;
618
619 CL cl_self;
620
621 using ElementWrapper = ElementGeometry<Element>;
622 using IntersectionWrapper = IntersectionGeometry<Intersection>;
623
624 TypeTree::applyToTreePair(p,lfs_e,VolumeConstraints<ElementWrapper,CL>(ElementWrapper(element),cl_self));
625
626 // iterate over intersections and call metaprogram
627 unsigned int intersection_index = 0;
628 for (const auto& intersection : intersections(es,element))
629 {
630
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);
634
635 switch (intersection_type) {
636
637 case IntersectionType::skeleton:
638 case IntersectionType::periodic:
639 {
640 auto idn = is.uniqueIndex(outside_element);
641
642 if(id > idn){
643 // bind local function space to element in neighbor
644 lfs_f.bind(outside_element);
645
646 CL cl_neighbor;
647
648 TypeTree::applyToTreePair(lfs_e,lfs_f,SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
649
650 if (!cl_neighbor.empty())
651 {
652 lfs_cache_f.update();
653 cg.import_local_transformation(cl_neighbor,lfs_cache_f);
654 }
655
656 }
657 break;
658 }
659
660 case IntersectionType::boundary:
661 TypeTree::applyToTreePair(p,lfs_e,BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
662 break;
663
664 case IntersectionType::processor:
665 TypeTree::applyToTree(lfs_e,ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
666 break;
667
668 }
669 ++intersection_index;
670 }
671
672 if (!cl_self.empty())
673 {
674 lfs_cache_e.update();
675 cg.import_local_transformation(cl_self,lfs_cache_e);
676 }
677
678 }
679
680 // print result
681 if(verbose){
682 std::cout << "constraints:" << std::endl;
683
684 std::cout << cg.size() << " constrained degrees of freedom" << std::endl;
685
686 for (const auto& col : cg)
687 {
688 std::cout << col.first << ": "; // col index
689 for (const auto& row : col.second)
690 std::cout << "(" << row.first << "," << row.second << ") "; // row index , value
691 std::cout << std::endl;
692 }
693 }
694 }
695 }; // end ConstraintsAssemblerHelper
696
697
698
699 // Disable constraints assembly for empty transformation
700 template<typename F, typename GFS>
701 struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, true>
702 {
703 static void assemble(const F& f, const GFS& gfs, EmptyTransformation& cg, const bool verbose)
704 {}
705 };
706
707 // Disable constraints assembly for empty transformation
708 template<typename F, typename GFS>
709 struct ConstraintsAssemblerHelper<F, GFS, EmptyTransformation, false>
710 {
711 static void assemble(const F& f, const GFS& gfs, EmptyTransformation& cg, const bool verbose)
712 {}
713 };
714
715
716
717 // Backwards compatibility shim
718 template<typename F, typename GFS, typename CG>
719 struct ConstraintsAssemblerHelper<F, GFS, CG, true>
720 {
721 static void
722 assemble(const F& f, const GFS& gfs, CG& cg, const bool verbose)
723 {
724 // type of transformed tree
725 typedef typename TypeTree::TransformTree<F,gf_to_constraints> Transformation;
726 typedef typename Transformation::Type P;
727 // transform tree
728 P p = Transformation::transform(f);
729 // call parameter based implementation
730 ConstraintsAssemblerHelper<P, GFS, CG, IsGridFunction<P>::value>::assemble(p,gfs,cg,verbose);
731 }
732 };
733#endif
734
736
748 template<typename GFS, typename CG>
749 void constraints(const GFS& gfs, CG& cg,
750 const bool verbose = false)
751 {
752 NoConstraintsParameters p;
753 ConstraintsAssemblerHelper<NoConstraintsParameters, GFS, CG, false>::assemble(p,gfs,cg,verbose);
754 }
755
757
774 template<typename P, typename GFS, typename CG>
775 void constraints(const P& p, const GFS& gfs, CG& cg,
776 const bool verbose = false)
777 {
778 // clear global constraints
779 cg.clear();
780 ConstraintsAssemblerHelper<P, GFS, CG, IsGridFunction<P>::value>::assemble(p,gfs,cg,verbose);
781 }
782
784
795 template<typename CG, typename XG>
796 void set_constrained_dofs(const CG& cg,
797 typename XG::ElementType x,
798 XG& xg)
799 {
800 for (const auto& col : cg)
801 xg[col.first] = x;
802 }
803
804
805#ifndef DOXYGEN
806
807 // Specialized version for unconstrained spaces
808 template<typename XG>
809 void set_constrained_dofs(const EmptyTransformation& cg,
810 typename XG::ElementType x,
811 XG& xg)
812 {}
813
814#endif // DOXYGEN
815
816
818
838 template<typename CG, typename XG, typename Cmp>
839 bool check_constrained_dofs(const CG& cg, typename XG::ElementType x,
840 XG& xg, const Cmp& cmp = Cmp())
841 {
842 for (const auto& col : cg)
843 if(cmp.ne(xg[col.first], x))
844 return false;
845 return true;
846 }
847
849
867 template<typename CG, typename XG>
868 bool check_constrained_dofs(const CG& cg, typename XG::ElementType x,
869 XG& xg)
870 {
871 return check_constrained_dofs(cg, x, xg,
872 FloatCmpOps<typename XG::ElementType>());
873 }
874
875
876#ifndef DOXYGEN
877
878 // Specialized version for unconstrained spaces
879 template<typename XG, typename Cmp>
880 bool check_constrained_dofs(const EmptyTransformation& cg, typename XG::ElementType x,
881 XG& xg, const Cmp& cmp = Cmp())
882 {
883 return true;
884 }
885
886 // Specialized version for unconstrained spaces
887 template<typename XG>
888 bool check_constrained_dofs(const EmptyTransformation& cg, typename XG::ElementType x,
889 XG& xg)
890 {
891 return true;
892 }
893
894#endif // DOXYGEN
895
896
898
903 template<typename CG, typename XG>
904 void constrain_residual (const CG& cg, XG& xg)
905 {
906 for (const auto& col : cg)
907 for (const auto& row : col.second)
908 xg[row.first] += row.second * xg[col.first];
909
910 // extra loop because constrained dofs might have contributions
911 // to constrained dofs
912 for (const auto& col : cg)
913 xg[col.first] = typename XG::ElementType(0);
914 }
915
916
917#ifndef DOXYGEN
918
919 // Specialized version for unconstrained spaces
920 template<typename XG>
921 void constrain_residual (const EmptyTransformation& cg, XG& xg)
922 {}
923
924#endif // DOXYGEN
925
929
935 template<typename CG, typename XG>
936 void copy_constrained_dofs (const CG& cg, const XG& xgin, XG& xgout)
937 {
938 for (const auto& col : cg)
939 xgout[col.first] = xgin[col.first];
940 }
941
942
943#ifndef DOXYGEN
944
945 // Specialized version for unconstrained spaces
946 template<typename XG>
947 void copy_constrained_dofs (const EmptyTransformation& cg, const XG& xgin, XG& xgout)
948 {}
949
950#endif // DOXYGEN
951
952
959 template<typename CG, typename XG>
960 void set_nonconstrained_dofs (const CG& cg, typename XG::ElementType x, XG& xg)
961 {
962 XG tmp(xg);
963 xg = x;
964 copy_constrained_dofs(cg,tmp,xg);
965 }
966
967
968#ifndef DOXYGEN
969
970 // Specialized version for unconstrained spaces
971 template<typename XG>
972 void set_nonconstrained_dofs (const EmptyTransformation& cg, typename XG::ElementType x, XG& xg)
973 {
974 xg = x;
975 }
976
977#endif // DOXYGEN
978
979
986 template<typename CG, typename XG>
987 void copy_nonconstrained_dofs (const CG& cg, const XG& xgin, XG& xgout)
988 {
989 XG tmp(xgin);
990 copy_constrained_dofs(cg,xgout,tmp);
991 xgout = tmp;
992 }
993
994
995#ifndef DOXYGEN
996
997 // Specialized version for unconstrained spaces
998 template<typename XG>
999 void copy_nonconstrained_dofs (const EmptyTransformation& cg, const XG& xgin, XG& xgout)
1000 {
1001 xgout = xgin;
1002 }
1003
1004#endif // DOXYGEN
1005
1006
1013 template<typename CG, typename XG>
1014 void set_shifted_dofs (const CG& cg, typename XG::ElementType x, XG& xg)
1015 {
1016 XG tmp(xg);
1017 tmp = x;
1018
1019 for (const auto& col : cg)
1020 {
1021 // this is our marker for Dirichlet constraints
1022 if (col.second.size() == 0)
1023 {
1024 tmp[col.first] = xg[col.first];
1025 }
1026 }
1027
1028 xg = tmp;
1029 }
1030
1031
1032#ifndef DOXYGEN
1033
1034 // Specialized version for unconstrained spaces
1035 template<typename XG>
1036 void set_shifted_dofs (const EmptyTransformation& cg, typename XG::ElementType x, XG& xg)
1037 {}
1038
1039#endif // DOXYGEN
1040
1042
1044 } // namespace PDELab
1045} // namespace Dune
1046
1047#endif // DUNE_PDELAB_CONSTRAINTS_COMMON_CONSTRAINTS_HH
Dune::FloatCmpOps
Class encapsulating a default epsilon.
Definition: float_cmp.hh:286
Dune::PDELab::NoConstraintsParameters
Definition: constraintsparameters.hh:293
exceptions.hh
A few common exception classes.
float_cmp.hh
Various ways to compare floating-point numbers.
Dune::PDELab::check_constrained_dofs
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
Dune::PDELab::set_nonconstrained_dofs
void set_nonconstrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:960
Dune::PDELab::constraints
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
Dune::PDELab::set_constrained_dofs
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
Dune::PDELab::constrain_residual
void constrain_residual(const CG &cg, XG &xg)
transform residual into transformed basis: r -> r~
Definition: constraints.hh:904
Dune::PDELab::set_shifted_dofs
void set_shifted_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
Definition: constraints.hh:1014
Dune::PDELab::copy_constrained_dofs
void copy_constrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:936
Dune::PDELab::copy_nonconstrained_dofs
void copy_nonconstrained_dofs(const CG &cg, const XG &xgin, XG &xgout)
Definition: constraints.hh:987
Dune::TypeTree::registerNodeTransformation
void registerNodeTransformation(SourceNode *, Transformation *, Tag *)
Register transformation descriptor to transform SourceNode with Transformation.
Dune::TypeTree::treePath
constexpr HybridTreePath< T... > treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:188
Dune::TypeTree::applyToTreePair
void applyToTreePair(Tree1 &&tree1, Tree2 &&tree2, Visitor &&visitor)
Apply visitor to a pair of TypeTrees.
Definition: pairtraversal.hh:107
Dune::TypeTree::applyToTree
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:213
Dune
Dune namespace.
Definition: alignedallocator.hh:14
Dune::stackobject_to_shared_ptr
shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:75
shared_ptr.hh
This file implements the class shared_ptr (a reference counting pointer), for those systems that don'...
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)