From c4557cc9e0bb2c999614358569dbac84c8959fd8 Mon Sep 17 00:00:00 2001 From: Quoc-Minh Ton-That Date: Thu, 18 Jul 2024 15:14:55 -0400 Subject: [PATCH] Upgrade bindings to simpler API --- bindings/pypbat/CMakeLists.txt | 1 - bindings/pypbat/fem/CMakeLists.txt | 101 -- bindings/pypbat/fem/Gradient.cpp | 6 +- bindings/pypbat/fem/HyperElasticPotential.cpp | 32 +- bindings/pypbat/fem/Jacobian.cpp | 18 +- bindings/pypbat/fem/Laplacian.cpp | 10 +- bindings/pypbat/fem/MassMatrix.cpp | 25 +- bindings/pypbat/fem/ShapeFunctions.cpp | 14 +- bindings/pypbat/fem/bindings.py | 1396 ----------------- python/examples/diffusion_smoothing.py | 7 +- python/examples/elasticity.py | 25 +- python/examples/elasticity_higher_order.py | 23 +- python/examples/heat.py | 15 +- python/examples/ipc.py | 18 +- python/examples/laplace.py | 5 +- python/examples/modes.py | 11 +- python/pbatoolkit/__init__.py | 3 + python/pbatoolkit/fem.py | 301 ---- source/pbat/fem/Jacobian.h | 24 +- source/pbat/fem/ShapeFunctions.h | 11 +- 20 files changed, 168 insertions(+), 1878 deletions(-) delete mode 100644 bindings/pypbat/fem/bindings.py delete mode 100644 python/pbatoolkit/fem.py diff --git a/bindings/pypbat/CMakeLists.txt b/bindings/pypbat/CMakeLists.txt index 61e7639..137a9a1 100644 --- a/bindings/pypbat/CMakeLists.txt +++ b/bindings/pypbat/CMakeLists.txt @@ -10,7 +10,6 @@ add_subdirectory(profiling) target_link_libraries(PhysicsBasedAnimationToolkit_Python PRIVATE - # PBAT_PyFem pybind11::headers Python::Module ) \ No newline at end of file diff --git a/bindings/pypbat/fem/CMakeLists.txt b/bindings/pypbat/fem/CMakeLists.txt index f6644fd..ad9813d 100644 --- a/bindings/pypbat/fem/CMakeLists.txt +++ b/bindings/pypbat/fem/CMakeLists.txt @@ -1,104 +1,3 @@ -# add_library(PBAT_PyFem) -# set_target_properties(PBAT_PyFem -# PROPERTIES -# FOLDER "PhysicsBasedAnimationToolkit/bindings/fem" -# WINDOWS_EXPORT_ALL_SYMBOLS ON -# ) - -# list(APPEND _pbat_fem_types -# "Gradient" -# "HyperElasticPotential" -# "Jacobian" -# "LaplacianMatrix" -# "LoadVector" -# "MassMatrix" -# "Mesh" -# "ShapeFunctions" -# ) - -# set(_pbat_python_autogen_dir "gen") - -# foreach(_pbat_fem_type IN ITEMS ${_pbat_fem_types}) -# set(_pbat_python_fem_target "PBAT_PyFem_${_pbat_fem_type}") -# add_library(${_pbat_python_fem_target}) -# set_target_properties(${_pbat_python_fem_target} -# PROPERTIES -# FOLDER "PhysicsBasedAnimationToolkit/bindings/fem/${_pbat_fem_type}" -# WINDOWS_EXPORT_ALL_SYMBOLS ON -# ) -# target_link_libraries(${_pbat_python_fem_target} -# PUBLIC -# pybind11::headers -# Python::Module -# PRIVATE -# PhysicsBasedAnimationToolkit_PhysicsBasedAnimationToolkit -# ) - -# execute_process( -# COMMAND ${Python_EXECUTABLE} bindings.py --cmake=1 --type=${_pbat_fem_type} -# WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}" -# OUTPUT_VARIABLE _pbat_fem_bindings_impl -# ERROR_VARIABLE _pbat_fem_bindings_error -# RESULT_VARIABLE _pbat_bindings_exit_code -# ECHO_ERROR_VARIABLE -# ) - -# if(NOT _pbat_bindings_exit_code EQUAL 0) -# message(FATAL_ERROR "Failed to generate FEM binding implementations.") -# endif() - -# target_sources(${_pbat_python_fem_target} -# PUBLIC -# FILE_SET api -# TYPE HEADERS -# BASE_DIRS ${CMAKE_CURRENT_LIST_DIR} -# ) - -# target_sources(${_pbat_python_fem_target} -# PUBLIC -# FILE_SET api -# FILES -# "${_pbat_python_autogen_dir}/${_pbat_fem_type}.h" -# ) - -# target_sources(${_pbat_python_fem_target} -# PRIVATE -# "${_pbat_python_autogen_dir}/${_pbat_fem_type}.cpp" -# ${_pbat_fem_bindings_impl} -# ) - -# target_link_libraries(PBAT_PyFem -# PRIVATE -# ${_pbat_python_fem_target} -# ) -# endforeach() - -# target_link_libraries(PBAT_PyFem -# PUBLIC -# pybind11::headers -# Python::Module -# PRIVATE -# PhysicsBasedAnimationToolkit_PhysicsBasedAnimationToolkit -# ) -# target_sources(PBAT_PyFem -# PUBLIC -# FILE_SET api -# TYPE HEADERS -# BASE_DIRS ${CMAKE_CURRENT_LIST_DIR} -# ) - -# target_sources(PBAT_PyFem -# PUBLIC -# FILE_SET api -# FILES -# "Fem.h" -# ) - -# target_sources(PBAT_PyFem -# PRIVATE -# "Fem.cpp" -# ) - target_sources(PhysicsBasedAnimationToolkit_Python PUBLIC FILE_SET api diff --git a/bindings/pypbat/fem/Gradient.cpp b/bindings/pypbat/fem/Gradient.cpp index f9862a6..e92fc76 100644 --- a/bindings/pypbat/fem/Gradient.cpp +++ b/bindings/pypbat/fem/Gradient.cpp @@ -47,9 +47,11 @@ void BindGradient(pybind11::module& m) pyb::init const&, int>(), pyb::arg("mesh"), pyb::arg("GNe"), - pyb::arg("quadrature_order") = 1) + pyb::arg("quadrature_order") = 1, + "Construct Gradient operator from mesh mesh, using precomputed shape function " + "gradients GNe at quadrature points given by quadrature rule of order quadrature_order") .def_readonly("dims", &Gradient::mDims) - .def_readonly("order", &Gradient::mOrder) + .def_readonly("order", &Gradient::mOrder, "Polynomial order of the gradient") .def_readonly("quadrature_order", &Gradient::mQuadratureOrder) .def_property_readonly("shape", &Gradient::Shape) .def("to_matrix", &Gradient::ToMatrix); diff --git a/bindings/pypbat/fem/HyperElasticPotential.cpp b/bindings/pypbat/fem/HyperElasticPotential.cpp index 59985ff..398bba1 100644 --- a/bindings/pypbat/fem/HyperElasticPotential.cpp +++ b/bindings/pypbat/fem/HyperElasticPotential.cpp @@ -182,7 +182,11 @@ void BindHyperElasticPotential(pybind11::module& m) pyb::arg("Y") = 1e6, pyb::arg("nu") = 0.45, pyb::arg("energy") = EHyperElasticEnergy::StableNeoHookean, - pyb::arg("quadrature_order") = 1) + pyb::arg("quadrature_order") = 1, + "Construct a HyperElasticPotential on mesh mesh, given precomputed jacobian " + "determinants detJe and shape function gradients GNe at mesh element quadrature points " + "given by quadrature rule of order quadrature_order. The corresponding energy has " + "Young's modulus Y and Poisson's ratio nu.") .def( pyb::init< Mesh const&, @@ -198,9 +202,13 @@ void BindHyperElasticPotential(pybind11::module& m) pyb::arg("Y"), pyb::arg("nu"), pyb::arg("energy") = EHyperElasticEnergy::StableNeoHookean, - pyb::arg("quadrature_order") = 1) + pyb::arg("quadrature_order") = 1, + "Construct a HyperElasticPotential on mesh mesh, given precomputed jacobian " + "determinants detJe and shape function gradients GNe at mesh element quadrature points " + "given by quadrature rule of order quadrature_order. The corresponding energy has " + "piecewise constant (at elements) Young's modulus Y and Poisson's ratio nu.") .def_readonly("dims", &HyperElasticPotential::mDims) - .def_readonly("order", &HyperElasticPotential::mOrder) + //.def_readonly("order", &HyperElasticPotential::mOrder) .def_readonly("quadrature_order", &HyperElasticPotential::mQuadratureOrder) .def("precompute_hessian_sparsity", &HyperElasticPotential::PrecomputeHessianSparsity) .def( @@ -215,22 +223,30 @@ void BindHyperElasticPotential(pybind11::module& m) .def_property( "mue", [](HyperElasticPotential const& M) { return M.mue(); }, - [](HyperElasticPotential& M, Eigen::Ref const& mue) { M.mue() = mue; }) + [](HyperElasticPotential& M, Eigen::Ref const& mue) { M.mue() = mue; }, + "Piecewise constant (per-element) vector of first Lame coefficients") .def_property( "lambdae", [](HyperElasticPotential const& M) { return M.mue(); }, [](HyperElasticPotential& M, Eigen::Ref const& lambdae) { M.lambdae() = lambdae; - }) + }, + "Piecewise constant (per-element) vector of second Lame coefficients") .def_property_readonly( "UE", - [](HyperElasticPotential const& M) { return M.ElementPotentials(); }) + [](HyperElasticPotential const& M) { return M.ElementPotentials(); }, + "|#elements| vector of element hyper elastic potentials") .def_property_readonly( "GE", - [](HyperElasticPotential const& M) { return M.ElementGradients(); }) + [](HyperElasticPotential const& M) { return M.ElementGradients(); }, + "|#element nodes * dims|x|#elements| matrix of element hyper elastic potential " + "gradients") .def_property_readonly( "HE", - [](HyperElasticPotential const& M) { return M.ElementHessians(); }) + [](HyperElasticPotential const& M) { return M.ElementHessians(); }, + "|#element nodes * dims|x|#elements nodes * dims * #elements| matrix of element hyper " + "elastic " + "potential hessians") .def_property_readonly("shape", &HyperElasticPotential::Shape) .def("to_matrix", &HyperElasticPotential::ToMatrix); } diff --git a/bindings/pypbat/fem/Jacobian.cpp b/bindings/pypbat/fem/Jacobian.cpp index f2caff1..3bdc494 100644 --- a/bindings/pypbat/fem/Jacobian.cpp +++ b/bindings/pypbat/fem/Jacobian.cpp @@ -26,7 +26,9 @@ void BindJacobian(pybind11::module& m) return detJe; }, pyb::arg("mesh"), - pyb::arg("quadrature_order") = 1); + pyb::arg("quadrature_order") = 1, + "|#quad.pts.|x|#elements| matrix of element jacobian determinants at element quadrature " + "points"); m.def( "inner_product_weights", @@ -41,7 +43,9 @@ void BindJacobian(pybind11::module& m) return I; }, pyb::arg("mesh"), - pyb::arg("quadrature_order") = 1); + pyb::arg("quadrature_order") = 1, + "|#quad.pts.|x|#elements| matrix of quadrature weights multiplied by jacobian determinants " + "at element quadrature points "); m.def( "inner_product_weights", @@ -57,10 +61,12 @@ void BindJacobian(pybind11::module& m) }, pyb::arg("mesh"), pyb::arg("detJe"), - pyb::arg("quadrature_order") = 1); + pyb::arg("quadrature_order") = 1, + "|#quad.pts.|x|#elements| matrix of quadrature weights multiplied by jacobian determinants " + "at element quadrature points"); m.def( - "inner_product_weights", + "reference_positions", [](Mesh const& M, Eigen::Ref const& E, Eigen::Ref const& X, @@ -76,7 +82,9 @@ void BindJacobian(pybind11::module& m) pyb::arg("E"), pyb::arg("X"), pyb::arg("max_iters") = 5, - pyb::arg("eps") = 1e-10); + pyb::arg("eps") = 1e-10, + "|#element dims| x |X.cols()| matrix of reference positions associated with domain points " + "X in corresponding elements E "); } } // namespace fem diff --git a/bindings/pypbat/fem/Laplacian.cpp b/bindings/pypbat/fem/Laplacian.cpp index f755991..134c030 100644 --- a/bindings/pypbat/fem/Laplacian.cpp +++ b/bindings/pypbat/fem/Laplacian.cpp @@ -64,7 +64,12 @@ void BindLaplacian(pybind11::module& m) pyb::arg("detJe"), pyb::arg("GNe"), pyb::arg("dims") = 1, - pyb::arg("quadrature_order") = 1) + pyb::arg("quadrature_order") = 1, + "Construct the symmetric part of the Laplacian operator on mesh mesh, using " + "precomputed jacobian determinants detJe and shape function gradients GNe evaluated at " + "quadrature points given by the quadrature rule of order quadrature_order. The " + "discretization is based on Galerkin projection. The dimensions dims can be set to " + "accommodate vector-valued functions.") .def_property( "dims", [](Laplacian const& L) { return L.dims(); }, @@ -76,7 +81,8 @@ void BindLaplacian(pybind11::module& m) [](Laplacian const& L) { return L.ElementLaplacians(); }, [](Laplacian& L, Eigen::Ref const& deltaE) { L.ElementLaplacians() = deltaE; - }) + }, + "|#element nodes|x|#element nodes * #elements| matrix of element Laplacians") .def_property_readonly("shape", &Laplacian::Shape) .def("to_matrix", &Laplacian::ToMatrix); } diff --git a/bindings/pypbat/fem/MassMatrix.cpp b/bindings/pypbat/fem/MassMatrix.cpp index 25bcd43..171999c 100644 --- a/bindings/pypbat/fem/MassMatrix.cpp +++ b/bindings/pypbat/fem/MassMatrix.cpp @@ -64,9 +64,14 @@ void BindMassMatrix(pybind11::module& m) pyb::init const&, Scalar, int, int>(), pyb::arg("mesh"), pyb::arg("detJe"), - pyb::arg("rho"), - pyb::arg("dims"), - pyb::arg("quadrature_order") = 1) + pyb::arg("rho") = 1., + pyb::arg("dims") = 1, + pyb::arg("quadrature_order") = 1, + "Construct the mass matrix operator on mesh mesh, using " + "precomputed jacobian determinants detJe evaluated at " + "quadrature points given by the quadrature rule of order quadrature_order. The " + "dimensions dims can be set to accommodate vector-valued functions. rho is a uniform " + "mass density.") .def( pyb::init< Mesh const&, @@ -76,9 +81,14 @@ void BindMassMatrix(pybind11::module& m) int>(), pyb::arg("mesh"), pyb::arg("detJe"), - pyb::arg("rho"), - pyb::arg("dims"), - pyb::arg("quadrature_order") = 1) + pyb::arg("rho") = 1., + pyb::arg("dims") = 1, + pyb::arg("quadrature_order") = 1, + "Construct the mass matrix operator on mesh mesh, using " + "precomputed jacobian determinants detJe evaluated at " + "quadrature points given by the quadrature rule of order quadrature_order. The " + "dimensions dims can be set to accommodate vector-valued functions. rho is a piecewise " + "constant (per element) mass density.") .def_property( "dims", [](MassMatrix const& L) { return L.dims(); }, @@ -90,7 +100,8 @@ void BindMassMatrix(pybind11::module& m) [](MassMatrix const& M) { return M.ElementMassMatrices(); }, [](MassMatrix& M, Eigen::Ref const& ME) { M.ElementMassMatrices() = ME; - }) + }, + "|#element nodes| x |#elements nodes * #elements| matrix of element mass matrices") .def_property_readonly("shape", &MassMatrix::Shape) .def("to_matrix", &MassMatrix::ToMatrix); } diff --git a/bindings/pypbat/fem/ShapeFunctions.cpp b/bindings/pypbat/fem/ShapeFunctions.cpp index 4ccd430..016e8a5 100644 --- a/bindings/pypbat/fem/ShapeFunctions.cpp +++ b/bindings/pypbat/fem/ShapeFunctions.cpp @@ -26,7 +26,9 @@ void BindShapeFunctions(pybind11::module& m) return GN; }, pyb::arg("mesh"), - pyb::arg("quadrature_order") = 1); + pyb::arg("quadrature_order") = 1, + "|#element nodes| x |#dims * #quad.pts. * #elements| matrix of shape functions at each " + "element quadrature point"); m.def( "shape_function_gradients_at", @@ -41,7 +43,9 @@ void BindShapeFunctions(pybind11::module& m) }, pyb::arg("mesh"), pyb::arg("E"), - pyb::arg("Xi")); + pyb::arg("Xi"), + "|#element nodes| x |E.size() * mesh.dims| nodal shape function gradients at reference " + "points Xi"); m.def( "shape_function_matrix", @@ -56,7 +60,8 @@ void BindShapeFunctions(pybind11::module& m) return N; }, pyb::arg("mesh"), - pyb::arg("quadrature_order") = 1); + pyb::arg("quadrature_order") = 1, + "|#elements * #quad.pts.| x |#nodes| shape function matrix"); m.def( "shape_functions_at", @@ -69,7 +74,8 @@ void BindShapeFunctions(pybind11::module& m) return N; }, pyb::arg("mesh"), - pyb::arg("Xi")); + pyb::arg("Xi"), + "|#element nodes| x |Xi.cols()| matrix of nodal shape functions at reference points Xi"); } } // namespace fem diff --git a/bindings/pypbat/fem/bindings.py b/bindings/pypbat/fem/bindings.py deleted file mode 100644 index 19836e9..0000000 --- a/bindings/pypbat/fem/bindings.py +++ /dev/null @@ -1,1396 +0,0 @@ -autogendir = "gen" - - -def mesh_types_of(max_order=3): - elements = ["Line", "Triangle", - "Quadrilateral", "Tetrahedron", "Hexahedron"] - ldims = [1, 2, 2, 3, 3] - udims = 3 - _mesh_types = [] - for e, element in enumerate(elements): - for order in range(1, max_order+1): - for dims in range(ldims[e], udims+1): - mesh_type = f"""pbat::fem::Mesh,{dims}>""" - mesh_type_py = f"""Mesh_{element.lower()}_Order_{ - order}_Dims_{dims}""" - includes = f"""#include \n#include \n""" - mesh_type_filename = f"M{element[:3].lower()}{order}_{dims}" - _mesh_types.append( - (mesh_type, mesh_type_py, includes, mesh_type_filename)) - return _mesh_types - - -def bind_gradient(mesh_types: list, max_qorder=6): - headers = [] - sources = [] - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - for qorder in range(1, max_qorder+1): - header = f"{autogendir}/Grad{qorder}{mesh_type_filename}.h" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_GRADIENT_{qorder}_{mesh_type_py}_H -#define PYPBAT_FEM_GRADIENT_{qorder}_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindGradient_{qorder}_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_GRADIENT_{qorder}_{mesh_type_py}_H - """ - file.write(code) - - source = f"{autogendir}/Grad{qorder}{mesh_type_filename}.cpp" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "Grad{qorder}{mesh_type_filename}.h" - -{mesh_includes} -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindGradient_{qorder}_{mesh_type_py}(pybind11::module& m) -{{ - namespace pyb = pybind11; - auto constexpr QuadratureOrder = {qorder}; - using MeshType = {mesh_type}; - using GradientMatrixType = pbat::fem::Gradient; - std::string const className = "GradientMatrix_QuadratureOrder_{qorder}_{mesh_type_py}"; - pyb::class_(m, className.data()) - .def( - pyb::init([](MeshType const& mesh, - Eigen::Ref const& GNe) {{ - return GradientMatrixType(mesh, GNe); - }}), - pyb::arg("mesh"), - pyb::arg("GNe")) - .def_property_readonly_static( - "dims", - [](pyb::object /*self*/) {{ return GradientMatrixType::kDims; }}) - .def_property_readonly_static( - "order", - [](pyb::object /*self*/) {{ return GradientMatrixType::kOrder; }}) - .def_property_readonly_static( - "quadrature_order", - [](pyb::object /*self*/) {{ return GradientMatrixType::kQuadratureOrder; }}) - .def("rows", &GradientMatrixType::OutputDimensions) - .def("cols", &GradientMatrixType::InputDimensions) - .def("to_matrix", &GradientMatrixType::ToMatrix); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat""" - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_GRADIENT_H -#define PYPBAT_FEM_GRADIENT_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindGradient(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_GRADIENT_H -""" - - includes = "\n".join([f"#include \"Grad{qorder}{mesh_type_filename}.h\"" for qorder in range( - 1, max_qorder+1) for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - bind_calls = "\n".join([f"BindGradient_{qorder}_{mesh_type_py}(m);" for qorder in range( - 1, max_qorder+1) for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - source = f""" -#include "Gradient.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindGradient(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/Gradient.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/Gradient.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/Gradient.h") - sources.append(f"{autogendir}/Gradient.cpp") - return (headers, sources) - - -def bind_hyper_elastic_potential(mesh_types: list, max_qorder: int = 6): - headers = [] - sources = [] - - psi_types = [ - ("SaintVenantKirchhoffEnergy", "StVk", "StVk"), ("StableNeoHookeanEnergy", "StableNeoHookean", "SNH")] - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - # Order 3 hexahedron element has too many DOFs to be stack allocated - # by Eigen in the fem::HessianWrtDofs calls, so skip it. - if mesh_type_py == "Mesh_hexahedron_Order_3_Dims_3": - continue - - for (psi_type, psi_type_py, psi_type_short) in psi_types: - for qorder in range(1, max_qorder+1): - header = f"""{ - autogendir}/HEP{psi_type_short}_{qorder}{mesh_type_filename}.h""" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_HYPER_ELASTIC_POTENTIAL_{psi_type_py}_{qorder}_{mesh_type_py}_H -#define PYPBAT_FEM_HYPER_ELASTIC_POTENTIAL_{psi_type_py}_{qorder}_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindHyperElasticPotential_{psi_type_py}_{qorder}_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_HYPER_ELASTIC_POTENTIAL_{psi_type_py}_{qorder}_{mesh_type_py}_H -""" - file.write(code) - - source = f"""{ - autogendir}/HEP{psi_type_short}_{qorder}{mesh_type_filename}.cpp""" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "HEP{psi_type_short}_{qorder}{mesh_type_filename}.h" - -{mesh_includes} -#include -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindHyperElasticPotential_{psi_type_py}_{qorder}_{mesh_type_py}(pybind11::module& m) -{{ - namespace pyb = pybind11; - auto constexpr QuadratureOrder = {qorder}; - using MeshType = {mesh_type}; - using ElementType = typename MeshType::ElementType; - using HyperElasticEnergy = pbat::physics::{psi_type}; - using ElasticPotentialType = pbat::fem:: - HyperElasticPotential; - std::string const className = - "HyperElasticPotential_{psi_type_py}_QuadratureOrder_{qorder}_Dims_" + std::to_string(MeshType::kDims) + - "_{mesh_type_py}"; - pyb::class_(m, className.data()) - .def( - pyb::init([](MeshType const& mesh, - Eigen::Ref const& detJe, - Eigen::Ref const& GNe, - Eigen::Ref const& Y, - Eigen::Ref const& nu) {{ - return ElasticPotentialType(mesh, detJe, GNe, Y, nu); - }}), - pyb::arg("mesh"), - pyb::arg("detJe"), - pyb::arg("GNe"), - pyb::arg("Y"), - pyb::arg("nu")) - .def( - pyb::init([](MeshType const& mesh, - Eigen::Ref const& detJe, - Eigen::Ref const& GNe, - Eigen::Ref const& x, - Eigen::Ref const& Y, - Eigen::Ref const& nu) {{ - return ElasticPotentialType(mesh, detJe, GNe, x, Y, nu); - }}), - pyb::arg("mesh"), - pyb::arg("detJe"), - pyb::arg("GNe"), - pyb::arg("x"), - pyb::arg("Y"), - pyb::arg("nu")) - .def_property_readonly_static( - "dims", - [](pyb::object /*self*/) {{ return ElasticPotentialType::kDims; }}) - .def_property_readonly_static( - "order", - [](pyb::object /*self*/) {{ return ElasticPotentialType::kOrder; }}) - .def_property_readonly_static( - "quadrature_order", - [](pyb::object /*self*/) {{ - return ElasticPotentialType::kQuadratureOrder; - }}) - .def( - "precompute_hessian_sparsity", - &ElasticPotentialType::PrecomputeHessianSparsity) - .def( - "compute_element_elasticity", - [](ElasticPotentialType& U, - Eigen::Ref const& x, - bool bWithGradient, - bool bWithHessian) {{ U.ComputeElementElasticity(x, bWithGradient, bWithHessian); }}, - pyb::arg("x"), - pyb::arg("grad") = true, - pyb::arg("hess") = true) - .def("to_matrix", &ElasticPotentialType::ToMatrix) - .def("to_vector", &ElasticPotentialType::ToVector) - .def("eval", &ElasticPotentialType::Eval) - .def("rows", &ElasticPotentialType::OutputDimensions) - .def("cols", &ElasticPotentialType::InputDimensions) - .def_readwrite("mue", &ElasticPotentialType::mue) - .def_readwrite("lambdae", &ElasticPotentialType::lambdae) - .def_readonly("He", &ElasticPotentialType::He) - .def_readonly("Ge", &ElasticPotentialType::Ge) - .def_readonly("Ue", &ElasticPotentialType::Ue); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_HYPER_ELASTIC_POTENTIAL_H -#define PYPBAT_FEM_HYPER_ELASTIC_POTENTIAL_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindHyperElasticPotential(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_HYPER_ELASTIC_POTENTIAL_H -""" - - includes = "\n".join([f"#include \"HEP{psi_type_short}_{qorder}{mesh_type_filename}.h\"" - for (psi_type, psi_type_py, psi_type_short) in psi_types - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types - if mesh_type_py != "Mesh_hexahedron_Order_3_Dims_3"]) - - bind_calls = "\n".join([f"BindHyperElasticPotential_{psi_type_py}_{qorder}_{mesh_type_py}(m);" - for (psi_type, psi_type_py, psi_type_short) in psi_types - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types - if mesh_type_py != "Mesh_hexahedron_Order_3_Dims_3"]) - - source = f""" -#include "HyperElasticPotential.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindHyperElasticPotential(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/HyperElasticPotential.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/HyperElasticPotential.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/HyperElasticPotential.h") - sources.append(f"{autogendir}/HyperElasticPotential.cpp") - return (headers, sources) - - -def bind_jacobian(mesh_types: list, max_qorder: int = 6): - headers = [] - sources = [] - - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - header = f"{autogendir}/Jac{mesh_type_filename}.h" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_JACOBIAN_{mesh_type_py}_H -#define PYPBAT_FEM_JACOBIAN_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindJacobian_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_JACOBIAN_{mesh_type_py}_H -""" - file.write(code) - - source = f"{autogendir}/Jac{mesh_type_filename}.cpp" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "Jac{mesh_type_filename}.h" - -#include -#include -#include -{mesh_includes} -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindJacobian_{mesh_type_py}(pybind11::module& m) -{{ - namespace pyb = pybind11; - using MeshType = {mesh_type}; - auto constexpr kMaxQuadratureOrder = {max_qorder}; - auto const throw_bad_quad_order = [&](int qorder) {{ - std::string const what = fmt::format( - "Invalid quadrature order={{}}, supported orders are [1,{{}}]", - qorder, - kMaxQuadratureOrder); - throw std::invalid_argument(what); - }}; - std::string const jacobianDeterminantsName = "jacobian_determinants_{mesh_type_py}"; - m.def( - jacobianDeterminantsName.data(), - [&](MeshType const& mesh, int qorder) -> MatrixX {{ - MatrixX R; - pbat::common::ForRange<1, kMaxQuadratureOrder + 1>([&]() {{ - if (qorder == QuadratureOrder) - {{ - R = pbat::fem::DeterminantOfJacobian(mesh); - }} - }}); - if (R.size() == 0) - throw_bad_quad_order(qorder); - return R; - }}, - pyb::arg("mesh"), - pyb::arg("quadrature_order")); - - std::string const referencePositionsName = "reference_positions_{mesh_type_py}"; - m.def( - referencePositionsName.data(), - [&](MeshType const& mesh, - Eigen::Ref const& E, - Eigen::Ref const& X, - int maxIterations, - Scalar eps) -> MatrixX {{ - return pbat::fem::ReferencePositions(mesh, E, X, maxIterations, eps); - }}, - pyb::arg("mesh"), - pyb::arg("E"), - pyb::arg("X"), - pyb::arg("max_iterations"), - pyb::arg("epsilon")); - - std::string const innerProductWeights = "inner_product_weights_{mesh_type_py}"; - m.def( - innerProductWeights.data(), - [&](MeshType const& mesh, int qorder) -> MatrixX {{ - MatrixX R; - pbat::common::ForRange<1, kMaxQuadratureOrder + 1>([&]() {{ - if (qorder == QuadratureOrder) - {{ - R = pbat::fem::InnerProductWeights(mesh); - }} - }}); - if (R.size() == 0) - throw_bad_quad_order(qorder); - return R; - }}, - pyb::arg("mesh"), - pyb::arg("quadrature_order")); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_JACOBIAN_H -#define PYPBAT_FEM_JACOBIAN_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindJacobian(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_JACOBIAN_H -""" - - includes = "\n".join([f"#include \"Jac{mesh_type_filename}.h\"" - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - bind_calls = "\n".join([f"BindJacobian_{mesh_type_py}(m);" - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - source = f""" -#include "Jacobian.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindJacobian(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/Jacobian.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/Jacobian.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/Jacobian.h") - sources.append(f"{autogendir}/Jacobian.cpp") - return (headers, sources) - - -def bind_laplacian_matrix(mesh_types: list, max_qorder: int = 6): - headers = [] - sources = [] - - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - for qorder in range(1, max_qorder+1): - header = f"{autogendir}/LM{qorder}{mesh_type_filename}.h" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_LAPLACIAN_MATRIX_{qorder}_{mesh_type_py}_H -#define PYPBAT_FEM_LAPLACIAN_MATRIX_{qorder}_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindLaplacianMatrix_{qorder}_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_LAPLACIAN_MATRIX_{qorder}_{mesh_type_py}_H -""" - file.write(code) - - source = f"{autogendir}/LM{qorder}{mesh_type_filename}.cpp" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "LM{qorder}{mesh_type_filename}.h" - -{mesh_includes} -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindLaplacianMatrix_{qorder}_{mesh_type_py}(pybind11::module& m) -{{ - namespace pyb = pybind11; - using MeshType = {mesh_type}; - auto constexpr QuadratureOrder = {qorder}; - using LaplacianMatrixType = - pbat::fem::SymmetricLaplacianMatrix; - std::string const className = "SymmetricLaplacianMatrix_QuadratureOrder_{qorder}_{mesh_type_py}"; - pyb::class_(m, className.data()) - .def( - pyb::init([](MeshType const& mesh, - Eigen::Ref const& detJe, - Eigen::Ref const& GNe) {{ - return LaplacianMatrixType(mesh, detJe, GNe); - }}), - pyb::arg("mesh"), - pyb::arg("detJe"), - pyb::arg("GNe")) - .def_property_readonly_static( - "order", - [](pyb::object /*self*/) {{ return LaplacianMatrixType::kOrder; }}) - .def_property_readonly_static( - "quadrature_order", - [](pyb::object /*self*/) {{ return LaplacianMatrixType::kQuadratureOrder; }}) - .def("to_matrix", &LaplacianMatrixType::ToMatrix) - .def("rows", &LaplacianMatrixType::OutputDimensions) - .def("cols", &LaplacianMatrixType::InputDimensions) - .def_readonly("deltae", &LaplacianMatrixType::deltaE) - .def_readwrite("dims", &LaplacianMatrixType::dims); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat - """ - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_LAPLACIAN_MATRIX_H -#define PYPBAT_FEM_LAPLACIAN_MATRIX_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindLaplacianMatrix(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_LAPLACIAN_MATRIX_H -""" - - includes = "\n".join([f"#include \"LM{qorder}{mesh_type_filename}.h\"" - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - bind_calls = "\n".join([f"BindLaplacianMatrix_{qorder}_{mesh_type_py}(m);" - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - source = f""" -#include "LaplacianMatrix.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindLaplacianMatrix(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/LaplacianMatrix.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/LaplacianMatrix.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/LaplacianMatrix.h") - sources.append(f"{autogendir}/LaplacianMatrix.cpp") - return (headers, sources) - - -def bind_load_vector(mesh_types: list, max_qorder: int = 3): - headers = [] - sources = [] - - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - for qorder in range(1, max_qorder+1): - header = f"{autogendir}/LV{qorder}{mesh_type_filename}.h" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_LOAD_VECTOR_{qorder}_{mesh_type_py}_H -#define PYPBAT_FEM_LOAD_VECTOR_{qorder}_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindLoadVector_{qorder}_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_LOAD_VECTOR_{qorder}_{mesh_type_py}_H -""" - file.write(code) - - source = f"{autogendir}/LV{qorder}{mesh_type_filename}.cpp" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "LV{qorder}{mesh_type_filename}.h" - -{mesh_includes} -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindLoadVector_{qorder}_{mesh_type_py}(pybind11::module& m) -{{ - namespace pyb = pybind11; - using MeshType = {mesh_type}; - auto constexpr QuadratureOrder = {qorder}; - using LoadVectorType = pbat::fem::LoadVector; - std::string const className = - "LoadVector_QuadratureOrder_{qorder}_{mesh_type_py}"; - pyb::class_(m, className.data()) - .def( - pyb::init([](MeshType const& mesh, - Eigen::Ref const& detJe, - Eigen::Ref const& fe, - int dims) {{ - return LoadVectorType(mesh, detJe, fe, dims); - }}), - pyb::arg("mesh"), - pyb::arg("detJe"), - pyb::arg("fe"), - pyb::arg("dims")) - .def_property_readonly_static( - "order", - [](pyb::object /*self*/) {{ return LoadVectorType::kOrder; }}) - .def_property_readonly_static( - "quadrature_order", - [](pyb::object /*self*/) {{ return LoadVectorType::kQuadratureOrder; }}) - .def_readonly("fe", &LoadVectorType::fe) - .def("to_vector", &LoadVectorType::ToVector) - .def( - "set_load", - [](LoadVectorType& f, Eigen::Ref const& fe) {{ - f.SetLoad(fe); - }}, - pyb::arg("fe")) - .def_readonly("N", &LoadVectorType::N) - .def_readwrite("dims", &LoadVectorType::dims); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_LOAD_VECTOR_H -#define PYPBAT_FEM_LOAD_VECTOR_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindLoadVector(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_LOAD_VECTOR_H -""" - - includes = "\n".join([f"#include \"LV{qorder}{mesh_type_filename}.h\"" - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - bind_calls = "\n".join([f"BindLoadVector_{qorder}_{mesh_type_py}(m);" - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - source = f""" -#include "LoadVector.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindLoadVector(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/LoadVector.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/LoadVector.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/LoadVector.h") - sources.append(f"{autogendir}/LoadVector.cpp") - return (headers, sources) - - -def bind_mass_matrix(mesh_types: list, max_qorder: int = 6): - headers = [] - sources = [] - - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - for qorder in range(1, max_qorder+1): - header = f"{autogendir}/MM{qorder}{mesh_type_filename}.h" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_MASS_MATRIX_{qorder}_{mesh_type_py}_H -#define PYPBAT_FEM_MASS_MATRIX_{qorder}_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindMassMatrix_{qorder}_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_MASS_MATRIX_{qorder}_{mesh_type_py}_H -""" - file.write(code) - - source = f"{autogendir}/MM{qorder}{mesh_type_filename}.cpp" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "MM{qorder}{mesh_type_filename}.h" - -{mesh_includes} -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindMassMatrix_{qorder}_{mesh_type_py}(pybind11::module& m) -{{ -namespace pyb = pybind11; -using MeshType = {mesh_type}; -auto constexpr kQuadratureOrder = {qorder}; -using MassMatrixType = pbat::fem::MassMatrix; -std::string const className = - "MassMatrix_QuadratureOrder_{qorder}_{mesh_type_py}"; -pyb::class_(m, className.data()) - .def( - pyb::init([](MeshType const& mesh, Eigen::Ref const& detJe) {{ - return MassMatrixType(mesh, detJe); - }}), - pyb::arg("mesh"), - pyb::arg("detJe")) - .def( - pyb::init([](MeshType const& mesh, - Eigen::Ref const& detJe, - Scalar rho) {{ return MassMatrixType(mesh, detJe, rho); }}), - pyb::arg("mesh"), - pyb::arg("detJe"), - pyb::arg("rho")) - .def( - pyb::init( - [](MeshType const& mesh, - Eigen::Ref const& detJe, - VectorX const& rhoe) {{ return MassMatrixType(mesh, detJe, rhoe); }}), - pyb::arg("mesh"), - pyb::arg("detJe"), - pyb::arg("rhoe")) - .def_property_readonly_static( - "order", - [](pyb::object /*self*/) {{ return MassMatrixType::kOrder; }}) - .def_property_readonly_static( - "quadrature_order", - [](pyb::object /*self*/) {{ return MassMatrixType::kQuadratureOrder; }}) - .def_readonly("Me", &MassMatrixType::Me) - .def("rows", &MassMatrixType::OutputDimensions) - .def("cols", &MassMatrixType::InputDimensions) - .def("to_matrix", &MassMatrixType::ToMatrix) - .def( - "compute_element_mass_matrices", - [](MassMatrixType& M, VectorX const& rhoe) {{ - M.ComputeElementMassMatrices(rhoe); - }}, - pyb::arg("rhoe")) - .def_readwrite("dims", &MassMatrixType::dims); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_MASS_MATRIX_H -#define PYPBAT_FEM_MASS_MATRIX_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindMassMatrix(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_MASS_MATRIX_H -""" - - includes = "\n".join([f"#include \"MM{qorder}{mesh_type_filename}.h\"" - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - bind_calls = "\n".join([f"BindMassMatrix_{qorder}_{mesh_type_py}(m);" - for qorder in range(1, max_qorder+1) - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - source = f""" -#include "MassMatrix.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindMassMatrix(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/MassMatrix.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/MassMatrix.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/MassMatrix.h") - sources.append(f"{autogendir}/MassMatrix.cpp") - return (headers, sources) - - -def bind_mesh(mesh_types: list): - headers = [] - sources = [] - - max_qorder = 8 - - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - header = f"{autogendir}/{mesh_type_filename}.h" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_MESH_{mesh_type_py}_H -#define PYPBAT_FEM_MESH_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindMesh_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_MESH_{mesh_type_py}_H -""" - file.write(code) - - source = f"{autogendir}/{mesh_type_filename}.cpp" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "{mesh_type_filename}.h" - -{mesh_includes} -#include -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindMesh_{mesh_type_py}(pybind11::module& m) -{{ - namespace pyb = pybind11; - using MeshType = {mesh_type}; - std::string const elementTypeName = "{str(mesh_type_py).split("_")[1].capitalize()}"; - std::string const className = "{mesh_type_py}"; - - auto constexpr kMaxQuadratureOrder = {max_qorder}; - auto const throw_bad_quad_order = [&](int qorder) {{ - std::string const what = fmt::format( - "Invalid quadrature order={{}}, supported orders are [1,{{}}]", - qorder, - kMaxQuadratureOrder); - throw std::invalid_argument(what); - }}; - - pyb::class_(m, className.data()) - .def(pyb::init<>()) - .def( - pyb:: - init const&, Eigen::Ref const&>(), - pyb::arg("V"), - pyb::arg("C")) - .def("quadrature_points", - [&](MeshType const& self, int qorder) -> MatrixX {{ - MatrixX R; - pbat::common::ForRange<1, kMaxQuadratureOrder + 1>([&]() {{ - if (qorder == QuadratureOrder) - {{ - R = self.QuadraturePoints(); - }} - }}); - if (R.size() == 0) - throw_bad_quad_order(qorder); - return R; - }}, - pyb::arg("quadrature_order")) - .def( - "quadrature_weights", - [&](MeshType const& self, int qorder) -> VectorX {{ - VectorX R; - pbat::common::ForRange<1, kMaxQuadratureOrder + 1>([&]() {{ - if (qorder == QuadratureOrder) - {{ - R = self.QuadratureWeights(); - }} - }}); - if (R.size() == 0) - throw_bad_quad_order(qorder); - return R; - }}, - pyb::arg("quadrature_order")) - .def_property_readonly_static( - "dims", - [](pyb::object /*self*/) {{ return MeshType::kDims; }}) - .def_property_readonly_static( - "order", - [](pyb::object /*self*/) {{ return MeshType::kOrder; }}) - .def_property_readonly_static( - "element", - [=](pyb::object /*self*/) {{ return elementTypeName; }}) - .def_readwrite("E", &MeshType::E) - .def_readwrite("X", &MeshType::X); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_MESH_H -#define PYPBAT_FEM_MESH_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindMesh(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_MESH_H -""" - - includes = "\n".join([f"#include \"{mesh_type_filename}.h\"" - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - bind_calls = "\n".join([f"BindMesh_{mesh_type_py}(m);" - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - source = f""" -#include "Mesh.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindMesh(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/Mesh.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/Mesh.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/Mesh.h") - sources.append(f"{autogendir}/Mesh.cpp") - return (headers, sources) - - -def bind_shape_functions(mesh_types: list, max_qorder: int = 6): - headers = [] - sources = [] - - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types: - header = f"{autogendir}/SF{mesh_type_filename}.h" - with open(header, 'w', encoding="utf-8") as file: - code = f""" -#ifndef PYPBAT_FEM_SHAPE_FUNCTIONS_{mesh_type_py}_H -#define PYPBAT_FEM_SHAPE_FUNCTIONS_{mesh_type_py}_H - -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindShapeFunctions_{mesh_type_py}(pybind11::module& m); - -}} // namespace fem -}} // namespace py -}} // namespace pbat - -#endif // PYPBAT_FEM_SHAPE_FUNCTIONS_{mesh_type_py}_H -""" - file.write(code) - - source = f"{autogendir}/SF{mesh_type_filename}.cpp" - with open(source, 'w', encoding="utf-8") as file: - code = f""" -#include "SF{mesh_type_filename}.h" - -#include -#include -#include -{mesh_includes} -#include -#include -#include - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindShapeFunctions_{mesh_type_py}(pybind11::module& m) -{{ - namespace pyb = pybind11; - using MeshType = {mesh_type}; - auto constexpr kMaxQuadratureOrder = {max_qorder}; - auto const throw_bad_quad_order = [&](int qorder) {{ - std::string const what = fmt::format( - "Invalid quadrature order={{}}, supported orders are [1,{{}}]", - qorder, - kMaxQuadratureOrder); - throw std::invalid_argument(what); - }}; - std::string const meshTypeName = "{mesh_type_py}"; - - std::string const integratedShapeFunctionsName = - "integrated_shape_functions_" + meshTypeName; - m.def( - integratedShapeFunctionsName.data(), - [&](MeshType const& mesh, - Eigen::Ref const& detJe, - int qorder) -> MatrixX {{ - MatrixX R; - pbat::common::ForRange<1, kMaxQuadratureOrder + 1>([&]() {{ - if (qorder == QuadratureOrder) - {{ - R = pbat::fem::IntegratedShapeFunctions( - mesh, - detJe); - }} - }}); - if (R.size() == 0) - throw_bad_quad_order(qorder); - return R; - }}, - pyb::arg("mesh"), - pyb::arg("detJe"), - pyb::arg("quadrature_order")); - std::string const shapeFunctionGradientsName = "shape_function_gradients_" + meshTypeName; - m.def( - shapeFunctionGradientsName.data(), - [&](MeshType const& mesh, int qorder) -> MatrixX {{ - MatrixX R; - pbat::common::ForRange<1, kMaxQuadratureOrder + 1>([&]() {{ - if (qorder == QuadratureOrder) - {{ - R = pbat::fem::ShapeFunctionGradients(mesh); - }} - }}); - if (R.size() == 0) - throw_bad_quad_order(qorder); - return R; - }}, - pyb::arg("mesh"), - pyb::arg("quadrature_order")); - - std::string const shapeFunctionsAtName = "shape_functions_at_" + meshTypeName; - m.def( - shapeFunctionsAtName.data(), - [&]([[maybe_unused]] MeshType const& mesh, Eigen::Ref const& Xi) -> MatrixX {{ - return pbat::fem::ShapeFunctionsAt(Xi); - }}, - pyb::arg("mesh"), - pyb::arg("Xi")); - - std::string const shapeFunctionGradientsAtName = "shape_function_gradients_at_" + meshTypeName; - m.def( - shapeFunctionGradientsAtName.data(), - [&](MeshType const& mesh, - Eigen::Ref const& E, - Eigen::Ref const& Xi) -> MatrixX {{ - return pbat::fem::ShapeFunctionGradientsAt(mesh, E, Xi); - }}, - pyb::arg("mesh"), - pyb::arg("E"), - pyb::arg("Xi")); - - std::string const shapeFunctionMatrixName = "shape_function_matrix_" + meshTypeName; - m.def( - shapeFunctionMatrixName.data(), - [&](MeshType const& mesh, int qorder) -> CSRMatrix {{ - CSRMatrix N; - pbat::common::ForRange<1, kMaxQuadratureOrder + 1>([&]() {{ - if (qorder == QuadratureOrder) - {{ - N = pbat::fem::ShapeFunctionMatrix(mesh); - }} - }}); - if (N.size() == 0) - throw_bad_quad_order(qorder); - return N; - }}, - pyb::arg("mesh"), - pyb::arg("quadrature_order")); -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - file.write(code) - - headers.append(header) - sources.append(source) - - header = """ -#ifndef PYPBAT_FEM_SHAPE_FUNCTIONS_H -#define PYPBAT_FEM_SHAPE_FUNCTIONS_H - -#include - -namespace pbat { -namespace py { -namespace fem { - -void BindShapeFunctions(pybind11::module& m); - -} // namespace fem -} // namespace py -} // namespace pbat - -#endif // PYPBAT_FEM_SHAPE_FUNCTIONS_H -""" - - includes = "\n".join([f"#include \"SF{mesh_type_filename}.h\"" - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - bind_calls = "\n".join([f"BindShapeFunctions_{mesh_type_py}(m);" - for mesh_type, mesh_type_py, mesh_includes, mesh_type_filename in mesh_types]) - - source = f""" -#include "ShapeFunctions.h" - -{includes} - -namespace pbat {{ -namespace py {{ -namespace fem {{ - -void BindShapeFunctions(pybind11::module& m) -{{ - {bind_calls} -}} - -}} // namespace fem -}} // namespace py -}} // namespace pbat -""" - - with open(f"{autogendir}/ShapeFunctions.h", 'w', encoding="utf-8") as fheader: - fheader.write(header) - with open(f"{autogendir}/ShapeFunctions.cpp", 'w', encoding="utf-8") as fsource: - fsource.write(source) - - headers.append(f"{autogendir}/ShapeFunctions.h") - sources.append(f"{autogendir}/ShapeFunctions.cpp") - return (headers, sources) - - -if __name__ == "__main__": - import argparse - parser = argparse.ArgumentParser( - prog="PBAT FEM Python Bindings Generator", - description="""Generates Python binding code for pbat::fem""", - ) - parser.add_argument( - "-c", - "--cmake", - help="Prints generated files to CMake via stdout", - type=bool, - dest="cmake") - parser.add_argument( - "-t", - "--type", - help="The type whose instantiations should be generated", - type=str, - dest="type") - args = parser.parse_args() - - import os - if not os.path.exists("gen"): - os.mkdir("gen") - - meshes = mesh_types_of(max_order=3) - headers, sources = None, None - if args.type == "Gradient": - header, sources = bind_gradient(meshes, max_qorder=5) - if args.type == "HyperElasticPotential": - headers, sources = bind_hyper_elastic_potential(meshes, max_qorder=6) - if args.type == "Jacobian": - headers, sources = bind_jacobian(meshes, max_qorder=6) - if args.type == "LaplacianMatrix": - headers, sources = bind_laplacian_matrix(meshes, max_qorder=4) - if args.type == "LoadVector": - headers, sources = bind_load_vector(meshes, max_qorder=3) - if args.type == "MassMatrix": - headers, sources = bind_mass_matrix(meshes, max_qorder=6) - if args.type == "Mesh": - headers, sources = bind_mesh(meshes) - if args.type == "ShapeFunctions": - headers, sources = bind_shape_functions(meshes, max_qorder=6) - - # headers = "\n".join([ - # "\n".join(gradient_headers), - # "\n".join(hyper_elastic_potential_headers), - # "\n".join(jacobian_headers), - # "\n".join(laplacian_matrix_headers), - # "\n".join(load_vector_headers), - # "\n".join(mass_matrix_headers), - # "\n".join(mesh_headers), - # "\n".join(shape_functions_headers) - # ]) - if args.cmake: - import sys - cmake_target_sources = ";".join(sources) - sys.stdout.write(cmake_target_sources) diff --git a/python/examples/diffusion_smoothing.py b/python/examples/diffusion_smoothing.py index 67cf50c..00ffc59 100644 --- a/python/examples/diffusion_smoothing.py +++ b/python/examples/diffusion_smoothing.py @@ -1,5 +1,4 @@ import pbatoolkit as pbat -import pbatoolkit.fem import pbatoolkit.math.linalg import igl import polyscope as ps @@ -22,15 +21,15 @@ V, F = V.astype(np.float64, order='c'), F.astype(np.int64, order='c') # Construct Galerkin laplacian, mass and gradient operators - mesh = pbat.fem.mesh( + mesh = pbat.fem.Mesh( V.T, F.T, element=pbat.fem.Element.Triangle, order=1) V, F = mesh.X.T, mesh.E.T detJeL = pbat.fem.jacobian_determinants(mesh, quadrature_order=1) GNeL = pbat.fem.shape_function_gradients(mesh, quadrature_order=1) - L = pbat.fem.laplacian_matrix( + L = pbat.fem.Laplacian( mesh, detJeL, GNeL, quadrature_order=1).to_matrix() detJeM = pbat.fem.jacobian_determinants(mesh, quadrature_order=2) - M = pbat.fem.mass_matrix(mesh, detJeM, dims=1, + M = pbat.fem.MassMatrix(mesh, detJeM, dims=1, quadrature_order=2).to_matrix() # Setup heat diffusion dt = 0.016 diff --git a/python/examples/elasticity.py b/python/examples/elasticity.py index 60e1c68..4ddcd53 100644 --- a/python/examples/elasticity.py +++ b/python/examples/elasticity.py @@ -1,5 +1,4 @@ import pbatoolkit as pbat -import pbatoolkit.fem import pbatoolkit.geometry import pbatoolkit.profiling import pbatoolkit.math.linalg @@ -29,14 +28,14 @@ # Construct FEM quantities for simulation imesh = meshio.read(args.input) V, C = imesh.points, imesh.cells_dict["tetra"] - mesh = pbat.fem.mesh( + mesh = pbat.fem.Mesh( V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=1) x = mesh.X.reshape(math.prod(mesh.X.shape), order='f') v = np.zeros(x.shape[0]) detJeM = pbat.fem.jacobian_determinants(mesh, quadrature_order=2) rho = args.rho - M = pbat.fem.mass_matrix(mesh, detJeM, rho=rho, - dims=3, quadrature_order=2).to_matrix() + M = pbat.fem.MassMatrix(mesh, detJeM, rho=rho, + dims=3, quadrature_order=2).to_matrix() Minv = pbat.math.linalg.ldlt(M) Minv.compute(M) # Could also lump the mass matrix like this @@ -48,18 +47,23 @@ # Construct load vector from gravity field detJeU = pbat.fem.jacobian_determinants(mesh, quadrature_order=1) GNeU = pbat.fem.shape_function_gradients(mesh, quadrature_order=1) + qgf = pbat.fem.inner_product_weights( + mesh, quadrature_order=1).flatten(order="F") + Qf = sp.sparse.diags_array([qgf], offsets=[0]) + Nf = pbat.fem.shape_function_matrix(mesh, quadrature_order=1) g = np.zeros(mesh.dims) g[-1] = -9.81 fe = np.tile(rho*g[:, np.newaxis], mesh.E.shape[1]) - f = pbat.fem.load_vector(mesh, detJeU, fe, quadrature_order=1).to_vector() + f = fe @ Qf @ Nf + f = f.reshape(math.prod(f.shape), order="F") a = Minv.solve(f).squeeze() # Create hyper elastic potential Y = np.full(mesh.E.shape[1], args.Y) nu = np.full(mesh.E.shape[1], args.nu) psi = pbat.fem.HyperElasticEnergy.StableNeoHookean - hep = pbat.fem.hyper_elastic_potential( - mesh, detJeU, GNeU, Y, nu, psi=psi, quadrature_order=1) + hep = pbat.fem.HyperElasticPotential( + mesh, detJeU, GNeU, Y, nu, energy=psi, quadrature_order=1) hep.precompute_hessian_sparsity() # Set Dirichlet boundary conditions @@ -80,7 +84,8 @@ # Setup linear solver Hdd = hep.to_matrix()[:, dofs].tocsr()[dofs, :] Mdd = M[:, dofs].tocsr()[dofs, :] - Addinv = pbat.math.linalg.ldlt(Hdd) + Addinv = pbat.math.linalg.ldlt( + Hdd, solver=pbat.math.linalg.SolverBackend.Eigen) Addinv.analyze(Hdd) ps.set_verbosity(0) @@ -134,8 +139,8 @@ def callback(): xtilde = x + dt*v + dt2*a xk = x for k in range(newton_maxiter): - hep.compute_element_elasticity(xk, grad=True, hess=True) - gradU, HU = hep.to_vector(), hep.to_matrix() + hep.compute_element_elasticity(xk, grad=True, hessian=True) + gradU, HU = hep.gradient(), hep.hessian() global bd, Add diff --git a/python/examples/elasticity_higher_order.py b/python/examples/elasticity_higher_order.py index 4d6ac4f..e5190f6 100644 --- a/python/examples/elasticity_higher_order.py +++ b/python/examples/elasticity_higher_order.py @@ -1,5 +1,4 @@ import pbatoolkit as pbat -import pbatoolkit.fem import pbatoolkit.geometry import pbatoolkit.profiling import pbatoolkit.math.linalg @@ -38,13 +37,13 @@ FR = FR.astype(np.int64, order='c') # Construct FEM quantities for simulation - mesh = pbat.fem.mesh( + mesh = pbat.fem.Mesh( V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=2) x = mesh.X.reshape(math.prod(mesh.X.shape), order='f') v = np.zeros(x.shape[0]) detJeM = pbat.fem.jacobian_determinants(mesh, quadrature_order=4) rho = args.rho - M = pbat.fem.mass_matrix(mesh, detJeM, rho=rho, + M = pbat.fem.MassMatrix(mesh, detJeM, rho=rho, dims=3, quadrature_order=4).to_matrix() Minv = pbat.math.linalg.ldlt(M) Minv.compute(M) @@ -53,8 +52,14 @@ g = np.zeros(mesh.dims) g[-1] = -9.81 detJef = pbat.fem.jacobian_determinants(mesh, quadrature_order=2) - fe = np.tile(rho*g[:, np.newaxis], mesh.E.shape[1]) - f = pbat.fem.load_vector(mesh, detJef, fe, quadrature_order=2).to_vector() + nquadpts = mesh.E.shape[1] * mesh.quadrature_weights(2).shape[0] + fe = np.tile(rho*g[:, np.newaxis], nquadpts) + qgf = pbat.fem.inner_product_weights( + mesh, quadrature_order=2).flatten(order="F") + Qf = sp.sparse.diags_array([qgf], offsets=[0]) + Nf = pbat.fem.shape_function_matrix(mesh, quadrature_order=2) + f = fe @ Qf @ Nf + f = f.reshape(math.prod(f.shape), order="F") a = Minv.solve(f).squeeze() # Create hyper elastic potential @@ -63,8 +68,8 @@ Y = np.full(mesh.E.shape[1], args.Y) nu = np.full(mesh.E.shape[1], args.nu) psi = pbat.fem.HyperElasticEnergy.StableNeoHookean - hep = pbat.fem.hyper_elastic_potential( - mesh, detJeU, GNeU, Y, nu, psi=psi, quadrature_order=4) + hep = pbat.fem.HyperElasticPotential( + mesh, detJeU, GNeU, Y, nu, energy=psi, quadrature_order=4) hep.precompute_hessian_sparsity() # Set Dirichlet boundary conditions @@ -118,8 +123,8 @@ def callback(): if animate or step: profiler.begin_frame("Physics") # 1 Newton step - hep.compute_element_elasticity(x, grad=True, hess=True) - gradU, HU = hep.to_vector(), hep.to_matrix() + hep.compute_element_elasticity(x, grad=True, hessian=True) + gradU, HU = hep.gradient(), hep.hessian() dt2 = dt**2 xtilde = x + dt*v + dt2*a A = M + dt2 * HU diff --git a/python/examples/heat.py b/python/examples/heat.py index 5b089fd..3177cec 100644 --- a/python/examples/heat.py +++ b/python/examples/heat.py @@ -1,5 +1,4 @@ import pbatoolkit as pbat -import pbatoolkit.fem import pbatoolkit.math.linalg import igl import polyscope as ps @@ -21,7 +20,7 @@ mesh = None if "tetra" in imesh.cells_dict.keys(): V, C = imesh.points, imesh.cells_dict["tetra"] - mesh = pbat.fem.mesh( + mesh = pbat.fem.Mesh( V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=1) if "triangle" in imesh.cells_dict.keys(): V, C = imesh.points, imesh.cells_dict["triangle"] @@ -30,7 +29,7 @@ V, C = mesh.X.T, mesh.E.T F = C - if mesh.element == "Tetrahedron": + if mesh.element == pbat.fem.Element.Tetrahedron: F = igl.boundary_facets(C) F[:, :2] = np.roll(F[:, :2], shift=1, axis=1) @@ -43,7 +42,7 @@ NM = pbat.fem.shape_function_matrix(mesh, quadrature_order=2) M = NM.T @ QM @ NM GNeL = pbat.fem.shape_function_gradients(mesh, quadrature_order=1) - G = pbat.fem.gradient_matrix( + G = pbat.fem.Gradient( mesh, GNeL, quadrature_order=1).to_matrix() qgL = pbat.fem.inner_product_weights( mesh, quadrature_order=1).flatten(order="F") @@ -110,17 +109,17 @@ def callback(): cn.set_color((0, 0, 0)) cn.set_radius(0.002) vm = ps.get_volume_mesh( - "model") if mesh.element == "Tetrahedron" else ps.get_surface_mesh("model") + "model") if mesh.element == pbat.fem.Element.Tetrahedron else ps.get_surface_mesh("model") vm.add_scalar_quantity("heat", u, cmap="reds") vm.add_scalar_quantity("distance", phi, cmap="reds", enabled=True) - grad_defined_on = "cells" if mesh.element == "Tetrahedron" else "faces" + grad_defined_on = "cells" if mesh.element == pbat.fem.Element.Tetrahedron else "faces" vm.add_vector_quantity("normalized heat grad", gradu, defined_on=grad_defined_on) vm.add_scalar_quantity("div unit gradient", divGu) - if mesh.element == "Tetrahedron": + if mesh.element == pbat.fem.Element.Tetrahedron: ps.register_volume_mesh("model", V, C) - if mesh.element == "Triangle": + if mesh.element == pbat.fem.Element.Triangle: ps.register_surface_mesh("model", V, C) ps.set_user_callback(callback) ps.show() diff --git a/python/examples/ipc.py b/python/examples/ipc.py index 547f50e..4754308 100644 --- a/python/examples/ipc.py +++ b/python/examples/ipc.py @@ -1,5 +1,4 @@ import pbatoolkit as pbat -import pbatoolkit.fem import pbatoolkit.geometry import pbatoolkit.profiling import pbatoolkit.math.linalg @@ -61,7 +60,7 @@ def combine(V: list, C: list): C.append(C2) V, C = combine(V, C) - mesh = pbat.fem.mesh( + mesh = pbat.fem.Mesh( V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=1) V, C = mesh.X.T, mesh.E.T @@ -74,7 +73,7 @@ def combine(V: list, C: list): detJeM = pbat.fem.jacobian_determinants(mesh, quadrature_order=2) rho = args.rho - M = pbat.fem.mass_matrix(mesh, detJeM, rho=rho, + M = pbat.fem.MassMatrix(mesh, detJeM, rho=rho, dims=3, quadrature_order=2).to_matrix() # Lump mass matrix lumpedm = M.sum(axis=0) @@ -85,18 +84,23 @@ def combine(V: list, C: list): # Construct load vector from gravity field detJeU = pbat.fem.jacobian_determinants(mesh, quadrature_order=1) GNeU = pbat.fem.shape_function_gradients(mesh, quadrature_order=1) + qgf = pbat.fem.inner_product_weights( + mesh, quadrature_order=1).flatten(order="F") + Qf = sp.sparse.diags_array([qgf], offsets=[0]) + Nf = pbat.fem.shape_function_matrix(mesh, quadrature_order=1) g = np.zeros(mesh.dims) g[-1] = -9.81 fe = np.tile(rho*g[:, np.newaxis], mesh.E.shape[1]) - f = pbat.fem.load_vector(mesh, detJeU, fe, quadrature_order=1).to_vector() + f = fe @ Qf @ Nf + f = f.reshape(math.prod(f.shape), order="F") a = Minv @ f # Create hyper elastic potential Y = np.full(mesh.E.shape[1], args.Y) nu = np.full(mesh.E.shape[1], args.nu) psi = pbat.fem.HyperElasticEnergy.StableNeoHookean - hep = pbat.fem.hyper_elastic_potential( - mesh, detJeU, GNeU, Y, nu, psi=psi, quadrature_order=1) + hep = pbat.fem.HyperElasticPotential( + mesh, detJeU, GNeU, Y, nu, energy=psi, quadrature_order=1) hep.precompute_hessian_sparsity() # Setup IPC contact handling @@ -183,7 +187,7 @@ def callback(): # Compute elasticity hep.compute_element_elasticity(xk, grad=True, hess=True) - U, gradU, HU = hep.eval(), hep.to_vector(), hep.to_matrix() + U, gradU, HU = hep.eval(), hep.gradient(), hep.hessian() # Compute adaptive barrier stiffness bboxdiag = ipctk.world_bbox_diagonal_length(BX) diff --git a/python/examples/laplace.py b/python/examples/laplace.py index a89a073..471417e 100644 --- a/python/examples/laplace.py +++ b/python/examples/laplace.py @@ -1,5 +1,4 @@ import pbatoolkit as pbat -import pbatoolkit.fem import pbatoolkit.math.linalg import pbatoolkit.geometry import igl @@ -11,14 +10,14 @@ def harmonic_field(V: np.ndarray, C: np.ndarray, order: int, eps: float = 0.1): # Construct order order mesh and its Laplacian - mesh = pbat.fem.mesh( + mesh = pbat.fem.Mesh( V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=order) quadrature_order = max(int(2*(order-1)), 1) detJeL = pbat.fem.jacobian_determinants( mesh, quadrature_order=quadrature_order) GNeL = pbat.fem.shape_function_gradients( mesh, quadrature_order=quadrature_order) - L = pbat.fem.laplacian_matrix( + L = pbat.fem.Laplacian( mesh, detJeL, GNeL, quadrature_order=quadrature_order).to_matrix() # Set Dirichlet boundary conditions at bottom and top of the model Xmin = mesh.X.min(axis=1) diff --git a/python/examples/modes.py b/python/examples/modes.py index 35d9ab7..db1390c 100644 --- a/python/examples/modes.py +++ b/python/examples/modes.py @@ -1,5 +1,4 @@ import pbatoolkit as pbat -import pbatoolkit.fem import numpy as np import scipy as sp import polyscope as ps @@ -32,22 +31,22 @@ def signal(w: float, v: np.ndarray, t: float, c: float, k: float): imesh = meshio.read(args.input) V, C = imesh.points, imesh.cells_dict["tetra"] - mesh = pbat.fem.mesh( + mesh = pbat.fem.Mesh( V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=1) x = mesh.X.reshape(mesh.X.shape[0]*mesh.X.shape[1], order='f') detJeM = pbat.fem.jacobian_determinants(mesh, quadrature_order=2) - M = pbat.fem.mass_matrix(mesh, detJeM, rho=args.rho, + M = pbat.fem.MassMatrix(mesh, detJeM, rho=args.rho, dims=3, quadrature_order=2).to_matrix() detJeU = pbat.fem.jacobian_determinants(mesh, quadrature_order=1) GNeU = pbat.fem.shape_function_gradients(mesh, quadrature_order=1) Y = np.full(mesh.E.shape[1], args.Y) nu = np.full(mesh.E.shape[1], args.nu) - hep = pbat.fem.hyper_elastic_potential( - mesh, detJeU, GNeU, Y, nu, psi=pbat.fem.HyperElasticEnergy.StableNeoHookean, quadrature_order=1) + hep = pbat.fem.HyperElasticPotential( + mesh, detJeU, GNeU, Y, nu, energy=pbat.fem.HyperElasticEnergy.StableNeoHookean, quadrature_order=1) hep.precompute_hessian_sparsity() hep.compute_element_elasticity(x) - U, gradU, HU = hep.eval(), hep.to_vector(), hep.to_matrix() + U, gradU, HU = hep.eval(), hep.gradient(), hep.hessian() sigma = -1e-5 leigs, Veigs = sp.sparse.linalg.eigsh(HU, k=args.modes, M=M, sigma=-1e-5, which='LM') Veigs = Veigs / sp.linalg.norm(Veigs, axis=0, keepdims=True) diff --git a/python/pbatoolkit/__init__.py b/python/pbatoolkit/__init__.py index e69de29..19f9711 100644 --- a/python/pbatoolkit/__init__.py +++ b/python/pbatoolkit/__init__.py @@ -0,0 +1,3 @@ +from ._pbat import fem as _fem + +fem = _fem \ No newline at end of file diff --git a/python/pbatoolkit/fem.py b/python/pbatoolkit/fem.py deleted file mode 100644 index 6ef64e1..0000000 --- a/python/pbatoolkit/fem.py +++ /dev/null @@ -1,301 +0,0 @@ -from ._pbat import fem as _fem -import numpy as np -from enum import Enum - - -class Element(Enum): - Line = 0 - Triangle = 1 - Quadrilateral = 2 - Tetrahedron = 3 - Hexahedron = 4 - - -def _mesh_type_name(mesh, element: str = None, order: int = None, dims: int = None): - if mesh is not None: - return f"Mesh_{mesh.element.lower()}_Order_{mesh.order}_Dims_{mesh.dims}" - class_name = f"Mesh_{element.lower()}_Order_{order}_Dims_{dims}" - return class_name - - -def mesh(V: np.ndarray, C: np.ndarray, element: Element, order: int = 1): - """Construct an FEM mesh from an input mesh - - Args: - V (np.ndarray): |#dims|x|#vertices| mesh vertex positions - C (np.ndarray): |#cell nodes|x|#cells| mesh cell indices into V - element (Element): Type of element - order (int, optional): Element shape function order. Defaults to 1. - - Returns: - The FEM mesh of the given order - """ - dims = V.shape[0] - class_name = _mesh_type_name(None, element.name, order, dims) - class_ = getattr(_fem, class_name) - return class_(V, C) - - -def gradient_matrix(mesh, GNe: np.ndarray, quadrature_order: int = 1): - """Computes the gradient operator acting on the FEM mesh's function space and returning gradients at all element quadrature points. - - Args: - mesh: The FEM mesh - GNe (np.ndarray): The shape function gradients evaluated at points of the specified - quadrature. - quadrature_order (int, optional): Polynomial quadrature to use. Defaults to 1. - - Returns: - The gradient operator - """ - mesh_name = _mesh_type_name(mesh) - class_name = f"GradientMatrix_QuadratureOrder_{quadrature_order}_{mesh_name}" - class_ = getattr(_fem, class_name) - return class_(mesh, GNe) - - -def jacobian_determinants(mesh, quadrature_order: int = 1) -> np.ndarray: - """Computes determinants of affine element jacobians, - i.e. the jacobians that map from the reference element to material space. - - Args: - mesh: The FEM mesh - quadrature_order (int, optional): Specifies the polynomial quadrature - to use, such that the jacobian determinants are computed at each quadrature - point. Defaults to 1. - - Returns: - np.ndarray: |#element quadrature points|x|#elements| array of jacobian determinants - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"jacobian_determinants_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, quadrature_order) - - -def inner_product_weights(mesh, quadrature_order: int = 1) -> np.ndarray: - """Computes weights such that element integrals may be computed on reference element quantities, i.e. weights weg = wg(g) * detJe(g,e) - - Args: - mesh: The FEM mesh - quadrature_order (int, optional): Specifies the polynomial quadrature - to use, such that the jacobian determinants and quadrature weights are computed at each quadrature - point. Defaults to 1. - - Returns: - np.ndarray: |#elements * #quadrature points| vector of inner product weights - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"inner_product_weights_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, quadrature_order) - - -def reference_positions(mesh, E: np.ndarray, X: np.ndarray, max_iterations: int = 5, epsilon: float = 1e-10): - """Computes reference positions of domain positions X in corresponding elements E, using Gauss-Newton. - - Args: - mesh: The FEM mesh - E (np.ndarray): 1D index array of element indices - X (np.ndarray): |#dims|x|E.shape[0]| matrix of domain positions in corresponding elements in E - max_iterations (int, optional): Max number of Gauss-Newton iterations. Defaults to 5. - epsilon (float, optional): Residual early out. Defaults to 1e-10. - - Returns: - np.ndarray: |#reference dims|x|E.shape[0]| The reference positions in elements E corresponding to domain positions X in the mesh. - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"reference_positions_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, E, X, max_iterations, epsilon) - - -def integrated_shape_functions(mesh, detJe: np.ndarray, quadrature_order: int = 1) -> np.ndarray: - """Integrates all element shape functions via polynomial quadrature rule - - Args: - mesh: The FEM mesh - detJe (np.ndarray): The jacobian determinants evaluated at points of the specified - quadrature. - quadrature_order (int, optional): Polynomial quadrature to use. Defaults to 1. - - Returns: - np.ndarray: |#element nodes| x |#elements| array of integrated shape functions, i.e. - IN[i,e] = \\int_{\\Omega^e} N_i^e(\\X) \\partial \\Omega^e, for i^{th} shape function - of element e. - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"integrated_shape_functions_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, detJe, quadrature_order) - - -def shape_function_matrix(mesh, quadrature_order: int = 1): - """Computes the matrix N of shape functions values evaluated at element quadrature points - such that N.T @ u evaluates the field u at all element quadrature points - - Args: - mesh: The FEM mesh - quadrature_order (int, optional): Polynomial quadrature rule to use in each element. Defaults to 1. - - Returns: - scipy.sparse.csr_matrix: |#elements * #quadrature points|x|#nodes| sparse matrix - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"shape_function_matrix_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, quadrature_order) - - -def shape_function_gradients(mesh, quadrature_order: int = 1) -> np.ndarray: - """Computes shape function gradients at all quadrature points. - Note that the mesh elements need to be linear transformations on the - reference elements for this method to work properly, even if elements - themselves support higher order shape functions. - - Args: - mesh: The FEM mesh - quadrature_order (int, optional): Polynomial quadrature rule to use. Defaults to 1. - - Returns: - np.ndarray: |#element nodes|x|#dims * #element quad. pts. * #elements| matrix of - shape functions, i.e. if offset=e*dims*nquad, then the block - GNe[i,offset+g*dims:offset+(g+1)*dims] gives the i^{th} shape function gradient of - element e. - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"shape_function_gradients_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, quadrature_order) - - -def shape_functions_at(mesh, Xi: np.ndarray): - """Computes shape functions at reference positions Xi (i.e. positions in element space) for the given mesh. - - Args: - mesh: The FEM mesh - Xi (np.ndarray): Positions in element reference space - - Returns: - np.ndarray: |#element nodes|x|Xi.shape[1]| matrix of nodal shape function values at reference positions Xi - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"shape_functions_at_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, Xi) - - -def shape_function_gradients_at(mesh, E: np.ndarray, Xi: np.ndarray): - """Computes shape function gradients at reference positions Xi (i.e. positions in element space) in elements E for the given mesh. - - Args: - mesh: The FEM mesh - E: Elements in which columns of Xi reside - Xi (np.ndarray): Positions in element reference space - - Returns: - np.ndarray: |#element nodes|x|Xi.shape[1] * mesh.dims| matrix of nodal shape function gradients at reference positions Xi - """ - mesh_name = _mesh_type_name(mesh) - function_name = f"shape_function_gradients_at_{mesh_name}" - function_ = getattr(_fem, function_name) - return function_(mesh, Xi) - - -def laplacian_matrix(mesh, detJe: np.ndarray, GNe: np.ndarray, dims: int = 1, quadrature_order: int = 1): - """Computes the input mesh's (negative semi-definite) symmetric part of the Laplacian matrix. - - Args: - mesh: The FEM mesh - detJe (np.ndarray): The jacobian determinants evaluated at points of the specified - quadrature. - GNe (np.ndarray): The shape function gradients evaluated at points of the specified - quadrature. - dims (int, optional): Image dimensionality of FEM function space. Defaults to 1. - quadrature_order (int, optional): Polynomial quadrature rule to use. Defaults to 1. - - Returns: - The negative semi-definite symmetric part of the FEM mesh's Laplacian matrix - """ - mesh_name = _mesh_type_name(mesh) - class_name = f"SymmetricLaplacianMatrix_QuadratureOrder_{quadrature_order}_{mesh_name}" - class_ = getattr(_fem, class_name) - L = class_(mesh, detJe, GNe) - L.dims = dims - return L - - -def mass_matrix(mesh, detJe: np.ndarray, rho: float = 1., dims: int = 3, quadrature_order: int = 2): - """Computes the input mesh's mass matrix - - Args: - mesh: The FEM mesh - detJe (np.ndarray): Element jacobian determinants at quadrature points - rho (float, optional): Uniform mass density (float) or per-element mass density (np.ndarray). Defaults to 1. - dims (int, optional): dims (int, optional): Image dimensionality of FEM function space. Defaults to 3. - quadrature_order (int, optional): Polynomial quadrature order to use for mass matrix computation. Defaults to 2. - - Returns: - The mass matrix operator - """ - mesh_name = _mesh_type_name(mesh) - class_name = f"MassMatrix_QuadratureOrder_{quadrature_order}_{mesh_name}" - class_ = getattr(_fem, class_name) - M = class_(mesh, detJe, rho) - M.dims = dims - return M - - -def load_vector(mesh, detJe: np.ndarray, fe: np.ndarray, quadrature_order: int = 1): - """Computes a load vector on the input mesh's - - Args: - mesh: The FEM mesh - detJe (np.ndarray): Element jacobian determinants at quadrature points - fe (np.ndarray): |#load dims|x|#elements| per-element load - quadrature_order (int, optional): Polynomial quadrature order to use for load vector computation. Defaults to 1. - - Returns: - The load vector - """ - mesh_name = _mesh_type_name(mesh) - dims = fe.shape[0] - class_name = f"LoadVector_QuadratureOrder_{quadrature_order}_{mesh_name}" - class_ = getattr(_fem, class_name) - return class_(mesh, detJe, fe, dims) - - -class HyperElasticEnergy(Enum): - StVk = 0 - StableNeoHookean = 1 - - -def hyper_elastic_potential( - mesh, - detJe: np.ndarray, - GNe: np.ndarray, - Y: np.ndarray, - nu: np.ndarray, - psi: HyperElasticEnergy = HyperElasticEnergy.StableNeoHookean, - quadrature_order: int = 1): - """Constructs the input mesh's hyper elastic potential, which can be used to evaluate - the potential, its gradient given some state vector (i.e. the DOFs). - - Args: - mesh: The FEM mesh - detJe (np.ndarray): Element jacobian determinants at quadrature points - GNe (np.ndarray): The shape function gradients evaluated at points of the specified - quadrature. - Y (np.ndarray): Array of per-element Young's modulus - nu (np.ndarray): Array of per-element Poisson's ratio - psi (HyperElasticEnergy, optional): The type of hyper elastic energy. Defaults to HyperElasticEnergy.StableNeohookean. - quadrature_order (int, optional): Polynomial quadrature order to use. Defaults to 1. - dims (int, optional): The problem's dimensionality. Defaults to 3. - - Returns: - A hyper elastic potential instance - """ - mesh_name = _mesh_type_name(mesh) - class_name = f"HyperElasticPotential_{psi.name}_QuadratureOrder_{quadrature_order}_Dims_{mesh.dims}_{mesh_name}" - class_ = getattr(_fem, class_name) - return class_(mesh, detJe, GNe, Y, nu) diff --git a/source/pbat/fem/Jacobian.h b/source/pbat/fem/Jacobian.h index a4b35e0..81fd5d7 100644 --- a/source/pbat/fem/Jacobian.h +++ b/source/pbat/fem/Jacobian.h @@ -39,6 +39,14 @@ template return detJ; } +/** + * @brief + * @tparam QuadratureOrder + * @tparam TMesh + * @param mesh + * @return |#quad.pts.|x|#elements| matrix of element jacobian determinants at element quadrature + * points + */ template MatrixX DeterminantOfJacobian(TMesh const& mesh) { @@ -105,7 +113,8 @@ MatrixX InnerProductWeights(TMesh const& mesh) * @tparam TDerivedDetJe * @param mesh * @param detJe - * @return MatrixX + * @return |#quad.pts.|x|#elements| matrix of quadrature weights multiplied by jacobian determinants + * at element quadrature points */ template MatrixX InnerProductWeights(TMesh const& mesh, Eigen::MatrixBase const& detJe) @@ -185,6 +194,19 @@ Vector ReferencePosition( return Xik; } +/** + * @brief + * @tparam TDerivedE + * @tparam TDerivedX + * @tparam TMesh + * @param mesh + * @param E + * @param X + * @param maxIterations + * @param eps + * @return |#element dims| x |X.cols()| matrix of reference positions associated with domain points + * X in corresponding elements E + */ template MatrixX ReferencePositions( TMesh const& mesh, diff --git a/source/pbat/fem/ShapeFunctions.h b/source/pbat/fem/ShapeFunctions.h index 368e60f..72760b6 100644 --- a/source/pbat/fem/ShapeFunctions.h +++ b/source/pbat/fem/ShapeFunctions.h @@ -42,7 +42,11 @@ ShapeFunctions() } /** - * + * @brief + * @tparam QuadratureOrder + * @tparam TMesh + * @param mesh + * @return |#elements * #quad.pts.| x |#nodes| shape function matrix */ template CSRMatrix ShapeFunctionMatrix(TMesh const& mesh) @@ -79,7 +83,7 @@ CSRMatrix ShapeFunctionMatrix(TMesh const& mesh) * @tparam TDerivedXi * @tparam TElement * @param Xi - * @return + * @return |#element nodes| x |Xi.cols()| matrix of nodal shape functions at reference points Xi */ template MatrixX ShapeFunctionsAt(Eigen::DenseBase const& Xi) @@ -275,7 +279,8 @@ MatrixX ShapeFunctionGradients(TMesh const& mesh) * @param mesh * @param E * @param Xi - * @return + * @return |#element nodes| x |E.size() * mesh.dims| nodal shape function gradients at reference + * points Xi */ template MatrixX ShapeFunctionGradientsAt(