diff --git a/.gitignore b/.gitignore
index f75b574..4440c8b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,8 +1,5 @@
/build/CMake*
*~
/build/*
/build-o3/
-/build*
-
-/src/external/pybind11
-/src/external/Eigen3
+/build*
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ac8e3b8..5e30318 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,130 +1,134 @@
# =============================================================================
# file CMakeLists.txt
#
# @author Till Junge <till.junge@epfl.ch>
#
# @date 08 Jan 2018
#
# @brief Main configuration file
#
# @section LICENCE
#
# Copyright © 2018 Till Junge
#
# µSpectre is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3, or (at
# your option) any later version.
#
# µSpectre is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with GNU Emacs; see the file COPYING. If not, write to the
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
# =============================================================================
cmake_minimum_required(VERSION 3.7.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)
add_compile_options(-Wall -Wextra -Weffc++)
enable_testing()
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake)
include(muspectreTools)
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
# using Clang
add_compile_options(-Wno-missing-braces)
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
# using GCC
string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type )
if ("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()
find_package(Boost COMPONENTS unit_test_framework REQUIRED)
add_external_package(Eigen3 VERSION 3.3 CONFIG)
add_external_package(pybind11 VERSION 2.2 CONFIG)
#find_package(Eigen3 3.3 REQUIRED NO_MODULE)
include_directories(
- tests
- src
+ ${CMAKE_SOURCE_DIR}/tests
+ ${CMAKE_SOURCE_DIR}/src
+ ${CMAKE_SOURCE_DIR}
)
#find_package(MPI REQUIRED)
find_package(FFTW REQUIRED)
#find_package(FFTWMPI REQUIRED)
#build tests
file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/test_*.cc")
add_executable(main_test_suite tests/main_test_suite.cc ${TEST_SRCS})
target_link_libraries(main_test_suite ${Boost_LIBRARIES} muSpectre)
add_test(main_test_suite main_test_suite --report_level=detailed --build_info=TRUE)
#add_test(main_test_suite main_test_suite --build_info=TRUE)
# compile the library
add_compile_options( -Werror)
add_subdirectory(
- "${CMAKE_SOURCE_DIR}/src/"
+ ${CMAKE_SOURCE_DIR}/src/
+ )
+add_subdirectory(
+ ${CMAKE_SOURCE_DIR}/language_bindings/
)
#compile executables
file( GLOB binaries "${CMAKE_SOURCE_DIR}/bin/*.cc")
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})
# compile 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})
# copy python test
configure_file(
"${CMAKE_SOURCE_DIR}/tests/python_binding_tests.py"
"${CMAKE_BINARY_DIR}/python_binding_tests.py"
COPYONLY)
add_test(python_binding_test python_binding_tests.py)
diff --git a/bin/small_case.py b/bin/small_case.py
index db5966f..9ff2693 100755
--- a/bin/small_case.py
+++ b/bin/small_case.py
@@ -1,71 +1,71 @@
#!/usr/bin/env python3
"""
file small_case.py
@author Till Junge <till.junge@epfl.ch>
@date 12 Jan 2018
@brief small case for debugging
@section LICENCE
Copyright © 2018 Till Junge
µSpectre is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation, either version 3, or (at
your option) any later version.
µSpectre is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with GNU Emacs; see the file COPYING. If not, write to the
Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.
"""
import sys
import os
import numpy as np
sys.path.append(os.path.join(os.getcwd(), "src/language_bindings/python"))
import pyMuSpectre as µ
resolution = [3, 3]
lengths = [3., 3.]
-formulation = µ.Formulation.small_strain
+formulation = µ.Formulation.finite_strain
rve = µ.SystemFactory(resolution,
lengths,
formulation)
hard = µ.material.MaterialHooke2d.make(
rve, "hard", 210e9, .33)
soft = µ.material.MaterialHooke2d.make(
rve, "soft", 70e9, .33)
for i, pixel in enumerate(rve):
if i < 3:
hard.add_pixel(pixel)
else:
soft.add_pixel(pixel)
print("{}, {}".format(i, tuple(pixel)))
rve.initialise()
tol = 1e-6
Del0 = np.array([[0, .1],
[0, 0]])
maxiter = 31
verbose = 2
r = µ.solvers.de_geus2d(rve, Del0, tol, tol, maxiter, verbose)
print(r) #print(r.eigen().T)
diff --git a/cmake/muspectreTools.cmake b/cmake/muspectreTools.cmake
index be7d8ea..db855d3 100644
--- a/cmake/muspectreTools.cmake
+++ b/cmake/muspectreTools.cmake
@@ -1,214 +1,214 @@
#==============================================================================
# file muspectreTools.cmake
#
# @author Nicolas Richart <nicolas.,richart@epfl.ch>
#
# @date 11 Jan 2018
#
# @brief some tool to help to do stuff with cmake in µSpectre
#
# @section LICENCE
#
# Copyright © 2018 Till Junge
#
# µSpectre is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3, or (at
# your option) any later version.
#
# µSpectre is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with GNU Emacs; see the file COPYING. If not, write to the
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
#
#[=[.rst:
µSpectreTools
-------------
This module provide some helper fuunctions for µSpectre
::
download_external_project(project_name
URL <url>
BACKEND <GIT|HG>
[TAG <tag>]
[THIRD_PARTY_SRC_DIR <dir>]
[NO_UPDATE]
)
Downloads the external project based on the uri
``THIRD_PARTY_SRC_DIR <dir>``
Specifies where to download the source
``NO_UPDATE``
Does not try to update existing download
::
mark_as_advanced_prefix(prefix)
Marks as advanced all variables whoes names start with prefix
::
add_external_dep(package
[IGNORE_SYSTEM]
[VERSION <version>]
[...])
Tries to find the external package and if not found but a local
${package}.cmake files exists, it includes it. The extra arguments are passed
to find_package()
``IGNORE_SYSTEM`` does not do the find_package
``VERSION`` Specifies the required minimum version
]=]
#
function(download_external_project project_name)
include(CMakeParseArguments)
set(_dep_flags
NO_UPDATE)
set(_dep_one_variables
URL
TAG
BACKEND
THIRD_PARTY_SRC_DIR
)
set(_dep_multi_variables)
cmake_parse_arguments(_dep_args
"${_dep_flags}"
"${_dep_one_variables}"
"${_dep_multi_variables}"
${ARGN}
)
if(NOT _dep_args_URL)
message(FATAL_ERROR "You have to provide a URL for the project ${project_name}")
endif()
if(NOT _dep_args_BACKEND)
message(FATAL_ERROR "You have to provide a backend to download ${project_name}")
endif()
if(_dep_args_TAG)
set(_ep_tag "${_dep_args_BACKEND}_TAG ${_dep_args_TAG}")
endif()
set(_src_dir ${PROJECT_SOURCE_DIR}/third-party/${project_name})
if (_dep_args_THIRD_PARTY_SRC_DIR)
set(_src_dir ${_dep_args_THIRD_PARTY_SRC_DIR}/${project_name})
endif()
if(EXISTS ${_src_dir} AND _dep_args_NO_UPDATE)
return()
endif()
set(_working_dir ${PROJECT_BINARY_DIR}/third-party/${project_name}-download)
file(WRITE ${_working_dir}/CMakeLists.txt
"
cmake_minimum_required(VERSION 3.1)
project(${project_name}-download NONE)
include(ExternalProject)
ExternalProject_Add(${project_name}
SOURCE_DIR ${_src_dir}
BINARY_DIR ${_working_dir}
${_dep_args_BACKEND}_REPOSITORY ${_dep_args_URL}
${_ep_tag}
CONFIGURE_COMMAND \"\"
BUILD_COMMAND \"\"
INSTALL_COMMAND \"\"
TEST_COMMAND \"\"
)
")
message(STATUS "Downloading ${project_name} ${_dep_args_GIT_TAG}")
execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" .
RESULT_VARIABLE _result
WORKING_DIRECTORY ${_working_dir}
OUTPUT_FILE ${_working_dir}/configure-out.log
ERROR_FILE ${_working_dir}/configure-error.log)
if(_result)
message(FATAL_ERROR "Something went wrong (${_result}) during the download"
" process of ${project_name} check the file"
" ${_working_dir}/configure-error.log for more details")
endif()
execute_process(COMMAND "${CMAKE_COMMAND}" --build .
RESULT_VARIABLE _result
WORKING_DIRECTORY ${_working_dir}
OUTPUT_FILE ${_working_dir}/build-out.log
ERROR_FILE ${_working_dir}/build-error.log)
if(_result)
message(FATAL_ERROR "Something went wrong (${_result}) during the download"
" process of ${project_name} check the file"
" ${_working_dir}/build-error.log for more details")
endif()
endfunction()
# ------------------------------------------------------------------------------
function(mark_as_advanced_prefix prefix)
get_property(_list DIRECTORY PROPERTY VARIABLES)
foreach(_var ${_list})
if("${_var}" MATCHES "^${prefix}.*")
mark_as_advanced(${_var})
endif()
endforeach()
endfunction()
# ------------------------------------------------------------------------------
function(add_external_package package)
include(CMakeParseArguments)
set(_cmake_includes ${PROJECT_SOURCE_DIR}/cmake)
- set(_${package}_external_dir ${PROJECT_SOURCE_DIR}/src/external)
+ set(_${package}_external_dir ${PROJECT_BINARY_DIR}/external)
set(_aep_flags
IGNORE_SYSTEM
)
set(_aep_one_variables
VERSION
)
set(_aep_multi_variables)
cmake_parse_arguments(_aep_args
"${_aep_flags}"
"${_aep_one_variables}"
"${_aep_multi_variables}"
${ARGN}
)
if(_aep_args_VERSION)
set(_${package}_version ${_aep_args_VERSION})
endif()
if(NOT EXISTS ${_cmake_includes}/${package}.cmake)
set(_required REQUIRED)
endif()
if(NOT _aep_args_IGNORE_SYSTEM)
find_package(${package} ${_${package}_version} ${_required} ${_aep_UNPARSED_ARGUMENTS} QUIET)
if(${package}_FOUND AND NOT ${package}_FOUND_EXTERNAL)
return()
endif()
endif()
if(EXISTS ${_cmake_includes}/${package}.cmake)
include(${_cmake_includes}/${package}.cmake)
endif()
endfunction()
diff --git a/src/external/cxxopts.hpp b/external/cxxopts.hpp
similarity index 100%
rename from src/external/cxxopts.hpp
rename to external/cxxopts.hpp
diff --git a/src/language_bindings/CMakeLists.txt b/language_bindings/CMakeLists.txt
similarity index 100%
rename from src/language_bindings/CMakeLists.txt
rename to language_bindings/CMakeLists.txt
diff --git a/src/language_bindings/python/CMakeLists.txt b/language_bindings/python/CMakeLists.txt
similarity index 100%
rename from src/language_bindings/python/CMakeLists.txt
rename to language_bindings/python/CMakeLists.txt
diff --git a/src/language_bindings/python/bind_py_common.cc b/language_bindings/python/bind_py_common.cc
similarity index 100%
rename from src/language_bindings/python/bind_py_common.cc
rename to language_bindings/python/bind_py_common.cc
diff --git a/src/language_bindings/python/bind_py_declarations.hh b/language_bindings/python/bind_py_declarations.hh
similarity index 100%
rename from src/language_bindings/python/bind_py_declarations.hh
rename to language_bindings/python/bind_py_declarations.hh
diff --git a/src/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc
similarity index 100%
rename from src/language_bindings/python/bind_py_material.cc
rename to language_bindings/python/bind_py_material.cc
diff --git a/src/language_bindings/python/bind_py_module.cc b/language_bindings/python/bind_py_module.cc
similarity index 100%
rename from src/language_bindings/python/bind_py_module.cc
rename to language_bindings/python/bind_py_module.cc
diff --git a/src/language_bindings/python/bind_py_solvers.cc b/language_bindings/python/bind_py_solvers.cc
similarity index 100%
rename from src/language_bindings/python/bind_py_solvers.cc
rename to language_bindings/python/bind_py_solvers.cc
diff --git a/src/language_bindings/python/bind_py_system.cc b/language_bindings/python/bind_py_system.cc
similarity index 100%
rename from src/language_bindings/python/bind_py_system.cc
rename to language_bindings/python/bind_py_system.cc
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0c8568a..b25dd39 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,66 +1,64 @@
# =============================================================================
# file CMakeLists.txt
#
# @author Till Junge <till.junge@epfl.ch>
#
# @date 08 Jan 2018
#
# @brief Configuration for libmuSpectre
#
# @section LICENCE
#
# Copyright © 2018 Till Junge
#
# µSpectre is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3, or (at
# your option) any later version.
#
# µSpectre is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with GNU Emacs; see the file COPYING. If not, write to the
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
# =============================================================================
include_directories(
common
materials
system
- external
)
set (muSpectre_SRC
)
add_subdirectory(common)
add_subdirectory(materials)
add_subdirectory(fft)
add_subdirectory(system)
add_subdirectory(solver)
find_package(FFTW REQUIRED)
# The following checks whether std::optional exists and replaces it by
# boost::optional if necessary
include(CheckCXXSourceCompiles)
CHECK_CXX_SOURCE_COMPILES(
"#include <experimental/optional>
int main() {
std::experimental::optional<double> A{};
}"
HAS_STD_OPTIONAL)
if( NOT HAS_STD_OPTIONAL)
add_definitions(-DNO_EXPERIMENTAL)
endif()
add_library(muSpectre ${muSpectre_SRC})
target_link_libraries(muSpectre Eigen3::Eigen ${FFTW_LIBRARIES})
-add_subdirectory(language_bindings)
diff --git a/src/common/field.hh b/src/common/field.hh
index 1ef3565..bf7daba 100644
--- a/src/common/field.hh
+++ b/src/common/field.hh
@@ -1,728 +1,778 @@
/**
* file field.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 07 Sep 2017
*
* @brief header-only implementation of a field for field collections
*
* @section LICENCE
*
* 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 <Eigen/Dense>
#include <string>
#include <sstream>
#include <utility>
#include <typeinfo>
#include <vector>
#include <algorithm>
#include <cmath>
#include <memory>
#include <type_traits>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
class FieldCollectionError: public std::runtime_error {
public:
explicit FieldCollectionError(const std::string& what)
:std::runtime_error(what){}
explicit FieldCollectionError(const char * what)
:std::runtime_error(what){}
};
class FieldError: public FieldCollectionError {
using Parent = FieldCollectionError;
public:
explicit FieldError(const std::string& what)
:Parent(what){}
explicit FieldError(const char * what)
:Parent(what){}
};
class FieldInterpretationError: public FieldError
{
public:
explicit FieldInterpretationError(const std::string & what)
:FieldError(what){}
explicit FieldInterpretationError(const char * what)
:FieldError(what){}
};
namespace internal {
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
class FieldBase
{
protected:
//! constructor
//! unique name (whithin Collection)
//! number of components
//! collection to which this field belongs (eg, material, system)
FieldBase(std::string unique_name,
size_t nb_components,
FieldCollection & collection);
public:
using collection_t = FieldCollection;
//! Copy constructor
FieldBase(const FieldBase &other) = delete;
//! Move constructor
- FieldBase(FieldBase &&other) noexcept = delete;
+ FieldBase(FieldBase &&other) = delete;
//! Destructor
- virtual ~FieldBase() noexcept = default;
+ virtual ~FieldBase() = default;
//! Copy assignment operator
FieldBase& operator=(const FieldBase &other) = delete;
//! Move assignment operator
- FieldBase& operator=(FieldBase &&other) noexcept = delete;
+ FieldBase& operator=(FieldBase &&other) = delete;
/* ---------------------------------------------------------------------- */
//!Identifying accessors
//! return field name
inline const std::string & get_name() const;
//! return field type
//inline const Field_t & get_type() const;
//! return my collection (for iterating)
inline const FieldCollection & get_collection() const;
//! return my collection (for iterating)
inline const size_t & get_nb_components() const;
//! return type_id of stored type
virtual const std::type_info & get_stored_typeid() const = 0;
virtual size_t size() const = 0;
//! initialise field to zero (do more complicated initialisations through
//! fully typed maps)
virtual void set_zero() = 0;
//! give access to collections
friend FieldCollection;
friend typename FieldCollection::Parent;
protected:
/* ---------------------------------------------------------------------- */
//! allocate memory etc
virtual void resize(size_t size) = 0;
const std::string name;
const size_t nb_components;
const FieldCollection & collection;
private:
};
+ /**
+ * dummy intermediate class to provide a run-time polymorphic
+ * typed field. Mainly for binding python
+ */
+ template <class FieldCollection, typename T>
+ class TypedFieldBase: public FieldBase<FieldCollection>
+ {
+ public:
+ using Parent = FieldBase<FieldCollection>;
+ using collection_t = typename Parent::collection_t;
+ using Scalar = T;
+ using Base = Parent;
+ using EigenRep = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic>;
+ using EigenMap = Eigen::Map<EigenRep>;
+ //! Default constructor
+ TypedFieldBase() = delete;
+
+ TypedFieldBase(std::string unique_name,
+ size_t nb_components,
+ FieldCollection& collection);
+
+ //! Copy constructor
+ TypedFieldBase(const TypedFieldBase &other) = delete;
+
+ //! Move constructor
+ TypedFieldBase(TypedFieldBase &&other) = delete;
+
+ //! Destructor
+ virtual ~TypedFieldBase() = default;
+
+ //! Copy assignment operator
+ TypedFieldBase& operator=(const TypedFieldBase &other) = delete;
+
+ //! Move assignment operator
+ TypedFieldBase& operator=(TypedFieldBase &&other) = delete;
+
+ //! return type_id of stored type
+ virtual const std::type_info & get_stored_typeid() const override final;
+
+
+ protected:
+ private:
+ };
+
/* ---------------------------------------------------------------------- */
//! declaraton for friending
template <class FieldCollection, typename T, Dim_t NbComponents, bool isConst>
class FieldMap;
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore=false>
- class TypedFieldBase: public FieldBase<FieldCollection>
+ class TypedSizedFieldBase: public TypedFieldBase<FieldCollection, T>
{
friend class FieldMap<FieldCollection, T, NbComponents, true>;
friend class FieldMap<FieldCollection, T, NbComponents, false>;
public:
constexpr static auto nb_components{NbComponents};
- using Parent = FieldBase<FieldCollection>;
+ using Parent = TypedFieldBase<FieldCollection, T>;
using collection_t = typename Parent::collection_t;
using Scalar = T;
- using Base = Parent;
+ using Base = typename Parent::Base;
//using storage_type = Eigen::Array<T, Eigen::Dynamic, NbComponents>;
using StoredType = Eigen::Array<T, NbComponents, 1>;
using StorageType = std::conditional_t
<ArrayStore,
std::vector<StoredType,
Eigen::aligned_allocator<StoredType>>,
std::vector<T,Eigen::aligned_allocator<T>>>;
using EigenRep = Eigen::Array<T, NbComponents, Eigen::Dynamic>;
using EigenMap = std::conditional_t<
ArrayStore,
Eigen::Map<EigenRep, Eigen::Aligned,
Eigen::OuterStride<sizeof(StoredType)/sizeof(T)>>,
Eigen::Map<EigenRep>>;
using ConstEigenMap = std::conditional_t<
ArrayStore,
Eigen::Map<const EigenRep, Eigen::Aligned,
Eigen::OuterStride<sizeof(StoredType)/sizeof(T)>>,
Eigen::Map<const EigenRep>>;
- TypedFieldBase(std::string unique_name,
+ TypedSizedFieldBase(std::string unique_name,
FieldCollection& collection);
- virtual ~TypedFieldBase() = default;
- //! return type_id of stored type
- virtual const std::type_info & get_stored_typeid() const override final;
+ virtual ~TypedSizedFieldBase() = default;
//! initialise field to zero (do more complicated initialisations through
//! fully typed maps)
inline void set_zero() override final;
template <bool isArrayStore = ArrayStore>
inline void push_back(const
std::enable_if_t<isArrayStore, StoredType> &
value);
template <bool componentStore = !ArrayStore>
inline std::enable_if_t<componentStore> push_back(const StoredType & value);
//! Number of stored arrays (i.e. total number of stored
//! scalars/NbComponents)
size_t size() const override final;
- static TypedFieldBase & check_ref(Parent & other);
- static const TypedFieldBase & check_ref(const Base & parent);
+ static TypedSizedFieldBase & check_ref(Base & other);
+ static const TypedSizedFieldBase & check_ref(const Base & parent);
inline T* data() {return this->get_ptr_to_entry(0);}
inline const T* data() const {return this->get_ptr_to_entry(0);}
inline EigenMap eigen();
inline ConstEigenMap eigen() const;
template<typename otherT>
- inline Real inner_product(const TypedFieldBase<FieldCollection, otherT, NbComponents,
+ inline Real inner_product(const TypedSizedFieldBase<FieldCollection, otherT, NbComponents,
ArrayStore> & other) const;
protected:
template <bool isArray=ArrayStore>
inline std::enable_if_t<isArray, T*> get_ptr_to_entry(const size_t&& index);
template <bool isArray=ArrayStore>
inline std::enable_if_t<isArray, const T*>
get_ptr_to_entry(const size_t&& index) const;
template <bool noArray = !ArrayStore>
inline T* get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index);
template <bool noArray = !ArrayStore>
inline const T*
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const;
inline virtual void resize(size_t size) override final;
StorageType values{};
};
} // internal
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
- class TensorField: public internal::TypedFieldBase<FieldCollection,
+ class TensorField: public internal::TypedSizedFieldBase<FieldCollection,
T,
ipow(dim,order)>
{
public:
- using Parent = internal::TypedFieldBase<FieldCollection,
+ using Parent = internal::TypedSizedFieldBase<FieldCollection,
T,
ipow(dim,order)>;
using Base = typename Parent::Base;
using Field_p = typename FieldCollection::Field_p;
using component_type = T;
using Scalar = typename Parent::Scalar;
//! Copy constructor
TensorField(const TensorField &other) = delete;
//! Move constructor
- TensorField(TensorField &&other) noexcept = delete;
+ TensorField(TensorField &&other) = delete;
//! Destructor
- virtual ~TensorField() noexcept = default;
+ virtual ~TensorField() = default;
//! Copy assignment operator
TensorField& operator=(const TensorField &other) = delete;
//! Move assignment operator
- TensorField& operator=(TensorField &&other) noexcept = delete;
+ TensorField& operator=(TensorField &&other) = delete;
//! accessors
inline Dim_t get_order() const;
inline Dim_t get_dim() const;
//! factory function
template<class FieldType, class CollectionType, typename... Args>
friend FieldType& make_field(std::string unique_name,
CollectionType & collection,
Args&&... args);
static TensorField & check_ref(Base & other) {
return static_cast<TensorField &>(Parent::check_ref(other));}
static const TensorField & check_ref(const Base & other) {
return static_cast<const TensorField &>(Parent::check_ref(other));}
/**
* 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) and the correct size (nbComponents).
*/
decltype(auto) get_map();
decltype(auto) get_map() const;
protected:
//! constructor protected!
TensorField(std::string unique_name, FieldCollection & collection);
private:
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol=NbRow>
- class MatrixField: public internal::TypedFieldBase<FieldCollection,
+ class MatrixField: public internal::TypedSizedFieldBase<FieldCollection,
T,
NbRow*NbCol>
{
public:
- using Parent = internal::TypedFieldBase<FieldCollection,
+ using Parent = internal::TypedSizedFieldBase<FieldCollection,
T,
NbRow*NbCol>;
using Base = typename Parent::Base;
using Field_p = std::unique_ptr<internal::FieldBase<FieldCollection>>;
using component_type = T;
//! Copy constructor
MatrixField(const MatrixField &other) = delete;
//! Move constructor
- MatrixField(MatrixField &&other) noexcept = delete;
+ MatrixField(MatrixField &&other) = delete;
//! Destructor
- virtual ~MatrixField() noexcept = default;
+ virtual ~MatrixField() = default;
//! Copy assignment operator
MatrixField& operator=(const MatrixField &other) = delete;
//! Move assignment operator
- MatrixField& operator=(MatrixField &&other) noexcept = delete;
+ MatrixField& operator=(MatrixField &&other) = delete;
//! accessors
inline Dim_t get_nb_row() const;
inline Dim_t get_nb_col() const;
//! factory function
template<class FieldType, class CollectionType, typename... Args>
friend FieldType& make_field(std::string unique_name,
CollectionType & collection,
Args&&... args);
static MatrixField & check_ref(Base & other) {
return static_cast<MatrixField &>(Parent::check_ref(other));}
static const MatrixField & check_ref(const Base & other) {
return static_cast<const MatrixField &>(Parent::check_ref(other));}
decltype(auto) get_map();
decltype(auto) get_map() const;
protected:
//! constructor protected!
MatrixField(std::string unique_name, FieldCollection & collection);
private:
};
/* ---------------------------------------------------------------------- */
//! convenience alias (
template <class FieldCollection, typename T>
using ScalarField = MatrixField<FieldCollection, T, 1, 1>;
/* ---------------------------------------------------------------------- */
// Implementations
namespace internal {
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
FieldBase<FieldCollection>::FieldBase(std::string unique_name,
size_t nb_components_,
FieldCollection & collection_)
:name(unique_name), nb_components(nb_components_),
collection(collection_) {}
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
inline const std::string & FieldBase<FieldCollection>::get_name() const {
return this->name;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
inline const FieldCollection & FieldBase<FieldCollection>::
get_collection() const {
return this->collection;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
inline const size_t & FieldBase<FieldCollection>::
get_nb_components() const {
return this->nb_components;
}
+ /* ---------------------------------------------------------------------- */
+ template <class FieldCollection, typename T>
+ TypedFieldBase<FieldCollection, T>::
+ TypedFieldBase(std::string unique_name, size_t nb_components,
+ FieldCollection & collection)
+ :Parent(unique_name, nb_components, collection)
+ {}
+
+ /* ---------------------------------------------------------------------- */
+ //! return type_id of stored type
+ template <class FieldCollection, typename T>
+ const std::type_info & TypedFieldBase<FieldCollection, T>::
+ get_stored_typeid() const {
+ return typeid(T);
+ }
+
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
- TypedFieldBase(std::string unique_name, FieldCollection & collection)
- :FieldBase<FieldCollection>(unique_name, NbComponents, collection){
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ TypedSizedFieldBase(std::string unique_name, FieldCollection & collection)
+ :Parent(unique_name, NbComponents, collection){
static_assert
((std::is_arithmetic<T>::value ||
std::is_same<T, Complex>::value),
- "Use TypedFieldBase for integer, real or complex scalars for T");
+ "Use TypedSizedFieldBase for integer, real or complex scalars for T");
static_assert(NbComponents > 0,
"Only fields with more than 0 components");
}
- /* ---------------------------------------------------------------------- */
- //! return type_id of stored type
- template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
- const std::type_info & TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
- get_stored_typeid() const {
- return typeid(T);
- }
-
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
void
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
set_zero() {
std::fill(this->values.begin(), this->values.end(), T{});
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
- size_t TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ size_t TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
size() const {
if (ArrayStore) {
return this->values.size();
} else {
return this->values.size()/NbComponents;
}
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
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<TypedFieldBase&>(other);
+ return static_cast<TypedSizedFieldBase&>(other);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
- const TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ const TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
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<const TypedFieldBase&>(other);
+ return static_cast<const TypedSizedFieldBase&>(other);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
- typename TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::EigenMap
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ typename TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::EigenMap
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
eigen() {
return EigenMap(this->data(), NbComponents, this->size());
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
- typename TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::ConstEigenMap
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ typename TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::ConstEigenMap
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
eigen() const {
return ConstEigenMap(this->data(), NbComponents, this->size());
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
template<typename otherT>
Real
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
- inner_product(const TypedFieldBase<FieldCollection, otherT, NbComponents,
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ inner_product(const TypedSizedFieldBase<FieldCollection, otherT, NbComponents,
ArrayStore> & other) const {
return (this->eigen() * other.eigen()).sum();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArray>
std::enable_if_t<isArray, T*>
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(const size_t&& index) {
static_assert (isArray == ArrayStore, "SFINAE");
return &this->values[std::move(index)](0, 0);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool noArray>
- T* TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ T* TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) {
static_assert (noArray != ArrayStore, "SFINAE");
return &this->values[NbComponents*std::move(index)];
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArray>
std::enable_if_t<isArray, const T*>
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(const size_t&& index) const {
static_assert (isArray == ArrayStore, "SFINAE");
return &this->values[std::move(index)](0, 0);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool noArray>
- const T* TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ const T* TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const {
static_assert (noArray != ArrayStore, "SFINAE");
return &this->values[NbComponents*std::move(index)];
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
- void TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ void TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
resize(size_t size) {
if (ArrayStore) {
this->values.resize(size);
} else {
this->values.resize(size*NbComponents);
}
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArrayStore>
void
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
push_back(const std::enable_if_t<isArrayStore,StoredType> & value) {
static_assert(isArrayStore == ArrayStore, "SFINAE");
this->values.push_back(value);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool componentStore>
std::enable_if_t<componentStore>
- TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
+ TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
push_back(const StoredType & value) {
static_assert(componentStore != ArrayStore, "SFINAE");
for (Dim_t i = 0; i < NbComponents; ++i) {
this->values.push_back(value(i));
}
}
} // internal
/* ---------------------------------------------------------------------- */
//! Factory function, guarantees that only fields get created that are
//! properly registered and linked to a collection.
template<class FieldType, class FieldCollection, typename... Args>
FieldType &
make_field(std::string unique_name,
FieldCollection & collection,
Args&&... args) {
auto ptr = std::unique_ptr<FieldType>{
new FieldType(unique_name, collection, args...)};
auto& retref = *ptr;
collection.register_field(std::move(ptr));
return retref;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
TensorField<FieldCollection, T, order, dim>::
TensorField(std::string unique_name, FieldCollection & collection)
:Parent(unique_name, collection) {}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
Dim_t TensorField<FieldCollection, T, order, dim>::
get_order() const {
return order;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
Dim_t TensorField<FieldCollection, T, order, dim>::
get_dim() const {
return dim;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
MatrixField<FieldCollection, T, NbRow, NbCol>::
MatrixField(std::string unique_name, FieldCollection & collection)
:Parent(unique_name, collection) {}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
get_nb_col() const {
return NbCol;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
get_nb_row() const {
return NbRow;
}
} // muSpectre
#include "common/field_map.hh"
namespace muSpectre {
namespace internal {
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, size_t order, Dim_t dim,
bool ConstMap>
struct tensor_map_type {
};
template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
struct tensor_map_type<FieldCollection, T, firstOrder, dim, ConstMap> {
using type = MatrixFieldMap<FieldCollection, T, dim, 1, ConstMap>;
};
template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
struct tensor_map_type<FieldCollection, T, secondOrder, dim, ConstMap> {
using type = MatrixFieldMap<FieldCollection, T, dim, dim, ConstMap>;
};
template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
struct tensor_map_type<FieldCollection, T, fourthOrder, dim, ConstMap> {
using type = T4MatrixFieldMap<FieldCollection, T, dim, ConstMap>;
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol,
bool ConstMap>
struct matrix_map_type {
using type =
MatrixFieldMap<FieldCollection, T, NbRow, NbCol, ConstMap>;
};
template <class FieldCollection, typename T, bool ConstMap>
struct matrix_map_type<FieldCollection, T, oneD, oneD, ConstMap> {
using type = ScalarFieldMap<FieldCollection, T, ConstMap>;
};
} // internal
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
decltype(auto) TensorField<FieldCollection, T, order, dim>::
get_map() {
constexpr bool map_constness{false};
using RawMap_t =
typename internal::tensor_map_type<FieldCollection, T, order, dim,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
decltype(auto) TensorField<FieldCollection, T, order, dim>::
get_map() const {
constexpr bool map_constness{true};
using RawMap_t =
typename internal::tensor_map_type<FieldCollection, T, order, dim,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
decltype(auto) MatrixField<FieldCollection, T, NbRow, NbCol>::
get_map() {
constexpr bool map_constness{false};
using RawMap_t =
typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
decltype(auto) MatrixField<FieldCollection, T, NbRow, NbCol>::
get_map() const {
constexpr bool map_constness{true};
using RawMap_t =
typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
map_constness>::type;
return RawMap_t(*this);
}
} // muSpectre
#endif /* FIELD_H */
diff --git a/src/common/field_collection_base.hh b/src/common/field_collection_base.hh
index 23599f8..c386359 100644
--- a/src/common/field_collection_base.hh
+++ b/src/common/field_collection_base.hh
@@ -1,162 +1,162 @@
/**
* file field_collection_base.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief Base class for field collections
*
* @section LICENCE
*
* 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 <map>
#include <vector>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
//! DimS spatial dimension (dimension of problem)
//! DimM material_dimension (dimension of constitutive law)
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
class FieldCollectionBase
{
public:
using Field = internal::FieldBase<FieldCollectionDerived>;
using Field_p = std::unique_ptr<Field>;
using Ccoord = Ccoord_t<DimS>;
//! Default constructor
FieldCollectionBase();
//! Copy constructor
FieldCollectionBase(const FieldCollectionBase &other) = delete;
//! Move constructor
- FieldCollectionBase(FieldCollectionBase &&other) = delete; //noexcept = default;
+ FieldCollectionBase(FieldCollectionBase &&other) = delete;
//! Destructor
- virtual ~FieldCollectionBase() noexcept = default;
+ virtual ~FieldCollectionBase() = default;
//! Copy assignment operator
FieldCollectionBase& operator=(const FieldCollectionBase &other) = delete;
//! Move assignment operator
- FieldCollectionBase& operator=(FieldCollectionBase &&other) = delete; //noexcept = default;
+ 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);
//! for return values of iterators
constexpr inline static Dim_t spatial_dim();
//! for return values of iterators
inline Dim_t get_spatial_dim() const;
//! retrieve field by unique_name
inline Field& operator[](std::string unique_name);
//! retrieve field by unique_name with bounds checking
inline Field& at(std::string unique_name);
//! 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_;}
protected:
std::map<const std::string, Field_p> fields{};
bool is_initialised{false};
const Uint id;
static Uint counter;
size_t size_{0};
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
Uint FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::counter{0};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::FieldCollectionBase()
:id(counter++){}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
void FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::register_field(Field_p &&field) {
auto&& search_it = this->fields.find(field->get_name());
auto&& does_exist = search_it != this->fields.end();
if (does_exist) {
std::stringstream err_str;
err_str << "a field named " << field->get_name()
<< "is already registered in this field collection. "
<< "Currently registered fields: ";
for (const auto& name_field_pair: this->fields) {
err_str << ", " << name_field_pair.first;
}
throw FieldCollectionError(err_str.str());
}
if (this->is_initialised) {
field->resize(this->size());
}
this->fields[field->get_name()] = std::move(field);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
constexpr Dim_t FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
spatial_dim() {
return DimS;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
Dim_t FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
get_spatial_dim() const {
return DimS;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
typename FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::Field&
FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
operator[](std::string unique_name) {
return *(this->fields[unique_name]);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
typename FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::Field&
FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
at(std::string unique_name) {
return *(this->fields.at(unique_name));
}
} // muSpectre
#endif /* FIELD_COLLECTION_BASE_H */
diff --git a/src/common/field_collection_global.hh b/src/common/field_collection_global.hh
index 57fe34a..09824fc 100644
--- a/src/common/field_collection_global.hh
+++ b/src/common/field_collection_global.hh
@@ -1,192 +1,192 @@
/**
* file field_collection_global.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief FieldCollection base-class for global fields
*
* @section LICENCE
*
* 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 {
/* ---------------------------------------------------------------------- */
//! DimS spatial dimension (dimension of problem)
//! DimM material_dimension (dimension of constitutive law)
template<Dim_t DimS, Dim_t DimM>
class GlobalFieldCollection:
public FieldCollectionBase<DimS, DimM, GlobalFieldCollection<DimS, DimM>>
{
public:
using Parent = FieldCollectionBase
<DimS, DimM, GlobalFieldCollection<DimS, DimM>>;
using Ccoord = typename Parent::Ccoord;
using Field_p = typename Parent::Field_p;
using iterator = typename CcoordOps::Pixels<DimS>::iterator;
//! Default constructor
GlobalFieldCollection();
//! Copy constructor
GlobalFieldCollection(const GlobalFieldCollection &other) = delete;
//! Move constructor
GlobalFieldCollection
- (GlobalFieldCollection &&other) noexcept = default;
+ (GlobalFieldCollection &&other) = default;
//! Destructor
- virtual ~GlobalFieldCollection() noexcept = default;
+ 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);
//! return the pixel sizes
inline const Ccoord & get_sizes() const;
//! returns the linear index corresponding to cell coordinates
template <class CcoordRef>
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();
inline iterator end();
static constexpr inline Dim_t spatial_dim() {return DimS;}
static constexpr inline Dim_t material_dim() {return DimM;}
protected:
//! number of discretisation cells in each of the DimS spatial directions
Ccoord sizes{};
Ccoord strides{};
CcoordOps::Pixels<DimS> pixels{};
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
GlobalFieldCollection<DimS, DimM>::GlobalFieldCollection()
:Parent()
{}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void GlobalFieldCollection<DimS, DimM>::
initialise(Ccoord sizes) {
if (this->is_initialised) {
throw std::runtime_error("double initialisation");
}
this->pixels = CcoordOps::Pixels<DimS>(sizes);
this->size_ = CcoordOps::get_size(sizes);
this->sizes = sizes;
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 the pixel sizes
template <Dim_t DimS, Dim_t DimM>
const typename GlobalFieldCollection<DimS, DimM>::Ccoord &
GlobalFieldCollection<DimS, DimM>::get_sizes() const {
return this->sizes;
}
//----------------------------------------------------------------------------//
//! returns the cell coordinates corresponding to a linear index
template <Dim_t DimS, Dim_t DimM>
typename GlobalFieldCollection<DimS, DimM>::Ccoord
GlobalFieldCollection<DimS, DimM>::get_ccoord(size_t index) const {
return CcoordOps::get_ccoord(this->get_sizes(), std::move(index));
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename GlobalFieldCollection<DimS, DimM>::iterator
GlobalFieldCollection<DimS, DimM>::begin() {
return this->pixels.begin();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename GlobalFieldCollection<DimS, DimM>::iterator
GlobalFieldCollection<DimS, DimM>::end() {
return this->pixels.end();
}
//----------------------------------------------------------------------------//
//! returns the linear index corresponding to cell coordinates
template <Dim_t DimS, Dim_t DimM>
template <class CcoordRef>
size_t
GlobalFieldCollection<DimS, DimM>::get_index(CcoordRef && ccoord) const {
static_assert(std::is_same<
Ccoord,
std::remove_const_t<
std::remove_reference_t<CcoordRef>>>::value,
"can only be called with values or references of Ccoord");
return CcoordOps::get_index(this->get_sizes(),
std::forward<CcoordRef>(ccoord));
}
} // muSpectre
#endif /* FIELD_COLLECTION_GLOBAL_H */
diff --git a/src/common/field_collection_local.hh b/src/common/field_collection_local.hh
index a4b15bf..3017afc 100644
--- a/src/common/field_collection_local.hh
+++ b/src/common/field_collection_local.hh
@@ -1,174 +1,174 @@
/**
* file field_collection_local.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief FieldCollection base-class for local fields
*
* @section LICENCE
*
* 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 {
/* ---------------------------------------------------------------------- */
//! DimS spatial dimension (dimension of problem)
//! DimM material_dimension (dimension of constitutive law)
template<Dim_t DimS, Dim_t DimM>
class LocalFieldCollection:
public FieldCollectionBase<DimS, DimM, LocalFieldCollection<DimS, DimM>>
{
public:
using Parent = FieldCollectionBase<DimS, DimM,
LocalFieldCollection<DimS, DimM>>;
using Ccoord = typename Parent::Ccoord;
using Field_p = typename Parent::Field_p;
using ccoords_container = std::vector<Ccoord>;
using iterator = typename ccoords_container::iterator;
//! Default constructor
LocalFieldCollection();
//! Copy constructor
LocalFieldCollection(const LocalFieldCollection &other) = delete;
//! Move constructor
- LocalFieldCollection(LocalFieldCollection &&other) noexcept = delete;
+ LocalFieldCollection(LocalFieldCollection &&other) = delete;
//! Destructor
- virtual ~LocalFieldCollection() noexcept = default;
+ virtual ~LocalFieldCollection() = default;
//! Copy assignment operator
LocalFieldCollection& operator=(const LocalFieldCollection &other) = delete;
//! Move assignment operator
- LocalFieldCollection& operator=(LocalFieldCollection &&other) noexcept = delete;
+ 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 <class CcoordRef>
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() {return this->ccoords.begin();}
inline iterator end() {return this->ccoords.end();}
protected:
//! container of pixel coords for non-global collections
ccoords_container ccoords{};
//! container of indices for non-global collections (slow!)
std::map<Ccoord, size_t> indices{};
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
LocalFieldCollection<DimS, DimM>::LocalFieldCollection()
:Parent(){}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void LocalFieldCollection<DimS, DimM>::
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);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void LocalFieldCollection<DimS, DimM>::
initialise() {
if (this->is_initialised) {
throw std::runtime_error("double initialisation");
}
this->size_ = this->ccoords.size();
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 <Dim_t DimS, Dim_t DimM>
template <class CcoordRef>
size_t
LocalFieldCollection<DimS, DimM>::get_index(CcoordRef && ccoord) const {
static_assert(std::is_same<
Ccoord,
std::remove_const_t<
std::remove_reference_t<CcoordRef>>>::value,
"can only be called with values or references of Ccoord");
return this->indices.at(std::forward<CcoordRef>(ccoord));
}
//----------------------------------------------------------------------------//
//! returns the cell coordinates corresponding to a linear index
template <Dim_t DimS, Dim_t DimM>
typename LocalFieldCollection<DimS, DimM>::Ccoord
LocalFieldCollection<DimS, DimM>::get_ccoord(size_t index) const {
return this->ccoords[std::move(index)];
}
} // muSpectre
#endif /* FIELD_COLLECTION_LOCAL_H */
diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh
index f2686b6..58dcdf4 100644
--- a/src/common/field_map_base.hh
+++ b/src/common/field_map_base.hh
@@ -1,614 +1,614 @@
/**
* file field_map.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 12 Sep 2017
*
* @brief Defined a strongly defines proxy that iterates efficiently over a field
*
* @section LICENCE
*
* 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_H
#define FIELD_MAP_H
#include "common/field.hh"
#include "field_collection_base.hh"
#include "common/common.hh"
#include <unsupported/Eigen/CXX11/Tensor>
#include <array>
#include <string>
#include <memory>
#include <type_traits>
namespace muSpectre {
namespace internal {
//----------------------------------------------------------------------------//
//! little helper to automate creation of const maps without duplication
template<class T, bool isConst>
struct const_corrector {
using type = typename T::reference;
};
template<class T>
struct const_corrector<T, true> {
using type = typename T::const_reference;
};
template<class T, bool isConst>
using const_corrector_t = typename const_corrector<T, isConst>::type;
//----------------------------------------------------------------------------//
template <class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
class FieldMap
{
public:
constexpr static auto nb_components{NbComponents};
- using TypedField_nc = TypedFieldBase
+ using TypedField_nc = TypedSizedFieldBase
<FieldCollection, T, NbComponents>;
using TypedField = std::conditional_t<ConstField,
const TypedField_nc,
TypedField_nc>;
- using Field = typename TypedField::Parent;
+ using Field = typename TypedField::Base;
using size_type = std::size_t;
using pointer = std::conditional_t<ConstField,
const T*,
T*>;
//! Default constructor
FieldMap() = delete;
template <bool isntConst=!ConstField>
FieldMap(std::enable_if_t<isntConst, Field &> field);
template <bool isConst=ConstField>
FieldMap(std::enable_if_t<isConst, const Field &> field);
template<class FC, typename T2, Dim_t NbC>
- FieldMap(TypedFieldBase<FC, T2, NbC> & field);
+ FieldMap(TypedSizedFieldBase<FC, T2, NbC> & field);
//! Copy constructor
FieldMap(const FieldMap &other) = default;
//! Move constructor
- FieldMap(FieldMap &&other) noexcept = default;
+ FieldMap(FieldMap &&other) = default;
//! Destructor
- virtual ~FieldMap() noexcept = default;
+ virtual ~FieldMap() = default;
//! Copy assignment operator
FieldMap& operator=(const FieldMap &other) = delete;
//! Move assignment operator
- FieldMap& operator=(FieldMap &&other) noexcept = delete;
+ 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<class TypedField>
struct is_compatible;
template <class FullyTypedFieldMap, bool ConstIter=false>
class iterator
{
static_assert(!((ConstIter==false) && (ConstField==true)),
"You can't have a non-const iterator over a const "
"field");
public:
using value_type =
const_corrector_t<FullyTypedFieldMap, ConstIter>;
using const_value_type =
const_corrector_t<FullyTypedFieldMap, true>;
using pointer = typename FullyTypedFieldMap::pointer;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
using Ccoord = typename FieldCollection::Ccoord;
using reference = typename FullyTypedFieldMap::reference;
using TypedRef = std::conditional_t<ConstIter,
const FullyTypedFieldMap &,
FullyTypedFieldMap>;
//! 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) noexcept = default;
+ iterator(iterator &&other) = default;
//! Destructor
- virtual ~iterator() noexcept = default;
+ virtual ~iterator() = default;
//! Copy assignment operator
iterator& operator=(const iterator &other) = default;
//! Move assignment operator
- iterator& operator=(iterator &&other) noexcept = default;
+ 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;
inline bool operator<=(const iterator & other) const;
inline bool operator>(const iterator & other) const;
inline bool operator>=(const iterator & other) const;
//! additions, subtractions and corresponding assignments
inline iterator operator+(difference_type diff) const;
inline iterator operator-(difference_type diff) const;
inline iterator& operator+=(difference_type diff);
inline iterator& operator-=(difference_type diff);
//! get pixel coordinates
inline Ccoord get_ccoord() const;
//! ostream operator (mainly for debug
friend std::ostream & operator<<(std::ostream & os,
const iterator& it) {
if (ConstIter) {
os << "const ";
}
os << "iterator on field '"
<< it.fieldmap.get_name()
<< "', entry " << it.index;
return os;
}
protected:
const FieldCollection & collection;
TypedRef fieldmap;
size_t index;
private:
};
protected:
inline pointer get_ptr_to_entry(size_t index);
inline const T* get_ptr_to_entry(size_t index) const;
const FieldCollection & collection;
TypedField & field;
private:
};
} // internal
namespace internal {
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <bool isntConst>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
FieldMap(std::enable_if_t<isntConst, Field &> field)
:collection(field.get_collection()), field(static_cast<TypedField&>(field)) {
static_assert(NbComponents>0,
"Only fields with more than 0 components allowed");
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <bool isConst>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
FieldMap(std::enable_if_t<isConst, const Field &> field)
:collection(field.get_collection()), field(static_cast<const TypedField&>(field)) {
static_assert(NbComponents>0,
"Only fields with more than 0 components allowed");
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <class FC, typename T2, Dim_t NbC>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
- FieldMap(TypedFieldBase<FC, T2, NbC> & field)
+ FieldMap(TypedSizedFieldBase<FC, T2, NbC> & field)
:collection(field.get_collection()), field(static_cast<TypedField&>(field)) {
static_assert(std::is_same<FC, FieldCollection>::value,
"The field does not have the expected FieldCollection type");
static_assert(std::is_same<T2, T>::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<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
void
FieldMap<FieldCollection, T, NbComponents, ConstField>::
check_compatibility() {
if (typeid(T).hash_code() !=
this->field.get_stored_typeid().hash_code()) {
std::string err{"Cannot create a Map of type '" +
this->info_string() +
"' for field '" + this->field.get_name() + "' of type '" +
this->field.get_stored_typeid().name() + "'"};
throw FieldInterpretationError
(err);
}
//check size compatibility
if (NbComponents != this->field.get_nb_components()) {
throw FieldInterpretationError
("Cannot create a Map of type '" +
this->info_string() +
"' for field '" + this->field.get_name() + "' with " +
std::to_string(this->field.get_nb_components()) + " components");
}
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
size_t
FieldMap<FieldCollection, T, NbComponents, ConstField>::
size() const {
return this->collection.size();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <class myField>
struct FieldMap<FieldCollection, T, NbComponents, ConstField>::is_compatible {
constexpr static bool explain() {
static_assert
(std::is_same<typename myField::collection_t, FieldCollection>::value,
"The field does not have the expected FieldCollection type");
static_assert
(std::is_same<typename myField::Scalar, T>::value,
"The // field does not have the expected Scalar type");
static_assert((TypedField::nb_components == NbComponents),
"The field does not have the expected number of components");
//The static asserts wouldn't pass in the incompatible case, so this is it
return true;
}
constexpr static bool value{std::is_base_of<TypedField, myField>::value};
};
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
const std::string &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_name() const {
return this->field.get_name();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
const FieldCollection &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_collection() const {
return this->collection;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
// Iterator implementations
//! constructor
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator
<FullyTypedFieldMap, ConstIter>::
iterator(TypedRef fieldmap, bool begin)
:collection(fieldmap.get_collection()), fieldmap(fieldmap),
index(begin ? 0 : fieldmap.field.size()) {}
/* ---------------------------------------------------------------------- */
//! constructor for random access
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator<FullyTypedFieldMap, ConstIter>::
iterator(TypedRef fieldmap, size_t index)
:collection(fieldmap.collection), fieldmap(fieldmap),
index(index) {}
/* ---------------------------------------------------------------------- */
//! pre-increment
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator<FullyTypedFieldMap, ConstIter>::
operator++() {
this->index++;
return *this;
}
/* ---------------------------------------------------------------------- */
//! post-increment
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator<FullyTypedFieldMap, ConstIter>::
operator++(int) {
iterator current = *this;
this->index++;
return current;
}
/* ---------------------------------------------------------------------- */
//! dereference
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
template iterator<FullyTypedFieldMap, ConstIter>::value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator*() {
return this->fieldmap.operator[](this->index);
}
/* ---------------------------------------------------------------------- */
//! dereference
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
template iterator<FullyTypedFieldMap, ConstIter>::const_value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator*() const {
return this->fieldmap.operator[](this->index);
}
/* ---------------------------------------------------------------------- */
//! member of pointer
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FullyTypedFieldMap::pointer
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator->() {
return this->fieldmap.ptr_to_val_t(this->index);
}
/* ---------------------------------------------------------------------- */
//! pre-decrement
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator--() {
this->index--;
return *this;
}
/* ---------------------------------------------------------------------- */
//! post-decrement
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator--(int) {
iterator current = *this;
this->index--;
return current;
}
/* ---------------------------------------------------------------------- */
//! Access subscripting
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
template iterator<FullyTypedFieldMap, ConstIter>::value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator[](difference_type diff) {
return this->fieldmap[this->index+diff];
}
/* ---------------------------------------------------------------------- */
//! Access subscripting
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template iterator<FullyTypedFieldMap, ConstIter>::const_value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator[](const difference_type diff) const {
return this->fieldmap[this->index+diff];
}
/* ---------------------------------------------------------------------- */
//! equality
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator==(const iterator & other) const {
return (this->index == other.index);
}
/* ---------------------------------------------------------------------- */
//! inquality
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator!=(const iterator & other) const {
return !(*this == other);
}
/* ---------------------------------------------------------------------- */
//! div. comparisons
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator<(const iterator & other) const {
return (this->index < other.index);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator<=(const iterator & other) const {
return (this->index <= other.index);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator>(const iterator & other) const {
return (this->index > other.index);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator>=(const iterator & other) const {
return (this->index >= other.index);
}
/* ---------------------------------------------------------------------- */
//! additions, subtractions and corresponding assignments
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator+(difference_type diff) const {
return iterator(this->fieldmap, this->index + diff);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator-(difference_type diff) const {
return iterator(this->fieldmap, this->index - diff);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator+=(difference_type diff) {
this->index += diff;
return *this;
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator-=(difference_type diff) {
this->index -= diff;
return *this;
}
/* ---------------------------------------------------------------------- */
//! get pixel coordinates
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldCollection::Ccoord
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
get_ccoord() const {
return this->collection.get_ccoord(this->index);
}
////----------------------------------------------------------------------------//
//template<class FieldCollection, typename T, Dim_t NbComponents, class FullyTypedFieldMap>
//std::ostream & operator <<
//(std::ostream &os,
// const typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
// template iterator<FullyTypedFieldMap, ConstIter> & it) {
// os << "iterator on field '"
// << it.field.get_name()
// << "', entry " << it.index;
// return os;
//}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::pointer
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_ptr_to_entry(size_t index) {
return this->field.get_ptr_to_entry(std::move(index));
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
const T*
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_ptr_to_entry(size_t index) const {
return this->field.get_ptr_to_entry(std::move(index));
}
} // internal
} // muSpectre
#endif /* FIELD_MAP_H */
diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh
index 9d21a04..6d8e9cf 100644
--- a/src/common/field_map_matrixlike.hh
+++ b/src/common/field_map_matrixlike.hh
@@ -1,358 +1,358 @@
/**
* file field_map_matrixlike.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 26 Sep 2017
*
* @brief Eigen-Matrix and -Array maps over strongly typed fields
*
* @section LICENCE
*
* 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 <Eigen/Dense>
#include <memory>
namespace muSpectre {
namespace internal {
/* ---------------------------------------------------------------------- */
enum class Map_t{Matrix, Array, T4Matrix};
template<Map_t map_type>
struct NameWrapper {
};
template<>
struct NameWrapper<Map_t::Array> {
static std::string field_info_root() {return "Array";}
};
template<>
struct NameWrapper<Map_t::Matrix> {
static std::string field_info_root() {return "Matrix";}
};
template<>
struct NameWrapper<Map_t::T4Matrix> {
static std::string field_info_root() {return "T4Matrix";}
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
class MatrixLikeFieldMap: public FieldMap<FieldCollection,
typename EigenArray::Scalar,
EigenArray::SizeAtCompileTime,
ConstField>
{
public:
using Parent = FieldMap<FieldCollection,
typename EigenArray::Scalar,
EigenArray::SizeAtCompileTime, ConstField>;
using T_t = EigenPlain;
using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
using value_type = EigenArray;
using const_reference = EigenConstArray;
using reference = std::conditional_t<ConstField,
const_reference,
value_type>; // since it's a resource handle
using size_type = typename Parent::size_type;
using pointer = std::unique_ptr<EigenArray>;
using TypedField = typename Parent::TypedField;
- using Field = typename TypedField::Parent;
+ using Field = typename TypedField::Base;
using const_iterator= typename Parent::template iterator<MatrixLikeFieldMap, true>;
using iterator = std::conditional_t<
ConstField,
const_iterator,
typename Parent::template iterator<MatrixLikeFieldMap>>;
using reverse_iterator = std::reverse_iterator<iterator>;
using const_reverse_iterator = std::reverse_iterator<const_iterator>;
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).
*/
template <bool isntConst=!ConstField>
MatrixLikeFieldMap(std::enable_if_t<isntConst, Field &> field);
template <bool isConst=ConstField>
MatrixLikeFieldMap(std::enable_if_t<isConst, const Field &> 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<class FC, typename T2, Dim_t NbC>
- MatrixLikeFieldMap(TypedFieldBase<FC, T2, NbC> & field);
+ MatrixLikeFieldMap(TypedSizedFieldBase<FC, T2, NbC> & field);
//! Copy constructor
MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = default;
//! Move constructor
- MatrixLikeFieldMap(MatrixLikeFieldMap &&other) noexcept = default;
+ MatrixLikeFieldMap(MatrixLikeFieldMap &&other) = default;
//! Destructor
- virtual ~MatrixLikeFieldMap() noexcept = default;
+ virtual ~MatrixLikeFieldMap() = default;
//! Copy assignment operator
MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete;
//! Move assignment operator
- MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) noexcept = delete;
+ MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) = delete;
//! Assign a matrixlike value to every entry
template <class Derived>
inline MatrixLikeFieldMap & operator=(const Eigen::EigenBase<Derived> & 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 head of field for ranges
inline iterator begin(){return iterator(*this);}
inline const_iterator cbegin() const {return const_iterator(*this);}
inline const_iterator begin() const {return this->cbegin();}
//! return an iterator to tail of field for ranges
inline iterator end(){return iterator(*this, false);};
inline const_iterator cend() const {return const_iterator(*this, false);}
inline const_iterator end() const {return this->cend();}
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;
private:
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
const std::string MatrixLikeFieldMap<FieldCollection, EigenArray,
EigenConstArray, EigenPlain, map_type,
ConstField>::
field_info_root{NameWrapper<map_type>::field_info_root()};
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template <bool isntConst>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
MatrixLikeFieldMap(std::enable_if_t<isntConst, Field &> field)
:Parent(field) {
static_assert((isntConst != ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template <bool isConst>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
MatrixLikeFieldMap(std::enable_if_t<isConst, const Field &> field)
:Parent(field) {
static_assert((isConst == ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template<class FC, typename T2, Dim_t NbC>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
- MatrixLikeFieldMap(TypedFieldBase<FC, T2, NbC> & field)
+ MatrixLikeFieldMap(TypedSizedFieldBase<FC, T2, NbC> & field)
:Parent(field) {
}
/* ---------------------------------------------------------------------- */
//! human-readable field map type
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
std::string
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
info_string() const {
std::stringstream info;
info << field_info_root << "("
<< typeid(typename EigenArray::value_type).name() << ", "
<< EigenArray::RowsAtCompileTime << "x"
<< EigenArray::ColsAtCompileTime << ")";
return info.str();
}
/* ---------------------------------------------------------------------- */
//! member access
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
operator[](size_type index) {
return reference(this->get_ptr_to_entry(index));
}
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
operator[](const Ccoord & ccoord) {
size_t index{};
index = this->collection.get_index(ccoord);
return reference(this->get_ptr_to_entry(std::move(index)));
}
/* ---------------------------------------------------------------------- */
//! member access
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::const_reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
operator[](size_type index) const {
return const_reference(this->get_ptr_to_entry(index));
}
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::const_reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
operator[](const Ccoord & ccoord) const{
size_t index{};
index = this->collection.get_index(ccoord);
return const_reference(this->get_ptr_to_entry(std::move(index)));
}
//----------------------------------------------------------------------------//
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template <class Derived>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField> &
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
operator=(const Eigen::EigenBase<Derived> & val) {
for (auto && entry: *this) {
entry = val;
}
return *this;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::T_t
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
mean() const {
T_t mean{T_t::Zero()};
for (auto && val: *this) {
mean += val;
}
mean *= 1./Real(this->size());
return mean;
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::pointer
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
ptr_to_val_t(size_type index) {
return std::make_unique<value_type>
(this->get_ptr_to_entry(std::move(index)));
}
} // internal
/* ---------------------------------------------------------------------- */
//! short-hand for an Eigen matrix map as iterate
template <class FieldCollection, typename T, Dim_t NbRows, Dim_t NbCols,
bool ConstField=false>
using MatrixFieldMap = internal::MatrixLikeFieldMap
<FieldCollection,
std::conditional_t<ConstField,
Eigen::Map<const Eigen::Matrix<T, NbRows, NbCols>>,
Eigen::Map<Eigen::Matrix<T, NbRows, NbCols>>>,
Eigen::Map<const Eigen::Matrix<T, NbRows, NbCols>>,
Eigen::Matrix<T, NbRows, NbCols>,
internal::Map_t::Matrix, ConstField>;
/* ---------------------------------------------------------------------- */
//! short-hand for an Eigen matrix map as iterate
template <class FieldCollection, typename T, Dim_t Dim,
bool MapConst=false>
using T4MatrixFieldMap = internal::MatrixLikeFieldMap
<FieldCollection,
T4MatMap<T, Dim, MapConst>,
T4MatMap<T, Dim, true>,
T4Mat<T, Dim>,
internal::Map_t::T4Matrix, MapConst>;
/* ---------------------------------------------------------------------- */
//! short-hand for an Eigen array map as iterate
template <class FieldCollection, typename T, Dim_t NbRows, Dim_t NbCols,
bool ConstField=false>
using ArrayFieldMap = internal::MatrixLikeFieldMap
<FieldCollection,
std::conditional_t<ConstField,
Eigen::Map<const Eigen::Array<T, NbRows, NbCols>>,
Eigen::Map<Eigen::Array<T, NbRows, NbCols>>>,
Eigen::Map<const Eigen::Array<T, NbRows, NbCols>>,
Eigen::Array<T, NbRows, NbCols>,
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 52707bb..914f88d 100644
--- a/src/common/field_map_scalar.hh
+++ b/src/common/field_map_scalar.hh
@@ -1,217 +1,217 @@
/**
* file field_map_scalar.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 26 Sep 2017
*
* @brief maps over scalar fields
*
* @section LICENCE
*
* 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 {
template <class FieldCollection, typename T, bool ConstField=false>
class ScalarFieldMap
: public internal::FieldMap<FieldCollection, T, 1, ConstField>
{
public:
using parent = internal::FieldMap<FieldCollection, T, 1, ConstField>;
using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
using value_type = T;
using const_reference = const value_type &;
using reference = std::conditional_t<ConstField,
const_reference,
value_type &>;
using size_type = typename parent::size_type;
using pointer = T*;
using Field = typename parent::Field;
using const_iterator= typename parent::template iterator<ScalarFieldMap, true>;
using iterator = std::conditional_t
<ConstField,
const_iterator,
typename parent::template iterator<ScalarFieldMap>>;
using reverse_iterator = std::reverse_iterator<iterator>;
using const_reverse_iterator = std::reverse_iterator<const_iterator>;
friend iterator;
//! Default constructor
ScalarFieldMap() = delete;
template <bool isntConst=!ConstField>
ScalarFieldMap(std::enable_if_t<isntConst, Field &> field);
template <bool isConst=ConstField>
ScalarFieldMap(std::enable_if_t<isConst, const Field &> field);
//! Copy constructor
ScalarFieldMap(const ScalarFieldMap &other) = default;
//! Move constructor
- ScalarFieldMap(ScalarFieldMap &&other) noexcept = default;
+ ScalarFieldMap(ScalarFieldMap &&other) = default;
//! Destructor
- virtual ~ScalarFieldMap() noexcept = default;
+ virtual ~ScalarFieldMap() = default;
//! Copy assignment operator
ScalarFieldMap& operator=(const ScalarFieldMap &other) = delete;
//! Move assignment operator
- ScalarFieldMap& operator=(ScalarFieldMap &&other) noexcept = delete;
+ 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 head of field for ranges
inline iterator begin(){return iterator(*this);}
inline const_iterator cbegin() const {return const_iterator(*this);}
inline const_iterator begin() const {return this->cbegin();}
//! return an iterator to tail of field for ranges
inline iterator end(){return iterator(*this, false);}
inline const_iterator cend() const {return const_iterator(*this, false);}
inline const_iterator end() const {return this->cend();}
inline 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;
private:
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
template <bool isntConst>
ScalarFieldMap<FieldCollection, T, ConstField>::
ScalarFieldMap(std::enable_if_t<isntConst, Field &> field)
:parent(field) {
static_assert((isntConst != ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
template <bool isConst>
ScalarFieldMap<FieldCollection, T, ConstField>::
ScalarFieldMap(std::enable_if_t<isConst, const Field &> field)
:parent(field) {
static_assert((isConst == ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
//! human-readable field map type
template<class FieldCollection, typename T, bool ConstField>
std::string
ScalarFieldMap<FieldCollection, T, ConstField>::info_string() const {
std::stringstream info;
info << "Scalar(" << typeid(T).name() << ")";
return info.str();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
const std::string ScalarFieldMap<FieldCollection, T, ConstField>::field_info_root{
"Scalar"};
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::reference
ScalarFieldMap<FieldCollection, T, ConstField>::operator[](size_type index) {
return this->get_ptr_to_entry(std::move(index))[0];
}
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::reference
ScalarFieldMap<FieldCollection, T, ConstField>::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 <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::const_reference
ScalarFieldMap<FieldCollection, T, ConstField>::
operator[](size_type index) const {
return this->get_ptr_to_entry(std::move(index))[0];
}
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::const_reference
ScalarFieldMap<FieldCollection, T, ConstField>::
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 <class FieldCollection, typename T, bool ConstField>
ScalarFieldMap<FieldCollection, T, ConstField> &
ScalarFieldMap<FieldCollection, T, ConstField>::operator=(T val) {
for (auto & scalar:*this) {
scalar = val;
}
return *this;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
T
ScalarFieldMap<FieldCollection, T, ConstField>::
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 95d1b2b..496e9e6 100644
--- a/src/common/field_map_tensor.hh
+++ b/src/common/field_map_tensor.hh
@@ -1,264 +1,264 @@
/**
* file field_map_tensor.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 26 Sep 2017
*
* @brief Defines an Eigen-Tensor map over strongly typed fields
*
* @section LICENCE
*
* 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 <sstream>
#include <memory>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField=false>
class TensorFieldMap: public
internal::FieldMap<FieldCollection,
T,
SizesByOrder<order, dim>::Sizes::total_size,
ConstField>
{
public:
using Parent = internal::FieldMap<FieldCollection, T,
TensorFieldMap::nb_components, ConstField>;
using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
using Sizes = typename SizesByOrder<order, dim>::Sizes;
using T_t = Eigen::TensorFixedSize<T, Sizes>;
using value_type = Eigen::TensorMap<T_t>;
using const_reference = Eigen::TensorMap<const T_t>;
using reference = std::conditional_t<ConstField,
const_reference,
value_type>; // since it's a resource handle
using size_type = typename Parent::size_type;
using pointer = std::unique_ptr<value_type>;
using TypedField = typename Parent::TypedField;
- using Field = typename TypedField::Parent;
+ using Field = typename TypedField::Base;
using const_iterator = typename Parent::template iterator<TensorFieldMap, true>;
using iterator = std::conditional_t<
ConstField,
const_iterator,
typename Parent::template iterator<TensorFieldMap>>;
using reverse_iterator = std::reverse_iterator<iterator>;
using const_reverse_iterator = std::reverse_iterator<const_iterator>;
friend iterator;
//! Default constructor
TensorFieldMap() = delete;
template <bool isntConst=!ConstField>
TensorFieldMap(std::enable_if_t<isntConst, Field &> field);
template <bool isConst=ConstField>
TensorFieldMap(std::enable_if_t<isConst, const Field &> field);
//! Copy constructor
TensorFieldMap(const TensorFieldMap &other) = default;
//! Move constructor
- TensorFieldMap(TensorFieldMap &&other) noexcept = default;
+ TensorFieldMap(TensorFieldMap &&other) = default;
//! Destructor
- virtual ~TensorFieldMap() noexcept = default;
+ 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) noexcept = delete;
+ 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);
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 head of field for ranges
inline iterator begin(){return iterator(*this);}
inline const_iterator cbegin() const {return const_iterator(*this);}
inline const_iterator begin() const {return this->cbegin();}
//! return an iterator to tail of field for ranges
inline iterator end(){return iterator(*this, false);}
inline const_iterator cend() const {return const_iterator(*this, false);}
inline const_iterator end() const {return this->cend();}
inline T_t mean() const;
protected:
//! for sad, legacy iterator use
inline pointer ptr_to_val_t(size_type index);
private:
};
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
template <bool isntConst>
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
TensorFieldMap(std::enable_if_t<isntConst, Field &> field)
:Parent(field) {
static_assert((isntConst != ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
template <bool isConst>
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
TensorFieldMap(std::enable_if_t<isConst, const Field &> field)
:Parent(field) {
static_assert((isConst == ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
//! human-readable field map type
template<class FieldCollection, typename T, Dim_t order, Dim_t dim, bool ConstField>
std::string
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::info_string() const {
std::stringstream info;
info << "Tensor(" << typeid(T).name() << ", " << order
<< "_o, " << dim << "_d)";
return info.str();
}
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator[](size_type index) {
auto && lambda = [this, &index](auto&&...tens_sizes) {
return reference(this->get_ptr_to_entry(index), tens_sizes...);
};
return call_sizes<order, dim>(lambda);
}
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
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<order, dim>(lambda);
}
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
const_reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
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<T*>(this->get_ptr_to_entry(index)),
tens_sizes...);
};
return call_sizes<order, dim>(lambda);
}
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
const_reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator[](const Ccoord & ccoord) const {
auto && index = this->collection.get_index(ccoord);
auto && lambda = [this, &index](auto&&...sizes) {
return const_reference(const_cast<T*>(this->get_ptr_to_entry(index)),
sizes...);
};
return call_sizes<order, dim>(lambda);
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
TensorFieldMap<FieldCollection, T, order, dim, ConstField> &
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator=(const T_t & val) {
for (auto && tens: *this) {
tens = val;
}
return *this;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::T_t
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
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<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::pointer
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
ptr_to_val_t(size_type index) {
auto && lambda = [this, &index](auto&&... tens_sizes) {
return std::make_unique<value_type>(this->get_ptr_to_entry(index),
tens_sizes...);
};
return call_sizes<order, dim>(lambda);
}
} // muSpectre
#endif /* FIELD_MAP_TENSOR_H */
diff --git a/src/fft/fft_engine_base.hh b/src/fft/fft_engine_base.hh
index ae9fb54..5888384 100644
--- a/src/fft/fft_engine_base.hh
+++ b/src/fft/fft_engine_base.hh
@@ -1,119 +1,119 @@
/**
* file fft_engine_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 Dec 2017
*
* @brief Interface for FFT engines
*
* @section LICENCE
*
* 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 FFT_ENGINE_BASE_H
#define FFT_ENGINE_BASE_H
#include "common/common.hh"
#include "common/field_collection.hh"
namespace muSpectre {
template <Dim_t DimS, Dim_t DimM>
class FFT_Engine_base
{
public:
constexpr static Dim_t sdim{DimS};
constexpr static Dim_t mdim{DimM};
using Ccoord = Ccoord_t<DimS>;
using Rcoord = std::array<Real, DimS>;
using GFieldCollection_t = FieldCollection<DimS, DimM, true>;
using LFieldCollection_t = FieldCollection<DimS, DimM, false>;
using Field_t = TensorField<GFieldCollection_t, Real, 2, DimM>;
using Workspace_t = TensorField<LFieldCollection_t, Complex, 2, DimM>;
using iterator = typename LFieldCollection_t::iterator;
//! Default constructor
FFT_Engine_base() = delete;
//! Constructor with system resolutions
FFT_Engine_base(Ccoord resolutions, Rcoord lengths);
//! Copy constructor
FFT_Engine_base(const FFT_Engine_base &other) = delete;
//! Move constructor
- FFT_Engine_base(FFT_Engine_base &&other) noexcept = default;
+ FFT_Engine_base(FFT_Engine_base &&other) = default;
//! Destructor
- virtual ~FFT_Engine_base() noexcept = default;
+ virtual ~FFT_Engine_base() = default;
//! Copy assignment operator
FFT_Engine_base& operator=(const FFT_Engine_base &other) = delete;
//! Move assignment operator
- FFT_Engine_base& operator=(FFT_Engine_base &&other) noexcept = default;
+ FFT_Engine_base& operator=(FFT_Engine_base &&other) = default;
// compute the plan, etc
virtual void initialise(FFT_PlanFlags /*plan_flags*/);
//! forward transform (dummy for interface)
virtual Workspace_t & fft(Field_t & /*field*/) = 0;
//! inverse transform (dummy for interface)
virtual void ifft(Field_t & /*field*/) const = 0;
/**
* iterators over only thos pixels that exist in frequency space
* (i.e. about half of all pixels, see rfft)
*/
inline iterator begin() {return this->work_space_container.begin();}
inline iterator end() {return this->work_space_container.end();}
//! nb of pixels (mostly for debugging)
size_t size() const;
size_t workspace_size() const;
//!
const Ccoord & get_resolutions() const {return this->resolutions;}
const Rcoord & get_lengths() const {return this->lengths;}
LFieldCollection_t & get_field_collection() {
return this->work_space_container;}
//! factor by which to multiply projection before inverse transform (this is
//! typically 1/nb_pixels for so-called unnormalized transforms (see,
//! e.g. http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
//! or https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html
//! . Rather than scaling the inverse transform (which would cost one more
//! loop), FFT engines provide this value so it can be used in the
//! projection operator (where no additional loop is required)
inline Real normalisation() const {return norm_factor;};
protected:
LFieldCollection_t work_space_container{};
const Ccoord resolutions;
const Rcoord lengths;
Workspace_t & work;
const Real norm_factor;
private:
};
} // muSpectre
#endif /* FFT_ENGINE_BASE_H */
diff --git a/src/fft/fft_utils.hh b/src/fft/fft_utils.hh
index 6832f85..3d44776 100644
--- a/src/fft/fft_utils.hh
+++ b/src/fft/fft_utils.hh
@@ -1,125 +1,125 @@
/**
* file fft_utils.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 06 Dec 2017
*
* @brief collection of functions used in the context of spectral operations
*
* @section LICENCE
*
* 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 FFT_UTILS_H
#define FFT_UTILS_H
#include "common/common.hh"
#include <Eigen/Dense>
#include <valarray>
#include <array>
namespace muSpectre {
/**
* compute fft frequencies (in time (or length) units of of sampling
* periods), see numpy's fftfreq function for reference
*/
std::valarray<Real> fft_freqs(size_t nb_samples);
/**
* compute fft frequencies in correct length or time units. Here,
* length refers to the total size of the domain over which the fft
* is taken (for instance the length of an edge of an RVE)
*/
std::valarray<Real> fft_freqs(size_t nb_samples, Real length);
/**
* Get fft_freqs for a grid
*/
template <size_t dim>
inline std::array<std::valarray<Real>, dim> fft_freqs(Ccoord_t<dim> sizes,
std::array<Real, dim> lengths) {
std::array<std::valarray<Real>, dim> retval{};
for (size_t i = 0; i < dim; ++i) {
retval[i] = std::move(fft_freqs(sizes[i], lengths[i]));
}
return retval;
}
/**
* simple class encapsulating the creation, and retrieval of
* wave vectors
*/
template <Dim_t dim>
class FFT_freqs
{
public:
using Vector = Eigen::Matrix<Real, dim, 1>;
//! Default constructor
FFT_freqs() = delete;
//! constructor with problem sizes
FFT_freqs(Ccoord_t<dim> sizes, std::array<Real, dim> lengths)
: freqs{fft_freqs(sizes, lengths)} {}
//! Copy constructor
FFT_freqs(const FFT_freqs &other) = delete;
//! Move constructor
FFT_freqs(FFT_freqs &&other) = default;
//! Destructor
- virtual ~FFT_freqs() noexcept = default;
+ virtual ~FFT_freqs() = default;
//! Copy assignment operator
FFT_freqs& operator=(const FFT_freqs &other) = delete;
//! Move assignment operator
FFT_freqs& operator=(FFT_freqs &&other) = default;
//! get unnormalised wave vector (in sampling units)
inline Vector get_xi(const Ccoord_t<dim> ccoord) const;
//! get normalised wave vector
inline Vector get_unit_xi(const Ccoord_t<dim> ccoord) const{
auto && xi = this->get_xi(std::move(ccoord));
return xi/xi.norm();
}
protected:
const std::array<std::valarray<Real>, dim> freqs;
private:
};
template <Dim_t dim>
typename FFT_freqs<dim>::Vector
FFT_freqs<dim>::get_xi(const Ccoord_t<dim> ccoord) const {
Vector retval{};
for (Dim_t i = 0; i < dim; ++i) {
retval(i) = this->freqs[i][ccoord[i]];
}
return retval;
}
} // muSpectre
#endif /* FFT_UTILS_H */
diff --git a/src/fft/projection_base.hh b/src/fft/projection_base.hh
index 7353750..068992c 100644
--- a/src/fft/projection_base.hh
+++ b/src/fft/projection_base.hh
@@ -1,106 +1,106 @@
/**
* file projection_base.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 03 Dec 2017
*
* @brief Base class for Projection operators
*
* @section LICENCE
*
* 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 PROJECTION_BASE_H
#define PROJECTION_BASE_H
#include "common/common.hh"
#include "common/field_collection.hh"
#include "common/field.hh"
#include "fft/fft_engine_base.hh"
#include <memory>
namespace muSpectre {
template<class Projection>
struct Projection_traits {
};
template <Dim_t DimS, Dim_t DimM>
class ProjectionBase
{
public:
using FFT_Engine = FFT_Engine_base<DimS, DimM>;
using FFT_Engine_ptr = std::unique_ptr<FFT_Engine>;
using Ccoord = typename FFT_Engine::Ccoord;
using Rcoord = typename FFT_Engine::Rcoord;
using GFieldCollection_t = typename FFT_Engine::GFieldCollection_t;
using LFieldCollection_t = typename FFT_Engine::LFieldCollection_t;
using Field_t = typename FFT_Engine::Field_t;
using iterator = typename FFT_Engine::iterator;
//! Default constructor
ProjectionBase() = delete;
//! Constructor with system sizes
ProjectionBase(FFT_Engine_ptr engine, Formulation form);
//! Copy constructor
ProjectionBase(const ProjectionBase &other) = delete;
//! Move constructor
- ProjectionBase(ProjectionBase &&other) noexcept = default;
+ ProjectionBase(ProjectionBase &&other) = default;
//! Destructor
- virtual ~ProjectionBase() noexcept = default;
+ virtual ~ProjectionBase() = default;
//! Copy assignment operator
ProjectionBase& operator=(const ProjectionBase &other) = delete;
//! Move assignment operator
- ProjectionBase& operator=(ProjectionBase &&other) noexcept = default;
+ ProjectionBase& operator=(ProjectionBase &&other) = default;
//! initialises the fft engine (plan the transform)
virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
//! apply the projection operator to a field
virtual void apply_projection(Field_t & field) = 0;
//!
const Ccoord & get_resolutions() const {
return this->fft_engine->get_resolutions();}
const Rcoord & get_lengths() const {
return this->fft_engine->get_lengths();}
const Formulation & get_formulation() const {return this->form;}
protected:
FFT_Engine_ptr fft_engine;
const Formulation form;
LFieldCollection_t & projection_container{};
private:
};
} // muSpectre
#endif /* PROJECTION_BASE_H */
diff --git a/src/fft/projection_finite_strain.hh b/src/fft/projection_finite_strain.hh
index 309eabd..ff3bec9 100644
--- a/src/fft/projection_finite_strain.hh
+++ b/src/fft/projection_finite_strain.hh
@@ -1,86 +1,86 @@
/**
* file projection_finite_strain.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Dec 2017
*
* @brief Class for standard finite-strain gradient projections see de Geus et
* al. (https://doi.org/10.1016/j.cma.2016.12.032) for derivation
*
* @section LICENCE
*
* 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 PROJECTION_FINITE_STRAIN_H
#define PROJECTION_FINITE_STRAIN_H
#include "fft/projection_default.hh"
#include "common/common.hh"
#include "common/field_collection.hh"
#include "common/field_map.hh"
namespace muSpectre {
template <Dim_t DimS, Dim_t DimM>
class ProjectionFiniteStrain: public ProjectionDefault<DimS, DimM>
{
public:
using Parent = ProjectionDefault<DimS, DimM>;
using FFT_Engine_ptr = typename Parent::FFT_Engine_ptr;
using Ccoord = typename Parent::Ccoord;
//using GFieldCollection_t = FieldCollection<DimS, DimM, true>;
using LFieldCollection_t = FieldCollection<DimS, DimM, false>;
//using Field_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
using Proj_t = TensorField<LFieldCollection_t, Real, fourthOrder, DimM>;
using Proj_map = T4MatrixFieldMap<LFieldCollection_t, Real, DimM>;
using Vector_map = MatrixFieldMap<LFieldCollection_t, Complex, DimM*DimM, 1>;
//! Default constructor
ProjectionFiniteStrain() = delete;
//! Constructor with fft_engine
ProjectionFiniteStrain(FFT_Engine_ptr engine);
//! Copy constructor
ProjectionFiniteStrain(const ProjectionFiniteStrain &other) = delete;
//! Move constructor
ProjectionFiniteStrain(ProjectionFiniteStrain &&other) = default;
//! Destructor
- virtual ~ProjectionFiniteStrain() noexcept = default;
+ virtual ~ProjectionFiniteStrain() = default;
//! Copy assignment operator
ProjectionFiniteStrain& operator=(const ProjectionFiniteStrain &other) = delete;
//! Move assignment operator
ProjectionFiniteStrain& operator=(ProjectionFiniteStrain &&other) = default;
//! initialises the fft engine (plan the transform)
virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final;
protected:
private:
};
} // muSpectre
#endif /* PROJECTION_FINITE_STRAIN_H */
diff --git a/src/fft/projection_finite_strain_fast.hh b/src/fft/projection_finite_strain_fast.hh
index a183343..fbe5f84 100644
--- a/src/fft/projection_finite_strain_fast.hh
+++ b/src/fft/projection_finite_strain_fast.hh
@@ -1,90 +1,90 @@
/**
* file projection_finite_strain_fast.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 12 Dec 2017
*
* @brief Faster alternative to ProjectionFinitestrain
*
* @section LICENCE
*
* 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 PROJECTION_FINITE_STRAIN_FAST_H
#define PROJECTION_FINITE_STRAIN_FAST_H
#include "fft/projection_base.hh"
#include "common/common.hh"
#include "common/field_collection.hh"
#include "common/field_map.hh"
namespace muSpectre {
template <Dim_t DimS, Dim_t DimM>
class ProjectionFiniteStrainFast: public ProjectionBase<DimS, DimM>
{
public:
using Parent = ProjectionBase<DimS, DimM>;
using FFT_Engine_ptr = typename Parent::FFT_Engine_ptr;
using Ccoord = typename Parent::Ccoord;
using GFieldCollection_t = FieldCollection<DimS, DimM, true>;
using LFieldCollection_t = FieldCollection<DimS, DimM, false>;
using Field_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
using Proj_t = TensorField<LFieldCollection_t, Real, firstOrder, DimM>;
using Proj_map = MatrixFieldMap<LFieldCollection_t, Real, DimM, 1>;
using Grad_map = MatrixFieldMap<LFieldCollection_t, Complex, DimM, DimM>;
//! Default constructor
ProjectionFiniteStrainFast() = delete;
//! Constructor with fft_engine
ProjectionFiniteStrainFast(FFT_Engine_ptr engine);
//! Copy constructor
ProjectionFiniteStrainFast(const ProjectionFiniteStrainFast &other) = delete;
//! Move constructor
- ProjectionFiniteStrainFast(ProjectionFiniteStrainFast &&other) noexcept = default;
+ ProjectionFiniteStrainFast(ProjectionFiniteStrainFast &&other) = default;
//! Destructor
- virtual ~ProjectionFiniteStrainFast() noexcept = default;
+ virtual ~ProjectionFiniteStrainFast() = default;
//! Copy assignment operator
ProjectionFiniteStrainFast& operator=(const ProjectionFiniteStrainFast &other) = delete;
//! Move assignment operator
ProjectionFiniteStrainFast& operator=(ProjectionFiniteStrainFast &&other) = default;
//! initialises the fft engine (plan the transform)
virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final;
//! apply the projection operator to a field
void apply_projection(Field_t & field) override final;
protected:
Proj_map xis;
private:
};
} // muSpectre
#endif /* PROJECTION_FINITE_STRAIN_FAST_H */
diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh
index ff29b20..1cb5ec1 100644
--- a/src/materials/material_base.hh
+++ b/src/materials/material_base.hh
@@ -1,142 +1,142 @@
/**
* file material_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 25 Oct 2017
*
* @brief Base class for materials (constitutive models)
*
* @section LICENCE
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef MATERIAL_BASE_H
#define MATERIAL_BASE_H
#include "common/common.hh"
#include "common/field.hh"
#include "common/field_collection.hh"
#include <string>
namespace muSpectre {
//! DimS spatial dimension (dimension of problem
//! DimM material_dimension (dimension of constitutive law)
template<Dim_t DimS, Dim_t DimM>
class MaterialBase
{
public:
//! typedefs for data handled by this interface
//! global field collection for system-wide fields, like stress, strain, etc
using GFieldCollection_t = FieldCollection<DimS, DimM, true>;
//! field collection for internal variables, such as eigen-strains,
//! plastic strains, damage variables, etc, but also for managing which
//! pixels the material is responsible for
using MFieldCollection_t = FieldCollection<DimS, DimM, false>;
using iterator = typename MFieldCollection_t::iterator;
using Field_t = internal::FieldBase<GFieldCollection_t>;
using StressField_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
using StrainField_t = StressField_t;
using TangentField_t = TensorField<GFieldCollection_t, Real, fourthOrder, DimM>;
using Ccoord = Ccoord_t<DimS>;
//! Default constructor
MaterialBase() = delete;
//! Construct by name
MaterialBase(std::string name);
//! Copy constructor
MaterialBase(const MaterialBase &other) = delete;
//! Move constructor
- MaterialBase(MaterialBase &&other) noexcept = delete;
+ MaterialBase(MaterialBase &&other) = delete;
//! Destructor
- virtual ~MaterialBase() noexcept = default;
+ virtual ~MaterialBase() = default;
//! Copy assignment operator
MaterialBase& operator=(const MaterialBase &other) = delete;
//! Move assignment operator
- MaterialBase& operator=(MaterialBase &&other) noexcept = delete;
+ MaterialBase& operator=(MaterialBase &&other) = delete;
/**
* take responsibility for a pixel identified by its cell coordinates
* WARNING: this won't work for materials with additional info per pixel
* (as, e.g. for eigenstrain), we need to pass more parameters. Materials
* of this tye need to overload add_pixel
*/
void add_pixel(const Ccoord & ccooord);
//! allocate memory, etc, but also: wipe history variables!
virtual void initialise(bool stiffness = false) = 0;
//! return the materil's name
const std::string & get_name() const;
//! for static inheritance stuff
constexpr static Dim_t sdim() {return DimS;}
constexpr static Dim_t mdim() {return DimM;}
//! computes stress
virtual void compute_stresses(const StrainField_t & F,
StressField_t & P,
Formulation form) = 0;
/**
* Convenience function to compute stresses, mostly for debugging and
* testing. Has runtime-cost associated with compatibility-checking and
* conversion of the Field_t arguments that can be avoided by using the
* version with strongly typed field references
*/
void compute_stresses(const Field_t & F,
Field_t & P,
Formulation form);
//! computes stress and tangent moduli
virtual void compute_stresses_tangent(const StrainField_t & F,
StressField_t & P,
TangentField_t & K,
Formulation form) = 0;
/**
* Convenience function to compute stresses and tangent moduli, mostly for
* debugging and testing. Has runtime-cost associated with
* compatibility-checking and conversion of the Field_t arguments that can
* be avoided by using the version with strongly typed field references
*/
void compute_stresses_tangent(const Field_t & F,
Field_t & P,
Field_t & K,
Formulation form);
inline iterator begin() {return this->internal_fields.begin();}
inline iterator end() {return this->internal_fields.end();}
protected:
//! members
const std::string name;
MFieldCollection_t internal_fields{};
private:
};
} // muSpectre
#endif /* MATERIAL_BASE_H */
diff --git a/src/materials/material_hyper_elastic1.hh b/src/materials/material_hyper_elastic1.hh
index 2fe91f2..30255c2 100644
--- a/src/materials/material_hyper_elastic1.hh
+++ b/src/materials/material_hyper_elastic1.hh
@@ -1,133 +1,133 @@
/**
* file material_hyper_elastic1.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 13 Nov 2017
*
* @brief Implementation for hyperelastic reference material like in de Geus
* 2017. This follows the simplest and likely not most efficient
* implementation (with exception of the Python law)
*
* @section LICENCE
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef MATERIAL_HYPER_ELASTIC1_H
#define MATERIAL_HYPER_ELASTIC1_H
#include "common/common.hh"
#include "materials/material_muSpectre_base.hh"
namespace muSpectre {
template<Dim_t DimS, Dim_t DimM>
class MaterialHyperElastic1;
template <Dim_t DimS, Dim_t DimM>
struct MaterialMuSpectre_traits<MaterialHyperElastic1<DimS, DimM>>:
public MaterialMuSpectre_traits<void> {
using Parent = MaterialMuSpectre_traits<void>;
using InternalVariables = typename Parent::DefaultInternalVariables;
};
//! DimS spatial dimension (dimension of problem
//! DimM material_dimension (dimension of constitutive law)
template<Dim_t DimS, Dim_t DimM>
class MaterialHyperElastic1:
public MaterialMuSpectre<MaterialHyperElastic1<DimS, DimM>, DimS, DimM>
{
public:
using Parent = MaterialMuSpectre<MaterialHyperElastic1, DimS, DimM>;
using NeedTangent = typename Parent::NeedTangent;
using GFieldCollection_t = typename Parent::GFieldCollection_t;
// declare what type of strain measure your law takes as input
constexpr static auto strain_measure{StrainMeasure::GreenLagrange};
// declare what type of stress measure your law yields as output
constexpr static auto stress_measure{StressMeasure::PK2};
// declare whether the derivative of stress with respect to strain is uniform
constexpr static bool uniform_stiffness = true;
// declare the type in which you wish to receive your strain measure
using StrainMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM, true>;
using StressMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM>;
using TangentMap_t = T4MatrixFieldMap<GFieldCollection_t, Real, DimM>;
using Strain_t = typename StrainMap_t::const_reference;
using Stress_t = typename StressMap_t::reference;
using Tangent_t = typename TangentMap_t::reference;
using Stiffness_t = Eigen::TensorFixedSize
<Real, Eigen::Sizes<DimM, DimM, DimM, DimM>>;
using InternalVariables = typename Parent::DefaultInternalVariables;
//! Default constructor
MaterialHyperElastic1() = delete;
//! Copy constructor
MaterialHyperElastic1(const MaterialHyperElastic1 &other) = delete;
//! Construct by name, Young's modulus and Poisson's ratio
MaterialHyperElastic1(std::string name, Real young, Real poisson);
//! Move constructor
- MaterialHyperElastic1(MaterialHyperElastic1 &&other) noexcept = delete;
+ MaterialHyperElastic1(MaterialHyperElastic1 &&other) = delete;
//! Destructor
- virtual ~MaterialHyperElastic1() noexcept = default;
+ virtual ~MaterialHyperElastic1() = default;
//! Copy assignment operator
MaterialHyperElastic1& operator=(const MaterialHyperElastic1 &other) = delete;
//! Move assignment operator
- MaterialHyperElastic1& operator=(MaterialHyperElastic1 &&other) noexcept = delete;
+ MaterialHyperElastic1& operator=(MaterialHyperElastic1 &&other) = delete;
template <class s_t>
inline decltype(auto) evaluate_stress(s_t && E);
template <class s_t>
inline decltype(auto) evaluate_stress_tangent(s_t && E);
const Tangent_t & get_stiffness() const;
protected:
const Real young, poisson, lambda, mu;
const Stiffness_t C;
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
template <class s_t>
decltype(auto)
MaterialHyperElastic1<DimS, DimM>::evaluate_stress(s_t && E) {
return E.trace()*lambda * Strain_t::Identity() + 2*mu*E;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
template <class s_t>
decltype(auto)
MaterialHyperElastic1<DimS, DimM>::evaluate_stress_tangent(s_t && E) {
return std::make_tuple(std::move(this->evaluate_stress(std::move(E))),
std::move(Tangent_t(const_cast<double*>(this->C.data()))));
}
} // muSpectre
#endif /* MATERIAL_HYPER_ELASTIC1_H */
diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh
index 550bb9d..49ace52 100644
--- a/src/materials/material_muSpectre_base.hh
+++ b/src/materials/material_muSpectre_base.hh
@@ -1,663 +1,663 @@
/**
* file material_muSpectre_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 25 Oct 2017
*
* @brief Base class for materials written for µSpectre specifically. These
* can take full advantage of the configuration-change utilities of
* µSpectre. The user can inherit from them to define new constitutive
* laws and is merely required to provide the methods for computing the
* second Piola-Kirchhoff stress and Tangent. This class uses the
* "curiously recurring template parameter" to avoid virtual calls.
*
* @section LICENCE
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef MATERIAL_MUSPECTRE_BASE_H
#define MATERIAL_MUSPECTRE_BASE_H
#include "common/common.hh"
#include "materials/material_base.hh"
#include "materials/materials_toolbox.hh"
#include "common/field_collection.hh"
#include "common/field.hh"
#include "common//utilities.hh"
#include <tuple>
#include <type_traits>
#include <iterator>
#include <stdexcept>
namespace muSpectre {
/**
* Forward declaration for factory function
*/
template <Dim_t DimS, Dim_t DimM>
class SystemBase;
template <class Material>
struct MaterialMuSpectre_traits {
};
template <class Material, Dim_t DimS, Dim_t DimM>
class MaterialMuSpectre;
template <>
struct MaterialMuSpectre_traits<void> {
using DefaultInternalVariables = std::tuple<>;
};
//! 'Material' is a CRTP
template <class Material, Dim_t DimS, Dim_t DimM>
class MaterialMuSpectre: public MaterialBase<DimS, DimM>
{
public:
using NeedTangent = MatTB::NeedTangent;
using Parent = MaterialBase<DimS, DimM>;
using GFieldCollection_t = typename Parent::GFieldCollection_t;
using MFieldCollection_t = typename Parent::MFieldCollection_t;
using StressField_t = typename Parent::StressField_t;
using StrainField_t = typename Parent::StrainField_t;
using TangentField_t = typename Parent::TangentField_t;
using DefaultInternalVariables = std::tuple<>;
using traits = MaterialMuSpectre_traits<Material>;
//! Default constructor
MaterialMuSpectre() = delete;
//! Construct by name
MaterialMuSpectre(std::string name);
//! Copy constructor
MaterialMuSpectre(const MaterialMuSpectre &other) = delete;
//! Move constructor
- MaterialMuSpectre(MaterialMuSpectre &&other) noexcept = delete;
+ MaterialMuSpectre(MaterialMuSpectre &&other) = delete;
//! Destructor
- virtual ~MaterialMuSpectre() noexcept = default;
+ virtual ~MaterialMuSpectre() = default;
//! Factory
template <class... ConstructorArgs>
static Material & make(SystemBase<DimS, DimM> & sys,
ConstructorArgs &&... args);
//! Copy assignment operator
MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete;
//! Move assignment operator
- MaterialMuSpectre& operator=(MaterialMuSpectre &&other) noexcept = delete;
+ MaterialMuSpectre& operator=(MaterialMuSpectre &&other) = delete;
//* allocate memory, etc
virtual void initialise(bool stiffness = false) override final;
using Parent::compute_stresses;
using Parent::compute_stresses_tangent;
//! computes stress
virtual void compute_stresses(const StrainField_t & F,
StressField_t & P,
Formulation form) override final;
//! computes stress and tangent modulus
virtual void compute_stresses_tangent(const StrainField_t & F,
StressField_t & P,
TangentField_t & K,
Formulation form) override final;
protected:
//! computes stress with the formulation available at compile time
//! __attribute__ required by g++-6 and g++-7 because of this bug:
//! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947
template <Formulation Form>
inline void compute_stresses_worker(const StrainField_t & F,
StressField_t & P)
__attribute__ ((visibility ("default")));
//! computes stress with the formulation available at compile time
//! __attribute__ required by g++-6 and g++-7 because of this bug:
//! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947
template <Formulation Form>
inline void compute_stresses_worker(const StrainField_t & F,
StressField_t & P,
TangentField_t & K)
__attribute__ ((visibility ("default")));
//! this iterable class is a default for simple laws that just take a strain
//! the iterable is just a templated wrapper to provide a range to iterate over
//! that does or does not include tangent moduli
template<NeedTangent need_tgt = NeedTangent::no>
class iterable_proxy;
/**
* inheriting classes with internal variables need to overload this function
*/
decltype(auto) get_internals() {
// the default material has no internal variables
return typename Material::InternalVariables{};}
decltype(auto) get_internals() const {
// the default material has no internal variables
return typename Material::InternalVariables{};}
typename traits::InternalVariables internal_variables{};
bool is_initialised{false};
private:
};
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
MaterialMuSpectre<Material, DimS, DimM>::
MaterialMuSpectre(std::string name)
:Parent(name) {
using stress_compatible = typename Material::StressMap_t::
template is_compatible<StressField_t>;
using strain_compatible = typename Material::StrainMap_t::
template is_compatible<StrainField_t>;
using tangent_compatible = typename Material::TangentMap_t::
template is_compatible<TangentField_t>;
static_assert((stress_compatible::value &&
stress_compatible::explain()),
"The material's declared stress map is not compatible "
"with the stress field. More info in previously shown "
"assert.");
static_assert((strain_compatible::value &&
strain_compatible::explain()),
"The material's declared strain map is not compatible "
"with the strain field. More info in previously shown "
"assert.");
static_assert((tangent_compatible::value &&
tangent_compatible::explain()),
"The material's declared tangent map is not compatible "
"with the tangent field. More info in previously shown "
"assert.");
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <class... ConstructorArgs>
Material & MaterialMuSpectre<Material, DimS, DimM>::
make(SystemBase<DimS, DimM> & sys,
ConstructorArgs && ... args) {
auto mat = std::make_unique<Material>(args...);
auto & mat_ref = *mat;
sys.add_material(std::move(mat));
return mat_ref;
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
initialise(bool /*stiffness*/) {
if (!this->is_initialised) {
this->internal_fields.initialise();
this->is_initialised = true;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses(const StrainField_t &F, StressField_t &P,
Formulation form) {
switch (form) {
case Formulation::finite_strain: {
this->template compute_stresses_worker<Formulation::finite_strain>(F, P);
break;
}
case Formulation::small_strain: {
this->template compute_stresses_worker<Formulation::small_strain>(F, P);
break;
}
default:
throw std::runtime_error("Unknown formulation");
break;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses_tangent(const StrainField_t & F, StressField_t & P,
TangentField_t & K,
Formulation form) {
switch (form) {
case Formulation::finite_strain: {
this->template compute_stresses_worker<Formulation::finite_strain>(F, P, K);
break;
}
case Formulation::small_strain: {
this->template compute_stresses_worker<Formulation::small_strain>(F, P, K);
break;
}
default:
throw std::runtime_error("Unknown formulation");
break;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <Formulation Form>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses_worker(const StrainField_t & F,
StressField_t & P,
TangentField_t & K){
/* These lambdas are executed for every integration point.
F contains the transformation gradient for finite strain calculations and
the infinitesimal strain tensor in small strain problems
The internal_variables tuple contains whatever internal variables
Material declared (e.g., eigenstrain, strain rate, etc.)
*/
using Strains_t = std::tuple<typename Material::Strain_t>;
using Stresses_t = std::tuple<typename Material::Stress_t,
typename Material::Tangent_t>;
auto constitutive_law_small_strain = [this]
(Strains_t Strains, Stresses_t Stresses, auto && internal_variables) {
constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
constexpr StrainMeasure expected_strain_m{
get_formulation_strain_type(Form, Material::strain_measure)};
auto & this_mat = static_cast<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && F = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
// return value contains a tuple of rvalue_refs to both stress and tangent moduli
Stresses =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress_tangent(std::move(strain),
internals...);},
internal_variables);
};
auto constitutive_law_finite_strain = [this]
(Strains_t Strains, Stresses_t Stresses, auto && internal_variables) {
constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
constexpr StrainMeasure expected_strain_m{
get_formulation_strain_type(Form, Material::strain_measure)};
auto & this_mat = static_cast<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && grad = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(grad);
// TODO: Figure this out: I can't std::move(internals...),
// because if there are no internals, compilation fails with "no
// matching function for call to ‘move()’'. These are tuples of
// lvalue references, so it shouldn't be too bad, but still
// irksome.
// return value contains a tuple of rvalue_refs to both stress
// and tangent moduli
auto stress_tgt =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress_tangent(std::move(strain),
internals...);},
internal_variables);
auto && stress = std::get<0>(stress_tgt);
auto && tangent = std::get<1>(stress_tgt);
Stresses = MatTB::PK1_stress<Material::stress_measure, Material::strain_measure>
(std::move(grad),
std::move(stress),
std::move(tangent));
};
iterable_proxy<NeedTangent::yes> fields{*this, F, P, K};
for (auto && arglist: fields) {
/**
* arglist is a tuple of three tuples containing only Lvalue
* references (see value_tye in the class definition of
* iterable_proxy::iterator). Tuples contain strains, stresses
* and internal variables, respectively,
*/
//auto && stress_tgt = std::get<0>(tuples);
//auto && inputs = std::get<1>(tuples);TODO:clean this
static_assert(std::is_same<typename Material::Strain_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<0>(arglist)))>>::value,
"Type mismatch for strain reference, check iterator "
"value_type");
static_assert(std::is_same<typename Material::Stress_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<1>(arglist)))>>::value,
"Type mismatch for stress reference, check iterator"
"value_type");
static_assert(std::is_same<typename Material::Tangent_t,
std::remove_reference_t<
decltype(std::get<1>(std::get<1>(arglist)))>>::value,
"Type mismatch for tangent reference, check iterator"
"value_type");
switch (Form) {
case Formulation::small_strain: {
apply(constitutive_law_small_strain, std::move(arglist));
break;
}
case Formulation::finite_strain: {
apply(constitutive_law_finite_strain, std::move(arglist));
break;
}
}
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <Formulation Form>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses_worker(const StrainField_t & F,
StressField_t & P){
/* These lambdas are executed for every integration point.
F contains the transformation gradient for finite strain calculations and
the infinitesimal strain tensor in small strain problems
The internal_variables tuple contains whatever internal variables
Material declared (e.g., eigenstrain, strain rate, etc.)
*/
using Strains_t = std::tuple<typename Material::Strain_t>;
using Stresses_t = std::tuple<typename Material::Stress_t>;
auto constitutive_law_small_strain = [this]
(Strains_t Strains, Stresses_t Stresses, auto && internal_variables) {
constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
constexpr StrainMeasure expected_strain_m{
get_formulation_strain_type(Form, Material::strain_measure)};
auto & this_mat = static_cast<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && F = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
// return value contains a tuple of rvalue_refs to both stress and tangent moduli
auto && sigma = std::get<0>(Stresses);
sigma =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress(std::move(strain),
internals...);},
internal_variables);
};
auto constitutive_law_finite_strain = [this]
(Strains_t Strains, Stresses_t && Stresses, auto && internal_variables) {
constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
constexpr StrainMeasure expected_strain_m{
get_formulation_strain_type(Form, Material::strain_measure)};
auto & this_mat = static_cast<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && F = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
// TODO: Figure this out: I can't std::move(internals...),
// because if there are no internals, compilation fails with "no
// matching function for call to ‘move()’'. These are tuples of
// lvalue references, so it shouldn't be too bad, but still
// irksome.
// return value contains a tuple of rvalue_refs to both stress
// and tangent moduli
auto && stress =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress(std::move(strain),
internals...);},
internal_variables);
auto && P = get<0>(Stresses);
P = MatTB::PK1_stress<Material::stress_measure, Material::strain_measure>
(F, stress);
};
iterable_proxy<NeedTangent::no> fields{*this, F, P};
for (auto && arglist: fields) {
/**
* arglist is a tuple of three tuples containing only Lvalue
* references (see value_tye in the class definition of
* iterable_proxy::iterator). Tuples contain strains, stresses
* and internal variables, respectively,
*/
//auto && stress_tgt = std::get<0>(tuples);
//auto && inputs = std::get<1>(tuples);TODO:clean this
static_assert(std::is_same<typename Material::Strain_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<0>(arglist)))>>::value,
"Type mismatch for strain reference, check iterator "
"value_type");
static_assert(std::is_same<typename Material::Stress_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<1>(arglist)))>>::value,
"Type mismatch for stress reference, check iterator"
"value_type");
switch (Form) {
case Formulation::small_strain: {
apply(constitutive_law_small_strain, std::move(arglist));
break;
}
case Formulation::finite_strain: {
apply(constitutive_law_finite_strain, std::move(arglist));
break;
}
}
}
}
/* ---------------------------------------------------------------------- */
//! this iterator class is a default for simple laws that just take a strain
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
class MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy {
public:
//! Default constructor
iterable_proxy() = delete;
using NeedTangent =
typename MaterialMuSpectre<Material, DimS, DimM>::NeedTangent;
/** Iterator uses the material's internal variables field
collection to iterate selectively over the global fields
(such as the transformation gradient F and first
Piola-Kirchhoff stress P.
**/
template<bool DoNeedTgt=(NeedTgt == NeedTangent::yes)>
iterable_proxy(const MaterialMuSpectre & mat,
const StrainField_t & F,
StressField_t & P,
std::enable_if_t<DoNeedTgt, TangentField_t> & K)
:material(mat), strain_field(F), stress_tup{P,K},
internals(material.internal_variables){};
template<bool DontNeedTgt=(NeedTgt == NeedTangent::no)>
iterable_proxy(const MaterialMuSpectre & mat,
const StrainField_t & F,
std::enable_if_t<DontNeedTgt, StressField_t> & P)
:material(mat), strain_field(F), stress_tup{P},
internals(material.internal_variables){};
using StrainMap_t = typename Material::StrainMap_t;
using StressMap_t = typename Material::StressMap_t;
using TangentMap_t = typename Material::TangentMap_t;
using Strain_t = typename Material::Strain_t;
using Stress_t = typename Material::Stress_t;
using Tangent_t = typename Material::Tangent_t;
using InternalVariables = typename Material::InternalVariables;
using StressFieldTup = std::conditional_t
<(NeedTgt == NeedTangent::yes),
std::tuple<StressField_t&, TangentField_t&>,
std::tuple<StressField_t&>>;
using StressMapTup = std::conditional_t
<(NeedTgt == NeedTangent::yes),
std::tuple<StressMap_t, TangentMap_t>,
std::tuple<StressMap_t>>;
using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes),
std::tuple<Stress_t, Tangent_t>,
std::tuple<Stress_t>>;
//! Copy constructor
iterable_proxy(const iterable_proxy &other) = default;
//! Move constructor
- iterable_proxy(iterable_proxy &&other) noexcept = default;
+ iterable_proxy(iterable_proxy &&other) = default;
//! Destructor
- virtual ~iterable_proxy() noexcept = default;
+ virtual ~iterable_proxy() = default;
//! Copy assignment operator
iterable_proxy& operator=(const iterable_proxy &other) = default;
//! Move assignment operator
iterable_proxy& operator=(iterable_proxy &&other) = default;
class iterator
{
public:
using InternalReferences = MatTB::ReferenceTuple_t<InternalVariables>;
using value_type =
std::tuple<std::tuple<Strain_t>, Stress_tTup, InternalReferences>;
using iterator_category = std::forward_iterator_tag;
//! Default constructor
iterator() = delete;
/** Iterator uses the material's internal variables field
collection to iterate selectively over the global fields
(such as the transformation gradient F and first
Piola-Kirchhoff stress P.
**/
iterator(const iterable_proxy & it, bool begin = true)
: it{it}, strain_map{it.strain_field},
stress_map {it.stress_tup},
index{begin ? 0:it.material.internal_fields.size()}{}
//! Copy constructor
iterator(const iterator &other) = default;
//! Move constructor
- iterator(iterator &&other) noexcept = default;
+ iterator(iterator &&other) = default;
//! Destructor
- virtual ~iterator() noexcept = default;
+ 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++();
//! dereference
inline value_type operator*();
//! inequality
inline bool operator!=(const iterator & other) const;
protected:
const iterable_proxy & it;
StrainMap_t strain_map;
StressMapTup stress_map;
size_t index;
private:
};
iterator begin() {return std::move(iterator(*this));}
iterator end() {return std::move(iterator(*this, false));}
protected:
const MaterialMuSpectre & material;
const StrainField_t & strain_field;
StressFieldTup stress_tup;
const InternalVariables & internals;
private:
};
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
bool
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
operator!=(const iterator & other) const {
return (this->index != other.index);
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
typename MaterialMuSpectre<Material, DimS, DimM>::
template iterable_proxy<NeedTgt>::
iterator &
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
operator++() {
this->index++;
return *this;
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgT>
typename MaterialMuSpectre<Material, DimS, DimM>::
template iterable_proxy<NeedTgT>::iterator::
value_type
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgT>::iterator::
operator*() {
const Ccoord_t<DimS> pixel{
this->it.material.internal_fields.get_ccoord(this->index)};
auto && strain = std::make_tuple(this->strain_map[pixel]);
auto && stresses =
apply([&pixel] (auto && ... stress_tgt) {
return std::make_tuple(stress_tgt[pixel]...);},
this->stress_map);
const auto & internal = this->it.material.get_internals();
const auto id{index};
auto && internals =
apply([id] (auto && ... internals) {
return std::make_tuple(internals[id]...);},
internal);
return std::make_tuple(std::move(strain),
std::move(stresses),
std::move(internals));
}
} // muSpectre
#endif /* MATERIAL_MUSPECTRE_BASE_H */
diff --git a/src/solver/solver_base.hh b/src/solver/solver_base.hh
index f4ee868..2071da7 100644
--- a/src/solver/solver_base.hh
+++ b/src/solver/solver_base.hh
@@ -1,98 +1,98 @@
/**
* file solver_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 18 Dec 2017
*
* @brief Base class for solvers
*
* @section LICENCE
*
* 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 SOLVER_BASE_H
#define SOLVER_BASE_H
#include "solver/solver_error.hh"
#include "common/common.hh"
#include "system/system_base.hh"
#include "common/tensor_algebra.hh"
#include <Eigen/Dense>
#include <vector>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM=DimS>
class SolverBase
{
public:
enum class TangentRequirement{NoNeed, NeedEffect, NeedTangents};
using Ccoord = Ccoord_t<DimS>;
using Collection_t = GlobalFieldCollection<DimS, DimM>;
//! Default constructor
SolverBase() = delete;
//! Constructor with domain resolutions
SolverBase(Ccoord resolutions): resolutions{resolutions} {}
//! Copy constructor
SolverBase(const SolverBase &other) = delete;
//! Move constructor
SolverBase(SolverBase &&other) = default;
//! Destructor
- virtual ~SolverBase() noexcept = default;
+ virtual ~SolverBase() = default;
//! Copy assignment operator
SolverBase& operator=(const SolverBase &other) = delete;
//! Move assignment operator
SolverBase& operator=(SolverBase &&other) = default;
//! Allocate fields used during the solution
void initialise() {this->collection.initialise(resolutions);}
bool need_tangents() const {
return (this->get_tangent_req() == TangentRequirement::NeedTangents);}
bool need_effect() const {
return (this->get_tangent_req() == TangentRequirement::NeedEffect);}
bool no_need_tangent() const {
return (this->get_tangent_req() == TangentRequirement::NoNeed);}
bool has_converged() const {return this->converged;}
protected:
virtual TangentRequirement get_tangent_req() const = 0;
Ccoord resolutions;
//! storage for internal fields to avoid reallocations between calls
Collection_t collection{};
bool converged{false};
private:
};
} // muSpectre
#endif /* SOLVER_BASE_H */
diff --git a/src/solver/solver_cg.hh b/src/solver/solver_cg.hh
index 7af6578..60264cb 100644
--- a/src/solver/solver_cg.hh
+++ b/src/solver/solver_cg.hh
@@ -1,96 +1,96 @@
/**
* file solver_cg.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Dec 2017
*
* @brief class for a simple implementation of a conjugate gradient
* solver. This follows algorithm 5.2 in Nocedal's Numerical
* Optimization (p 112)
*
* @section LICENCE
*
* 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 SOLVER_CG_H
#define SOLVER_CG_H
#include "solver/solver_base.hh"
#include "common/field.hh"
#include <functional>
namespace muSpectre {
template <Dim_t DimS, Dim_t DimM=DimS>
class SolverCG: public SolverBase<DimS, DimM>
{
public:
using Parent = SolverBase<DimS, DimM>;
using Ccoord = typename Parent::Ccoord;
using Tg_req_t = typename Parent::TangentRequirement;
// cg only needs to handle fields that look like strain and stress
using Field_t = TensorField<
typename Parent::Collection_t, Real, secondOrder, DimM>;
using Fun_t = std::function<void(const Field_t& in, Field_t & out)>;
constexpr static Tg_req_t tangent_requirement{Tg_req_t::NeedEffect};
//! Default constructor
SolverCG() = delete;
//! Constructor with domain resolutions, etc,
SolverCG(Ccoord resolutions, Real tol, Uint maxiter=0, bool verbose =false);
//! Copy constructor
SolverCG(const SolverCG &other) = delete;
//! Move constructor
- SolverCG(SolverCG &&other) noexcept = default;
+ SolverCG(SolverCG &&other) = default;
//! Destructor
- virtual ~SolverCG() noexcept = default;
+ virtual ~SolverCG() = default;
//! Copy assignment operator
SolverCG& operator=(const SolverCG &other) = delete;
//! Move assignment operator
SolverCG& operator=(SolverCG &&other) = default;
//! actual solver
void solve(const Fun_t & tangent_effect,
const Field_t & rhs,
Field_t & x);
protected:
Tg_req_t get_tangent_req() const override final;
Field_t & r_k; // residual
Field_t & p_k; // search direction
Field_t & Ap_k; // effect of tangent on search direction
const Real tol;
const Uint maxiter;
const bool verbose;
private:
};
} // muSpectre
#endif /* SOLVER_CG_H */
diff --git a/src/system/system_base.hh b/src/system/system_base.hh
index b37a549..4c63ca4 100644
--- a/src/system/system_base.hh
+++ b/src/system/system_base.hh
@@ -1,166 +1,166 @@
/**
* file system_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 Nov 2017
*
* @brief Base class representing a unit cell system with single
* projection operator
*
* @section LICENCE
*
* 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 SYSTEM_BASE_H
#define SYSTEM_BASE_H
#include "common/common.hh"
#include "common/ccoord_operations.hh"
#include "common/field.hh"
#include "common/utilities.hh"
#include "materials/material_base.hh"
#include "fft/projection_base.hh"
#include <vector>
#include <memory>
#include <tuple>
#include <functional>
namespace muSpectre {
//! DimS spatial dimension (dimension of problem
//! DimM material_dimension (dimension of constitutive law)
template <Dim_t DimS, Dim_t DimM>
class SystemBase
{
public:
using Ccoord = Ccoord_t<DimS>;
using Rcoord = Rcoord_t<DimS>;
using FieldCollection_t = FieldCollection<DimS, DimM>;
using Collection_ptr = std::unique_ptr<FieldCollection_t>;
using Material_t = MaterialBase<DimS, DimM>;
using Material_ptr = std::unique_ptr<Material_t>;
using Projection_ptr = std::unique_ptr<ProjectionBase<DimS, DimM>>;
using StrainField_t =
TensorField<FieldCollection_t, Real, secondOrder, DimM>;
using StressField_t =
TensorField<FieldCollection_t, Real, secondOrder, DimM>;
using TangentField_t =
TensorField<FieldCollection_t, Real, fourthOrder, DimM>;
using FullResponse_t =
std::tuple<const StressField_t&, const TangentField_t&>;
using iterator = typename CcoordOps::Pixels<DimS>::iterator;
//! Default constructor
SystemBase() = delete;
//! constructor using sizes and resolution
SystemBase(Projection_ptr projection);
//! Copy constructor
SystemBase(const SystemBase &other) = delete;
//! Move constructor
- SystemBase(SystemBase &&other) noexcept = default;
+ SystemBase(SystemBase &&other) = default;
//! Destructor
- virtual ~SystemBase() noexcept = default;
+ virtual ~SystemBase() = default;
//! Copy assignment operator
SystemBase& operator=(const SystemBase &other) = delete;
//! Move assignment operator
SystemBase& operator=(SystemBase &&other) = default;
/**
* Materials can only be moved. This is to assure exclusive
* ownership of any material by this system
*/
Material_t & add_material(Material_ptr mat);
/**
* evaluates all materials
*/
FullResponse_t evaluate_stress_tangent(StrainField_t & F);
/**
* evaluate directional stiffness (i.e. G:K:δF or G:K:δε)
*/
StressField_t & directional_stiffness(const TangentField_t & K,
const StrainField_t & delF,
StressField_t & delP);
/**
* Convenience function circumventing the neeed to use the
* underlying projection
*/
StressField_t & project(StressField_t & field);
StrainField_t & get_strain();
const StressField_t & get_stress() const;
/**
* general initialisation; initialises the projection and
* fft_engine (i.e. infrastructure) but not the materials. These
* need to be initialised separately
*/
void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
/**
* initialise materials (including resetting any history variables)
*/
void initialise_materials(bool stiffness=false);
iterator begin();
iterator end();
size_t size() const {return pixels.size();}
const Ccoord & get_resolutions() const {return this->resolutions;}
const Rcoord & get_lengths() const {return this->lengths;}
/**
* formulation is hard set by the choice of the projection class
*/
const Formulation & get_formulation() const {
return this->projection->get_formulation();}
protected:
//! make sure that every pixel is assigned to one and only one material
void check_material_coverage();
const Ccoord & resolutions;
CcoordOps::Pixels<DimS> pixels;
const Rcoord & lengths;
Collection_ptr fields;
StrainField_t & F;
StressField_t & P;
//! Tangent field might not even be required; so this is an
//! optional ref_wrapper instead of a ref
optional<std::reference_wrapper<TangentField_t>> K{};
std::vector<Material_ptr> materials{};
Projection_ptr projection;
bool is_initialised{false};
const Formulation form;
private:
};
} // muSpectre
#endif /* SYSTEM_BASE_H */
diff --git a/tests/test_field_collections.cc b/tests/test_field_collections.cc
index 6377ae6..a4547de 100644
--- a/tests/test_field_collections.cc
+++ b/tests/test_field_collections.cc
@@ -1,221 +1,223 @@
/**
* file test_field_collections.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Sep 2017
*
* @brief Test the FieldCollection classes which provide fast optimized iterators
* over run-time typed fields
*
* @section LICENCE
*
* 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_header.hh"
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(field_collection_tests);
BOOST_AUTO_TEST_CASE(simple) {
const Dim_t sdim = 2;
const Dim_t mdim = 2;
using FC_t = FieldCollection<sdim, mdim>;
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<FC_t, Real, order, F::mdim()>;
auto && myfield = make_field<TF_t>("TensorField real 2", F::fc);
using TensorMap = TensorFieldMap<FC_t, Real, order,F::mdim()>;
using MatrixMap = MatrixFieldMap<FC_t, Real, F::mdim(), F::mdim()>;
using ArrayMap = ArrayFieldMap<FC_t, Real, F::mdim(), F::mdim()>;
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<FC_t, T_type, order, F::mdim()>;
using T_TFM2_t = TensorFieldMap<FC_t, T_type, 2, F::mdim()*F::mdim()>; //! dangerous
using T4_Map_t = T4MatrixFieldMap<FC_t, Real, F::mdim()>;
// impossible maptypes for Real tensor fields
using T_SFM_t = ScalarFieldMap<FC_t, T_type>;
using T_MFM_t = MatrixFieldMap<FC_t, T_type, 1, 1>;
using T_AFM_t = ArrayFieldMap<FC_t, T_type, 1, 1>;
using T_MFMw1_t = MatrixFieldMap<FC_t, Int, 1, 2>;
using T_MFMw2_t = MatrixFieldMap<FC_t, Real, 1, 2>;
using T_MFMw3_t = MatrixFieldMap<FC_t, Complex, 1, 2>;
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<FC_t, S_type>;
using S_TFM1_t = TensorFieldMap<FC_t, S_type, 1, 1>;
using S_TFM2_t = TensorFieldMap<FC_t, S_type, 2, 1>;
using S_MFM_t = MatrixFieldMap<FC_t, S_type, 1, 1>;
using S_AFM_t = ArrayFieldMap<FC_t, S_type, 1, 1>;
using S4_Map_t = T4MatrixFieldMap<FC_t, S_type, 1>;
// impossible maptypes for integer scalar fields
using S_MFMw1_t = MatrixFieldMap<FC_t, Int, 1, 2>;
using S_MFMw2_t = MatrixFieldMap<FC_t, Real, 1, 2>;
using S_MFMw3_t = MatrixFieldMap<FC_t, Complex, 1, 2>;
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<FC_t, M_type, F::sdim(), F::mdim()>;
using M_AFM_t = ArrayFieldMap<FC_t, M_type, F::sdim(), F::mdim()>;
// impossible maptypes for complex matrix fields
using M_SFM_t = ScalarFieldMap<FC_t, M_type>;
using M_MFMw1_t = MatrixFieldMap<FC_t, Int, 1, 2>;
using M_MFMw2_t = MatrixFieldMap<FC_t, Real, 1, 2>;
using M_MFMw3_t = MatrixFieldMap<FC_t, Complex, 1, 2>;
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, 2, true>,
FC_multi_fixture<2, 3, true>,
FC_multi_fixture<3, 3, true>>;
using mult_collections_f = boost::mpl::list<FC_multi_fixture<2, 2, false>,
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<F::sdim()> size;
for (auto && s: size) {
s = 3;
}
BOOST_CHECK_NO_THROW(F::fc.initialise(size));
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca, F, mult_collections_f, F) {
testGoodies::RandRange<Int> rng;
for (int i = 0; i < 7; ++i) {
Ccoord_t<F::sdim()> 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) {
+ 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<Int> rng;
- using ftype = internal::TypedFieldBase<decltype(F::fc), Real, mdim*mdim*mdim*mdim>;
+ using ftype = internal::TypedSizedFieldBase<
+ decltype(F::fc), Real, mdim*mdim*mdim*mdim>;
using stype = Eigen::Array<Real, mdim*mdim*mdim*mdim, 1>;
auto & field = reinterpret_cast<ftype&>(F::fc["Tensorfield Real o4"]);
field.push_back(stype());
for (int i = 0; i < nb_pix; ++i) {
Ccoord_t<F::sdim()> 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