Implementation challenges in Spline based FEM (Isogeometric analysis) in Caribou

[bhagath555]

Dear Caribou Community,

At the University of Luxembourg, we have successfully implemented Linear and Non-Linear Isogeometric Analysis. The formulation is working well, and we have validated it with simple standard geometries. However, there are some challenges that we need to address to make the code generalize and optimize. I try will explain what are the challenges we face in detail. Please provide your thoughts and suggestions on possible solutions.

Source code is available here : Caribou_IGA

Note : I have used the words splines and NURBS inter changeably. NURBS are nothing but generalized splines.

First challenge

Our first challenge is the spline geometry reader. With the inspired from available VTKReader, we have developed a NURBSReader that can read spline data from a TXT file and creates the spline topology. This module works well in C++; however, we are currently unable to generate the Python bindings for it. I kindly request you to review the code and suggest any improvements.

Also, I haven't seen any available example in python where VTKReader is used. Most of the time, meshio module is used for FEM geometry processing. VTKReader is not even added in unit-tests. It would be nice if you comment on this too.

Second challenge:

There are fundamental differences between the implementation of Finite Element Analysis and Isogeometric Analysis. In the context of the Caribou code structure, IGA elements (spline/NURBS) do not have predefined nodes per element or Gauss points per element. The number of nodes per element depends on the degree of the element (as shown in the figure below). Once the degree is known, the rest of the computation is common for all cases.

I believe that with the help of the VTK module, Caribou FEM can determine the element type at compile time, allowing templates instantiation and the memory management to occur during compilation. We have also developed a module called CoreNurbs from scratch, which is analogous to the VTK module but not as powerful.

Currently, we can only obtain the following element information during runtime:

1. Degree of the element.
2. Number of nodes per element.
3. Number of Gauss nodes per element.

This information is known after reading the data from the TXT file.

In this type of scenario, I am considering two possible solutions, but I am not sure how to implement them.

1. Improve the TXT file reader module (CoreNurbs) to read the geometry information at compile time and assign it to spline elements.
2. Update the element information at runtime and perform memory assignment during runtime.

Please let me know if you require any additional information regarding splines or IGA. I would appreciate your thoughts on the possible solutions to over come this challenges.

@Ziemnono, @Sidaty1 @jnbrunet Please look into it.

Thank you,

Bhagath Mamindlapelly

[jnbrunet]

Hi Bhagath,

Concerning your first challenge, the VTKLoader python bindings can be found here. However it does seem like no unittest were done, or I can't find them anymore. In any cases, I'm pretty sure we tested it while developing it, hence it should be very close to a working solution.

Now, for your second challenge, elements in the Geometry module were made so that it is your element class that manages its nodes and Gauss points. You are therefore free to implement this as you wish. You could take advantage of the Eigen::Dynamic constant. It could be something that looks like (not tested at all here, just writing this very quickly):

namespace caribou::geometry {
template<unsigned int _Dimension>
struct SplineCurve;
template<unsigned int _Dimension>
struct traits<SplineCurve> {
    static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 1;
    static constexpr UNSIGNED_INTEGER_TYPE Dimension = _Dimension;
    static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = Eigen::Dynamic;
    static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = Eigen::Dynamic;
};
template<unsigned int _Dimension>
struct SplineCurve: public Element<SplineCurve<_Dimension>> {
    // Types
    using Base = Element<SplineCurve<_Dimension>>;
    using LocalCoordinates = typename Base::LocalCoordinates;
    using WorldCoordinates = typename Base::WorldCoordinates;
    using GaussNode = typename Base::GaussNode;

    /** Default empty constructor is restricted */
    SplineCurve() = deleted;

    /** Constructor from an Eigen matrix containing the positions of the Spline's nodes and Gauss points */
    template<typename EigenType, REQUIRES(EigenType::ColumnsAtCompileTime == _Dimension)>
    SplineCurve(const Eigen::EigenBase<EigenType> & nodes, const std::vector<GaussNode> & gauss_nodes) : p_gauss_nodes(gauss_nodes) {
        p_nodes.resize(nodes.rows(), _Dimension);
        for (unsigned int node_id; node_id < nodes.rows(); ++node_id) {
            for (unsigned int d; d< _Dimension; ++d) {
                p_nodes(node_id, d) = nodes(node_id, d);
            }
        }
    }
    
private:
    // Implementations
    friend struct Element<Derived>;
    inline auto get_number_of_nodes() const {return p_nodes.rows();}
    inline auto get_number_of_gauss_nodes() const {return p_gauss_nodes.size();}
    inline auto get_node(const UNSIGNED_INTEGER_TYPE & index) const {return WorldCoordinates(p_nodes.row(index));};
    inline auto get_nodes() const -> const auto & {return p_nodes;};
    inline auto get_center() const {return Base::world_coordinates(LocalCoordinates(0));};
    inline auto get_contains_local(const LocalCoordinates & xi, const FLOATING_POINT_TYPE & eps) const -> bool {
        /* ... */
    }
    inline auto get_L(const LocalCoordinates & xi) const -> Vector<Eigen::Dynamic> {
        /* ... */
    };

