dune-geometry 2.10
Loading...
Searching...
No Matches
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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
6#define DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
7
8// This file is part of DUNE, a Distributed and Unified Numerics Environment
9// This file is copyright (C) 2005 Jorrit Fahlke <jorrit@jorrit.de>
10// This file is licensed under version 2 of the GNU General Public License,
11// with a special "runtime exception." See COPYING at the top of the source
12// tree for the full licence.
13
243#include <algorithm>
244
245#include <dune/common/fvector.hh>
246#include <dune/common/math.hh>
247
250#include <dune/geometry/type.hh>
251
252#include "base.cc"
253
254namespace Dune {
255
256 namespace RefinementImp {
257
264 namespace Simplex {
265
266 // //////////////////
267 //
269 //
270
272
279 template<int dimension>
280 int pointIndex(const FieldVector<int, dimension> &point)
281 {
282 int index = 0;
283 for(int i = 0; i < dimension; ++i)
284 index += Dune::binomial(dimension-i + point[i]-1, dimension-i);
285 return index;
286 }
287
292 template<int n>
293 FieldVector<int, n> getPermutation(int m)
294 {
295 FieldVector<int, n> perm;
296 for(int i = 0; i < n; ++i)
297 perm[i] = i;
298
299 int base = 1;
300 for(int i = 1; i <= n; ++i)
301 base *= i;
302
303 for(int i = n; i > 0; --i) {
304 base /= i;
305 int d = m / base;
306 m %= base;
307 int t = perm[i-1]; perm[i-1] = perm[i-1-d]; perm[i-1-d] = t;
308 }
309 return perm;
310 }
311
312 // map between the reference simplex and some arbitrary kuhn simplex (denoted by it's permutation)
320 template<int dimension, class CoordType>
321 FieldVector<CoordType, dimension>
323 FieldVector<CoordType, dimension> point,
325 const FieldVector<int, dimension> &kuhn)
326 {
327 for(int i = dimension - 1; i > 0; --i)
328 point[kuhn[i-1]] += point[kuhn[i]];
329 return point;
330 }
331
339 template<int dimension, class CoordType>
340 FieldVector<CoordType, dimension>
342 FieldVector<CoordType, dimension> point,
344 const FieldVector<int, dimension> &kuhn)
345 {
346 for(int i = 0; i < dimension - 1; ++i)
347 point[kuhn[i]] -= point[kuhn[i+1]];
348 return point;
349 }
350
351
353
354 // /////////////////////////////////////////
355 //
356 // refinement implementation for simplices
357 //
358
359 template<int dimension_, class CoordType>
361 {
362 public:
363 constexpr static int dimension = dimension_;
364 typedef CoordType ctype;
365
366 template<int codimension>
367 struct Codim;
369 typedef FieldVector<CoordType, dimension> CoordVector;
371 typedef FieldVector<int, dimension+1> IndexVector;
372
373 static int nVertices(int nIntervals);
374 static VertexIterator vBegin(int nIntervals);
375 static VertexIterator vEnd(int nIntervals);
376
377 static int nElements(int nIntervals);
378 static ElementIterator eBegin(int nIntervals);
379 static ElementIterator eEnd(int nIntervals);
380 };
381
382 template<int dimension, class CoordType>
383 template<int codimension>
384 struct RefinementImp<dimension, CoordType>::Codim
385 {
386 class SubEntityIterator;
387 // We don't need the caching, but the uncached MultiLinearGeometry has bug FS#1209
389 };
390
391 template<int dimension, class CoordType>
392 int
394 nVertices(int nIntervals)
395 {
396 return Dune::binomial(dimension + nIntervals, (int)dimension);
397 }
398
399 template<int dimension, class CoordType>
402 vBegin(int nIntervals)
403 {
404 return VertexIterator(nIntervals);
405 }
406
407 template<int dimension, class CoordType>
410 vEnd(int nIntervals)
411 {
412 return VertexIterator(nIntervals, true);
413 }
414
415 template<int dimension, class CoordType>
416 int
418 nElements(int nIntervals)
419 {
420 return Dune::power(nIntervals, int(dimension));
421 }
422
423 template<int dimension, class CoordType>
426 eBegin(int nIntervals)
427 {
428 return ElementIterator(nIntervals);
429 }
430
431 template<int dimension, class CoordType>
434 eEnd(int nIntervals)
435 {
436 return ElementIterator(nIntervals, true);
437 }
438
439 // //////////////
440 //
441 // The iterator
442 //
443
444 template<int dimension, class CoordType, int codimension>
446
447 // vertices
448
449 template<int dimension, class CoordType>
450 class RefinementIteratorSpecial<dimension, CoordType, dimension>
451 {
452 public:
454 typedef typename Refinement::CoordVector CoordVector;
455 typedef typename Refinement::template Codim<dimension>::Geometry Geometry;
457
458 RefinementIteratorSpecial(int nIntervals, bool end = false);
459
460 void increment();
461 bool equals(const This &other) const;
462
463 CoordVector coords() const;
464 Geometry geometry () const;
465
466 int index() const;
467 protected:
468 typedef FieldVector<int, dimension> Vertex;
469
470 int size;
472 };
473
474 template<int dimension, class CoordType>
476 RefinementIteratorSpecial(int nIntervals, bool end)
477 : size(nIntervals)
478 {
479 vertex[0] = (end) ? size + 1 : 0;
480 for(int i = 1; i < dimension; ++ i)
481 vertex[i] = 0;
482 }
483
484 template<int dimension, class CoordType>
485 void
488 {
489 assert(vertex[0] <= size);
490 for(int i = dimension - 1; i >= 0; --i) {
491 ++vertex[i];
492 if(i == 0 || vertex[i] <= vertex[i-1])
493 break;
494 else
495 vertex[i] = 0;
496 }
497 }
498
499 template<int dimension, class CoordType>
500 bool
502 equals(const This &other) const
503 {
504 return size == other.size && vertex == other.vertex;
505 }
506
507 template<int dimension, class CoordType>
510 coords() const
511 {
512 Vertex ref = kuhnToReference(vertex, getPermutation<dimension>(0));
513
514 CoordVector coords;
515 for(int i = 0; i < dimension; ++i)
516 coords[i] = CoordType(ref[i]) / size;
517 return coords;
518 }
519
520 template<int dimension, class CoordType>
523 {
524 std::vector<CoordVector> corners(1);
525 corners[0] = (CoordVector)vertex;
526 return Geometry(GeometryTypes::vertex, corners);
527 }
528
529 template<int dimension, class CoordType>
530 int
536
537 // elements
538
539 template<int dimension, class CoordType>
540 class RefinementIteratorSpecial<dimension, CoordType, 0>
541 {
542 public:
546 typedef typename Refinement::template Codim<0>::Geometry Geometry;
548
549 RefinementIteratorSpecial(int nIntervals, bool end = false);
550
551 void increment();
552 bool equals(const This &other) const;
553
554 IndexVector vertexIndices() const;
555 int index() const;
556 CoordVector coords() const;
557
558 Geometry geometry () const;
559
560 private:
561 CoordVector global(const CoordVector &local) const;
562
563 protected:
564 typedef FieldVector<int, dimension> Vertex;
565 constexpr static int nKuhnIntervals = Dune::factorial(dimension);
566
569 int size;
571 };
572
573 template<int dimension, class CoordType>
575 RefinementIteratorSpecial(int nIntervals, bool end)
576 : kuhnIndex(0), size(nIntervals), index_(0)
577 {
578 for(int i = 0; i < dimension; ++i)
579 origin[i] = 0;
580 if(end) {
581 index_ = Refinement::nElements(nIntervals);
582 origin[0] = size;
583 }
584 }
585
586 template<int dimension, class CoordType>
587 void
590 {
591 assert(origin[0] < size);
592
593 ++index_;
594
595 while(1) {
596 ++kuhnIndex;
597 if(kuhnIndex == nKuhnIntervals) {
598 kuhnIndex = 0;
599 // increment origin
600 for(int i = dimension - 1; i >= 0; --i) {
601 ++origin[i];
602 if(i == 0 || origin[i] <= origin[i-1])
603 break;
604 else
605 origin[i] = 0;
606 }
607 }
608
609 // test whether the current simplex has any corner outside the kuhn0 simplex
610 FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
611 Vertex corner = origin;
612 bool outside = false;
613 for(int i = 0; i < dimension; ++i) {
614 // next corner
615 ++corner[perm[i]];
616 if(perm[i] > 0)
617 if(corner[perm[i]] > corner[perm[i]-1]) {
618 outside = true;
619 break;
620 }
621 }
622 if(!outside)
623 return;
624 }
625 }
626
627 template<int dimension, class CoordType>
628 bool
630 equals(const This &other) const
631 {
632 return size == other.size && index_ == other.index_;
633 }
634
635 template<int dimension, class CoordType>
638 vertexIndices() const
639 {
640 IndexVector indices;
641 FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
642 Vertex vertex = origin;
643 indices[0] = pointIndex(vertex);
644 for(int i = 0; i < dimension; ++i) {
645 ++vertex[perm[i]];
646 indices[i+1] = pointIndex(vertex);
647 }
648 if (kuhnIndex%2 == 1)
649 for(int i = 0; i < (dimension+1)/2; ++i) {
650 int t = indices[i];
651 indices[i] = indices[dimension-i];
652 indices[dimension-i] = t;
653 }
654 return indices;
655 }
656
657 template<int dimension, class CoordType>
658 int
660 index() const
661 {
662 return index_;
663 }
664
665 template<int dimension, class CoordType>
668 coords() const
669 {
671 ::simplex().position(0,0));
672 }
673
674 template<int dimension, class CoordType>
677 {
678 std::vector<CoordVector> corners(dimension+1);
680 for(int i = 0; i <= dimension; ++i)
681 corners[i] = global(refelem.position(i, dimension));
682 return Geometry(refelem.type(), corners);
683 }
684
685 template<int dimension, class CoordType>
688 global(const CoordVector &local) const {
689 CoordVector v =
690 referenceToKuhn(local, getPermutation<dimension>(kuhnIndex));
691 v += origin;
692 v /= (typename CoordVector::value_type)size;
693 return kuhnToReference(v, getPermutation<dimension>(0));
694 }
695
696 // common
697
698 template<int dimension, class CoordType>
699 template<int codimension>
700 class RefinementImp<dimension, CoordType>::Codim<codimension>::SubEntityIterator
701 : public ForwardIteratorFacade<typename RefinementImp<dimension, CoordType>::template Codim<codimension>::SubEntityIterator, int>,
702 public RefinementIteratorSpecial<dimension, CoordType, codimension>
703 {
704 public:
706
707 SubEntityIterator(int nIntervals, bool end = false);
708 };
709
710#ifndef DOXYGEN
711
712 template<int dimension, class CoordType>
713 template<int codimension>
715 SubEntityIterator(int nIntervals, bool end)
716 : RefinementIteratorSpecial<dimension, CoordType, codimension>(nIntervals, end)
717 {}
718
719#endif
720
721 } // namespace Simplex
722
723 } // namespace RefinementImp
724
725
726 namespace RefinementImp {
727
728 // ///////////////////////
729 //
730 // The refinement traits
731 //
732
733#ifndef DOXYGEN
734 template<unsigned topologyId, class CoordType, unsigned coerceToId,
735 int dim>
736 struct Traits<
737 topologyId, CoordType, coerceToId, dim,
738 typename std::enable_if<
739 ((GeometryTypes::simplex(dim).id() >> 1) ==
740 (topologyId >> 1) &&
741 (GeometryTypes::simplex(dim).id() >> 1) ==
742 (coerceToId >> 1)
743 )>::type
744 >
745 {
746 typedef Simplex::RefinementImp<dim, CoordType> Imp;
747 };
748#endif
749
750
751 } // namespace RefinementImp
752
753} // namespace Dune
754
755#endif //DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
A unique label for each type of element that can occur in a grid.
This file contains the parts independent of a particular Refinement implementation.
STL namespace.
Definition affinegeometry.hh:21
int pointIndex(const FieldVector< int, dimension > &point)
calculate the index of a given gridpoint within a Kuhn0 simplex
Definition simplex.cc:280
FieldVector< int, n > getPermutation(int m)
Calculate permutation from it's index.
Definition simplex.cc:293
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:322
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:341
Class providing access to the singletons of the reference elements.
Definition referenceelements.hh:128
static const ReferenceElement & simplex()
get simplex reference elements
Definition referenceelements.hh:162
Static tag representing a codimension.
Definition dimension.hh:24
Implement a MultiLinearGeometry with additional caching.
Definition multilineargeometry.hh:526
Codim< dimension >::SubEntityIterator VertexIterator
Definition simplex.cc:368
FieldVector< int, dimension+1 > IndexVector
Definition simplex.cc:371
CoordType ctype
Definition simplex.cc:364
static int nVertices(int nIntervals)
Definition simplex.cc:394
static int nElements(int nIntervals)
Definition simplex.cc:418
static ElementIterator eEnd(int nIntervals)
Definition simplex.cc:434
static VertexIterator vEnd(int nIntervals)
Definition simplex.cc:410
Codim< 0 >::SubEntityIterator ElementIterator
Definition simplex.cc:370
static VertexIterator vBegin(int nIntervals)
Definition simplex.cc:402
static ElementIterator eBegin(int nIntervals)
Definition simplex.cc:426
FieldVector< CoordType, dimension > CoordVector
Definition simplex.cc:369
static constexpr int dimension
Definition simplex.cc:363
Dune::CachedMultiLinearGeometry< CoordType, dimension-codimension, dimension > Geometry
Definition simplex.cc:388
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:453
Refinement::template Codim< dimension >::Geometry Geometry
Definition simplex.cc:455
RefinementIteratorSpecial< dimension, CoordType, dimension > This
Definition simplex.cc:456
FieldVector< int, dimension > Vertex
Definition simplex.cc:564
Refinement::template Codim< 0 >::Geometry Geometry
Definition simplex.cc:546
RefinementIteratorSpecial< dimension, CoordType, 0 > This
Definition simplex.cc:547
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:543
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:705