gsIOUtils.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsDebug.h>
#include <gsHSplines/gsHDomainIterator.h>
#include <gsHSplines/gsHTensorBasis.h>
#include <gsIO/gsBase64.h>

namespace gismo {

template<class T>
void makeMesh(const gsBasis<T>& basis, gsMesh<T> & mesh, int n = 0)
{
    const unsigned d = basis.dim();

    typedef typename gsMesh<T>::VertexHandle Vertex;
    typename gsBasis<T>::domainIter domIter = basis.makeDomainIterator();

    // variables for iterating over a cube (element is a cube)
    const gsVector<unsigned> zeros = gsVector<unsigned>::Zero(d);
    const gsVector<unsigned> ones  = gsVector<unsigned>::Ones(d);
    gsVector<unsigned> cur;

    // maps integer representation of a vertex into pointer to the
    // vertex coordinates
    std::vector<Vertex> map(1ULL<<d);

    // neighbour[i] are integer representations of certain neighbours of
    // vertex i (i counts in lexicographics order over all vertices)
    std::vector<std::vector<unsigned> > neighbour(1ULL<<d,
                                                  std::vector<unsigned>() );

    cur.setZero(d);
    int counter = 0;
    do
    {
        // set neighbour
        for (unsigned dim = 0; dim < d; dim++)
        {
            if (cur(dim) == 0)
            {
                const unsigned tmp =  counter | (1<< dim) ;
                neighbour[counter].push_back(tmp);
            }
        }
        counter++;

    } while (nextCubePoint<gsVector<unsigned> >(cur, zeros, ones));

    gsVector<T> vertex(d);

    for (; domIter->good(); domIter->next())
    {
        const gsVector<T>& low = domIter->lowerCorner();
        const gsVector<T>& upp = domIter->upperCorner();
        const T vol = domIter->volume();

        vertex.setZero();
        cur.setZero();
        counter = 0;

        // add points to the mesh
        do
        {
            // get appropriate coordinate of a point
            for (unsigned dim = 0; dim < d; dim++)
            {
                vertex(dim) = ( cur(dim) ?  upp(dim) : low(dim) );
            }

            Vertex v = mesh.addVertex(vertex);
            v->data  = vol;
            map[counter++] = v;

        } while (nextCubePoint<gsVector<unsigned> >(cur, zeros, ones));


        // add edges to the mesh (connect points)
        for (size_t index = 0; index != neighbour.size(); index++)
        {
            const std::vector<unsigned> & v = neighbour[index];

            for (size_t ngh = 0; ngh != v.size(); ngh++)
            {
                // add more vertices (n) for better physical resolution
                mesh.addLine( map[index], map[v[ngh]], n );
                //mesh.addEdge( map[index], map[v[ngh]] );
            }
        }

        // idea: instead of edges add the faces to the mesh
        // mesh->addFace( mesh.vertices().back(),
        //                *(mesh.vertices().end()-3),
        //                *(mesh.vertices.end()-4),
        //                *(mesh.vertices.end()-2)
        //     );

    }
}

namespace internal
{

template <short_t d, typename T>
void makeHierarchicalMesh(const gsHTensorBasis<d, T>& basis,
                          std::vector<gsMesh<T> >& meshes,
                          int n = 0)
{
    // prepare meshes
    meshes.clear();
    for (unsigned i = 0; i < basis.maxLevel() + 1; i++)
    {
        meshes.push_back(gsMesh<T>());
    }

    // variables for iterating over a cube (element is a cube)
    const gsVector<unsigned> zeros = gsVector<unsigned>::Zero(d);
    const gsVector<unsigned> ones  = gsVector<unsigned>::Ones(d);
    gsVector<unsigned> cur;

    // neighbour[i] are integer representations of certain neighbours of
    // vertex i (i counts in lexicographics order over all vertices)
    std::vector<std::vector<unsigned> > neighbour(1 << d,
                                                  std::vector<unsigned>() );


    cur.setZero(d);
    int counter = 0;
    do
    {
        // set neighbour
        for (unsigned dim = 0; dim < d; dim++)
        {
            if (cur(dim) == 0)
            {
                const unsigned tmp =  counter | (1<< dim) ;
                neighbour[counter].push_back(tmp);
            }
        }
        counter++;

    } while (nextCubePoint<gsVector<unsigned> >(cur, zeros, ones));


    // maps integer representation of a vertex into pointer to the
    // vertex coordinates
    typedef typename gsMesh<T>::VertexHandle Vertex;
    std::vector<Vertex> map(1 << d);

    gsVector<T> vertex(d);

    gsHDomainIterator<T, d> domIter(basis);

    for (; domIter.good(); domIter.next())
    {

        int level = domIter.getLevel();

        const gsVector<T>& low = domIter.lowerCorner();
        const gsVector<T>& upp = domIter.upperCorner();


        vertex.setZero();
        cur.setZero();
        counter = 0;

        // add points to the mesh
        do
        {
            // get appropriate coordinate of a point
            for (unsigned dim = 0; dim < d; dim++)
            {
                vertex(dim) = ( cur(dim) ?  upp(dim) : low(dim) );
            }

            meshes[level].addVertex(vertex);
            map[counter++] = meshes[level].vertices().back();

        } while (nextCubePoint<gsVector<unsigned> >(cur, zeros, ones));


        // add edges to the mesh (connect points)
        for (size_t index = 0; index != neighbour.size(); index++)
        {
            const std::vector<unsigned> & v = neighbour[index];

            for (size_t ngh = 0; ngh != v.size(); ngh++)
            {
                // add more vertices (n) for better physical resolution
                meshes[level].addLine( map[index], map[v[ngh]], n );
            }
        }
    }
}

} // end namespace internal



template <typename T>
bool makeHierarchicalMesh(const gsBasis<T>& basis,
                          std::vector<gsMesh<T> >& meshes,
                          int n = 0)
{

    const gsHTensorBasis<1, T>* hBasis1 =
            dynamic_cast<const gsHTensorBasis<1, T>* >(&basis);

    if (hBasis1 != NULL)
    {
        internal::makeHierarchicalMesh<1, T>(*hBasis1, meshes, n);
        return true;
    }


    const gsHTensorBasis<2, T>* hBasis2 =
            dynamic_cast<const gsHTensorBasis<2, T>* >(&basis);

    if (hBasis2 != NULL)
    {
        internal::makeHierarchicalMesh<2, T>(*hBasis2, meshes, n);
        return true;
    }


    const gsHTensorBasis<3, T>* hBasis3 =
            dynamic_cast<const gsHTensorBasis<3, T>* >(&basis);

    if (hBasis3 != NULL)
    {
        internal::makeHierarchicalMesh<3, T>(*hBasis3, meshes, n);
        return true;
    }

    return false;
}

}  // namespace gismo