    inline auto get_dL(const LocalCoordinates & /*xi*/) const -> Matrix<Eigen::Dynamic, CanonicalDimension> {
        /* ... */
    };

protected:
    Matrix<Eigen::Dynamic, Dimension> p_nodes;
    std::vector<GaussNode> p_gauss_nodes;
}
}

Finally, I think it is great that you forked Caribou into a public repository. However, I think that your work would gain more visibility if you were creating a pull request to bring it into this repository. I am pretty sure that many could benefit from what you have done, and it would also be easier for us to help you. Just something to consider :)

Let me know if something wasn't clear or if I missed something

Jean-Nicolas

[bhagath555]

Thank you @jnbrunet for the response.

• I have tried to imitate VTKLoader for NURBSReader's python bindings.
• Regarding number of control points, The solution seems quite interesting.

I'll try to implement what you suggested and update the process.

Thanks again ;)
Bhagath

[bhagath555]

Hello @jnbrunet , I tried implementing the ideas you suggested it worked well. However, there are few challenges that i am struggling to resolve. Could you please go them and suggest how can I tackle.

1. Implementation of managing the varying nodes per element is implemented successfully, but there is a issue in python binding of a data type in CaribouSplineTopology.
   To handle different degree spline geometries, I have changed d_indices data type from Data<sofa::type::vector<sofa::type::fixed_array<PointID, 9>>> to Data<sofa::type::vector<sofa::type::vector<PointID>>>. This change worked well in C++ test, but not in python.
   I suspect there is a binding issue in using sofa::type::vector<sofa::type::vector<>>, because the same code worked when I replaced it with fixed_array<> for a particular problem.

2. As you mentioned in earlier, I have followed VTKReader code and came up with this code, but it is not working. I am accessing NURBS Reader function as below

#include <pybind11/pybind11.h>
#include <Caribou/Topology/config.h>
#include <Caribou/Topology/IO/NURBSReader.h>
#ifdef CARIBOU_WITH_VTK
#include <Caribou/Topology/IO/VTKReader.h>
#endif
namespace caribou::topology::io::bindings {

#ifdef CARIBOU_WITH_VTK
template<UNSIGNED_INTEGER_TYPE Dimension>
void add_reader(pybind11::module & m) {
	std::string name = "VTKReader" + std::to_string(Dimension) + "D";
	pybind11::class_<VTKReader<Dimension>> c(m, name.c_str());

	c.def_static("Read", &VTKReader<Dimension>::Read);
	c.def("__str__", [](const VTKReader<Dimension> & self) {
    	std::stringstream ss;
    	ss << self;
    	return ss.str();
	});
}
#endif

template<UNSIGNED_INTEGER_TYPE Dimension>
void add_nurbs_reader(pybind11::module & m) {
	std::string name = "NURBSReader" + std::to_string(Dimension) + "D";
	pybind11::class_<NURBSReader<Dimension>> c(m, name.c_str());

	c.def_static("Read", &NURBSReader<Dimension>::Read);
	c.def("__str__", [](const NURBSReader<Dimension> & self) {
    	std::stringstream ss;
    	ss << self;
    	return ss.str();
	});
}

void create_IO(pybind11::module & m) {
	pybind11::module io = m.def_submodule("IO");

#ifdef CARIBOU_WITH_VTK
	add_reader<1>(m);
	add_reader<2>(m);
	add_reader<3>(m);

	io.def("VTKReader", [](const std::string & filepath, unsigned int dimension) {
    	if (dimension == 1) {
        	return pybind11::cast(VTKReader<1>::Read(filepath));
    	} else if (dimension == 2) {
        	return pybind11::cast(VTKReader<2>::Read(filepath));
    	} else if (dimension == 3) {
        	return pybind11::cast(VTKReader<3>::Read(filepath));
    	} else {
        	throw std::runtime_error("Trying to create a VTKReader with a dimension that is not 1, 2 or 3.");
    	}
	}, pybind11::arg("filepath"), pybind11::arg("dimension") = 3);
#endif
}

} // namespace caribou::topology::io::bindings

Thank you