geometry_example.cpp

In this example we show how to read a geometry from a file, and how to analyze its properties and perform some basic operations.

First we include the library and use the gismo namespace.

#include <iostream>
#include <gismo.h>
 
using namespace gismo;

Then we define our command line options. For example, we use the option -f to set the path to the file that contains our geometry.

    std::string input("surfaces/simple.xml");
    std::string output("");
 
    gsCmdLine cmd("Tutorial on gsGeometry class.");
    cmd.addString("f", "filename", "G+Smo input geometry file.", input);
    cmd.addString("o", "output", "Name of the output file", output);
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }

Then we read the input from the file that contains the geometry. In this example, we use the function gsFileData::getFirst() to read the first gsGeometry contained in the input xml file.

    gsFileData<> fileData(input);
 
    gsGeometry<>::uPtr pGeom;
    if (fileData.has< gsGeometry<> >())
    {
        pGeom = fileData.getFirst< gsGeometry<> >();
    }
    else
    {
        gsInfo << "Input file doesn't have a geometry inside." << "\n";
        return -1;
    }
 
 
    if (!pGeom)
    {
        gsInfo << "Didn't find any geometry." << "\n";
        return -1;
    }

We can then video-print the the geometry, together with its basis and coefficients:

    gsInfo << "The file contains: \n" << *pGeom << "\n";
 
    // G+Smo geometries contains basis and coefficients
    const gsBasis<>& basis = pGeom->basis();
    gsInfo << "\nBasis: \n" << basis << "\n";
 
    const gsMatrix<>& coefs = pGeom->coefs();
    gsInfo << "\nCoefficients: \n" << coefs << "\n" << "\n";

as well its properties. Among others, we show here how to print the dimenstion of the parameter space, the dimension of the geometry, and the geometry support.

    gsInfo << "Dimension of the parameter space: " << pGeom->parDim() << "\n"
              << "Dimension of the geometry: " << pGeom->geoDim() << "\n";
 
    // support of the geometry, this is the same as gsBasis::support
    // (dim x 2 matrix, the parametric domain)
    gsMatrix<> support = pGeom->support();
    gsInfo << "Support: \n"
              << support << "\n" << "\n";

Moreover, we can evaluate the geometry and its 1st and 2nd derivaties at \(N\) given points \(\mathbf{u}_i\in\mathbb{R}^n\) for \(i = 1, \ldots, N\), and print the results. Let \(f: \mathbb{R}^n \to \mathbb{R}^m\), and \(f^{(j)}\) be the \(j\)-th component of function \(f\), for \(j = 1, \ldots, m\).

The function evaluation format is:

