maia::grid::hilbert Namespace Reference

Namespaces
namespace	detail_

Functions
template<MInt nDim, typename FloatType , typename IdType >
IdType	index (const FloatType *const x, const IdType level)
	Return Hilbert index for given location and level in 2D or 3D. More...
template<MInt nDim, typename FloatType , typename IdType >
MBool	coordinatesToTreeId (IdType &treeId, const FloatType *const x, const IdType level, const FloatType *const centerOfGravity, FloatType const lengthOnLevel0)
template<MInt nDim, typename FloatType , typename IdType >
void	treeIdToCoordinates (FloatType *const x, const IdType treeId, const IdType level, const FloatType *const centerOfGravity, const FloatType lengthOnLevel0)

Function Documentation

coordinatesToTreeId()

template<MInt nDim, typename FloatType , typename IdType >

MBool maia::grid::hilbert::coordinatesToTreeId	(	IdType &	treeId,
		const FloatType *const	x,
		const IdType	level,
		const FloatType *const	centerOfGravity,
		FloatType const	lengthOnLevel0
	)

Author

Lennart

Date

Definition at line 188 of file hilbert.h.

                                                          {
  // Check that level is at least 0
  ASSERT(level >= 0, "Level must be at least zero");
 
  // Check for integer overflow ("- 1" to account for signed integers, cast to
  // avoid "unsigned/signed comparison" warning)
  ASSERT(nDim * level <= CHAR_BIT * static_cast<IdType>(sizeof(IdType)) - 1, "Integer overflow: level too large");
 
  MBool exists = true;
  treeId = 0;
  for(IdType k = 0; k < nDim; k++) {
    if(x[k] < centerOfGravity[k] - 0.5 * lengthOnLevel0 || x[k] > centerOfGravity[k] + 0.5 * lengthOnLevel0) {
      exists = false;
    }
  }
  if(!exists) {
    treeId = std::numeric_limits<IdType>::is_signed ? -1 : ~0;
  } else {
    IdType bitCount = 0;
    IdType coord[3];
    constexpr IdType id2 = (IdType)2;
    for(IdType k = 0; k < nDim; k++) {
      coord[k] = (IdType)((x[k] - centerOfGravity[k] + 0.5 * lengthOnLevel0) * FPOW2(level) / lengthOnLevel0);
    }
    for(IdType lvl = level - 1; lvl >= 0; lvl--) {
      for(IdType k = 0; k < nDim; k++) {
        IdType bit = (coord[k] / (IdType)IPOW2(lvl)) % id2;
        treeId |= bit << bitCount;
        bitCount++;
      }
    }
  }
  return exists;
}

index()

template<MInt nDim, typename FloatType , typename IdType >

IdType maia::grid::hilbert::index	(	const FloatType *const	x,
		const IdType	level
	)

Author

Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de, Peter Philippen

Date

2014-12-08

Template Parameters

nDim	Dimension for which the Hilbert index should be generated.
FloatType	The floating point type for the coordinates.
IdType	The integer type that should be used for the index.

Parameters

[in]	x	Pointer to nDim coordinates for which the index is calculated. The coordinates must be within the unit square (2D, [0,1]x[0,1])/unit cube (3D, [0,1]x[0,1]x[0,1]).
[in]	level	Refinement level (not Hilbert level) for which the index should be calculated.

Returns

The calculated Hilbert index.

The algorithms used here are taken from

Peter Philippen. Investigation of two different domain decompositioning
methods for a cartesian flow solver. Studienarbeit, 2009.

The original implementation is by Peter, it was only simplified (no recursion!) and its robustness against bad input parameters increased. Also, the original implementation did not have any comments.

This is the basic idea:

• determine the quadrant/octant of the unit hypercube in which the point is located
• calculate its Hilbert index and multiply with the level weight
• mirror/scale/translate quadrant/octant to new unit hypercube
• repeat from beginning

Definition at line 165 of file hilbert.h.

                                                           {
  // Check that level is at least 0
  ASSERT(level >= 0, "Level must be at least zero");
 
  // Check for integer overflow ("- 1" to account for signed integers, cast to
  // avoid "unsigned/signed comparison" warning)
  ASSERT(nDim * level <= CHAR_BIT * static_cast<IdType>(sizeof(IdType)) - 1, "Integer overflow: level too large");
 
  // Assert that coordinates are in unit hypercube
  ASSERT(x[0] >= 0.0 && x[0] <= 1.0, "Coordinate 0 out of bounds.");
  ASSERT(x[1] >= 0.0 && x[1] <= 1.0, "Coordinate 1 out of bounds.");
  IF_CONSTEXPR(nDim == 3) { ASSERT(x[2] >= 0.0 && x[2] <= 1.0, "Coordinate 2 out of bounds."); }
 
  // Call template function
  return detail_::Impl<nDim, FloatType, IdType>::f(x, level);
}

treeIdToCoordinates()

template<MInt nDim, typename FloatType , typename IdType >

void maia::grid::hilbert::treeIdToCoordinates	(	FloatType *const	x,
		const IdType	treeId,
		const IdType	level,
		const FloatType *const	centerOfGravity,
		const FloatType	lengthOnLevel0
	)

Author

Lennart

Date

Definition at line 232 of file hilbert.h.

                                                         {
  // Check that level is at least 0
  ASSERT(level >= 0, "Level must be at least zero");
 
  // Check for integer overflow ("- 1" to account for signed integers, cast to
  // avoid "unsigned/signed comparison" warning)
  ASSERT(nDim * level <= CHAR_BIT * static_cast<IdType>(sizeof(IdType)) - 1, "Integer overflow: level too large");
 
  IdType bitCount = 0;
  constexpr IdType id1 = (IdType)1;
  for(IdType j = 0; j < nDim; j++) {
    x[j] = centerOfGravity[j];
  }
  for(IdType lvl = 1; lvl <= level; lvl++) {
    for(IdType j = 0; j < nDim; j++) {
      IdType tmpBit = (treeId >> bitCount) & id1;
      x[j] += (tmpBit ? 1.0 : -1.0) * lengthOnLevel0 * FFPOW2(lvl + 1);
      bitCount++;
    }
  }
}