diff --git a/doc/dev-docs/source/ConstitutiveLaws.rst b/doc/dev-docs/source/ConstitutiveLaws.rst new file mode 100644 index 0000000..309b52b --- /dev/null +++ b/doc/dev-docs/source/ConstitutiveLaws.rst @@ -0,0 +1,15 @@ + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + +.. _constitutive_laws: + + +Constitutive Laws +~~~~~~~~~~~~~~~~~ + +.. toctree:: + :maxdepth: 2 + + MaterialLinearElasticGeneric diff --git a/doc/dev-docs/source/MaterialLinearElasticGeneric.rst b/doc/dev-docs/source/MaterialLinearElasticGeneric.rst new file mode 100644 index 0000000..d41db25 --- /dev/null +++ b/doc/dev-docs/source/MaterialLinearElasticGeneric.rst @@ -0,0 +1,31 @@ +Generic Linear Elastic Material +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The generic linear elastic material is implemented in the class :cpp:class:`muSpectre::MaterialLinearElasticGeneric`, and is defined solely by the elastic stiffness tensor :math:`\mathbb{C}`, which has to be specified in `Voigt notation `_. The constitutive relation between the Cauchy stress :math:`\boldsymbol{\sigma}` and the small strain tensor :math:`\boldsymbol{\varepsilon}` is given by + +.. math:: + :nowrap: + + \begin{align} + \boldsymbol{\sigma} &= \mathbb{C}:\boldsymbol{\varepsilon}\\ + \sigma_{ij} &= C_{ijkl}\,\varepsilon_{kl} + \end{align} + +This implementation is convenient, as it covers all possible linear elastic behaviours, but it is by far not as efficient as :cpp:class:`muSpectre::MaterialLinearElastic1` for isotropic linear elasticity. + +This law can be used in both small strain and finite strain calculations. + +The following snippet shows how to use this law in python to implement isotropic linear elasticity: + +.. code-block:: python + + C = np.array([[2 * mu + lam, lam, lam, 0, 0, 0], + [ lam, 2 * mu + lam, lam, 0, 0, 0], + [ lam, lam, 2 * mu + lam, 0, 0, 0], + [ 0, 0, 0, mu, 0, 0], + [ 0, 0, 0, 0, mu, 0], + [ 0, 0, 0, 0, 0, mu]]) + + mat = muSpectre.material.MaterialLinearElasticGeneric_3d.make( + cell, "material", C) + diff --git a/doc/dev-docs/source/NewMaterial.rst b/doc/dev-docs/source/NewMaterial.rst index c093398..600a281 100644 --- a/doc/dev-docs/source/NewMaterial.rst +++ b/doc/dev-docs/source/NewMaterial.rst @@ -1,191 +1,206 @@ Writing a New Constitutive Law ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The abstraction for a constitutive law in *µ*\Spectre** is the ``Material``, and new such materials must inherit from the class :cpp:class:`muSpectre::MaterialBase`. Most often, however, it will be most convenient to inherit from the derived class :cpp:class:`muSpectre::MaterialMuSpectre`, as it implements a lot of the machinery that is commonly used by constitutive laws. This section describes how to implement a new constitutive law with internal variables (sometimes also called state variables). The example material implemented here is :cpp:class:`MaterialTutorial`, an objective linear elastic law with a distribution of eigenstrains as internal variables. The constitutive law is defined by the relationship between material (or second Piola-Kirchhoff) stress :math:`\mathbf{S}` and Green-Lagrange strain :math:`\mathbf{E}` +The abstraction for a constitutive law in *µ*\Spectre** is the ``Material``, and +new such materials must inherit from the class +:cpp:class:`muSpectre::MaterialBase`. Most often, however, it will be most +convenient to inherit from the derived class +:cpp:class:`muSpectre::MaterialMuSpectre`, as it implements a lot of the +machinery that is commonly used by constitutive laws. This section describes how +to implement a new constitutive law with internal variables (sometimes also +called state variables). The example material implemented here is +:cpp:class:`MaterialTutorial`, an objective linear elastic law with a +distribution of eigenstrains as internal variables. The constitutive law is +defined by the relationship between material (or second Piola-Kirchhoff) stress +:math:`\mathbf{S}` and Green-Lagrange strain :math:`\mathbf{E}` .. math:: :nowrap: \begin{align} \mathbf{S} &= \mathbb C:\left(\mathbf{E}-\mathbf{e}\right), \\ S_{ij} &= C_{ijkl}\left(E_{kl}-e_{kl}\right), \\ \end{align} -where :math:`\mathbb C` is the elastic stiffness tensor and :math:`\mathbf e` is the local eigenstrain. Note that the implementation sketched out here is the most convenient to quickly get started with using *µ*\Spectre**, but not the most efficient one. For a most efficient implementation, refer to the implementation of :cpp:class:`muSpectre::MaterialLinearElastic2`. +where :math:`\mathbb C` is the elastic stiffness tensor and :math:`\mathbf e` is +the local eigenstrain. Note that the implementation sketched out here is the +most convenient to quickly get started with using *µ*\Spectre**, but not the +most efficient one. For a most efficient implementation, refer to the +implementation of :cpp:class:`muSpectre::MaterialLinearElastic2`. The :cpp:class:`muSpectre::MaterialMuSpectre` class *************************************************** The class :cpp:class:`muSpectre::MaterialMuSpectre` is defined in ``material_muSpectre_base.hh`` and takes three template parameters; #. ``class Material`` is a `CRTP `_ and names the material inheriting from it. The reason for this construction is that we want to avoid virtual method calls from :cpp:class:`muSpectre::MaterialMuSpectre` to its derived classes. Rather, we want :cpp:class:`muSpectre::MaterialMuSpectre` to be able to call methods of the inheriting class directly without runtime overhead. #. ``Dim_t DimS`` defines the number of spatial dimensions of the problem, i.e., whether we are dealing with a two- or three-dimensional grid of pixels/voxels. #. ``Dim_t DimM`` defines the number of dimensions of our material description. This value will typically be the same as ``DimS``, but in cases like generalised plane strain, we can for instance have a three three-dimensional material response in a two-dimensional pixel grid. The main job of :cpp:class:`muSpectre::MaterialMuSpectre` is to #. loop over all pixels to which this material has been assigned, transform the global gradient :math:`\mathbf{F}` (or small strain tensor :math:`\boldsymbol\varepsilon`) into the new material's required strain measure (e.g., the Green-Lagrange strain tensor :math:`\mathbf{E}`), #. for each pixel evaluate the constitutive law by calling its ``evaluate_stress`` (computes the stress response) or ``evaluate_stress_tangent`` (computes both stress and consistent tangent) method with the local strain and internal variables, and finally #. transform the stress (and possibly tangent) response from the material's stress measure into first Piola-Kirchhoff stress :math:`\mathbf{P}` (or Cauchy stress :math:`\boldsymbol\sigma` in small strain). In order to use these facilities, the new material needs to inherit from :cpp:class:`muSpectre::MaterialMuSpectre` (where we calculation of the response) and specialise the type :cpp:class:`muSpectre::MaterialMuSpectre_traits` (where we tell :cpp:class:`muSpectre::MaterialMuSpectre` how to use the new material). These two steps are described here for our example material. Specialising the :cpp:class:`muSpectre::MaterialMuSpectre_traits` structure *************************************************************************** This structure is templated by the new material (in this case :cpp:class:`MaterialTutorial`) and needs to specify #. the types used to communicate per-pixel strains, stresses and stiffness tensors to the material (i.e., whether you want to get maps to `Eigen` matrices or raw pointers, or ...). Here we will use the convenient :cpp:type:`muSpectre::MatrixFieldMap` for strains and stresses, and :cpp:type:`muSpectre::T4MatrixFieldMap` for the stiffness. Look through the classes deriving from :cpp:type:`muSpectre::FieldMap` for all available options. #. the strain measure that is expected (e.g., gradient, Green-Lagrange strain, left Cauchy-Green strain, etc.). Here we will use Green-Lagrange strain. The options are defined by the enum :cpp:enum:`muSpectre::StrainMeasure`. #. the stress measure that is computed by the law (e.g., Cauchy, first Piola-Kirchhoff, etc,). Here, it will be first Piola-Kirchhoff stress. The available options are defined by the enum :cpp:enum:`muSpectre::StressMeasure`. Our traits look like this (assuming we are in the namespace ``muSpectre``:: template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename GlobalFieldCollection; //! expected map type for strain fields using StrainMap_t = MatrixFieldMap; //! expected map type for stress fields using StressMap_t = MatrixFieldMap; //! expected map type for tangent stiffness fields using TangentMap_t = T4MatrixFieldMap; //! declare what type of strain measure your law takes as input constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; //! declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::PK2}; //! local field_collections used for internals using LFieldColl_t = LocalFieldCollection; //! local strain type using LStrainMap_t = MatrixFieldMap; //! elasticity with eigenstrain using InternalVariables = std::tuple; }; Implementing the new material ***************************** The new law needs to implement the methods ``add_pixel``, ``get_internals``, ``evaluate_stress``, and ``evaluate_stress_tangent``. Below is a commented example header:: template class MaterialTutorial: public MaterialMuSpectre, DimS, DimM> { public: //! traits of this material using traits = MaterialMuSpectre_traits; //! Type of container used for storing eigenstrain using InternalVariables = typename traits::InternalVariables; //! Construct by name, Young's modulus and Poisson's ratio MaterialTutorial(std::string name, Real young, Real poisson); /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor) */ template inline decltype(auto) evaluate_stress(s_t && E, eigen_s_t && E_eig); /** * evaluates both second Piola-Kirchhoff stress and stiffness given * the Green-Lagrange strain (or Cauchy stress and stiffness if * called with a small strain tensor) */ template inline decltype(auto) evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig); /** * return the internals tuple (needed by `muSpectre::MaterialMuSpectre`) */ InternalVariables & get_internals() { return this->internal_variables;}; /** * overload add_pixel to write into eigenstrain */ void add_pixel(const Ccoord_t & pixel, const Eigen::Matrix & E_eig); protected: //! stiffness tensor T4Mat C; //! storage for eigenstrain using Field_t = TensorField, Real, secondOrder, DimM>; Field_t & eigen_field; //!< field of eigenstrains //! tuple for iterable eigen_field InternalVariables internal_variables; private: }; A possible implementation for the constructor would be:: template MaterialTutorial::MaterialTutorial(std::string name, Real young, Real poisson) :MaterialMuSpectre(name) { // Lamé parameters Real lambda{young*poisson/((1+poisson)*(1-2*poisson))}; Real mu{young/(2*(1+poisson))}; // Kronecker delta Eigen::Matrix del{Eigen::Matrix::Identity()}; // fill the stiffness tensor this->C.setZero(); for (Dim_t i = 0; i < DimM; ++i) { for (Dim_t j = 0; j < DimM; ++j) { for (Dim_t k = 0; k < DimM; ++k) { for (Dim_t l = 0; l < DimM; ++l) { get(this->C, i, j, k, l) += (lambda * del(i,j)*del(k,l) + mu * (del(i,k)*del(j,l) + del(i,l)*del(j,k))); } } } } } as an exercise, you could check how :cpp:class:`muSpectre::MaterialLinearElastic1` uses *µ*\Spectre**'s materials toolbox (in namespace ``MatTB``) to compute :math:`\mathbb C` in a much more convenient fashion. The evaluation of the stress could be (here, we make use of the ``Matrices`` namespace that defines common tensor algebra operations):: template template decltype(auto) MaterialTutorial:: evaluate_stress(s_t && E, eigen_s_t && E_eig) { return Matrices::tens_mult(this->C, E-E_eig); } The remaining two methods are straight-forward:: template template decltype(auto) MaterialTutorial:: evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig) { return return std::make_tuple (evaluate_stress(E, E_eig), this->C); } template InternalVariables & MaterialTutorial::get_internals() { return this->internal_variables; } Note that the methods ``evaluate_stress`` and ``evaluate_stress_tangent`` need to be in the header, as both their input parameter types and output type depend on the compile-time context. diff --git a/doc/dev-docs/source/conf.py b/doc/dev-docs/source/conf.py index 12bbc54..e4ece92 100644 --- a/doc/dev-docs/source/conf.py +++ b/doc/dev-docs/source/conf.py @@ -1,224 +1,226 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- # # µSpectre documentation build configuration file, created by # sphinx-quickstart on Thu Feb 1 12:01:24 2018. # # This file is execfile()d with the current directory set to its # containing dir. # # Note that not all possible configuration values are present in this # autogenerated file. # # All configuration values have a default; values that are commented out # serve to show the default. # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # import os import subprocess # import sys # sys.path.insert(0, os.path.abspath('.')) # -- General configuration ------------------------------------------------ # If your documentation needs a minimal Sphinx version, state it here. # # needs_sphinx = '1.0' # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = ['sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.todo', 'sphinx.ext.coverage', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', 'breathe'] read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True' if os.environ.get('READTHEDOCS', None) is not None: print("${READTHEDOCS} = " + os.environ.get('READTHEDOCS', None)) if read_the_docs_build: muspectre_path = "." else: muspectre_path = "@CMAKE_CURRENT_BINARY_DIR@" os.makedirs("@CMAKE_CURRENT_BINARY_DIR@/_static", exist_ok=True) print("muspectre_path = '{}'".format(muspectre_path)) subprocess.call('ls; pwd', shell=True) subprocess.call("cd {} && doxygen".format(muspectre_path), shell=True) breathe_projects = {"µSpectre": os.path.join(muspectre_path, "doxygenxml")} breathe_default_project = "µSpectre" # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: # # source_suffix = ['.rst', '.md'] source_suffix = '.rst' # The master toctree document. master_doc = 'index' # General information about the project. project = 'µSpectre' copyright = '2018, Till Junge' author = 'Till Junge' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the # built documents. # # The short X.Y version. version = 'v0.1' # The full version, including alpha/beta/rc tags. release = 'v0.1 α' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. language = None # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This patterns also effect to html_static_path and html_extra_path exclude_patterns = ["CmakeLists.txt"] # The name of the Pygments (syntax highlighting) style to use. pygments_style = 'sphinx' # If true, `todo` and `todoList` produce output, else they produce nothing. todo_include_todos = True # -- Options for HTML output ---------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # on_rtd = os.environ.get('READTHEDOCS') == 'True' if on_rtd: html_theme = 'default' else: html_theme = 'classic' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. # # html_theme_options = {} # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". html_static_path = ['_static'] # Custom sidebar templates, must be a dictionary that maps document names # to template names. # # This is required for the alabaster theme # refs: http://alabaster.readthedocs.io/en/latest/installation.html#sidebars html_sidebars = { '**': [ 'relations.html', # needs 'show_related': True theme option to display 'searchbox.html', ] } # -- Options for HTMLHelp output ------------------------------------------ # Output file base name for HTML help builder. htmlhelp_basename = 'Spectredoc' # -- Options for LaTeX output --------------------------------------------- latex_elements = { # The paper size ('letterpaper' or 'a4paper'). # # 'papersize': 'letterpaper', # The font size ('10pt', '11pt' or '12pt'). # # 'pointsize': '10pt', # Additional stuff for the LaTeX preamble. # - # 'preamble': '', + 'preamble': r''' +\usepackage{amsmath} +''', # Latex figure (float) alignment # # 'figure_align': 'htbp', } # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, # author, documentclass [howto, manual, or own class]). latex_documents = [ (master_doc, 'Spectre.tex', 'µSpectre Documentation', 'Till Junge', 'manual'), ] # -- Options for manual page output --------------------------------------- # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). man_pages = [ (master_doc, 'spectre', 'µSpectre Documentation', [author], 1) ] # -- Options for Texinfo output ------------------------------------------- # Grouping the document tree into Texinfo files. List of tuples # (source start file, target name, title, author, # dir menu entry, description, category) texinfo_documents = [ (master_doc, 'Spectre', 'µSpectre Documentation', author, 'Spectre', 'One line description of project.', 'Miscellaneous'), ] # -- Options for Epub output ---------------------------------------------- # Bibliographic Dublin Core info. epub_title = project epub_author = author epub_publisher = author epub_copyright = copyright # The unique identifier of the text. This can be a ISBN number # or the project homepage. # # epub_identifier = '' # A unique identification for the text. # # epub_uid = '' # A list of files that should not be packed into the epub file. epub_exclude_files = ['search.html'] # Example configuration for intersphinx: refer to the Python standard library. intersphinx_mapping = {'https://docs.python.org/': None} primary_domain = 'cpp' highlight_language = 'cpp' diff --git a/doc/dev-docs/source/index.rst b/doc/dev-docs/source/index.rst index b1dd156..aea28cf 100644 --- a/doc/dev-docs/source/index.rst +++ b/doc/dev-docs/source/index.rst @@ -1,41 +1,42 @@ .. image:: ../logo_flat.png *µ*\Spectre, FFT-based Homogenisation without Linear Reference Medium ===================================================================== .. toctree:: :maxdepth: 2 :caption: Contents: Summary Tutorials CodingConvention + ConstitutiveLaws Reference License *µ*\Spectre is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 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 Lesser General Public License along with µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Indices and tables ================== * :ref:`genindex` * :ref:`modindex` * :ref:`search` diff --git a/language_bindings/python/CMakeLists.txt b/language_bindings/python/CMakeLists.txt index 3b2e483..f3d7472 100644 --- a/language_bindings/python/CMakeLists.txt +++ b/language_bindings/python/CMakeLists.txt @@ -1,87 +1,88 @@ #============================================================================== # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief configuration for python binding using pybind11 # # @section LICENSE # # Copyright © 2018 Till Junge # # µSpectre is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public License as # published by the Free Software Foundation, either version 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 Lesser General Public License # along with µSpectre; see the file COPYING. If not, write to the # Free Software Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. # # Additional permission under GNU GPL version 3 section 7 # # If you modify this Program, or any covered work, by linking or combining it # with proprietary FFT implementations or numerical libraries, containing parts # covered by the terms of those libraries' licenses, the licensors of this # Program grant you additional permission to convey the resulting work. # ============================================================================= # FIXME! The user should have a choice to configure this path. execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-m" "site" "--user-site" RESULT_VARIABLE _PYTHON_SUCCESS OUTPUT_VARIABLE PYTHON_USER_SITE ERROR_VARIABLE _PYTHON_ERROR_VALUE) if(NOT _PYTHON_SUCCESS MATCHES 0) message(FATAL_ERROR "Python config failure:\n${_PYTHON_ERROR_VALUE}") endif() string(REGEX REPLACE "\n" "" PYTHON_USER_SITE ${PYTHON_USER_SITE}) set (PY_BINDING_SRCS ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_module.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_common.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_cell.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material.cc + ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic_generic.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_solvers.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_fftengine.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_projections.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_field_collection.cc ) if (${USE_FFTWMPI}) add_definitions(-DWITH_FFTWMPI) endif(${USE_FFTWMPI}) if (${USE_PFFT}) add_definitions(-DWITH_PFFT) endif(${USE_PFFT}) find_package(PythonLibsNew ${MUSPECTRE_PYTHON_MAJOR_VERSION} MODULE REQUIRED) # On OS X, add -Wno-deprecated-declarations IF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") add_compile_options(-Wno-deprecated-declarations) ENDIF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") pybind11_add_module(pyMuSpectreLib ${PY_BINDING_SRCS}) target_link_libraries(pyMuSpectreLib PRIVATE muSpectre) # Want to rename the output, so that the python module is called muSpectre set_target_properties(pyMuSpectreLib PROPERTIES OUTPUT_NAME _muSpectre) target_include_directories(pyMuSpectreLib PUBLIC ${PYTHON_INCLUDE_DIRS}) add_custom_target(pyMuSpectre ALL SOURCES muSpectre/__init__.py muSpectre/fft.py) add_custom_command(TARGET pyMuSpectre POST_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/language_bindings/python/muSpectre $/muSpectre) install(TARGETS pyMuSpectreLib LIBRARY DESTINATION ${PYTHON_USER_SITE}) install(FILES muSpectre/__init__.py muSpectre/fft.py DESTINATION ${PYTHON_USER_SITE}/muSpectre) diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc index c8a1c49..4001d6e 100644 --- a/language_bindings/python/bind_py_material.cc +++ b/language_bindings/python/bind_py_material.cc @@ -1,183 +1,188 @@ /** * @file bind_py_material.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for µSpectre's materials * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/common.hh" #include "materials/material_linear_elastic1.hh" #include "materials/material_linear_elastic2.hh" #include "materials/material_linear_elastic3.hh" #include "materials/material_linear_elastic4.hh" #include "cell/cell_base.hh" #include "common/field_collection.hh" #include #include #include #include #include using namespace muSpectre; // NOLINT // TODO(junge): figure this out namespace py = pybind11; using namespace pybind11::literals; // NOLINT: recommended usage +/* ---------------------------------------------------------------------- */ +template +void add_material_linear_elastic_generic_helper(py::module & mod); + /** * python binding for the optionally objective form of Hooke's law */ template void add_material_linear_elastic1_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic1_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElastic1; using Sys_t = CellBase; py::class_>(mod, name.c_str()) .def_static("make", [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { return Mat_t::make(sys, n, e, p); }, "cell"_a, "name"_a, "Young"_a, "Poisson"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [](Mat_t & mat, Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a) .def("size", &Mat_t::size); } template void add_material_linear_elastic2_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic2_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElastic2; using Sys_t = CellBase; py::class_>(mod, name.c_str()) .def_static("make", [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { return Mat_t::make(sys, n, e, p); }, "cell"_a, "name"_a, "Young"_a, "Poisson"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [](Mat_t & mat, Ccoord_t pix, py::EigenDRef & eig) { Eigen::Matrix eig_strain{eig}; mat.add_pixel(pix, eig_strain); }, "pixel"_a, "eigenstrain"_a) .def("size", &Mat_t::size); } template void add_material_linear_elastic3_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic3_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElastic3; using Sys_t = CellBase; py::class_>(mod, name.c_str()) .def(py::init(), "name"_a) .def_static("make", [](Sys_t & sys, std::string n) -> Mat_t & { return Mat_t::make(sys, n); }, "cell"_a, "name"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [](Mat_t & mat, Ccoord_t pix, Real Young, Real Poisson) { mat.add_pixel(pix, Young, Poisson); }, "pixel"_a, "Young"_a, "Poisson"_a) .def("size", &Mat_t::size); } template void add_material_linear_elastic4_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic4_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElastic4; using Sys_t = CellBase; py::class_>(mod, name.c_str()) .def(py::init(), "name"_a) .def_static("make", [](Sys_t & sys, std::string n) -> Mat_t & { return Mat_t::make(sys, n); }, "cell"_a, "name"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [](Mat_t & mat, Ccoord_t pix, Real Young, Real Poisson) { mat.add_pixel(pix, Young, Poisson); }, "pixel"_a, "Young"_a, "Poisson"_a) .def("size", &Mat_t::size); } template void add_material_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "Material_" << dim << 'd'; const std::string name{name_stream.str()}; using Mat_t = MaterialBase; using FC_t = LocalFieldCollection; using FCBase_t = FieldCollectionBase; py::class_(mod, name.c_str()) .def_property_readonly("collection", [](Mat_t & material) -> FCBase_t & { return material.get_collection(); }, "returns the field collection containing internal " "fields of this material", py::return_value_policy::reference_internal); add_material_linear_elastic1_helper(mod); add_material_linear_elastic2_helper(mod); add_material_linear_elastic3_helper(mod); add_material_linear_elastic4_helper(mod); + add_material_linear_elastic_generic_helper(mod); } void add_material(py::module & mod) { auto material{mod.def_submodule("material")}; material.doc() = "bindings for constitutive laws"; add_material_helper(material); add_material_helper(material); } diff --git a/language_bindings/python/bind_py_material_linear_elastic_generic.cc b/language_bindings/python/bind_py_material_linear_elastic_generic.cc new file mode 100644 index 0000000..9c0cbe3 --- /dev/null +++ b/language_bindings/python/bind_py_material_linear_elastic_generic.cc @@ -0,0 +1,86 @@ +/** + * @file bind_py_material_linear_elastic_generic.cc + * + * @author Till Junge + * + * @date 20 Dec 2018 + * + * @brief bindings for the generic linear elastic law defined by its stiffness + * tensor + * + * Copyright © 2018 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 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 Lesser General Public License + * along with µSpectre; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + * + * Additional permission under GNU GPL version 3 section 7 + * + * If you modify this Program, or any covered work, by linking or combining it + * with proprietary FFT implementations or numerical libraries, containing parts + * covered by the terms of those libraries' licenses, the licensors of this + * Program grant you additional permission to convey the resulting work. + */ + +#include "materials/material_linear_elastic_generic.hh" +#include "cell/cell_base.hh" + +#include +#include +#include + + +#include +#include + +using namespace muSpectre; // NOLINT // TODO(junge): figure this out +namespace py = pybind11; +using namespace pybind11::literals; // NOLINT: recommended usage + +/** + * python binding for the generic linear elastic material + */ +template +void add_material_linear_elastic_generic_helper(py::module & mod) { +std::stringstream name_stream{}; +name_stream << "MaterialLinearElasticGeneric_" << Dim << "d"; +const auto name{name_stream.str()}; + +using Mat_t = MaterialLinearElasticGeneric; +using Cell_t = CellBase; + +py::class_>(mod, name.c_str()) + .def_static("make", + [](Cell_t & cell, std::string name, + const py::EigenDRef< + Eigen::Matrix> & + elastic_tensor) -> Mat_t & { + return Mat_t::make(cell, name, elastic_tensor); + }, + "cell"_a, "name"_a, "elastic_tensor"_a, + py::return_value_policy::reference, py::keep_alive<1, 0>(), + "Factory function returning a MaterialLinearElastic instance. " + "The elastic tensor has to be specified in Voigt notation.") + .def("add_pixel", + [](Mat_t & mat, Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a, + "Register a new pixel to this material. subsequent evaluations of the " + "stress and tangent in the cell will use this constitutive law for " + "this " + "particular pixel") + .def("size", &Mat_t::size); +} + + +template void add_material_linear_elastic_generic_helper(py::module &); +template void add_material_linear_elastic_generic_helper(py::module &); diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py index 636c99b..d87a950 100755 --- a/tests/python_binding_tests.py +++ b/tests/python_binding_tests.py @@ -1,201 +1,202 @@ #!/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 Lesser 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 Lesser General Public License along with µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Additional permission under GNU GPL version 3 section 7 If you modify this Program, or any covered work, by linking or combining it with proprietary FFT implementations or numerical libraries, containing parts covered by the terms of those libraries' licenses, the licensors of this Program grant you additional permission to convey the resulting work. """ 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_material_linear_elastic_generic_test import MaterialLinearElasticGeneric_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) self.mat3 = µ.material.MaterialLinearElastic2_2d.make( self.cell2, "eigen2", 120e9, .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_real_array("Eigenstrain") error = np.linalg.norm(loc_eigenstrain-glo_eigenstrain) self.assertEqual(error, 0) def test_globalisation_constant(self): for i, pixel in enumerate(self.cell2): if i%2 == 0: self.mat2.add_pixel(pixel, np.ones((2,2))) else: self.mat3.add_pixel(pixel, np.ones((2,2))) glo_eigenstrain = self.cell2.get_globalised_internal_real_array("Eigenstrain") error = np.linalg.norm(glo_eigenstrain-1) 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() diff --git a/tests/python_material_linear_elastic_generic_test.py b/tests/python_material_linear_elastic_generic_test.py new file mode 100644 index 0000000..7b8754c --- /dev/null +++ b/tests/python_material_linear_elastic_generic_test.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python3 +# -*- coding:utf-8 -*- +""" +@file python_material_linear_elastic_generic.py + +@author Till Junge + +@date 20 Dec 2018 + +@brief tests the python bindings of the generic linear elastic material + +Copyright © 2018 Till Junge + +µSpectre is free software; you can redistribute it and/or +modify it under the terms of the GNU General Lesser 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 Lesser General Public License +along with µSpectre; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 59 Temple Place - Suite 330, +Boston, MA 02111-1307, USA. + +Additional permission under GNU GPL version 3 section 7 + +If you modify this Program, or any covered work, by linking or combining it +with proprietary FFT implementations or numerical libraries, containing parts +covered by the terms of those libraries' licenses, the licensors of this +Program grant you additional permission to convey the resulting work. +""" + + +import unittest +import numpy as np + +from python_test_imports import µ + +class MaterialLinearElasticGeneric_Check(unittest.TestCase): + def setUp(self): + self.resolution = [5,7,5] + self.dim = len(self.resolution) + self.lengths = [5.2, 8.3, 2.7] + self.formulation = µ.Formulation.small_strain + self.cell1 = µ.Cell(self.resolution, + self.lengths, + self.formulation) + self.Young = 210e9 + self.Poisson = .33 + self.mat1 = µ.material.MaterialLinearElastic1_3d.make( + self.cell1, "material", self.Young, self.Poisson) + self.matO1 = µ.material.MaterialLinearElastic1_3d.make( + self.cell1, "material", 2* self.Young, self.Poisson) + + E, nu = self.Young, self.Poisson + lam, mu = E*nu/((1+nu)*(1-2*nu)), E/(2*(1+nu)) + + C = np.array([[2 * mu + lam, lam, lam, 0, 0, 0], + [ lam, 2 * mu + lam, lam, 0, 0, 0], + [ lam, lam, 2 * mu + lam, 0, 0, 0], + [ 0, 0, 0, mu, 0, 0], + [ 0, 0, 0, 0, mu, 0], + [ 0, 0, 0, 0, 0, mu]]) + + self.cell2 = µ.Cell(self.resolution, + self.lengths, + self.formulation) + self.mat2 = µ.material.MaterialLinearElasticGeneric_3d.make( + self.cell2, "material", C) + self.matO2 = µ.material.MaterialLinearElastic1_3d.make( + self.cell2, "material", 2* self.Young, self.Poisson) + + def test_equivalence(self): + sym = lambda x: .5*(x + x.T) + Del0 = sym((np.random.random((self.dim, self.dim))-.5)/10) + for pixel in self.cell1: + if pixel[0] == 0: + self.matO1.add_pixel(pixel) + self.matO2.add_pixel(pixel) + else: + self.mat1.add_pixel(pixel) + self.mat2.add_pixel(pixel) + + tol = 1e-6 + equil_tol = tol + maxiter = 100 + verbose = 0 + + solver1 = µ.solvers.SolverCG(self.cell1, tol, maxiter, verbose) + solver2 = µ.solvers.SolverCG(self.cell2, tol, maxiter, verbose) + + + r1 = µ.solvers.de_geus(self.cell1, Del0, + solver1, tol, equil_tol, verbose) + + + r2 = µ.solvers.de_geus(self.cell2, Del0, + solver2, tol, equil_tol, verbose) + + error = (np.linalg.norm(r1.stress - r2.stress) / + np.linalg.norm(r1.stress + r2.stress)) + + self.assertLess(error, 1e-13) +