dune-geometry 3.0-git
simplex.cc
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
4#define DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
5
6// This file is part of DUNE, a Distributed and Unified Numerics Environment
7// This file is copyright (C) 2005 Jorrit Fahlke <jorrit@jorrit.de>
8// This file is licensed under version 2 of the GNU General Public License,
9// with a special "runtime exception." See COPYING at the top of the source
10// tree for the full licence.
11
249#include <algorithm>
250
251#include <dune/common/fvector.hh>
252
256
257#include "base.cc"
258
259namespace Dune {
260
261 namespace RefinementImp {
262
269 namespace Simplex {
270
271 // //////////////////
272 //
274 //
275
277
282 inline int factorial(int n)
283 {
284 int prod = 1;
285 for(int i = 1; i <= n; ++i)
286 prod *= i;
287 return prod;
288 }
289
294 inline int binomial(int upper, int lower)
295 {
296 lower = std::min( lower, upper - lower );
297 if(lower < 0)
298 return 0;
299 int prod = 1;
300 for(int i = upper - lower; i < upper; ++i)
301 prod *= (i+1);
302 return prod / factorial(lower);
303 }
304
311 template<int dimension>
312 int pointIndex(const FieldVector<int, dimension> &point)
313 {
314 int index = 0;
315 for(int i = 0; i < dimension; ++i)
316 index += binomial(dimension-i + point[i]-1, dimension-i);
317 return index;
318 }
319
324 template<int n>
325 FieldVector<int, n> getPermutation(int m)
326 {
327 FieldVector<int, n> perm;
328 for(int i = 0; i < n; ++i)
329 perm[i] = i;
330
331 int base = 1;
332 for(int i = 1; i <= n; ++i)
333 base *= i;
334
335 for(int i = n; i > 0; --i) {
336 base /= i;
337 int d = m / base;
338 m %= base;
339 int t = perm[i-1]; perm[i-1] = perm[i-1-d]; perm[i-1-d] = t;
340 }
341 return perm;
342 }
343
344#if 0
345 Has to be checked
346 // calculate the index of a permutation
347 template<int n>
348 int getPermIndex(const FieldVector<int, n>& test) // O(n^2)
349 {
350 int m = 0;
351 FieldVector<int, n> perm;
352 for(int i = 0; i < n; ++i)
353 perm[i] = i;
354
355 int base = 1;
356 for(int i = 1; i <= n; ++i)
357 base *= i;
358
359 for(int i = n; i > 0; --i) {
360 base /= i;
361 int d;
362 for(d = 0; d < i; ++d)
363 if(test[i-1] == perm[i-1-d])
364 break;
365 m += d * base;
366 int d = m / base;
367 m %= base;
368 perm[i-1-d] = perm[i-1];
369 }
370 }
371#endif
372
373 // map between the reference simplex and some arbitrary kuhn simplex (denoted by it's permutation)
381 template<int dimension, class CoordType>
382 FieldVector<CoordType, dimension>
384 FieldVector<CoordType, dimension> point,
386 const FieldVector<int, dimension> &kuhn)
387 {
388 for(int i = dimension - 1; i > 0; --i)
389 point[kuhn[i-1]] += point[kuhn[i]];
390 return point;
391 }
392
400 template<int dimension, class CoordType>
401 FieldVector<CoordType, dimension>
403 FieldVector<CoordType, dimension> point,
405 const FieldVector<int, dimension> &kuhn)
406 {
407 for(int i = 0; i < dimension - 1; ++i)
408 point[kuhn[i]] -= point[kuhn[i+1]];
409 return point;
410 }
411
412
414
415 // /////////////////////////////////////////
416 //
417 // refinement implementation for simplices
418 //
419
420 template<int dimension_, class CoordType>
422 {
423 public:
424 enum { dimension = dimension_ };
425 typedef CoordType ctype;
426
427 template<int codimension>
428 struct Codim;
430 typedef FieldVector<CoordType, dimension> CoordVector;
432 typedef FieldVector<int, dimension+1> IndexVector;
433
434 static int nVertices(int level);
435 static VertexIterator vBegin(int level);
436 static VertexIterator vEnd(int level);
437
438 static int nElements(int level);
439 static ElementIterator eBegin(int level);
440 static ElementIterator eEnd(int level);
441 };
442
443 template<int dimension, class CoordType>
444 template<int codimension>
445 struct RefinementImp<dimension, CoordType>::Codim
446 {
447 class SubEntityIterator;
448 // We don't need the caching, but the uncached MultiLinearGeometry has bug FS#1209
450 };
451
452 template<int dimension, class CoordType>
453 int
455 nVertices(int level)
456 {
457 return binomial(dimension + (1 << level), dimension);
458 }
459
460 template<int dimension, class CoordType>
463 vBegin(int level)
464 {
465 return VertexIterator(level);
466 }
467
468 template<int dimension, class CoordType>
471 vEnd(int level)
472 {
473 return VertexIterator(level, true);
474 }
475
476 template<int dimension, class CoordType>
477 int
479 nElements(int level)
480 {
481 return 1 << (level * dimension);
482 }
483
484 template<int dimension, class CoordType>
487 eBegin(int level)
488 {
489 return ElementIterator(level);
490 }
491
492 template<int dimension, class CoordType>
495 eEnd(int level)
496 {
497 return ElementIterator(level, true);
498 }
499
500 // //////////////
501 //
502 // The iterator
503 //
504
505 template<int dimension, class CoordType, int codimension>
507
508 // vertices
509
510 template<int dimension, class CoordType>
511 class RefinementIteratorSpecial<dimension, CoordType, dimension>
512 {
513 public:
515 typedef typename Refinement::CoordVector CoordVector;
516 typedef typename Refinement::template Codim<dimension>::Geometry Geometry;
518
519 RefinementIteratorSpecial(int level, bool end = false);
520
521 void increment();
522 bool equals(const This &other) const;
523
524 CoordVector coords() const;
525 Geometry geometry () const;
526
527 int index() const;
528 protected:
529 typedef FieldVector<int, dimension> Vertex;
530
531 int size;
533 };
534
535 template<int dimension, class CoordType>
537 RefinementIteratorSpecial(int level, bool end)
538 : size(1<<level)
539 {
540 vertex[0] = (end) ? size + 1 : 0;
541 for(int i = 1; i < dimension; ++ i)
542 vertex[i] = 0;
543 }
544
545 template<int dimension, class CoordType>
546 void
549 {
550 assert(vertex[0] <= size);
551 for(int i = dimension - 1; i >= 0; --i) {
552 ++vertex[i];
553 if(i == 0 || vertex[i] <= vertex[i-1])
554 break;
555 else
556 vertex[i] = 0;
557 }
558 }
559
560 template<int dimension, class CoordType>
561 bool
563 equals(const This &other) const
564 {
565 return size == other.size && vertex == other.vertex;
566 }
567
568 template<int dimension, class CoordType>
571 coords() const
572 {
573 Vertex ref = kuhnToReference(vertex, getPermutation<dimension>(0));
574
575 CoordVector coords;
576 for(int i = 0; i < dimension; ++i)
577 coords[i] = CoordType(ref[i]) / size;
578 return coords;
579 }
580
581 template<int dimension, class CoordType>
584 {
585 std::vector<CoordVector> corners(1);
586 corners[0] = (CoordVector)vertex;
587 return Geometry(GeometryType(0), corners);
588 }
589
590 template<int dimension, class CoordType>
591 int
597
598 // elements
599
600 template<int dimension, class CoordType>
601 class RefinementIteratorSpecial<dimension, CoordType, 0>
602 {
603 public:
607 typedef typename Refinement::template Codim<0>::Geometry Geometry;
609
610 RefinementIteratorSpecial(int level, bool end = false);
611
612 void increment();
613 bool equals(const This &other) const;
614
615 IndexVector vertexIndices() const;
616 int index() const;
617 CoordVector coords() const;
618
619 Geometry geometry () const;
620
621 private:
622 CoordVector global(const CoordVector &local) const;
623
624 protected:
625 typedef FieldVector<int, dimension> Vertex;
626 enum { nKuhnSimplices = Factorial<dimension>::factorial };
627
630 int size;
632 };
633
634 template<int dimension, class CoordType>
636 RefinementIteratorSpecial(int level, bool end)
637 : kuhnIndex(0), size(1<<level), index_(0)
638 {
639 for(int i = 0; i < dimension; ++i)
640 origin[i] = 0;
641 if(end) {
642 index_ = Refinement::nElements(level);
643 origin[0] = size;
644 }
645 }
646
647 template<int dimension, class CoordType>
648 void
651 {
652 assert(origin[0] < size);
653
654 ++index_;
655
656 while(1) {
657 ++kuhnIndex;
658 if(kuhnIndex == nKuhnSimplices) {
659 kuhnIndex = 0;
660 // increment origin
661 for(int i = dimension - 1; i >= 0; --i) {
662 ++origin[i];
663 if(i == 0 || origin[i] <= origin[i-1])
664 break;
665 else
666 origin[i] = 0;
667 }
668 }
669
670 // test whether the current simplex has any corner outside the kuhn0 simplex
671 FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
672 Vertex corner = origin;
673 bool outside = false;
674 for(int i = 0; i < dimension; ++i) {
675 // next corner
676 ++corner[perm[i]];
677 if(perm[i] > 0)
678 if(corner[perm[i]] > corner[perm[i]-1]) {
679 outside = true;
680 break;
681 }
682 }
683 if(!outside)
684 return;
685 }
686 }
687
688 template<int dimension, class CoordType>
689 bool
691 equals(const This &other) const
692 {
693 return size == other.size && index_ == other.index_;
694 }
695
696 template<int dimension, class CoordType>
699 vertexIndices() const
700 {
701 IndexVector indices;
702 FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
703 Vertex vertex = origin;
704 indices[0] = pointIndex(vertex);
705 for(int i = 0; i < dimension; ++i) {
706 ++vertex[perm[i]];
707 indices[i+1] = pointIndex(vertex);
708 }
709 if (kuhnIndex%2 == 1)
710 for(int i = 0; i < (dimension+1)/2; ++i) {
711 int t = indices[i];
712 indices[i] = indices[dimension-i];
713 indices[dimension-i] = t;
714 }
715 return indices;
716 }
717
718 template<int dimension, class CoordType>
719 int
721 index() const
722 {
723 return index_;
724 }
725
726 template<int dimension, class CoordType>
729 coords() const
730 {
732 ::simplex().position(0,0));
733 }
734
735 template<int dimension, class CoordType>
738 {
739 std::vector<CoordVector> corners(dimension+1);
740 CoordVector v;
743 for(int i = 0; i <= dimension; ++i)
744 corners[i] = global(refelem.position(i, dimension));
745 return Geometry(refelem.type(), corners);
746 }
747
748 template<int dimension, class CoordType>
751 global(const CoordVector &local) const {
752 CoordVector v =
753 referenceToKuhn(local, getPermutation<dimension>(kuhnIndex));
754 v += origin;
755 v /= (typename CoordVector::value_type)size;
756 return kuhnToReference(v, getPermutation<dimension>(0));
757 }
758
759 // common
760
761 template<int dimension, class CoordType>
762 template<int codimension>
763 class RefinementImp<dimension, CoordType>::Codim<codimension>::SubEntityIterator
764 : public ForwardIteratorFacade<typename RefinementImp<dimension, CoordType>::template Codim<codimension>::SubEntityIterator, int>,
765 public RefinementIteratorSpecial<dimension, CoordType, codimension>
766 {
767 public:
769
770 SubEntityIterator(int level, bool end = false);
771 };
772
773#ifndef DOXYGEN
774
775 template<int dimension, class CoordType>
776 template<int codimension>
778 SubEntityIterator(int level, bool end)
779 : RefinementIteratorSpecial<dimension, CoordType, codimension>(level, end)
780 {}
781
782#endif
783
784 } // namespace Simplex
785
786 } // namespace RefinementImp
787
788
789 namespace RefinementImp {
790
791 // ///////////////////////
792 //
793 // The refinement traits
794 //
795
796#ifndef DOXYGEN
797 template<unsigned topologyId, class CoordType, unsigned coerceToId,
798 int dim>
799 struct Traits<
800 topologyId, CoordType, coerceToId, dim,
801 typename std::enable_if<
802 ((GenericGeometry::SimplexTopology<dim>::type::id >> 1) ==
803 (topologyId >> 1) &&
804 (GenericGeometry::SimplexTopology<dim>::type::id >> 1) ==
805 (coerceToId >> 1)
806 )>::type
807 >
808 {
809 typedef Simplex::RefinementImp<dim, CoordType> Imp;
810 };
811#endif
812
813
814 } // namespace RefinementImp
815
816} // namespace Dune
817
818#endif //DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
This file contains the parts independent of a particular Refinement implementation.
Definition affinegeometry.hh:19
int pointIndex(const FieldVector< int, dimension > &point)
calculate the index of a given gridpoint within a Kuhn0 simplex
Definition simplex.cc:312
FieldVector< int, n > getPermutation(int m)
Calculate permutation from it's index.
Definition simplex.cc:325
int factorial(int n)
Calculate n!
Definition simplex.cc:282
int binomial(int upper, int lower)
calculate
Definition simplex.cc:294
FieldVector< CoordType, dimension > referenceToKuhn(FieldVector< CoordType, dimension > point, const FieldVector< int, dimension > &kuhn)
Map from the reference simplex to some Kuhn simplex.
Definition simplex.cc:383
FieldVector< CoordType, dimension > kuhnToReference(FieldVector< CoordType, dimension > point, const FieldVector< int, dimension > &kuhn)
Map from some Kuhn simplex to the reference simplex.
Definition simplex.cc:402
This class provides access to geometric and topological properties of a reference element.
Definition referenceelements.hh:53
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition referenceelements.hh:148
const GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition referenceelements.hh:130
Class providing access to the singletons of the reference elements.
Definition referenceelements.hh:455
static const ReferenceElement< ctype, dim > & simplex()
get simplex reference elements
Definition referenceelements.hh:466
Static tag representing a codimension.
Definition dimension.hh:22
Implement a MultiLinearGeometry with additional caching.
Definition multilineargeometry.hh:485
Codim< dimension >::SubEntityIterator VertexIterator
Definition simplex.cc:429
FieldVector< int, dimension+1 > IndexVector
Definition simplex.cc:432
CoordType ctype
Definition simplex.cc:425
static VertexIterator vEnd(int level)
Definition simplex.cc:471
static int nVertices(int level)
Definition simplex.cc:455
static int nElements(int level)
Definition simplex.cc:479
Codim< 0 >::SubEntityIterator ElementIterator
Definition simplex.cc:431
static ElementIterator eBegin(int level)
Definition simplex.cc:487
@ dimension
Definition simplex.cc:424
FieldVector< CoordType, dimension > CoordVector
Definition simplex.cc:430
static VertexIterator vBegin(int level)
Definition simplex.cc:463
static ElementIterator eEnd(int level)
Definition simplex.cc:495
Dune::CachedMultiLinearGeometry< CoordType, dimension-codimension, dimension > Geometry
Definition simplex.cc:449
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:514
Refinement::template Codim< dimension >::Geometry Geometry
Definition simplex.cc:516
RefinementIteratorSpecial< dimension, CoordType, dimension > This
Definition simplex.cc:517
FieldVector< int, dimension > Vertex
Definition simplex.cc:625
Refinement::template Codim< 0 >::Geometry Geometry
Definition simplex.cc:607
RefinementIteratorSpecial< dimension, CoordType, 0 > This
Definition simplex.cc:608
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:604
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:768
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:25