fixed timedependent simulation

Author: teseoch

cmake/PolyfemPythonDownloadExternal.cmake

@@ -27,7 +27,7 @@ endfunction()
 function(polyfem_python_download_polyfem)
     polyfem_python_download_project(polyfem
         GIT_REPOSITORY https://github.com/polyfem/polyfem.git
-        GIT_TAG         41005711400c3e5f0957f70d9f071b1a45d2d123
+        GIT_TAG         c679ec3b52e8be63d4e092fd2f315432045a73f5
     )
 endfunction()
 

setup.py

@@ -24,10 +24,11 @@ def run(self):
             raise RuntimeError(
                 "CMake must be installed to build the following extensions: , ".join(e.name for e in self.extensions))
 
-        # if platform.system() == "Windows":
+        # self.debug = True
+
         cmake_version = LooseVersion(re.search(r'version\s*([\d.]+)', out.decode()).group(1))
         if cmake_version < '3.1.0':
-            raise RuntimeError("CMake >= 3.1.0 is required on Windows")
+            raise RuntimeError("CMake >= 3.1.0 is required")
 
         for ext in self.extensions:
             self.build_extension(ext)

src/binding.cpp

@@ -14,6 +14,7 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/eigen.h>
 #include <pybind11/functional.h>
+#include <pybind11/stl.h>
 
 namespace py = pybind11;
 
@@ -111,6 +112,11 @@ PYBIND11_MODULE(polyfempy, m) {
 		else
 			polyfem::to_geogram_mesh_3d(V, F, M);
 		s.load_mesh(M, [](const polyfem::RowVectorNd&){ return -1; }, true);
+
+		double boundary_id_threshold = s.args["boundary_id_threshold"];
+		if(boundary_id_threshold <= 0)
+			boundary_id_threshold = s.mesh->is_volume() ? 1e-2 : 1e-7;
+		s.mesh->compute_boundary_ids(boundary_id_threshold);
 	},
 	"Loads a mesh from vertices and connectivity",
 	py::arg("vertices"), py::arg("connectivity"))
@@ -162,7 +168,10 @@ PYBIND11_MODULE(polyfempy, m) {
 		s.assemble_rhs();
 		s.assemble_stiffness_mat();
 
+		s.solve_export_to_file = false;
+		s.solution_frames.clear();
 		s.solve_problem();
+		s.solve_export_to_file = true;
 	},
 	"solve the pde")
 
@@ -239,7 +248,80 @@ PYBIND11_MODULE(polyfempy, m) {
 		return py::make_tuple(fun, tfun);
 	},
 	"returns the von mises stresses and stress tensor averaged around a vertex on a densly sampled mesh, use 'vismesh_rel_area' to control density",
-	py::arg("boundary_only") = bool(false));
+	py::arg("boundary_only") = bool(false))
+
+
+
+
+	////////////////////////////////////////////////////////////////////////////////////////////
+	.def("get_sampled_points_frames", 	[](polyfem::State &s) {
+		assert(!s.solution_frames.empty());
+
+		std::vector<Eigen::MatrixXd> pts;
+
+		for(const auto &sol : s.solution_frames){
+			pts.push_back(sol.points);
+		}
+
+
+		return pts;
+	},
+	"returns the points frames for a time dependent problem on a densly sampled mesh, use 'vismesh_rel_area' to control density")
+
+	.def("get_sampled_connectivity_frames", 	[](polyfem::State &s) {
+		assert(!s.solution_frames.empty());
+
+		std::vector<Eigen::MatrixXi> tets;
+
+		for(const auto &sol : s.solution_frames)
+			tets.push_back(sol.connectivity);
+
+
+		return tets;
+	},
+	"returns the connectivity frames for a time dependent problem on a densly sampled mesh, use 'vismesh_rel_area' to control density")
+
+
+	.def("get_sampled_solution_frames", 	[](polyfem::State &s) {
+		assert(!s.solution_frames.empty());
+
+		std::vector<Eigen::MatrixXd> fun;
+
+		for(const auto &sol : s.solution_frames){
+			fun.push_back(sol.solution);
+		}
+
+
+		return fun;
+	},
+	"returns the solution frames for a time dependent problem on a densly sampled mesh, use 'vismesh_rel_area' to control density")
+
+	.def("get_sampled_mises_frames", 	[](polyfem::State &s) {
+		assert(!s.solution_frames.empty());
+
+		std::vector<Eigen::MatrixXd> mises;
+
+		for(const auto &sol : s.solution_frames)
+			mises.push_back(sol.scalar_value);
+
+		return mises;
+	},
+	"returns the von mises stresses frames on a densly sampled mesh, use 'vismesh_rel_area' to control density")
+
+	.def("get_sampled_mises_avg_frames", 	[](polyfem::State &s) {
+		assert(!s.solution_frames.empty());
+
+		std::vector<Eigen::MatrixXd> mises;
+
+		for(const auto &sol : s.solution_frames)
+			mises.push_back(sol.scalar_value_avg);
+
+		return mises;
+	},
+	"returns the von mises stresses per frame averaged around a vertex on a densly sampled mesh, use 'vismesh_rel_area' to control density");
+
+
+
 
 	solver.doc() = "Polyfem solver";
 

tests/gravity.py

@@ -0,0 +1,72 @@
+import unittest
+
+import numpy as np
+import polyfempy as pf
+from utils import plot
+
+
+class Gravity(unittest.TestCase):
+    def test_run(self):
+        n_pts = 3
+
+        extend = np.linspace(0, 1, n_pts)
+        x, y = np.meshgrid(extend, extend, sparse=False, indexing='xy')
+        pts = np.column_stack((x.ravel(), y.ravel()))
+
+
+        # Create connectivity
+        faces = np.ndarray([(n_pts-1)**2, 4],dtype=np.int32)
+
+        index = 0
+        for i in range(n_pts-1):
+            for j in range(n_pts-1):
+                faces[index, :] = np.array([
+                    j + i * n_pts,
+                    j+1 + i * n_pts,
+                    j+1 + (i+1) * n_pts,
+                    j + (i+1) * n_pts
+                ])
+                index = index + 1
+
+        # create settings
+        settings = pf.Settings()
+        settings.discr_order = 1
+
+        settings.set_material_params("E", 2100)
+        settings.set_material_params("nu", 0.3)
+
+        # We use a linear material model
+        settings.tensor_formulation = pf.TensorFormulations.LinearElasticity
+
+        problem = pf.Gravity()
+        settings.set_problem(problem)
+
+        solver = pf.Solver()
+        solver.settings(str(settings))
+
+        # This time we are using pts and faces instead of loading from the disk
+        solver.set_mesh(pts, faces)
+
+        solver.solve()
+
+        #Animation frame, last one
+        frame = -1
+
+        # Get the solution
+        pts = solver.get_sampled_points_frames()
+        tets = solver.get_sampled_connectivity_frames()
+        disp = solver.get_sampled_solution_frames()
+
+
+        # diplace the mesh
+        vertices = pts[frame] + disp[frame]
+
+        # and get the stresses
+        mises = solver.get_sampled_mises_avg_frames()
+
+        # finally plot
+        plot(vertices, tets[frame], mises[frame])
+
+
+if __name__ == '__main__':
+    unittest.main()