gsInterpolateMap.cpp

Annotated source file

Here is the full file examples/gsInterpolateMap.cpp. Clicking on a function or class name will lead you to its reference documentation.

#include <gismo.h>

using namespace gismo;

int main(int argc, char* argv[])
{
    index_t p = 3; // Degree
    index_t k = 3; // Number of interor knots (more knots --> better approximation)
    index_t s = 1000; // samples for plotting
    std::string func_name_x("x");
    std::string func_name_y("y");
    std::string func_name_z("z");
    index_t d = 3; //dimension

    gsCmdLine cmd("Input is coordinate functions X(x,y,z), Y(x,y,z). Z(x,y,z), or"
                  "X(x,y), Y(x,y). Z(x,y), or X(x), Y(x). Z(x). This is controlled by"
                  "the dimension parameter d. The parameters take values in the interval [0,1]");

    cmd.addInt("p", "degree", "this is the degree", p);
    cmd.addInt("k", "knots", "This is the number of interior knots", k);
    cmd.addInt("d", "dim", "this is the parametric dimension", d);
    cmd.addInt("s", "samples", "this is samples for paraview", s);
    cmd.addString("X", "f1", "The X-coordinate of the function", func_name_x);
    cmd.addString("Y", "f2", "The Y-coordinate of the function", func_name_y);
    cmd.addString("Z", "f3", "The Z-coordinate of the function", func_name_z);

    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }

    // Define a function R^d --> R^2
    gsFunctionExpr<> func(func_name_x, func_name_y, func_name_z, d);

    // Define a B-spline space of dimension 2
    // start, end, numInt, multiplicities at end-points
    gsKnotVector<> KV(0.0, 1.0, k, p+1);
    gsBasis<>::uPtr tBasis;

    switch (d) // static dispatch..
    {
    case 1:
        tBasis = gsBSplineBasis<>::make(KV);
        break;
    case 2:
        tBasis = memory::make_unique(new gsTensorBSplineBasis<2>(KV,KV));
        break;
    case 3:
        tBasis = memory::make_unique(new gsTensorBSplineBasis<3>(KV,KV,KV));
        break;
    default:
    {
        gsWarn<<"Dimension must be 1, 2 or 3.";
        return 0;
    }
    };

    gsInfo <<"We are going to interpolate "<< func <<"\n";

    // Support: 2x2 matrix, first column is the lower left
    // point of the parametric domain (eg. 0,0)
    // second column is the upper right corner of the support (eg. 1,1)
    // gsMatrix<> support = tBasis->support();

    // Points to interpolate at (Greville points):
    gsMatrix<> intGrid = tBasis->anchors();
    gsInfo <<"Int. grid dim: "<< intGrid.dim() <<"\n";

    // Evaluate f at the Greville points
    gsMatrix<> fValues = func.eval(intGrid);
    gsInfo <<"Function values dim: "<< fValues.dim() <<"\n";

    // Returns a geometry with basis = tBasis
    // and coefficients being
    // computed as the interpolant of \a funct
    gsGeometry<>::uPtr interpolant = tBasis->interpolateAtAnchors(fValues);

    gsInfo << "Result :"<< *interpolant <<"\n";

    // Save the result as an XML file
    gsFileData<> fd;
    gsMultiPatch<> mp;
    mp.addPatch(give(interpolant));
    fd << mp;
    std::string remark("Made by interpolation from function: ( ");
    for (index_t i = 0; i<3; ++i)
    {
        remark += func.expression(i);
        remark += ( i==2 ? " )" : " , ");
    }
    fd.addComment( remark );
    fd.save("interpolant_spline");
    gsInfo <<"Result saved as interpolant_spline.xml \n";

    // Produce Paraview file for the interpolant
    //gsWriteParaview( *interpolant, support, "interpolant", s );

    return 0;
}