Skip to content

Commit 0feffec

Browse files
th-skamth-skam
committed
Add bindings for computeStress
1 parent 43e7896 commit 0feffec

4 files changed

Lines changed: 143 additions & 0 deletions

File tree

‎bindings/Modules/src/SofaPython3/SofaDeformable/Binding_CauchyStress.cpp‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
/******************************************************************************
2+
* SofaPython3 plugin *
3+
* (c) 2021 CNRS, University of Lille, INRIA *
4+
* *
5+
* This program is free software; you can redistribute it and/or modify it *
6+
* under the terms of the GNU Lesser General Public License as published by *
7+
* the Free Software Foundation; either version 2.1 of the License, or (at *
8+
* your option) any later version. *
9+
* *
10+
* This program is distributed in the hope that it will be useful, but WITHOUT *
11+
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
12+
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
13+
* for more details. *
14+
* *
15+
* You should have received a copy of the GNU Lesser General Public License *
16+
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
17+
*******************************************************************************
18+
* Contact information: contact@sofa-framework.org *
19+
******************************************************************************/
20+
#include <SofaPython3/SofaDeformable/Binding_CauchyStress.h>
21+
22+
#include <pybind11/numpy.h>
23+
24+
#include <sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h>
25+
#include <sofa/core/objectmodel/Base.h>
26+
#include <sofa/defaulttype/VecTypes.h>
27+
#include <sofa/type/VoigtNotation.h>
28+
29+
#include <string>
30+
31+
// Bindings for CauchyStressEvaluator: the Cauchy stress a force field or a material reports for given
32+
// deformation gradients.
33+
34+
namespace sofapython3
35+
{
36+
namespace py = pybind11;
37+
38+
using RealArray = py::array_t<SReal, py::array::c_style | py::array::forcecast>;
39+
40+
template <class DataTypes>
41+
static py::array_t<SReal> computeCauchyStressFor(sofa::core::objectmodel::Base* component,
42+
const RealArray& deformationGradientsArray)
43+
{
44+
using Evaluator = sofa::component::solidmechanics::fem::elastic::CauchyStressEvaluator<DataTypes>;
45+
constexpr py::ssize_t dim = DataTypes::spatial_dimensions;
46+
47+
auto* evaluator = dynamic_cast<Evaluator*>(component);
48+
if (!evaluator)
49+
throw py::type_error("computeCauchyStress: '" + component->getName()
50+
+ "' is not a CauchyStressEvaluator<" + DataTypes::Name() + ">");
51+
52+
const auto deformationGradients = deformationGradientsArray.unchecked<4>(); // (nbElements, Q, dim, dim)
53+
const py::ssize_t nbElements = deformationGradients.shape(0);
54+
const py::ssize_t Q = deformationGradients.shape(1);
55+
56+
py::array_t<SReal> stress({nbElements, Q, dim, dim});
57+
auto stressView = stress.mutable_unchecked<4>();
58+
59+
typename Evaluator::DeformationGradient F;
60+
for (py::ssize_t e = 0; e < nbElements; ++e)
61+
{
62+
for (py::ssize_t q = 0; q < Q; ++q)
63+
{
64+
for (auto i = 0; i < dim; ++i)
65+
for (auto j = 0; j < dim; ++j)
66+
F(i, j) = deformationGradients(e, q, i, j);
67+
68+
const auto stressVoigt = evaluator->computeStress(F, static_cast<sofa::Size>(e));
69+
for (sofa::Size v = 0; v < sofa::type::NumberOfIndependentElements<dim>; ++v)
70+
{
71+
const auto [i, j] = sofa::type::toTensorIndices<dim>(v);
72+
stressView(e, q, i, j) = stressVoigt[v];
73+
stressView(e, q, j, i) = stressVoigt[v];
74+
}
75+
}
76+
}
77+
return stress;
78+
}
79+
80+
static py::array_t<SReal> computeCauchyStress(sofa::core::objectmodel::Base* component,
81+
RealArray deformationGradients)
82+
{
83+
using namespace sofa::defaulttype;
84+
85+
if (!component)
86+
throw py::value_error("computeCauchyStress: component is None");
87+
88+
if (deformationGradients.ndim() != 4 || deformationGradients.shape(2) != deformationGradients.shape(3))
89+
throw py::value_error("computeCauchyStress: deformationGradients must be a 4D array "
90+
"(nb_elements, nb_quadrature_points, spatial_dimension, spatial_dimension)");
91+
92+
switch (deformationGradients.shape(3))
93+
{
94+
case 1: return computeCauchyStressFor<Vec1Types>(component, deformationGradients);
95+
case 2: return computeCauchyStressFor<Vec2Types>(component, deformationGradients);
96+
case 3: return computeCauchyStressFor<Vec3Types>(component, deformationGradients);
97+
default:
98+
throw py::value_error("computeCauchyStress: unsupported dimension "
99+
+ std::to_string(deformationGradients.shape(3)));
100+
}
101+
}
102+
103+
void moduleAddCauchyStress(py::module& m)
104+
{
105+
m.def("computeCauchyStress", &computeCauchyStress,
106+
py::arg("component"), py::arg("deformationGradients"),
107+
"Cauchy stress reported by a CauchyStressEvaluator (a force field or a material) for deformation "
108+
"gradients of shape (nb_elements, nb_quadrature_points, spatial_dimension, spatial_dimension), "
109+
"element e being the component's element e: returns the full symmetric stress tensors, same shape.");
110+
}
111+
112+
} // namespace sofapython3