\[ \left[ \begin{array}{ccccc} f^{(1)}(\mathbf{u}_1) & f^{(1)}(\mathbf{u}_2) & \ldots & f^{(1)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots\\ f^{(m)}(\mathbf{u}_1) & f^{(m)}(\mathbf{u}_2) & \ldots & f^{(m)}(\mathbf{u}_N)\\ \end{array} \right]. \]

The first derivative format is:

\[ \left[ \begin{array}{ccccc} \partial_{1}f^{(1)}(\mathbf{u}_1) & \partial_{1}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{1}f^{(1)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots\\ \partial_{n}f^{(1)}(\mathbf{u}_1) & \partial_{n}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{n}f^{(1)}(\mathbf{u}_N)\\ \partial_{1}f^{(2)}(\mathbf{u}_1) & \partial_{1}f^{(2)}(\mathbf{u}_2) & \ldots & \partial_{1}f^{(2)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots\\ \partial_{n}f^{(2)}(\mathbf{u}_1) & \partial_{n}f^{(2)}(\mathbf{u}_2) & \ldots & \partial_{n}f^{(2)}(\mathbf{u}_N)\\ \partial_{1}f^{(m)}(\mathbf{u}_1) & \partial_{1}f^{(m)}(\mathbf{u}_2) & \ldots & \partial_{1}f^{(m)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots\\ \partial_{n}f^{(m)}(\mathbf{u}_1) & \partial_{n}f^{(m)}(\mathbf{u}_2) & \ldots & \partial_{n}f^{(m)}(\mathbf{u}_N) \end{array} \right]. \]

The second derivative format is:

\[ \left[ \begin{array}{ccccc} \partial_{1}\partial_{1}f^{(1)}(\mathbf{u}_1) & \partial_{1}\partial_{1}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{1}\partial_{1}f^{(1)}(\mathbf{u}_N)\\ \partial_{2}\partial_{2}f^{(1)}(\mathbf{u}_1) & \partial_{2}\partial_{2}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{2}\partial_{2}f^{(1)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots \\ \partial_{n}\partial_{n}f^{(1)}(\mathbf{u}_1) & \partial_{n}\partial_{n}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{n}\partial_{n}f^{(1)}(\mathbf{u}_N)\\ \partial_{1}\partial_{2}f^{(1)}(\mathbf{u}_1) & \partial_{1}\partial_{2}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{1}\partial_{2}f^{(1)}(\mathbf{u}_N)\\ \partial_{1}\partial_{3}f^{(1)}(\mathbf{u}_1) & \partial_{1}\partial_{3}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{1}\partial_{3}f^{(1)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots \\ \partial_{1}\partial_{n}f^{(1)}(\mathbf{u}_1) & \partial_{1}\partial_{n}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{1}\partial_{n}f^{(1)}(\mathbf{u}_N)\\ \partial_{2}\partial_{3}f^{(1)}(\mathbf{u}_1) & \partial_{2}\partial_{3}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{2}\partial_{3}f^{(1)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots\\ \partial_{2}\partial_{n}f^{(1)}(\mathbf{u}_1) & \partial_{2}\partial_{n}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{2}\partial_{n}f^{(1)}(\mathbf{u}_N)\\ \partial_{3}\partial_{4}f^{(1)}(\mathbf{u}_1) & \partial_{3}\partial_{4}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{3}\partial_{4}f^{(1)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots \\ \partial_{n-1}\partial_{n}f^{(1)}(\mathbf{u}_1) & \partial_{n-1}\partial_{n}f^{(1)}(\mathbf{u}_2) & \ldots & \partial_{n-1}\partial_{n}f^{(1)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots \\ \partial_{1}\partial_{1}f^{(m)}(\mathbf{u}_1) & \partial_{1}\partial_{1}f^{(m)}(\mathbf{u}_2) & \ldots & \partial_{1}\partial_{1}f^{(m)}(\mathbf{u}_N)\\ \vdots & \vdots & & \vdots \\ \partial_{n-1}\partial_{n}f^{(m)}(\mathbf{u}_1) & \partial_{n-1}\partial_{n}f^{(m)}(\mathbf{u}_2) & \ldots & \partial_{n-1}\partial_{n}f^{(m)}(\mathbf{u}_N)\\ \end{array} \right]. \]

    gsMatrix<> u = 0.3 * support.col(0) + 0.7 * support.col(1);
    gsInfo << "u " << size(u) << ": \n" << u << "\n" << "\n";
 
    // geoDim x 1 matrix
    gsMatrix<> value = pGeom->eval(u);
    // geoDim * parDim x 1 matrix (columns represent gradients)
    gsMatrix<> der1 = pGeom->deriv(u);
    // [geoDim * (parDim + parDim * (parDim - 1) / 2)] x 1 matrix
    gsMatrix<> der2 = pGeom->deriv2(u);
 
    gsInfo << "Value at u " << size(value) << ": \n"
              << value
              << "\n\nDerivative at u " << size(der1) << ": \n"
              << der1
              << "\n\nSecond derivative at u " << size(der2) << ": \n"
              << der2
              << "\n" << "\n";

Finally, we can export tesselations (meshes) as files for visualization in Paraview

    // mesh holds the control net of a geometry
    // mesh is a set of vertices and lines (connections between vertices)
    gsMesh<> mesh;
    pGeom->controlNet(mesh);

and print the files in Paraview.

        std::string out = output + "Geometry";
        gsInfo << "Writing the geometry to a paraview file: " << out
                  << "\n\n";
 
        gsWriteParaview(*pGeom, out);
 
        out = output + "Basis";
        gsInfo << "Writing the basis to a paraview file: " << out
                  << "\n\n";
 
        gsWriteParaview(basis, out);
 
        out = output + "ContolNet";
        gsInfo << "Writing the control net to a paraview file: " << out
                  << "\n" << "\n";
 
        gsWriteParaview(mesh, out);

Annotated source file

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

 
#include <iostream>
#include <gismo.h>
 
using namespace gismo;
 
 
// Returns the string with the size of a matrix.
template <typename T>
std::string size(const gsMatrix<T>& matrix)
{
    std::string result = "(" + util::to_string(matrix.rows()) + " x " +
        util::to_string(matrix.cols()) + ")";
 
    return result;
}
 
 
int main(int argc, char* argv[])
{
    std::string input("surfaces/simple.xml");
    std::string output("");
 
    gsCmdLine cmd("Tutorial on gsGeometry class.");
    cmd.addString("f", "filename", "G+Smo input geometry file.", input);
    cmd.addString("o", "output", "Name of the output file", output);
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
    
    
    // ======================================================================
    // reading the geometry
    // ======================================================================
    gsFileData<> fileData(input);
 
    gsGeometry<>::uPtr pGeom;
    if (fileData.has< gsGeometry<> >())
    {
        pGeom = fileData.getFirst< gsGeometry<> >();
    }
    else
    {
        gsInfo << "Input file doesn't have a geometry inside." << "\n";
        return -1;
    }
 
 
    if (!pGeom)
    {
        gsInfo << "Didn't find any geometry." << "\n";
        return -1;
    }
 
    // ======================================================================
    // printing some information about the basis
    // ======================================================================
 
 
    // ----------------------------------------------------------------------
    // printing the geometry
    // ----------------------------------------------------------------------
 
    gsInfo << "The file contains: \n" << *pGeom << "\n";
 
    // G+Smo geometries contains basis and coefficients
    const gsBasis<>& basis = pGeom->basis();
    gsInfo << "\nBasis: \n" << basis << "\n";
 
    const gsMatrix<>& coefs = pGeom->coefs();
    gsInfo << "\nCoefficients: \n" << coefs << "\n" << "\n";
 
 
    // ----------------------------------------------------------------------
    // printing some properties about the basis
    // ----------------------------------------------------------------------
 
    gsInfo << "Dimension of the parameter space: " << pGeom->parDim() << "\n"
              << "Dimension of the geometry: " << pGeom->geoDim() << "\n";
 
    // support of the geometry, this is the same as gsBasis::support
    // (dim x 2 matrix, the parametric domain)
    gsMatrix<> support = pGeom->support();
    gsInfo << "Support: \n"
              << support << "\n" << "\n";
 
 
    // ======================================================================
    // evaluation
    // ======================================================================
 
 
    // ----------------------------------------------------------------------
    // values, 1st derivatives, and 2nd derivatives
    // ----------------------------------------------------------------------
 
    gsMatrix<> u = 0.3 * support.col(0) + 0.7 * support.col(1);
    gsInfo << "u " << size(u) << ": \n" << u << "\n" << "\n";
 
    // geoDim x 1 matrix
    gsMatrix<> value = pGeom->eval(u);
    // geoDim * parDim x 1 matrix (columns represent gradients)
    gsMatrix<> der1 = pGeom->deriv(u);
    // [geoDim * (parDim + parDim * (parDim - 1) / 2)] x 1 matrix
    gsMatrix<> der2 = pGeom->deriv2(u);
 
    gsInfo << "Value at u " << size(value) << ": \n"
              << value
              << "\n\nDerivative at u " << size(der1) << ": \n"
              << der1
              << "\n\nSecond derivative at u " << size(der2) << ": \n"
              << der2
              << "\n" << "\n";
 
    gsInfo << "\nFor more information about evaluation "
              << "(and order of derivatives) look at doxygen documentation."
              << "\n" << "\n";
 
 
    // ======================================================================
    // contol net
    // ======================================================================
 
    // mesh holds the control net of a geometry
    // mesh is a set of vertices and lines (connections between vertices)
    gsMesh<> mesh;
    pGeom->controlNet(mesh);
 
    // ======================================================================
    // writing to paraview
    // ======================================================================
 
    if (!output.empty())
    {
        std::string out = output + "Geometry";
        gsInfo << "Writing the geometry to a paraview file: " << out
                  << "\n\n";
 
        gsWriteParaview(*pGeom, out);
 
        out = output + "Basis";
        gsInfo << "Writing the basis to a paraview file: " << out
                  << "\n\n";
 
        gsWriteParaview(basis, out);
 
        out = output + "ContolNet";
        gsInfo << "Writing the control net to a paraview file: " << out
                  << "\n" << "\n";
 
        gsWriteParaview(mesh, out);
    }
    else
        gsInfo << "Done. No output created, re-run with --output <fn> to get a ParaView "
                  "file containing the solution.\n";
 
    return 0;
}