Computing Volumes in a DEC Mesh

Home    Publications    Blog    About Me   

18 Feb 2014

One of the more interesting blocks of code I wrote for my DEC implementation was generating the primal and dual volumes for computing the diagonal Hodge star. Through the pyDEC codebase I found a link to this which provides some nice insights for how to compute volumes for simplices.

Dual volumes are slightly more interesting because, on top of the dual mesh not being simplicial, the cells cross primal boundaries. The above method doesn’t directly work, so we have to tesselate the dual cells into simplices and accumulating the volumes of the simplices. One convenient and fairly nice tesselation can be generated by taking the sequences of circumcenters in the simplicial complex.

Let $n$ be the dimension of the simplicial complex and $m$ the dimension of the primal simplex $S^k$, dual to the cell we want to compute the volume of. Using $<$ be the partial ordering of simplices, consider all sequences $\{S^k_i\}_{i=m\ldots n}$ where $S^k_i$ is an $i$-simplex, $S^k_m=S^k$, and $\forall i, S^m_i < S^m_{i+1}$. If we let each set of circumcenters define a dual partial-simplex, $\bar S_j^k = \{\mathcal{C}(S_i^k)\}$, we get that the volume of the dual cell is $\sum_j vol(\bar S_j^k).$

We can order the evaluation of these dual cells to fill a single, statically sized, array of circumcenters and compute the volumes of the first $n-m$ entries of the array for the dual of a $m$-simplex.

template <int N,int TopDim,int EmbedDim> 
void genDualVolume(Simplex<N> & s, std::vector<Vec<EmbedDim> > & clist) {
    size_t denom = factorial(TopDim - N);
    clist[N] = s.center;
    if(N == TopDim) {
        s.dualVolume = 1;
    } else {
        int M = EmbedDim;
        Matrix m(M,clist.size()-N-1);
        auto&& origin = clist.back();
        for(size_t i=N; i < TopDim; ++i) {
            m.col(i-N) = clist[i] - origin;
        }
        s.dualVolume += std::sqrt((m.transpose()*m).determinant())/denom;
    }
    for(Simplex<N-1>& c: boundary(s)) {
        genDualVolume<N-1,TopDim,EmbedDim>(c, clist);
    }
}
template <int TopDim, int EmbedDim>
void populateDualVolumes() {
    //Initialize every simplex with their circumcenter
    std::vector<VecNf> clist(TopDim+1,Vec<EmbedDim>::Zero());
    for(Simplex<N>& s: m_simplices) {
        genDualVolume<TopDim,TopDim, EmbedDim>>(s, clist);
    }
}
comments powered by Disqus