Loading...
Searching...
No Matches
simplex.cc
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
#define DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
// This file is part of DUNE, a Distributed and Unified Numerics Environment
// This file is copyright (C) 2005 Jorrit Fahlke <jorrit@jorrit.de>
// This file is licensed under version 2 of the GNU General Public License,
// with a special "runtime exception." See COPYING at the top of the source
// tree for the full licence.
#include <algorithm>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
#include "base.cc"
namespace Dune {
namespace RefinementImp {
namespace Simplex {
// //////////////////
//
//
template<int dimension>
int pointIndex(const FieldVector<int, dimension> &point)
{
int index = 0;
for(int i = 0; i < dimension; ++i)
index += Dune::binomial(dimension-i + point[i]-1, dimension-i);
return index;
}
template<int n>
FieldVector<int, n> getPermutation(int m)
{
FieldVector<int, n> perm;
for(int i = 0; i < n; ++i)
perm[i] = i;
int base = 1;
for(int i = 1; i <= n; ++i)
base *= i;
for(int i = n; i > 0; --i) {
base /= i;
int d = m / base;
m %= base;
int t = perm[i-1]; perm[i-1] = perm[i-1-d]; perm[i-1-d] = t;
}
return perm;
}
// map between the reference simplex and some arbitrary kuhn simplex (denoted by it's permutation)
template<int dimension, class CoordType>
FieldVector<CoordType, dimension>
referenceToKuhn(
FieldVector<CoordType, dimension> point,
const FieldVector<int, dimension> &kuhn)
{
for(int i = dimension - 1; i > 0; --i)
point[kuhn[i-1]] += point[kuhn[i]];
return point;
}
template<int dimension, class CoordType>
FieldVector<CoordType, dimension>
kuhnToReference(
FieldVector<CoordType, dimension> point,
const FieldVector<int, dimension> &kuhn)
{
for(int i = 0; i < dimension - 1; ++i)
point[kuhn[i]] -= point[kuhn[i+1]];
return point;
}
// /////////////////////////////////////////
//
// refinement implementation for simplices
//
template<int dimension_, class CoordType>
class RefinementImp
{
public:
constexpr static int dimension = dimension_;
typedef CoordType ctype;
template<int codimension>
struct Codim;
typedef typename Codim<dimension>::SubEntityIterator VertexIterator;
typedef FieldVector<CoordType, dimension> CoordVector;
typedef typename Codim<0>::SubEntityIterator ElementIterator;
typedef FieldVector<int, dimension+1> IndexVector;
static int nVertices(int nIntervals);
static VertexIterator vBegin(int nIntervals);
static VertexIterator vEnd(int nIntervals);
static int nElements(int nIntervals);
static ElementIterator eBegin(int nIntervals);
static ElementIterator eEnd(int nIntervals);
};
template<int dimension, class CoordType>
template<int codimension>
struct RefinementImp<dimension, CoordType>::Codim
{
class SubEntityIterator;
// We don't need the caching, but the uncached MultiLinearGeometry has bug FS#1209
typedef Dune::CachedMultiLinearGeometry<CoordType,dimension-codimension,dimension> Geometry;
};
template<int dimension, class CoordType>
int
RefinementImp<dimension, CoordType>::
nVertices(int nIntervals)
{
return Dune::binomial(dimension + nIntervals, (int)dimension);
}
template<int dimension, class CoordType>
typename RefinementImp<dimension, CoordType>::VertexIterator
RefinementImp<dimension, CoordType>::
vBegin(int nIntervals)
{
return VertexIterator(nIntervals);
}
template<int dimension, class CoordType>
typename RefinementImp<dimension, CoordType>::VertexIterator
RefinementImp<dimension, CoordType>::
vEnd(int nIntervals)
{
return VertexIterator(nIntervals, true);
}
template<int dimension, class CoordType>
int
RefinementImp<dimension, CoordType>::
nElements(int nIntervals)
{
return Dune::power(nIntervals, int(dimension));
}
template<int dimension, class CoordType>
typename RefinementImp<dimension, CoordType>::ElementIterator
RefinementImp<dimension, CoordType>::
eBegin(int nIntervals)
{
return ElementIterator(nIntervals);
}
template<int dimension, class CoordType>
typename RefinementImp<dimension, CoordType>::ElementIterator
RefinementImp<dimension, CoordType>::
eEnd(int nIntervals)
{
return ElementIterator(nIntervals, true);
}
// //////////////
//
// The iterator
//
template<int dimension, class CoordType, int codimension>
class RefinementIteratorSpecial;
// vertices
template<int dimension, class CoordType>
class RefinementIteratorSpecial<dimension, CoordType, dimension>
{
public:
typedef RefinementImp<dimension, CoordType> Refinement;
typedef typename Refinement::CoordVector CoordVector;
typedef typename Refinement::template Codim<dimension>::Geometry Geometry;
typedef RefinementIteratorSpecial<dimension, CoordType, dimension> This;
RefinementIteratorSpecial(int nIntervals, bool end = false);
void increment();
bool equals(const This &other) const;
CoordVector coords() const;
Geometry geometry () const;
int index() const;
protected:
typedef FieldVector<int, dimension> Vertex;
int size;
Vertex vertex;
};
template<int dimension, class CoordType>
RefinementIteratorSpecial<dimension, CoordType, dimension>::
RefinementIteratorSpecial(int nIntervals, bool end)
: size(nIntervals)
{
vertex[0] = (end) ? size + 1 : 0;
for(int i = 1; i < dimension; ++ i)
vertex[i] = 0;
}
template<int dimension, class CoordType>
void
RefinementIteratorSpecial<dimension, CoordType, dimension>::
increment()
{
assert(vertex[0] <= size);
for(int i = dimension - 1; i >= 0; --i) {
++vertex[i];
if(i == 0 || vertex[i] <= vertex[i-1])
break;
else
vertex[i] = 0;
}
}
template<int dimension, class CoordType>
bool
RefinementIteratorSpecial<dimension, CoordType, dimension>::
equals(const This &other) const
{
return size == other.size && vertex == other.vertex;
}
template<int dimension, class CoordType>
typename RefinementIteratorSpecial<dimension, CoordType, dimension>::CoordVector
RefinementIteratorSpecial<dimension, CoordType, dimension>::
coords() const
{
Vertex ref = kuhnToReference(vertex, getPermutation<dimension>(0));
CoordVector coords;
for(int i = 0; i < dimension; ++i)
coords[i] = CoordType(ref[i]) / size;
return coords;
}
template<int dimension, class CoordType>
typename RefinementIteratorSpecial<dimension, CoordType, dimension>::Geometry
RefinementIteratorSpecial<dimension, CoordType, dimension>::geometry () const
{
std::vector<CoordVector> corners(1);
corners[0] = (CoordVector)vertex;
return Geometry(GeometryTypes::vertex, corners);
}
template<int dimension, class CoordType>
int
RefinementIteratorSpecial<dimension, CoordType, dimension>::
index() const
{
return pointIndex(vertex);
}
// elements
template<int dimension, class CoordType>
class RefinementIteratorSpecial<dimension, CoordType, 0>
{
public:
typedef RefinementImp<dimension, CoordType> Refinement;
typedef typename Refinement::IndexVector IndexVector;
typedef typename Refinement::CoordVector CoordVector;
typedef typename Refinement::template Codim<0>::Geometry Geometry;
typedef RefinementIteratorSpecial<dimension, CoordType, 0> This;
RefinementIteratorSpecial(int nIntervals, bool end = false);
void increment();
bool equals(const This &other) const;
IndexVector vertexIndices() const;
int index() const;
CoordVector coords() const;
Geometry geometry () const;
private:
CoordVector global(const CoordVector &local) const;
protected:
typedef FieldVector<int, dimension> Vertex;
constexpr static int nKuhnIntervals = Dune::factorial(dimension);
Vertex origin;
int kuhnIndex;
int size;
int index_;
};
template<int dimension, class CoordType>
RefinementIteratorSpecial<dimension, CoordType, 0>::
RefinementIteratorSpecial(int nIntervals, bool end)
: kuhnIndex(0), size(nIntervals), index_(0)
{
for(int i = 0; i < dimension; ++i)
origin[i] = 0;
if(end) {
index_ = Refinement::nElements(nIntervals);
origin[0] = size;
}
}
template<int dimension, class CoordType>
void
RefinementIteratorSpecial<dimension, CoordType, 0>::
increment()
{
assert(origin[0] < size);
++index_;
while(1) {
++kuhnIndex;
if(kuhnIndex == nKuhnIntervals) {
kuhnIndex = 0;
// increment origin
for(int i = dimension - 1; i >= 0; --i) {
++origin[i];
if(i == 0 || origin[i] <= origin[i-1])
break;
else
origin[i] = 0;
}
}
// test whether the current simplex has any corner outside the kuhn0 simplex
FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
Vertex corner = origin;
bool outside = false;
for(int i = 0; i < dimension; ++i) {
// next corner
++corner[perm[i]];
if(perm[i] > 0)
if(corner[perm[i]] > corner[perm[i]-1]) {
outside = true;
break;
}
}
if(!outside)
return;
}
}
template<int dimension, class CoordType>
bool
RefinementIteratorSpecial<dimension, CoordType, 0>::
equals(const This &other) const
{
return size == other.size && index_ == other.index_;
}
template<int dimension, class CoordType>
typename RefinementIteratorSpecial<dimension, CoordType, 0>::IndexVector
RefinementIteratorSpecial<dimension, CoordType, 0>::
vertexIndices() const
{
IndexVector indices;
FieldVector<int, dimension> perm = getPermutation<dimension>(kuhnIndex);
Vertex vertex = origin;
indices[0] = pointIndex(vertex);
for(int i = 0; i < dimension; ++i) {
++vertex[perm[i]];
indices[i+1] = pointIndex(vertex);
}
if (kuhnIndex%2 == 1)
for(int i = 0; i < (dimension+1)/2; ++i) {
int t = indices[i];
indices[i] = indices[dimension-i];
indices[dimension-i] = t;
}
return indices;
}
template<int dimension, class CoordType>
int
RefinementIteratorSpecial<dimension, CoordType, 0>::
index() const
{
return index_;
}
template<int dimension, class CoordType>
typename RefinementIteratorSpecial<dimension, CoordType, 0>::CoordVector
RefinementIteratorSpecial<dimension, CoordType, 0>::
coords() const
{
return global(ReferenceElements<CoordType, dimension>
::simplex().position(0,0));
}
template<int dimension, class CoordType>
typename RefinementIteratorSpecial<dimension, CoordType, 0>::Geometry
RefinementIteratorSpecial<dimension, CoordType, 0>::geometry () const
{
std::vector<CoordVector> corners(dimension+1);
auto refelem = ReferenceElements<CoordType, dimension>::simplex();
for(int i = 0; i <= dimension; ++i)
corners[i] = global(refelem.position(i, dimension));
return Geometry(refelem.type(), corners);
}
template<int dimension, class CoordType>
typename RefinementIteratorSpecial<dimension, CoordType, 0>::CoordVector
RefinementIteratorSpecial<dimension, CoordType, 0>::
global(const CoordVector &local) const {
CoordVector v =
referenceToKuhn(local, getPermutation<dimension>(kuhnIndex));
v += origin;
v /= (typename CoordVector::value_type)size;
return kuhnToReference(v, getPermutation<dimension>(0));
}
// common
template<int dimension, class CoordType>
template<int codimension>
class RefinementImp<dimension, CoordType>::Codim<codimension>::SubEntityIterator
: public ForwardIteratorFacade<typename RefinementImp<dimension, CoordType>::template Codim<codimension>::SubEntityIterator, int>,
public RefinementIteratorSpecial<dimension, CoordType, codimension>
{
public:
typedef RefinementImp<dimension, CoordType> Refinement;
SubEntityIterator(int nIntervals, bool end = false);
};
#ifndef DOXYGEN
template<int dimension, class CoordType>
template<int codimension>
RefinementImp<dimension, CoordType>::Codim<codimension>::SubEntityIterator::
SubEntityIterator(int nIntervals, bool end)
: RefinementIteratorSpecial<dimension, CoordType, codimension>(nIntervals, end)
{}
#endif
} // namespace Simplex
} // namespace RefinementImp
namespace RefinementImp {
// ///////////////////////
//
// The refinement traits
//
#ifndef DOXYGEN
template<unsigned topologyId, class CoordType, unsigned coerceToId,
int dim>
struct Traits<
topologyId, CoordType, coerceToId, dim,
typename std::enable_if<
((GeometryTypes::simplex(dim).id() >> 1) ==
(topologyId >> 1) &&
(GeometryTypes::simplex(dim).id() >> 1) ==
(coerceToId >> 1)
)>::type
>
{
typedef Simplex::RefinementImp<dim, CoordType> Imp;
};
#endif
} // namespace RefinementImp
} // namespace Dune
#endif //DUNE_GRID_COMMON_REFINEMENT_SIMPLEX_CC
base.cc
This file contains the parts independent of a particular Refinement implementation.
multilineargeometry.hh
referenceelements.hh
type.hh
A unique label for each type of element that can occur in a grid.
Dune::GeometryTypes::vertex
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition type.hh:492
std
STL namespace.
Dune
Definition affinegeometry.hh:22
Dune::RefinementImp
This namespace contains the implementation of Refinement.
Definition base.cc:29
Dune::RefinementImp::Simplex
This namespace contains the Refinement implementation for simplices (triangles, tetrahedrons....
Definition simplex.cc:264
Dune::RefinementImp::Simplex::pointIndex
int pointIndex(const FieldVector< int, dimension > &point)
calculate the index of a given gridpoint within a Kuhn0 simplex
Definition simplex.cc:280
Dune::RefinementImp::Simplex::getPermutation
FieldVector< int, n > getPermutation(int m)
Calculate permutation from it's index.
Definition simplex.cc:293
Dune::RefinementImp::Simplex::referenceToKuhn
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
Dune::RefinementImp::Simplex::kuhnToReference
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
Dune::Geo::ReferenceElements
Class providing access to the singletons of the reference elements.
Definition referenceelements.hh:128
Dune::Codim
Static tag representing a codimension.
Definition dimension.hh:24
Dune::CachedMultiLinearGeometry
Implement a MultiLinearGeometry with additional caching.
Definition multilineargeometry.hh:526
Dune::RefinementImp::Simplex::RefinementImp
Definition simplex.cc:361
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::VertexIterator
Codim< dimension >::SubEntityIterator VertexIterator
Definition simplex.cc:368
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::IndexVector
FieldVector< int, dimension+1 > IndexVector
Definition simplex.cc:371
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::ctype
CoordType ctype
Definition simplex.cc:364
Dune::RefinementImp::Simplex::RefinementImp::nVertices
static int nVertices(int nIntervals)
Definition simplex.cc:394
Dune::RefinementImp::Simplex::RefinementImp::nElements
static int nElements(int nIntervals)
Definition simplex.cc:418
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::Codim< codimension >::SubEntityIterator
SubEntityIterator(int nIntervals, bool end=false)
Dune::RefinementImp::Simplex::RefinementImp::eEnd
static ElementIterator eEnd(int nIntervals)
Definition simplex.cc:434
Dune::RefinementImp::Simplex::RefinementImp::vEnd
static VertexIterator vEnd(int nIntervals)
Definition simplex.cc:410
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::Codim< codimension >::Refinement
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:705
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::ElementIterator
Codim< 0 >::SubEntityIterator ElementIterator
Definition simplex.cc:370
Dune::RefinementImp::Simplex::RefinementImp::vBegin
static VertexIterator vBegin(int nIntervals)
Definition simplex.cc:402
Dune::RefinementImp::Simplex::RefinementImp::eBegin
static ElementIterator eBegin(int nIntervals)
Definition simplex.cc:426
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::CoordVector
FieldVector< CoordType, dimension > CoordVector
Definition simplex.cc:369
Dune::RefinementImp::Simplex::RefinementImp< dimension, CoordType >::dimension
static constexpr int dimension
Definition simplex.cc:363
Dune::RefinementImp::Simplex::RefinementImp::Codim
Definition simplex.cc:385
Dune::RefinementImp::Simplex::RefinementIteratorSpecial
Definition simplex.cc:445
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::CoordVector
Refinement::CoordVector CoordVector
Definition simplex.cc:454
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::Refinement
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:453
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::vertex
Vertex vertex
Definition simplex.cc:471
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::index
int index() const
Definition simplex.cc:532
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::increment
void increment()
Definition simplex.cc:487
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::RefinementIteratorSpecial
RefinementIteratorSpecial(int nIntervals, bool end=false)
Definition simplex.cc:476
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::equals
bool equals(const This &other) const
Definition simplex.cc:502
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::coords
CoordVector coords() const
Definition simplex.cc:510
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::Geometry
Refinement::template Codim< dimension >::Geometry Geometry
Definition simplex.cc:455
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::geometry
Geometry geometry() const
Definition simplex.cc:522
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::This
RefinementIteratorSpecial< dimension, CoordType, dimension > This
Definition simplex.cc:456
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::Vertex
FieldVector< int, dimension > Vertex
Definition simplex.cc:468
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, dimension >::size
int size
Definition simplex.cc:470
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::vertexIndices
IndexVector vertexIndices() const
Definition simplex.cc:638
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::kuhnIndex
int kuhnIndex
Definition simplex.cc:568
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::increment
void increment()
Definition simplex.cc:589
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::Vertex
FieldVector< int, dimension > Vertex
Definition simplex.cc:564
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::RefinementIteratorSpecial
RefinementIteratorSpecial(int nIntervals, bool end=false)
Definition simplex.cc:575
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::origin
Vertex origin
Definition simplex.cc:567
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::geometry
Geometry geometry() const
Definition simplex.cc:676
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::index_
int index_
Definition simplex.cc:570
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::index
int index() const
Definition simplex.cc:660
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::nKuhnIntervals
static constexpr int nKuhnIntervals
Definition simplex.cc:565
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::Geometry
Refinement::template Codim< 0 >::Geometry Geometry
Definition simplex.cc:546
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::coords
CoordVector coords() const
Definition simplex.cc:668
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::IndexVector
Refinement::IndexVector IndexVector
Definition simplex.cc:544
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::CoordVector
Refinement::CoordVector CoordVector
Definition simplex.cc:545
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::This
RefinementIteratorSpecial< dimension, CoordType, 0 > This
Definition simplex.cc:547
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::Refinement
RefinementImp< dimension, CoordType > Refinement
Definition simplex.cc:543
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::size
int size
Definition simplex.cc:569
Dune::RefinementImp::Simplex::RefinementIteratorSpecial< dimension, CoordType, 0 >::equals
bool equals(const This &other) const
Definition simplex.cc:630
Dune::ReferenceElements::simplex
static const ReferenceElement & simplex()
get simplex reference elements
Definition referenceelements.hh:162