diff --git a/Jenkinsfile b/Jenkinsfile index 969bb63..0cafdd5 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -1,218 +1,218 @@ pipeline { parameters {string(defaultValue: '', description: 'api-token', name: 'API_TOKEN') string(defaultValue: '', description: 'Token for readthedocs', name: 'RTD_TOKEN') string(defaultValue: '', description: 'buildable phid', name: 'TARGET_PHID') string(defaultValue: 'docker_debian_testing', description: 'docker file to use', name: 'DOCKERFILE') string(defaultValue: '', description: 'Commit id', name: 'COMMIT_ID') } agent any environment { OMPI_MCA_plm = 'isolated' OMPI_MCA_btl = 'tcp,self' } options { disableConcurrentBuilds() } stages { stage ('wipe build') { when { anyOf{ changeset glob: "**/*.cmake" changeset glob: "**/CMakeLists.txt" } } steps { sh ' rm -rf build_*' } } stage ('configure') { parallel { stage ('docker_debian_testing') { agent { dockerfile { filename 'docker_debian_testing' dir 'dockerfiles' } } steps { configure('docker_debian_testing') } } stage ('docker_debian_stable') { agent { dockerfile { filename 'docker_debian_stable' dir 'dockerfiles' } } steps { configure('docker_debian_stable') } } } } stage ('build') { parallel { stage ('docker_debian_testing') { agent { dockerfile { filename 'docker_debian_testing' dir 'dockerfiles' } } steps { build('docker_debian_testing') } } stage ('docker_debian_stable') { agent { dockerfile { filename 'docker_debian_stable' dir 'dockerfiles' } } steps { build('docker_debian_stable') } } } } stage ('Warnings gcc') { steps { warnings(consoleParsers: [[parserName: 'GNU Make + GNU C Compiler (gcc)']]) } } stage ('test') { parallel { //////////////////////////////////////////////////// stage ('docker_debian_testing') { agent { dockerfile { filename 'docker_debian_testing' dir 'dockerfiles' } } steps { run_test('docker_debian_testing', 'g++') run_test('docker_debian_testing', 'clang++') } post { always { collect_test_results('docker_debian_testing', 'g++') collect_test_results('docker_debian_testing', 'clang++') } } } //////////////////////////////////////////////////// stage ('docker_debian_stable') { agent { dockerfile { filename 'docker_debian_stable' dir 'dockerfiles' } } steps { run_test('docker_debian_stable', 'g++') run_test('docker_debian_stable', 'clang++') } post { always { collect_test_results('docker_debian_stable', 'g++') collect_test_results('docker_debian_stable', 'clang++') } } } } } } post { always { createartifact() } success { send_fail_pass('pass') trigger_readthedocs() } failure { send_fail_pass('fail') } } } def configure(container_name) { def BUILD_DIR = "build_${container_name}" for (CXX_COMPILER in ["g++", "clang++"]) { sh """ mkdir -p ${BUILD_DIR}_${CXX_COMPILER} cd ${BUILD_DIR}_${CXX_COMPILER} CXX=${CXX_COMPILER} cmake -DCMAKE_BUILD_TYPE:STRING=Release -DRUNNING_IN_CI=ON .. """ } } def build(container_name) { def BUILD_DIR = "build_${container_name}" for (CXX_COMPILER in ["g++", "clang++"]) { sh "make -C ${BUILD_DIR}_${CXX_COMPILER}" } } def run_test(container_name, cxx_compiler) { def BUILD_DIR = "build_${container_name}" sh "cd ${BUILD_DIR}_${cxx_compiler} && ctest || true" } def send_fail_pass(state) { sh """ set +x curl https://c4science.ch/api/harbormaster.sendmessage \ -d api.token=${API_TOKEN} \ -d buildTargetPHID=${TARGET_PHID} \ -d type=${state} """ } def trigger_readthedocs() { - sh """ + sh """ set +x curl -X POST \ -d "token=${RTD_TOKEN}" \ -d "branches=master" \ https://readthedocs.org/api/v2/webhook/muspectre/26537/ curl -X POST \ -d "token=${RTD_TOKEN}" \ https://readthedocs.org/api/v2/webhook/muspectre/26537/ """ } def collect_test_results(container_name, cxx_compiler) { def BUILD_DIR = "build_${container_name}" junit "${BUILD_DIR}_${cxx_compiler}/test_results*.xml" } def createartifact() { sh """ set +x curl https://c4science.ch/api/harbormaster.createartifact \ -d api.token=${API_TOKEN} \ -d buildTargetPHID=${TARGET_PHID} \ -d artifactKey="Jenkins URI" \ -d artifactType=uri \ -d artifactData[uri]=${BUILD_URL} \ -d artifactData[name]="View Jenkins result" \ -d artifactData[ui.external]=1 """ } diff --git a/doc/dev-docs/source/conf.py b/doc/dev-docs/source/conf.py index b441fdf..12bbc54 100644 --- a/doc/dev-docs/source/conf.py +++ b/doc/dev-docs/source/conf.py @@ -1,223 +1,224 @@ #!/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' -print("${READTHEDOCS} = " + os.environ.get('READTHEDOCS', None)) +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': '', # 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/language_bindings/python/CMakeLists.txt b/language_bindings/python/CMakeLists.txt index c6087b5..7677306 100644 --- a/language_bindings/python/CMakeLists.txt +++ b/language_bindings/python/CMakeLists.txt @@ -1,75 +1,76 @@ #============================================================================== # 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 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. # # 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_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) 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 DESTINATION ${PYTHON_USER_SITE}/muSpectre) diff --git a/language_bindings/python/bind_py_cell.cc b/language_bindings/python/bind_py_cell.cc index 82afac3..ddb649b 100644 --- a/language_bindings/python/bind_py_cell.cc +++ b/language_bindings/python/bind_py_cell.cc @@ -1,210 +1,210 @@ /** * @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_field(in_name); auto & output = cell.get_managed_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_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("size", &sys_t::size) + .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"); 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/language_bindings/python/bind_py_declarations.hh b/language_bindings/python/bind_py_declarations.hh index b10f2fa..8ef2ff4 100644 --- a/language_bindings/python/bind_py_declarations.hh +++ b/language_bindings/python/bind_py_declarations.hh @@ -1,41 +1,42 @@ /** * @file bind_py_common.hh * * @author Till Junge * * @date 12 Jan 2018 * * @brief header for python bindings for the common part of µSpectre * * 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. */ #ifndef BIND_PY_DECLARATIONS_H #define BIND_PY_DECLARATIONS_H #include namespace py = pybind11; void add_common(py::module & mod); void add_cell(py::module & mod); void add_material(py::module & mod); void add_solvers(py::module & mod); void add_fft_engines(py::module & mod); void add_projections(py::module & submodule); +void add_field_collections(py::module & submodule); #endif /* BIND_PY_DECLARATIONS_H */ diff --git a/language_bindings/python/bind_py_field_collection.cc b/language_bindings/python/bind_py_field_collection.cc new file mode 100644 index 0000000..3e27dae --- /dev/null +++ b/language_bindings/python/bind_py_field_collection.cc @@ -0,0 +1,222 @@ +/** + * file bind_py_field_collection.cc + * + * @author Till Junge + * + * @date 05 Jul 2018 + * + * @brief Python bindings for µSpectre field collections + * + * @section LICENSE + * + * 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/field.hh" +#include "common/field_collection.hh" + +#include +#include +#include + +#include +#include + +using namespace muSpectre; +namespace py = pybind11; +using namespace pybind11::literals; + +template +void add_field_collection(py::module & mod) { + std::stringstream name_stream{}; + name_stream << "_" << (FieldCollectionDerived::Global ? "Global" : "Local") + << "FieldCollection_" << Dim << 'd'; + const auto name {name_stream.str()}; + using FC_t = FieldCollectionBase; + py::class_(mod, name.c_str()) + .def("get_real_field", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedField_t & { + return field_collection.template get_typed_field(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def("get_int_field", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedField_t & { + return field_collection.template get_typed_field(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def("get_uint_field", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedField_t & { + return field_collection.template get_typed_field(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def("get_complex_field", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedField_t & { + return field_collection.template get_typed_field(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def("get_real_statefield", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedStateField_t & { + return field_collection.template get_typed_statefield(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def("get_int_statefield", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedStateField_t & { + return field_collection.template get_typed_statefield(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def("get_uint_statefield", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedStateField_t & { + return field_collection.template get_typed_statefield(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def("get_complex_statefield", + [](FC_t & field_collection, std::string unique_name) + -> typename FC_t::template TypedStateField_t & { + return field_collection.template get_typed_statefield(unique_name); + }, + "unique_name"_a, + py::return_value_policy::reference_internal) + .def_property_readonly("field_names", &FC_t::get_field_names, + "returns the names of all fields in this collection") + .def_property_readonly("statefield_names", &FC_t::get_statefield_names, + "returns the names of all state fields in this " + "collection"); +} + +template +void add_field(py::module & mod, std::string dtype_name) { + using Field_t = TypedField; + std::stringstream name_stream{}; + name_stream << (FieldCollection::Global ? "Global" : "Local") + << "Field" << dtype_name << "_" << FieldCollection::spatial_dim(); + std::string name{name_stream.str()}; + using Ref_t = py::EigenDRef>; + py::class_(mod, name.c_str()) + .def_property("array", [](Field_t & field) {return field.eigen();}, + [](Field_t & field, Ref_t mat) {field.eigen() = mat;}, + "array of stored data") + .def_property_readonly("array", + [](const Field_t & field) {return field.eigen();}, + "array of stored data") + .def_property("vector", + [](Field_t& field) {return field.eigenvec();}, + [](Field_t & field, Ref_t mat) {field.eigen() = mat;}, + "flattened array of stored data") + .def_property_readonly("vector", + [](const Field_t& field) {return field.eigenvec();}, + "flattened array of stored data"); +} + +template +void add_field_helper(py::module & mod) { + std::stringstream name_stream{}; + name_stream << (FieldCollection::Global ? "Global" : "Local") + << "Field" << "_" << Dim; + std::string name{name_stream.str()}; + using Field_t = internal::FieldBase; + py::class_(mod, name.c_str()) + .def_property_readonly("name", &Field_t::get_name, "field name") + .def_property_readonly("collection", &Field_t::get_collection, + "Collection containing this field") + .def_property_readonly("nb_components", &Field_t::get_nb_components, + "number of scalars stored per pixel in this field") + .def_property_readonly("stored_type", + [](const Field_t & field) { + return field.get_stored_typeid().name(); + }, + "fundamental type of scalars stored in this field") + .def_property_readonly("size", &Field_t::size, "number of pixels in this field") + .def("set_zero", &Field_t::set_zero, + "Set all components in the field to zero"); + + add_field(mod, "Real"); + add_field(mod, "Int"); +} + +template +void add_statefield(py::module & mod, std::string dtype_name) { + using StateField_t = TypedStateField; + std::stringstream name_stream{}; + name_stream << (FieldCollection::Global ? "Global" : "Local") + << "StateField" << dtype_name + << "_" << FieldCollection::spatial_dim(); + std::string name{name_stream.str()}; + py::class_(mod, name.c_str()) + .def("get_current_field", &StateField_t::get_current_field, + "returns the current field value", + py::return_value_policy::reference_internal) + .def("get_old_field", &StateField_t::get_old_field, + "nb_steps_ago"_a = 1, + "returns the value this field held 'nb_steps_ago' steps ago", + py::return_value_policy::reference_internal); +} + +template +void add_statefield_helper(py::module & mod) { + std::stringstream name_stream{}; + name_stream << (FieldCollection::Global ? "Global" : "Local") + << "StateField" << "_" << Dim; + std::string name{name_stream.str()}; + using StateField_t = StateFieldBase; + py::class_(mod, name.c_str()) + .def_property_readonly("prefix", &StateField_t::get_prefix, "state field prefix") + .def_property_readonly("collection", &StateField_t::get_collection, + "Collection containing this field") + .def_property_readonly("nb_memory", &StateField_t::get_nb_memory, + "number of old states stored") + .def_property_readonly("stored_type", + [](const StateField_t & field) { + return field.get_stored_typeid().name(); + }, + "fundamental type of scalars stored in this field"); + + add_statefield(mod, "Real"); + add_statefield(mod, "Int"); +} + +template +void add_field_collection_helper(py::module & mod) { + add_field_helper>(mod); + add_field_helper>(mod); + + add_statefield_helper>(mod); + add_statefield_helper>(mod); + + add_field_collection>(mod); + add_field_collection>(mod); +} + +void add_field_collections (py::module & mod) { + add_field_collection_helper< twoD>(mod); + add_field_collection_helper(mod); +} diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc index e3a1e1b..01fb87f 100644 --- a/language_bindings/python/bind_py_material.cc +++ b/language_bindings/python/bind_py_material.cc @@ -1,163 +1,180 @@ /** * @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 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 "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; namespace py = pybind11; using namespace pybind11::literals; /** * 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()) + 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()) + 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()) + 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()) + 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); } 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_module.cc b/language_bindings/python/bind_py_module.cc index 031386c..7302b5d 100644 --- a/language_bindings/python/bind_py_module.cc +++ b/language_bindings/python/bind_py_module.cc @@ -1,43 +1,44 @@ /** * @file bind_py_module.cc * * @author Till Junge * * @date 12 Jan 2018 * * @brief Python bindings for µSpectre * * 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 "bind_py_declarations.hh" #include using namespace pybind11::literals; namespace py=pybind11; PYBIND11_MODULE(_muSpectre, mod) { mod.doc() = "Python bindings to the µSpectre library"; add_common(mod); add_cell(mod); add_material(mod); add_solvers(mod); add_fft_engines(mod); + add_field_collections(mod); } diff --git a/src/common/field.hh b/src/common/field.hh index 16b9df2..af4a289 100644 --- a/src/common/field.hh +++ b/src/common/field.hh @@ -1,662 +1,662 @@ /** * @file field.hh * * @author Till Junge * * @date 07 Sep 2017 * * @brief header-only implementation of a field for field collections * * Copyright © 2017 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. */ #ifndef FIELD_H #define FIELD_H #include "common/T4_map_proxy.hh" #include "field_typed.hh" #include #include #include #include #include #include #include #include #include #include namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ //! declaraton for friending template class FieldMap; /* ---------------------------------------------------------------------- */ /** * A `TypedSizedFieldBase` is the base class for fields that contain a * statically known number of scalars of a statically known type per pixel * in a `FieldCollection`. The actual data for all pixels is * stored in `TypedSizeFieldBase::values`. * `TypedSizedFieldBase` is the base class for `MatrixField` and * `TensorField`. */ template class TypedSizedFieldBase: public TypedField { friend class FieldMap; friend class FieldMap; public: //! for compatibility checks constexpr static auto nb_components{NbComponents}; using Parent = TypedField; //!< base class using Scalar = T; //!< for type checking using Base = typename Parent::Base; //!< root base class //! type stored if ArrayStore is true using Stored_t = Eigen::Array; //! storage container using Storage_t = typename Parent::Storage_t; //! Plain type that is being mapped (Eigen lingo) using EigenRep_t = Eigen::Array; //! maps returned when iterating over field using EigenMap_t = Eigen::Map; //! maps returned when iterating over field using ConstEigenMap_t = Eigen::Map; //! constructor TypedSizedFieldBase(std::string unique_name, FieldCollection& collection); virtual ~TypedSizedFieldBase() = default; //! add a new value at the end of the field inline void push_back(const Stored_t & value); //! add a new scalar value at the end of the field template inline std::enable_if_t push_back(const T & value); /** * returns an upcasted reference to a field, or throws an * exception if the field is incompatible */ static TypedSizedFieldBase & check_ref(Base & other); /** * returns an upcasted reference to a field, or throws an * exception if the field is incompatible */ static const TypedSizedFieldBase & check_ref(const Base & other); //! return a map representing the entire field as a single `Eigen::Array` inline EigenMap_t eigen(); //! return a map representing the entire field as a single `Eigen::Array` inline ConstEigenMap_t eigen() const; /** * return a map representing the entire field as a single * dynamically sized `Eigen::Array` (for python bindings) */ inline typename Parent::EigenMap_t dyn_eigen() {return Parent::eigen();} //! inner product between compatible fields template inline Real inner_product(const TypedSizedFieldBase< FieldCollection, T2, NbComponents> & other) const; protected: //! returns a raw pointer to the entry, for `Eigen::Map` inline T* get_ptr_to_entry(const size_t&& index); //! returns a raw pointer to the entry, for `Eigen::Map` inline const T* get_ptr_to_entry(const size_t&& index) const; }; } // internal /* ---------------------------------------------------------------------- */ /** * The `TensorField` is a subclass of `muSpectre::internal::TypedSizedFieldBase` * that represents tensorial fields, i.e. arbitrary-dimensional arrays with * identical number of rows/columns (that typically correspond to the spatial * cartesian dimensions). It is defined by the stored scalar type @a T, the * tensorial order @a order (often also called degree or rank) and the * number of spatial dimensions @a dim. */ template class TensorField: public internal::TypedSizedFieldBase { public: //! base class using Parent = internal::TypedSizedFieldBase; using Base = typename Parent::Base; //!< root base class //! polymorphic base class using Field_p = typename FieldCollection::Field_p; using Scalar = typename Parent::Scalar; //!< for type checking //! Copy constructor TensorField(const TensorField &other) = delete; //! Move constructor TensorField(TensorField &&other) = delete; //! Destructor virtual ~TensorField() = default; //! Copy assignment operator TensorField& operator=(const TensorField &other) = delete; //! Move assignment operator TensorField& operator=(TensorField &&other) = delete; //! return the order of the stored tensor inline Dim_t get_order() const; //! return the dimension of the stored tensor inline Dim_t get_dim() const; //! factory function template friend FieldType& make_field(std::string unique_name, CollectionType & collection, Args&&... args); //! return a reference or throw an exception if `other` is incompatible static TensorField & check_ref(Base & other) { return static_cast(Parent::check_ref(other));} //! return a reference or throw an exception if `other` is incompatible static const TensorField & check_ref(const Base & other) { return static_cast(Parent::check_ref(other));} /** * Convenience functions to return a map onto this field. A map allows * iteration over all pixels. The map's iterator returns an object that * represents the underlying mathematical structure of the field and * implements common linear algebra operations on it. * Specifically, this function returns * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial * order @a order is unity. * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial * order @a order is 2. * - A `T4MatrixFieldMap` if the tensorial order is 4. */ decltype(auto) get_map(); /** * Convenience functions to return a map onto this field. A map allows * iteration over all pixels. The map's iterator returns an object that * represents the underlying mathematical structure of the field and * implements common linear algebra operations on it. * Specifically, this function returns * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial * order @a order is unity. * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial * order @a order is 2. * - A `T4MatrixFieldMap` if the tensorial order is 4. */ decltype(auto) get_const_map(); /** * Convenience functions to return a map onto this field. A map allows * iteration over all pixels. The map's iterator returns an object that * represents the underlying mathematical structure of the field and * implements common linear algebra operations on it. * Specifically, this function returns * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial * order @a order is unity. * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial * order @a order is 2. * - A `T4MatrixFieldMap` if the tensorial order is 4. */ decltype(auto) get_map() const; protected: //! constructor protected! TensorField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ /** * The `MatrixField` is subclass of `muSpectre::internal::TypedSizedFieldBase` * that represents matrix fields, i.e. a two dimensional arrays, defined by * the stored scalar type @a T and the number of rows @a NbRow and columns * @a NbCol of the matrix. */ template class MatrixField: public internal::TypedSizedFieldBase { public: //! base class using Parent = internal::TypedSizedFieldBase; using Base = typename Parent::Base; //!< root base class //! polymorphic base field ptr to store using Field_p = std::unique_ptr>; //! Copy constructor MatrixField(const MatrixField &other) = delete; //! Move constructor MatrixField(MatrixField &&other) = delete; //! Destructor virtual ~MatrixField() = default; //! Copy assignment operator MatrixField& operator=(const MatrixField &other) = delete; //! Move assignment operator MatrixField& operator=(MatrixField &&other) = delete; //! returns the number of rows inline Dim_t get_nb_row() const; //! returns the number of columns inline Dim_t get_nb_col() const; //! factory function template friend FieldType& make_field(std::string unique_name, CollectionType & collection, Args&&... args); //! returns a `MatrixField` reference if `other` is a compatible field static MatrixField & check_ref(Base & other) { return static_cast(Parent::check_ref(other));} //! returns a `MatrixField` reference if `other` is a compatible field static const MatrixField & check_ref(const Base & other) { return static_cast(Parent::check_ref(other));} /** * Convenience functions to return a map onto this field. A map allows * iteration over all pixels. The map's iterator returns an object that * represents the underlying mathematical structure of the field and * implements common linear algebra operations on it. * Specifically, this function returns * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity. * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns * otherwise. */ decltype(auto) get_map(); /** * Convenience functions to return a map onto this field. A map allows * iteration over all pixels. The map's iterator returns an object that * represents the underlying mathematical structure of the field and * implements common linear algebra operations on it. * Specifically, this function returns * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity. * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns * otherwise. */ decltype(auto) get_const_map(); /** * Convenience functions to return a map onto this field. A map allows * iteration over all pixels. The map's iterator returns an object that * represents the underlying mathematical structure of the field and * implements common linear algebra operations on it. * Specifically, this function returns * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity. * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns * otherwise. */ decltype(auto) get_map() const; protected: //! constructor protected! MatrixField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ //! convenience alias ( template using ScalarField = MatrixField; /* ---------------------------------------------------------------------- */ // Implementations /* ---------------------------------------------------------------------- */ namespace internal { /* ---------------------------------------------------------------------- */ template TypedSizedFieldBase:: TypedSizedFieldBase(std::string unique_name, FieldCollection & collection) :Parent(unique_name, collection, NbComponents){ static_assert ((std::is_arithmetic::value || std::is_same::value), "Use TypedSizedFieldBase for integer, real or complex scalars for T"); static_assert(NbComponents > 0, "Only fields with more than 0 components"); } /* ---------------------------------------------------------------------- */ template TypedSizedFieldBase & TypedSizedFieldBase:: check_ref(Base & other) { if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) { std::string err ="Cannot create a Reference of requested type " +( "for field '" + other.get_name() + "' of type '" + other.get_stored_typeid().name() + "'"); throw std::runtime_error (err); } //check size compatibility if (NbComponents != other.get_nb_components()) { throw std::runtime_error ("Cannot create a Reference to a field with " + std::to_string(NbComponents) + " components " + "for field '" + other.get_name() + "' with " + std::to_string(other.get_nb_components()) + " components"); } return static_cast(other); } /* ---------------------------------------------------------------------- */ template const TypedSizedFieldBase & TypedSizedFieldBase:: check_ref(const Base & other) { if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) { std::stringstream err_str{}; err_str << "Cannot create a Reference of requested type " << "for field '" << other.get_name() << "' of type '" << other.get_stored_typeid().name() << "'"; throw std::runtime_error (err_str.str()); } //check size compatibility if (NbComponents != other.get_nb_components()) { throw std::runtime_error ("Cannot create a Reference to a field with " + std::to_string(NbComponents) + " components " + "for field '" + other.get_name() + "' with " + std::to_string(other.get_nb_components()) + " components"); } return static_cast(other); } /* ---------------------------------------------------------------------- */ template auto TypedSizedFieldBase:: eigen() -> EigenMap_t{ return EigenMap_t(this->data(), NbComponents, this->size()); } /* ---------------------------------------------------------------------- */ template auto TypedSizedFieldBase:: eigen() const -> ConstEigenMap_t{ return ConstEigenMap_t(this->data(), NbComponents, this->size()); } /* ---------------------------------------------------------------------- */ template template Real TypedSizedFieldBase:: inner_product(const TypedSizedFieldBase & other) const { return (this->eigen() * other.eigen()).sum(); } /* ---------------------------------------------------------------------- */ template T* TypedSizedFieldBase:: get_ptr_to_entry(const size_t&& index) { return this->data_ptr + NbComponents*std::move(index); } /* ---------------------------------------------------------------------- */ template const T* TypedSizedFieldBase:: get_ptr_to_entry(const size_t&& index) const { return this->data_ptr + NbComponents*std::move(index); } /* ---------------------------------------------------------------------- */ template void TypedSizedFieldBase:: push_back(const Stored_t & value) { static_assert (not FieldCollection::Global, "You can only push_back data into local field " "collections"); for (Dim_t i = 0; i < NbComponents; ++i) { this->values.push_back(value(i)); } ++this->current_size; this->data_ptr = &this->values.front(); } /* ---------------------------------------------------------------------- */ template template std::enable_if_t TypedSizedFieldBase:: push_back(const T & value) { static_assert(scalar_store, "SFINAE"); this->values.push_back(value); ++this->current_size; this->data_ptr = &this->values.front(); } } // internal /* ---------------------------------------------------------------------- */ //! Factory function, guarantees that only fields get created that are //! properly registered and linked to a collection. template inline FieldType & make_field(std::string unique_name, FieldCollection & collection, Args&&... args) { std::unique_ptr ptr{ new FieldType(unique_name, collection, args...)}; auto& retref{*ptr}; collection.register_field(std::move(ptr)); return retref; } /* ---------------------------------------------------------------------- */ template TensorField:: TensorField(std::string unique_name, FieldCollection & collection) :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_order() const { return order; } /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_dim() const { return dim; } /* ---------------------------------------------------------------------- */ template MatrixField:: MatrixField(std::string unique_name, FieldCollection & collection) :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_col() const { return NbCol; } /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_row() const { return NbRow; } } // muSpectre #include "common/field_map.hh" namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ /** * defines the default mapped type obtained when calling * `muSpectre::TensorField::get_map()` */ template struct tensor_map_type { }; /// specialisation for vectors template struct tensor_map_type { //! use this type using type = MatrixFieldMap; }; /// specialisation to second-order tensors (matrices) template struct tensor_map_type { //! use this type using type = MatrixFieldMap; }; /// specialisation to fourth-order tensors template struct tensor_map_type { //! use this type using type = T4MatrixFieldMap; }; /* ---------------------------------------------------------------------- */ /** * defines the default mapped type obtained when calling * `muSpectre::MatrixField::get_map()` */ template struct matrix_map_type { //! mapping type using type = MatrixFieldMap; }; //! specialisation to scalar fields template struct matrix_map_type { //! mapping type using type = ScalarFieldMap; }; } // internal /* ---------------------------------------------------------------------- */ template - decltype(auto) TensorField:: - get_map() { + auto TensorField:: + get_map() -> decltype(auto) { constexpr bool map_constness{false}; using RawMap_t = typename internal::tensor_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template - decltype(auto) TensorField:: - get_const_map() { + auto TensorField:: + get_const_map() -> decltype(auto) { constexpr bool map_constness{true}; using RawMap_t = typename internal::tensor_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template - decltype(auto) TensorField:: - get_map() const { + auto TensorField:: + get_map() const -> decltype(auto) { constexpr bool map_constness{true}; using RawMap_t = typename internal::tensor_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template - decltype(auto) MatrixField:: - get_map() { + auto MatrixField:: + get_map() -> decltype(auto) { constexpr bool map_constness{false}; using RawMap_t = typename internal::matrix_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template - decltype(auto) MatrixField:: - get_const_map() { + auto MatrixField:: + get_const_map() -> decltype(auto) { constexpr bool map_constness{true}; using RawMap_t = typename internal::matrix_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template - decltype(auto) MatrixField:: - get_map() const { + auto MatrixField:: + get_map() const -> decltype(auto) { constexpr bool map_constness{true}; using RawMap_t = typename internal::matrix_map_type::type; return RawMap_t(*this); } } // muSpectre #endif /* FIELD_H */ diff --git a/src/common/field_base.hh b/src/common/field_base.hh index 2c29501..95f825e 100644 --- a/src/common/field_base.hh +++ b/src/common/field_base.hh @@ -1,211 +1,210 @@ /** * file field_base.hh * * @author Till Junge * * @date 10 Apr 2018 * * @brief Virtual base class for fields * - * @section LICENSE - * * 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. */ #ifndef FIELD_BASE_H #define FIELD_BASE_H #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ /** * base class for field collection-related exceptions */ class FieldCollectionError: public std::runtime_error { public: //! constructor explicit FieldCollectionError(const std::string& what) :std::runtime_error(what){} //! constructor explicit FieldCollectionError(const char * what) :std::runtime_error(what){} }; /// base class for field-related exceptions class FieldError: public FieldCollectionError { using Parent = FieldCollectionError; public: //! constructor explicit FieldError(const std::string& what) :Parent(what){} //! constructor explicit FieldError(const char * what) :Parent(what){} }; /** * Thrown when a associating a field map to and incompatible field * is attempted */ class FieldInterpretationError: public FieldError { public: //! constructor explicit FieldInterpretationError(const std::string & what) :FieldError(what){} //! constructor explicit FieldInterpretationError(const char * what) :FieldError(what){} }; namespace internal{ /* ---------------------------------------------------------------------- */ /** * Virtual base class for all fields. A field represents * meta-information for the per-pixel storage for a scalar, vector * or tensor quantity and is therefore the abstract class defining * the field. It is used for type and size checking at runtime and * for storage of polymorphic pointers to fully typed and sized * fields. `FieldBase` (and its children) are templated with a * specific `FieldCollection` (derived from * `muSpectre::FieldCollectionBase`). A `FieldCollection` stores * multiple fields that all apply to the same set of * pixels. Addressing and managing the data for all pixels is * handled by the `FieldCollection`. Note that `FieldBase` does * not know anything about about mathematical operations on the * data or how to iterate over all pixels. Mapping the raw data * onto for instance Eigen maps and iterating over those is * handled by the `FieldMap`. */ template class FieldBase { protected: //! constructor //! unique name (whithin Collection) //! number of components //! collection to which this field belongs (eg, material, cell) FieldBase(std::string unique_name, size_t nb_components, FieldCollection & collection); public: using collection_t = FieldCollection; //!< for type checks //! Copy constructor FieldBase(const FieldBase &other) = delete; //! Move constructor FieldBase(FieldBase &&other) = delete; //! Destructor virtual ~FieldBase() = default; //! Copy assignment operator FieldBase& operator=(const FieldBase &other) = delete; //! Move assignment operator FieldBase& operator=(FieldBase &&other) = delete; /* ---------------------------------------------------------------------- */ //!Identifying accessors //! return field name inline const std::string & get_name() const; //! return field type //inline const Field_t & get_type() const; //! return my collection (for iterating) inline const FieldCollection & get_collection() const; - //! return my collection (for iterating) + //! return number of components (e.g., dimensions) of this field inline const size_t & get_nb_components() const; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const = 0; //! number of pixels in the field virtual size_t size() const = 0; //! add a pad region to the end of the field buffer; required for //! using this as e.g. an FFT workspace virtual void set_pad_size(size_t pad_size_) = 0; //! pad region size virtual size_t get_pad_size() const {return this->pad_size;}; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) virtual void set_zero() = 0; //! give access to collections friend FieldCollection; //! give access to collection's base class - friend typename FieldCollection::Parent; + using FParent_t = typename FieldCollection::Parent; + friend FParent_t; protected: /* ---------------------------------------------------------------------- */ //! allocate memory etc virtual void resize(size_t size) = 0; const std::string name; //!< the field's unique name const size_t nb_components; //!< number of components per entry //! reference to the collection this field belongs to const FieldCollection & collection; size_t pad_size; //!< size of padding region at end of buffer private: }; /* ---------------------------------------------------------------------- */ // Implementations /* ---------------------------------------------------------------------- */ template FieldBase::FieldBase(std::string unique_name, size_t nb_components_, FieldCollection & collection_) :name(unique_name), nb_components(nb_components_), collection(collection_), pad_size{0} {} /* ---------------------------------------------------------------------- */ template inline const std::string & FieldBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template inline const FieldCollection & FieldBase:: get_collection() const { return this->collection; } /* ---------------------------------------------------------------------- */ template inline const size_t & FieldBase:: get_nb_components() const { return this->nb_components; } } // internal } // muSpectre #endif /* FIELD_BASE_H */ diff --git a/src/common/field_collection_base.hh b/src/common/field_collection_base.hh index ddb4c95..1e8e333 100644 --- a/src/common/field_collection_base.hh +++ b/src/common/field_collection_base.hh @@ -1,181 +1,351 @@ /** * @file field_collection_base.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief Base class for field collections * * Copyright © 2017 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. */ #ifndef FIELD_COLLECTION_BASE_H #define FIELD_COLLECTION_BASE_H #include "common/common.hh" #include "common/field.hh" +#include "common/statefield.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ /** `FieldCollectionBase` is the base class for collections of fields. All * fields in a field collection have the same number of pixels. The field * collection is templated with @a DimS is the spatial dimension (i.e. * whether the simulation domain is one, two or three-dimensional). * All fields within a field collection have a unique string identifier. * A `FieldCollectionBase` is therefore comparable to a dictionary of fields * that live on the same grid. * `FieldCollectionBase` has the specialisations `GlobalFieldCollection` and * `LocalFieldCollection`. */ template class FieldCollectionBase { public: //! polymorphic base type to store - using Field = internal::FieldBase; - using Field_p = std::unique_ptr; //!< stored type + using Field_t = internal::FieldBase; + template + using TypedField_t = TypedField; + + using Field_p = std::unique_ptr; //!< stored type + using StateField_t = StateFieldBase; + template + using TypedStateField_t = TypedStateField; + + using StateField_p = std::unique_ptr; using Ccoord = Ccoord_t; //!< cell coordinates type //! Default constructor FieldCollectionBase(); //! Copy constructor FieldCollectionBase(const FieldCollectionBase &other) = delete; //! Move constructor FieldCollectionBase(FieldCollectionBase &&other) = delete; //! Destructor virtual ~FieldCollectionBase() = default; //! Copy assignment operator FieldCollectionBase& operator=(const FieldCollectionBase &other) = delete; //! Move assignment operator FieldCollectionBase& operator=(FieldCollectionBase &&other) = delete; //! Register a new field (fields need to be in heap, so I want to keep them //! as shared pointers void register_field(Field_p&& field); + //! Register a new field (fields need to be in heap, so I want to keep them + //! as shared pointers + void register_statefield(StateField_p&& field); + //! for return values of iterators constexpr inline static Dim_t spatial_dim(); //! for return values of iterators inline Dim_t get_spatial_dim() const; + //! return names of all stored fields + std::vector get_field_names() const { + std::vector names{}; + for (auto & tup: this->fields) { + names.push_back(std::get<0>(tup)); + } + return names; + } + + //! return names of all state fields + std::vector get_statefield_names() const { + std::vector names{}; + for (auto & tup: this->statefields) { + names.push_back(std::get<0>(tup)); + } + return names; + } + //! retrieve field by unique_name - inline Field& operator[](std::string unique_name); + inline Field_t& operator[](std::string unique_name); //! retrieve field by unique_name with bounds checking - inline Field& at(std::string unique_name); + inline Field_t& at(std::string unique_name); + + //! retrieve typed field by unique_name + template + inline TypedField_t & get_typed_field(std::string unique_name); + + //! retrieve state field by unique_prefix with bounds checking + template + inline TypedStateField_t& get_typed_statefield(std::string unique_prefix); + + //! retrieve state field by unique_prefix with bounds checking + inline StateField_t& get_statefield(std::string unique_prefix) { + return *(this->statefields.at(unique_prefix)); + } + + //! retrieve state field by unique_prefix with bounds checking + inline const StateField_t& get_statefield(std::string unique_prefix) const { + return *(this->statefields.at(unique_prefix)); + } + + /** + * retrieve current value of typed state field by unique_prefix with + * bounds checking + */ + template + inline TypedField_t& get_current(std::string unique_prefix); + + /** + * retrieve old value of typed state field by unique_prefix with + * bounds checking + */ + template + inline const TypedField_t & get_old(std::string unique_prefix, + size_t nb_steps_ago = 1) const; //! returns size of collection, this refers to the number of pixels handled //! by the collection, not the number of fields inline size_t size() const {return this->size_;} //! check whether a field is present bool check_field_exists(std::string unique_name); protected: std::map fields{}; //!< contains the field ptrs + //! contains ptrs to state fields + std::map statefields{}; bool is_initialised{false}; //!< to handle double initialisation correctly const Uint id; //!< unique identifier static Uint counter; //!< used to assign unique identifiers size_t size_{0}; //!< holds the number of pixels after initialisation private: }; /* ---------------------------------------------------------------------- */ template Uint FieldCollectionBase::counter{0}; /* ---------------------------------------------------------------------- */ template FieldCollectionBase::FieldCollectionBase() :id(counter++){} /* ---------------------------------------------------------------------- */ template - void FieldCollectionBase::register_field(Field_p &&field) { + void FieldCollectionBase:: + register_field(Field_p &&field) { auto&& search_it = this->fields.find(field->get_name()); auto&& does_exist = search_it != this->fields.end(); if (does_exist) { std::stringstream err_str; err_str << "a field named " << field->get_name() << "is already registered in this field collection. " << "Currently registered fields: "; for (const auto& name_field_pair: this->fields) { err_str << ", " << name_field_pair.first; } throw FieldCollectionError(err_str.str()); } if (this->is_initialised) { field->resize(this->size()); } this->fields[field->get_name()] = std::move(field); } + /* ---------------------------------------------------------------------- */ + template + void FieldCollectionBase:: + register_statefield(StateField_p&& field) { + auto&& search_it = this->statefields.find(field->get_prefix()); + auto&& does_exist = search_it != this->statefields.end(); + if (does_exist) { + std::stringstream err_str; + err_str << "a state field named " << field->get_prefix() + << "is already registered in this field collection. " + << "Currently registered fields: "; + for (const auto& name_field_pair: this->statefields) { + err_str << ", " << name_field_pair.first; + } + throw FieldCollectionError(err_str.str()); + } + this->statefields[field->get_prefix()] = std::move(field); + } + /* ---------------------------------------------------------------------- */ template constexpr Dim_t FieldCollectionBase:: spatial_dim() { return DimS; } /* ---------------------------------------------------------------------- */ template Dim_t FieldCollectionBase:: get_spatial_dim() const { return DimS; } /* ---------------------------------------------------------------------- */ template - typename FieldCollectionBase::Field& + auto FieldCollectionBase:: - operator[](std::string unique_name) { + operator[](std::string unique_name) -> Field_t & { return *(this->fields[unique_name]); } /* ---------------------------------------------------------------------- */ template - typename FieldCollectionBase::Field& + auto FieldCollectionBase:: - at(std::string unique_name) { + at(std::string unique_name) -> Field_t & { return *(this->fields.at(unique_name)); } /* ---------------------------------------------------------------------- */ template bool FieldCollectionBase:: check_field_exists(std::string unique_name) { return this->fields.find(unique_name) != this->fields.end(); } + //! retrieve typed field by unique_name + template + template + auto FieldCollectionBase:: + get_typed_field(std::string unique_name) -> TypedField_t & { + auto & unqualified_field{this->at(unique_name)}; + if (unqualified_field.get_stored_typeid().hash_code() != + typeid(T).hash_code()) { + std::stringstream err{}; + err << "Field '" << unique_name << "' is of type " + << unqualified_field.get_stored_typeid().name() + << ", but should be of type " << typeid(T).name() << std::endl; + throw FieldCollectionError(err.str()); + } + return static_cast &>(unqualified_field); + } + + /* ---------------------------------------------------------------------- */ + //! retrieve state field by unique_prefix with bounds checking + template + template + auto FieldCollectionBase:: + get_typed_statefield(std::string unique_prefix) + -> TypedStateField_t & { + auto & unqualified_statefield{this->get_statefield(unique_prefix)}; + if (unqualified_statefield.get_stored_typeid().hash_code() != + typeid(T).hash_code()) { + std::stringstream err{}; + err << "Statefield '" << unique_prefix << "' is of type " + << unqualified_statefield.get_stored_typeid().name() + << ", but should be of type " << typeid(T).name() << std::endl; + throw FieldCollectionError(err.str()); + } + return static_cast &>(unqualified_statefield); + } + + + /* ---------------------------------------------------------------------- */ + template + template + auto + FieldCollectionBase:: + get_current(std::string unique_prefix) -> TypedField_t & { + auto & unqualified_statefield = this->get_statefield(unique_prefix); + + //! check for correct underlying fundamental type + if (unqualified_statefield.get_stored_typeid().hash_code() != + typeid(T).hash_code()) { + std::stringstream err{}; + err << "StateField '" << unique_prefix << "' is of type " + << unqualified_statefield.get_stored_typeid().name() + << ", but should be of type " << typeid(T).name() << std::endl; + throw FieldCollectionError(err.str()); + } + + using Typed_t = TypedStateField; + auto & typed_field{static_cast(unqualified_statefield)}; + return typed_field.get_current_field(); + } + + /* ---------------------------------------------------------------------- */ + /* ---------------------------------------------------------------------- */ + template + template + auto + FieldCollectionBase:: + get_old(std::string unique_prefix, size_t nb_steps_ago) const + -> const TypedField_t & { + auto & unqualified_statefield = this->get_statefield(unique_prefix); + + //! check for correct underlying fundamental type + if (unqualified_statefield.get_stored_typeid().hash_code() != + typeid(T).hash_code()) { + std::stringstream err{}; + err << "StateField '" << unique_prefix << "' is of type " + << unqualified_statefield.get_stored_typeid().name() + << ", but should be of type " << typeid(T).name() << std::endl; + throw FieldCollectionError(err.str()); + } + + using Typed_t = TypedStateField; + auto & typed_field{static_cast(unqualified_statefield)}; + return typed_field.get_old_field(nb_steps_ago); + } } // muSpectre #endif /* FIELD_COLLECTION_BASE_H */ diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh index 2b5d370..da1388b 100644 --- a/src/common/field_map_base.hh +++ b/src/common/field_map_base.hh @@ -1,643 +1,645 @@ /** * @file field_map.hh * * @author Till Junge * * @date 12 Sep 2017 * * @brief Defined a strongly defines proxy that iterates efficiently over a field * * Copyright © 2017 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. */ #ifndef FIELD_MAP_BASE_H #define FIELD_MAP_BASE_H #include "common/field.hh" #include "field_collection_base.hh" #include "common/common.hh" #include #include #include #include namespace muSpectre { namespace internal { //----------------------------------------------------------------------------// //! little helper to automate creation of const maps without duplication template struct const_corrector { //! non-const type using type = typename T::reference; }; //! specialisation for constant case template struct const_corrector { //! const type using type = typename T::const_reference; }; //! convenience alias template using const_corrector_t = typename const_corrector::type; //----------------------------------------------------------------------------// /** * `FieldMap` provides two central mechanisms: * - Map a field (that knows only about the size of the underlying object, * onto the mathematical object (reprensented by the respective Eigen class) * that provides linear algebra functionality. * - Provide an iterator that allows to iterate over all pixels. * A field is represented by `FieldBase` or a derived class. * `FieldMap` has the specialisations `MatrixLikeFieldMap`, * `ScalarFieldMap` and `TensorFieldMap`. */ template class FieldMap { public: + //! Fundamental type stored + using Scalar = T; //! number of scalars per entry constexpr static auto nb_components{NbComponents}; using TypedField_nc = TypedSizedFieldBase ; //!< non-constant version of field //! field type as seen from iterator using TypedField = std::conditional_t; using Field = typename TypedField::Base; //!< iterated field type //! const-correct field type using Field_c = std::conditional_t; using size_type = std::size_t; //!< stl conformance using pointer = std::conditional_t; //!< stl conformance //! Default constructor FieldMap() = delete; //! constructor FieldMap(Field_c & field); //! constructor with run-time cost (for python and debugging) template FieldMap(TypedSizedFieldBase & field); //! Copy constructor FieldMap(const FieldMap &other) = default; //! Move constructor FieldMap(FieldMap &&other) = default; //! Destructor virtual ~FieldMap() = default; //! Copy assignment operator FieldMap& operator=(const FieldMap &other) = delete; //! Move assignment operator FieldMap& operator=(FieldMap &&other) = delete; //! give human-readable field map type virtual std::string info_string() const = 0; //! return field name inline const std::string & get_name() const; //! return my collection (for iterating) inline const FieldCollection & get_collection() const; //! member access needs to be implemented by inheriting classes //inline value_type operator[](size_t index); //inline value_type operator[](Ccoord ccord); //! check compatibility (must be called by all inheriting classes at the //! end of their constructors inline void check_compatibility(); //! convenience call to collection's size method inline size_t size() const; //! compile-time compatibility check template struct is_compatible; /** * iterates over all pixels in the `muSpectre::FieldCollection` * and dereferences to an Eigen map to the currently used field. */ template class iterator { static_assert(!((ConstIter==false) && (ConstField==true)), "You can't have a non-const iterator over a const " "field"); public: //! stl conformance using value_type = const_corrector_t; //! stl conformance using const_value_type = const_corrector_t; //! stl conformance using pointer = typename FullyTypedFieldMap::pointer; //! stl conformance using difference_type = std::ptrdiff_t; //! stl conformance using iterator_category = std::random_access_iterator_tag; //! cell coordinates type using Ccoord = typename FieldCollection::Ccoord; //! stl conformance using reference = typename FullyTypedFieldMap::reference; //! fully typed reference as seen by the iterator using TypedRef = std::conditional_t; //! Default constructor iterator() = delete; //! constructor inline iterator(TypedRef fieldmap, bool begin=true); //! constructor for random access inline iterator(TypedRef fieldmap, size_t index); //! Copy constructor iterator(const iterator &other)= default; //! Move constructor iterator(iterator &&other) = default; //! Destructor virtual ~iterator() = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) = default; //! pre-increment inline iterator & operator++(); //! post-increment inline iterator operator++(int); //! dereference inline value_type operator*(); //! dereference inline const_value_type operator*() const; //! member of pointer inline pointer operator->(); //! pre-decrement inline iterator & operator--(); //! post-decrement inline iterator operator--(int); //! access subscripting inline value_type operator[](difference_type diff); //! access subscripting inline const_value_type operator[](const difference_type diff) const; //! equality inline bool operator==(const iterator & other) const; //! inequality inline bool operator!=(const iterator & other) const; //! div. comparisons inline bool operator<(const iterator & other) const; //! div. comparisons inline bool operator<=(const iterator & other) const; //! div. comparisons inline bool operator>(const iterator & other) const; //! div. comparisons inline bool operator>=(const iterator & other) const; //! additions, subtractions and corresponding assignments inline iterator operator+(difference_type diff) const; //! additions, subtractions and corresponding assignments inline iterator operator-(difference_type diff) const; //! additions, subtractions and corresponding assignments inline iterator& operator+=(difference_type diff); //! additions, subtractions and corresponding assignments inline iterator& operator-=(difference_type diff); //! get pixel coordinates inline Ccoord get_ccoord() const; //! ostream operator (mainly for debug friend std::ostream & operator<<(std::ostream & os, const iterator& it) { if (ConstIter) { os << "const "; } os << "iterator on field '" << it.fieldmap.get_name() << "', entry " << it.index; return os; } protected: const FieldCollection & collection; //!< collection of the field TypedRef fieldmap; //!< ref to the field itself size_t index; //!< index of currently pointed-to pixel private: }; TypedField & get_field() {return this->field;} protected: //! raw pointer to entry (for Eigen Map) inline pointer get_ptr_to_entry(size_t index); //! raw pointer to entry (for Eigen Map) inline const T* get_ptr_to_entry(size_t index) const; const FieldCollection & collection; //!< collection holding Field TypedField & field; //!< mapped Field private: }; } // internal namespace internal { /* ---------------------------------------------------------------------- */ template FieldMap:: FieldMap(Field_c& field) :collection(field.get_collection()), field(static_cast(field)) { static_assert(NbComponents>0, "Only fields with more than 0 components allowed"); } /* ---------------------------------------------------------------------- */ template template FieldMap:: FieldMap(TypedSizedFieldBase & field) :collection(field.get_collection()), field(static_cast(field)) { static_assert(std::is_same::value, "The field does not have the expected FieldCollection type"); static_assert(std::is_same::value, "The field does not have the expected Scalar type"); static_assert((NbC == NbComponents), "The field does not have the expected number of components"); } /* ---------------------------------------------------------------------- */ template void FieldMap:: check_compatibility() { if (typeid(T).hash_code() != this->field.get_stored_typeid().hash_code()) { std::string err{"Cannot create a Map of type '" + this->info_string() + "' for field '" + this->field.get_name() + "' of type '" + this->field.get_stored_typeid().name() + "'"}; throw FieldInterpretationError (err); } //check size compatibility if (NbComponents != this->field.get_nb_components()) { throw FieldInterpretationError ("Cannot create a Map of type '" + this->info_string() + "' for field '" + this->field.get_name() + "' with " + std::to_string(this->field.get_nb_components()) + " components"); } } /* ---------------------------------------------------------------------- */ template size_t FieldMap:: size() const { return this->collection.size(); } /* ---------------------------------------------------------------------- */ template template struct FieldMap::is_compatible { //! creates a more readable compile error constexpr static bool explain() { static_assert (std::is_same::value, "The field does not have the expected FieldCollection type"); static_assert (std::is_same::value, "The // field does not have the expected Scalar type"); static_assert((TypedField::nb_components == NbComponents), "The field does not have the expected number of components"); //The static asserts wouldn't pass in the incompatible case, so this is it return true; } //! evaluated compatibility constexpr static bool value{std::is_base_of::value}; }; /* ---------------------------------------------------------------------- */ template const std::string & FieldMap:: get_name() const { return this->field.get_name(); } /* ---------------------------------------------------------------------- */ template const FieldCollection & FieldMap:: get_collection() const { return this->collection; } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ // Iterator implementations //! constructor template template FieldMap::iterator :: iterator(TypedRef fieldmap, bool begin) :collection(fieldmap.get_collection()), fieldmap(fieldmap), index(begin ? 0 : fieldmap.field.size()) {} /* ---------------------------------------------------------------------- */ //! constructor for random access template template FieldMap::iterator:: iterator(TypedRef fieldmap, size_t index) :collection(fieldmap.collection), fieldmap(fieldmap), index(index) {} /* ---------------------------------------------------------------------- */ //! pre-increment template template typename FieldMap::template iterator & FieldMap::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ //! post-increment template template typename FieldMap::template iterator FieldMap::iterator:: operator++(int) { iterator current = *this; this->index++; return current; } /* ---------------------------------------------------------------------- */ //! dereference template template typename FieldMap:: template iterator::value_type FieldMap:: iterator:: operator*() { return this->fieldmap.operator[](this->index); } /* ---------------------------------------------------------------------- */ //! dereference template template typename FieldMap:: template iterator::const_value_type FieldMap:: iterator:: operator*() const { return this->fieldmap.operator[](this->index); } /* ---------------------------------------------------------------------- */ //! member of pointer template template typename FullyTypedFieldMap::pointer FieldMap:: iterator:: operator->() { return this->fieldmap.ptr_to_val_t(this->index); } /* ---------------------------------------------------------------------- */ //! pre-decrement template template typename FieldMap::template iterator & FieldMap:: iterator:: operator--() { this->index--; return *this; } /* ---------------------------------------------------------------------- */ //! post-decrement template template typename FieldMap::template iterator FieldMap:: iterator:: operator--(int) { iterator current = *this; this->index--; return current; } /* ---------------------------------------------------------------------- */ //! Access subscripting template template typename FieldMap:: template iterator::value_type FieldMap:: iterator:: operator[](difference_type diff) { return this->fieldmap[this->index+diff]; } /* ---------------------------------------------------------------------- */ //! Access subscripting template template typename FieldMap::template iterator::const_value_type FieldMap:: iterator:: operator[](const difference_type diff) const { return this->fieldmap[this->index+diff]; } /* ---------------------------------------------------------------------- */ //! equality template template bool FieldMap:: iterator:: operator==(const iterator & other) const { return (this->index == other.index); } /* ---------------------------------------------------------------------- */ //! inquality template template bool FieldMap:: iterator:: operator!=(const iterator & other) const { return !(*this == other); } /* ---------------------------------------------------------------------- */ //! div. comparisons template template bool FieldMap:: iterator:: operator<(const iterator & other) const { return (this->index < other.index); } template template bool FieldMap:: iterator:: operator<=(const iterator & other) const { return (this->index <= other.index); } template template bool FieldMap:: iterator:: operator>(const iterator & other) const { return (this->index > other.index); } template template bool FieldMap:: iterator:: operator>=(const iterator & other) const { return (this->index >= other.index); } /* ---------------------------------------------------------------------- */ //! additions, subtractions and corresponding assignments template template typename FieldMap::template iterator FieldMap:: iterator:: operator+(difference_type diff) const { return iterator(this->fieldmap, this->index + diff); } template template typename FieldMap::template iterator FieldMap:: iterator:: operator-(difference_type diff) const { return iterator(this->fieldmap, this->index - diff); } template template typename FieldMap::template iterator & FieldMap:: iterator:: operator+=(difference_type diff) { this->index += diff; return *this; } template template typename FieldMap::template iterator & FieldMap:: iterator:: operator-=(difference_type diff) { this->index -= diff; return *this; } /* ---------------------------------------------------------------------- */ //! get pixel coordinates template template typename FieldCollection::Ccoord FieldMap:: iterator:: get_ccoord() const { return this->collection.get_ccoord(this->index); } ////----------------------------------------------------------------------------// //template //std::ostream & operator << //(std::ostream &os, // const typename FieldMap:: // template iterator & it) { // os << "iterator on field '" // << it.field.get_name() // << "', entry " << it.index; // return os; //} /* ---------------------------------------------------------------------- */ template typename FieldMap::pointer FieldMap:: get_ptr_to_entry(size_t index) { return this->field.get_ptr_to_entry(std::move(index)); } /* ---------------------------------------------------------------------- */ template const T* FieldMap:: get_ptr_to_entry(size_t index) const { return this->field.get_ptr_to_entry(std::move(index)); } } // internal } // muSpectre #endif /* FIELD_MAP_BASE_H */ diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh index 76b9c54..e69d18b 100644 --- a/src/common/field_map_matrixlike.hh +++ b/src/common/field_map_matrixlike.hh @@ -1,383 +1,376 @@ /** * @file field_map_matrixlike.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief Eigen-Matrix and -Array maps over strongly typed fields * * Copyright © 2017 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. */ #ifndef FIELD_MAP_MATRIXLIKE_H #define FIELD_MAP_MATRIXLIKE_H #include "common/field_map_base.hh" #include "common/T4_map_proxy.hh" #include #include namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ /** * lists all matrix-like types consideres by * `muSpectre::internal::MatrixLikeFieldMap` */ enum class Map_t{ Matrix, //!< for wrapping `Eigen::Matrix` Array, //!< for wrapping `Eigen::Array` T4Matrix //!< for wrapping `Eigen::T4Matrix` }; /** * traits structure to define the name shown when a * `muSpectre::MatrixLikeFieldMap` output into an ostream */ template struct NameWrapper { }; /// specialisation for `muSpectre::ArrayFieldMap` template<> struct NameWrapper { //! string to use for printing static std::string field_info_root() {return "Array";} }; /// specialisation for `muSpectre::MatrixFieldMap` template<> struct NameWrapper { //! string to use for printing static std::string field_info_root() {return "Matrix";} }; /// specialisation for `muSpectre::T4MatrixFieldMap` template<> struct NameWrapper { //! string to use for printing static std::string field_info_root() {return "T4Matrix";} }; /* ---------------------------------------------------------------------- */ /*! * A `MatrixLikeFieldMap` is the base class for maps of matrices, arrays and * fourth-order tensors mapped onto matrices. * * It should never be necessary to call directly any of the * constructors if this class, but rather use the template aliases: * - `muSpectre::ArrayFieldMap`: iterate in the form of `Eigen::Array<...>`. * - `muSpectre::MatrixFieldMap`: iterate in the form of * `Eigen::Matrix<...>`. * - `muSpectre::T4MatrixFieldMap`: iterate in the form of `muSpectre::T4MatMap`. */ template class MatrixLikeFieldMap: public FieldMap { public: //! base class using Parent = FieldMap; //! sister class with all params equal, but ConstField guaranteed true using ConstMap = MatrixLikeFieldMap; using T_t = EigenPlain; //!< plain Eigen type to map //! cell coordinates type using Ccoord = Ccoord_t; using value_type = EigenArray; //!< stl conformance using const_reference = EigenConstArray; //!< stl conformance //! stl conformance using reference = std::conditional_t; // since it's a resource handle using size_type = typename Parent::size_type; //!< stl conformance using pointer = std::unique_ptr; //!< stl conformance using TypedField = typename Parent::TypedField; //!< stl conformance using Field = typename Parent::Field; //!< stl conformance using Field_c = typename Parent::Field_c; //!< stl conformance //! stl conformance using const_iterator= typename Parent::template iterator; //! stl conformance using iterator = std::conditional_t< ConstField, const_iterator, typename Parent::template iterator>; using reverse_iterator = std::reverse_iterator; //!< stl conformance //! stl conformance using const_reverse_iterator = std::reverse_iterator; //! stl conformance friend iterator; //! Default constructor MatrixLikeFieldMap() = delete; /** * Constructor using a (non-typed) field. Compatibility is enforced at * runtime. This should not be a performance concern, as this constructor * will not be called in anny inner loops (if used as intended). */ MatrixLikeFieldMap(Field_c & field); /** * Constructor using a typed field. Compatibility is enforced * statically. It is not always possible to call this constructor, as the * type of the field might not be known at compile time. */ template MatrixLikeFieldMap(TypedSizedFieldBase & field); //! Copy constructor MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = default; //! Move constructor MatrixLikeFieldMap(MatrixLikeFieldMap &&other) = default; //! Destructor virtual ~MatrixLikeFieldMap() = default; //! Copy assignment operator MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete; //! Move assignment operator MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) = delete; //! Assign a matrixlike value to every entry template inline MatrixLikeFieldMap & operator=(const Eigen::EigenBase & val); //! give human-readable field map type inline std::string info_string() const override final; //! member access inline reference operator[](size_type index); //! member access inline reference operator[](const Ccoord& ccoord); //! member access inline const_reference operator[](size_type index) const; //! member access inline const_reference operator[](const Ccoord& ccoord) const; //! return an iterator to head of field for ranges inline iterator begin(){return iterator(*this);} //! return an iterator to head of field for ranges inline const_iterator cbegin() const {return const_iterator(*this);} //! return an iterator to head of field for ranges inline const_iterator begin() const {return this->cbegin();} //! return an iterator to tail of field for ranges inline iterator end(){return iterator(*this, false);}; //! return an iterator to tail of field for ranges inline const_iterator cend() const {return const_iterator(*this, false);} //! return an iterator to tail of field for ranges inline const_iterator end() const {return this->cend();} //! evaluate the average of the field inline T_t mean() const; protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); - const static std::string field_info_root; //!< for printing and debug + const static std::string field_info_root() { + return NameWrapper::field_info_root();} //!< for printing and debug private: }; - /* ---------------------------------------------------------------------- */ - template - const std::string MatrixLikeFieldMap:: - field_info_root{NameWrapper::field_info_root()}; - /* ---------------------------------------------------------------------- */ template MatrixLikeFieldMap:: MatrixLikeFieldMap(Field_c & field) :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ template template MatrixLikeFieldMap:: MatrixLikeFieldMap(TypedSizedFieldBase & field) :Parent(field) { } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string MatrixLikeFieldMap:: info_string() const { std::stringstream info; - info << field_info_root << "(" + info << this->field_info_root() << "(" << typeid(typename EigenArray::value_type).name() << ", " << EigenArray::RowsAtCompileTime << "x" << EigenArray::ColsAtCompileTime << ")"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template typename MatrixLikeFieldMap::reference MatrixLikeFieldMap:: operator[](size_type index) { return reference(this->get_ptr_to_entry(index)); } template typename MatrixLikeFieldMap::reference MatrixLikeFieldMap:: operator[](const Ccoord & ccoord) { size_t index{}; index = this->collection.get_index(ccoord); return reference(this->get_ptr_to_entry(std::move(index))); } /* ---------------------------------------------------------------------- */ //! member access template typename MatrixLikeFieldMap::const_reference MatrixLikeFieldMap:: operator[](size_type index) const { return const_reference(this->get_ptr_to_entry(index)); } template typename MatrixLikeFieldMap::const_reference MatrixLikeFieldMap:: operator[](const Ccoord & ccoord) const{ size_t index{}; index = this->collection.get_index(ccoord); return const_reference(this->get_ptr_to_entry(std::move(index))); } //----------------------------------------------------------------------------// template template MatrixLikeFieldMap & MatrixLikeFieldMap:: operator=(const Eigen::EigenBase & val) { for (auto && entry: *this) { entry = val; } return *this; } /* ---------------------------------------------------------------------- */ template typename MatrixLikeFieldMap::T_t MatrixLikeFieldMap:: mean() const { T_t mean{T_t::Zero()}; for (auto && val: *this) { mean += val; } mean *= 1./Real(this->size()); return mean; } /* ---------------------------------------------------------------------- */ template typename MatrixLikeFieldMap::pointer MatrixLikeFieldMap:: ptr_to_val_t(size_type index) { return std::make_unique (this->get_ptr_to_entry(std::move(index))); } } // internal /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate template using MatrixFieldMap = internal::MatrixLikeFieldMap >, Eigen::Map>, Eigen::Matrix, internal::Map_t::Matrix, ConstField>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate template using T4MatrixFieldMap = internal::MatrixLikeFieldMap , T4MatMap, T4Mat, internal::Map_t::T4Matrix, MapConst>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen array map as iterate template using ArrayFieldMap = internal::MatrixLikeFieldMap >, Eigen::Map>, Eigen::Array, internal::Map_t::Array, ConstField>; } // muSpectre #endif /* FIELD_MAP_MATRIXLIKE_H */ diff --git a/src/common/field_typed.hh b/src/common/field_typed.hh index b14a5f5..e7f3687 100644 --- a/src/common/field_typed.hh +++ b/src/common/field_typed.hh @@ -1,327 +1,342 @@ /** * file field_typed.hh * * @author Till Junge * * @date 10 Apr 2018 * * @brief Typed Field for dynamically sized fields and base class for fields * of tensors, matrices, etc * - * @section LICENSE - * * 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. */ #ifndef FIELD_TYPED_H #define FIELD_TYPED_H #include "field_base.hh" #include namespace muSpectre { /** * Dummy intermediate class to provide a run-time polymorphic * typed field. Mainly for binding Python. TypedField specifies methods * that return typed Eigen maps and vectors in addition to pointers to the * raw data. */ template class TypedField: public internal::FieldBase { public: using Parent = internal::FieldBase; //!< base class //! for type checks when mapping this field using collection_t = typename Parent::collection_t; using Scalar = T; //!< for type checks using Base = Parent; //!< for uniformity of interface //! Plain Eigen type to map using EigenRep_t = Eigen::Array; - //! map returned when iterating over field + using EigenVecRep_t = Eigen::Matrix; + //! map returned when accessing entire field using EigenMap_t = Eigen::Map; + //! map returned when accessing entire const field + using EigenMapConst_t = Eigen::Map; //! Plain eigen vector to map - using EigenVec_t = Eigen::Map; - //! vector map returned when iterating over field - using EigenVecConst_t = Eigen::Map; + using EigenVec_t = Eigen::Map; + //! vector map returned when accessing entire field + using EigenVecConst_t = Eigen::Map; /** * type stored (unfortunately, we can't statically size the second - * dimension due to an Eigen bug,i.e., creating a row vector + * dimension due to an Eigen bug, i.e., creating a row vector * reference to a column vector does not raise an error :( */ using Stored_t = Eigen::Array; //! storage container using Storage_t = std::vector>; //! Default constructor TypedField() = delete; //! constructor TypedField(std::string unique_name, FieldCollection& collection, size_t nb_components); /** * constructor for field proxies which piggy-back on existing * memory. These cannot be registered in field collections and * should only be used for transient temporaries */ TypedField(std::string unique_name, FieldCollection& collection, Eigen::Ref> vec, size_t nb_components); //! Copy constructor TypedField(const TypedField &other) = delete; //! Move constructor TypedField(TypedField &&other) = default; //! Destructor virtual ~TypedField() = default; //! Copy assignment operator TypedField& operator=(const TypedField &other) = delete; //! Move assignment operator TypedField& operator=(TypedField &&other) = delete; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const override final; virtual size_t size() const override final; //! add a pad region to the end of the field buffer; required for //! using this as e.g. an FFT workspace void set_pad_size(size_t pad_size_) override final; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) virtual void set_zero() override final; //! add a new value at the end of the field inline void push_back(const Stored_t & value); //! raw pointer to content (e.g., for Eigen maps) virtual T* data() {return this->get_ptr_to_entry(0);} //! raw pointer to content (e.g., for Eigen maps) virtual const T* data() const {return this->get_ptr_to_entry(0);} //! return a map representing the entire field as a single `Eigen::Array` EigenMap_t eigen(); + //! return a map representing the entire field as a single `Eigen::Array` + EigenMapConst_t eigen() const ; //! return a map representing the entire field as a single Eigen vector EigenVec_t eigenvec(); //! return a map representing the entire field as a single Eigen vector EigenVecConst_t eigenvec() const; //! return a map representing the entire field as a single Eigen vector EigenVecConst_t const_eigenvec() const; + + protected: //! returns a raw pointer to the entry, for `Eigen::Map` inline T* get_ptr_to_entry(const size_t&& index); //! returns a raw pointer to the entry, for `Eigen::Map` inline const T* get_ptr_to_entry(const size_t&& index) const; //! set the storage size of this field inline virtual void resize(size_t size) override final; //! The actual storage container Storage_t values{}; /** * an unregistered typed field can be mapped onto an array of * existing values */ optional>> alt_values{}; /** * maintains a tally of the current size, as it cannot be reliably * determined from either `values` or `alt_values` alone. */ size_t current_size; /** * in order to accomodate both registered fields (who own and * manage their data) and unregistered temporary field proxies * (piggy-backing on a chunk of existing memory as e.g., a numpy * array) *efficiently*, the `get_ptr_to_entry` methods need to be * branchless. this means that we cannot decide on the fly whether * to return pointers pointing into values or into alt_values, we * need to maintain an (shudder) raw data pointer that is set * either at construction (for unregistered fields) or at any * resize event (which may invalidate existing pointers). For the * coder, this means that they need to be absolutely vigilant that * *any* operation on the values vector that invalidates iterators * needs to be followed by an update of data_ptr, or we will get * super annoying memory bugs. */ T* data_ptr{}; private: }; /* ---------------------------------------------------------------------- */ /* Implementations */ /* ---------------------------------------------------------------------- */ template TypedField:: TypedField(std::string unique_name, FieldCollection & collection, size_t nb_components): Parent(unique_name, nb_components, collection), current_size{0} {} /* ---------------------------------------------------------------------- */ template TypedField:: TypedField(std::string unique_name, FieldCollection & collection, Eigen::Ref> vec, size_t nb_components): Parent(unique_name, nb_components, collection), alt_values{vec}, current_size{vec.size()/nb_components}, data_ptr{vec.data()} { if (vec.size()%nb_components) { std::stringstream err{}; err << "The vector you supplied has a size of " << vec.size() << ", which is not a multiple of the number of components (" << nb_components << ")"; throw FieldError(err.str()); } if (current_size != collection.size()) { std::stringstream err{}; err << "The vector you supplied has the size for " << current_size << " pixels with " << nb_components << "components each, but the " << "field collection has " << collection.size() << " pixels."; throw FieldError(err.str()); } } /* ---------------------------------------------------------------------- */ //! return type_id of stored type template const std::type_info & TypedField:: get_stored_typeid() const { return typeid(T); } /* ---------------------------------------------------------------------- */ template auto TypedField::eigen() -> EigenMap_t { return EigenMap_t(this->data(), this->get_nb_components(), this->size()); } + /* ---------------------------------------------------------------------- */ + template + auto TypedField::eigen() const -> EigenMapConst_t { + return EigenMapConst_t(this->data(), this->get_nb_components(), this->size()); + } + /* ---------------------------------------------------------------------- */ template auto TypedField::eigenvec() -> EigenVec_t { - return EigenVec_t(this->data(), this->get_nb_components() * this->size()); + return EigenVec_t(this->data(), + this->get_nb_components() * this->size(), + 1); } /* ---------------------------------------------------------------------- */ template auto TypedField:: eigenvec() const -> EigenVecConst_t { - return EigenVecConst_t(this->data(), this->get_nb_components() * this->size()); + return EigenVecConst_t(this->data(), + this->get_nb_components() * this->size(), + 1); } /* ---------------------------------------------------------------------- */ template auto TypedField:: const_eigenvec() const -> EigenVecConst_t { return EigenVecConst_t(this->data(), this->get_nb_components() * this->size()); } /* ---------------------------------------------------------------------- */ template void TypedField::resize(size_t size) { if (this->alt_values) { throw FieldError("Field proxies can't resize."); } this->current_size = size; this->values.resize(size*this->get_nb_components() + this->pad_size); this->data_ptr = &this->values.front(); } /* ---------------------------------------------------------------------- */ template void TypedField::set_zero() { std::fill(this->values.begin(), this->values.end(), T{}); } /* ---------------------------------------------------------------------- */ template size_t TypedField:: size() const { return this->current_size; } /* ---------------------------------------------------------------------- */ template void TypedField:: set_pad_size(size_t pad_size) { if (this->alt_values) { throw FieldError("You can't set the pad size of a field proxy."); } this->pad_size = pad_size; this->resize(this->size()); this->data_ptr = &this->values.front(); } /* ---------------------------------------------------------------------- */ template T* TypedField:: get_ptr_to_entry(const size_t&& index) { return this->data_ptr + this->get_nb_components()*std::move(index); } /* ---------------------------------------------------------------------- */ template const T* TypedField:: get_ptr_to_entry(const size_t&& index) const { return this->data_ptr+this->get_nb_components()*std::move(index); } /* ---------------------------------------------------------------------- */ template void TypedField:: push_back(const Stored_t & value) { static_assert (not FieldCollection::Global, "You can only push_back data into local field collections"); if (value.cols() != 1) { std::stringstream err{}; err << "Expected a column vector, but received and array with " << value.cols() <<" colums."; throw FieldError(err.str()); } if (value.rows() != static_cast(this->get_nb_components())) { std::stringstream err{}; err << "Expected a column vector of length " << this->get_nb_components() << ", but received one of length " << value.rows() <<"."; throw FieldError(err.str()); } for (size_t i = 0; i < this->get_nb_components(); ++i) { this->values.push_back(value(i)); } ++this->current_size; this->data_ptr = &this->values.front(); } } // muSpectre #endif /* FIELD_TYPED_H */ diff --git a/src/common/statefield.hh b/src/common/statefield.hh index 19b11ff..460e272 100644 --- a/src/common/statefield.hh +++ b/src/common/statefield.hh @@ -1,546 +1,661 @@ /** * file statefield.hh * * @author Till Junge * * @date 28 Feb 2018 * * @brief A state field is an abstraction of a field that can hold * current, as well as a chosen number of previous values. This is * useful for instance for internal state variables in plastic laws, * where a current, new, or trial state is computed based on its * previous state, and at convergence, this new state gets cycled into * the old, the old into the old-1 etc. The state field abstraction * helps doing this safely (i.e. only const references to the old * states are available, while the current state can be assigned * to/modified), and efficiently (i.e., no need to copy values from * new to old, we just cycle the labels). This file implements the * state field as well as state maps using the Field, FieldCollection * and FieldMap abstractions of µSpectre * * 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. */ #ifndef STATEFIELD_H #define STATEFIELD_H #include "common/field.hh" #include "common/utilities.hh" #include #include #include namespace muSpectre { /** * Base class for state fields, useful for storing polymorphic references */ - template + template class StateFieldBase { public: + //! get naming prefix + const std::string & get_prefix() const {return this->prefix;} + + //! get a ref to the `StateField` 's field collection + const FieldCollection & get_collection() const { + return this->collection;} + + virtual ~StateFieldBase() = default; + + /** + * returns number of old states that are stored + */ + size_t get_nb_memory() const {return this->nb_memory;} + + //! return type_id of stored type + virtual const std::type_info & get_stored_typeid() const = 0; + + /** + * cycle the fields (current becomes old, old becomes older, + * oldest becomes current) + */ + virtual void cycle() = 0; + protected: + //! constructor + StateFieldBase(std::string unique_prefix, + const FieldCollection & collection, + size_t nb_memory=1): + prefix{unique_prefix}, + nb_memory{nb_memory}, + collection{collection} {} + + /** + * the unique prefix is used as the first part of the unique name + * of the subfields belonging to this state field + */ + std::string prefix; + /** + * number of old states to store, defaults to 1 + */ + const size_t nb_memory; + //! reference to the collection this statefield belongs to + const FieldCollection & collection; + }; + + /* ---------------------------------------------------------------------- */ + template + class TypedStateField: + public StateFieldBase { + public: + //! Parent class + using Parent = StateFieldBase; + //! Typed field + using TypedField_t = TypedField; + + //! returns a TypedField ref to the current value of this state field + virtual TypedField_t & get_current_field() = 0; + + //! returns a const TypedField ref to an old value of this state field + virtual const TypedField_t & get_old_field(size_t nb_steps_ago=1) const = 0; + + //! return type_id of stored type + const std::type_info & get_stored_typeid() const override final { + return typeid(T); + }; + + virtual ~TypedStateField() = default; + + protected: + //! constructor + TypedStateField(const std::string & unique_prefix, + const FieldCollection & collection, + size_t nb_memory): + Parent{unique_prefix, collection, nb_memory} + {} + }; + + + + /* ---------------------------------------------------------------------- */ + template + class TypedSizedStateField: public TypedStateField { + public: + //! Parent class + using Parent = TypedStateField; //! the current (historically accurate) ordering of the fields using index_t = std::array; //! get the current ordering of the fields inline const index_t & get_indices() const {return this->indices;} + //! destructor + virtual ~TypedSizedStateField() = default; + protected: //! constructor - StateFieldBase(index_t indices):indices{indices}{}; + TypedSizedStateField(std::string unique_prefix, + const FieldCollection& collection, + index_t indices): + Parent{unique_prefix, collection, nb_memory}, + indices{indices}{}; index_t indices; ///< these are cycled through }; //! early declaration template class StateFieldMap; namespace internal { template inline decltype(auto) build_fields_helper(std::string prefix, typename Field::Base::collection_t & collection, std::index_sequence) { auto get_field{[&prefix, &collection](size_t i) -> Field& { std::stringstream name_stream{}; name_stream << prefix << ", sub_field index " << i; return make_field(name_stream.str(), collection); }}; return std::tie(get_field(I)...); } /* ---------------------------------------------------------------------- */ template - decltype(auto) build_indices(std::index_sequence) { - return std::array{I...}; + inline decltype(auto) build_indices(std::index_sequence) { + return std::array{(size-I)%size...}; } } // internal /** * A statefield is an abstraction around a Field that can hold a * current and `nb_memory` previous values. There are useful for * history variables, for instance. */ template - class StateField: public StateFieldBase + class StateField: + public TypedSizedStateField { public: - //! Base class for all state fields of same memory - using Base = StateFieldBase; //! the underlying field's collection type - using FieldCollection = typename Field_t::Base::collection_t; + using FieldCollection_t = typename Field_t::Base::collection_t; + //! base type for fields + using Scalar = typename Field_t::Scalar; + //! Base class for all state fields of same memory + using Base = TypedSizedStateField; /** * storage of field refs (can't be a `std::array`, because arrays * of refs are explicitely forbidden */ using Fields_t = tuple_array; + //! Typed field + using TypedField_t = TypedField; //! Default constructor StateField() = delete; - /** - * Constructor. @param unique_prefix is used to create the names - * of the fields that this abstraction creates in the background - * @param collection is the field collection in which the - * subfields will be stored - */ - inline StateField(std::string unique_prefix, FieldCollection & collection) - : Base{internal::build_indices(std::make_index_sequence{})}, - prefix{unique_prefix}, fields{internal::build_fields_helper - (unique_prefix, collection, std::make_index_sequence{})} - {} - //! Copy constructor StateField(const StateField &other) = delete; //! Move constructor StateField(StateField &&other) = delete; //! Destructor virtual ~StateField() = default; //! Copy assignment operator StateField& operator=(const StateField &other) = delete; //! Move assignment operator StateField& operator=(StateField &&other) = delete; //! get (modifiable) current field - inline StateField& current() { + inline Field_t& current() { return this->fields[this->indices[0]]; } //! get (constant) previous field template - inline const StateField& old() { + inline const Field_t& old() { static_assert(nb_steps_ago <= nb_memory, "you can't go that far inte the past"); static_assert(nb_steps_ago > 0, "Did you mean to call current()?"); - return this->fields[this->indices[nb_steps_ago]]; + return this->fields[this->indices.at(nb_steps_ago)]; + } + + //! returns a TypedField ref to the current value of this state field + TypedField_t & get_current_field() override final { + return this->current(); + } + + //! returns a const TypedField ref to an old value of this state field + const TypedField_t & + get_old_field(size_t nb_steps_ago=1) const override final { + return this->fields[this->indices.at(nb_steps_ago)]; } - // //! factory function - // template - // friend StateFieldType& make_statefield(std::string unique_prefix, - // CollectionType & collection, - // Args&&... args); + //! factory function + template + friend StateFieldType& make_statefield(const std::string & unique_prefix, + CollectionType & collection); //! returns a `StateField` reference if `other is a compatible state field inline static StateField& check_ref(Base& other) { // the following triggers and exception if the fields are incompatible Field_t::check_ref(other.fields[0]); return static_cast (other); } //! returns a const `StateField` reference if `other` is a compatible state field inline static const StateField& check_ref(const Base& other) { // the following triggers and exception if the fields are incompatible Field_t::check_ref(other.fields[0]); return static_cast (other); } - //! get naming prefix - const std::string & get_prefix() const {return this->prefix;} - - //! get a ref to the `StateField` 's field collection - const FieldCollection & get_collection() const { - return std::get<0>(this->fields).get_collection();} - //! get a ref to the `StateField` 's fields Fields_t & get_fields() { return this->fields; } /** * Pure convenience functions to get a MatrixFieldMap of * appropriate dimensions mapped to this field. You can also * create other types of maps, as long as they have the right * fundamental type (T), the correct size (nbComponents), and * memory (nb_memory). */ inline decltype(auto) get_map() { using FieldMap = decltype(std::get<0>(this->fields).get_map()); return StateFieldMap(*this); } /** * Pure convenience functions to get a MatrixFieldMap of * appropriate dimensions mapped to this field. You can also * create other types of maps, as long as they have the right * fundamental type (T), the correct size (nbComponents), and * memory (nb_memory). */ inline decltype(auto) get_const_map() { using FieldMap = decltype(std::get<0>(this->fields).get_const_map()); return StateFieldMap(*this); } /** * cycle the fields (current becomes old, old becomes older, * oldest becomes current) */ - inline void cycle() { + inline void cycle() override final { for (auto & val: this->indices) { val = (val+1)%(nb_memory+1); } } protected: /** - * the unique prefix is used as the first part of the unique name - * of the subfields belonging to this state field + * Constructor. @param unique_prefix is used to create the names + * of the fields that this abstraction creates in the background + * @param collection is the field collection in which the + * subfields will be stored */ - std::string prefix; + inline StateField(const std::string & unique_prefix, + FieldCollection_t & collection) + : Base{unique_prefix, collection, + internal::build_indices + (std::make_index_sequence{})}, + fields{internal::build_fields_helper + (unique_prefix, collection, std::make_index_sequence{})} + {} + Fields_t fields; //!< container for the states private: }; namespace internal { template inline decltype(auto) build_maps_helper(Fields & fields, std::index_sequence) { return std::array{FieldMap(std::get(fields))...}; } } // internal + /* ---------------------------------------------------------------------- */ + template + inline StateFieldType & + make_statefield(const std::string & unique_prefix, + CollectionType & collection) { + std::unique_ptr ptr { + new StateFieldType(unique_prefix, collection)}; + auto & retref{*ptr}; + collection.register_statefield(std::move(ptr)); + return retref; + } + /** * extends the StateField <-> Field equivalence to StateFieldMap <-> FieldMap */ template class StateFieldMap { public: /** * iterates over all pixels in the `muSpectre::FieldCollection` and * dereferences to a proxy giving access to the appropriate iterates * of the underlying `FieldMap` type. */ class iterator; //! stl conformance using reference = typename iterator::reference; //! stl conformance using value_type = typename iterator::value_type; //! stl conformance using size_type = typename iterator::size_type; //! field collection type where this state field can be stored - using FieldCollection= typename FieldMap::Field::collection_t; - //! base class (for polymorphic references) - using StateFieldBase_t = StateFieldBase; + using FieldCollection_t= typename FieldMap::Field::collection_t; + + //! Fundamental type stored + using Scalar = typename FieldMap::Scalar; + //! base class (must be at least sized) + using TypedSizedStateField_t = TypedSizedStateField; //! for traits access using FieldMap_t = FieldMap; //! for traits access using ConstFieldMap_t = typename FieldMap::ConstMap; //! Default constructor StateFieldMap() = delete; //! constructor using a StateField template StateFieldMap(StateField & statefield) :collection{statefield.get_collection()}, statefield{statefield}, maps{internal::build_maps_helper (statefield.get_fields(), std::make_index_sequence{})}, const_maps{internal::build_maps_helper (statefield.get_fields(), std::make_index_sequence{})} { - static_assert(std::is_base_of::value, + static_assert(std::is_base_of::value, "Not the right type of StateField ref"); } //! Copy constructor StateFieldMap(const StateFieldMap &other) = delete; //! Move constructor StateFieldMap(StateFieldMap &&other) = default; //! Destructor virtual ~StateFieldMap() = default; //! Copy assignment operator StateFieldMap& operator=(const StateFieldMap &other) = delete; //! Move assignment operator StateFieldMap& operator=(StateFieldMap &&other) = delete; //! access the wrapper to a given pixel directly value_type operator[](size_type index) { return *iterator(*this, index); } /** * return a ref to the current field map. useful for instance for * initialisations of `StateField` instances */ FieldMap& current() { return this->maps[this->statefield.get_indices()[0]]; } //! stl conformance iterator begin() { return iterator(*this, 0);} //! stl conformance iterator end() { return iterator(*this, this->collection.size());} protected: - const FieldCollection & collection; //!< collection holding the field - StateFieldBase_t & statefield; //!< ref to the field itself + const FieldCollection_t & collection; //!< collection holding the field + TypedSizedStateField_t & statefield; //!< ref to the field itself std::array maps;//!< refs to the addressable maps; //! const refs to the addressable maps; std::array const_maps; private: }; /** * Iterator class used by the `StateFieldMap` */ template class StateFieldMap::iterator { public: class StateWrapper; using Ccoord = typename FieldMap::Ccoord; //!< cell coordinates type using value_type = StateWrapper; //!< stl conformance using const_value_type = value_type; //!< stl conformance using pointer_type = value_type*; //!< stl conformance using difference_type = std::ptrdiff_t; //!< stl conformance using size_type = size_t; //!< stl conformance using iterator_category = std::random_access_iterator_tag; //!< stl conformance using reference = StateWrapper; //!< stl conformance //! Default constructor iterator() = delete; //! constructor iterator(StateFieldMap& map, size_t index = 0) :index{index}, map{map} {}; //! Copy constructor iterator(const iterator &other) = default; //! Move constructor iterator(iterator &&other) = default; //! Destructor virtual ~iterator() = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) = default; //! pre-increment inline iterator & operator++() { this->index++; return *this;} //! post-increment inline iterator operator++(int) { iterator curr{*this}; this->index++; return curr;} //! dereference inline value_type operator*() { return value_type(*this);} //! pre-decrement inline iterator & operator--() { this->index--; return *this;} //! post-decrement inline iterator operator--(int) { iterator curr{*this}; this->index--; return curr;} //! access subscripting inline value_type operator[](difference_type diff) { return value_type{iterator{this->map, this->index+diff}};} //! equality inline bool operator==(const iterator & other) const { return this->index == other.index; } //! inequality inline bool operator!=(const iterator & other) const { return this->index != other.index;} //! div. comparisons inline bool operator<(const iterator & other) const { return this->index < other.index; } //! div. comparisons inline bool operator<=(const iterator & other) const { return this->index <= other.index; } //! div. comparisons inline bool operator>(const iterator & other) const { return this->index > other.index; } //! div. comparisons inline bool operator>=(const iterator & other) const { return this->index >= other.index; } //! additions, subtractions and corresponding assignments inline iterator operator+(difference_type diff) const { return iterator{this->map, this-index + diff}; } //! additions, subtractions and corresponding assignments inline iterator operator-(difference_type diff) const { return iterator{this->map, this-index - diff};} //! additions, subtractions and corresponding assignments inline iterator& operator+=(difference_type diff) { this->index += diff; return *this;} //! additions, subtractions and corresponding assignments inline iterator& operator-=(difference_type diff) { this->index -= diff; return *this; } //! get pixel coordinates inline Ccoord get_ccoord() const { return this->map.collection.get_ccoord(this->index); } //! access the index inline const size_t & get_index() const {return this->index;} protected: size_t index; //!< current pixel this iterator refers to StateFieldMap& map; //!< map over with `this` iterates private: }; namespace internal { //! FieldMap is an `Eigen::Map` or `Eigen::TensorMap` here template inline decltype(auto) build_old_vals_helper(iterator& it, maps_t & maps, indices_t & indices, std::index_sequence) { return tuple_array(std::forward_as_tuple(maps[indices[I+1]][it.get_index()]...)); } template inline decltype(auto) build_old_vals(iterator& it, maps_t & maps, indices_t & indices) { return tuple_array{build_old_vals_helper (it, maps, indices, std::make_index_sequence{})}; } } // internal /** * Light-weight resource-handle representing the current and old * values of a field at a given pixel identified by an iterator * pointing to it */ template class StateFieldMap::iterator::StateWrapper { public: //! short-hand using iterator = typename StateFieldMap::iterator; //! short-hand using Ccoord = typename iterator::Ccoord; //! short-hand using Map = typename FieldMap::reference; //! short-hand using ConstMap = typename FieldMap::const_reference; //! Default constructor StateWrapper() = delete; //! Copy constructor StateWrapper(const StateWrapper &other) = default; //! Move constructor StateWrapper(StateWrapper &&other) = default; //! construct with `StateFieldMap::iterator` StateWrapper(iterator & it) :it{it}, current_val{it.map.maps[it.map.statefield.get_indices()[0]][it.index]}, old_vals(internal::build_old_vals (it, it.map.const_maps, it.map.statefield.get_indices())) { } //! Destructor virtual ~StateWrapper() = default; //! Copy assignment operator StateWrapper& operator=(const StateWrapper &other) = default; //! Move assignment operator StateWrapper& operator=(StateWrapper &&other) = default; //! returns reference to the currectly mapped value inline Map& current() { return this->current_val; } //! recurnts reference the the value that was current `nb_steps_ago` ago template inline const ConstMap & old() const{ static_assert (nb_steps_ago <= nb_memory, "You have not stored that time step"); static_assert (nb_steps_ago > 0, "Did you mean to access the current value? If so, use " "current()"); - return std::get(this->old_vals); + return std::get(this->old_vals); } //! read the coordinates of the current pixel inline Ccoord get_ccoord() const { return this->it.get_ccoord(); } protected: iterator& it; //!< ref to the iterator that dereferences to `this` Map current_val; //!< current value tuple_array old_vals; //!< all stored old values private: }; } // muSpectre #endif /* STATEFIELD_H */ diff --git a/src/common/utilities.hh b/src/common/utilities.hh index b95c29c..4fabaf5 100644 --- a/src/common/utilities.hh +++ b/src/common/utilities.hh @@ -1,268 +1,338 @@ /** * @file utilities.hh * * @author Till Junge * * @date 17 Nov 2017 * * @brief additions to the standard name space to anticipate C++17 features * * Copyright © 2017 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. */ #ifndef UTILITIES_H #define UTILITIES_H #include #include #ifdef NO_EXPERIMENTAL # if defined(__INTEL_COMPILER) //# pragma warning ( disable : 383 ) # elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu # pragma clang diagnostic push # pragma clang diagnostic ignored "-Weffc++" # elif (defined(__GNUC__) || defined(__GNUG__)) # pragma GCC diagnostic push # pragma GCC diagnostic ignored "-Weffc++" # endif # include # if defined(__INTEL_COMPILER) //# pragma warning ( disable : 383 ) # elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu # pragma clang diagnostic pop # pragma clang diagnostic ignored "-Weffc++" # elif (defined(__GNUC__) || defined(__GNUG__)) # pragma GCC diagnostic pop # pragma GCC diagnostic ignored "-Weffc++" # endif #else # include #endif namespace std_replacement { namespace detail { template struct is_reference_wrapper : std::false_type {}; template struct is_reference_wrapper> : std::true_type {}; //! from cppreference template auto INVOKE(T Base::*pmf, Derived&& ref, Args&&... args) noexcept(noexcept((std::forward(ref).*pmf)(std::forward(args)...))) -> std::enable_if_t::value && std::is_base_of>::value, decltype((std::forward(ref).*pmf)(std::forward(args)...))> { return (std::forward(ref).*pmf)(std::forward(args)...); } //! from cppreference template auto INVOKE(T Base::*pmf, RefWrap&& ref, Args&&... args) noexcept(noexcept((ref.get().*pmf)(std::forward(args)...))) -> std::enable_if_t::value && is_reference_wrapper>::value, decltype((ref.get().*pmf)(std::forward(args)...))> { return (ref.get().*pmf)(std::forward(args)...); } //! from cppreference template auto INVOKE(T Base::*pmf, Pointer&& ptr, Args&&... args) noexcept(noexcept(((*std::forward(ptr)).*pmf)(std::forward(args)...))) -> std::enable_if_t::value && !is_reference_wrapper>::value && !std::is_base_of>::value, decltype(((*std::forward(ptr)).*pmf)(std::forward(args)...))> { return ((*std::forward(ptr)).*pmf)(std::forward(args)...); } //! from cppreference template auto INVOKE(T Base::*pmd, Derived&& ref) noexcept(noexcept(std::forward(ref).*pmd)) -> std::enable_if_t::value && std::is_base_of>::value, decltype(std::forward(ref).*pmd)> { return std::forward(ref).*pmd; } //! from cppreference template auto INVOKE(T Base::*pmd, RefWrap&& ref) noexcept(noexcept(ref.get().*pmd)) -> std::enable_if_t::value && is_reference_wrapper>::value, decltype(ref.get().*pmd)> { return ref.get().*pmd; } //! from cppreference template auto INVOKE(T Base::*pmd, Pointer&& ptr) noexcept(noexcept((*std::forward(ptr)).*pmd)) -> std::enable_if_t::value && !is_reference_wrapper>::value && !std::is_base_of>::value, decltype((*std::forward(ptr)).*pmd)> { return (*std::forward(ptr)).*pmd; } //! from cppreference template auto INVOKE(F&& f, Args&&... args) noexcept(noexcept(std::forward(f)(std::forward(args)...))) -> std::enable_if_t>::value, decltype(std::forward(f)(std::forward(args)...))> { return std::forward(f)(std::forward(args)...); } } // namespace detail //! from cppreference template< class F, class... ArgTypes > auto invoke(F&& f, ArgTypes&&... args) // exception specification for QoI noexcept(noexcept(detail::INVOKE(std::forward(f), std::forward(args)...))) -> decltype(detail::INVOKE(std::forward(f), std::forward(args)...)) { return detail::INVOKE(std::forward(f), std::forward(args)...); } namespace detail { //! from cppreference template constexpr decltype(auto) apply_impl(F &&f, Tuple &&t, std::index_sequence) { return std_replacement::invoke(std::forward(f), std::get(std::forward(t))...); } } // namespace detail //! from cppreference template constexpr decltype(auto) apply(F &&f, Tuple &&t) { return detail::apply_impl (std::forward(f), std::forward(t), std::make_index_sequence>::value>{}); } } //namespace std_replacement namespace muSpectre { namespace internal { /** * helper struct template to compute the type of a tuple with a * given number of entries of the same type */ template struct tuple_array_helper { //! underlying tuple using type = typename tuple_array_helper::type; }; /** * helper struct template to compute the type of a tuple with a * given number of entries of the same type */ template< typename T, typename... tail> struct tuple_array_helper<0, T, tail...> { //! underlying tuple using type = std::tuple; }; + /** + * helper struct for runtime index access to + * tuples. RecursionLevel indicates how much more we can recurse + * down + */ + template + struct Accessor { + using Stored_t = typename TupArr::Stored_t; + + inline static Stored_t + get(const size_t & index, TupArr & container) { + if (index == Index) { + return std::get(container); + } else { + return Accessor::get(index, container); + } + } + inline static const Stored_t + get(const size_t & index, const TupArr & container) { + if (index == Index) { + return std::get(container); + } else { + return Accessor::get(index, container); + } + } + }; + + /** + * specialisation for recursion end + */ + template + struct Accessor { + using Stored_t = typename TupArr::Stored_t; + + inline static Stored_t + get(const size_t & index, TupArr & container) { + if (index == Index) { + return std::get(container); + } else { + std::stringstream err{}; + err << "Index " << index << "is out of range."; + throw std::runtime_error(err.str()); + } + } + + inline static const Stored_t + get(const size_t & index, const TupArr & container) { + if (index == Index) { + return std::get(container); + } else { + std::stringstream err{}; + err << "Index " << index << "is out of range."; + throw std::runtime_error(err.str()); + } + } + }; + /** * helper struct that provides the tuple_array. */ template struct tuple_array_provider { //! tuple type that can be used (almost) like an `std::array` class type: public tuple_array_helper::type { public: //! short-hand using Parent = typename tuple_array_helper::type; + using Stored_t = T; + constexpr static size_t Size{size}; //! constructor inline type(Parent && parent):Parent{parent}{}; + + //! element access + T operator[] (const size_t & index) { + return Accessor::get(index, *this); + } + + //! element access + const T operator[](const size_t & index) const { + return Accessor::get(index, *this); + } + protected: }; }; } // internal /** * This is a convenience structure to create a tuple of `nb_elem` * entries of type `T`. It is named tuple_array, because it is * somewhat similar to an `std::array`. The reason for * this structure is that the `std::array` is not allowed by the * standard to store references (8.3.2 References, paragraph 5: * "There shall be no references to references, no arrays of * references, and no pointers to references.") use this, if you * want to have a statically known number of references to store, * and you wish to do so efficiently. */ template using tuple_array = typename internal::tuple_array_provider::type; using std_replacement::apply; /** * emulation `std::optional` (a C++17 feature) */ template #ifdef NO_EXPERIMENTAL using optional = typename boost::optional; #else using optional = typename std::experimental::optional; #endif /* ---------------------------------------------------------------------- */ /** * conversion helper from `boost::tuple` to `std::tuple` */ template auto asStdTuple(BoostTuple&& boostTuple, std::index_sequence) { return std::tuple>::type...> (boost::get(std::forward(boostTuple))...); } /** * conversion from `boost::tuple` to `std::tuple` */ template auto asStdTuple(BoostTuple&& boostTuple) { return asStdTuple(std::forward(boostTuple), std::make_index_sequence>::value>()); } } // muSpectre #endif /* UTILITIES_H */ diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh index 9d1365c..b05865f 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,161 +1,164 @@ /** * @file material_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials (constitutive models) * * Copyright © 2017 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. */ #ifndef MATERIAL_BASE_H #define MATERIAL_BASE_H #include "common/common.hh" #include "common/field.hh" #include "common/field_collection.hh" #include namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) and /** * @a DimM is the material dimension (i.e., the dimension of constitutive * law; even for e.g. two-dimensional problems the constitutive law could * live in three-dimensional space for e.g. plane strain or stress problems) */ template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for cell-wide fields, like stress, strain, etc using GFieldCollection_t = GlobalFieldCollection; //! field collection for internal variables, such as eigen-strains, //! plastic strains, damage variables, etc, but also for managing which //! pixels the material is responsible for using MFieldCollection_t = LocalFieldCollection; using iterator = typename MFieldCollection_t::iterator; //!< pixel iterator //! polymorphic base class for fields only to be used for debugging using Field_t = internal::FieldBase; //! Full type for stress fields using StressField_t = TensorField; //! Full type for strain fields using StrainField_t = StressField_t; //! Full type for tangent stiffness fields fields using TangentField_t = TensorField; using Ccoord = Ccoord_t; //!< cell coordinates type //! Default constructor MaterialBase() = delete; //! Construct by name MaterialBase(std::string name); //! Copy constructor MaterialBase(const MaterialBase &other) = delete; //! Move constructor MaterialBase(MaterialBase &&other) = delete; //! Destructor virtual ~MaterialBase() = default; //! Copy assignment operator MaterialBase& operator=(const MaterialBase &other) = delete; //! Move assignment operator MaterialBase& operator=(MaterialBase &&other) = delete; /** * take responsibility for a pixel identified by its cell coordinates * WARNING: this won't work for materials with additional info per pixel * (as, e.g. for eigenstrain), we need to pass more parameters. Materials * of this tye need to overload add_pixel */ virtual void add_pixel(const Ccoord & ccooord); //! allocate memory, etc, but also: wipe history variables! virtual void initialise() = 0; /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, the virtual * base implementation does nothing, but materials with history * variables need to implement this */ virtual void save_history_variables() {}; //! return the material's name const std::string & get_name() const; //! spatial dimension for static inheritance constexpr static Dim_t sdim() {return DimS;} //! material dimension for static inheritance constexpr static Dim_t mdim() {return DimM;} //! computes stress virtual void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form) = 0; /** * Convenience function to compute stresses, mostly for debugging and * testing. Has runtime-cost associated with compatibility-checking and * conversion of the Field_t arguments that can be avoided by using the * version with strongly typed field references */ void compute_stresses(const Field_t & F, Field_t & P, Formulation form); //! computes stress and tangent moduli virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) = 0; /** * Convenience function to compute stresses and tangent moduli, mostly for * debugging and testing. Has runtime-cost associated with * compatibility-checking and conversion of the Field_t arguments that can * be avoided by using the version with strongly typed field references */ void compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, Formulation form); //! iterator to first pixel handled by this material inline iterator begin() {return this->internal_fields.begin();} //! iterator past the last pixel handled by this material inline iterator end() {return this->internal_fields.end();} //! number of pixels assigned to this material inline size_t size() const {return this->internal_fields.size();} + + //! gives access to internal fields + inline MFieldCollection_t& get_collection() {return this->internal_fields;} protected: const std::string name; //!< material's name (for output and debugging) MFieldCollection_t internal_fields{};//!< storage for internal variables private: }; } // muSpectre #endif /* MATERIAL_BASE_H */ diff --git a/src/materials/material_hyper_elasto_plastic1.cc b/src/materials/material_hyper_elasto_plastic1.cc index 9f4ecb7..0a7ec51 100644 --- a/src/materials/material_hyper_elasto_plastic1.cc +++ b/src/materials/material_hyper_elasto_plastic1.cc @@ -1,77 +1,82 @@ /** * @file material_hyper_elasto_plastic1.cc * * @author Till Junge * * @date 21 Feb 2018 * * @brief implementation for MaterialHyperElastoPlastic1 * * 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 "materials/material_hyper_elasto_plastic1.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialHyperElastoPlastic1:: MaterialHyperElastoPlastic1(std::string name, Real young, Real poisson, Real tau_y0, Real H) : Parent{name}, - plast_flow_field("cumulated plastic flow εₚ", this->internal_fields), - F_prev_field("Previous placement gradient Fᵗ", this->internal_fields), - be_prev_field("Previous left Cauchy-Green deformation bₑᵗ", - this->internal_fields), + plast_flow_field{make_statefield>> + ("cumulated plastic flow εₚ", this->internal_fields)}, + F_prev_field{make_statefield>> + ("Previous placement gradient Fᵗ", this->internal_fields)}, + be_prev_field{make_statefield>> + ("Previous left Cauchy-Green deformation bₑᵗ", + this->internal_fields)}, young{young}, poisson{poisson}, lambda{Hooke::compute_lambda(young, poisson)}, mu{Hooke::compute_mu(young, poisson)}, K{Hooke::compute_K(young, poisson)}, tau_y0{tau_y0}, H{H}, // the factor .5 comes from equation (18) in Geers 2003 // (https://doi.org/10.1016/j.cma.2003.07.014) C{0.5*Hooke::compute_C_T4(lambda, mu)}, internal_variables{F_prev_field.get_map(), be_prev_field.get_map(), plast_flow_field.get_map()} {} /* ---------------------------------------------------------------------- */ template void MaterialHyperElastoPlastic1::save_history_variables() { this->plast_flow_field.cycle(); this->F_prev_field.cycle(); this->be_prev_field.cycle(); } /* ---------------------------------------------------------------------- */ template void MaterialHyperElastoPlastic1::initialise() { Parent::initialise(); this->F_prev_field.get_map().current() = Eigen::Matrix::Identity(); this->be_prev_field.get_map().current() = Eigen::Matrix::Identity(); this->save_history_variables(); } template class MaterialHyperElastoPlastic1< twoD, twoD>; template class MaterialHyperElastoPlastic1< twoD, threeD>; template class MaterialHyperElastoPlastic1; } // muSpectre diff --git a/src/materials/material_hyper_elasto_plastic1.hh b/src/materials/material_hyper_elasto_plastic1.hh index a7ce609..8a6ce26 100644 --- a/src/materials/material_hyper_elasto_plastic1.hh +++ b/src/materials/material_hyper_elasto_plastic1.hh @@ -1,307 +1,308 @@ /** * @file material_hyper_elasto_plastic1.hh * * @author Till Junge * * @date 20 Feb 2018 * * @brief Material for logarithmic hyperelasto-plasticity, as defined in de * Geus 2017 (https://doi.org/10.1016/j.cma.2016.12.032) and further * explained in Geers 2003 (https://doi.org/10.1016/j.cma.2003.07.014) * * 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. */ #ifndef MATERIAL_HYPER_ELASTO_PLASTIC1_H #define MATERIAL_HYPER_ELASTO_PLASTIC1_H #include "materials/material_muSpectre_base.hh" #include "materials/materials_toolbox.hh" #include "common/eigen_tools.hh" #include "common/statefield.hh" #include namespace muSpectre { template class MaterialHyperElastoPlastic1; /** * traits for hyper-elastoplastic material */ template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! 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::Gradient}; //! declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::Kirchhoff}; //! local field collection used for internals using LFieldColl_t = LocalFieldCollection; //! storage type for plastic flow measure (εₚ in the papers) using LScalarMap_t = StateFieldMap>; /** * storage type for for previous gradient Fᵗ and elastic left * Cauchy-Green deformation tensor bₑᵗ */ using LStrainMap_t = StateFieldMap< MatrixFieldMap>; /** * format in which to receive internals (previous gradient Fᵗ, * previous elastic lef Cauchy-Green deformation tensor bₑᵗ, and * the plastic flow measure εₚ */ using InternalVariables = std::tuple; }; /** * Material implementation for hyper-elastoplastic constitutive law */ template class MaterialHyperElastoPlastic1: public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre , DimS, DimM>; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! shortcut to traits using traits = MaterialMuSpectre_traits; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! type in which the previous strain state is referenced using StrainStRef_t = typename traits::LStrainMap_t::reference; //! type in which the previous plastic flow is referenced using FlowStRef_t = typename traits::LScalarMap_t::reference; //! Default constructor MaterialHyperElastoPlastic1() = delete; //! Constructor with name and material properties MaterialHyperElastoPlastic1(std::string name, Real young, Real poisson, Real tau_y0, Real H); //! Copy constructor MaterialHyperElastoPlastic1(const MaterialHyperElastoPlastic1 &other) = delete; //! Move constructor MaterialHyperElastoPlastic1(MaterialHyperElastoPlastic1 &&other) = delete; //! Destructor virtual ~MaterialHyperElastoPlastic1() = default; //! Copy assignment operator MaterialHyperElastoPlastic1& operator=(const MaterialHyperElastoPlastic1 &other) = delete; //! Move assignment operator MaterialHyperElastoPlastic1& operator=(MaterialHyperElastoPlastic1 &&other) = delete; /** * evaluates Kirchhoff stress given the current placement gradient * Fₜ, the previous Gradient Fₜ₋₁ and the cumulated plastic flow * εₚ */ template inline decltype(auto) evaluate_stress(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t plast_flow); /** * evaluates Kirchhoff stress and stiffness given the current placement gradient * Fₜ, the previous Gradient Fₜ₋₁ and the cumulated plastic flow * εₚ */ template inline decltype(auto) evaluate_stress_tangent(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t plast_flow); /** * The statefields need to be cycled at the end of each load increment */ virtual void save_history_variables() override; /** * set the previous gradients to identity */ virtual void initialise() override final; /** * return the internals tuple */ typename traits::InternalVariables & get_internals() { return this->internal_variables;}; protected: /** * worker function computing stresses and internal variables */ template inline decltype(auto) stress_n_internals_worker(grad_t && F, StrainStRef_t& F_prev, StrainStRef_t& be_prev, FlowStRef_t& plast_flow); //! Local FieldCollection type for field storage using LColl_t = LocalFieldCollection; //! storage for cumulated plastic flow εₚ - StateField> plast_flow_field; + StateField> & plast_flow_field; //! storage for previous gradient Fᵗ - StateField> F_prev_field; + StateField> & F_prev_field; //! storage for elastic left Cauchy-Green deformation tensor bₑᵗ - StateField> be_prev_field; + StateField> & be_prev_field; // material properties const Real young; //!< Young's modulus const Real poisson; //!< Poisson's ratio const Real lambda; //!< first Lamé constant const Real mu; //!< second Lamé constant (shear modulus) const Real K; //!< Bulk modulus const Real tau_y0; //!< initial yield stress const Real H; //!< hardening modulus const T4Mat C; //!< stiffness tensor //! Field maps and state field maps over internal fields typename traits::InternalVariables internal_variables; private: }; //----------------------------------------------------------------------------// template template - decltype(auto) + auto MaterialHyperElastoPlastic1:: stress_n_internals_worker(grad_t && F, StrainStRef_t& F_prev, - StrainStRef_t& be_prev, FlowStRef_t& eps_p) { + StrainStRef_t& be_prev, FlowStRef_t& eps_p) + -> decltype(auto) { // the notation in this function follows Geers 2003 // (https://doi.org/10.1016/j.cma.2003.07.014). // computation of trial state using Mat_t = Eigen::Matrix; auto && f{F*F_prev.old().inverse()}; Mat_t be_star{f*be_prev.old()*f.transpose()}; Mat_t ln_be_star{logm(std::move(be_star))}; Mat_t tau_star{.5*Hooke::evaluate_stress(this->lambda, this->mu, ln_be_star)}; // deviatoric part of Kirchhoff stress Mat_t tau_d_star{tau_star - tau_star.trace()/DimM*tau_star.Identity()}; auto && tau_eq_star{std::sqrt(3*.5*(tau_d_star.array()* tau_d_star.transpose().array()).sum())}; Mat_t N_star{3*.5*tau_d_star/tau_eq_star}; // this is eq (27), and the std::min enforces the Kuhn-Tucker relation (16) Real phi_star{std::max(tau_eq_star - this->tau_y0 - this->H * eps_p.old(), 0.)}; // return mapping Real Del_gamma{phi_star/(this->H + 3 * this->mu)}; auto && tau{tau_star - 2*Del_gamma*this->mu*N_star}; /////auto && tau_eq{tau_eq_star - 3*this->mu*Del_gamma}; // update the previous values to the new ones F_prev.current() = F; ln_be_star -= 2*Del_gamma*N_star; be_prev.current() = expm(std::move(ln_be_star)); eps_p.current() += Del_gamma; // transmit info whether this is a plastic step or not bool is_plastic{phi_star > 0}; return std::tuple (std::move(tau), std::move(tau_eq_star), std::move(Del_gamma), std::move(N_star), std::move(is_plastic)); } //----------------------------------------------------------------------------// template template - decltype(auto) + auto MaterialHyperElastoPlastic1:: evaluate_stress(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, - FlowStRef_t eps_p) { + FlowStRef_t eps_p) -> decltype(auto) { auto retval(std::move(std::get<0>(this->stress_n_internals_worker (std::forward(F), F_prev, be_prev, eps_p)))); return retval; } //----------------------------------------------------------------------------// template template - decltype(auto) + auto MaterialHyperElastoPlastic1:: evaluate_stress_tangent(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, - FlowStRef_t eps_p) { + FlowStRef_t eps_p) -> decltype(auto) { //! after the stress computation, all internals are up to date auto && vals{this->stress_n_internals_worker (std::forward(F), F_prev, be_prev, eps_p)}; auto && tau {std::get<0>(vals)}; auto && tau_eq_star{std::get<1>(vals)}; auto && Del_gamma {std::get<2>(vals)}; auto && N_star {std::get<3>(vals)}; auto && is_plastic {std::get<4>(vals)}; if (is_plastic) { auto && a0 = Del_gamma*this->mu/tau_eq_star; auto && a1 = this->mu/(this->H + 3*this->mu); return std::make_tuple(std::move(tau), T4Mat{ ((this->K/2. - this->mu/3 + a0*this->mu)*Matrices::Itrac() + (1 - 3*a0) * this->mu*Matrices::Isymm() + 2*this->mu * (a0 - a1)*Matrices::outer(N_star, N_star))}); } else { return std::make_tuple(std::move(tau), T4Mat{this->C}); } } } // muSpectre #endif /* MATERIAL_HYPER_ELASTO_PLASTIC1_H */ diff --git a/src/materials/material_linear_elastic1.hh b/src/materials/material_linear_elastic1.hh index f963853..c5d822c 100644 --- a/src/materials/material_linear_elastic1.hh +++ b/src/materials/material_linear_elastic1.hh @@ -1,189 +1,191 @@ /** * @file material_linear_elastic1.hh * * @author Till Junge * * @date 13 Nov 2017 * * @brief Implementation for linear elastic reference material like in de Geus * 2017. This follows the simplest and likely not most efficient * implementation (with exception of the Python law) * * Copyright © 2017 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. */ #ifndef MATERIAL_LINEAR_ELASTIC1_H #define MATERIAL_LINEAR_ELASTIC1_H #include "common/common.hh" #include "materials/material_muSpectre_base.hh" #include "materials/materials_toolbox.hh" namespace muSpectre { template class MaterialLinearElastic1; /** * traits for objective linear elasticity */ template struct MaterialMuSpectre_traits> { using Parent = MaterialMuSpectre_traits;//!< base for elasticity //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! 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}; //! elasticity without internal variables using InternalVariables = std::tuple<>; }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) /** * implements objective linear elasticity */ template class MaterialLinearElastic1: public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! global field collection using Stiffness_t = Eigen::TensorFixedSize >; //! traits of this material using traits = MaterialMuSpectre_traits; //! this law does not have any internal variables using InternalVariables = typename traits::InternalVariables; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! Default constructor MaterialLinearElastic1() = delete; //! Copy constructor MaterialLinearElastic1(const MaterialLinearElastic1 &other) = delete; //! Construct by name, Young's modulus and Poisson's ratio MaterialLinearElastic1(std::string name, Real young, Real poisson); //! Move constructor MaterialLinearElastic1(MaterialLinearElastic1 &&other) = delete; //! Destructor virtual ~MaterialLinearElastic1() = default; //! Copy assignment operator MaterialLinearElastic1& operator=(const MaterialLinearElastic1 &other) = delete; //! Move assignment operator MaterialLinearElastic1& operator=(MaterialLinearElastic1 &&other) = delete; /** * 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); /** * 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); /** * return the empty internals tuple */ InternalVariables & get_internals() { return this->internal_variables;}; protected: const Real young; //!< Young's modulus const Real poisson;//!< Poisson's ratio const Real lambda; //!< first Lamé constant const Real mu; //!< second Lamé constant (shear modulus) const Stiffness_t C; //!< stiffness tensor //! empty tuple InternalVariables internal_variables{}; private: }; /* ---------------------------------------------------------------------- */ template template - decltype(auto) - MaterialLinearElastic1::evaluate_stress(s_t && E) { + auto + MaterialLinearElastic1::evaluate_stress(s_t && E) + -> decltype(auto) { return Hooke::evaluate_stress(this->lambda, this->mu, std::move(E)); } /* ---------------------------------------------------------------------- */ template template - decltype(auto) - MaterialLinearElastic1::evaluate_stress_tangent(s_t && E) { + auto + MaterialLinearElastic1::evaluate_stress_tangent(s_t && E) + -> decltype(auto) { using Tangent_t = typename traits::TangentMap_t::reference; // using mat = Eigen::Matrix; // mat ecopy{E}; // std::cout << "E" << std::endl << ecopy << std::endl; // std::cout << "P1" << std::endl << mat{ // std::get<0>(Hooke::evaluate_stress(this->lambda, this->mu, // Tangent_t(const_cast(this->C.data())), // std::move(E)))} << std::endl; return Hooke::evaluate_stress(this->lambda, this->mu, Tangent_t(const_cast(this->C.data())), std::move(E)); } } // muSpectre #endif /* MATERIAL_LINEAR_ELASTIC1_H */ diff --git a/src/materials/material_linear_elastic2.hh b/src/materials/material_linear_elastic2.hh index 58a9b7e..f90c405 100644 --- a/src/materials/material_linear_elastic2.hh +++ b/src/materials/material_linear_elastic2.hh @@ -1,199 +1,199 @@ /** * @file material_linear_elastic2.hh * * @author Till Junge * * @date 03 Feb 2018 * * @brief linear elastic material with imposed eigenstrain and its * type traits. Uses the MaterialMuSpectre facilities to keep it * simple * * 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. */ #ifndef MATERIAL_LINEAR_ELASTIC_EIGENSTRAIN_H #define MATERIAL_LINEAR_ELASTIC_EIGENSTRAIN_H #include "materials/material_linear_elastic1.hh" #include "common/field.hh" #include namespace muSpectre { template class MaterialLinearElastic2; /** * traits for objective linear elasticity with eigenstrain */ template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! 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; }; /** * implements objective linear elasticity with an eigenstrain per pixel */ template class MaterialLinearElastic2: public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre; //! type for stiffness tensor construction using Stiffness_t = Eigen::TensorFixedSize >; //! traits of this material using traits = MaterialMuSpectre_traits; //! Type of container used for storing eigenstrain using InternalVariables = typename traits::InternalVariables; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! reference to any type that casts to a matrix using StrainTensor = Eigen::Ref>; //! Default constructor MaterialLinearElastic2() = delete; //! Construct by name, Young's modulus and Poisson's ratio MaterialLinearElastic2(std::string name, Real young, Real poisson); //! Copy constructor MaterialLinearElastic2(const MaterialLinearElastic2 &other) = delete; //! Move constructor MaterialLinearElastic2(MaterialLinearElastic2 &&other) = delete; //! Destructor virtual ~MaterialLinearElastic2() = default; //! Copy assignment operator MaterialLinearElastic2& operator=(const MaterialLinearElastic2 &other) = delete; //! Move assignment operator MaterialLinearElastic2& operator=(MaterialLinearElastic2 &&other) = delete; /** * 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 */ InternalVariables & get_internals() { return this->internal_variables;}; /** * overload add_pixel to write into eigenstrain */ void add_pixel(const Ccoord_t & pixel) override final; /** * overload add_pixel to write into eigenstrain */ void add_pixel(const Ccoord_t & pixel, const StrainTensor & E_eig); protected: //! linear material without eigenstrain used to compute response MaterialLinearElastic1 material; //! storage for eigenstrain using Field_t = TensorField, Real, secondOrder, DimM>; Field_t & eigen_field; //!< field holding the eigen strain per pixel //! tuple for iterable eigen_field InternalVariables internal_variables; private: }; /* ---------------------------------------------------------------------- */ template template - decltype(auto) + auto MaterialLinearElastic2:: - evaluate_stress(s_t && E, eigen_s_t && E_eig) { + evaluate_stress(s_t && E, eigen_s_t && E_eig) -> decltype(auto) { return this->material.evaluate_stress(E-E_eig); } /* ---------------------------------------------------------------------- */ template template - decltype(auto) + auto MaterialLinearElastic2:: - evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig) { + evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig) -> decltype(auto) { // using mat = Eigen::Matrix; // mat ecopy{E}; // mat eig_copy{E_eig}; // mat ediff{ecopy-eig_copy}; // std::cout << "eidff - (E-E_eig)" << std::endl << ediff-(E-E_eig) << std::endl; // std::cout << "P1 " << std::endl << mat{std::get<0>(this->material.evaluate_stress_tangent(E-E_eig))} << "" << std::endl; // std::cout << "P2" << std::endl << mat{std::get<0>(this->material.evaluate_stress_tangent(std::move(ediff)))} << std::endl; return this->material.evaluate_stress_tangent(E-E_eig); } } // muSpectre #endif /* MATERIAL_LINEAR_ELASTIC_EIGENSTRAIN_H */ diff --git a/src/materials/material_linear_elastic3.hh b/src/materials/material_linear_elastic3.hh index 3b27cba..c724058 100644 --- a/src/materials/material_linear_elastic3.hh +++ b/src/materials/material_linear_elastic3.hh @@ -1,196 +1,196 @@ /** * @file material_linear_elastic3.hh * * @author Richard Leute * * @date 20 Feb 2018 * * @brief linear elastic material with distribution of stiffness properties. * Uses the MaterialMuSpectre facilities to keep it simple. * * 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. */ #ifndef MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_H #define MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_H #include "materials/material_linear_elastic1.hh" #include "common/field.hh" #include "common/tensor_algebra.hh" #include namespace muSpectre { template class MaterialLinearElastic3; /** * traits for objective linear elasticity with eigenstrain */ template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! 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 stiffness tensor type using LStiffnessMap_t = T4MatrixFieldMap; //! elasticity without internal variables using InternalVariables = std::tuple; }; /** * implements objective linear elasticity with an eigenstrain per pixel */ template class MaterialLinearElastic3: public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! global field collection using Stiffness_t = Eigen::TensorFixedSize >; //! traits of this material using traits = MaterialMuSpectre_traits; //! Type of container used for storing eigenstrain using InternalVariables = typename traits::InternalVariables; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! Default constructor MaterialLinearElastic3() = delete; //! Construct by name MaterialLinearElastic3(std::string name); //! Copy constructor MaterialLinearElastic3(const MaterialLinearElastic3 &other) = delete; //! Move constructor MaterialLinearElastic3(MaterialLinearElastic3 &&other) = delete; //! Destructor virtual ~MaterialLinearElastic3() = default; //! Copy assignment operator MaterialLinearElastic3& operator=(const MaterialLinearElastic3 &other) = delete; //! Move assignment operator MaterialLinearElastic3& operator=(MaterialLinearElastic3 &&other) = delete; /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor) * and the local stiffness tensor. */ template inline decltype(auto) evaluate_stress(s_t && E, stiffness_t && C); /** * 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) and the local stiffness tensor. */ template inline decltype(auto) evaluate_stress_tangent(s_t && E, stiffness_t && C); /** * return the empty internals tuple */ InternalVariables & get_internals() { return this->internal_variables;}; /** * overload add_pixel to write into loacal stiffness tensor */ void add_pixel(const Ccoord_t & pixel) override final; /** * overload add_pixel to write into local stiffness tensor */ void add_pixel(const Ccoord_t & pixel, const Real & Young, const Real & PoissonRatio); protected: //! storage for stiffness tensor using Field_t = TensorField, Real, fourthOrder, DimM>; Field_t & C_field; //!< field of stiffness tensors //! tuple for iterable eigen_field InternalVariables internal_variables; private: }; /* ---------------------------------------------------------------------- */ template template - decltype(auto) + auto MaterialLinearElastic3:: - evaluate_stress(s_t && E, stiffness_t && C) { + evaluate_stress(s_t && E, stiffness_t && C) -> decltype(auto) { return Matrices::tensmult(C, E); } /* ---------------------------------------------------------------------- */ template template - decltype(auto) + auto MaterialLinearElastic3:: - evaluate_stress_tangent(s_t && E, stiffness_t && C) { + evaluate_stress_tangent(s_t && E, stiffness_t && C) -> decltype(auto) { return std::make_tuple (evaluate_stress(E, C), C); } } // muSpectre #endif /* MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_H */ diff --git a/src/materials/material_linear_elastic4.hh b/src/materials/material_linear_elastic4.hh index 6c562a3..a5aee8b 100644 --- a/src/materials/material_linear_elastic4.hh +++ b/src/materials/material_linear_elastic4.hh @@ -1,207 +1,208 @@ /** * @file material_linear_elastic4.hh * * @author Richard Leute * * @date 15 March 2018 * * @brief linear elastic material with distribution of stiffness properties. * In difference to material_linear_elastic3 two Lame constants are * stored per pixel instead of the whole elastic matrix C. * Uses the MaterialMuSpectre facilities to keep it simple. * * 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. */ #ifndef MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_2_H #define MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_2_H #include "materials/material_linear_elastic1.hh" #include "common/field.hh" #include "common/tensor_algebra.hh" #include namespace muSpectre { template class MaterialLinearElastic4; /** * traits for objective linear elasticity with eigenstrain */ template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! 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 Lame constant type using LLameConstantMap_t = ScalarFieldMap; //! elasticity without internal variables using InternalVariables = std::tuple; }; /** * implements objective linear elasticity with an eigenstrain per pixel */ template class MaterialLinearElastic4: public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! global field collection using Stiffness_t = Eigen::TensorFixedSize >; //! traits of this material using traits = MaterialMuSpectre_traits; //! Type of container used for storing eigenstrain using InternalVariables = typename traits::InternalVariables; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! Default constructor MaterialLinearElastic4() = delete; //! Construct by name MaterialLinearElastic4(std::string name); //! Copy constructor MaterialLinearElastic4(const MaterialLinearElastic4 &other) = delete; //! Move constructor MaterialLinearElastic4(MaterialLinearElastic4 &&other) = delete; //! Destructor virtual ~MaterialLinearElastic4() = default; //! Copy assignment operator MaterialLinearElastic4& operator=(const MaterialLinearElastic4 &other) = delete; //! Move assignment operator MaterialLinearElastic4& operator=(MaterialLinearElastic4 &&other) = delete; /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor), the first * Lame constant (lambda) and the second Lame constant (shear modulus/mu). */ template inline decltype(auto) evaluate_stress(s_t && E, const Real & lambda, const Real & mu); /** * 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), the first Lame constant (lambda) and * the second Lame constant (shear modulus/mu). */ template inline decltype(auto) evaluate_stress_tangent(s_t && E, const Real & lambda, const Real & mu); /** * return the empty internals tuple */ InternalVariables & get_internals() { return this->internal_variables;}; /** * overload add_pixel to write into loacal stiffness tensor */ void add_pixel(const Ccoord_t & pixel) override final; /** * overload add_pixel to write into local stiffness tensor */ void add_pixel(const Ccoord_t & pixel, const Real & Poisson_ratio, const Real & Youngs_modulus); protected: //! storage for first Lame constant 'lambda' //! and second Lame constant(shear modulus) 'mu' using Field_t = MatrixField, Real, oneD, oneD>; Field_t & lambda_field; Field_t & mu_field; //! tuple for iterable eigen_field InternalVariables internal_variables; private: }; /* ---------------------------------------------------------------------- */ template template - decltype(auto) + auto MaterialLinearElastic4:: - evaluate_stress(s_t && E, const Real & lambda, const Real & mu) { + evaluate_stress(s_t && E, const Real & lambda, const Real & mu) + -> decltype(auto) { auto C = Hooke::compute_C_T4(lambda, mu); return Matrices::tensmult(C, E); } /* ---------------------------------------------------------------------- */ template template - decltype(auto) + auto MaterialLinearElastic4:: evaluate_stress_tangent(s_t && E, const Real & lambda, - const Real & mu) + const Real & mu) -> decltype(auto) { auto C = Hooke::compute_C_T4(lambda, mu); return std::make_tuple(Matrices::tensmult(C, E), C); } } // muSpectre #endif /* MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_2_H */ diff --git a/src/solver/solver_base.cc b/src/solver/solver_base.cc index ff3b889..7b19cad 100644 --- a/src/solver/solver_base.cc +++ b/src/solver/solver_base.cc @@ -1,65 +1,63 @@ /** * file solver_base.cc * * @author Till Junge * * @date 24 Apr 2018 * * @brief implementation of SolverBase * - * @section LICENSE - * * 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 "solver/solver_base.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ SolverBase::SolverBase(Cell & cell, Real tol, Uint maxiter, bool verbose): cell(cell), tol{tol}, maxiter{maxiter}, verbose{verbose} {} /* ---------------------------------------------------------------------- */ bool SolverBase::has_converged() const { return this->converged; } /* ---------------------------------------------------------------------- */ void SolverBase::reset_counter() { this->counter = 0; this->converged = false; } /* ---------------------------------------------------------------------- */ Uint SolverBase::get_counter() const { return this->counter; } /* ---------------------------------------------------------------------- */ Real SolverBase::get_tol() const { return this->tol; } /* ---------------------------------------------------------------------- */ Uint SolverBase::get_maxiter() const { return this->maxiter; } } // muSpectre diff --git a/src/solver/solver_base.hh b/src/solver/solver_base.hh index 331ecbf..49ac024 100644 --- a/src/solver/solver_base.hh +++ b/src/solver/solver_base.hh @@ -1,120 +1,118 @@ /** * file solver_base.hh * * @author Till Junge * * @date 24 Apr 2018 * * @brief Base class for iterative solvers for linear systems of equations * - * @section LICENSE - * * 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. */ #ifndef SOLVER_BASE_H #define SOLVER_BASE_H #include "solver/solver_common.hh" #include "cell/cell_base.hh" #include namespace muSpectre { /** * Virtual base class for solvers. An implementation of this interface * can be used with the solution strategies in solvers.hh */ class SolverBase { public: //! underlying vector type using Vector_t = Eigen::Matrix; //! Input vector for solvers using Vector_ref = Eigen::Ref; //! Input vector for solvers using ConstVector_ref = Eigen::Ref; //! Output vector for solvers using Vector_map = Eigen::Map; //! Default constructor SolverBase() = delete; /** * Constructor takes a Cell, tolerance, max number of iterations * and verbosity flag as input */ SolverBase(Cell & cell, Real tol, Uint maxiter, bool verbose=false); //! Copy constructor SolverBase(const SolverBase &other) = delete; //! Move constructor SolverBase(SolverBase &&other) = default; //! Destructor virtual ~SolverBase() = default; //! Copy assignment operator SolverBase& operator=(const SolverBase &other) = delete; //! Move assignment operator SolverBase& operator=(SolverBase &&other) = default; //! Allocate fields used during the solution virtual void initialise() = 0; //! returns whether the solver has converged bool has_converged() const ; //! reset the iteration counter to zero void reset_counter(); //! get the count of how many solve steps have been executed since //! construction of most recent counter reset Uint get_counter() const; //! returns the max number of iterations Uint get_maxiter() const; //! returns the resolution tolerance Real get_tol() const; //! returns the solver's name (i.e. 'CG', 'GMRES', etc) virtual std::string get_name() const = 0; //! run the solve operation virtual Vector_map solve(const ConstVector_ref rhs) = 0; protected: Cell & cell; //!< reference to the problem's cell Real tol; //!< convergence tolerance Uint maxiter; //!< maximum allowed number of iterations bool verbose; //!< whether to write information to the stdout Uint counter{0}; //!< iteration counter bool converged{false}; //!< whether the solver has converged private: }; } // muSpectre #endif /* SOLVER_BASE_H */ diff --git a/src/solver/solver_cg.cc b/src/solver/solver_cg.cc index 2b38b29..e17417b 100644 --- a/src/solver/solver_cg.cc +++ b/src/solver/solver_cg.cc @@ -1,109 +1,107 @@ /** * file solver_cg.cc * * @author Till Junge * * @date 24 Apr 2018 * * @brief implements SolverCG * - * @section LICENSE - * * 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 "solver/solver_cg.hh" #include "common/communicator.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ SolverCG::SolverCG(Cell & cell, Real tol, Uint maxiter, bool verbose): Parent(cell, tol, maxiter, verbose), r_k(cell.get_nb_dof()), p_k(cell.get_nb_dof()), Ap_k(cell.get_nb_dof()), x_k(cell.get_nb_dof()) {} /* ---------------------------------------------------------------------- */ auto SolverCG::solve(const ConstVector_ref rhs) -> Vector_map { this->x_k.setZero(); const Communicator & comm = this->cell.get_communicator(); // Following implementation of algorithm 5.2 in Nocedal's // Numerical Optimization (p. 112) //initialisation of algorithm this->r_k = (this->cell.evaluate_projected_directional_stiffness(this->x_k) - rhs); this->p_k = -this->r_k; this->converged = false; Real rdr = comm.sum(this->r_k.dot(this->r_k)); Real rhs_norm2 = comm.sum(rhs.squaredNorm()); Real tol2 = ipow(this->tol, 2) * rhs_norm2; size_t count_width{}; // for output formatting in verbose case if (this->verbose) { count_width = size_t(std::log10(this->maxiter))+1; } for (Uint i = 0; i < this->maxiter && (rdr > tol2 || i == 0); ++i, ++this->counter) { this->Ap_k = this->cell.evaluate_projected_directional_stiffness(this->p_k); Real alpha = rdr/comm.sum(this->p_k.dot(this->Ap_k)); this->x_k += alpha * this->p_k; this->r_k += alpha * this->Ap_k; Real new_rdr = comm.sum(this->r_k.dot(this->r_k)); Real beta = new_rdr/rdr; rdr = new_rdr; if (this->verbose && comm.rank() == 0) { std::cout << " at CG step " << std::setw(count_width) << i << ": |r|/|b| = " << std::setw(15) << std::sqrt(rdr/rhs_norm2) << ", cg_tol = " << this->tol << std::endl; } this->p_k = - this->r_k + beta * this->p_k; } if (rdr < tol2) { this->converged = true; } else { std::stringstream err{}; err << " After " << this->counter << " steps, the solver " << " FAILED with |r|/|b| = " << std::setw(15) << std::sqrt(rdr/rhs_norm2) << ", cg_tol = " << this->tol << std::endl; throw ConvergenceError("Conjugate gradient has not converged." + err.str()); } return Vector_map(this->x_k.data(), this->x_k.size()); } } // muSpectre diff --git a/src/solver/solver_cg.hh b/src/solver/solver_cg.hh index 3a803ea..bc1cc18 100644 --- a/src/solver/solver_cg.hh +++ b/src/solver/solver_cg.hh @@ -1,105 +1,103 @@ /** * file solver_cg.hh * * @author Till Junge * * @date 24 Apr 2018 * * @brief class fo a simple implementation of a conjugate gradient solver. * This follows algorithm 5.2 in Nocedal's Numerical Optimization * (p 112) * - * @section LICENSE - * * 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. */ #ifndef SOLVER_CG_H #define SOLVER_CG_H #include "solver/solver_base.hh" namespace muSpectre { /** * implements the `muSpectre::SolverBase` interface using a * conjugate gradient solver. This particular class is useful for * trouble shooting, as it can be made very verbose, but for * production runs, it is probably better to use * `muSpectre::SolverCGEigen`. */ class SolverCG: public SolverBase { public: using Parent = SolverBase; //!< standard short-hand for base class //! for storage of fields using Vector_t = Parent::Vector_t; //! Input vector for solvers using Vector_ref = Parent::Vector_ref; //! Input vector for solvers using ConstVector_ref = Parent::ConstVector_ref; //! Output vector for solvers using Vector_map = Parent::Vector_map; //! Default constructor SolverCG() = delete; //! Copy constructor SolverCG(const SolverCG &other) = delete; /** * Constructor takes a Cell, tolerance, max number of iterations * and verbosity flag as input */ SolverCG(Cell & cell, Real tol, Uint maxiter, bool verbose=false); //! Move constructor SolverCG(SolverCG &&other) = default; //! Destructor virtual ~SolverCG() = default; //! Copy assignment operator SolverCG& operator=(const SolverCG &other) = delete; //! Move assignment operator SolverCG& operator=(SolverCG &&other) = default; //! initialisation does not need to do anything in this case void initialise() override final {}; //! returns the solver's name std::string get_name() const override final {return "CG";} //! the actual solver Vector_map solve(const ConstVector_ref rhs) override final; protected: Vector_t r_k; //!< residual Vector_t p_k; //!< search direction Vector_t Ap_k; //!< directional stiffness Vector_t x_k; //!< current solution private: }; } // muSpectre #endif /* SOLVER_CG_H */ diff --git a/src/solver/solver_common.cc b/src/solver/solver_common.cc index 85c0cbc..88372cd 100644 --- a/src/solver/solver_common.cc +++ b/src/solver/solver_common.cc @@ -1,42 +1,40 @@ /** * file solver_common.cc * * @author Till Junge * * @date 15 May 2018 * * @brief implementation for solver utilities * - * @section LICENSE - * * 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 "solver/solver_common.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ bool check_symmetry(const Eigen::Ref& eps, Real rel_tol){ return (rel_tol >= (eps-eps.transpose()).matrix().norm()/eps.matrix().norm() || rel_tol >= eps.matrix().norm()); } } // muSpectre diff --git a/src/solver/solver_eigen.cc b/src/solver/solver_eigen.cc index 494068d..4872dd7 100644 --- a/src/solver/solver_eigen.cc +++ b/src/solver/solver_eigen.cc @@ -1,93 +1,91 @@ /** * file solver_eigen.cc * * @author Till Junge * * @date 15 May 2018 * * @brief Implementations for bindings to Eigen's iterative solvers * - * @section LICENSE - * * 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 "solver/solver_eigen.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template SolverEigen::SolverEigen(Cell& cell, Real tol, Uint maxiter, bool verbose): Parent(cell, tol, maxiter, verbose), adaptor{cell.get_adaptor()}, solver{}, result{} {} /* ---------------------------------------------------------------------- */ template void SolverEigen::initialise() { this->solver.setTolerance(this->get_tol()); this->solver.setMaxIterations(this->get_maxiter()); this->solver.compute(this->adaptor); } /* ---------------------------------------------------------------------- */ template auto SolverEigen::solve(const ConstVector_ref rhs) -> Vector_map { // for crtp auto & this_solver = static_cast (*this); this->result = this->solver.solve(rhs); this->counter += this->solver.iterations(); if (this->solver.info() != Eigen::Success) { std::stringstream err {}; err << this_solver.get_name() << " has not converged," << " After " << this->solver.iterations() << " steps, the solver " << " FAILED with |r|/|b| = " << std::setw(15) << this->solver.error() << ", cg_tol = " << this->tol << std::endl; throw ConvergenceError(err.str()); } if (this->verbose) { std::cout << " After " << this->solver.iterations() << " " << this_solver.get_name() << " steps, |r|/|b| = " << std::setw(15) << this->solver.error() << ", cg_tol = " << this->tol << std::endl; } return Vector_map(this->result.data(), this->result.size()); } /* ---------------------------------------------------------------------- */ template class SolverEigen; template class SolverEigen; template class SolverEigen; template class SolverEigen; template class SolverEigen; } // muSpectre diff --git a/src/solver/solver_eigen.hh b/src/solver/solver_eigen.hh index c4a6cdc..ed982a5 100644 --- a/src/solver/solver_eigen.hh +++ b/src/solver/solver_eigen.hh @@ -1,215 +1,213 @@ /** * file solver_eigen.hh * * @author Till Junge * * @date 15 May 2018 * * @brief Bindings to Eigen's iterative solvers * - * @section LICENSE - * * 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. */ #ifndef SOLVER_EIGEN_H #define SOLVER_EIGEN_H #include "solver/solver_base.hh" #include "cell/cell_base.hh" #include #include namespace muSpectre { template class SolverEigen; class SolverCGEigen; class SolverGMRESEigen; class SolverBiCGSTABEigen; class SolverDGMRESEigen; class SolverMINRESEigen; namespace internal { template struct Solver_traits { }; //! traits for the Eigen conjugate gradient solver template<> struct Solver_traits { //! Eigen Iterative Solver using Solver = Eigen::ConjugateGradient; }; //! traits for the Eigen GMRES solver template<> struct Solver_traits { //! Eigen Iterative Solver using Solver = Eigen::GMRES; }; //! traits for the Eigen BiCGSTAB solver template<> struct Solver_traits { //! Eigen Iterative Solver using Solver = Eigen::BiCGSTAB; }; //! traits for the Eigen DGMRES solver template<> struct Solver_traits { //! Eigen Iterative Solver using Solver = Eigen::DGMRES; }; //! traits for the Eigen MINRES solver template<> struct Solver_traits { //! Eigen Iterative Solver using Solver = Eigen::MINRES; }; } // internal /** * base class for iterative solvers from Eigen */ template class SolverEigen: public SolverBase { public: using Parent = SolverBase; //!< base class //! traits obtained from CRTP using Solver = typename internal::Solver_traits::Solver; //! Input vectors for solver using ConstVector_ref = Parent::ConstVector_ref; //! Output vector for solver using Vector_map = Parent::Vector_map; //! storage for output vector using Vector_t = Parent::Vector_t; //! Default constructor SolverEigen() = delete; //! Constructor with domain resolutions, etc, SolverEigen(Cell& cell, Real tol, Uint maxiter=0, bool verbose =false); //! Copy constructor SolverEigen(const SolverEigen &other) = delete; //! Move constructor SolverEigen(SolverEigen &&other) = default; //! Destructor virtual ~SolverEigen() = default; //! Copy assignment operator SolverEigen& operator=(const SolverEigen &other) = delete; //! Move assignment operator SolverEigen& operator=(SolverEigen &&other) = default; //! Allocate fields used during the solution void initialise() override final; //! executes the solver Vector_map solve(const ConstVector_ref rhs) override final; protected: Cell::Adaptor adaptor; //!< cell handle Solver solver; //!< Eigen's Iterative solver Vector_t result; //!< storage for result }; /** * Binding to Eigen's conjugate gradient solver */ class SolverCGEigen: public SolverEigen { public: using SolverEigen::SolverEigen; std::string get_name() const override final {return "CG";} }; /** * Binding to Eigen's GMRES solver */ class SolverGMRESEigen: public SolverEigen { public: using SolverEigen::SolverEigen; std::string get_name() const override final {return "GMRES";} }; /** * Binding to Eigen's BiCGSTAB solver */ class SolverBiCGSTABEigen: public SolverEigen { public: using SolverEigen::SolverEigen; //! Solver's name std::string get_name() const override final {return "BiCGSTAB";} }; /** * Binding to Eigen's DGMRES solver */ class SolverDGMRESEigen: public SolverEigen { public: using SolverEigen::SolverEigen; //! Solver's name std::string get_name() const override final {return "DGMRES";} }; /** * Binding to Eigen's MINRES solver */ class SolverMINRESEigen: public SolverEigen { public: using SolverEigen::SolverEigen; //! Solver's name std::string get_name() const override final {return "MINRES";} }; } // muSpectre #endif /* SOLVER_EIGEN_H */ diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc index 53f0b0b..c4e19bc 100644 --- a/src/solver/solvers.cc +++ b/src/solver/solvers.cc @@ -1,404 +1,402 @@ /** * file solvers.cc * * @author Till Junge * * @date 24 Apr 2018 * * @brief implementation of dynamic newton-cg solver * - * @section LICENSE - * * 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 "solvers.hh" #include #include namespace muSpectre { //----------------------------------------------------------------------------// std::vector newton_cg(Cell & cell, const LoadSteps_t & load_steps, SolverBase & solver, Real newton_tol, Real equil_tol, Dim_t verbose) { const Communicator & comm = cell.get_communicator(); using Vector_t = Eigen::Matrix; using Matrix_t = Eigen::Matrix; // Corresponds to symbol δF or δε Vector_t incrF(cell.get_nb_dof()); // field to store the rhs for cg calculations Vector_t rhs(cell.get_nb_dof()); solver.initialise(); size_t count_width{}; const auto form{cell.get_formulation()}; std::string strain_symb{}; if (verbose > 0 && comm.rank() == 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "Newton-" << solver.get_name() << " for "; switch (form) { case Formulation::small_strain: { strain_symb = "ε"; std::cout << "small"; break; } case Formulation::finite_strain: { strain_symb = "F"; std::cout << "finite"; break; } default: throw SolverError("unknown formulation"); break; } std::cout << " strain with" << std::endl << "newton_tol = " << newton_tol << ", cg_tol = " << solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ" << strain_symb << " =" <(tup)}; auto && grad{std::get<1>(tup)}; std::cout << "Step " << counter + 1 << ":" << std::endl << grad << std::endl; } count_width = size_t(std::log10(solver.get_maxiter()))+1; } auto shape{cell.get_strain_shape()}; switch (form) { case Formulation::finite_strain: { cell.set_uniform_strain(Matrix_t::Identity(shape[0], shape[1])); for (const auto & delF: load_steps) { if (not ((delF.rows() == shape[0]) and (delF.cols() == shape[1]))) { std::stringstream err{}; err << "Load increments need to be given in " << shape[0] << "×" << shape[1] << " matrices, but I got a " << delF.rows() << "×" << delF.cols() << " matrix:" << std::endl << delF; throw SolverError(err.str()); } } break; } case Formulation::small_strain: { cell.set_uniform_strain(Matrix_t::Zero(shape[0], shape[1])); for (const auto & delF: load_steps) { if (not ((delF.rows() == shape[0]) and (delF.cols() == shape[1]))) { std::stringstream err{}; err << "Load increments need to be given in " << shape[0] << "×" << shape[1] << " matrices, but I got a " << delF.rows() << "×" << delF.cols() << " matrix:" << std::endl << delF; throw SolverError(err.str()); } if (not check_symmetry(delF)) { throw SolverError("all Δε must be symmetric!"); } } break; } default: throw SolverError("Unknown strain measure"); break; } // initialise return value std::vector ret_val{}; // storage for the previous mean strain (to compute ΔF or Δε) Matrix_t previous_macro_strain{load_steps.back().Zero(shape[0], shape[1])}; auto F{cell.get_strain_vector()}; //! incremental loop for (const auto & macro_strain: load_steps) { using StrainMap_t = RawFieldMap>; for (auto && strain: StrainMap_t(F, shape[0], shape[1])) { strain += macro_strain - previous_macro_strain; } std::string message{"Has not converged"}; Real incr_norm{2*newton_tol}, grad_norm{1}; Real stress_norm{2*equil_tol}; bool has_converged{false}; auto convergence_test = [&incr_norm, &grad_norm, &newton_tol, &stress_norm, &equil_tol, &message, &has_converged] () { bool incr_test = incr_norm/grad_norm <= newton_tol; bool stress_test = stress_norm < equil_tol; if (incr_test) { message = "Residual tolerance reached"; } else if (stress_test) { message = "Reached stress divergence tolerance"; } has_converged = incr_test || stress_test; return has_converged; }; Uint newt_iter{0}; for (; newt_iter < solver.get_maxiter() && !has_converged; ++newt_iter) { // obtain material response auto res_tup{cell.evaluate_stress_tangent()}; auto & P{std::get<0>(res_tup)}; rhs = -P; cell.apply_projection(rhs); stress_norm = std::sqrt(comm.sum(rhs.squaredNorm())); if (convergence_test()) { break; } //! this is a potentially avoidable copy TODO: check this out incrF = solver.solve(rhs); F += incrF; incr_norm = std::sqrt(comm.sum(incrF.squaredNorm())); grad_norm = std::sqrt(comm.sum(F.squaredNorm())); if ((verbose > 0) and (comm.rank() == 0)) { std::cout << "at Newton step " << std::setw(count_width) << newt_iter << ", |δ" << strain_symb << "|/|Δ" << strain_symb << "| = " << std::setw(17) << incr_norm/grad_norm << ", tol = " << newton_tol << std::endl; if (verbose-1>1) { std::cout << "<" << strain_symb << "> =" << std::endl << StrainMap_t(F, shape[0], shape[1]).mean() << std::endl; } } convergence_test(); } // update previous macroscopic strain previous_macro_strain = macro_strain; // store results ret_val.emplace_back(OptimizeResult{F, cell.get_stress_vector(), convergence_test(), Int(convergence_test()), message, newt_iter, solver.get_counter()}); // store history variables for next load increment cell.save_history_variables(); } return ret_val; } //----------------------------------------------------------------------------// std::vector de_geus(Cell & cell, const LoadSteps_t & load_steps, SolverBase & solver, Real newton_tol, Real equil_tol, Dim_t verbose) { const Communicator & comm = cell.get_communicator(); using Vector_t = Eigen::Matrix; using Matrix_t = Eigen::Matrix; // Corresponds to symbol δF or δε Vector_t incrF(cell.get_nb_dof()); // Corresponds to symbol ΔF or Δε Vector_t DeltaF(cell.get_nb_dof()); // field to store the rhs for cg calculations Vector_t rhs(cell.get_nb_dof()); solver.initialise(); size_t count_width{}; const auto form{cell.get_formulation()}; std::string strain_symb{}; if (verbose > 0 && comm.rank() == 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "de Geus-" << solver.get_name() << " for "; switch (form) { case Formulation::small_strain: { strain_symb = "ε"; std::cout << "small"; break; } case Formulation::finite_strain: { strain_symb = "F"; std::cout << "finite"; break; } default: throw SolverError("unknown formulation"); break; } std::cout << " strain with" << std::endl << "newton_tol = " << newton_tol << ", cg_tol = " << solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ" << strain_symb << " =" <(tup)}; auto && grad{std::get<1>(tup)}; std::cout << "Step " << counter + 1 << ":" << std::endl << grad << std::endl; } count_width = size_t(std::log10(solver.get_maxiter()))+1; } auto shape{cell.get_strain_shape()}; switch (form) { case Formulation::finite_strain: { cell.set_uniform_strain(Matrix_t::Identity(shape[0], shape[1])); for (const auto & delF: load_steps) { auto rows = delF.rows(); auto cols = delF.cols(); if (not ((rows == shape[0]) and (cols == shape[1]))) { std::stringstream err{}; err << "Load increments need to be given in " << shape[0] << "×" << shape[1] << " matrices, but I got a " << delF.rows() << "×" << delF.cols() << " matrix:" << std::endl << delF; throw SolverError(err.str()); } } break; } case Formulation::small_strain: { cell.set_uniform_strain(Matrix_t::Zero(shape[0], shape[1])); for (const auto & delF: load_steps) { if (not ((delF.rows() == shape[0]) and (delF.cols() == shape[1]))) { std::stringstream err{}; err << "Load increments need to be given in " << shape[0] << "×" << shape[1] << " matrices, but I got a " << delF.rows() << "×" << delF.cols() << " matrix:" << std::endl << delF; throw SolverError(err.str()); } if (not check_symmetry(delF)) { throw SolverError("all Δε must be symmetric!"); } } break; } default: throw SolverError("Unknown strain measure"); break; } // initialise return value std::vector ret_val{}; // storage for the previous mean strain (to compute ΔF or Δε) Matrix_t previous_macro_strain{load_steps.back().Zero(shape[0], shape[1])}; auto F{cell.get_strain_vector()}; //! incremental loop for (const auto & macro_strain: load_steps) { using StrainMap_t = RawFieldMap>; std::string message{"Has not converged"}; Real incr_norm{2*newton_tol}, grad_norm{1}; Real stress_norm{2*equil_tol}; bool has_converged{false}; auto convergence_test = [&incr_norm, &grad_norm, &newton_tol, &stress_norm, &equil_tol, &message, &has_converged] () { bool incr_test = incr_norm/grad_norm <= newton_tol; bool stress_test = stress_norm < equil_tol; if (incr_test) { message = "Residual tolerance reached"; } else if (stress_test) { message = "Reached stress divergence tolerance"; } has_converged = incr_test || stress_test; return has_converged; }; Uint newt_iter{0}; for (; newt_iter < solver.get_maxiter() && !has_converged; ++newt_iter) { // obtain material response auto res_tup{cell.evaluate_stress_tangent()}; auto & P{std::get<0>(res_tup)}; if (newt_iter == 0) { for (auto && strain: StrainMap_t(DeltaF, shape[0], shape[1])) { strain = macro_strain -previous_macro_strain; } rhs = - cell.evaluate_projected_directional_stiffness(DeltaF); stress_norm = std::sqrt(comm.sum(rhs.matrix().squaredNorm())); if (convergence_test()) { break; } incrF = solver.solve(rhs); F += DeltaF; } else { rhs = -P; cell.apply_projection(rhs); stress_norm = std::sqrt(comm.sum(rhs.matrix().squaredNorm())); if (convergence_test()) { break; } incrF = solver.solve(rhs); } F += incrF; incr_norm = std::sqrt(comm.sum(incrF.squaredNorm())); grad_norm = std::sqrt(comm.sum(F.squaredNorm())); if ((verbose > 0) and (comm.rank() == 0)) { std::cout << "at Newton step " << std::setw(count_width) << newt_iter << ", |δ" << strain_symb << "|/|Δ" << strain_symb << "| = " << std::setw(17) << incr_norm/grad_norm << ", tol = " << newton_tol << std::endl; if (verbose-1>1) { std::cout << "<" << strain_symb << "> =" << std::endl << StrainMap_t(F, shape[0], shape[1]).mean() << std::endl; } } convergence_test(); } // update previous macroscopic strain previous_macro_strain = macro_strain; // store results ret_val.emplace_back(OptimizeResult{F, cell.get_stress_vector(), convergence_test(), Int(convergence_test()), message, newt_iter, solver.get_counter()}); // store history variables for next load increment cell.save_history_variables(); } return ret_val; } } // muSpectre diff --git a/src/solver/solvers.hh b/src/solver/solvers.hh index 391c7b6..64345be 100644 --- a/src/solver/solvers.hh +++ b/src/solver/solvers.hh @@ -1,100 +1,98 @@ /** * file solvers.hh * * @author Till Junge * * @date 24 Apr 2018 * * @brief Free functions for solving rve problems * - * @section LICENSE - * * 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. */ #ifndef SOLVERS_H #define SOLVERS_H #include "solver/solver_base.hh" #include #include #include namespace muSpectre { using LoadSteps_t = std::vector; /** * Uses the Newton-conjugate Gradient method to find the static * equilibrium of a cell given a series of mean applied strains */ std::vector newton_cg(Cell & cell, const LoadSteps_t & load_steps, SolverBase & solver, Real newton_tol, Real equil_tol, Dim_t verbose = 0); /** * Uses the Newton-conjugate Gradient method to find the static * equilibrium of a cell given a mean applied strain */ OptimizeResult newton_cg(Cell & cell, const Eigen::Ref load_step, SolverBase & solver, Real newton_tol, Real equil_tol, Dim_t verbose = 0) { LoadSteps_t load_steps{load_step}; return newton_cg(cell, load_steps, solver, newton_tol, equil_tol, verbose).front(); } /* ---------------------------------------------------------------------- */ /** * Uses the method proposed by de Geus method to find the static * equilibrium of a cell given a series of mean applied strains */ std::vector de_geus(Cell & cell, const LoadSteps_t & load_steps, SolverBase & solver, Real newton_tol, Real equil_tol, Dim_t verbose = 0); /* ---------------------------------------------------------------------- */ /** * Uses the method proposed by de Geus method to find the static * equilibrium of a cell given a mean applied strain */ OptimizeResult de_geus(Cell & cell, const Eigen::Ref load_step, SolverBase & solver, Real newton_tol, Real equil_tol, Dim_t verbose = 0){ return de_geus(cell, LoadSteps_t{load_step}, solver, newton_tol, equil_tol, verbose)[0]; } } // muSpectre #endif /* SOLVERS_H */ diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py index 7a8c2e1..afe329b 100755 --- a/tests/python_binding_tests.py +++ b/tests/python_binding_tests.py @@ -1,172 +1,173 @@ #!/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_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_field_tests.py b/tests/python_field_tests.py new file mode 100644 index 0000000..2850301 --- /dev/null +++ b/tests/python_field_tests.py @@ -0,0 +1,71 @@ +""" +file python_field_tests.py + +@author Till Junge + +@date 06 Jul 2018 + +@brief tests the python bindings for fieldcollections, fields, and statefields + +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 µ + +class FieldCollection_Check(unittest.TestCase): + """Because field collections do not have python-accessible + constructors, this test creates a problem with a material with + statefields + + """ + def setUp(self): + self.resolution = [3, 3] + self.lengths = [1.58, 5.87] + self.formulation = µ.Formulation.finite_strain + self.cell = µ.Cell(self.resolution, + self.lengths, + self.formulation) + self.dim = len(self.lengths) + self.mat = µ.material.MaterialLinearElastic2_2d.make( + self.cell, "material", 210e9, .33) + + def test_fields(self): + eigen_strain = np.array([[.01, .02], + [.03, -.01]]) + for i, pixel in enumerate(self.cell): + self.mat.add_pixel(pixel, i/self.cell.size*eigen_strain) + + self.cell.initialise() + dir(µ.material.Material_2d) + self.assertTrue(isinstance(self.mat, µ.material.Material_2d)) + collection = self.mat.collection + field_name = collection.field_names[0] + self.assertRaises(Exception, collection.get_complex_field, field_name) + self.assertRaises(Exception, collection.get_int_field, field_name) + self.assertRaises(Exception, collection.get_uint_field, field_name) + eigen_strain_field = collection.get_real_field(field_name) + + print(eigen_strain_field.array.T) + for i, row in enumerate(eigen_strain_field.array.T): + error = np.linalg.norm(i/self.cell.size*eigen_strain - + row.reshape(eigen_strain.shape).T) + self.assertEqual(0, error) + diff --git a/tests/test_field_collections_1.cc b/tests/test_field_collections_1.cc index 44a5988..2706a95 100644 --- a/tests/test_field_collections_1.cc +++ b/tests/test_field_collections_1.cc @@ -1,221 +1,221 @@ /** * @file test_field_collections_1.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test the FieldCollection classes which provide fast optimized iterators * over run-time typed fields * * Copyright © 2017 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 "test_field_collections.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_collection_tests); BOOST_AUTO_TEST_CASE(simple) { constexpr Dim_t sdim = 2; using FC_t = GlobalFieldCollection; FC_t fc; BOOST_CHECK_EQUAL(FC_t::spatial_dim(), sdim); BOOST_CHECK_EQUAL(fc.get_spatial_dim(), sdim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, test_collections, F) { BOOST_CHECK_EQUAL(F::FC_t::spatial_dim(), F::sdim()); BOOST_CHECK_EQUAL(F::fc.get_spatial_dim(), F::sdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_field2_test, F, test_collections, F) { const auto order{2}; using FC_t = typename F::FC_t; using TF_t = TensorField; auto && myfield = make_field("TensorField real 2", F::fc); using TensorMap = TensorFieldMap; using MatrixMap = MatrixFieldMap; using ArrayMap = ArrayFieldMap; TensorMap TFM(myfield); MatrixMap MFM(myfield); ArrayMap AFM(myfield); BOOST_CHECK_EQUAL(TFM.info_string(), "Tensor(d, "+ std::to_string(order) + "_o, " + std::to_string(F::mdim()) + "_d)"); BOOST_CHECK_EQUAL(MFM.info_string(), "Matrix(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); BOOST_CHECK_EQUAL(AFM.info_string(), "Array(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(multi_field_test, F, mult_collections, F) { using FC_t = typename F::FC_t; // possible maptypes for Real tensor fields using T_type = Real; using T_TFM1_t = TensorFieldMap; using T_TFM2_t = TensorFieldMap; //! dangerous using T4_Map_t = T4MatrixFieldMap; // impossible maptypes for Real tensor fields using T_SFM_t = ScalarFieldMap; using T_MFM_t = MatrixFieldMap; using T_AFM_t = ArrayFieldMap; using T_MFMw1_t = MatrixFieldMap; using T_MFMw2_t = MatrixFieldMap; using T_MFMw3_t = MatrixFieldMap; const std::string T_name{"Tensorfield Real o4"}; const std::string T_name_w{"TensorField Real o4 wrongname"}; BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(T_TFM1_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T_TFM2_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T4_Map_t(F::fc.at(T_name))); BOOST_CHECK_THROW(T4_Map_t(F::fc.at(T_name_w)), std::out_of_range); BOOST_CHECK_THROW(T_MFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_AFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw1_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw3_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name_w)), std::out_of_range); // possible maptypes for integer scalar fields using S_type = Int; using S_SFM_t = ScalarFieldMap; using S_TFM1_t = TensorFieldMap; using S_TFM2_t = TensorFieldMap; using S_MFM_t = MatrixFieldMap; using S_AFM_t = ArrayFieldMap; using S4_Map_t = T4MatrixFieldMap; // impossible maptypes for integer scalar fields using S_MFMw1_t = MatrixFieldMap; using S_MFMw2_t = MatrixFieldMap; using S_MFMw3_t = MatrixFieldMap; const std::string S_name{"integer Scalar"}; const std::string S_name_w{"integer Scalar wrongname"}; BOOST_CHECK_NO_THROW(S_SFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM1_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM2_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_MFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_AFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S4_Map_t(F::fc.at(S_name))); BOOST_CHECK_THROW(S_MFMw1_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(T4_Map_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw3_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_SFM_t(F::fc.at(S_name_w)), std::out_of_range); // possible maptypes for complex matrix fields using M_type = Complex; using M_MFM_t = MatrixFieldMap; using M_AFM_t = ArrayFieldMap; // impossible maptypes for complex matrix fields using M_SFM_t = ScalarFieldMap; using M_MFMw1_t = MatrixFieldMap; using M_MFMw2_t = MatrixFieldMap; using M_MFMw3_t = MatrixFieldMap; const std::string M_name{"Matrixfield Complex sdim x mdim"}; const std::string M_name_w{"Matrixfield Complex sdim x mdim wrongname"}; BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(M_MFM_t(F::fc.at(M_name))); BOOST_CHECK_NO_THROW(M_AFM_t(F::fc.at(M_name))); BOOST_CHECK_THROW(M_MFMw1_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw3_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name_w)), std::out_of_range); } /* ---------------------------------------------------------------------- */ //! Check whether fields can be initialized using mult_collections_t = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>>; using mult_collections_f = boost::mpl::list, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_glob, F, mult_collections_t, F) { Ccoord_t size; Ccoord_t loc{}; for (auto && s: size) { s = 3; } BOOST_CHECK_NO_THROW(F::fc.initialise(size, loc)); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca, F, mult_collections_f, F) { testGoodies::RandRange rng; for (int i = 0; i < 7; ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca_with_push_back, F, mult_collections_f, F) { - constexpr auto mdim = F::mdim(); - constexpr int nb_pix = 7; - testGoodies::RandRange rng; + constexpr auto mdim{F::mdim()}; + constexpr int nb_pix{7}; + testGoodies::RandRange rng{}; using ftype = internal::TypedSizedFieldBase< decltype(F::fc), Real, mdim*mdim*mdim*mdim>; using stype = Eigen::Array; - auto & field = reinterpret_cast(F::fc["Tensorfield Real o4"]); + auto & field = static_cast(F::fc["Tensorfield Real o4"]); field.push_back(stype()); for (int i = 0; i < nb_pix; ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_THROW(F::fc.initialise(), FieldCollectionError); for (int i = 0; i < nb_pix-1; ++i) { field.push_back(stype()); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_field_collections_3.cc b/tests/test_field_collections_3.cc index 2ba1ada..02321fe 100644 --- a/tests/test_field_collections_3.cc +++ b/tests/test_field_collections_3.cc @@ -1,143 +1,152 @@ /** * @file test_field_collections_3.cc * * @author Till Junge * * @date 19 Dec 2017 * * @brief Continuation of tests from test_field_collection_2.cc, split for faster * compilation * * Copyright © 2017 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 #include "test_field_collections.hh" #include "tests/test_goodies.hh" #include "common/tensor_algebra.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_collection_tests); /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(assignment_test, Fix, iter_collections, Fix) { auto t4map = Fix::t4_field.get_map(); auto t2map = Fix::t2_field.get_map(); auto scmap = Fix::sc_field.get_map(); auto m2map = Fix::m2_field.get_map(); const auto t4map_val{Matrices::Isymm()}; t4map = t4map_val; const auto t2map_val{Matrices::I2()}; t2map = t2map_val; const Int scmap_val{1}; scmap = scmap_val; Eigen::Matrix m2map_val; m2map_val.setRandom(); m2map = m2map_val; const size_t nb_pts{Fix::fc.size()}; testGoodies::RandRange rnd{}; BOOST_CHECK_EQUAL((t4map[rnd.randval(0, nb_pts-1)] - t4map_val).norm(), 0.); BOOST_CHECK_EQUAL((t2map[rnd.randval(0, nb_pts-1)] - t2map_val).norm(), 0.); BOOST_CHECK_EQUAL((scmap[rnd.randval(0, nb_pts-1)] - scmap_val), 0.); BOOST_CHECK_EQUAL((m2map[rnd.randval(0, nb_pts-1)] - m2map_val).norm(), 0.); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Eigentest, Fix, iter_collections, Fix) { auto t4eigen = Fix::t4_field.eigen(); auto t2eigen = Fix::t2_field.eigen(); BOOST_CHECK_EQUAL(t4eigen.rows(), ipow(Fix::mdim(), 4)); BOOST_CHECK_EQUAL(t4eigen.cols(), Fix::t4_field.size()); using T2_t = typename Eigen::Matrix; T2_t test_mat; test_mat.setRandom(); Eigen::Map> test_map(test_mat.data()); t2eigen.col(0) = test_map; BOOST_CHECK_EQUAL((Fix::t2_field.get_map()[0] - test_mat).norm(), 0.); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(field_proxy_test, Fix, iter_collections, Fix) { Eigen::VectorXd t4values{Fix::t4_field.eigenvec()}; using FieldProxy_t = TypedField; //! create a field proxy FieldProxy_t proxy("proxy to 'Tensorfield Real o4'", Fix::fc, t4values, Fix::t4_field.get_nb_components()); Eigen::VectorXd wrong_size_not_multiple{ Eigen::VectorXd::Zero(t4values.size()+1)}; BOOST_CHECK_THROW(FieldProxy_t("size not a multiple of nb_components", Fix::fc, wrong_size_not_multiple, Fix::t4_field.get_nb_components()), FieldError); Eigen::VectorXd wrong_size_but_multiple{ Eigen::VectorXd::Zero(t4values.size()+ Fix::t4_field.get_nb_components())}; BOOST_CHECK_THROW(FieldProxy_t("size wrong multiple of nb_components", Fix::fc, wrong_size_but_multiple, Fix::t4_field.get_nb_components()), FieldError); using Tensor4Map = T4MatrixFieldMap; Tensor4Map ref_map{Fix::t4_field}; Tensor4Map proxy_map{proxy}; for (auto tup: akantu::zip(ref_map, proxy_map)) { auto & ref = std::get<0>(tup); auto & prox = std::get<1>(tup); BOOST_CHECK_EQUAL((ref-prox).norm(), 0); } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(field_proxy_of_existing_field, Fix, iter_collections, Fix ) { Eigen::Ref t4values{Fix::t4_field.eigenvec()}; using FieldProxy_t = TypedField; //! create a field proxy FieldProxy_t proxy("proxy to 'Tensorfield Real o4'", Fix::fc, t4values, Fix::t4_field.get_nb_components()); using Tensor4Map = T4MatrixFieldMap; Tensor4Map ref_map{Fix::t4_field}; Tensor4Map proxy_map{proxy}; for (auto tup: akantu::zip(ref_map, proxy_map)) { auto & ref = std::get<0>(tup); auto & prox = std::get<1>(tup); prox += prox.Identity(); BOOST_CHECK_EQUAL((ref-prox).norm(), 0); } } + /* ---------------------------------------------------------------------- */ + BOOST_FIXTURE_TEST_CASE_TEMPLATE(typed_field_getter, Fix, + mult_collections, Fix) { + constexpr auto mdim{Fix::mdim()}; + auto & fc{Fix::fc}; + auto & field = fc.template get_typed_field("Tensorfield Real o4"); + BOOST_CHECK_EQUAL(field.get_nb_components(), ipow(mdim, fourthOrder)); + } + BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_material_hyper_elasto_plastic1.cc b/tests/test_material_hyper_elasto_plastic1.cc index 4cc0ef5..f749892 100644 --- a/tests/test_material_hyper_elasto_plastic1.cc +++ b/tests/test_material_hyper_elasto_plastic1.cc @@ -1,343 +1,343 @@ /** * @file test_material_hyper_elasto_plastic1.cc * * @author Till Junge * * @date 25 Feb 2018 * * @brief Tests for the large-strain Simo-type plastic law implemented * using MaterialMuSpectre * * 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 "boost/mpl/list.hpp" #include "materials/material_hyper_elasto_plastic1.hh" #include "materials/materials_toolbox.hh" #include "tests.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_hyper_elasto_plastic_1); template struct MaterialFixture { using Mat = Mat_t; constexpr static Real K{.833}; // bulk modulus constexpr static Real mu{.386}; // shear modulus constexpr static Real H{.004}; // hardening modulus constexpr static Real tau_y0{.003}; // initial yield stress constexpr static Real young{ MatTB::convert_elastic_modulus(K, mu)}; constexpr static Real poisson{ MatTB::convert_elastic_modulus(K, mu)}; MaterialFixture():mat("Name", young, poisson, tau_y0, H){}; constexpr static Dim_t sdim{Mat_t::sdim()}; constexpr static Dim_t mdim{Mat_t::mdim()}; Mat_t mat; }; using mats = boost::mpl::list>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mats, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_stress, Fix, mats, Fix) { // This test uses precomputed reference values (computed using // elasto-plasticity.py) for the 3d case only // need higher tol because of printout precision of reference solutions constexpr Real hi_tol{1e-8}; constexpr Dim_t mdim{Fix::mdim}, sdim{Fix::sdim}; constexpr bool has_precomputed_values{(mdim == sdim) && (mdim == threeD)}; constexpr bool verbose{false}; using Strain_t = Eigen::Matrix; using traits = MaterialMuSpectre_traits>; using LColl_t = typename traits::LFieldColl_t; using StrainStField_t = StateField< TensorField>; using FlowStField_t = StateField< ScalarField>; // using StrainStRef_t = typename traits::LStrainMap_t::reference; // using ScalarStRef_t = typename traits::LScalarMap_t::reference; // create statefields LColl_t coll{}; coll.add_pixel({0}); coll.initialise(); - StrainStField_t F_("previous gradient", coll); - StrainStField_t be_("previous elastic strain", coll); - FlowStField_t eps_("plastic flow", coll); + auto & F_{make_statefield("previous gradient", coll)}; + auto & be_{make_statefield("previous elastic strain", coll)}; + auto & eps_{make_statefield("plastic flow", coll)}; auto F_prev{F_.get_map()}; F_prev[0].current() = Strain_t::Identity(); auto be_prev{be_.get_map()}; be_prev[0].current() = Strain_t::Identity(); auto eps_prev{eps_.get_map()}; eps_prev[0].current() = 0; // elastic deformation Strain_t F{Strain_t::Identity()}; F(0, 1) = 1e-5; F_.cycle(); be_.cycle(); eps_.cycle(); Strain_t stress{Fix::mat.evaluate_stress(F, F_prev[0], be_prev[0], eps_prev[0])}; if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.92999522e-11, 3.86000000e-06, 0.00000000e+00, 3.86000000e-06, -1.93000510e-11, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -2.95741950e-17; Real error{(tau_ref-stress).norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00000000e+00, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00; error = (be_ref-be_prev[0].current()).norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0}; error = ep_ref-eps_prev[0].current(); BOOST_CHECK_LT(error, hi_tol); } if (verbose) { std::cout << "τ =" << std::endl << stress << std::endl; std::cout << "F =" << std::endl << F << std::endl; std::cout << "Fₜ =" << std::endl << F_prev[0].current() << std::endl; std::cout << "bₑ =" << std::endl << be_prev[0].current() << std::endl; std::cout << "εₚ =" << std::endl << eps_prev[0].current() << std::endl; } F_.cycle(); be_.cycle(); eps_.cycle(); // plastic deformation F(0, 1) = .2; stress = Fix::mat.evaluate_stress(F, F_prev[0], be_prev[0], eps_prev[0]); if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.98151335e-04, 1.98151335e-03, 0.00000000e+00, 1.98151335e-03, -1.98151335e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.60615155e-16; Real error{(tau_ref-stress).norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00052666, 0.00513348, 0., 0.00513348, 0.99949996, 0., 0., 0., 1.; error = (be_ref-be_prev[0].current()).norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0.11229988}; error = ep_ref-eps_prev[0].current(); BOOST_CHECK_LT(error, hi_tol); } if (verbose) { std::cout << "Post Cycle" << std::endl; std::cout << "τ =" << std::endl << stress << std::endl << "F =" << std::endl << F << std::endl << "Fₜ =" << std::endl << F_prev[0].current() << std::endl << "bₑ =" << std::endl << be_prev[0].current() << std::endl << "εₚ =" << std::endl << eps_prev[0].current() << std::endl; } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_stiffness, Fix, mats, Fix) { // This test uses precomputed reference values (computed using // elasto-plasticity.py) for the 3d case only // need higher tol because of printout precision of reference solutions constexpr Real hi_tol{1e-8}; constexpr Dim_t mdim{Fix::mdim}, sdim{Fix::sdim}; constexpr bool has_precomputed_values{(mdim == sdim) && (mdim == threeD)}; constexpr bool verbose{has_precomputed_values && false}; using Strain_t = Eigen::Matrix; using Stiffness_t = T4Mat; using traits = MaterialMuSpectre_traits>; using LColl_t = typename traits::LFieldColl_t; using StrainStField_t = StateField< TensorField>; using FlowStField_t = StateField< ScalarField>; // using StrainStRef_t = typename traits::LStrainMap_t::reference; // using ScalarStRef_t = typename traits::LScalarMap_t::reference; // create statefields LColl_t coll{}; coll.add_pixel({0}); coll.initialise(); - StrainStField_t F_("previous gradient", coll); - StrainStField_t be_("previous elastic strain", coll); - FlowStField_t eps_("plastic flow", coll); + auto & F_{make_statefield("previous gradient", coll)}; + auto & be_{make_statefield("previous elastic strain", coll)}; + auto & eps_{make_statefield("plastic flow", coll)}; auto F_prev{F_.get_map()}; F_prev[0].current() = Strain_t::Identity(); auto be_prev{be_.get_map()}; be_prev[0].current() = Strain_t::Identity(); auto eps_prev{eps_.get_map()}; eps_prev[0].current() = 0; // elastic deformation Strain_t F{Strain_t::Identity()}; F(0, 1) = 1e-5; F_.cycle(); be_.cycle(); eps_.cycle(); Strain_t stress{}; Stiffness_t stiffness{}; std::tie(stress, stiffness) = Fix::mat.evaluate_stress_tangent(F, F_prev[0], be_prev[0], eps_prev[0]); if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.92999522e-11, 3.86000000e-06, 0.00000000e+00, 3.86000000e-06, -1.93000510e-11, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -2.95741950e-17; Real error{(tau_ref-stress).norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00000000e+00, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00; error = (be_ref-be_prev[0].current()).norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0}; error = ep_ref-eps_prev[0].current(); BOOST_CHECK_LT(error, hi_tol); Stiffness_t C4_ref{}; C4_ref << 0.67383333, 0., 0., 0., 0.28783333, 0., 0., 0., 0.28783333, 0., 0.193, 0., 0.193, 0., 0., 0., 0., 0., 0., 0., 0.193, 0., 0., 0., 0.193, 0., 0., 0., 0.193, 0., 0.193, 0., 0., 0., 0., 0., 0.28783333, 0., 0., 0., 0.67383333, 0., 0., 0., 0.28783333, 0., 0., 0., 0., 0., 0.193, 0., 0.193, 0., 0., 0., 0.193, 0., 0., 0., 0.193, 0., 0., 0., 0., 0., 0., 0., 0.193, 0., 0.193, 0., 0.28783333, 0., 0., 0., 0.28783333, 0., 0., 0., 0.67383333; error = (C4_ref - stiffness).norm(); BOOST_CHECK_LT(error, hi_tol); } if (verbose) { std::cout << "C₄ =" << std::endl << stiffness << std::endl; } F_.cycle(); be_.cycle(); eps_.cycle(); // plastic deformation F(0, 1) = .2; std::tie(stress, stiffness) = Fix::mat.evaluate_stress_tangent(F, F_prev[0], be_prev[0], eps_prev[0]); if (has_precomputed_values) { Strain_t tau_ref{}; tau_ref << 1.98151335e-04, 1.98151335e-03, 0.00000000e+00, 1.98151335e-03, -1.98151335e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.60615155e-16; Real error{(tau_ref-stress).norm()}; BOOST_CHECK_LT(error, hi_tol); Strain_t be_ref{}; be_ref << 1.00052666, 0.00513348, 0., 0.00513348, 0.99949996, 0., 0., 0., 1.; error = (be_ref-be_prev[0].current()).norm(); BOOST_CHECK_LT(error, hi_tol); Real ep_ref{0.11229988}; error = ep_ref-eps_prev[0].current(); BOOST_CHECK_LT(error, hi_tol); Stiffness_t C4_ref{}; C4_ref << +4.23106224e-01, -4.27959704e-04, 0.00000000e+00, -4.27959704e-04, 4.13218286e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.13175490e-01, -4.27959704e-04, 7.07167743e-04, 0.00000000e+00, 7.07167743e-04, 4.27959704e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.79121029e-18, +0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, -4.27959704e-04, 7.07167743e-04, 0.00000000e+00, 7.07167743e-04, 4.27959704e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.79121029e-18, +4.13218286e-01, 4.27959704e-04, 0.00000000e+00, 4.27959704e-04, 4.23106224e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.13175490e-01, +0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, +0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, +0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, +4.13175490e-01, 2.79121029e-18, 0.00000000e+00, 2.79121029e-18, 4.13175490e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.23149020e-01; error = (C4_ref - stiffness).norm(); BOOST_CHECK_LT(error, hi_tol); } if (verbose) { std::cout << "Post Cycle" << std::endl; std::cout << "C₄ =" << std::endl << stiffness << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_statefields.cc b/tests/test_statefields.cc index 3d7658c..119cf07 100644 --- a/tests/test_statefields.cc +++ b/tests/test_statefields.cc @@ -1,208 +1,298 @@ /** * file test_statefields.cc * * @author Till Junge * * @date 01 Mar 2018 * * @brief Test the StateField abstraction and the associated maps * * @section LICENSE * * 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/field.hh" #include "common/field_collection.hh" #include "common/statefield.hh" #include "common/ccoord_operations.hh" #include "tests.hh" #include #include namespace muSpectre { template struct SF_Fixture { using FC_t = std::conditional_t, LocalFieldCollection>; using Field_t = TensorField; using ScalField_t = ScalarField; constexpr static size_t nb_mem{2}; constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; constexpr static bool global{Global}; - + constexpr static size_t get_nb_mem() {return nb_mem;} + constexpr static Dim_t get_sdim () {return sdim;} + constexpr static Dim_t get_mdim () {return mdim;} + constexpr static bool get_global() {return global;} SF_Fixture() - :fc{}, sf("prefix", fc), scalar_f("scalar", fc), self{*this} {} + :fc{}, + sf{make_statefield>("prefix", fc)}, + scalar_f{make_statefield>("scalar", fc)}, + self{*this} {} FC_t fc; - StateField sf; - StateField scalar_f; + StateField & sf; + StateField & scalar_f; SF_Fixture & self; }; using typelist = boost::mpl::list, SF_Fixture< twoD, threeD, false>, SF_Fixture, SF_Fixture< twoD, twoD, true>, SF_Fixture< twoD, threeD, true>, SF_Fixture>; + BOOST_AUTO_TEST_SUITE(statefield); + + BOOST_AUTO_TEST_CASE(old_values_test) { + constexpr Dim_t Dim{twoD}; + constexpr size_t NbMem{2}; + constexpr bool verbose{false}; + using FC_t = LocalFieldCollection; + FC_t fc{}; + using Field_t = ScalarField; + auto & statefield{make_statefield>("name", fc)}; + fc.add_pixel({}); + fc.initialise(); + for (size_t i{0}; i < NbMem+1; ++i) { + statefield.current().eigen() = i+1; + if (verbose) { + std::cout << "current = " << statefield.current().eigen() << std::endl + << "old 1 = " << statefield.old().eigen() << std::endl + << "old 2 = " << statefield.template old<2>().eigen() + << std::endl + << "indices = " << statefield.get_indices() << std::endl + << std::endl; + } + statefield.cycle(); + } + BOOST_CHECK_EQUAL(statefield.current().eigen()(0), 1); + BOOST_CHECK_EQUAL(statefield.old().eigen()(0), 3); + BOOST_CHECK_EQUAL(statefield.template old<2>().eigen()(0), 2); + + } + BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, Fix, typelist, Fix) { - BOOST_CHECK_EQUAL("prefix", Fix::sf.get_prefix()); + const std::string ref{"prefix"}; + const std::string & fix{Fix::sf.get_prefix()}; + BOOST_CHECK_EQUAL(ref, fix); } namespace internal { template struct init{ static void run(Fixture_t & fix) { constexpr Dim_t dim{std::remove_reference_t::sdim}; fix.fc.initialise(CcoordOps::get_cube(3), CcoordOps::get_cube(0)); } }; template struct init{ static void run(Fixture_t & fix) { constexpr Dim_t dim{std::remove_reference_t::sdim}; CcoordOps::Pixels pixels(CcoordOps::get_cube(3), CcoordOps::get_cube(0)); for (auto && pix: pixels) { fix.fc.add_pixel(pix); } fix.fc.initialise(); } }; } // internal BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_iteration, Fix, typelist, Fix) { internal::init::run(Fix::self); constexpr Dim_t mdim{Fix::mdim}; constexpr bool verbose{false}; using StateFMap = StateFieldMap< MatrixFieldMap, Fix::nb_mem>; StateFMap matrix_map(Fix::sf); for (size_t i = 0; i < Fix::nb_mem+1; ++i) { for (auto && wrapper: matrix_map) { wrapper.current() += (i+1)*wrapper.current().Identity(); if (verbose) { std::cout << "pixel " << wrapper.get_ccoord() << ", memory cycle " << i << std::endl; std::cout << wrapper.current() << std::endl; std::cout << wrapper.old() << std::endl; std::cout << wrapper.template old<2>() << std::endl << std::endl; } } Fix::sf.cycle(); } for (auto && wrapper: matrix_map) { auto I{wrapper.current().Identity()}; Real error{(wrapper.current() - I).norm()}; BOOST_CHECK_LT(error, tol); error = (wrapper.old() - 3*I).norm(); BOOST_CHECK_LT(error, tol); error = (wrapper.template old<2>() - 2* I).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_default_map, Fix, typelist, Fix) { internal::init::run(Fix::self); constexpr bool verbose{false}; auto matrix_map{Fix::sf.get_map()}; for (size_t i = 0; i < Fix::nb_mem+1; ++i) { for (auto && wrapper: matrix_map) { wrapper.current() += (i+1)*wrapper.current().Identity(); if (verbose) { std::cout << "pixel " << wrapper.get_ccoord() << ", memory cycle " << i << std::endl; std::cout << wrapper.current() << std::endl; std::cout << wrapper.old() << std::endl; std::cout << wrapper.template old<2>() << std::endl << std::endl; } } Fix::sf.cycle(); } auto matrix_const_map{Fix::sf.get_const_map()}; for (auto && wrapper: matrix_const_map) { auto I{wrapper.current().Identity()}; Real error{(wrapper.current() - I).norm()}; BOOST_CHECK_LT(error, tol); error = (wrapper.old() - 3*I).norm(); BOOST_CHECK_LT(error, tol); error = (wrapper.template old<2>() - 2* I).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_scalar_map, Fix, typelist, Fix) { internal::init::run(Fix::self); constexpr bool verbose{false}; auto scalar_map{Fix::scalar_f.get_map()}; for (size_t i = 0; i < Fix::nb_mem+1; ++i) { for (auto && wrapper: scalar_map) { wrapper.current() += (i+1); if (verbose) { std::cout << "pixel " << wrapper.get_ccoord() << ", memory cycle " << i << std::endl; std::cout << wrapper.current() << std::endl; std::cout << wrapper.old() << std::endl; std::cout << wrapper.template old<2>() << std::endl << std::endl; } } Fix::scalar_f.cycle(); } auto scalar_const_map{Fix::scalar_f.get_const_map()}; BOOST_CHECK_EQUAL(scalar_const_map[0].current(), scalar_const_map[1].current()); for (auto wrapper: scalar_const_map) { Real error{wrapper.current() - 1}; BOOST_CHECK_LT(error, tol); error = wrapper.old() - 3; BOOST_CHECK_LT(error, tol); error = wrapper.template old<2>() - 2; BOOST_CHECK_LT(error, tol); } } + + /* ---------------------------------------------------------------------- */ + BOOST_FIXTURE_TEST_CASE_TEMPLATE(Polymorphic_access_by_name, Fix, typelist, Fix) { + internal::init::run(Fix::self); + + //constexpr bool verbose{true}; + auto & tensor_field = Fix::fc.get_statefield("prefix"); + BOOST_CHECK_EQUAL(tensor_field.get_nb_memory(), Fix::get_nb_mem()); + + auto & field = Fix::fc.template get_current("prefix"); + BOOST_CHECK_EQUAL(field.get_nb_components(), + ipow(Fix::get_mdim(), secondOrder)); + BOOST_CHECK_THROW(Fix::fc.template get_current("prefix"), std::runtime_error); + auto & old_field = Fix::fc.template get_old("prefix"); + BOOST_CHECK_EQUAL(old_field.get_nb_components(), + field.get_nb_components()); + BOOST_CHECK_THROW(Fix::fc.template get_old("prefix", Fix::get_nb_mem()+1), + std::out_of_range); + + auto & statefield{Fix::fc.get_statefield("prefix")}; + auto & typed_statefield{Fix::fc.template get_typed_statefield("prefix")}; + auto map{ArrayFieldMap + + (typed_statefield.get_current_field())}; + for (auto arr: map) { + arr.setConstant(1); + } + + + + Eigen::ArrayXXd field_copy{field.eigen()}; + statefield.cycle(); + auto & alt_old_field{typed_statefield.get_old_field()}; + + Real err{(field_copy - alt_old_field.eigen()).matrix().norm()/ + field_copy.matrix().norm()}; + BOOST_CHECK_LT(err, tol); + if (not(err