diff --git a/CMakeLists.txt b/CMakeLists.txt index 165a9e5..d709555 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,222 +1,243 @@ # ============================================================================= # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief Main configuration file # # @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. # ============================================================================= cmake_minimum_required(VERSION 3.0.0) project(µSpectre) set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) set(BUILD_SHARED_LIBS ON) set(MUSPECTRE_PYTHON_MAJOR_VERSION 3) add_compile_options(-Wall -Wextra -Weffc++) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) set(MAKE_DOC_TARGET "OFF" CACHE BOOL "If on, a target dev_doc (which builds the documentation) is added") set(MAKE_TESTS "ON" CACHE BOOL "If on, several ctest targets will be built automatically") set(MAKE_EXAMPLES "ON" CACHE BOOL "If on, the executables in the bin folder will be compiled") set(MAKE_BENCHMARKS "ON" CACHE BOOL "If on, the benchmarks will be compiled") set(MPI_PARALLEL "OFF" CACHE BOOL "If on, MPI-parallel solvers become available") set(SPLIT_CELL "OFF" CACHE BOOL "if on, split cell and the cgal intersection calculator is accsessible ") set(RUNNING_IN_CI "OFF" CACHE INTERNAL "changes output format for tests") if(${MAKE_TESTS}) enable_testing() find_package(Boost COMPONENTS unit_test_framework REQUIRED) endif(${MAKE_TESTS}) if(${MPI_PARALLEL}) add_definitions(-DWITH_MPI) find_package(MPI) if (NOT ${MPI_FOUND}) message(SEND_ERROR "You chose MPI but CMake cannot find the MPI package") endif(NOT ${MPI_FOUND}) endif(${MPI_PARALLEL}) include(muspectreTools) +string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type ) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") # using Clang add_compile_options(-Wno-missing-braces) + if ("debug" STREQUAL "${build_type}") + add_compile_options(-O0) + endif() elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") # using GCC add_compile_options(-Wno-non-virtual-dtor) - string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type ) - if ("release" STREQUAL "${build_type}" ) + add_compile_options(-march=native) + if (("relwithdebinfo" STREQUAL "${build_type}") OR ("release" STREQUAL "${build_type}" )) add_compile_options(-march=native) endif() if ("debug" STREQUAL "${build_type}" ) add_compile_options(-O0) endif() elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") # using Intel C++ elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") # using Visual Studio C++ endif() # Do not trust old gcc. the std::optional has memory bugs if(${CMAKE_COMPILER_IS_GNUCC}) if(${CMAKE_CXX_COMPILER_VERSION} VERSION_LESS 6.0.0) add_definitions(-DNO_EXPERIMENTAL) endif() endif() add_external_package(Eigen3 VERSION 3.3.0 CONFIG) add_external_package(pybind11 VERSION 2.2 CONFIG) if(${SPLIT_CELL}) add_external_package(CGAL VERSION 4.12 CONFIG) endif() include_directories( ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR} ) if(APPLE) include_directories(${CMAKE_INSTALL_PREFIX}/include ${Boost_INCLUDE_DIRS}) endif() #build tests (these are before we add -Werror to the compile options) if (${MAKE_TESTS}) ############################################################################## - # build tests + # build library tests file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/test_*.cc") if(${SPLIT_CELL}) file( GLOB SPLIT_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/split_test_*.cc") list (APPEND TEST_SRCS ${SPLIT_TEST_SRCS}) endif(${SPLIT_CELL}) add_executable(main_test_suite tests/main_test_suite.cc ${TEST_SRCS}) target_link_libraries(main_test_suite ${Boost_LIBRARIES} muSpectre mpfr) + muSpectre_add_test(main_test_suite TYPE BOOST main_test_suite --report_level=detailed) + # build header tests + file( GLOB HEADER_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/header_test_*.cc") + foreach(header_test ${HEADER_TEST_SRCS}) + get_filename_component(header_test_name ${header_test} NAME_WE) + string(SUBSTRING ${header_test_name} 12 -1 test_name) + list(APPEND header_tests ${test_name}) + add_executable(${test_name} tests/main_test_suite.cc ${header_test}) + target_link_libraries(${test_name} ${Boost_LIBRARIES} Eigen3::Eigen) + target_include_directories(${test_name} INTERFACE ${muSpectre_INCLUDES}) + muSpectre_add_test(${test_name} TYPE BOOST ${test_name} + --report_level=detailed) + endforeach(header_test ${HEADER_TEST_SRCS}) + + add_custom_target(header_tests) + add_dependencies(header_tests ${header_tests}) + ############################################################################## # copy python test file( GLOB PY_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/python_*.py") foreach(pytest ${PY_TEST_SRCS}) get_filename_component(pytest_name ${pytest} NAME) configure_file( ${pytest} "${CMAKE_BINARY_DIR}/${pytest_name}" COPYONLY) endforeach(pytest ${PY_TEST_SRCS}) find_package(PythonInterp ${MUSPECTRE_PYTHON_MAJOR_VERSION} REQUIRED) muSpectre_add_test(python_binding_test TYPE PYTHON python_binding_tests.py) if(${MPI_PARALLEL}) ############################################################################ # add MPI tests file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/mpi_test_*.cc") add_executable(mpi_main_test_suite tests/mpi_main_test_suite.cc ${TEST_SRCS}) target_link_libraries(mpi_main_test_suite ${Boost_LIBRARIES} muSpectre) muSpectre_add_test(mpi_main_test_suite1 TYPE BOOST MPI_NB_PROCS 1 mpi_main_test_suite --report_level=detailed) muSpectre_add_test(mpi_main_test_suite2 TYPE BOOST MPI_NB_PROCS 2 mpi_main_test_suite --report_level=detailed) muSpectre_add_test(python_mpi_binding_test1 TYPE PYTHON MPI_NB_PROCS 1 python_mpi_binding_tests.py) muSpectre_add_test(python_mpi_binding_test2 TYPE PYTHON MPI_NB_PROCS 2 python_mpi_binding_tests.py) endif(${MPI_PARALLEL}) endif(${MAKE_TESTS}) ################################################################################ # compile the library add_compile_options( -Werror) add_subdirectory( ${CMAKE_SOURCE_DIR}/src/ ) add_subdirectory( ${CMAKE_SOURCE_DIR}/language_bindings/ ) if (${MAKE_DOC_TARGET}) add_subdirectory( ${CMAKE_SOURCE_DIR}/doc/ ) endif() ################################################################################ if (${MAKE_EXAMPLES}) #compile executables set(binaries ${CMAKE_SOURCE_DIR}/bin/demonstrator1.cc ${CMAKE_SOURCE_DIR}/bin/demonstrator_dynamic_solve.cc ${CMAKE_SOURCE_DIR}/bin/demonstrator2.cc ${CMAKE_SOURCE_DIR}/bin/hyper-elasticity.cc ${CMAKE_SOURCE_DIR}/bin/small_case.cc) if (${MPI_PARALLEL}) set (binaries ${binaries} ${CMAKE_SOURCE_DIR}/bin/demonstrator_mpi.cc ) endif (${MPI_PARALLEL}) foreach(binaryfile ${binaries}) get_filename_component(binaryname ${binaryfile} NAME_WE) add_executable(${binaryname} ${binaryfile}) target_link_libraries(${binaryname} ${Boost_LIBRARIES} muSpectre) endforeach(binaryfile ${binaries}) #or copy them file (GLOB pybins "${CMAKE_SOURCE_DIR}/bin/*.py") foreach(pybin ${pybins}) get_filename_component(binaryname ${pybin} NAME_WE) configure_file( ${pybin} "${CMAKE_BINARY_DIR}/${binaryname}.py" COPYONLY) endforeach(pybin ${pybins}) endif (${MAKE_EXAMPLES}) ################################################################################ # compile benchmarks if(${MAKE_BENCHMARKS}) file(GLOB benchmarks "${CMAKE_SOURCE_DIR}/benchmarks/benchmark*cc") foreach(benchmark ${benchmarks}) get_filename_component(benchmark_name ${benchmark} NAME_WE) add_executable(${benchmark_name} ${benchmark}) target_link_libraries(${benchmark_name} ${BOOST_LIBRARIES} muSpectre) endforeach(benchmark ${benchmark}) endif(${MAKE_BENCHMARKS}) diff --git a/src/common/ccoord_operations.hh b/src/common/ccoord_operations.hh index 00ca0c2..f2f9906 100644 --- a/src/common/ccoord_operations.hh +++ b/src/common/ccoord_operations.hh @@ -1,325 +1,324 @@ /** * @file ccoord_operations.hh * * @author Till Junge * * @date 29 Sep 2017 * * @brief common operations on pixel addressing * * 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 #include #include #include #include "common/common.hh" #include "common/iterators.hh" #ifndef CCOORD_OPERATIONS_H #define CCOORD_OPERATIONS_H namespace muSpectre { namespace CcoordOps { namespace internal { //! simple helper returning the first argument and ignoring the second template constexpr T ret(T val, size_t /*dummy*/) {return val;} //! helper to build cubes template constexpr std::array cube_fun(T val, std::index_sequence) { return std::array{ret(val, I)...}; } //! computes hermitian size according to FFTW template constexpr Ccoord_t herm(const Ccoord_t & full_sizes, std::index_sequence) { return Ccoord_t{full_sizes[I]..., full_sizes.back()/2+1}; } //! compute the stride in a direction of a row-major grid template constexpr Dim_t stride(const Ccoord_t & sizes, const size_t index) { static_assert(Dim > 0, "only for positive numbers of dimensions"); auto const diff{Dim - 1 - Dim_t(index)}; Dim_t ret_val{1}; for (Dim_t i{0}; i < diff; ++i) { ret_val *= sizes[Dim-1-i]; } return ret_val; } //! get all strides from a row-major grid (helper function) template constexpr Ccoord_t compute_strides(const Ccoord_t & sizes, std::index_sequence) { return Ccoord_t{stride(sizes, I)...}; } } // internal //----------------------------------------------------------------------------// //! returns a grid of equal resolutions in each direction template constexpr std::array get_cube(T size) { return internal::cube_fun(size, std::make_index_sequence{}); } /* ---------------------------------------------------------------------- */ //! returns the hermition grid to correcsponding to a full grid template constexpr Ccoord_t get_hermitian_sizes(Ccoord_t full_sizes) { return internal::herm(full_sizes, std::make_index_sequence{}); } - /* ---------------------------------------------------------------------- */ //! return physical vector of a cell of cubic pixels template Eigen::Matrix get_vector(const Ccoord_t & ccoord, Real pix_size = 1.) { Eigen::Matrix retval; for (size_t i = 0; i < dim; ++i) { retval[i] = pix_size * ccoord[i]; } return retval; } /* ---------------------------------------------------------------------- */ //! return physical vector of a cell of general pixels template Eigen::Matrix get_vector(const Ccoord_t & ccoord, Eigen::Matrix pix_size) { Eigen::Matrix retval = pix_size; for (size_t i = 0; i < dim; ++i) { retval[i] *= ccoord[i]; } return retval; } /* ---------------------------------------------------------------------- */ //! return physical vector of a cell of general pixels template Eigen::Matrix get_vector(const Ccoord_t & ccoord, const std::array& pix_size) { Eigen::Matrix retval{}; for (size_t i = 0; i < dim; ++i) { retval[i] = pix_size[i]*ccoord[i]; } return retval; } /* ---------------------------------------------------------------------- */ //! get all strides from a row-major grid template constexpr Ccoord_t get_default_strides(const Ccoord_t & sizes) { return internal::compute_strides(sizes, std::make_index_sequence{}); } //----------------------------------------------------------------------------// //! get the i-th pixel in a grid of size sizes template constexpr Ccoord_t get_ccoord(const Ccoord_t & resolutions, const Ccoord_t & locations, Dim_t index) { Ccoord_t retval{{0}}; Dim_t factor{1}; for (Dim_t i = dim-1; i >=0; --i) { retval[i] = index/factor%resolutions[i] + locations[i]; if (i != 0 ) { factor *= resolutions[i]; } } return retval; } //----------------------------------------------------------------------------// //! get the i-th pixel in a grid of size sizes template constexpr Ccoord_t get_ccoord(const Ccoord_t & resolutions, const Ccoord_t & locations, Dim_t index, std::index_sequence) { Ccoord_t ccoord{get_ccoord(resolutions, locations, index)}; return Ccoord_t({ccoord[I]...}); } //----------------------------------------------------------------------------// //! get the linear index of a pixel in a given grid template constexpr Dim_t get_index(const Ccoord_t & sizes, const Ccoord_t & locations, const Ccoord_t & ccoord) { Dim_t retval{0}; Dim_t factor{1}; for (Dim_t i = dim-1; i >=0; --i) { retval += (ccoord[i]-locations[i])*factor; if (i != 0) { factor *= sizes[i]; } } return retval; } //----------------------------------------------------------------------------// //! get the linear index of a pixel given a set of strides template constexpr Dim_t get_index_from_strides(const Ccoord_t & strides, const Ccoord_t & ccoord) { Dim_t retval{0}; for (const auto & tup: akantu::zip(strides, ccoord)) { const auto & stride = std::get<0>(tup); const auto & ccord_ = std::get<1>(tup); retval += stride * ccord_; } return retval; } //----------------------------------------------------------------------------// //! get the number of pixels in a grid template constexpr size_t get_size(const Ccoord_t& sizes) { Dim_t retval{1}; for (size_t i = 0; i < dim; ++i) { retval *= sizes[i]; } return retval; } //----------------------------------------------------------------------------// //! get the number of pixels in a grid given its strides template constexpr size_t get_size_from_strides(const Ccoord_t& sizes, const Ccoord_t& strides) { return sizes[0]*strides[0]; } /* ---------------------------------------------------------------------- */ /** * centralises iterating over square (or cubic) discretisation * grids. The optional parameter pack `dmap` can be used to * specify the order of the axes in which to iterate over the * dimensions (i.e., dmap = 0, 1, 2 is rowmajor, and 0, 2, 1 would * be a custom order in which the second and third dimension are * transposed */ template class Pixels { public: //! constructor Pixels(const Ccoord_t & resolutions=Ccoord_t{}, const Ccoord_t & locations=Ccoord_t{}) :resolutions{resolutions}, locations{locations}{}; //! copy constructor Pixels(const Pixels & other) = default; //! assignment operator Pixels & operator=(const Pixels & other) = default; virtual ~Pixels() = default; /** * iterators over `Pixels` dereferences to cell coordinates */ class iterator { public: using value_type = Ccoord_t; //!< stl conformance using const_value_type = const value_type; //!< stl conformance using pointer = value_type*; //!< stl conformance using difference_type = std::ptrdiff_t; //!< stl conformance using iterator_category = std::forward_iterator_tag;//!resolutions);} protected: Ccoord_t resolutions; //!< resolutions of this domain Ccoord_t locations; //!< locations of this domain }; /* ---------------------------------------------------------------------- */ template Pixels::iterator::iterator(const Pixels & pixels, bool begin) :pixels{pixels}, index{begin? 0: get_size(pixels.resolutions)} {} /* ---------------------------------------------------------------------- */ template typename Pixels::iterator::value_type Pixels::iterator::operator*() const { return get_ccoord(pixels.resolutions, pixels.locations, this->index, std::conditional_t, std::index_sequence>{}); } /* ---------------------------------------------------------------------- */ template bool Pixels::iterator::operator!=(const iterator &other) const { return (this->index != other.index) || (&this->pixels != &other.pixels); } /* ---------------------------------------------------------------------- */ template bool Pixels::iterator::operator==(const iterator &other) const { return !(*this!= other); } /* ---------------------------------------------------------------------- */ template typename Pixels::iterator& Pixels::iterator::operator++() { ++this->index; return *this; } } // CcoordOps } // muSpectre #endif /* CCOORD_OPERATIONS_H */ diff --git a/src/common/field.hh b/src/common/field.hh index af4a289..4f45c3c 100644 --- a/src/common/field.hh +++ b/src/common/field.hh @@ -1,662 +1,678 @@ /** * @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 "common/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); + template + inline void push_back(const Eigen::DenseBase & 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: + 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 - { + ipow(dim,order)> { 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(); + inline 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(); + inline 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; + inline decltype(auto) get_map() const; + + /** + * creates a `TensorField` same size and type as this, but all + * entries are zero. Convenience function + */ + inline TensorField & get_zeros_like(std::string unique_name) 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(); + inline 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(); + inline 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; + inline decltype(auto) get_map() const; + /** + * creates a `MatrixField` same size and type as this, but all + * entries are zero. Convenience function + */ + inline MatrixField & get_zeros_like(std::string unique_name) 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 + template void TypedSizedFieldBase:: - push_back(const Stored_t & value) { + push_back(const Eigen::DenseBase & value) { + static_assert(Derived::SizeAtCompileTime == NbComponents, + "You provided an array with the wrong number of entries."); + static_assert((Derived::RowsAtCompileTime == 1) or + (Derived::ColsAtCompileTime == 1), + "You have not provided a column or row vector."); static_assert (not FieldCollection::Global, "You can only push_back data into local field " - "collections"); + "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 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 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 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 + auto TensorField:: + get_zeros_like(std::string unique_name) const -> TensorField & { + return make_field(unique_name, + this->collection); + } + /* ---------------------------------------------------------------------- */ template 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 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 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); } + /* ---------------------------------------------------------------------- */ + template + auto MatrixField:: + get_zeros_like(std::string unique_name) const -> MatrixField & { + return make_field(unique_name, + this->collection); + } + + } // muSpectre #endif /* FIELD_H */ diff --git a/src/common/field_base.hh b/src/common/field_base.hh index 95f825e..752bb14 100644 --- a/src/common/field_base.hh +++ b/src/common/field_base.hh @@ -1,210 +1,210 @@ /** * file field_base.hh * * @author Till Junge * * @date 10 Apr 2018 * * @brief Virtual base class for fields * * 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 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 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; + 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 1e8e333..7c85e6c 100644 --- a/src/common/field_collection_base.hh +++ b/src/common/field_collection_base.hh @@ -1,351 +1,355 @@ /** * @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_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_t& operator[](std::string unique_name); //! retrieve field by unique_name with bounds checking 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) { 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. " + err_str << "a field named '" << field->get_name() + << "' is already registered in this field collection. " << "Currently registered fields: "; + std::string prelude{""}; for (const auto& name_field_pair: this->fields) { - err_str << ", " << name_field_pair.first; + err_str << prelude << '\'' << name_field_pair.first << '\''; + prelude = ", "; } 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. " + err_str << "a state field named '" << field->get_prefix() + << "' is already registered in this field collection. " << "Currently registered fields: "; + std::string prelude{""}; for (const auto& name_field_pair: this->statefields) { - err_str << ", " << name_field_pair.first; + err_str << prelude << '\'' << name_field_pair.first << '\''; + prelude = ", "; } 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 auto FieldCollectionBase:: operator[](std::string unique_name) -> Field_t & { return *(this->fields[unique_name]); } /* ---------------------------------------------------------------------- */ template auto FieldCollectionBase:: 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_collection_global.hh b/src/common/field_collection_global.hh index 7252c9b..76058c1 100644 --- a/src/common/field_collection_global.hh +++ b/src/common/field_collection_global.hh @@ -1,208 +1,220 @@ /** * @file field_collection_global.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief FieldCollection base-class for global 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_COLLECTION_GLOBAL_H #define FIELD_COLLECTION_GLOBAL_H #include "common/field_collection_base.hh" #include "common/ccoord_operations.hh" namespace muSpectre { - /* ---------------------------------------------------------------------- */ + /** + * forward declaration + */ + template + class LocalFieldCollection; + /** `GlobalFieldCollection` derives from `FieldCollectionBase` and stores * global fields that live throughout the whole computational domain, i.e. * are defined for each pixel. */ - template + template class GlobalFieldCollection: public FieldCollectionBase> { public: //! for compile time check constexpr static bool Global{true}; using Parent = FieldCollectionBase >; //!< base class + //! helpful for functions that fill global fields from local fields + using LocalFieldCollection_t = LocalFieldCollection; + //! helpful for functions that fill global fields from local fields + using GlobalFieldCollection_t = GlobalFieldCollection; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type using Field_p = typename Parent::Field_p; //!< spatial coordinates type //! iterator over all pixels contained it the collection using iterator = typename CcoordOps::Pixels::iterator; //! Default constructor GlobalFieldCollection(); //! Copy constructor GlobalFieldCollection(const GlobalFieldCollection &other) = delete; //! Move constructor GlobalFieldCollection (GlobalFieldCollection &&other) = default; //! Destructor virtual ~GlobalFieldCollection() = default; //! Copy assignment operator GlobalFieldCollection& operator=(const GlobalFieldCollection &other) = delete; //! Move assignment operator GlobalFieldCollection& operator=(GlobalFieldCollection &&other) = default; /** allocate memory, etc. At this point, the collection is informed aboud the size and shape of the domain (through the sizes parameter). The job of initialise is to make sure that all fields are either of size 0, in which case they need to be allocated, or are of the same size as the product of 'sizes' (if standard strides apply) any field of a different size is wrong. TODO: check whether it makes sense to put a runtime check here **/ inline void initialise(Ccoord sizes, Ccoord locations); //! return subdomain resolutions inline const Ccoord & get_sizes() const; //! return subdomain locations inline const Ccoord & get_locations() const; //! returns the linear index corresponding to cell coordinates template inline size_t get_index(CcoordRef && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; - inline iterator begin(); //!< returns iterator to first pixel - inline iterator end(); //!< returns iterator past the last pixel + inline iterator begin() const; //!< returns iterator to first pixel + inline iterator end() const; //!< returns iterator past the last pixel //! return spatial dimension (template parameter) static constexpr inline Dim_t spatial_dim() {return DimS;} + + //! return globalness at compile time + static constexpr inline bool is_global() {return Global;} protected: //! number of discretisation cells in each of the DimS spatial directions Ccoord sizes{}; Ccoord locations{}; CcoordOps::Pixels pixels{}; //!< helper to iterate over the grid private: }; /* ---------------------------------------------------------------------- */ template GlobalFieldCollection::GlobalFieldCollection() :Parent() {} /* ---------------------------------------------------------------------- */ template void GlobalFieldCollection:: initialise(Ccoord sizes, Ccoord locations) { if (this->is_initialised) { throw std::runtime_error("double initialisation"); } this->pixels = CcoordOps::Pixels(sizes, locations); this->size_ = CcoordOps::get_size(sizes); this->sizes = sizes; this->locations = locations; std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { auto && field = *item.second; const auto field_size = field.size(); if (field_size == 0) { field.resize(this->size()); } else if (field_size != this->size()) { std::stringstream err_stream; err_stream << "Field '" << field.get_name() << "' contains " << field_size << " entries, but the field collection " << "has " << this->size() << " pixels"; throw FieldCollectionError(err_stream.str()); } }); this->is_initialised = true; } //----------------------------------------------------------------------------// //! return subdomain resolutions template const typename GlobalFieldCollection::Ccoord & GlobalFieldCollection::get_sizes() const { return this->sizes; } //----------------------------------------------------------------------------// //! return subdomain locations template const typename GlobalFieldCollection::Ccoord & GlobalFieldCollection::get_locations() const { return this->locations; } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename GlobalFieldCollection::Ccoord GlobalFieldCollection::get_ccoord(size_t index) const { return CcoordOps::get_ccoord(this->get_sizes(), this->get_locations(), std::move(index)); } /* ---------------------------------------------------------------------- */ template typename GlobalFieldCollection::iterator - GlobalFieldCollection::begin() { + GlobalFieldCollection::begin() const { return this->pixels.begin(); } /* ---------------------------------------------------------------------- */ template typename GlobalFieldCollection::iterator - GlobalFieldCollection::end() { + GlobalFieldCollection::end() const { return this->pixels.end(); } //----------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template template size_t GlobalFieldCollection::get_index(CcoordRef && ccoord) const { static_assert(std::is_same< Ccoord, std::remove_const_t< std::remove_reference_t>>::value, "can only be called with values or references of Ccoord"); return CcoordOps::get_index(this->get_sizes(), this->get_locations(), std::forward(ccoord)); } } // muSpectre #endif /* FIELD_COLLECTION_GLOBAL_H */ diff --git a/src/common/field_collection_local.hh b/src/common/field_collection_local.hh index e8d22b3..a2085c0 100644 --- a/src/common/field_collection_local.hh +++ b/src/common/field_collection_local.hh @@ -1,184 +1,203 @@ /** * @file field_collection_local.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief FieldCollection base-class for local 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_COLLECTION_LOCAL_H #define FIELD_COLLECTION_LOCAL_H #include "common/field_collection_base.hh" namespace muSpectre { - /* ---------------------------------------------------------------------- */ + /** + * forward declaration + */ + template + class GlobalFieldCollection; + /** `LocalFieldCollection` derives from `FieldCollectionBase` and stores * local fields, i.e. fields that are only defined for a subset of all pixels * in the computational domain. The coordinates of these active pixels are * explicitly stored by this field collection. * `LocalFieldCollection::add_pixel` allows to add individual pixels to the * field collection. */ - template + template class LocalFieldCollection: public FieldCollectionBase> { public: //! for compile time check constexpr static bool Global{false}; //! base class using Parent = FieldCollectionBase>; + //! helpful for functions that fill local fields from global fields + using GlobalFieldCollection_t = GlobalFieldCollection; + //! helpful for functions that fill local fields from global fields + using LocalFieldCollection_t = LocalFieldCollection; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type using Field_p = typename Parent::Field_p; //!< field pointer using ccoords_container = std::vector; //!< list of pixels //! iterator over managed pixels using iterator = typename ccoords_container::iterator; + //! const iterator over managed pixels + using const_iterator = typename ccoords_container::const_iterator; //! Default constructor LocalFieldCollection(); //! Copy constructor LocalFieldCollection(const LocalFieldCollection &other) = delete; //! Move constructor LocalFieldCollection(LocalFieldCollection &&other) = delete; //! Destructor virtual ~LocalFieldCollection() = default; //! Copy assignment operator LocalFieldCollection& operator=(const LocalFieldCollection &other) = delete; //! Move assignment operator LocalFieldCollection& operator=(LocalFieldCollection &&other) = delete; //! add a pixel/voxel to the field collection inline void add_pixel(const Ccoord & local_ccoord); /** allocate memory, etc. at this point, the field collection knows how many entries it should have from the size of the coords containes (which grows by one every time add_pixel is called. The job of initialise is to make sure that all fields are either of size 0, in which case they need to be allocated, or are of the same size as the product of 'sizes' any field of a different size is wrong TODO: check whether it makes sense to put a runtime check here **/ inline void initialise(); //! returns the linear index corresponding to cell coordinates template inline size_t get_index(CcoordRef && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; //! iterator to first pixel inline iterator begin() {return this->ccoords.begin();} //! iterator past last pixel inline iterator end() {return this->ccoords.end();} + //! const iterator to first pixel + inline const_iterator begin() const {return this->ccoords.cbegin();} + //! const iterator past last pixel + inline const_iterator end() const {return this->ccoords.cend();} + + + //! return globalness at compile time + static constexpr inline bool is_global() {return Global;} protected: //! container of pixel coords for non-global collections ccoords_container ccoords{}; //! container of indices for non-global collections (slow!) std::map indices{}; // std::vector assigned_ratio{}; private: }; /* ---------------------------------------------------------------------- */ template LocalFieldCollection::LocalFieldCollection() :Parent(){} /* ---------------------------------------------------------------------- */ template void LocalFieldCollection:: add_pixel(const Ccoord & local_ccoord) { if (this->is_initialised) { throw FieldCollectionError ("once a field collection has been initialised, you can't add new " "pixels."); } this->indices[local_ccoord] = this->ccoords.size(); this->ccoords.push_back(local_ccoord); this->size_++; } /* ---------------------------------------------------------------------- */ template void LocalFieldCollection:: initialise() { if (this->is_initialised) { throw std::runtime_error("double initialisation"); } std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { auto && field = *item.second; const auto field_size = field.size(); if (field_size == 0) { field.resize(this->size()); } else if (field_size != this->size()) { std::stringstream err_stream; err_stream << "Field '" << field.get_name() << "' contains " << field_size << " entries, but the field collection " << "has " << this->size() << " pixels"; throw FieldCollectionError(err_stream.str()); } }); this->is_initialised = true; } //----------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template template size_t LocalFieldCollection::get_index(CcoordRef && ccoord) const { static_assert(std::is_same< Ccoord, std::remove_const_t< std::remove_reference_t>>::value, "can only be called with values or references of Ccoord"); return this->indices.at(std::forward(ccoord)); } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename LocalFieldCollection::Ccoord LocalFieldCollection::get_ccoord(size_t index) const { return this->ccoords[std::move(index)]; } } // muSpectre #endif /* FIELD_COLLECTION_LOCAL_H */ diff --git a/src/common/field_helpers.hh b/src/common/field_helpers.hh new file mode 100644 index 0000000..7bad374 --- /dev/null +++ b/src/common/field_helpers.hh @@ -0,0 +1,58 @@ +/** + * file field_helpers.hh + * + * @author Till Junge + * + * @date 30 Aug 2018 + * + * @brief helper functions that needed to be sequestered to avoid circular + * inclusions + * + * @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_HELPERS_H +#define FIELD_HELPERS_H + +#include + +namespace muSpectre { + + /** + * 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; + } + + + +} // muSpectre + +#endif /* FIELD_HELPERS_H */ diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh index da1388b..66c49bd 100644 --- a/src/common/field_map_base.hh +++ b/src/common/field_map_base.hh @@ -1,645 +1,870 @@ /** * @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 { - //----------------------------------------------------------------------------// + /** + * Forward-declaration + */ + template + class TypedSizedFieldBase; + + //! 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 { + static_assert((NbComponents != 0), + "Fields with now components make no sense."); + /* + * Eigen::Dynamic is equal to -1, and is a legal value, hence + * the following peculiar check + */ + static_assert((NbComponents > -2), + "Fields with a negative number of components make no sense."); 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 + //! non-constant version of field + using TypedField_nc = std::conditional_t<(NbComponents >= 1), + TypedSizedFieldBase, + TypedField>; //! field type as seen from iterator - using TypedField = std::conditional_t; - using Field = typename TypedField::Base; //!< iterated field type + using TypedField_t = std::conditional_t; + using Field = typename TypedField_nc::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; - } + class iterator; - 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: - }; + /** + * Simple iterable proxy wrapper object around a FieldMap. When + * iterated over, rather than dereferencing to the reference + * type of iterator, it dereferences to a tuple of the pixel, + * and the reference type of iterator + */ + template + class enumerator; - TypedField & get_field() {return this->field;} + TypedField_t & 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 + TypedField_t & field; //!< mapped Field private: }; + + + /** + * iterates over all pixels in the `muSpectre::FieldCollection` + * and dereferences to an Eigen map to the currently used field. + */ + template + template + class FieldMap::iterator + { + static_assert(!((ConstIter==false) && (ConstField==true)), + "You can't have a non-const iterator over a const " + "field"); + public: + //! for use by enumerator + using FullyTypedFieldMap_t = FullyTypedFieldMap; + //! 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); + + //! 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 debugging) + 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: + //! Copy constructor + iterator(const iterator &other) = default; + + const FieldCollection & collection; //!< collection of the field + TypedRef fieldmap; //!< ref to the field itself + size_t index; //!< index of currently pointed-to pixel + private: + }; + + + /* ---------------------------------------------------------------------- */ + template + template + class FieldMap::enumerator + { + public: + //! fully typed reference as seen by the iterator + using TypedRef =typename Iterator::TypedRef; + + //! Default constructor + enumerator() = delete; + + //! constructor with field mapped + enumerator(TypedRef& field_map): field_map{field_map} {} + + /** + * similar to iterators of the field map, but dereferences to a + * tuple containing the cell coordinates and teh corresponding + * entry + */ + class iterator; + + iterator begin() {return iterator(this->field_map);} + iterator end() {return iterator(this->field_map, false);} + + protected: + TypedRef & field_map; + }; + + /* ---------------------------------------------------------------------- */ + template + template + class FieldMap:: + enumerator::iterator { + public: + //! cell coordinates type + using Ccoord = typename FieldCollection::Ccoord; + //! stl conformance + using value_type = std::tuple; + //! stl conformance + using const_value_type = std::tuple; + //! stl conformance + using difference_type = std::ptrdiff_t; + //! stl conformance + using iterator_category = std::random_access_iterator_tag; + //! stl conformance + using reference = std::tuple< + Ccoord, + typename SimpleIterator::FullyTypedFieldMap_t::reference>; + + //! Default constructor + iterator() = delete; + + //! constructor for begin/end + iterator(TypedRef fieldmap, bool begin=true): it{fieldmap, begin} {} + + //! constructor for random access + iterator(TypedRef fieldmap, size_t index): it{fieldmap, index} {} + + //! constructor from iterator + iterator(const SimpleIterator & it): it{it} {} + + //! 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 + iterator & operator++() { + ++(this->it); + return *this; + } + + //! post-increment + iterator operator++(int) { + iterator current = *this; + ++(this->it); + return current; + } + + //! dereference + value_type operator*() { + return value_type{it.get_ccoord(), *it}; + } + + //! dereference + const_value_type operator*() const { + return const_value_type{it.get_ccoord(), *it}; + }; + + //! pre-decrement + iterator & operator--() { + --(this->it); + return *this; + } + + //! post-decrement + iterator operator--(int) { + iterator current = *this; + --(this->it); + return current; + } + + //! access subscripting + value_type operator[](difference_type diff) { + SimpleIterator accessed{this->it + diff}; + return *accessed; + } + + //! access subscripting + const_value_type operator[](const difference_type diff) const { + SimpleIterator accessed{this->it + diff}; + return *accessed; + } + + //! equality + bool operator==(const iterator & other) const { + return this->it == other.it; + } + + //! inequality + bool operator!=(const iterator & other) const { + return this->it != other.it; + } + + //! div. comparisons + bool operator<(const iterator & other) const { + return this->it < other.it; + } + + //! div. comparisons + bool operator<=(const iterator & other) const { + return this->it <= other.it; + } + + //! div. comparisons + bool operator>(const iterator & other) const { + return this->it > other.it; + } + + //! div. comparisons + bool operator>=(const iterator & other) const { + return this->it >= other.it; + } + + //! additions, subtractions and corresponding assignments + iterator operator+(difference_type diff) const { + return iterator{this->it + diff}; + } + + //! additions, subtractions and corresponding assignments + iterator operator-(difference_type diff) const { + return iterator{this->it - diff}; + } + + //! additions, subtractions and corresponding assignments + iterator& operator+=(difference_type diff) { + this->it += diff; + } + + //! additions, subtractions and corresponding assignments + iterator& operator-=(difference_type diff) { + this->it -= diff; + } + + protected: + SimpleIterator it; + private: + }; + + + } // internal namespace internal { /* ---------------------------------------------------------------------- */ template FieldMap:: FieldMap(Field_c& field) - :collection(field.get_collection()), field(static_cast(field)) { - static_assert(NbComponents>0, + :collection(field.get_collection()), field(static_cast(field)) { + static_assert((NbComponents > 0) or (NbComponents == Eigen::Dynamic), "Only fields with more than 0 components allowed"); } /* ---------------------------------------------------------------------- */ template template FieldMap:: FieldMap(TypedSizedFieldBase & field) - :collection(field.get_collection()), field(static_cast(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() != + 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() + "'"}; + 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()) { + if ((NbComponents != Dim_t(this->field.get_nb_components())) and + (NbComponents != Eigen::Dynamic)) { throw FieldInterpretationError ("Cannot create a Map of type '" + - this->info_string() + - "' for field '" + this->field.get_name() + "' with " + + 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), + static_assert((TypedField_t::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}; + 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_dynamic.hh b/src/common/field_map_dynamic.hh new file mode 100644 index 0000000..a875519 --- /dev/null +++ b/src/common/field_map_dynamic.hh @@ -0,0 +1,217 @@ +/** + * @file field_map_dynamic.hh + * + * @author Till Junge + * + * @date 24 Jul 2018 + * + * @brief Field map for dynamically sized maps. for use in non-critical + * applications (i.e., i/o, postprocessing, etc, but *not* in a hot + * loop + * + * 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_MAP_DYNAMIC_H +#define FIELD_MAP_DYNAMIC_H + +#include "common/field_map_base.hh" + +namespace muSpectre { + + /** + * Maps onto any `muSpectre::TypedField` and lets you iterate over + * it in the form of `Eigen::Map`. This is significantly slower than the statically sized field + * maps and should only be used in non-critical contexts. + */ + template + class TypedFieldMap: public internal::FieldMap + { + public: + //! base class + using Parent = internal::FieldMap; + //! sister class with all params equal, but ConstField guaranteed true + using ConstMap = TypedFieldMap; + //! cell coordinates type + using Ccoord = Ccoord_t; + //! plain Eigen type + using Arr_t = Eigen::Array; + using value_type = Eigen::Map; //!< stl conformance + using const_reference = Eigen::Map; //!< 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 + + //! polymorphic base field type (for debug and python) + using Field = typename Parent::Field; + //! polymorphic base field type (for debug and python) + using Field_c = typename Parent::Field_c; + //! stl conformance + using const_iterator = typename Parent::template iterator; + //! stl conformance + using iterator = std::conditional_t< + ConstField, + const_iterator, + typename Parent::template iterator>; + //! stl conformance + using reverse_iterator = std::reverse_iterator; + //! stl conformance + using const_reverse_iterator = std::reverse_iterator; + //! enumerator over a constant scalar field + using const_enumerator = typename Parent::template enumerator; + //! enumerator over a scalar field + using enumerator = std::conditional_t< + ConstField, + const_enumerator, + typename Parent::template enumerator>; + //! give access to the protected fields + friend iterator; + + + //! Default constructor + TypedFieldMap() = delete; + + //! Constructor + TypedFieldMap(Field_c & field); + + //! Copy constructor + TypedFieldMap(const TypedFieldMap & other) = delete; + + //! Move constructor + TypedFieldMap(TypedFieldMap &&other) = default; + + //! Destructor + virtual ~TypedFieldMap() = default; + + //! Copy assignment operator + TypedFieldMap& operator=(const TypedFieldMap &other) = delete; + + //! Assign a matrixlike value to every entry + template + inline TypedFieldMap & operator=(const Eigen::EigenBase & val); + + //! Move assignment operator + TypedFieldMap& operator=(TypedFieldMap &&other) = default; + + //! 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 first entry of field + inline iterator begin(){return iterator(*this);} + //! return an iterator to first entry of field + inline const_iterator cbegin() const {return const_iterator(*this);} + //! return an iterator to first entry of field + inline const_iterator begin() const {return this->cbegin();} + //! return an iterator past the last entry of field + inline iterator end(){return iterator(*this, false);} + //! return an iterator past the last entry of field + inline const_iterator cend() const {return const_iterator(*this, false);} + //! return an iterator past the last entry of field + inline const_iterator end() const {return this->cend();} + + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + enumerator enumerate() {return enumerator(*this);} + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + const_enumerator enumerate() const {return const_enumerator(*this);} + + //! evaluate the average of the field + inline Arr_t mean() const; + + protected: + //! for sad, legacy iterator use + inline pointer ptr_to_val_t(size_type index); + private: + }; + + //----------------------------------------------------------------------------// + template + TypedFieldMap:: + TypedFieldMap(Field_c & field): + Parent(field) { + this->check_compatibility(); + } + + + //----------------------------------------------------------------------------// + template + std::string TypedFieldMap:: + info_string() const { + std::stringstream info; + info << "Dynamic(" << typeid(T).name() << ", " + << this->field.get_nb_components() + << ")"; + return info.str(); + } + + //----------------------------------------------------------------------------// + template + auto TypedFieldMap:: + operator[](size_type index) -> reference { + return reference{this->get_ptr_to_entry(index), + Dim_t(this->field.get_nb_components())}; + } + + //----------------------------------------------------------------------------// + template + auto TypedFieldMap:: + operator[](const Ccoord & ccoord) -> reference { + size_t index{this->collection.get_index(ccoord)}; + return (*this)[index]; + } + + //----------------------------------------------------------------------------// + template + auto TypedFieldMap:: + operator[](size_type index) const -> const_reference { + return const_reference{this->get_ptr_to_entry(index), + Dim_t(this->field.get_nb_components())}; + } + + //----------------------------------------------------------------------------// + template + auto TypedFieldMap:: + operator[](const Ccoord & ccoord) const -> const_reference { + size_t index{this->collection.get_index(ccoord)}; + return (*this)[index]; + } + +} // muSpectre + +#endif /* FIELD_MAP_DYNAMIC_H */ diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh index e69d18b..7f730d4 100644 --- a/src/common/field_map_matrixlike.hh +++ b/src/common/field_map_matrixlike.hh @@ -1,376 +1,392 @@ /** * @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 + //! enumerator over a constant scalar field + using const_enumerator = typename Parent::template enumerator; + //! enumerator over a scalar field + using enumerator = std::conditional_t< + ConstField, + const_enumerator, + typename Parent::template enumerator>; + //! give access to the protected fields 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; + MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = delete; - //! Move constructor + //! Move constructorxo 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);} + inline iterator begin(){return std::move(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();} + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + enumerator enumerate() {return enumerator(*this);} + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + const_enumerator enumerate() const {return const_enumerator(*this);} + //! 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() { return NameWrapper::field_info_root();} //!< for printing and debug private: }; /* ---------------------------------------------------------------------- */ 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 << 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))); + size_t && index{this->collection.get_index(ccoord)}; + return reference(this->get_ptr_to_entry(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))); + size_t && index{this->collection.get_index(ccoord)}; + return const_reference(this->get_ptr_to_entry(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_map_scalar.hh b/src/common/field_map_scalar.hh index 4ebaf7f..70baffa 100644 --- a/src/common/field_map_scalar.hh +++ b/src/common/field_map_scalar.hh @@ -1,217 +1,235 @@ /** * @file field_map_scalar.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief maps over scalar 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_SCALAR_H #define FIELD_MAP_SCALAR_H #include "common/field_map_base.hh" namespace muSpectre { /** * Implements maps on scalar fields (i.e. material properties, * temperatures, etc). Maps onto a `muSpectre::internal::TypedSizedFieldBase` * and lets you iterate over it in the form of the bare type of the field. */ template class ScalarFieldMap : public internal::FieldMap { public: //! base class using Parent = internal::FieldMap; //! sister class with all params equal, but ConstField guaranteed true using ConstMap = ScalarFieldMap; //! cell coordinates type using Ccoord = Ccoord_t; using value_type = T; //!< stl conformance using const_reference = const value_type &; //!< stl conformance //! stl conformance using reference = std::conditional_t; using size_type = typename Parent::size_type; //!< stl conformance using pointer = T*; //!< 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; //! iterator over a scalar field using iterator = std::conditional_t >; using reverse_iterator = std::reverse_iterator; //!< stl conformance //! stl conformance using const_reverse_iterator = std::reverse_iterator; + //! enumerator over a constant scalar field + using const_enumerator = typename Parent::template enumerator; + //! enumerator over a scalar field + using enumerator = std::conditional_t< + ConstField, + const_enumerator, + typename Parent::template enumerator>; //! give access to the protected fields friend iterator; //! Default constructor ScalarFieldMap() = delete; //! constructor ScalarFieldMap(Field_c & field); //! Copy constructor ScalarFieldMap(const ScalarFieldMap &other) = default; //! Move constructor ScalarFieldMap(ScalarFieldMap &&other) = default; //! Destructor virtual ~ScalarFieldMap() = default; //! Copy assignment operator ScalarFieldMap& operator=(const ScalarFieldMap &other) = delete; //! Move assignment operator ScalarFieldMap& operator=(ScalarFieldMap &&other) = delete; //! Assign a value to every entry ScalarFieldMap& operator=(T val); //! give human-readable field map type inline std::string info_string() const override final; //! member access inline reference operator[](size_type index); inline reference operator[](const Ccoord& ccoord); inline const_reference operator[] (size_type index) const; inline const_reference operator[] (const Ccoord& ccoord) const; //! return an iterator to the first pixel of the field - inline iterator begin(){return iterator(*this);} + iterator begin(){return iterator(*this);} //! return an iterator to the first pixel of the field - inline const_iterator cbegin() const {return const_iterator(*this);} + const_iterator cbegin() const {return const_iterator(*this);} //! return an iterator to the first pixel of the field - inline const_iterator begin() const {return this->cbegin();} + const_iterator begin() const {return this->cbegin();} //! return an iterator to tail of field for ranges - inline iterator end(){return iterator(*this, false);} + 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);} + 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();} + const_iterator end() const {return this->cend();} + + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + enumerator enumerate() {return enumerator(*this);} + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + const_enumerator enumerate() const {return const_enumerator(*this);} //! evaluate the average of the field inline T mean() const; protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); //! type identifier for printing and debugging const static std::string field_info_root; private: }; /* ---------------------------------------------------------------------- */ template ScalarFieldMap:: ScalarFieldMap(Field_c & field) :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string ScalarFieldMap::info_string() const { std::stringstream info; info << "Scalar(" << typeid(T).name() << ")"; return info.str(); } /* ---------------------------------------------------------------------- */ template const std::string ScalarFieldMap::field_info_root{ "Scalar"}; /* ---------------------------------------------------------------------- */ //! member access template typename ScalarFieldMap::reference ScalarFieldMap::operator[](size_type index) { return this->get_ptr_to_entry(std::move(index))[0]; } /* ---------------------------------------------------------------------- */ //! member access template typename ScalarFieldMap::reference ScalarFieldMap::operator[](const Ccoord& ccoord) { auto && index = this->collection.get_index(std::move(ccoord)); return this->get_ptr_to_entry(std::move(index))[0]; } /* ---------------------------------------------------------------------- */ //! member access template typename ScalarFieldMap::const_reference ScalarFieldMap:: operator[](size_type index) const { return this->get_ptr_to_entry(std::move(index))[0]; } /* ---------------------------------------------------------------------- */ //! member access template typename ScalarFieldMap::const_reference ScalarFieldMap:: operator[](const Ccoord& ccoord) const { auto && index = this->collection.get_index(std::move(ccoord)); return this->get_ptr_to_entry(std::move(index))[0]; } /* ---------------------------------------------------------------------- */ //! Assign a value to every entry template ScalarFieldMap & ScalarFieldMap::operator=(T val) { for (auto & scalar:*this) { scalar = val; } return *this; } /* ---------------------------------------------------------------------- */ template T ScalarFieldMap:: mean() const { T mean{0}; for (auto && val: *this) { mean += val; } mean /= Real(this->size()); return mean; } } // muSpectre #endif /* FIELD_MAP_SCALAR_H */ diff --git a/src/common/field_map_tensor.hh b/src/common/field_map_tensor.hh index 6898dda..b823444 100644 --- a/src/common/field_map_tensor.hh +++ b/src/common/field_map_tensor.hh @@ -1,269 +1,287 @@ /** * @file field_map_tensor.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief Defines an Eigen-Tensor map 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_TENSOR_H #define FIELD_MAP_TENSOR_H #include "common/eigen_tools.hh" #include "common/field_map_base.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ /** * Maps onto a `muSpectre::internal::TypedSizedFieldBase` and lets * you iterate over it in the form of `Eigen::TensorMap>` */ template class TensorFieldMap: public internal::FieldMap::Sizes::total_size, ConstField> { public: //! base class using Parent = internal::FieldMap; //! sister class with all params equal, but ConstField guaranteed true using ConstMap = TensorFieldMap; //! cell coordinates type using Ccoord = Ccoord_t; //! static tensor size using Sizes = typename SizesByOrder::Sizes; //! plain Eigen type using T_t = Eigen::TensorFixedSize; using value_type = Eigen::TensorMap; //!< stl conformance using const_reference = Eigen::TensorMap; //!< 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; //!< field to map + //! polymorphic base field type (for debug and python) using Field = typename Parent::Field; //! polymorphic base field type (for debug and python) using Field_c = typename Parent::Field_c; //! stl conformance using const_iterator = typename Parent::template iterator; //! stl conformance using iterator = std::conditional_t< ConstField, const_iterator, typename Parent::template iterator>; //! stl conformance using reverse_iterator = std::reverse_iterator; //! stl conformance using const_reverse_iterator = std::reverse_iterator; - //! stl conformance + //! enumerator over a constant scalar field + using const_enumerator = typename Parent::template enumerator; + //! enumerator over a scalar field + using enumerator = std::conditional_t< + ConstField, + const_enumerator, + typename Parent::template enumerator>; + //! give access to the protected fields friend iterator; //! Default constructor TensorFieldMap() = delete; //! constructor TensorFieldMap(Field_c & field); //! Copy constructor - TensorFieldMap(const TensorFieldMap &other) = default; + TensorFieldMap(const TensorFieldMap &other) = delete; //! Move constructor TensorFieldMap(TensorFieldMap &&other) = default; //! Destructor virtual ~TensorFieldMap() = default; //! Copy assignment operator TensorFieldMap& operator=(const TensorFieldMap &other) = delete; //! Assign a matrixlike value to every entry inline TensorFieldMap & operator=(const T_t & val); //! Move assignment operator TensorFieldMap& operator=(TensorFieldMap &&other) = delete; //! 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 first entry of field inline iterator begin(){return iterator(*this);} //! return an iterator to first entry of field inline const_iterator cbegin() const {return const_iterator(*this);} //! return an iterator to first entry of field inline const_iterator begin() const {return this->cbegin();} //! return an iterator past the last entry of field inline iterator end(){return iterator(*this, false);} //! return an iterator past the last entry of field inline const_iterator cend() const {return const_iterator(*this, false);} //! return an iterator past the last entry of field inline const_iterator end() const {return this->cend();} + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + enumerator enumerate() {return enumerator(*this);} + /** + * return an iterable proxy to this field that can be iterated + * in Ccoord-value tuples + */ + const_enumerator enumerate() const {return const_enumerator(*this);} + //! 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); private: }; /* ---------------------------------------------------------------------- */ template TensorFieldMap:: TensorFieldMap(Field_c & field) :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string TensorFieldMap::info_string() const { std::stringstream info; info << "Tensor(" << typeid(T).name() << ", " << order << "_o, " << dim << "_d)"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template typename TensorFieldMap::reference TensorFieldMap:: operator[](size_type index) { auto && lambda = [this, &index](auto&&...tens_sizes) { return reference(this->get_ptr_to_entry(index), tens_sizes...); }; return call_sizes(lambda); } template typename TensorFieldMap::reference TensorFieldMap:: operator[](const Ccoord & ccoord) { auto && index = this->collection.get_index(ccoord); auto && lambda = [this, &index](auto&&...sizes) { return reference(this->get_ptr_to_entry(index), sizes...); }; return call_sizes(lambda); } template typename TensorFieldMap:: const_reference TensorFieldMap:: operator[](size_type index) const { // Warning: due to a inconsistency in Eigen's API, tensor maps // cannot be constructed from a const ptr, hence this nasty const // cast :( auto && lambda = [this, &index](auto&&...tens_sizes) { return const_reference(const_cast(this->get_ptr_to_entry(index)), tens_sizes...); }; return call_sizes(lambda); } template typename TensorFieldMap:: const_reference TensorFieldMap:: operator[](const Ccoord & ccoord) const { auto && index = this->collection.get_index(ccoord); auto && lambda = [this, &index](auto&&...sizes) { return const_reference(const_cast(this->get_ptr_to_entry(index)), sizes...); }; return call_sizes(lambda); } /* ---------------------------------------------------------------------- */ template TensorFieldMap & TensorFieldMap:: operator=(const T_t & val) { for (auto && tens: *this) { tens = val; } return *this; } /* ---------------------------------------------------------------------- */ template typename TensorFieldMap::T_t TensorFieldMap:: mean() const { T_t mean{T_t::Zero()}; for (auto && val: *this) { mean += val; } mean *= 1./Real(this->size()); return mean; } /* ---------------------------------------------------------------------- */ //! for sad, legacy iterator use. Don't use unless you have to. template typename TensorFieldMap::pointer TensorFieldMap:: ptr_to_val_t(size_type index) { auto && lambda = [this, &index](auto&&... tens_sizes) { return std::make_unique(this->get_ptr_to_entry(index), tens_sizes...); }; return call_sizes(lambda); } } // muSpectre #endif /* FIELD_MAP_TENSOR_H */ diff --git a/src/common/field_typed.hh b/src/common/field_typed.hh index e7f3687..55f7d4f 100644 --- a/src/common/field_typed.hh +++ b/src/common/field_typed.hh @@ -1,342 +1,507 @@ /** * file field_typed.hh * * @author Till Junge * * @date 10 Apr 2018 * - * @brief Typed Field for dynamically sized fields and base class for fields + * @brief Typed Field for dynamically sized fields and base class for fields * of tensors, matrices, etc * * 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 "common/field_base.hh" +#include "common/field_helpers.hh" #include namespace muSpectre { + /** + * forward-declaration + */ + template + class TypedFieldMap; + + namespace internal { + + /* ---------------------------------------------------------------------- */ + //! declaraton for friending + template + class FieldMap; + + } // internal + + /** * 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 { + friend class internal::FieldMap; + friend class internal::FieldMap; + + static constexpr bool Global{FieldCollection::is_global()}; + public: using Parent = internal::FieldBase; //!< base class //! for type checks when mapping this field using collection_t = typename Parent::collection_t; + + //! for filling global fields from local fields and vice-versa + using LocalField_t = + std::conditional_t, + TypedField>; + //! for filling global fields from local fields and vice-versa + using GlobalField_t = + std::conditional_t>; + using Scalar = T; //!< for type checks using Base = Parent; //!< for uniformity of interface //! Plain Eigen type to map using EigenRep_t = Eigen::Array; 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 accessing entire field using EigenVecConst_t = Eigen::Map; + //! associated non-const field map + using FieldMap_t = TypedFieldMap; + //! associated const field map + using ConstFieldMap_t = TypedFieldMap; /** * type stored (unfortunately, we can't statically size the second * 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); + template + inline void push_back(const Eigen::DenseBase & 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; - + /** + * Convenience function to return a map onto this field. A map + * allows iteration over all pixels. The map's iterator returns a + * dynamically sized `Eigen::Map` the data associated with a + * pixel. + */ + inline FieldMap_t get_map(); + + /** + * Convenience function to return a map onto this field. A map + * allows iteration over all pixels. The map's iterator returns a + * dynamically sized `Eigen::Map` the data associated with a + * pixel. + */ + inline ConstFieldMap_t get_map() const; + + /** + * Convenience function to return a map onto this field. A map + * allows iteration over all pixels. The map's iterator returns a + * dynamically sized `Eigen::Map` the data associated with a + * pixel. + */ + inline ConstFieldMap_t get_const_map() const; + + + /** + * creates a `TypedField` same size and type as this, but all + * entries are zero. Convenience function + */ + inline TypedField & get_zeros_like(std::string unique_name) const; + + /** + * Fill the content of the local field into the global field + * (obviously only for pixels that actually are present in the + * local field) + */ + template + inline std::enable_if_t + fill_from_local(const LocalField_t & local); + + /** + * For pixels that are present in the local field, fill them with + * the content of the global field at those pixels + */ + template + inline std::enable_if_t + fill_from_global(const GlobalField_t & global); 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: }; +} // muSpectre + +#include "common/field_map_dynamic.hh" + +namespace muSpectre { + /* ---------------------------------------------------------------------- */ /* 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(), 1); } /* ---------------------------------------------------------------------- */ template auto TypedField:: eigenvec() const -> EigenVecConst_t { 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 + auto TypedField::get_map() -> FieldMap_t { + return FieldMap_t(*this); + } + + /* ---------------------------------------------------------------------- */ + template + auto TypedField::get_map() const -> ConstFieldMap_t { + return ConstFieldMap_t(*this); + } + + /* ---------------------------------------------------------------------- */ + template + auto TypedField::get_const_map() const -> ConstFieldMap_t { + return ConstFieldMap_t(*this); + } + + /* ---------------------------------------------------------------------- */ + template + auto TypedField:: + get_zeros_like(std::string unique_name) const -> TypedField& { + return make_field(unique_name, + this->collection, + this->nb_components); + } + + /* ---------------------------------------------------------------------- */ + template + template + std::enable_if_t TypedField:: + fill_from_local(const LocalField_t & local) { + static_assert(IsGlobal == Global, "SFINAE parameter, do not touch"); + if (not (local.get_nb_components() == this->get_nb_components())) { + std::stringstream err_str{}; + err_str << "Fields not compatible: You are trying to write a local " + << local.get_nb_components() << "-component field into a global " + << this->get_nb_components() << "-component field."; + throw std::runtime_error(err_str.str()); + } + auto this_map{this->get_map()}; + for (const auto && key_val: local.get_map().enumerate()) { + const auto & key{std::get<0>(key_val)}; + const auto & value{std::get<1>(key_val)}; + this_map[key] = value; + } + } + + /* ---------------------------------------------------------------------- */ + template + template + std::enable_if_t TypedField:: + fill_from_global(const GlobalField_t & global) { + static_assert(IsLocal == not Global, "SFINAE parameter, do not touch"); + if (not (global.get_nb_components() == this->get_nb_components())) { + std::stringstream err_str{}; + err_str << "Fields not compatible: You are trying to write a global " + << global.get_nb_components() << "-component field into a local " + << this->get_nb_components() << "-component field."; + throw std::runtime_error(err_str.str()); + } + + auto global_map{global.get_map()}; + + for (auto && key_val: this->get_map().enumerate()) { + const auto & key{std::get<0>(key_val)}; + auto & value{std::get<1>(key_val)}; + value = global_map[key]; + } + } + /* ---------------------------------------------------------------------- */ 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 + template void TypedField:: - push_back(const Stored_t & value) { + push_back(const Eigen::DenseBase & 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 460e272..9d2101b 100644 --- a/src/common/statefield.hh +++ b/src/common/statefield.hh @@ -1,661 +1,669 @@ /** * 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_helpers.hh" #include "common/field.hh" #include "common/utilities.hh" #include #include #include namespace muSpectre { + /** + * Forward-declaration + */ + template + class TypedField; + + /** * Base class for state fields, useful for storing polymorphic references */ 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 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 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 TypedSizedStateField { public: //! the underlying field's collection type 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; //! 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 Field_t& current() { return this->fields[this->indices[0]]; } //! get (constant) previous field template 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.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(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 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() override final { for (auto & val: this->indices) { val = (val+1)%(nb_memory+1); } } protected: /** * 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(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_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, "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_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); } //! 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/tests/test_ccoord_operations.cc b/tests/header_test_ccoord_operations.cc similarity index 100% rename from tests/test_ccoord_operations.cc rename to tests/header_test_ccoord_operations.cc diff --git a/tests/test_eigen_tools.cc b/tests/header_test_eigen_tools.cc similarity index 97% rename from tests/test_eigen_tools.cc rename to tests/header_test_eigen_tools.cc index dc8481c..7685fee 100644 --- a/tests/test_eigen_tools.cc +++ b/tests/header_test_eigen_tools.cc @@ -1,59 +1,59 @@ /** - * @file test_eigen_tools.cc + * @file header_test_eigen_tools.cc * * @author Till Junge * * @date 07 Mar 2018 * * @brief test the eigen_tools * * 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/eigen_tools.hh" #include "tests.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(eigen_tools); BOOST_AUTO_TEST_CASE(exponential_test) { using Mat_t = Eigen::Matrix; Mat_t input{}; input << 0, .25*pi, 0, .25*pi, 0, 0, 0, 0, 1; Mat_t output{}; output << 1.32460909, 0.86867096, 0, 0.86867096, 1.32460909, 0, 0, 0, 2.71828183; auto my_output{expm(input)}; Real error{(my_output-output).norm()}; BOOST_CHECK_LT(error, 1e-8); if (error >= 1e-8) { std::cout << "input:" << std::endl << input << std::endl; std::cout << "output:" << std::endl << output << std::endl; std::cout << "my_output:" << std::endl << my_output << std::endl; } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/header_test_field_collections.cc b/tests/header_test_field_collections.cc new file mode 100644 index 0000000..62e1b5d --- /dev/null +++ b/tests/header_test_field_collections.cc @@ -0,0 +1,693 @@ +/** + * @file header_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" +#include "common/field_map_dynamic.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{}; + using ftype = internal::TypedSizedFieldBase< + decltype(F::fc), Real, mdim*mdim*mdim*mdim>; + using stype = Eigen::Array; + 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_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) { + using FC_t = typename F::Parent::FC_t; + using Tensor4Map = TensorFieldMap; + Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; + TypedFieldMap dyn_map{F::fc["Tensorfield Real o4"]}; + F::fc["Tensorfield Real o4"].set_zero(); + + for (auto && tens:T4map) { + BOOST_CHECK_EQUAL(Real(Eigen::Tensor(tens.abs().sum().eval())()), 0); + } + for (auto && tens: T4map) { + tens.setRandom(); + } + + for (auto && tup: akantu::zip(T4map, dyn_map) ) { + auto & tens = std::get<0>(tup); + auto & dyn = std::get<1>(tup); + constexpr Dim_t nb_comp{ipow(F::mdim(), order)}; + Eigen::Map> tens_arr(tens.data()); + Real error{(dyn-tens_arr).matrix().norm()}; + BOOST_CHECK_EQUAL(error, 0); + } + + using Tensor2Map = TensorFieldMap; + using MSqMap = MatrixFieldMap; + using ASqMap = ArrayFieldMap; + using A2Map = ArrayFieldMap; + using WrongMap = ArrayFieldMap; + Tensor2Map T2map{F::fc["Tensorfield Real o2"]}; + MSqMap Mmap{F::fc["Tensorfield Real o2"]}; + ASqMap Amap{F::fc["Tensorfield Real o2"]}; + A2Map DynMap{F::fc["Dynamically sized Field"]}; + auto & fc_ref{F::fc}; + BOOST_CHECK_THROW(WrongMap{fc_ref["Dynamically sized Field"]}, + FieldInterpretationError); + auto t2_it = T2map.begin(); + auto t2_it_end = T2map.end(); + auto m_it = Mmap.begin(); + auto a_it = Amap.begin(); + for (; t2_it != t2_it_end; ++t2_it, ++m_it, ++a_it) { + t2_it->setRandom(); + auto && m = *m_it; + bool comp = (m == a_it->matrix()); + BOOST_CHECK(comp); + } + + size_t counter{0}; + for (auto val: DynMap) { + ++counter; + val += val.Ones()*counter; + } + + counter = 0; + for (auto val: DynMap) { + ++counter; + val -= val.Ones()*counter; + auto error {val.matrix().norm()}; + BOOST_CHECK_LT(error, tol); + } + + + using ScalarMap = ScalarFieldMap; + ScalarMap s_map{F::fc["integer Scalar"]}; + for (Uint i = 0; i < s_map.size(); ++i) { + s_map[i] = i; + } + + counter = 0; + for (const auto& val: s_map) { + BOOST_CHECK_EQUAL(counter++, val); + } + + } + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(ccoord_indexing_test, F, glob_iter_colls, F) { + using FC_t = typename F::Parent::FC_t; + using ScalarMap = ScalarFieldMap; + ScalarMap s_map{F::fc["integer Scalar"]}; + for (Uint i = 0; i < s_map.size(); ++i) { + s_map[i] = i; + } + + + for (size_t i = 0; i < CcoordOps::get_size(F::fc.get_sizes()); ++i) { + BOOST_CHECK_EQUAL(CcoordOps::get_index(F::fc.get_sizes(), + F::fc.get_locations(), + CcoordOps::get_ccoord + (F::fc.get_sizes(), + F::fc.get_locations(), + i)), i); + } + } + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(iterator_methods_test, F, iter_collections, F) { + using FC_t = typename F::Parent::FC_t; + using Tensor4Map = TensorFieldMap; + Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; + using it_t = typename Tensor4Map::iterator; + std::ptrdiff_t diff{3}; // arbitrary, as long as it is smaller than the container size + + // check constructors + auto itstart = T4map.begin(); // standard way of obtaining iterator + auto itend = T4map.end(); // ditto + + it_t it1{T4map}; + it_t it2{T4map, false}; + it_t it3{T4map, size_t(diff)}; + BOOST_CHECK(itstart == itstart); + BOOST_CHECK(itstart != itend); + BOOST_CHECK_EQUAL(itstart, it1); + BOOST_CHECK_EQUAL(itend, it2); + + // check ostream operator + std::stringstream response; + response << it3; + BOOST_CHECK_EQUAL + (response.str(), + std::string + ("iterator on field 'Tensorfield Real o4', entry ") + + std::to_string(diff)); + + // check move and assigment constructor (and operator+) + it_t it_repl{T4map}; + it_t itmove = std::move(T4map.begin()); + it_t it4 = it1+diff; + it_t it7 = it4 -diff; + //BOOST_CHECK_EQUAL(itcopy, it1); + BOOST_CHECK_EQUAL(itmove, it1); + BOOST_CHECK_EQUAL(it4, it3); + BOOST_CHECK_EQUAL(it7, it1); + + // check increments/decrements + BOOST_CHECK_EQUAL(it1++, it_repl); // post-increment + BOOST_CHECK_EQUAL(it1, it_repl+1); + BOOST_CHECK_EQUAL(--it1, it_repl); // pre -decrement + BOOST_CHECK_EQUAL(++it1, it_repl+1); // pre -increment + BOOST_CHECK_EQUAL(it1--, it_repl+1); // post-decrement + BOOST_CHECK_EQUAL(it1, it_repl); + + + // dereference and member-of-pointer check + Eigen::Tensor Tens = *it1; + Eigen::Tensor Tens2 = *itstart; + Eigen::Tensor check = (Tens==Tens2).all(); + BOOST_CHECK_EQUAL(bool(check()), true); + + BOOST_CHECK_NO_THROW(itstart->setZero()); + + //check access subscripting + auto T3a = *it3; + auto T3b = itstart[diff]; + BOOST_CHECK(bool(Eigen::Tensor((T3a==T3b).all())())); + + // div. comparisons + BOOST_CHECK_LT(itstart, itend); + BOOST_CHECK(!(itend < itstart)); + BOOST_CHECK(!(itstart < itstart)); + + BOOST_CHECK_LE(itstart, itend); + BOOST_CHECK_LE(itstart, itstart); + BOOST_CHECK(!(itend <= itstart)); + + BOOST_CHECK_GT(itend, itstart); + BOOST_CHECK(!(itend>itend)); + BOOST_CHECK(!(itstart>itend)); + + BOOST_CHECK_GE(itend, itstart); + BOOST_CHECK_GE(itend, itend); + BOOST_CHECK(!(itstart >= itend)); + + // check assignment increment/decrement + BOOST_CHECK_EQUAL(it1+=diff, it3); + BOOST_CHECK_EQUAL(it1-=diff, itstart); + + // check cell coordinates + using Ccoord = Ccoord_t; + Ccoord a{itstart.get_ccoord()}; + Ccoord b{0}; + + // Weirdly, boost::has_left_shift::value is false for Ccoords, even though the operator is implemented :( + //BOOST_CHECK_EQUAL(a, b); + bool check2 = (a==b); + BOOST_CHECK(check2); + } + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_tensor_iter_test, F, iter_collections, F) { + using FC_t = typename F::Parent::FC_t; + using Tensor4Map = TensorFieldMap; + Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; + + using T_t = typename Tensor4Map::T_t; + Eigen::TensorMap Tens2(T4map[0].data(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim()); + + for (auto it = T4map.cbegin(); it != T4map.cend(); ++it) { + // maps to const tensors can't be initialised with a const pointer this sucks + auto&& tens = *it; + auto&& ptr = tens.data(); + + static_assert(std::is_pointer>::value, "should be getting a pointer"); + //static_assert(std::is_const>::value, "should be const"); + + // If Tensor were written well, above static_assert should pass, and the + // following check shouldn't. If it get triggered, it means that a newer + // version of Eigen now does have const-correct + // TensorMap. This means that const-correct field maps + // are then also possible for tensors + BOOST_CHECK(!std::is_const>::value); + + } + } + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_matrix_iter_test, F, iter_collections, F) { + using FC_t = typename F::Parent::FC_t; + using MatrixMap = MatrixFieldMap; + MatrixMap Mmap{F::fc["Matrixfield Complex sdim x mdim"]}; + + for (auto it = Mmap.cbegin(); it != Mmap.cend(); ++it) { + // maps to const tensors can't be initialised with a const pointer this sucks + auto&& mat = *it; + auto&& ptr = mat.data(); + + static_assert(std::is_pointer>::value, "should be getting a pointer"); + //static_assert(std::is_const>::value, "should be const"); + + // If Matrix were written well, above static_assert should pass, and the + // following check shouldn't. If it get triggered, it means that a newer + // version of Eigen now does have const-correct + // MatrixMap. This means that const-correct field maps + // are then also possible for matrices + BOOST_CHECK(!std::is_const>::value); + + } + } + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_scalar_iter_test, F, iter_collections, F) { + using FC_t = typename F::Parent::FC_t; + using ScalarMap = ScalarFieldMap; + ScalarMap Smap{F::fc["integer Scalar"]}; + + for (auto it = Smap.cbegin(); it != Smap.cend(); ++it) { + auto&& scal = *it; + static_assert(std::is_const>::value, + "referred type should be const"); + static_assert(std::is_lvalue_reference::value, + "Should have returned an lvalue ref"); + + } + } + /* ---------------------------------------------------------------------- */ + 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()}; + auto dymap{Fix::dyn_field.get_map()}; + + auto t4map_c{Fix::t4_field.get_const_map()}; + auto t2map_c{Fix::t2_field.get_const_map()}; + auto scmap_c{Fix::sc_field.get_const_map()}; + auto m2map_c{Fix::m2_field.get_const_map()}; + auto dymap_c{Fix::dyn_field.get_const_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_FIXTURE_TEST_CASE_TEMPLATE(enumeration, 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()}; + auto dymap{Fix::dyn_field.get_map()}; + + for (auto && tup: akantu::zip(scmap.get_collection(), scmap, + scmap.enumerate())) { + + const auto & ccoord_ref = std::get<0>(tup); + const auto & val_ref = std::get<1>(tup); + const auto & key_val = std::get<2>(tup); + const auto & ccoord = std::get<0>(key_val); + const auto & val = std::get<1>(key_val); + + for (auto && ccoords: akantu::zip(ccoord_ref, ccoord)) { + const auto & ref{std::get<0>(ccoords)}; + const auto & val{std::get<1>(ccoords)}; + BOOST_CHECK_EQUAL(ref, val); + } + + const auto error{std::abs(val-val_ref)}; + BOOST_CHECK_EQUAL(error, 0); + } + + for (auto && tup: akantu::zip(t4map.get_collection(), t4map, + t4map.enumerate())) { + + const auto & ccoord_ref = std::get<0>(tup); + const auto & val_ref = std::get<1>(tup); + const auto & key_val = std::get<2>(tup); + const auto & ccoord = std::get<0>(key_val); + const auto & val = std::get<1>(key_val); + + for (auto && ccoords: akantu::zip(ccoord_ref, ccoord)) { + const auto & ref{std::get<0>(ccoords)}; + const auto & val{std::get<1>(ccoords)}; + BOOST_CHECK_EQUAL(ref, val); + } + + const auto error{(val-val_ref).norm()}; + BOOST_CHECK_EQUAL(error, 0); + } + + for (auto && tup: akantu::zip(t2map.get_collection(), t2map, + t2map.enumerate())) { + + const auto & ccoord_ref = std::get<0>(tup); + const auto & val_ref = std::get<1>(tup); + const auto & key_val = std::get<2>(tup); + const auto & ccoord = std::get<0>(key_val); + const auto & val = std::get<1>(key_val); + + for (auto && ccoords: akantu::zip(ccoord_ref, ccoord)) { + const auto & ref{std::get<0>(ccoords)}; + const auto & val{std::get<1>(ccoords)}; + BOOST_CHECK_EQUAL(ref, val); + } + + const auto error{(val-val_ref).norm()}; + BOOST_CHECK_EQUAL(error, 0); + } + + for (auto && tup: akantu::zip(m2map.get_collection(), m2map, + m2map.enumerate())) { + + const auto & ccoord_ref = std::get<0>(tup); + const auto & val_ref = std::get<1>(tup); + const auto & key_val = std::get<2>(tup); + const auto & ccoord = std::get<0>(key_val); + const auto & val = std::get<1>(key_val); + + for (auto && ccoords: akantu::zip(ccoord_ref, ccoord)) { + const auto & ref{std::get<0>(ccoords)}; + const auto & val{std::get<1>(ccoords)}; + BOOST_CHECK_EQUAL(ref, val); + } + + const auto error{(val-val_ref).norm()}; + BOOST_CHECK_EQUAL(error, 0); + } + + for (auto && tup: akantu::zip(dymap.get_collection(), dymap, + dymap.enumerate())) { + + const auto & ccoord_ref = std::get<0>(tup); + const auto & val_ref = std::get<1>(tup); + const auto & key_val = std::get<2>(tup); + const auto & ccoord = std::get<0>(key_val); + const auto & val = std::get<1>(key_val); + + for (auto && ccoords: akantu::zip(ccoord_ref, ccoord)) { + const auto & ref{std::get<0>(ccoords)}; + const auto & val{std::get<1>(ccoords)}; + BOOST_CHECK_EQUAL(ref, val); + } + + const auto error{(val-val_ref).matrix().norm()}; + BOOST_CHECK_EQUAL(error, 0); + } + } + + + + BOOST_AUTO_TEST_SUITE_END(); + +} // muSpectre diff --git a/tests/test_fields.cc b/tests/header_test_fields.cc similarity index 67% rename from tests/test_fields.cc rename to tests/header_test_fields.cc index 0f985c3..2c4d0b4 100644 --- a/tests/test_fields.cc +++ b/tests/header_test_fields.cc @@ -1,190 +1,264 @@ /** - * @file test_fields.cc + * @file header_test_fields.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test Fields that are used in FieldCollections * * 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 "tests.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common/ccoord_operations.hh" #include #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_test); template struct FieldFixture { constexpr static bool IsGlobal{Global}; constexpr static Dim_t Order{secondOrder}; constexpr static Dim_t SDim{twoD}; constexpr static Dim_t MDim{threeD}; constexpr static Dim_t NbComponents{ipow(MDim, Order)}; using FieldColl_t = std::conditional_t, LocalFieldCollection>; using TField_t = TensorField; + using MField_t = MatrixField; using DField_t = TypedField; FieldFixture() : tensor_field{make_field("TensorField", this->fc)}, + matrix_field{make_field("MatrixField", this->fc)}, dynamic_field1{ make_field("Dynamically sized field with correct number of" " components", this->fc, ipow(MDim, Order))}, dynamic_field2{ make_field("Dynamically sized field with incorrect number" - "of components", this->fc, NbComponents+1)} + " of components", this->fc, NbComponents+1)} {} ~FieldFixture() = default; FieldColl_t fc{}; TField_t & tensor_field; + MField_t & matrix_field; DField_t & dynamic_field1; DField_t & dynamic_field2; }; using field_fixtures = boost::mpl::list, FieldFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE(size_check_global, FieldFixture) { // check that fields are initialised with empty vector BOOST_CHECK_EQUAL(tensor_field.size(), 0); BOOST_CHECK_EQUAL(dynamic_field1.size(), 0); BOOST_CHECK_EQUAL(dynamic_field2.size(), 0); // check that returned size is correct Dim_t len{2}; auto sizes{CcoordOps::get_cube(len)}; fc.initialise(sizes, {}); auto nb_pixels{CcoordOps::get_size(sizes)}; BOOST_CHECK_EQUAL(tensor_field.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field1.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field2.size(), nb_pixels); constexpr Dim_t pad_size{3}; tensor_field.set_pad_size(pad_size); dynamic_field1.set_pad_size(pad_size); dynamic_field2.set_pad_size(pad_size); // check that setting pad size won't change logical size BOOST_CHECK_EQUAL(tensor_field.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field1.size(), nb_pixels); BOOST_CHECK_EQUAL(dynamic_field2.size(), nb_pixels); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE(size_check_local, FieldFixture) { // check that fields are initialised with empty vector BOOST_CHECK_EQUAL(tensor_field.size(), 0); BOOST_CHECK_EQUAL(dynamic_field1.size(), 0); BOOST_CHECK_EQUAL(dynamic_field2.size(), 0); // check that returned size is correct Dim_t nb_pixels{3}; Eigen::Array new_elem; Eigen::Array wrong_elem; for (Dim_t i{0}; i; FC_t fc; using TF_t = TensorField; auto & field{make_field("TensorField 1", fc)}; // check that fields are initialised with empty vector BOOST_CHECK_EQUAL(field.size(), 0); Dim_t len{2}; fc.initialise(CcoordOps::get_cube(len), {}); // check that returned size is correct BOOST_CHECK_EQUAL(field.size(), ipow(len, sdim)); // check that setting pad size won't change logical size field.set_pad_size(24); BOOST_CHECK_EQUAL(field.size(), ipow(len, sdim)); } BOOST_AUTO_TEST_CASE(dynamic_field_creation) { constexpr Dim_t sdim{threeD}; Dim_t nb_components{2}; using FC_t = GlobalFieldCollection; FC_t fc{}; make_field>("Dynamic Field", fc, nb_components); } + BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_zeros_like, Fix, field_fixtures, Fix) { + auto & t_clone{Fix::tensor_field.get_zeros_like("tensor clone")}; + static_assert(std::is_same< + std::remove_reference_t, + typename Fix::TField_t>::value, "wrong overloaded function"); + + auto & m_clone{Fix::matrix_field.get_zeros_like("matrix clone")}; + static_assert(std::is_same< + std::remove_reference_t, + typename Fix::MField_t>::value, "wrong overloaded function"); + using FieldColl_t = typename Fix::FieldColl_t; + using T = typename Fix::TField_t::Scalar; + TypedField & t_ref{t_clone}; + + auto & typed_clone{t_ref.get_zeros_like("dynamically sized clone")}; + static_assert(std::is_same< + std::remove_reference_t, + TypedField>::value, + "Field type incorrectly deduced"); + BOOST_CHECK_EQUAL(typed_clone.get_nb_components(), t_clone.get_nb_components()); + + auto & dyn_clone{Fix::dynamic_field1.get_zeros_like("dynamic clone")}; + static_assert(std::is_same::value, + "mismatch"); + BOOST_CHECK_EQUAL(typed_clone.get_nb_components(), dyn_clone.get_nb_components()); + } + + /* ---------------------------------------------------------------------- */ + BOOST_AUTO_TEST_CASE(fill_global_local) { + FieldFixture global; + FieldFixture local; + constexpr Dim_t len{2}; + constexpr auto sizes{CcoordOps::get_cube::SDim>(len)}; + + global.fc.initialise(sizes,{}); + + local.fc.add_pixel({1, 1}); + local.fc.add_pixel({0, 1}); + local.fc.initialise(); + + // fill the local matrix field and then transfer it to the global field + for (auto mat: local.matrix_field.get_map()) { + mat.setRandom(); + } + global.matrix_field.fill_from_local(local.matrix_field); + + for (const auto & ccoord: local.fc) { + const auto & a{local.matrix_field.get_map()[ccoord]}; + const auto & b{global.matrix_field.get_map()[ccoord]}; + const Real error{(a -b).norm()}; + BOOST_CHECK_EQUAL(error, 0.); + } + + // fill the global tensor field and then transfer it to the global field + for (auto mat: global.tensor_field.get_map()) { + mat.setRandom(); + } + local.tensor_field.fill_from_global(global.tensor_field); + for (const auto & ccoord: local.fc) { + const auto & a{local.matrix_field.get_map()[ccoord]}; + const auto & b{global.matrix_field.get_map()[ccoord]}; + const Real error{(a -b).norm()}; + BOOST_CHECK_EQUAL(error, 0.); + } + + BOOST_CHECK_THROW(local.tensor_field.fill_from_global(global.matrix_field), + std::runtime_error); + BOOST_CHECK_THROW(global.tensor_field.fill_from_local(local.matrix_field), + std::runtime_error); + } + BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_raw_field_map.cc b/tests/header_test_raw_field_map.cc similarity index 98% rename from tests/test_raw_field_map.cc rename to tests/header_test_raw_field_map.cc index 6eaf611..6fd75ab 100644 --- a/tests/test_raw_field_map.cc +++ b/tests/header_test_raw_field_map.cc @@ -1,93 +1,93 @@ /** - * file test_raw_field_map.cc + * file header_test_raw_field_map.cc * * @author Till Junge * * @date 17 Apr 2018 * * @brief tests for the raw field map type * * @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 "test_field_collections.hh" #include "common/field_map.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(raw_field_map_tests); BOOST_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) { using FC_t = typename F::Parent::FC_t; using MSqMap = MatrixFieldMap; MSqMap Mmap{F::fc["Tensorfield Real o2"]}; auto m_it = Mmap.begin(); auto m_it_end = Mmap.end(); RawFieldMap>> raw_map{Mmap.get_field().eigenvec()}; for (auto && mat: Mmap) { mat.setRandom(); } for (auto tup: akantu::zip(Mmap, raw_map)) { auto & mat_A = std::get<0>(tup); auto & mat_B = std::get<1>(tup); BOOST_CHECK_EQUAL((mat_A-mat_B).norm(), 0.); } Mmap.get_field().eigenvec().setZero(); for (auto && mat: raw_map) { mat.setIdentity(); } for (auto && mat: Mmap) { BOOST_CHECK_EQUAL((mat-mat.Identity()).norm(), 0.); } } BOOST_AUTO_TEST_CASE(Const_correctness_test) { Eigen::VectorXd vec1(12); vec1.setRandom(); RawFieldMap> map1{vec1}; static_assert(not map1.IsConst, "should not have been const"); RawFieldMap> cmap1{vec1}; static_assert(cmap1.IsConst, "should have been const"); const Eigen::VectorXd vec2{vec1}; RawFieldMap> cmap2{vec2}; } BOOST_AUTO_TEST_CASE(incompatible_size_check) { Eigen::VectorXd vec1(11); using RawFieldMap_t = RawFieldMap>; BOOST_CHECK_THROW(RawFieldMap_t {vec1}, std::runtime_error); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_statefields.cc b/tests/header_test_statefields.cc similarity index 99% rename from tests/test_statefields.cc rename to tests/header_test_statefields.cc index 119cf07..5963427 100644 --- a/tests/test_statefields.cc +++ b/tests/header_test_statefields.cc @@ -1,298 +1,298 @@ /** - * file test_statefields.cc + * file header_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{make_statefield>("prefix", fc)}, scalar_f{make_statefield>("scalar", fc)}, self{*this} {} FC_t fc; 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) { 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 * * @date 20 Nov 2017 * * @brief Test the fourth-order map on second-order tensor implementation * * 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 #include #include #include "common/common.hh" #include "tests.hh" #include "common/T4_map_proxy.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(T4map_tests); /** * Test fixture for construction of T4Map for the time being, symmetry is not * exploited */ template struct T4_fixture { T4_fixture():matrix{}, tensor(matrix.data()){} EIGEN_MAKE_ALIGNED_OPERATOR_NEW; using M4 = T4Mat; using T4 = T4MatMap; constexpr static Dim_t dim{Dim}; M4 matrix; T4 tensor; }; using fix_collection = boost::mpl::list, T4_fixture>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, fix_collection, F) { BOOST_CHECK_EQUAL(F::tensor.cols(), F::dim*F::dim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(write_access_test, F, fix_collection, F) { auto & t4 = F::tensor; constexpr Dim_t dim{F::dim}; Eigen::TensorFixedSize> t4c; Eigen::Map t4c_map(t4c.data()); for (Dim_t i = 0; i < F::dim; ++i) { for (Dim_t j = 0; j < F::dim; ++j) { for (Dim_t k = 0; k < F::dim; ++k) { for (Dim_t l = 0; l < F::dim; ++l) { get(t4,i,j,k,l) = 1000*(i+1) + 100*(j+1) + 10*(k+1) + l+1; t4c(i,j,k,l) = 1000*(i+1) + 100*(j+1) + 10*(k+1) + l+1; } } } } for (Dim_t i = 0; i < ipow(dim,4); ++i) { BOOST_CHECK_EQUAL(F::matrix.data()[i], t4c.data()[i]); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(assign_matrix_test, F, fix_collection, F) { decltype(F::matrix) matrix; matrix.setRandom(); F::tensor = matrix; for (Dim_t i = 0; i < ipow(F::dim,4); ++i) { BOOST_CHECK_EQUAL(F::matrix.data()[i], matrix.data()[i]); } } BOOST_AUTO_TEST_CASE(Return_ref_from_const_test) { constexpr Dim_t dim{2}; using T = int; using M4 = Eigen::Matrix; using M4c = const Eigen::Matrix; using T4 = T4MatMap; using T4c = T4MatMap; M4 mat; mat.setRandom(); M4c cmat{mat}; T4 tensor{mat.data()}; T4c ctensor{mat.data()}; T a = get(tensor,0,0,0,1); T b = get(ctensor,0,0,0,1); BOOST_CHECK_EQUAL(a, b); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_tensor_algebra.cc b/tests/header_test_tensor_algebra.cc similarity index 99% rename from tests/test_tensor_algebra.cc rename to tests/header_test_tensor_algebra.cc index 21c36e3..d584207 100644 --- a/tests/test_tensor_algebra.cc +++ b/tests/header_test_tensor_algebra.cc @@ -1,292 +1,292 @@ /** - * @file test_tensor_algebra.cc + * @file header_test_tensor_algebra.cc * * @author Till Junge * * @date 05 Nov 2017 * * @brief Tests for the tensor algebra functions * * 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 #include "common/tensor_algebra.hh" #include "tests.hh" #include "tests/test_goodies.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(tensor_algebra) auto TerrNorm = [](auto && t){ return Eigen::Tensor(t.abs().sum())(); }; /* ---------------------------------------------------------------------- */ BOOST_AUTO_TEST_CASE(tensor_outer_product_test) { constexpr Dim_t dim{2}; Eigen::TensorFixedSize> A, B; // use prime numbers so that every multiple is uniquely identifiable A.setValues({{1, 2}, {3 ,7}}); B.setValues({{11, 13}, {17 ,19}}); Eigen::TensorFixedSize> Res1, Res2, Res3; for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t k = 0; k < dim; ++k) { for (Dim_t l = 0; l < dim; ++l) { Res1(i, j, k, l) = A(i, j)*B(k, l); Res2(i, j, k, l) = A(i, k)*B(j, l); Res3(i, j, k, l) = A(i, l)*B(j, k); } } } } Real error = TerrNorm(Res1 - Tensors::outer(A,B)); BOOST_CHECK_LT(error, tol); error = TerrNorm(Res2 - Tensors::outer_under(A, B)); BOOST_CHECK_LT(error, tol); error = TerrNorm(Res3 - Tensors::outer_over(A, B)); if (error > tol) { std::cout << "reference:" << std::endl << Res3 << std::endl; std::cout << "result:" << std::endl << Tensors::outer_over(A, B) << std::endl; std::cout << "A:" << std::endl << A << std::endl; std::cout << "B" << std::endl << B << std::endl; decltype(Res3) tmp = Tensors::outer_over(A, B); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t k = 0; k < dim; ++k) { for (Dim_t l = 0; l < dim; ++l) { std::cout << "for (" << i << ", " << j << ", " << k << ", " << l << "), ref: " << std::setw(3)<< Res3(i,j,k,l) << ", res: " << std::setw(3)<< tmp(i,j,k,l) << std::endl; } } } } } BOOST_CHECK_LT(error, tol); error = TerrNorm(Res3 - Tensors::outer(A, B)); BOOST_CHECK_GT(error, tol); }; BOOST_FIXTURE_TEST_CASE_TEMPLATE(outer_products, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using T2 = Tensors::Tens2_t; using M2 = Matrices::Tens2_t; using Map2 = Eigen::Map; using T4 = Tensors::Tens4_t; using M4 = Matrices::Tens4_t; using Map4 = Eigen::Map; T2 A, B; T4 RT; A.setRandom(); B.setRandom(); Map2 Amap(A.data()); Map2 Bmap(B.data()); M2 C, D; M4 RM; C = Amap; D = Bmap; auto error =[](const T4& A, const M4& B) {return (B - Map4(A.data())).norm();}; // Check outer product RT = Tensors::outer(A, B); RM = Matrices::outer(C, D); BOOST_CHECK_LT(error(RT, RM), tol); // Check outer_under product RT = Tensors::outer_under(A, B); RM = Matrices::outer_under(C, D); BOOST_CHECK_LT(error(RT, RM), tol); // Check outer_over product RT = Tensors::outer_over(A, B); RM = Matrices::outer_over(C, D); BOOST_CHECK_LT(error(RT, RM), tol); } BOOST_AUTO_TEST_CASE(tensor_multiplication) { constexpr Dim_t dim{2}; using Strain_t = Eigen::TensorFixedSize>; Strain_t A, B; A.setValues({{1, 2}, {3 ,7}}); B.setValues({{11, 13}, {17 ,19}}); Strain_t FF1 = A*B; // element-wise multiplication std::array, 1> prod_dims { Eigen::IndexPair{1, 0}}; Strain_t FF2 = A.contract(B, prod_dims); // correct option 1 Strain_t FF3; using Mat_t = Eigen::Map>; // following only works for evaluated tensors (which already have data()) Mat_t(FF3.data()) = Mat_t(A.data()) * Mat_t(B.data()); Strain_t ref; ref.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t a = 0; a < dim; ++a) { ref(i,j) += A(i, a) * B(a, j); } } } using Strain_tw = Eigen::TensorFixedSize>; Strain_tw C; C.setConstant(100); //static_assert(!std::is_convertible::value, // "Tensors not size-protected"); if (std::is_convertible::value) { //std::cout << "this is not good, should I abandon Tensors?"; } // this test seems useless. I use to detect if Eigen changed the // default tensor product Real error = TerrNorm(FF1-ref); if (error < tol) { std::cout << "A =" << std::endl << A << std::endl; std::cout << "B =" << std::endl << B << std::endl; std::cout << "FF1 =" << std::endl << FF1 << std::endl; std::cout << "ref =" << std::endl << ref << std::endl; } BOOST_CHECK_GT(error, tol); error = TerrNorm(FF2-ref); if (error > tol) { std::cout << "A =" << std::endl << A << std::endl; std::cout << "B =" << std::endl << B << std::endl; std::cout << "FF2 =" << std::endl << FF2 << std::endl; std::cout << "ref =" << std::endl << ref << std::endl; } BOOST_CHECK_LT(error, tol); error = TerrNorm(FF3-ref); if (error > tol) { std::cout << "A =" << std::endl << A << std::endl; std::cout << "B =" << std::endl << B << std::endl; std::cout << "FF3 =" << std::endl << FF3 << std::endl; std::cout << "ref =" << std::endl << ref << std::endl; } BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensmult, Fix, testGoodies::dimlist, Fix) { using namespace Matrices; constexpr Dim_t dim{Fix::dim}; using T4 = Tens4_t; using T2 = Tens2_t; using V2 = Eigen::Matrix; T4 C; C.setRandom(); T2 E; E.setRandom(); Eigen::Map Ev(E.data()); T2 R = tensmult(C, E); auto error = (Eigen::Map(R.data())-C*Ev).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tracer, Fix, testGoodies::dimlist, Fix) { using namespace Matrices; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto tracer = Itrac(); T2 F; F.setRandom(); auto Ftrac = tensmult(tracer, F); auto error = (Ftrac-F.trace()*F.Identity()).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_identity, Fix, testGoodies::dimlist, Fix) { using namespace Matrices; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto ident = Iiden(); T2 F; F.setRandom(); auto Fiden = tensmult(ident, F); auto error = (Fiden-F).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_transposer, Fix, testGoodies::dimlist, Fix) { using namespace Matrices; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto trnst = Itrns(); T2 F; F.setRandom(); auto Ftrns = tensmult(trnst, F); auto error = (Ftrns-F.transpose()).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_symmetriser, Fix, testGoodies::dimlist, Fix) { using namespace Matrices; constexpr Dim_t dim{Fix::dim}; using T2 = Tens2_t; auto symmt = Isymm(); T2 F; F.setRandom(); auto Fsymm = tensmult(symmt, F); auto error = (Fsymm-.5*(F+F.transpose())).norm(); BOOST_CHECK_LT(error, tol); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_field_collections_1.cc b/tests/test_field_collections_1.cc deleted file mode 100644 index 2706a95..0000000 --- a/tests/test_field_collections_1.cc +++ /dev/null @@ -1,221 +0,0 @@ -/** - * @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{}; - using ftype = internal::TypedSizedFieldBase< - decltype(F::fc), Real, mdim*mdim*mdim*mdim>; - using stype = Eigen::Array; - 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_2.cc b/tests/test_field_collections_2.cc deleted file mode 100644 index c3c4201..0000000 --- a/tests/test_field_collections_2.cc +++ /dev/null @@ -1,275 +0,0 @@ -/** - * @file test_field_collections_2.cc - * - * @author Till Junge - * - * @date 23 Nov 2017 - * - * @brief Continuation of tests from test_field_collection.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 "test_field_collections.hh" -namespace muSpectre { - BOOST_AUTO_TEST_SUITE(field_collection_tests); - - BOOST_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) { - using FC_t = typename F::Parent::FC_t; - using Tensor4Map = TensorFieldMap; - Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; - F::fc["Tensorfield Real o4"].set_zero(); - for (auto && tens:T4map) { - BOOST_CHECK_EQUAL(Real(Eigen::Tensor(tens.abs().sum().eval())()), 0); - } - - using Tensor2Map = TensorFieldMap; - using MSqMap = MatrixFieldMap; - using ASqMap = ArrayFieldMap; - using A2Map = ArrayFieldMap; - using WrongMap = ArrayFieldMap; - Tensor2Map T2map{F::fc["Tensorfield Real o2"]}; - MSqMap Mmap{F::fc["Tensorfield Real o2"]}; - ASqMap Amap{F::fc["Tensorfield Real o2"]}; - A2Map DynMap{F::fc["Dynamically sized Field"]}; - auto & fc_ref{F::fc}; - BOOST_CHECK_THROW(WrongMap{fc_ref["Dynamically sized Field"]}, - FieldInterpretationError); - auto t2_it = T2map.begin(); - auto t2_it_end = T2map.end(); - auto m_it = Mmap.begin(); - auto a_it = Amap.begin(); - for (; t2_it != t2_it_end; ++t2_it, ++m_it, ++a_it) { - t2_it->setRandom(); - auto && m = *m_it; - bool comp = (m == a_it->matrix()); - BOOST_CHECK(comp); - } - - size_t counter{0}; - for (auto val: DynMap) { - ++counter; - val += val.Ones()*counter; - } - - counter = 0; - for (auto val: DynMap) { - ++counter; - val -= val.Ones()*counter; - auto error {val.matrix().norm()}; - BOOST_CHECK_LT(error, tol); - } - - - using ScalarMap = ScalarFieldMap; - ScalarMap s_map{F::fc["integer Scalar"]}; - for (Uint i = 0; i < s_map.size(); ++i) { - s_map[i] = i; - } - - counter = 0; - for (const auto& val: s_map) { - BOOST_CHECK_EQUAL(counter++, val); - } - - } - - BOOST_FIXTURE_TEST_CASE_TEMPLATE(ccoord_indexing_test, F, glob_iter_colls, F) { - using FC_t = typename F::Parent::FC_t; - using ScalarMap = ScalarFieldMap; - ScalarMap s_map{F::fc["integer Scalar"]}; - for (Uint i = 0; i < s_map.size(); ++i) { - s_map[i] = i; - } - - - for (size_t i = 0; i < CcoordOps::get_size(F::fc.get_sizes()); ++i) { - BOOST_CHECK_EQUAL(CcoordOps::get_index(F::fc.get_sizes(), - F::fc.get_locations(), - CcoordOps::get_ccoord(F::fc.get_sizes(), - F::fc.get_locations(), - i)), i); - } - - } - - BOOST_FIXTURE_TEST_CASE_TEMPLATE(iterator_methods_test, F, iter_collections, F) { - using FC_t = typename F::Parent::FC_t; - using Tensor4Map = TensorFieldMap; - Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; - using it_t = typename Tensor4Map::iterator; - std::ptrdiff_t diff{3}; // arbitrary, as long as it is smaller than the container size - - // check constructors - auto itstart = T4map.begin(); // standard way of obtaining iterator - auto itend = T4map.end(); // ditto - - it_t it1{T4map}; - it_t it2{T4map, false}; - it_t it3{T4map, size_t(diff)}; - BOOST_CHECK(itstart == itstart); - BOOST_CHECK(itstart != itend); - BOOST_CHECK_EQUAL(itstart, it1); - BOOST_CHECK_EQUAL(itend, it2); - - // check ostream operator - std::stringstream response; - response << it3; - BOOST_CHECK_EQUAL - (response.str(), - std::string - ("iterator on field 'Tensorfield Real o4', entry ") + - std::to_string(diff)); - - // check copy, move, and assigment constructor (and operator+) - it_t itcopy = it1; - it_t itmove = std::move(T4map.begin()); - it_t it4 = it1+diff; - it_t it5(it2); - it_t tmp(it2); - it_t it6(std::move(tmp)); - it_t it7 = it4 -diff; - BOOST_CHECK_EQUAL(itcopy, it1); - BOOST_CHECK_EQUAL(itmove, it1); - BOOST_CHECK_EQUAL(it4, it3); - BOOST_CHECK_EQUAL(it5, it2); - BOOST_CHECK_EQUAL(it6, it5); - BOOST_CHECK_EQUAL(it7, it1); - - // check increments/decrements - BOOST_CHECK_EQUAL(it1++, itcopy); // post-increment - BOOST_CHECK_EQUAL(it1, itcopy+1); - BOOST_CHECK_EQUAL(--it1, itcopy); // pre -decrement - BOOST_CHECK_EQUAL(++it1, itcopy+1); // pre -increment - BOOST_CHECK_EQUAL(it1--, itcopy+1); // post-decrement - BOOST_CHECK_EQUAL(it1, itcopy); - - - // dereference and member-of-pointer check - Eigen::Tensor Tens = *it1; - Eigen::Tensor Tens2 = *itstart; - Eigen::Tensor check = (Tens==Tens2).all(); - BOOST_CHECK_EQUAL(bool(check()), true); - - BOOST_CHECK_NO_THROW(itstart->setZero()); - - //check access subscripting - auto T3a = *it3; - auto T3b = itstart[diff]; - BOOST_CHECK(bool(Eigen::Tensor((T3a==T3b).all())())); - - // div. comparisons - BOOST_CHECK_LT(itstart, itend); - BOOST_CHECK(!(itend < itstart)); - BOOST_CHECK(!(itstart < itstart)); - - BOOST_CHECK_LE(itstart, itend); - BOOST_CHECK_LE(itstart, itstart); - BOOST_CHECK(!(itend <= itstart)); - - BOOST_CHECK_GT(itend, itstart); - BOOST_CHECK(!(itend>itend)); - BOOST_CHECK(!(itstart>itend)); - - BOOST_CHECK_GE(itend, itstart); - BOOST_CHECK_GE(itend, itend); - BOOST_CHECK(!(itstart >= itend)); - - // check assignment increment/decrement - BOOST_CHECK_EQUAL(it1+=diff, it3); - BOOST_CHECK_EQUAL(it1-=diff, itstart); - - // check cell coordinates - using Ccoord = Ccoord_t; - Ccoord a{itstart.get_ccoord()}; - Ccoord b{0}; - - // Weirdly, boost::has_left_shift::value is false for Ccoords, even though the operator is implemented :( - //BOOST_CHECK_EQUAL(a, b); - bool check2 = (a==b); - BOOST_CHECK(check2); - } - - BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_tensor_iter_test, F, iter_collections, F) { - using FC_t = typename F::Parent::FC_t; - using Tensor4Map = TensorFieldMap; - Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; - - using T_t = typename Tensor4Map::T_t; - Eigen::TensorMap Tens2(T4map[0].data(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim()); - - for (auto it = T4map.cbegin(); it != T4map.cend(); ++it) { - // maps to const tensors can't be initialised with a const pointer this sucks - auto&& tens = *it; - auto&& ptr = tens.data(); - - static_assert(std::is_pointer>::value, "should be getting a pointer"); - //static_assert(std::is_const>::value, "should be const"); - - // If Tensor were written well, above static_assert should pass, and the - // following check shouldn't. If it get triggered, it means that a newer - // version of Eigen now does have const-correct - // TensorMap. This means that const-correct field maps - // are then also possible for tensors - BOOST_CHECK(!std::is_const>::value); - - } - } - - BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_matrix_iter_test, F, iter_collections, F) { - using FC_t = typename F::Parent::FC_t; - using MatrixMap = MatrixFieldMap; - MatrixMap Mmap{F::fc["Matrixfield Complex sdim x mdim"]}; - - for (auto it = Mmap.cbegin(); it != Mmap.cend(); ++it) { - // maps to const tensors can't be initialised with a const pointer this sucks - auto&& mat = *it; - auto&& ptr = mat.data(); - - static_assert(std::is_pointer>::value, "should be getting a pointer"); - //static_assert(std::is_const>::value, "should be const"); - - // If Matrix were written well, above static_assert should pass, and the - // following check shouldn't. If it get triggered, it means that a newer - // version of Eigen now does have const-correct - // MatrixMap. This means that const-correct field maps - // are then also possible for matrices - BOOST_CHECK(!std::is_const>::value); - - } - } - - BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_scalar_iter_test, F, iter_collections, F) { - using FC_t = typename F::Parent::FC_t; - using ScalarMap = ScalarFieldMap; - ScalarMap Smap{F::fc["integer Scalar"]}; - - for (auto it = Smap.cbegin(); it != Smap.cend(); ++it) { - auto&& scal = *it; - static_assert(std::is_const>::value, - "referred type should be const"); - static_assert(std::is_lvalue_reference::value, - "Should have returned an lvalue ref"); - - } - } - - BOOST_AUTO_TEST_SUITE_END(); - -} // muSpectre diff --git a/tests/test_field_collections_3.cc b/tests/test_field_collections_3.cc deleted file mode 100644 index 02321fe..0000000 --- a/tests/test_field_collections_3.cc +++ /dev/null @@ -1,152 +0,0 @@ -/** - * @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