diff --git a/bindings/Modules/CMakeLists.txt b/bindings/Modules/CMakeLists.txt index bdee373d..e47f5307 100644 --- a/bindings/Modules/CMakeLists.txt +++ b/bindings/Modules/CMakeLists.txt @@ -6,6 +6,7 @@ set(MODULEBINDINGS_MODULE_LIST SofaLinearSolver SofaLinearSystem SofaConstraintSolver + SofaFEM ) sofa_find_package(Sofa.GL QUIET) diff --git a/bindings/Modules/src/SofaPython3/SofaFEM/Binding_Quadrature.cpp b/bindings/Modules/src/SofaPython3/SofaFEM/Binding_Quadrature.cpp new file mode 100644 index 00000000..f7135584 --- /dev/null +++ b/bindings/Modules/src/SofaPython3/SofaFEM/Binding_Quadrature.cpp @@ -0,0 +1,291 @@ +/****************************************************************************** +* SofaPython3 plugin * +* (c) 2021 CNRS, University of Lille, INRIA * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +// Bindings for sofa::fem::FiniteElement's reference-element functions. +// Exposes the element kernel (quadrature rule, shape functions, reference gradients) so that an +// integration loop can be carried out using SofaPython3's API. + +namespace sofapython3 +{ +namespace py = pybind11; + +using RealArray = py::array_t; +// Topology connectivity arrives as whatever integer type the SOFA container uses; forcecast makes +// the binding indifferent to it. +using IndexArray = py::array_t; + +// FiniteElement class dispatch: (element name, ambient dimension) +template +struct ElementTag { using FiniteElement = sofa::fem::FiniteElement; }; + +template +static py::tuple withElement(const std::string& element, py::ssize_t dim, Function&& f) +{ + using namespace sofa::defaulttype; + namespace geo = sofa::geometry; + + if (element == "Edge" && dim == 1) return f(ElementTag{}); + if (element == "Edge" && dim == 2) return f(ElementTag{}); + if (element == "Edge" && dim == 3) return f(ElementTag{}); + if (element == "Triangle" && dim == 2) return f(ElementTag{}); + if (element == "Triangle" && dim == 3) return f(ElementTag{}); + if (element == "Quad" && dim == 2) return f(ElementTag{}); + if (element == "Quad" && dim == 3) return f(ElementTag{}); + if (element == "Tetrahedron" && dim == 3) return f(ElementTag{}); + if (element == "Hexahedron" && dim == 3) return f(ElementTag{}); + + throw py::value_error("SofaFEM: unsupported element '" + element + "' for dimension " + std::to_string(dim)); +} + +// Reference-space quadrature points' data for one FiniteElement specialization: +// - weights +// - shape functions' values N_a +// - shape functions' gradients dN_a/dxi +template +static py::tuple quadratureDataFor(sofa::Size degree) +{ + constexpr py::ssize_t nbNodes = FE::NumberOfNodesInElement; + constexpr py::ssize_t topoDim = FE::TopologicalDimension; + + std::span rule; + try + { + rule = FE::quadratureRule(degree); + } + catch (const std::exception& error) + { + throw py::value_error(error.what()); + } + const py::ssize_t Q = static_cast(rule.size()); + + using Real = typename FE::Real; + + py::array_t weights({Q}); + py::array_t shapeFunctions({Q, nbNodes}); + py::array_t shapeFunctionGrads({Q, nbNodes, topoDim}); + auto weightsView = weights.template mutable_unchecked<1>(); + auto functionsView = shapeFunctions.template mutable_unchecked<2>(); + auto gradientsView = shapeFunctionGrads.template mutable_unchecked<3>(); + + py::ssize_t q = 0; + for (const auto& [referencePoint, weight] : rule) + { + weightsView(q) = weight; + + const auto N = FE::shapeFunctions(referencePoint); + for (auto a = 0; a < nbNodes; ++a) + functionsView(q, a) = N[a]; + + const auto gradient = FE::gradientShapeFunctions(referencePoint); // Mat + for (auto a = 0; a < nbNodes; ++a) + for (auto j = 0; j < topoDim; ++j) + gradientsView(q, a, j) = gradient[a][j]; + + ++q; + } + return py::make_tuple(weights, shapeFunctions, shapeFunctionGrads); +} + +static py::tuple quadratureData(const std::string& element, py::ssize_t dim, sofa::Size degree) +{ + return withElement(element, dim, [&](auto tag) + { + return quadratureDataFor(degree); + }); +} + +// Physical-space shape function gradients dN_a/dx and the integration measure per quadrature point +// - physical space gradients dN_a/dxi +// - measures +template +static py::tuple elementMappingFor(const RealArray& nodeCoordinatesArray, const RealArray& referenceGradientsArray) +{ + using Real = typename FE::Real; + using Coord = typename FE::Coord; + using Helper = typename FE::Helper; + constexpr py::ssize_t nbNodes = FE::NumberOfNodesInElement; + constexpr py::ssize_t spatialDim = FE::spatial_dimensions; + constexpr py::ssize_t topoDim = FE::TopologicalDimension; + + const auto nodeCoordinates = nodeCoordinatesArray.unchecked<2>(); // (nbNodes, spatialDim) + const auto referenceGrads = referenceGradientsArray.unchecked<3>(); // (Q, nbNodes, topoDim) + const py::ssize_t Q = referenceGrads.shape(0); + + std::array elementNodes; + for (auto a = 0; a < nbNodes; ++a) + for (auto d = 0; d < spatialDim; ++d) + elementNodes[a][d] = nodeCoordinates(a, d); + + py::array_t physicalGradients({Q, nbNodes, spatialDim}); + py::array_t measures({Q}); + auto gradientsView = physicalGradients.template mutable_unchecked<3>(); + auto measuresView = measures.template mutable_unchecked<1>(); + + for (py::ssize_t q = 0; q < Q; ++q) + { + sofa::type::Mat referenceGradient; + for (auto a = 0; a < nbNodes; ++a) + for (auto j = 0; j < topoDim; ++j) + referenceGradient[a][j] = referenceGrads(q, a, j); + + const auto jacobian = Helper::jacobianFromReferenceToPhysical(elementNodes, referenceGradient); + measuresView(q) = sofa::type::absGeneralizedDeterminant(jacobian); // |det J|, or sqrt(det(J^T J)) if embedded + const auto inverseJacobian = sofa::type::inverse(jacobian); // inverse, or left pseudo-inverse if embedded + + for (auto a = 0; a < nbNodes; ++a) + { + const auto physicalGradient = inverseJacobian.transposed() * referenceGradient[a]; // dN_a/dx + for (auto d = 0; d < spatialDim; ++d) + gradientsView(q, a, d) = physicalGradient[d]; + } + } + return py::make_tuple(physicalGradients, measures); +} + +static py::tuple elementMapping(const std::string& element, RealArray nodeCoordinates, RealArray referenceGradients) +{ + if (nodeCoordinates.ndim() != 2) + throw py::value_error("element_mapping: node_coordinates must be a 2D array (nodes_per_element, spatial_dimension)"); + + return withElement(element, nodeCoordinates.shape(1), [&](auto tag) + { + return elementMappingFor(nodeCoordinates, referenceGradients); + }); +} + +// Same mapping, but for every element of a mesh in one call. A Python integration loop that calls +// element_mapping per element pays a pybind crossing and two array allocations per element, which +// dominates the cost of walking a large mesh; here the whole sweep is one crossing and one pair of +// allocations. The per-quadrature-point block is identical to elementMappingFor. +template +static py::tuple elementMappingBatchFor(const RealArray& nodeCoordinatesArray, + const IndexArray& nodeIndicesArray, + const RealArray& referenceGradientsArray) +{ + using Real = typename FE::Real; + using Coord = typename FE::Coord; + using Helper = typename FE::Helper; + constexpr py::ssize_t nbNodes = FE::NumberOfNodesInElement; + constexpr py::ssize_t spatialDim = FE::spatial_dimensions; + constexpr py::ssize_t topoDim = FE::TopologicalDimension; + + const auto nodeCoordinates = nodeCoordinatesArray.unchecked<2>(); // (nbMeshNodes, spatialDim) + const auto nodeIndices = nodeIndicesArray.unchecked<2>(); // (nbElements, nbNodes) + const auto referenceGrads = referenceGradientsArray.unchecked<3>(); // (Q, nbNodes, topoDim) + + const py::ssize_t nbElements = nodeIndices.shape(0); + const py::ssize_t Q = referenceGrads.shape(0); + const py::ssize_t nbMeshNodes = nodeCoordinates.shape(0); + + if (nodeIndices.shape(1) != nbNodes) + throw py::value_error("element_mapping_batch: node_indices has the wrong number of nodes per element"); + + // The reference gradients do not depend on the element, so lift them out of the element loop. + std::vector> referenceGradient(Q); + for (py::ssize_t q = 0; q < Q; ++q) + for (auto a = 0; a < nbNodes; ++a) + for (auto j = 0; j < topoDim; ++j) + referenceGradient[q][a][j] = referenceGrads(q, a, j); + + py::array_t physicalGradients({nbElements, Q, nbNodes, spatialDim}); + py::array_t measures({nbElements, Q}); + auto gradientsView = physicalGradients.template mutable_unchecked<4>(); + auto measuresView = measures.template mutable_unchecked<2>(); + + std::array elementNodes; + for (py::ssize_t e = 0; e < nbElements; ++e) + { + for (auto a = 0; a < nbNodes; ++a) + { + const auto node = static_cast(nodeIndices(e, a)); + if (node < 0 || node >= nbMeshNodes) + throw py::value_error("element_mapping_batch: node index out of range"); + for (auto d = 0; d < spatialDim; ++d) + elementNodes[a][d] = nodeCoordinates(node, d); + } + + for (py::ssize_t q = 0; q < Q; ++q) + { + const auto jacobian = Helper::jacobianFromReferenceToPhysical(elementNodes, referenceGradient[q]); + measuresView(e, q) = sofa::type::absGeneralizedDeterminant(jacobian); + const auto inverseJacobian = sofa::type::inverse(jacobian); + + for (auto a = 0; a < nbNodes; ++a) + { + const auto physicalGradient = inverseJacobian.transposed() * referenceGradient[q][a]; + for (auto d = 0; d < spatialDim; ++d) + gradientsView(e, q, a, d) = physicalGradient[d]; + } + } + } + return py::make_tuple(physicalGradients, measures); +} + +static py::tuple elementMappingBatch(const std::string& element, RealArray nodeCoordinates, + IndexArray nodeIndices, RealArray referenceGradients) +{ + if (nodeCoordinates.ndim() != 2) + throw py::value_error("element_mapping_batch: nodes must be a 2D array (nb_mesh_nodes, spatial_dimension)"); + if (nodeIndices.ndim() != 2) + throw py::value_error("element_mapping_batch: node_indices must be a 2D array (nb_elements, nodes_per_element)"); + + return withElement(element, nodeCoordinates.shape(1), [&](auto tag) + { + return elementMappingBatchFor( + nodeCoordinates, nodeIndices, referenceGradients); + }); +} + +void moduleAddQuadrature(py::module& m) +{ + m.def("quadrature_data", &quadratureData, + py::arg("element"), py::arg("dim"), py::arg("degree"), + "Reference-space quadrature data for the element at the given degree: " + "returns (quadrature weights, shape functions, shape function gradients)."); + + m.def("element_mapping", &elementMapping, + py::arg("element"), py::arg("node_coordinates"), py::arg("reference_gradients"), + "Reference->physical mapping for one element (reuses SOFA's jacobianFromReferenceToPhysical, " + "inverse and absGeneralizedDeterminant): returns (physical shape-function gradients dN_a/dx, " + "integration measures) per quadrature point; handles square and embedded (rectangular Jacobian) elements."); + + m.def("element_mapping_batch", &elementMappingBatch, + py::arg("element"), py::arg("nodes"), py::arg("node_indices"), py::arg("reference_gradients"), + "Reference->physical mapping for every element of a mesh in one call, so a Python " + "integration loop pays one pybind crossing instead of one per element: returns " + "(physical shape-function gradients dN_a/dx of shape (nb_elements, nb_quadrature_points, " + "nodes_per_element, spatial_dimension), integration measures of shape (nb_elements, " + "nb_quadrature_points))."); +} + +} // namespace sofapython3 diff --git a/bindings/Modules/src/SofaPython3/SofaFEM/Binding_Quadrature.h b/bindings/Modules/src/SofaPython3/SofaFEM/Binding_Quadrature.h new file mode 100644 index 00000000..5359aa2c --- /dev/null +++ b/bindings/Modules/src/SofaPython3/SofaFEM/Binding_Quadrature.h @@ -0,0 +1,27 @@ +/****************************************************************************** +* SofaPython3 plugin * +* (c) 2021 CNRS, University of Lille, INRIA * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +namespace sofapython3 +{ +void moduleAddQuadrature(pybind11::module& m); +} diff --git a/bindings/Modules/src/SofaPython3/SofaFEM/CMakeLists.txt b/bindings/Modules/src/SofaPython3/SofaFEM/CMakeLists.txt new file mode 100644 index 00000000..f6111c32 --- /dev/null +++ b/bindings/Modules/src/SofaPython3/SofaFEM/CMakeLists.txt @@ -0,0 +1,26 @@ +project(Bindings.Modules.SofaFEM) + +set(SOURCE_FILES + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Quadrature.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Module_SofaFEM.cpp +) + +set(HEADER_FILES + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Quadrature.h +) + +if (NOT TARGET SofaPython3::Plugin) + find_package(SofaPython3 REQUIRED COMPONENTS Plugin Bindings.Sofa) +endif() + +sofa_find_package(Sofa.FEM REQUIRED) + +SP3_add_python_module( + TARGET ${PROJECT_NAME} + PACKAGE Bindings.Modules + MODULE SofaFEM + DESTINATION Sofa + SOURCES ${SOURCE_FILES} + HEADERS ${HEADER_FILES} + DEPENDS Sofa.FEM SofaPython3::Plugin SofaPython3::Bindings.Sofa.Core +) diff --git a/bindings/Modules/src/SofaPython3/SofaFEM/Module_SofaFEM.cpp b/bindings/Modules/src/SofaPython3/SofaFEM/Module_SofaFEM.cpp new file mode 100644 index 00000000..ff000ca2 --- /dev/null +++ b/bindings/Modules/src/SofaPython3/SofaFEM/Module_SofaFEM.cpp @@ -0,0 +1,36 @@ +/****************************************************************************** +* SofaPython3 plugin * +* (c) 2021 CNRS, University of Lille, INRIA * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include + +#include + +namespace py { using namespace pybind11; } + +namespace sofapython3 +{ + +PYBIND11_MODULE(SofaFEM, m) +{ + m.doc() = "Bindings for the sofa::fem finite-element machinery"; + + moduleAddQuadrature(m); +} + +} // namespace sofapython3 diff --git a/examples/source-terms/example-quadrature-rules.py b/examples/source-terms/example-quadrature-rules.py new file mode 100644 index 00000000..0a837b91 --- /dev/null +++ b/examples/source-terms/example-quadrature-rules.py @@ -0,0 +1,174 @@ +"""Solves a linear elastic problem with a prescribed load on the RHS. The load +is assembled from a manufactured body-force density. The assembly is done by +integrating the density using quadrature information provided by the SofaFEM +quadrature bindings. + +The resulting displacement should match the manufactured solution: +ue(x,y) = 0.03 * (sin(2 pi x) sin(2 pi y), sin(2 pi x) sin(2 pi y)) on [0,1]^2, +with zero Dirichlet BCs on all boundaries. + +The mesh is coarse enough to reveal a difference in accuracy between the +quadrature degrees in use. + +White: reference solution ue +Red & Green: solved using a load assembled with quadrature order 1 & 2, respectively. +""" +import numpy as np +import Sofa +import Sofa.Core +from Sofa import SofaFEM + +AMPLITUDE = 0.03 +WAVENUMBER = 2 * np.pi + + +def manufactured_body_force_density(positions, young_modulus, poisson_ratio): + """The RHS side f of the linear elasticity PDE: + f = Ak²[(λ+3μ) sin(kx)sin(ky) − (λ+μ) cos(kx)cos(ky)], with A=AMPLITUDE, k=WAVENUMBER + Acts as a source term, exciting the displacement field to yield ue, the manufactured solution""" + mu = young_modulus / (2 * (1 + poisson_ratio)) + lam = young_modulus * poisson_ratio / (1 - poisson_ratio ** 2) + + x, y = positions[:, 0], positions[:, 1] + k = WAVENUMBER + s = AMPLITUDE * np.sin(k * x) * np.sin(k * y) + c = AMPLITUDE * np.cos(k * x) * np.cos(k * y) + f = k ** 2 * ((lam + 3 * mu) * s - (lam + mu) * c) + return np.column_stack([f, f]) + + +def assemble_nodal_force(node_coordinates, node_indices, quadrature_degree, nodal_density): + """Assembles a per-node force vector by integrating a nodal source density: + F_a = sum_q w_q * measure_q * N_a(q) * r(q), with r(q) = sum_a N_a(q) * density_a""" + + # First get the quadrature data and mapping information for the corresponding topological element + weights, shape_functions, reference_gradients = SofaFEM.quadrature_data("Triangle", 2, quadrature_degree) + _, measures = SofaFEM.element_mapping_batch("Triangle", node_coordinates, node_indices, reference_gradients) + + element_density = nodal_density[node_indices] # density at each element's local nodes, (E, nbNodes, dof) + nb_nodes_per_element = node_indices.shape[1] + + nodal_force = np.zeros_like(nodal_density) + for weight, measure, shape_function in zip(weights, measures.T, shape_functions): + # r(q), per element: interpolate the local nodal densities with the shape functions + density_at_q = sum(shape_function[a] * element_density[:, a] for a in range(nb_nodes_per_element)) + + # scatter w_q * measure_q * N_a(q) * r(q) onto each local node a + for a in range(nb_nodes_per_element): + contribution = weight * measure * shape_function[a] # (E,) + np.add.at(nodal_force, node_indices[:, a], contribution[:, None] * density_at_q) + return nodal_force + + +class AssembleSourceTerms(Sofa.Core.Controller): + """Assembles the manufactured body-force load for and assigns it to the Data of a ConstantForceField.""" + def __init__(self, node, quadrature_degree, young_modulus, poisson_ratio, *args, **kwargs): + Sofa.Core.Controller.__init__(self, *args, **kwargs) + self.node = node + self.quadrature_degree = quadrature_degree + self.young_modulus = young_modulus + self.poisson_ratio = poisson_ratio + + def onSimulationInitDoneEvent(self, _): + node_coordinates = self.node.dofs.rest_position.array() + node_indices = np.array(self.node.Triangle_topo.triangles.array()) + + nodal_density = manufactured_body_force_density(node_coordinates, self.young_modulus, self.poisson_ratio) + nodal_force = assemble_nodal_force(node_coordinates, node_indices, self.quadrature_degree, nodal_density) + + self.node.addObject("ConstantForceField", name="source", + indices=list(range(len(nodal_force))), forces=nodal_force.tolist()).init() + + +def createScene(root): + root.dt = 1 + root.gravity = [0, 0, 0] + + root.addObject("RequiredPlugin", pluginName=[ + "Sofa.Component.AnimationLoop", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.Engine.Select", + "Sofa.Component.LinearSolver.Direct", + "Sofa.Component.MechanicalLoad", + "Sofa.Component.ODESolver.Backward", + "Sofa.Component.SolidMechanics.FEM.Elastic", + "Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Dynamic", + "Sofa.Component.Topology.Container.Grid", + "Sofa.Component.Topology.Mapping", + "Sofa.Component.Visual", + "Sofa.GL.Component.Rendering2D", + ]) + + root.addObject("DefaultAnimationLoop") + root.addObject("VisualStyle", displayFlags="showBehaviorModels showVisualModels showWireframe") + + with root.addChild("Labels") as labels: + labels.addObject("OglLabel", name="Reference", x="10", y="810", fontsize="16", color="1 1 1 1") + labels.addObject("OglLabel", name="Deg1", x="10", y="832", fontsize="16", color="1 0 0 1") + labels.addObject("OglLabel", name="Deg2", x="10", y="854", fontsize="16", color="0 1 0 1") + + quadratureDegrees = root.addChild("QuadratureDegrees") + quadratureDegrees.addObject("RegularGridTopology", name="grid", n="5 5 1", min="0 0 0", max="1 1 0") + + # Reference Mesh: nodes are placed exactly on the manufactured solution, i.e. at grid position + ue(grid position). + with quadratureDegrees.addChild("ManufacturedSolution") as manufacturedSolution: + # Compute the exact displacement from the analytical expression + grid_coordinates = np.linspace(0, 1, 5) + grid_x, grid_y = np.meshgrid(grid_coordinates, grid_coordinates) + displacement = AMPLITUDE * np.sin(WAVENUMBER * grid_x) * np.sin(WAVENUMBER * grid_y) + referencePositions = np.column_stack([(grid_x + displacement).ravel(), (grid_y + displacement).ravel()]) + + # Assign the computed displacement the dofs + manufacturedSolution.addObject("MechanicalObject", name="dofs", template="Vec2", showObject="true", drawMode="1", showObjectScale="0.005", showColor="1 1 1 0.7", + position=referencePositions.tolist()) + manufacturedSolution.addObject("TriangleSetTopologyContainer", name="Triangle_topo") + manufacturedSolution.addObject("TriangleSetTopologyModifier", name="Modifier") + manufacturedSolution.addObject("TriangleSetGeometryAlgorithms", template="Vec2", name="GeomAlgo", drawEdges="1", drawColorEdges="1 1 1 1") + manufacturedSolution.addObject("Quad2TriangleTopologicalMapping", input="@../grid", output="@Triangle_topo", swapping="True") + + boundaryBox = ("-0.125 -0.125 -0.125 0.125 1.125 0.125 " + "0.875 -0.125 -0.125 1.125 1.125 0.125 " + "-0.125 -0.125 -0.125 1.125 0.125 0.125 " + "-0.125 0.875 -0.125 1.125 1.125 0.125") + + youngModulus = 100000 + poissonRatio = 0.3 + + def addQuadratureDegreeNode(name, quadrature_degree, showColor, edgeColor): + # Starts at rest, loaded by the manufactured body force, integrated with the given quadrature degree. + with quadratureDegrees.addChild(name) as node: + node.addObject("NewtonRaphsonSolver", name="newton", maxNbIterationsNewton="1", absoluteResidualStoppingThreshold="1e-8") + node.addObject("StaticSolver", newtonSolver="@newton") + node.addObject("SparseLDLSolver", name="linear_solver", template="CompressedRowSparseMatrixd") + + node.addObject("MechanicalObject", name="dofs", template="Vec2", src="@../grid", showObject="true", drawMode="1", showObjectScale="0.005", showColor=showColor) + node.addObject("BoxROI", name="boundary", template="Vec2", box=boundaryBox) + node.addObject("FixedProjectiveConstraint", name="dirichlet", indices="@boundary.indices") + + node.addObject("TriangleSetTopologyContainer", name="Triangle_topo") + node.addObject("TriangleSetTopologyModifier", name="Modifier") + node.addObject("TriangleSetGeometryAlgorithms", template="Vec2", name="GeomAlgo", drawEdges="1", drawColorEdges=edgeColor) + node.addObject("Quad2TriangleTopologicalMapping", input="@../grid", output="@Triangle_topo", swapping="True") + + node.addObject("LinearSmallStrainFEMForceField", name="FEM", template="Vec2,Triangle", + youngModulus=youngModulus, poissonRatio=poissonRatio, topology="@Triangle_topo") + + node.addObject(AssembleSourceTerms( + node=node, quadrature_degree=quadrature_degree, + young_modulus=youngModulus, poisson_ratio=poissonRatio, name="assembleSourceTerm")) + + addQuadratureDegreeNode("QuadratureDegree1", 1, showColor="1 0 0 0.9", edgeColor="1 0 0 1") + addQuadratureDegreeNode("QuadratureDegree2", 2, showColor="0 1 0 1", edgeColor="0 1 0 1") + + return root + + +def main(): + root = Sofa.Core.Node("root") + createScene(root) + Sofa.Simulation.init(root) + + +if __name__ == '__main__': + main()