gsFlowUtils.h

Go to the documentation of this file.

 
#pragma once
 
#include <gismo.h>
 
namespace gismo
{
 
#ifdef GISMO_WITH_PARDISO
template <class T>
inline void pardisoSetup(typename gsSparseSolver<T>::PardisoLU& solver)
{
    solver.setParam(7, 15);
    solver.setParam(9, 13);
    solver.setParam(12, 0);
}
#endif
 
 
inline void startAnimationFile(std::ofstream& file)
{
    file << "<?xml version=\"1.0\"?>\n";
    file << "<VTKFile type=\"Collection\" version=\"0.1\">";
    file << "<Collection>\n";
}
 
 
inline void endAnimationFile(std::ofstream& file)
{
    file << "</Collection>\n";
    file << "</VTKFile>\n";
    file.close();
}
 
 
inline void gsWriteOutput(std::ofstream& file, const std::string output, bool fileOutput, bool dispInTerminal)
{
    if (fileOutput)
        file << output;
 
    if (dispInTerminal)
        gsInfo << output;
}
 
 
inline void gsWriteOutputLine(std::ofstream& file, const std::string line, bool fileOutput, bool dispInTerminal)
{
    gsWriteOutput(file, line, fileOutput, dispInTerminal);
 
    if (fileOutput)
        file << std::endl;
 
    if (dispInTerminal)
            gsInfo << std::endl;
}
 
 
template<class T>
gsMatrix<T> createVectorOfUniqueIndices(const gsMatrix<T>& mat)
{
    std::map<T, bool> map;
    std::vector<T> vecTmp;
 
    for (index_t j = 0; j < mat.cols(); j++)
    {
        for (index_t i = 0; i < mat.rows(); i++)
        {
            if(map.count(mat(i,j)) == 0)
            {
                map[mat(i,j)] = true;
                vecTmp.push_back(mat(i,j));
            }
        }
    }
 
    gsMatrix<T> vec(vecTmp.size(), 1);
    for (size_t i = 0; i < vecTmp.size(); i++)
        vec(i) = vecTmp[i];
 
 
    return vec;
}
 
 
template<class T, int MatOrder>
gsVector<index_t> getNnzVectorPerOuter(const gsSparseMatrix<T, MatOrder>& mat)
{
    gsVector<index_t> nnzPerOuter(mat.outerSize());
    nnzPerOuter.setZero();
 
    for (index_t outer = 0; outer < mat.outerSize(); outer++)
        for (typename gsSparseMatrix<T, MatOrder>::InnerIterator it(mat, outer); it; ++it)
            ++nnzPerOuter[outer];
 
    return nnzPerOuter;
}
 
 
template<class T, int MatOrder>
void diagInvMatrix_into(const gsSparseMatrix<T, MatOrder>& mat, gsSparseMatrix<T, MatOrder>& diagInv, int repeat, bool lumping = false)
{
    GISMO_ENSURE(mat.nonZeros() != 0, "diagInvMatrix_into(): The matrix is empty!");
 
    int varDofs = mat.rows();
 
    diagInv.resize(repeat*varDofs, repeat*varDofs);
    diagInv.reserve(gsVector<int>::Constant(diagInv.cols(), 1));
 
    const gsSparseMatrix<T, MatOrder>* matPtr = &mat;
    gsSparseMatrix<T, MatOrder> lumped(varDofs, varDofs);
 
    if (lumping)
    {
        lumped.reserve(gsVector<int>::Constant(varDofs, 1));
 
        for (int j = 0; j < varDofs; j++)
            lumped.insert(j, j) = math::abs(mat.at(j, j)); // abs value because of "lumping" diag(A) in SIMPLE-type prec., does not change lumped mass matrix in IgA
 
        for (int j = 0; j < varDofs; j++)
        {
            for (typename gsSparseMatrix<T, MatOrder>::InnerIterator it(mat, j); it; ++it)
            {
                int i = it.row();
 
                if (i != j)
                    lumped.coeffRef(i, i) += math::abs(it.value());
            }
        }
 
        matPtr = &lumped;
    }
 
    for (int i = 0; i < varDofs; i++)
    {
        T tmp = 1 / matPtr->coeff(i, i);
 
        for (int s = 0; s < repeat; s++)
            diagInv.coeffRef(i + s * varDofs, i + s * varDofs) = tmp;
    }
}
 
 
template<class T>
gsTensorBSpline<2, T> BSplineRectangle(int deg, const T llx, const T lly, const T a, const T b, int numSep = 0)
{
    gsKnotVector<T> kv(0, 1, numSep, deg + 1, deg); // first, last, num_inter, mult_end, mult_inter
 
    int n = kv.size() - deg - 1;
    gsMatrix<T> coef(n*n, 2);
 
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            coef.row(i + j*n) << llx + (i*a) / (deg*(numSep + 1)), lly + (j*b) / (deg*(numSep + 1));
 
    return gsTensorBSpline<2, T>(kv, kv, coef);
}
 
 
template<class T>
gsTensorBSpline<3, T> BSplineBlock(int deg, const T llx, const T lly, const T llz, const T a, const T b, const T c, int numSep = 0)
{
    gsKnotVector<T> kv(0, 1, numSep, deg + 1, deg);
 
    int n = kv.size() - deg - 1;
    gsMatrix<T> coef(n*n*n, 3);
 
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                coef.row(i + j*n + k*n*n) << llx + (i*a) / (deg*(numSep + 1)), lly + (j*b) / (deg*(numSep + 1)), llz + (k*c) / (deg*(numSep + 1));
 
    return gsTensorBSpline<3, T>(kv, kv, kv, coef);
}
 
 
template<class T>
gsMultiPatch<T> BSplineCavity2D(int deg, const T a, const T b, const int np = 1, int numSep = 0)
{
    gsMultiPatch<T> mp;
 
    T aPatch = a / np;
    T bPatch = b / np;
 
    for (int j = 0; j < np; j++)
        for (int i = 0; i < np; i++)
            mp.addPatch(BSplineRectangle(deg, i*aPatch, j*bPatch, aPatch, bPatch, numSep));
 
    mp.computeTopology();
 
    return mp;
}
 
 
template<class T>
gsMultiPatch<T> BSplineCavity3D(int deg, const T a, const T b, const T c, int numSep = 0)
{
    gsMultiPatch<T> mp;
    mp.addPatch(BSplineBlock(deg, 0.0, 0.0, 0.0, a, b, c, numSep));
    mp.computeTopology();
    return mp;
}
 
 
template<class T>
gsMultiPatch<T> BSplineStep2D(int deg, const T a, const T b, const T a_in, T h = 0, bool periodic = false)
{
    gsMultiPatch<T> mp;
 
    if (h == 0)
        h = b / 2;
 
    mp.addPatch(BSplineRectangle(deg, 0.0, 0.0, a, h));
    mp.addPatch(BSplineRectangle(deg, 0.0, h, a, b - h));
    mp.addPatch(BSplineRectangle(deg, -a_in, h, a_in, b - h));
 
    mp.addInterface(0, boundary::north, 1, boundary::south);
    mp.addInterface(1, boundary::west, 2, boundary::east);
 
    if(periodic)
        mp.addInterface(0, boundary::south, 1, boundary::north);
 
    mp.addAutoBoundaries();
 
    return mp;
}
 
 
template<class T>
gsMultiPatch<T> BSplineStep3D(int deg, const T a, const T b, const T c, const T a_in, T h = 0, bool periodic = false)
{
    gsMultiPatch<T> mp;
 
    if (h == 0)
        h = b / 2;
 
    mp.addPatch(BSplineBlock<T>(deg, 0.0, 0.0, 0.0, a, h, c));
    mp.addPatch(BSplineBlock<T>(deg, 0.0, h, 0.0, a, b - h, c));
    mp.addPatch(BSplineBlock<T>(deg, -a_in, h, 0.0, a_in, b - h, c));
 
    mp.addInterface(0, boundary::north, 1, boundary::south);
    mp.addInterface(2, boundary::east, 1, boundary::west);
 
    if (periodic)
    {
        mp.addInterface(0, boundary::front, (size_t)0, boundary::back);
        mp.addInterface(1, boundary::front, 1, boundary::back);
        mp.addInterface(2, boundary::front, 2, boundary::back);
    }
    
    mp.addAutoBoundaries();
 
    return mp;
}
 
 
template<class T>
void addBCs(gsBoundaryConditions<T>& bcInfo, std::vector<std::pair<int, boxSide> >& bndIn, std::vector<std::pair<int, boxSide> >& bndWall, gsFunctionExpr<T> Uin, gsFunctionExpr<T> Uwall)
{
    for (size_t i = 0; i < bndIn.size(); i++)
        bcInfo.addCondition(bndIn[i].first, bndIn[i].second, condition_type::dirichlet, Uin);
 
    for (size_t i = 0; i < bndWall.size(); i++)
        bcInfo.addCondition(bndWall[i].first, bndWall[i].second, condition_type::dirichlet, Uwall);
}
 
 
template<class T>
void defineBCs_cavity2D(gsBoundaryConditions<T>& bcInfo, const int np, std::vector<std::pair<int, boxSide> >& bndWall, std::string lidVel = "1")
{
    gsFunctionExpr<T> Uwall("0", "0", 2);
    gsFunctionExpr<T> Ulid(lidVel, "0", 2);
 
    for (int i = 1; i <= np; i++)
    {
        bcInfo.addCondition(np*np - i, boundary::north, condition_type::dirichlet, Ulid, 0);
        bndWall.push_back(std::make_pair(np*np - i, boundary::north));
    }
 
    for (int i = 0; i < np; i++)
    {
        bcInfo.addCondition(i, boundary::south, condition_type::dirichlet, Uwall, 0);
        bndWall.push_back(std::make_pair(i, boundary::south));
    }
 
    for (int i = 0; i < np; i++)
    {
        bcInfo.addCondition(i * np, boundary::west, condition_type::dirichlet, Uwall, 0);
        bcInfo.addCondition((i + 1)*np - 1, boundary::east, condition_type::dirichlet, Uwall, 0);
        bndWall.push_back(std::make_pair(i * np, boundary::west));
        bndWall.push_back(std::make_pair((i + 1)*np - 1, boundary::east));
    }
}
 
 
template<class T>
void defineBCs_cavity3D(gsBoundaryConditions<T>& bcInfo, std::vector<std::pair<int, boxSide> >& bndWall, std::string lidVel = "1")
{
    gsFunctionExpr<T> Uwall("0", "0", "0", 3);
    gsFunctionExpr<T> Ulid(lidVel, "0", "0", 3);
 
    bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, Ulid, 0);
    bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, Uwall, 0);
    bcInfo.addCondition(0, boundary::west, condition_type::dirichlet, Uwall, 0);
    bcInfo.addCondition(0, boundary::east, condition_type::dirichlet, Uwall, 0);
    bcInfo.addCondition(0, boundary::front, condition_type::dirichlet, Uwall, 0);
    bcInfo.addCondition(0, boundary::back, condition_type::dirichlet, Uwall, 0);
    bndWall.push_back(std::make_pair(0, boundary::north));
    bndWall.push_back(std::make_pair(0, boundary::south));
    bndWall.push_back(std::make_pair(0, boundary::west));
    bndWall.push_back(std::make_pair(0, boundary::east));
    bndWall.push_back(std::make_pair(0, boundary::front));
    bndWall.push_back(std::make_pair(0, boundary::back));
}
 
 
template<class T>
void defineBCs_cavity(gsBoundaryConditions<T>& bcInfo, std::vector<std::pair<int, boxSide> >& bndWall, int dim, const int np = 1, std::string lidVel = "1")
{
    switch (dim)
    {
    case 2:
        defineBCs_cavity2D(bcInfo, np, bndWall, lidVel);
        break;
    case 3:
        defineBCs_cavity3D(bcInfo, bndWall, lidVel);
        break;
    default:
        GISMO_ERROR("Wrong dimension!");
        break;
    }
}
 
 
template<class T>
void defineBCs_step2D(gsBoundaryConditions<T>& bcInfo, std::vector<std::pair<int, boxSide> >& bndIn, std::vector<std::pair<int, boxSide> >& bndOut, std::vector<std::pair<int, boxSide> >& bndWall, bool periodic = false, std::string inVel = "default")
{
    if (inVel == "default")
        inVel = "(-4*(y-1.5)^2 + 1)";
 
    gsFunctionExpr<T> Uin, Uwall;
    Uin = gsFunctionExpr<T>(inVel, "0", 2);
    Uwall = gsFunctionExpr<T>("0", "0", 2);
 
    if (!periodic)
    {
        bndIn.push_back(std::make_pair(2, boundary::west));
        bndWall.push_back(std::make_pair(0, boundary::west));
        bndWall.push_back(std::make_pair(0, boundary::south));
        bndWall.push_back(std::make_pair(1, boundary::north));
        bndWall.push_back(std::make_pair(2, boundary::south));
        bndWall.push_back(std::make_pair(2, boundary::north));
        bndOut.push_back(std::make_pair(0, boundary::east));
        bndOut.push_back(std::make_pair(1, boundary::east));
    }
    else
    {
        bndIn.push_back(std::make_pair(2, boundary::west));
        bndWall.push_back(std::make_pair(0, boundary::west));
        bndWall.push_back(std::make_pair(2, boundary::south));
        bndWall.push_back(std::make_pair(2, boundary::north));
        bndOut.push_back(std::make_pair(0, boundary::east));
        bndOut.push_back(std::make_pair(1, boundary::east));
    }
 
    addBCs(bcInfo, bndIn, bndWall, Uin, Uwall);
}
 
 
template<class T>
void defineBCs_step3D(gsBoundaryConditions<T>& bcInfo, std::vector<std::pair<int, boxSide> >& bndIn, std::vector<std::pair<int, boxSide> >& bndOut, std::vector<std::pair<int, boxSide> >& bndWall, bool periodic = false, std::string inVel = "default")
{
    gsFunctionExpr<T> Uin, Uwall;
 
    if (!periodic)
    {
        if (inVel == "default")
            inVel = "(-4*(y-1.5)^2 + 1)*(-(z-1)^2 + 1)";
 
        Uin = gsFunctionExpr<T>(inVel, "0", "0", 3);
        Uwall = gsFunctionExpr<T>("0", "0", "0", 3);
 
        bndIn.push_back(std::make_pair(2, boundary::west));
        bndWall.push_back(std::make_pair(0, boundary::west));
        bndWall.push_back(std::make_pair(0, boundary::south));
        bndWall.push_back(std::make_pair(0, boundary::front));
        bndWall.push_back(std::make_pair(0, boundary::back));
        bndWall.push_back(std::make_pair(1, boundary::north));
        bndWall.push_back(std::make_pair(1, boundary::front));
        bndWall.push_back(std::make_pair(1, boundary::back));
        bndWall.push_back(std::make_pair(2, boundary::south));
        bndWall.push_back(std::make_pair(2, boundary::north));
        bndWall.push_back(std::make_pair(2, boundary::front));
        bndWall.push_back(std::make_pair(2, boundary::back));
        bndOut.push_back(std::make_pair(0, boundary::east));
        bndOut.push_back(std::make_pair(1, boundary::east));
    }
    else
    {
        if (inVel == "default")
            inVel = "(-4*(y-1.5)^2 + 1)";
 
        Uin = gsFunctionExpr<T>(inVel, "0", "0", 3);
        Uwall = gsFunctionExpr<T>("0", "0", "0", 3);
 
        bndIn.push_back(std::make_pair(2, boundary::west));
        bndWall.push_back(std::make_pair(0, boundary::west));
        bndWall.push_back(std::make_pair(0, boundary::south));
        bndWall.push_back(std::make_pair(1, boundary::north));
        bndWall.push_back(std::make_pair(2, boundary::south));
        bndWall.push_back(std::make_pair(2, boundary::north));
        bndOut.push_back(std::make_pair(0, boundary::east));
        bndOut.push_back(std::make_pair(1, boundary::east));
    }
 
    addBCs(bcInfo, bndIn, bndWall, Uin, Uwall);
}
 
 
template<class T>
void defineBCs_step(gsBoundaryConditions<T>& bcInfo, std::vector<std::pair<int, boxSide> >& bndIn, std::vector<std::pair<int, boxSide> >& bndOut, std::vector<std::pair<int, boxSide> >& bndWall, int dim, bool periodic = false, std::string inVel = "default")
{
    switch (dim)
    {
    case 2:
        defineBCs_step2D(bcInfo, bndIn, bndOut, bndWall, periodic, inVel);
        break;
    case 3:
        defineBCs_step3D(bcInfo, bndIn, bndOut, bndWall, periodic, inVel);
        break;
    default:
        GISMO_ERROR("Wrong dimension!");
        break;
    }
}
 
 
template<class T>
void defineBCs_profile2D(gsBoundaryConditions<T>& bcInfo, std::vector<std::pair<int, boxSide> >& bndIn, std::vector<std::pair<int, boxSide> >& bndOut, std::vector<std::pair<int, boxSide> >& bndWall, T inVelX, T inVelY)
{
    gsFunctionExpr<T> Uin = gsFunctionExpr<T>(util::to_string(inVelX), util::to_string(inVelY), 2);
    gsFunctionExpr<T> Uwall = gsFunctionExpr<T>("0", "0", 2);
 
    bndIn.push_back(std::make_pair(0, boundary::west));
    bndWall.push_back(std::make_pair(1, boundary::north));
    bndWall.push_back(std::make_pair(1, boundary::south));
    bndOut.push_back(std::make_pair(2, boundary::east));
 
    addBCs(bcInfo, bndIn, bndWall, Uin, Uwall);
}
 
 
template< int d, class T>
void refineFirstKnotSpan(gsMultiBasis<T>& basis, int numRefine, int patch, int dir)
{
    const gsTensorBSplineBasis<d, T>*  patchBasis = dynamic_cast<const gsTensorBSplineBasis<d, T>*>(&(basis.basis(patch)));
    gsMatrix<T> box(d, 2);
 
    for (int i = 0; i < numRefine; i++)
    {
        box.setZero();
        T knot = patchBasis->knot(dir, patchBasis->degree(dir) + 1);
        box(dir, 1) = knot;
        basis.refine(patch, box);
    }
}
 
 
template< int d, class T>
void refineLastKnotSpan(gsMultiBasis<T>& basis, int numRefine, int patch, int dir)
{
    const gsTensorBSplineBasis<d, T>*  patchBasis = dynamic_cast<const gsTensorBSplineBasis<d, T>*>(&(basis.basis(patch)));
    gsMatrix<T> box(d, 2);
 
    for (int i = 0; i < numRefine; i++)
    {
        box.setZero();
        int sizeKnots = patchBasis->knots(dir).size() - 1;
        T lastKnot = patchBasis->knot(dir, sizeKnots);
        T knot = patchBasis->knot(dir, sizeKnots - (patchBasis->degree(dir) + 1));
        box(dir, 0) = knot;
        box(dir, 1) = lastKnot;
        basis.refine(patch, box);
    }
}
 
 
template<class T>
void refineBasis_cavity(gsMultiBasis<T>& basis, int numRefine, int numRefineLocal, int dim, int numSep = 0)
{
    switch (dim)
    {
        case 2:
        {
            int npDir = math::sqrt(basis.nPieces());
            int npRef = std::log2(npDir);
            int sepRef = std::log2(numSep + 1);
 
            for (int i = 0; i < numRefine - npRef - sepRef; ++i)
                basis.uniformRefine();
 
            if (numRefineLocal && npDir == 1)
            {
                for (int dir = 0; dir < dim; dir++)
                {
                    refineFirstKnotSpan<2, T>(basis, numRefineLocal, 0, dir);
                    refineLastKnotSpan<2, T>(basis, numRefineLocal, 0, dir);
                }
            }
 
            break;
        }
        case 3:
        {
            for (int i = 0; i < numRefine; ++i)
                basis.uniformRefine();
 
            for (int dir = 0; dir < dim; dir++)
            {
                refineFirstKnotSpan<3, T>(basis, numRefineLocal, 0, dir);
                refineLastKnotSpan<3, T>(basis, numRefineLocal, 0, dir);
            }
 
            break;
        }
        default:
            GISMO_ERROR("Wrong dimension!");
            break;
    }
}
 
 
template<int d, class T>
void refineLocal_step(gsMultiBasis<T>& basis, int numRefineWalls, int numRefineCorner)
{
    // refinement near upper and bottom walls
    refineFirstKnotSpan<d, T>(basis, numRefineWalls, 0, 1);
    refineLastKnotSpan<d, T>(basis, numRefineWalls, 1, 1);
    refineLastKnotSpan<d, T>(basis, numRefineWalls, 2, 1);
 
    // refinement near the corner 
    refineFirstKnotSpan<d, T>(basis, numRefineCorner, 0, 0);
    refineFirstKnotSpan<d, T>(basis, numRefineCorner, 1, 0);
    refineLastKnotSpan<d, T>(basis, numRefineCorner, 2, 0);
    refineLastKnotSpan<d, T>(basis, numRefineCorner, 0, 1);
    refineFirstKnotSpan<d, T>(basis, numRefineCorner, 1, 1);
    refineFirstKnotSpan<d, T>(basis, numRefineCorner, 2, 1);
}
 
 
template<class T>
void refineBasis_step(gsMultiBasis<T>& basis, int numRefine, int numRefineWalls, int numRefineCorner, int numRefineU, real_t addRefPart, int dim, real_t a, real_t b, real_t c = 0.0)
{
    for (int i = 0; i < numRefine; ++i)
        basis.uniformRefine();
 
    gsMatrix<T> box(dim, 2);
 
    int uRefine = math::floor(std::log2(a / b)) + 1 + numRefineU;
    box.setZero();
    box(0, 1) = 1; 
    for (int i = 0; i < uRefine; i++)
        for (int p = 0; p < 2; p++)
            basis.refine(p, box);
 
    if (dim == 3)
    {
        int wRefine = math::floor(std::log2(c / b)) + 1;
        box.setZero();
        box(2, 1) = 1;
        for (int i = 0; i < wRefine; i++)
            for (int p = 0; p < 3; p++)
                basis.refine(p, box);
    }
 
    box.setZero();
    box(0,1) = addRefPart;
    for (int p = 0; p < 2; p++)
        basis.refine(p, box);
 
    switch (dim)
    {
    case 2:
        refineLocal_step<2, T>(basis, numRefineWalls, numRefineCorner);
        break;
    case 3:
        refineLocal_step<3, T>(basis, numRefineWalls, numRefineCorner);
        break;
    default:
        GISMO_ERROR("Wrong dimension!");
        break;
    }
}
 
 
template<class T>
void refineBasis_profile2D(gsMultiBasis<T>& basis, int numRefine, int numRefineBlade, int numRefineLead)
{
    for (int i = 0; i < numRefine; ++i)
        basis.uniformRefine();
 
    if (numRefineBlade)
    {
        for (int p = 0; p < basis.nPieces(); p++)
        {
            refineFirstKnotSpan<2, T>(basis, numRefineBlade, p, 1);
            refineLastKnotSpan<2, T>(basis, numRefineBlade, p, 1);
        }
    }
 
    if (numRefineLead)
    {
        refineLastKnotSpan<2, T>(basis, numRefineLead, 0, 0);
        refineFirstKnotSpan<2, T>(basis, numRefineLead, 1, 0);
    }
}
 
 
template<class T>
void plotQuantityFromSolution(std::string quantity, const gsField<T>& solField, std::string const & fn, unsigned npts)
{
    gsParaviewCollection pwCollection(fn);
    std::string fileName;
 
    for (unsigned int p = 0; p < solField.patches().nPatches(); p++)
    {
        fileName = fn + util::to_string(p);
 
        gsMatrix<T> geoVals, quantityVals;
        gsVector<unsigned> np;
 
        gsGeometry<T>* patch = &solField.patches().patch(p);
        const int parDim = patch->domainDim();
        const int tarDim = patch->targetDim();
 
        gsMatrix<T> ab = patch->support();
        gsVector<T> a = ab.col(0);
        gsVector<T> b = ab.col(1);
        np = uniformSampleCount(a, b, npts);
        gsMatrix<T> pts = gsPointGrid(a, b, np);
 
        // evaluate geometry at pts
        geoVals = patch->eval(pts);
 
        if (3 - parDim > 0)
        {
            np.conservativeResize(3);
            np.bottomRows(3 - parDim).setOnes();
        }
        if (3 - tarDim > 0)
        {
            geoVals.conservativeResize(3, geoVals.cols());
            geoVals.bottomRows(3 - tarDim).setZero();
        }
 
        // evaluate quantity at pts
        if (quantity == "divergence")
            solField.function(p).div_into(pts, quantityVals); // incorrect?? (derivatives in parametric space?)
        else
            GISMO_ERROR("Function plotQuantityFromSolution for " + quantity + " not implemented.");
 
        // paraview export
        gsWriteParaviewTPgrid(geoVals, quantityVals, np.template cast<index_t>(), fileName);
 
        pwCollection.addPart(fileName + ".vts");
    }
    pwCollection.save();
}
 
 
// -----------------------------------------------
 
 
template<class T>
struct gsVectorHash
{
    std::size_t operator()(const gsVector<T>& vec) const
    {
        std::size_t seed = vec.size();
        for (int i = 0; i < vec.size(); i++)
            seed ^= std::hash<T>()(vec(i)) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
 
        return seed;
    }
};
 
 
inline index_t mapTensorID(index_t i, index_t j, index_t k, size_t size1, size_t size2)
{
    index_t result = (k*size1*size2)+(j*size1)+i;
    return result;
}
 
 
} //namespace gismo