diff --git a/language_bindings/python/bind_py_cell.cc b/language_bindings/python/bind_py_cell.cc index a09bc44..8eba510 100644 --- a/language_bindings/python/bind_py_cell.cc +++ b/language_bindings/python/bind_py_cell.cc @@ -1,210 +1,227 @@ /** * @file bind_py_cell.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief Python bindings for the cell factory function * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "common/ccoord_operations.hh" #include "cell/cell_factory.hh" #include "cell/cell_base.hh" #ifdef WITH_FFTWMPI #include "fft/fftwmpi_engine.hh" #endif #ifdef WITH_PFFT #include "fft/pfft_engine.hh" #endif #include #include #include "pybind11/eigen.h" #include #include using namespace muSpectre; namespace py=pybind11; using namespace pybind11::literals; /** * cell factory for specific FFT engine */ #ifdef WITH_MPI template void add_parallel_cell_factory_helper(py::module & mod, const char *name) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def (name, [](Ccoord res, Rcoord lens, Formulation form, size_t comm) { return make_parallel_cell , FFTEngine> (std::move(res), std::move(lens), std::move(form), std::move(Communicator(MPI_Comm(comm)))); }, "resolutions"_a, "lengths"_a=CcoordOps::get_cube(1.), "formulation"_a=Formulation::finite_strain, "communicator"_a=size_t(MPI_COMM_SELF)); } #endif /** * the cell factory is only bound for default template parameters */ template void add_cell_factory_helper(py::module & mod) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def ("CellFactory", [](Ccoord res, Rcoord lens, Formulation form) { return make_cell(std::move(res), std::move(lens), std::move(form)); }, "resolutions"_a, "lengths"_a=CcoordOps::get_cube(1.), "formulation"_a=Formulation::finite_strain); #ifdef WITH_FFTWMPI add_parallel_cell_factory_helper>( mod, "FFTWMPICellFactory"); #endif #ifdef WITH_PFFT add_parallel_cell_factory_helper>( mod, "PFFTCellFactory"); #endif } void add_cell_factory(py::module & mod) { add_cell_factory_helper(mod); add_cell_factory_helper(mod); } /** * CellBase for which the material and spatial dimension are identical */ template void add_cell_base_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "CellBase" << dim << 'd'; const std::string name = name_stream.str(); using sys_t = CellBase; py::class_(mod, name.c_str()) .def("__len__", &sys_t::size) .def("__iter__", [](sys_t & s) { return py::make_iterator(s.begin(), s.end()); }) .def("initialise", &sys_t::initialise, "flags"_a=FFT_PlanFlags::estimate) .def("directional_stiffness", [](sys_t& cell, py::EigenDRef& v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim*dim)) { std::stringstream err{}; err << "need array of shape (" << dim*dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; constexpr bool create_tangent{true}; auto & K = cell.get_tangent(create_tangent); auto & input = cell.get_managed_T2_field(in_name); auto & output = cell.get_managed_T2_field(out_name); input.eigen() = v; cell.directional_stiffness(K, input, output); return output.eigen(); }, "δF"_a) .def("project", [](sys_t& cell, py::EigenDRef& v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim*dim)) { std::stringstream err{}; err << "need array of shape (" << dim*dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string in_name{"temp input for projection"}; auto & input = cell.get_managed_T2_field(in_name); input.eigen() = v; cell.project(input); return input.eigen(); }, "field"_a) .def("get_strain",[](sys_t & s) { return Eigen::ArrayXXd(s.get_strain().eigen()); }) .def("get_stress",[](sys_t & s) { return Eigen::ArrayXXd(s.get_stress().eigen()); }) .def_property_readonly("size", &sys_t::size) .def("evaluate_stress_tangent", [](sys_t& cell, py::EigenDRef& v ) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim*dim)) { std::stringstream err{}; err << "need array of shape (" << dim*dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } auto & strain{cell.get_strain()}; strain.eigen() = v; cell.evaluate_stress_tangent(strain); }, "strain"_a) .def("get_projection", &sys_t::get_projection) .def("get_subdomain_resolutions", &sys_t::get_subdomain_resolutions) .def("get_subdomain_locations", &sys_t::get_subdomain_locations) .def("get_domain_resolutions", &sys_t::get_domain_resolutions) .def("get_domain_lengths", &sys_t::get_domain_resolutions); } void add_cell_base(py::module & mod) { - py::class_(mod, "Cell"); + py::class_(mod, "Cell") + .def("get_globalised_internal_array", + &Cell::get_globalised_internal_array, + "unique_name"_a, + "Convenience function to copy local (internal) fields of " + "materials into a global field. At least one of the materials in " + "the cell needs to contain an internal field named " + "`unique_name`. If multiple materials contain such a field, they " + "all need to be of same scalar type and same number of " + "components. This does not work for split pixel cells or " + "laminate pixel cells, as they can have multiple entries for the " + "same pixel. Pixels for which no field named `unique_name` " + "exists get an array of zeros." + "\n" + "Parameters:\n" + "unique_name: fieldname to fill the global field with. At " + "least one material must have such a field, or an " + "Exception is raised."); add_cell_base_helper (mod); add_cell_base_helper(mod); } void add_cell(py::module & mod) { add_cell_factory(mod); auto cell{mod.def_submodule("cell")}; cell.doc() = "bindings for cells and cell factories"; add_cell_base(cell); } diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py index afe329b..ad971ab 100755 --- a/tests/python_binding_tests.py +++ b/tests/python_binding_tests.py @@ -1,173 +1,182 @@ #!/usr/bin/env python3 """ file python_binding_tests.py @author Till Junge @date 09 Jan 2018 @brief Unit tests for python bindings @section LICENCE Copyright © 2018 Till Junge µSpectre is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3, or (at your option) any later version. µSpectre 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 General Public License for more details. You should have received a copy of the GNU General Public License along with GNU Emacs; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. """ import unittest import numpy as np from python_test_imports import µ from python_fft_tests import FFT_Check from python_projection_tests import * from python_material_linear_elastic3_test import MaterialLinearElastic3_Check from python_material_linear_elastic4_test import MaterialLinearElastic4_Check from python_field_tests import FieldCollection_Check class CellCheck(unittest.TestCase): def test_Construction(self): """ Simple check for cell constructors """ resolution = [5,7] lengths = [5.2, 8.3] formulation = µ.Formulation.small_strain try: sys = µ.Cell(resolution, lengths, formulation) mat = µ.material.MaterialLinearElastic1_2d.make(sys, "material", 210e9, .33) + except Exception as err: print(err) raise err class MaterialLinearElastic1_2dCheck(unittest.TestCase): def setUp(self): self.resolution = [5,7] self.lengths = [5.2, 8.3] self.formulation = µ.Formulation.small_strain self.sys = µ.Cell(self.resolution, self.lengths, self.formulation) self.mat = µ.material.MaterialLinearElastic1_2d.make( self.sys, "material", 210e9, .33) def test_add_material(self): self.mat.add_pixel([2,1]) class SolverCheck(unittest.TestCase): def setUp(self): self.resolution = [3, 3]#[5,7] self.lengths = [3., 3.]#[5.2, 8.3] self.formulation = µ.Formulation.finite_strain self.sys = µ.Cell(self.resolution, self.lengths, self.formulation) self.hard = µ.material.MaterialLinearElastic1_2d.make( self.sys, "hard", 210e9, .33) self.soft = µ.material.MaterialLinearElastic1_2d.make( self.sys, "soft", 70e9, .33) def test_solve(self): for i, pixel in enumerate(self.sys): if i < 3: self.hard.add_pixel(pixel) else: self.soft.add_pixel(pixel) self.sys.initialise() tol = 1e-6 Del0 = np.array([[0, .1], [0, 0]]) maxiter = 100 verbose = 0 solver=µ.solvers.SolverCG(self.sys, tol, maxiter, verbose) r = µ.solvers.de_geus(self.sys, Del0, solver,tol, verbose) #print(r) class EigenStrainCheck(unittest.TestCase): def setUp(self): self.resolution = [3, 3]#[5,7] self.lengths = [3., 3.]#[5.2, 8.3] self.formulation = µ.Formulation.small_strain self.cell1 = µ.Cell(self.resolution, self.lengths, self.formulation) self.cell2 = µ.Cell(self.resolution, self.lengths, self.formulation) self.mat1 = µ.material.MaterialLinearElastic1_2d.make( self.cell1, "simple", 210e9, .33) self.mat2 = µ.material.MaterialLinearElastic2_2d.make( self.cell2, "eigen", 210e9, .33) + def test_globalisation(self): + for pixel in self.cell2: + self.mat2.add_pixel(pixel, np.random.rand(2,2)) + loc_eigenstrain = self.mat2.collection.get_real_field("Eigenstrain").array + glo_eigenstrain = self.cell2.get_globalised_internal_array("Eigenstrain") + error = np.linalg.norm(loc_eigenstrain-glo_eigenstrain) + self.assertEqual(error, 0) + def test_solve(self): verbose_test = False if verbose_test: print("start test_solve") grad = np.array([[1.1, .2], [ .3, 1.5]]) gl_strain = -0.5*(grad.T.dot(grad) - np.eye(2)) gl_strain = -0.5*(grad.T + grad - 2*np.eye(2)) grad = -gl_strain if verbose_test: print("grad =\n{}\ngl_strain =\n{}".format(grad, gl_strain)) for i, pixel in enumerate(self.cell1): self.mat1.add_pixel(pixel) self.mat2.add_pixel(pixel, gl_strain) self.cell1.initialise() self.cell2.initialise() tol = 1e-6 Del0_1 = grad Del0_2 = np.zeros_like(grad) maxiter = 2 verbose = 0 def solve(cell, grad): solver=µ.solvers.SolverCG(cell, tol, maxiter, verbose) r = µ.solvers.newton_cg(cell, grad, solver, tol, tol, verbose) return r results = [solve(cell, del0) for (cell, del0) in zip((self.cell1, self.cell2), (Del0_1, Del0_2))] P1 = results[0].stress P2 = results[1].stress error = np.linalg.norm(P1-P2)/np.linalg.norm(.5*(P1+P2)) if verbose_test: print("cell 1, no eigenstrain") print("P1:\n{}".format(P1[:,0])) print("F1:\n{}".format(results[0].grad[:,0])) print("cell 2, with eigenstrain") print("P2:\n{}".format(P2[:,0])) print("F2:\n{}".format(results[1].grad[:,0])) print("end test_solve") self.assertLess(error, tol) if __name__ == '__main__': unittest.main()