‎bindings/Modules/src/SofaPython3/SofaDeformable/Binding_CauchyStress.h‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
/******************************************************************************
2+
* SofaPython3 plugin *
3+
* (c) 2021 CNRS, University of Lille, INRIA *
4+
* *
5+
* This program is free software; you can redistribute it and/or modify it *
6+
* under the terms of the GNU Lesser General Public License as published by *
7+
* the Free Software Foundation; either version 2.1 of the License, or (at *
8+
* your option) any later version. *
9+
* *
10+
* This program is distributed in the hope that it will be useful, but WITHOUT *
11+
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
12+
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
13+
* for more details. *
14+
* *
15+
* You should have received a copy of the GNU Lesser General Public License *
16+
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
17+
*******************************************************************************
18+
* Contact information: contact@sofa-framework.org *
19+
******************************************************************************/
20+
#pragma once
21+
22+
#include <pybind11/pybind11.h>
23+
24+
namespace sofapython3
25+
{
26+
void moduleAddCauchyStress(pybind11::module& m);
27+
}

‎bindings/Modules/src/SofaPython3/SofaDeformable/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,15 @@
11
project(Bindings.Modules.SofaDeformable)
22

33
set(SOURCE_FILES
4+
${CMAKE_CURRENT_SOURCE_DIR}/Binding_CauchyStress.cpp
45
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LameParameters.cpp
56
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSpring.cpp
67
${CMAKE_CURRENT_SOURCE_DIR}/Binding_SpringForceField.cpp
78
${CMAKE_CURRENT_SOURCE_DIR}/Module_SofaDeformable.cpp
89
)
910

1011
set(HEADER_FILES
12+
${CMAKE_CURRENT_SOURCE_DIR}/Binding_CauchyStress.h
1113
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LameParameters.h
1214
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSpring.h
1315
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSpring_doc.h

‎bindings/Modules/src/SofaPython3/SofaDeformable/Module_SofaDeformable.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
#include <SofaPython3/SofaDeformable/Binding_LinearSpring.h>
2323
#include <SofaPython3/SofaDeformable/Binding_SpringForceField.h>
2424
#include <SofaPython3/SofaDeformable/Binding_LameParameters.h>
25+
#include <SofaPython3/SofaDeformable/Binding_CauchyStress.h>
2526

2627

2728
namespace py { using namespace pybind11; }
@@ -36,6 +37,7 @@ PYBIND11_MODULE(SofaDeformable, m)
3637
moduleAddLinearSpring(m);
3738
moduleAddSpringForceField(m);
3839
moduleAddLameParameters(m);
40+
moduleAddCauchyStress(m);
3941
}
4042

4143
} // namespace sofapython3

0 commit comments

Comments
 (0)