#include <cstring>
#include <algorithm>
#include <vector>
+#include <limits>
#include <iostream>
/** A simplicial complex is a std::vector of simplices such that each face is
* also part of the complex. Every simplex has dimension at most MAXDIM. The
* indices of simplices resp. their faces are of type IT. To each simplex a
- * value is assigend, which is of type VT. */
-template<size_t MAXDIM, class IT=uint32_t, class VT=double>
+ * value is assigend, which is of type VT. When a SimplicialComplex is
+ * instantiated, a single (-1) dimensional simplex is automatically created.
+ * Each 0-dimensional simplex automatically has this simplex as its face.
+ * Consequently, the innner class SimplexOrder gives the extended boundary
+ * matrix. */
+template<int MAXDIM, class IT=uint32_t, class VT=double>
class SimplicialComplex {
public:
/** A simplex of the complex. */
struct Simplex {
/** Dimension of the simplex. */
- size_t dim;
+ int dim;
/** The indices of the faces of the simplex. */
index_type faces[MAXDIM+1];
/** The value of the simplex. */
value_type value;
/** Create a new simplex with dimension 'dim', (dim+1)-faces and
- * its value. */
- static Simplex create(size_t dim, index_type* faces, value_type value) {
- assert(dim <= MAXDIM);
+ * its value. If simpley is 0-dimensional, its face is
+ * automatically set to one (-1)-dimensional simplex. */
+ static Simplex create(int dim, index_type* faces, value_type value) {
+ assert(0 <= dim && dim <= MAXDIM);
Simplex s;
s.dim = dim;
s.value = value;
- memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type));
+
+ if (dim > 0)
+ memcpy(s.faces, faces, faceCountByDim(dim)*sizeof(index_type));
+ else
+ s.faces[0] = 0;
+
+ return s;
+ }
+
+ /** Create a (-1)-dimensional simplex. It has the lowest possible value. */
+ static Simplex create_minusonedim_simplex() {
+ Simplex s;
+
+ s.dim = -1;
+ s.faces[0] = 0;
+ s.value = std::numeric_limits<VT>::has_infinity
+ ? -std::numeric_limits<VT>::infinity()
+ : std::numeric_limits<VT>::min();
return s;
}
/** Get number of faces. */
- size_t face_count() const {
- return face_count_bydim(dim);
+ size_t faceCount() const {
+ return faceCountByDim(dim);
}
/** Get number of faces of a dim-dimensional simplex. */
- static size_t face_count_bydim(size_t dim) {
- if (dim == 0)
- return 0;
+ static size_t faceCountByDim(int dim) {
+ assert(-1 <= dim && dim <= MAXDIM);
return dim + 1;
}
};
assert(size() == c.size());
for (unsigned i=0; i < size(); ++i)
- for (unsigned f=0; f < getSimplex(i).face_count(); ++f)
+ for (unsigned f=0; f < getSimplex(i).faceCount(); ++f)
if (revorder[getSimplex(i).faces[f]] >= i)
return false;
BoundaryMatrix mat(size());
for (unsigned c=0; c < size(); ++c)
- for(unsigned r=0; r < getSimplex(c).face_count(); ++r)
+ for(unsigned r=0; r < getSimplex(c).faceCount(); ++r)
mat.set(revorder[getSimplex(c).faces[r]], c, true);
return mat;
};
public:
+ SimplicialComplex() {
+ // Add the one minus-one dimensional simplex
+ addSimplex(Simplex::create_minusonedim_simplex());
+ }
+
/** Return number of simplices. */
size_t size() const {
return simplices.size();
/** Add a simplex to the complex. The dimension of the faces must be
* dim-1, and they must already be part of the complex. */
- void addSimplex(size_t dim, index_type* faces, value_type value) {
+ void addSimplex(int dim, index_type* faces, value_type value) {
addSimplex(Simplex::create(dim, faces, value));
}
* dim-1, and they must already be part of the complex. */
void addSimplex(Simplex s) {
// Check requirements for faces
- for (unsigned i=0; i < s.face_count(); ++i) {
+ for (unsigned i=0; i < s.faceCount(); ++i) {
// Faces are already in complex.
assert(s.faces[i] < size());
// Faces have dimension dim-1
/** Return true iff for each simplex i with dimension dim it holds that
* the faces of i are contained in the complex and have dimension dim-1. */
bool isComplex() const {
- typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
for (unsigned i=0; i < size(); ++i) {
const Simplex &s = simplices[i];
- for (unsigned f=0; f < s.face_count(); ++f) {
+ for (unsigned f=0; f < s.faceCount(); ++f) {
if (s.faces[f] >= size())
return false;
typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
for (; it != simplices.end(); ++it)
- for (unsigned f=0; f < it->face_count(); ++f)
+ for (unsigned f=0; f < it->faceCount(); ++f)
if (simplices[it->faces[f]].value > it->value)
return false;