diff --git a/CMakeLists.txt b/CMakeLists.txt index ccc7ea3..7b56fb4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,55 +1,55 @@ cmake_minimum_required(VERSION 3.7.0) project(µSpectre_prototype) set(CMAKE_CXX_STANDARD 14) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -add_compile_options(-Wall -Wextra -Werror) +add_compile_options(-Wall -Wextra )#-Werror) enable_testing() set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") # using Clang 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) add_compile_options(-mavx) add_compile_options(-mfma) add_compile_options(-msse) 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) find_package(Eigen3 3.3 REQUIRED NO_MODULE) include_directories( tests src ${EIGEN3_INCLUDE_DIRS} ) #find_package(MPI REQUIRED) find_package(FFTW REQUIRED) #find_package(FFTWMPI REQUIRED) add_subdirectory( "${CMAKE_SOURCE_DIR}/src/" ) #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} Eigen3::Eigen materials common) +target_link_libraries(main_test_suite ${Boost_LIBRARIES} Eigen3::Eigen materials common systems) add_test(main_test_suite main_test_suite) diff --git a/helpers/fft_notes.ipynb b/helpers/fft_notes.ipynb new file mode 100644 index 0000000..8b2b328 --- /dev/null +++ b/helpers/fft_notes.ipynb @@ -0,0 +1,123 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "shape =(5, 3)\n", + "in_arr = np.arange(np.prod(shape)).reshape(shape)\n", + "print(in_arr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "out_arr = np.fft.fftn(in_arr, s=shape)\n", + "print(out_arr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "shape2 = (2,5,3)\n", + "in_arr2 = np.zeros(shape2)\n", + "in_arr2[0,:,:] = in_arr\n", + "in_arr2[1,:,:] = 2.1*in_arr\n", + "print(in_arr2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "out_arr2 = np.fft.fftn(in_arr2, shape)\n", + "print(out_arr2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "out_arr2[0,:,:] == out_arr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "name": "python3" + }, + "name": "fft_notes.ipynb" + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/helpers/material_hyper_elastic.ipynb b/helpers/material_hyper_elastic.ipynb index f684992..4cb108c 100644 --- a/helpers/material_hyper_elastic.ipynb +++ b/helpers/material_hyper_elastic.ipynb @@ -1,318 +1,460 @@ { "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ - "import numpy as np" + "import numpy as np\n", + "import itertools" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "#ndim = 2#3 # number of dimensions\n", "N = 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# tensor operations/products: np.einsum enables index notation, avoiding loops\n", "# e.g. ddot42 performs $C_ij = A_ijkl B_lk$ for the entire grid\n", "trans2 = lambda A2 : np.einsum('ij... ->ji... ',A2 )\n", "ddot42 = lambda A4,B2: np.einsum('ijkl...,lk... ->ij... ',A4,B2)\n", "ddot44 = lambda A4,B4: np.einsum('ijkl...,lkmn...->ijmn...',A4,B4)\n", "dot22 = lambda A2,B2: np.einsum('ij... ,jk... ->ik... ',A2,B2)\n", "dot24 = lambda A2,B4: np.einsum('ij... ,jkmn...->ikmn...',A2,B4)\n", "dot42 = lambda A4,B2: np.einsum('ijkl...,lm... ->ijkm...',A4,B2)\n", "dyad22 = lambda A2,B2: np.einsum('ij... ,kl... ->ijkl...',A2,B2)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# material parameters + function to convert to grid of scalars\n", "phase = .1\n", "param = lambda M0,M1: M0*np.ones([N,N,N])*(1.-phase)+M1*np.ones([N,N,N])*phase\n", "Kv = param(0.833,8.33) # bulk modulus [grid of scalars]\n", "mu = param(0.386,3.86) # shear modulus [grid of scalars]\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# bulk modulus K = lambda - 2*µ/3\n", "# shear modulus µ = lamé's µ\n", "Young = 9*Kv*mu/(3*Kv + mu)\n", "Poisson = (3*Kv-2*mu)/(2*(3*Kv+ mu))\n", "print(\"Young = {}\".format(Young.flatten()[0]))\n", "print(\"Poisson = {}\".format(Poisson.flatten()[0]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def mat(dim):\n", " if dim == 2:\n", " return np.array([[0, 0],\n", " [1, 1],\n", " [0, 1],\n", " [1, 0]])\n", " elif dim == 3:\n", " return np.array([[0, 0],\n", " [1, 1],\n", " [2, 2],\n", " [1, 2],\n", " [0, 2],\n", " [0, 1],\n", " [2, 1],\n", " [2, 0],\n", " [1, 0]])\n", " else:\n", " raise Exception(\"unsupported dimension\")\n", "\n", "def vsize(dim):\n", " return dim * (dim - 1) // 2 + dim\n", "def V4(arr, sym=True):\n", "\n", " dim = arr.shape[0]\n", " if sym:\n", " siz = vsize(dim)\n", " else:\n", " siz = arr.shape[0]**2\n", "\n", " vout = np.zeros([siz, siz])\n", " ids = mat(dim)\n", " for I in range(siz):\n", " i, j = ids[I]\n", " for J in range(siz):\n", " k, l = ids[J]\n", " vout[I, J] = arr [i, j, k, l]\n", " pass\n", " pass\n", " return vout\n", "\n", "def print_eigen_input(arr):\n", " dim = len(arr.shape)\n", " assert dim in {1, 2}, \"only 1d and 2d arrays\"\n", "\n", " if dim == 1:\n", " return \"<<\\n\" + \", \".join((str(i) for i in arr.tolist()))\n", " else:\n", " return \"<<\\n\" + \",\\n\".join(\n", " \", \".join((str(i) for i in row))\n", " for row in arr.tolist())\n", " pass\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def constitutive(F):\n", " ndim = F.shape[0]\n", " # identity tensor [single tensor]\n", " i = np.eye(ndim)\n", " # identity tensors [grid of tensors]\n", " I = i\n", "\n", " I4 = np.einsum('il,jk',i,i)\n", "\n", " I4rt = np.einsum('ik,jl',i,i)\n", "\n", " I4s = (I4+I4rt)/2.\n", "\n", " II = dyad22(I,I)\n", "\n", " C4 = Kv*II+2.*mu*(I4s-1./3.*II)\n", " S = ddot42(C4,.5*(dot22(trans2(F),F)-I))\n", " P = dot22(F,S)\n", " K4 = dot24(S,I4)+ddot44(ddot44(I4rt,dot42(dot24(F,C4),trans2(F))),I4rt)\n", " return P,K4,C4,S\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "\n", + "Fs = [np.array([[1, .1],\n", + " [.2, 1]]),\n", + " np.array([[1, .1, 0],\n", + " [.2, 1, 0],\n", + " [0, 0, 1]]),\n", + " np.array([[1, .1, .3],\n", + " [.2, 1, .4],\n", + " [.5, .6, 1]])]\n", + "\n", + "for F in Fs:\n", + " P, K, C, S = constitutive(F)\n", + " print(\"F =\\n{}\".format(print_eigen_input(F )))\n", + " print(\"P =\\n{}\".format(print_eigen_input(P )))\n", + " print(\"K =\\n{}\".format(print_eigen_input(V4(K, sym=False))))\n", + "\n" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def illustrate(F):\n", " ndim = F.shape[0]\n", " # identity tensor [single tensor]\n", " i = np.eye(ndim)\n", " # identity tensors [grid of tensors]\n", " I = i\n", "\n", " I4 = np.einsum('il,jk',i,i)\n", "\n", " I4rt = np.einsum('ik,jl',i,i)\n", "\n", " I4s = (I4+I4rt)/2.\n", "\n", " II = dyad22(I,I)\n", "\n", " C4 = Kv*II+2.*mu*(I4s-1./3.*II)\n", " S = ddot42(C4,.5*(dot22(trans2(F),F)-I))\n", " P = dot22(F,S)\n", " def printa(name, var):\n", " print(\"{} =\\n{}\".format(name, V4(var, sym=False)))\n", " C_FT = dot42(C4, trans2(F))\n", " printa('C_FT', C_FT)\n", " F_C_FT = dot24(F, C_FT)\n", " printa('F_C_FT', F_C_FT)\n", " IRT_F_C_FT = ddot44(I4rt, F_C_FT)\n", " printa('IRT_F_C_FT', IRT_F_C_FT)\n", " IRT_F_C_FT_IRT = ddot44(IRT_F_C_FT, I4rt)\n", " printa('IRT_F_C_FT_IRT', IRT_F_C_FT_IRT)\n", " S_I = dot24(S,I4)\n", " printa('S_I', S_I)\n", " my_K4 = S_I + IRT_F_C_FT_IRT\n", " K4 = dot24(S,I4)+ddot44(ddot44(I4rt,dot42(dot24(F,C4),trans2(F))),I4rt)\n", " print(\"max Error = {}\".format(abs(my_K4 - K4).max()))\n", "\n", "illustrate(Fs[0])" ] }, + { + "cell_type": "markdown", + "metadata": { + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "source": [ + "Computation of Projection operator" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ + "def Ghat(ndim, N):\n", + " delta = lambda i,j: np.float(i==j) # Dirac delta function\n", + " freq = np.arange(-(N-1)/2.,+(N+1)/2.) # coordinate axis -> freq. axis\n", + " Ns = [N for _ in range(ndim)]\n", + " Ghat4 = np.zeros([ndim,ndim,ndim,ndim]+Ns) # zero initialize\n", + " # - compute\n", + " for i,j,l,m in itertools.product(range(ndim),repeat=4):\n", + " for xyz in itertools.product(range(N), repeat=ndim):\n", + " q = np.array([freq[coord] for coord in xyz]) # frequency vector\n", + " # zero freq. -> mean\n", + " if not q.dot(q) == 0:\n", + " if len(xyz) == 2:\n", + " Ghat4[i,j,l,m, xyz[0], xyz[1]] = delta(i,m)*q[j]*q[l]/(q.dot(q))\n", + " elif len(xyz) == 3:\n", + " Ghat4[i,j,l,m, xyz[0], xyz[1], xyz[2]] = delta(i,m)*q[j]*q[l]/(q.dot(q))\n", + " else:\n", + " print(xyz)\n", + " raise Exception(\"'{}' is an illegal dim\".format(len(xyz)))\n", + " return Ghat4\n", "\n", - "Fs = [np.array([[1, .1],\n", - " [.2, 1]]),\n", - " np.array([[1, .1, 0],\n", - " [.2, 1, 0],\n", - " [0, 0, 1]]),\n", - " np.array([[1, .1, .3],\n", - " [.2, 1, .4],\n", - " [.5, .6, 1]])]\n", - "\n", - "for F in Fs:\n", - " P, K, C, S = constitutive(F)\n", - " print(\"F =\\n{}\".format(print_eigen_input(F )))\n", - " print(\"P =\\n{}\".format(print_eigen_input(P )))\n", - " print(\"K =\\n{}\".format(print_eigen_input(V4(K, sym=False))))\n" + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "G = Ghat(2, 55).flatten()\n", + "print(\"empty-ratio: {}\".format((G==0).sum()/G.size))\n", + "print(\"size: {}:\".format(G.size))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "for n in (8, 9):\n", + " print(print_eigen_input(np.fft.fftfreq(n, 15.2/n)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "np.fft.fftfreq??" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "import scipy as sc\n", + "import scipy.sparse.linalg as sp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "autoscroll": "json-false", + "ein.tags": [ + "worksheet-0" + ], + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "sp.cg??" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "name": "python3" }, "name": "material_hyper_elastic.ipynb" }, "nbformat": 4, "nbformat_minor": 2 } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3baf61f..851311e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,8 +1,10 @@ include_directories( common materials + system ) add_subdirectory(common) add_subdirectory(materials) +add_subdirectory(system) diff --git a/src/materials/material.hh b/src/materials/material.hh index d382114..731bfe4 100644 --- a/src/materials/material.hh +++ b/src/materials/material.hh @@ -1,232 +1,232 @@ /** * file material.hh * * @author Till Junge * * @date 01 May 2017 * * @brief Base class for materials (constitutive models) * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include #include "common/common.hh" #include "common/voigt_conversion.hh" #ifndef MATERIAL_H #define MATERIAL_H namespace muSpectre { - //! s_dim spatial dimension (dimension of problem - //! m_dim material_dimension (dimension of constitutive law) + //! DimS spatial dimension (dimension of problem + //! DimM material_dimension (dimension of constitutive law) template class MaterialBase { public: //! Gradients are second-order tensors per pixel using SecondArray = Eigen::Tensor; using SecondArrayMap = Eigen::TensorMap; //! Tangent moduli are fourth-order tensors per pixel using FourthArray = Eigen::Tensor; using FourthArrayMap = Eigen::TensorMap; // commonly used representations od data; using Index = Index_t; using VoigtConv = VoigtConversion; const static size_t nb_voigt_sym = vsize(DimM); const static size_t nb_voigt_asy = vsize(DimM); using SecondTens = Eigen::Matrix; using SecondMap = Eigen::Map; using SecondSymVoigt = Eigen::Map>; using SecondVoigt = Eigen::Map>; using SecondTuple = std::tuple; using FourthTens = Eigen::TensorFixedSize>; using FourthMap = Eigen::TensorMap; using FourthSymVoigt = Eigen::Map>; using FourthVoigt = Eigen::Map>; using FourthTuple = std::tuple; //! Default constructor MaterialBase(); //! take responsibility for a pixel void add_pixel(const Index && local_index); //! allocate memory, etc virtual void initialize() = 0; virtual void compute_First_Piola_Kirchhoff_stress(const SecondArrayMap & F, SecondArrayMap & P, FourthArrayMap & K) = 0; // Anything that has a data() method that yields a Real* to the first elem template void compute_First_Piola_Kirchhoff_stress(const SecondArrF & F, SecondArrP & P, FourthArrK & K); //! convenience function returning a suitable map view to a raw array inline static SecondArrayMap get_second_array_map(Real * data, Index sizes); inline static FourthArrayMap get_fourth_array_map(Real * data, Index sizes); inline static const SecondArrayMap get_second_array_map(const Real * data, Index sizes) { return get_second_array_map(const_cast(data), sizes); } inline static const FourthArrayMap get_fourth_array_map(const Real * data, Index sizes) { return get_fourth_array_map(const_cast(data), sizes); } protected: inline SecondMap get_view(size_t id, SecondArrayMap & arr); inline FourthMap get_view(size_t id, FourthArrayMap & arr); inline const SecondMap get_const_view(size_t id, const SecondArrayMap & arr); inline const FourthMap get_const_view(size_t id, const FourthArrayMap & arr); //! container of pixel indices using this material std::vector indices; private: }; namespace internal { //----------------------------------------------------------------------------// template inline auto * getter (tens_t & tens, arr_t1 & leading_zeros, arr_t2 & index, const std::index_sequence&, const std::index_sequence&) { return &tens(leading_zeros[J]..., index[I]...); } //----------------------------------------------------------------------------// template inline auto * get(Eigen::TensorMap> & tens, Index_t & index) { const Dim_t diff = order - dim; Index_t leading_zeros{0}; return getter(tens, leading_zeros, index, std::make_index_sequence{}, std::make_index_sequence{}); } //----------------------------------------------------------------------------// template inline auto * const_getter (const tens_t & tens, arr_t1 & leading_zeros, arr_t2 & index, const std::index_sequence&, const std::index_sequence&) { return &tens(leading_zeros[J]..., index[I]...); } //----------------------------------------------------------------------------// template inline auto * const_get(const Eigen::TensorMap> & tens, Index_t & index) { const Dim_t diff = order - dim; Index_t leading_zeros{0}; return const_getter(tens, leading_zeros, index, std::make_index_sequence{}, std::make_index_sequence{}); } //----------------------------------------------------------------------------// template inline auto get_map(Real * data, Index_t& leading_dims, Index_t & sizes, const std::index_sequence&, const std::index_sequence&) { return Eigen::TensorMap>(data, leading_dims[I]..., sizes[J]...); } } // internal template inline typename MaterialBase::SecondArrayMap MaterialBase::get_second_array_map(Real * data, Index sizes) { const size_t order{2}; Index_t leading_dims; leading_dims.fill(DimM); return internal::get_map(data, leading_dims, sizes, std::make_index_sequence{}, std::make_index_sequence{}); } template inline typename MaterialBase::FourthArrayMap MaterialBase::get_fourth_array_map(Real * data, Index sizes) { const size_t order{4}; Index_t leading_dims; leading_dims.fill(DimM); return internal::get_map(data, leading_dims, sizes, std::make_index_sequence{}, std::make_index_sequence{}); } //----------------------------------------------------------------------------// template inline typename MaterialBase::SecondMap MaterialBase::get_view(size_t id, SecondArrayMap & arr) { auto * data_ptr = internal::get(arr, this->indices[id]); return SecondMap(data_ptr); } //----------------------------------------------------------------------------// template inline const typename MaterialBase::SecondMap MaterialBase::get_const_view(size_t id, const SecondArrayMap & arr) { auto * data_ptr = internal::const_get(arr, this->indices[id]); return SecondMap(const_cast(data_ptr)); } //----------------------------------------------------------------------------// template inline typename MaterialBase::FourthMap MaterialBase::get_view(size_t id, FourthArrayMap & arr) { auto * data_ptr = internal::get(arr, this->indices[id]); return FourthMap(data_ptr, DimM, DimM, DimM, DimM); } //----------------------------------------------------------------------------// template inline const typename MaterialBase::FourthMap MaterialBase::get_const_view(size_t id, const FourthArrayMap & arr) { auto * data_ptr = internal::const_get(arr, this->indices[id]); return FourthMap(const_cast(data_ptr), DimM, DimM, DimM, DimM); } } // muSpectre #endif /* MATERIAL_H */ diff --git a/src/system/CMakeLists.txt b/src/system/CMakeLists.txt new file mode 100644 index 0000000..838b37b --- /dev/null +++ b/src/system/CMakeLists.txt @@ -0,0 +1,8 @@ +set (systems_SRC + system_base.cc + fftw_engine_c2c.cc + ) +find_package(FFTW REQUIRED) + +add_library(systems ${systems_SRC}) +target_link_libraries(systems Eigen3::Eigen common ${FFTW_LIBRARIES}) diff --git a/src/system/fft_engine.hh b/src/system/fft_engine.hh new file mode 100644 index 0000000..4b09f0e --- /dev/null +++ b/src/system/fft_engine.hh @@ -0,0 +1,63 @@ +/** + * file fft_engine.hh + * + * @author Till Junge + * + * @date 11 May 2017 + * + * @brief abstract base class defining interface for fft_systems + * + * @section LICENCE + * + * Copyright (C) 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 "common/common.hh" +#include +#include + + +#ifndef FFT_ENGINE_H +#define FFT_ENGINE_H + +namespace muSpectre { + enum class FFT_PlanFlags {estimate, measure, patient}; + + template + class FFT_Engine { + public: + using Index = Index_t; + FFT_Engine(Index nb_pixels) + :nb_pixels(nb_pixels), + tot_nb_pixels(std::accumulate(nb_pixels.begin(), nb_pixels.end(), + 1, std::multiplies())){}; + + virtual void init(FFT_PlanFlags flags) = 0; + + virtual Real * convolve(const Real * const Ghat, const Real * const arr) = 0; + + protected: + const Index nb_pixels; + size_t tot_nb_pixels; + + }; + +} // muSpectre + + +#endif /* FFT_ENGINE_H */ diff --git a/src/system/fftw_engine_c2c.cc b/src/system/fftw_engine_c2c.cc new file mode 100644 index 0000000..ca30b8c --- /dev/null +++ b/src/system/fftw_engine_c2c.cc @@ -0,0 +1,139 @@ +/** + * file fftw_engine_c2c.cc + * + * @author Till Junge + * + * @date 11 May 2017 + * + * @brief fft_engine usinge FFTW + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +#include +#include +#include "system/fftw_engine_c2c.hh" +#include + +namespace muSpectre { + + template + FFTW_EngineC2C::FFTW_EngineC2C(Index nb_pixels) + :FFT_Engine(nb_pixels) + { + this->i_work_space = + fftw_alloc_complex(this->tot_nb_pixels* DimM*DimM); + this->o_work_space = + fftw_alloc_complex(this->tot_nb_pixels* DimM*DimM); + } + + //----------------------------------------------------------------------------// + template + FFTW_EngineC2C::~FFTW_EngineC2C() + { + fftw_destroy_plan(this->plan_fft); + fftw_destroy_plan(this->plan_ifft); + fftw_free(this->i_work_space); + fftw_free(this->o_work_space); + fftw_cleanup(); + } + + //----------------------------------------------------------------------------// + template + void FFTW_EngineC2C::init(FFT_PlanFlags inflags) { + const int & rank = DimS; + std::array narr; + const int * const n = &narr[0]; + std::copy(this->nb_pixels.begin(), this->nb_pixels.end(), narr.begin()); + int howmany = DimM*DimM; + fftw_complex * in = reinterpret_cast(this->i_work_space); + const int * const inembed = n; + int istride = howmany; + int idist = 1; + fftw_complex * out = in; + const int * const onembed = n; + int ostride = howmany; + int odist = idist; + int sign = FFTW_FORWARD; + unsigned int flags; + if (inflags == FFT_PlanFlags::estimate) { + flags = FFTW_ESTIMATE; + } + if (inflags == FFT_PlanFlags::measure) { + flags = FFTW_MEASURE; + } + if (inflags == FFT_PlanFlags::patient) { + flags = FFTW_PATIENT; + } + this->plan_fft = fftw_plan_many_dft(rank, n, howmany, in, inembed, istride, + idist, out, onembed, ostride, odist, sign, + flags); + + sign = FFTW_BACKWARD; + in = reinterpret_cast(this->o_work_space); + + this->plan_ifft = fftw_plan_many_dft(rank, n, howmany, in, inembed, istride, + idist, out, onembed, ostride, odist, sign, + flags); + } + + //----------------------------------------------------------------------------// + template + Real * FFTW_EngineC2C::convolve(const Real * const Ghat, + const Real * const arr) { + + std::copy(arr, arr+this->tot_nb_pixels, + reinterpret_cast*>(this->i_work_space)); + + fftw_execute(this->plan_fft); + + //----------------------------------------------------------------------------// + const std::array,2> prod_dims {Eigen::IndexPair{2, 1}, + Eigen::IndexPair{3, 0}}; + for (size_t i = 0; i < this->tot_nb_pixels; ++i) { + + using T4shape = Eigen::Sizes; + using T2shape = Eigen::Sizes; + using T4type = Eigen::TensorFixedSize; + using T2type = Eigen::TensorFixedSize, T2shape>; + using T4map = Eigen::TensorMap; + using T2map = Eigen::TensorMap; + const Real* const G_ptr = &Ghat[i*DimM*DimM*DimM*DimM]; + const T4map G(const_cast(G_ptr),DimM, DimM, DimM,DimM); + std::complex * P_ptr = + &reinterpret_cast*>(&this->i_work_space)[i*DimM*DimM]; + const T2map P(P_ptr, DimM, DimM); + + std::complex* Pout_ptr = + &reinterpret_cast*>(&this->o_work_space)[i*DimM*DimM]; + T2map Pout(Pout_ptr, DimM, DimM); + auto Gc = G.template cast< std::complex >(); + Pout = Gc.contract(P, prod_dims); + } + //----------------------------------------------------------------------------// + fftw_execute(this->plan_ifft); + + return reinterpret_cast(this->o_work_space); + } + + template class FFTW_EngineC2C<1, 1>; + template class FFTW_EngineC2C<2, 2>; + template class FFTW_EngineC2C<3, 3>; +} // muSpectre diff --git a/src/system/fftw_engine_c2c.hh b/src/system/fftw_engine_c2c.hh new file mode 100644 index 0000000..a255514 --- /dev/null +++ b/src/system/fftw_engine_c2c.hh @@ -0,0 +1,58 @@ +/** + * file fftw_engine_c2c.hh + * + * @author Till Junge + * + * @date 11 May 2017 + * + * @brief implementaion of fft_engine using redundant c2c fftw transform + * + * @section LICENCE + * + * Copyright (C) 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 "system/fft_engine.hh" +#include + +#ifndef FFTW_ENGINE_C2C_H +#define FFTW_ENGINE_C2C_H + + +namespace muSpectre { + + template + class FFTW_EngineC2C: public FFT_Engine { + public: + using Index = Index_t; + FFTW_EngineC2C(Index nb_pixels); + virtual ~FFTW_EngineC2C(); + virtual void init(FFT_PlanFlags flags) override; + virtual Real * convolve(const Real * const Ghat, const Real * const arr) override; + + protected: + void * i_work_space, * o_work_space; + fftw_plan plan_fft; + fftw_plan plan_ifft; + + }; + +} // muSpectre + + +#endif /* FFTW_ENGINE_C2C_H */ diff --git a/src/system/system_base.cc b/src/system/system_base.cc new file mode 100644 index 0000000..e6427e1 --- /dev/null +++ b/src/system/system_base.cc @@ -0,0 +1,186 @@ +/** + * file system_base.cc + * + * @author Till Junge + * + * @date 09 May 2017 + * + * @brief Implementation for system base class + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#include +#include +#include +#include "system/system_base.hh" +#include "iostream" + + +namespace muSpectre { + + namespace Projection { + + //----------------------------------------------------------------------------// + //! computes fft frequecies. see function fftfreq in numpy's helper.py + vector fft_freqs(Uint nb_samples, Real length) { + vector retval(1, nb_samples); + int N = (nb_samples-1)/2 + 1; // needs to be signed int for neg freqs + for (int i = 0; i < N; ++i) { + retval(0, i) = i; + } + for (int i = N; i < static_cast(nb_samples); ++i) { + retval(0, i) = -(static_cast(nb_samples)/2)+i-N; + + } + + return retval/length; + } + + //----------------------------------------------------------------------------// + template + FreqStruc::FreqStruc(const std::array& sizes, + const Index_t& nb_pixels) + :nb_pixels(nb_pixels), sizes(sizes) + { + for (size_t i = 0; i < Dim; ++i) { + this->xi[i] = Projection::fft_freqs(this->nb_pixels[i], this->sizes[i]); + } + } + + //----------------------------------------------------------------------------// + template + inline typename FreqStruc::iterator FreqStruc::begin() { + return iterator(*this, {0}); + } + + //----------------------------------------------------------------------------// + template + inline typename FreqStruc::iterator FreqStruc::end() { + return iterator(*this, this->nb_pixels); + } + + //----------------------------------------------------------------------------// + template + FreqStruc::iterator::iterator(const FreqStruc & parent, + Index_t index) + :parent(parent), + index(std::accumulate(index.begin(), index.end(), 1, std::multiplies())){ } + + //----------------------------------------------------------------------------// + template + inline Index_t FreqStruc::iterator::get_index() { + // implements row-major choice. might have to be templated at some point + static_assert(((Dim == twoD) || (Dim == threeD)), "only 2d or 3d"); + auto & stride = this->parent.nb_pixels[1]; + return Index_t{static_cast(this->index/stride), + static_cast(this->index%stride)}; + } + //----------------------------------------------------------------------------// + template<> + inline Index_t FreqStruc::iterator::get_index() { + // implements row-major choice. might have to be templated at some point + auto && kstride = this->parent.nb_pixels[2]; + Dim_t && k = this->index%kstride; + auto && jstride = this->parent.nb_pixels[1]; + Dim_t && j = (this->index/kstride)%jstride; + Dim_t && i = (this->index/kstride)/jstride; + return Index_t{i, j, k}; + } + + //----------------------------------------------------------------------------// + template + inline Eigen::Matrix FreqStruc::iterator::operator*() { + Eigen::Matrix retval; + auto && indices = this->get_index(); + if (Dim == twoD) { + retval << parent.xi[0](indices[0]), parent.xi[1](indices[1]); + } else if (Dim == threeD) { + retval << + parent.xi[0](indices[0]), + parent.xi[1](indices[1]), + parent.xi[2](indices[2]); + } + return retval; + } + + //----------------------------------------------------------------------------// + template + inline typename FreqStruc::iterator& FreqStruc::iterator::operator++() { + this->index++; + return *this; + } + + } // Projection + + //----------------------------------------------------------------------------// + template + SystemBase::SystemBase(std::array sizes, Index nb_pixels, + Formulation form) + : sizes{sizes}, nb_pixels{nb_pixels}, + nb_pixel(std::accumulate(nb_pixels.begin(), nb_pixels.end(), 1, std::multiplies())), + form{form} + { + for (const auto & size: this->nb_pixels) { + if (size%2 != 1) { + throw std::runtime_error("Sorry, for now can only handle odd grid sizes"); + } + } + + if (form != Formulation::finite_strain) { + throw std::runtime_error("Sorry, only finite strain for the time being"); + } + //resize and allocate + this->Ghats = T4Vector(DimM, DimM, DimM, DimM, this->size()); + this->Ghats.setZero(); + this->build_ghats(); + } + + //----------------------------------------------------------------------------// + template + void SystemBase::build_ghats() { + //! row-major storage for pixels + Projection::FreqStruc freqs(this->sizes, this->nb_pixels); + auto it = freqs.begin(); + //drop the first frequency, as it is zero + ++it; + for (size_t pix_id = 1; pix_id < this->nb_pixel; ++pix_id, ++it) { + auto xi = *it; + xi/=xi.norm(); + using T4shape = Eigen::Sizes; + using T4 = Eigen::TensorFixedSize; + Eigen::TensorMap Ghat_kk(&this->Ghats(0, 0, 0, 0, pix_id), + DimS, DimS, DimS, DimS); + for (Dim_t im = 0; im < DimS; ++im) { + for (Dim_t j = 0; j < DimS; ++j) { + for (Dim_t l = 0; l < DimS; ++l) { + Ghat_kk(im, j, l, im) = xi(j)*xi(l); + } + } + } + } + } + + template class SystemBase<2, 2>; + template class SystemBase<2, 3>; + template class SystemBase<3, 3>; + + +} // muSpectre diff --git a/src/system/system_base.hh b/src/system/system_base.hh new file mode 100644 index 0000000..6b2c1c9 --- /dev/null +++ b/src/system/system_base.hh @@ -0,0 +1,119 @@ +/** + * file system_base.hh + * + * @author Till Junge + * + * @date 09 May 2017 + * + * @brief Base class representing a periodic cell and it's solution + * + * @section LICENCE + * + * Copyright (C) 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 "common/common.hh" +#include +#include + + +#ifndef SYSTEM_BASE_H +#define SYSTEM_BASE_H + + +namespace muSpectre { + namespace Projection { + + using vector = Eigen::Matrix; + + + //----------------------------------------------------------------------------// + //! computes fft frequecies. see function fftfreq in numpy's helper.py + vector fft_freqs(Uint nb_samples, Real length); + + template + class FreqStruc { + public: + FreqStruc(const std::array& sizes, const Index_t& nb_pixels); + + class iterator { + public: + iterator(const FreqStruc & parent, Index_t index); + inline Eigen::Matrix operator*(); + inline iterator& operator++ (); + inline bool operator!=(const iterator & other) const + {return this->index !=other.index;} + protected: + inline Index_t get_index(); + const FreqStruc & parent; + size_t index; + }; + + inline iterator begin(); + inline iterator end(); + + protected: + const Index_t& nb_pixels; + const std::array& sizes; + //! scaled frequencies (k_i/L_i) + std::array xi; + }; + + } // Projection + + //! DimS spatial dimension (dimension of problem + //! DimM material_dimension (dimension of constitutive law) + template + class SystemBase + { + public: + using Index = Index_t; + using T4Vector = Eigen::Tensor; + using T2Vector = Eigen::Tensor; + using T4VecMap = Eigen::TensorMap; + using T2VecMap = Eigen::TensorMap; + + enum class Formulation{finite_strain, small_strain}; + + SystemBase(std::array sizes, Index nb_pixels, + Formulation = Formulation::finite_strain); + virtual ~SystemBase() = default; + + inline size_t size() {return this->nb_pixel;} + protected: + void build_ghats(); + T4Vector get_scaled_freqs; + + std::array sizes; + Index nb_pixels; + size_t nb_pixel; + Formulation form; + //! projection operator + T4Vector Ghats; + //! tangent stiffness tensors by pixel; + //T4VecMap Ks; + //! gradient or strain (depending on finite/small strain formulation) + //T2VecMap grad; + //! Cauchy or First Piola-Kirchhoff stress + //T2VecMap stress; + }; + + +} // muSpectre + +#endif /* SYSTEM_BASE_H */ diff --git a/tests/test_fftw_c2c.cc b/tests/test_fftw_c2c.cc new file mode 100644 index 0000000..b60267c --- /dev/null +++ b/tests/test_fftw_c2c.cc @@ -0,0 +1,68 @@ +/** + * file test_fftw_c2c.cc + * + * @author Till Junge + * + * @date 11 May 2017 + * + * @brief test fftw redundant class + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#include +#include +#include "complex" + +#include "system/fftw_engine_c2c.hh" +#include "tests.hh" + +namespace muSpectre { + template + struct fftw_c2c_fixture: public FFTW_EngineC2C { + fftw_c2c_fixture() + :FFTW_EngineC2C(std_nb_pix){} + + const static typename FFTW_EngineC2C::Index std_nb_pix; + }; + + template<> const Index_t<1> fftw_c2c_fixture<1,1>::std_nb_pix{3}; + template<> const Index_t<2> fftw_c2c_fixture<2,2>::std_nb_pix{5,3}; + template<> const Index_t<3> fftw_c2c_fixture<3,3>::std_nb_pix{7,5,3}; + + using fix2 = fftw_c2c_fixture<2, 2>; + BOOST_FIXTURE_TEST_CASE(numpy_comparison, fix2) { + const size_t dim = 2; + init(FFT_PlanFlags::estimate); + Real in_arr[tot_nb_pixels*dim*dim]; + for (size_t i = 0; i < tot_nb_pixels; ++i) { + for (size_t d = 0; d < dim*dim; ++d) { + in_arr[i*dim*dim+d] = i; + } + + in_arr[i] = i; + } + + Real * out_arr = convolve(nullptr, in_arr); + std::complex* cptr = reinterpret_cast*>(out_arr); + + } + +} // muSpectre diff --git a/tests/test_material_hyperelastic.cc b/tests/test_material_hyperelastic.cc index b22bb91..0453096 100644 --- a/tests/test_material_hyperelastic.cc +++ b/tests/test_material_hyperelastic.cc @@ -1,216 +1,216 @@ /** * file test_material_hyperelastic.cc * * @author Till Junge * * @date 01 May 2017 * * @brief Test the hyper-elastic material used in Section 4 of * de Geus et al. / Comput. Methods in Appl. Mech. Engrg. * 318, (2017), 412–430 https://doi.org/10.1016/j.cma.2016.12.032 * * @section LICENCE * * Copyright (C) 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 "materials/material_hyper_elastic.hh" #include #include #include +#include "tests.hh" namespace muSpectre { - const Real tol = 1e-14*100; //it's in percent // params from "material_hyper_elastic.ipynb"; const Real E = 1.9058231542461002; const Real nu = 0.2993067590987868; template struct mat_fixture: public MaterialHyperElastic{ mat_fixture(): MaterialHyperElastic(E, nu){ }; const static Dim_t DimS = s_dim; const static Dim_t DimM = m_dim; /** inputs with expected outputs computed with * "material_hyper_elastic.ipynb" in folder "helpers"*/ using parent = MaterialHyperElastic; using SecondTens = typename parent::SecondTens; using FourthVoigt = Eigen::Matrix; const static SecondTens F, P; const static FourthVoigt K; }; //----------------------------------------------------------------------------// template <> const mat_fixture<2, 2>::SecondTens mat_fixture<2, 2>::F = (mat_fixture<2, 2>::SecondTens() << 1.0, 0.1, 0.2, 1.0).finished(); template <> const mat_fixture<2, 2>::SecondTens mat_fixture<2, 2>::P = (mat_fixture<2, 2>::SecondTens()<< 0.07868216666666673, 0.22348781666666673, 0.2313560333333334, 0.07868216666666672).finished(); template <> const mat_fixture<2, 2>::FourthVoigt mat_fixture<2, 2>::K = (mat_fixture<2, 2>::FourthVoigt()<< 2.624580833333334, 1.1084346666666667, 0.40273666666666674, 0.5854533333333334, 1.1084346666666667, 2.6245808333333334, 0.40273666666666674, 0.5854533333333334, 0.5854533333333334, 0.5854533333333334, 0.7552753333333334, 0.8925028333333335, 0.40273666666666674, 0.40273666666666674, 0.7936838333333334, 0.7552753333333334).finished(); //----------------------------------------------------------------------------// template <> const mat_fixture<2, 3>::SecondTens mat_fixture<2, 3>::F = (mat_fixture<2, 3>::SecondTens()<< 1.0, 0.1, 0.0, 0.2, 1.0, 0.0, 0.0, 0.0, 1.0).finished(); template <> const mat_fixture<2, 3>::SecondTens mat_fixture<2, 3>::P = (mat_fixture<2, 3>::SecondTens()<< 0.07868216666666673, 0.22348781666666673, 0.0, 0.2313560333333334, 0.07868216666666672, 0.0, 0.0, 0.0, 0.027344166666666694).finished(); template <> const mat_fixture<2, 3>::FourthVoigt mat_fixture<2, 3>::K = (mat_fixture<2, 3>::FourthVoigt()<<2.624580833333334, 1.1084346666666667, 1.0937666666666668, 0.0, 0.0, 0.40273666666666674, 0.0, 0.0, 0.5854533333333334, 1.1084346666666667, 2.6245808333333334, 1.0937666666666668, 0.0, 0.0, 0.40273666666666674, 0.0, 0.0, 0.5854533333333334, 1.0937666666666668, 1.0937666666666668, 2.5879108333333334, 0.0, 0.0, 0.10937666666666668, 0.0, 0.0, 0.21875333333333336, 0.0, 0.0, 0.0, 0.7334, 0.07334, 0.0, 0.7680781666666667, 0.22002000000000005, 0.0, 0.0, 0.0, 0.0, 0.14668, 0.7334, 0.0, 0.22002000000000005, 0.7900801666666668, 0.0, 0.5854533333333334, 0.5854533333333334, 0.21875333333333336, 0.0, 0.0, 0.7552753333333334, 0.0, 0.0, 0.8925028333333335, 0.0, 0.0, 0.0, 0.7900801666666668, 0.22002, 0.0, 0.7334, 0.14668, 0.0, 0.0, 0.0, 0.0, 0.22002, 0.7680781666666667, 0.0, 0.07334, 0.7334, 0.0, 0.40273666666666674, 0.40273666666666674, 0.10937666666666668, 0.0, 0.0, 0.7936838333333334, 0.0, 0.0, 0.7552753333333334).finished(); //----------------------------------------------------------------------------// template <> const mat_fixture<3, 3>::SecondTens mat_fixture<3, 3>::F = (mat_fixture<3, 3>::SecondTens()<< 1.0, 0.1, 0.3, 0.2, 1.0, 0.4, 0.5, 0.6, 1.0 ).finished(); template <> const mat_fixture<3, 3>::SecondTens mat_fixture<3, 3>::P = (mat_fixture<3, 3>::SecondTens()<< 0.9479714333333337, 0.7435627833333334, 0.92523635, 0.8402667666666668, 1.1591906333333337, 1.1568859333333334, 1.264590916666667, 1.4368351000000001, 1.4569510333333335).finished(); template <> const mat_fixture<3, 3>::FourthVoigt mat_fixture<3, 3>::K = (mat_fixture<3,3>::FourthVoigt() << 3.3442565, 1.1084346666666667, 1.2037766666666667, 0.4815106666666667, 1.193542, 0.6227566666666668, 0.69293, 1.5443073333333335, 0.6734613333333335, 1.1084346666666667, 3.4762685000000006, 1.2697826666666667, 1.4862686666666667, 0.35746600000000006, 0.4907446666666667, 1.90304, 0.6348913333333334, 0.8054733333333335, 1.2037766666666667, 1.2697826666666667, 3.6889545000000004, 1.5376066666666668, 1.178874, 0.2413886666666667, 1.851702, 1.5589753333333336, 0.3654333333333334, 0.69293, 1.90304, 1.851702, 0.9959040000000001, 0.270218, 0.7403540000000001, 2.6075758333333336, 0.9881900000000001, 0.49795200000000006, 1.5443073333333335, 0.6348913333333334, 1.5589753333333336, 0.3654333333333334, 0.8974650000000001, 0.4947283333333334, 0.9881900000000001, 2.3479155, 0.9894566666666668, 0.6734613333333335, 0.8054733333333335, 0.3654333333333334, 0.7915653333333335, 0.358986, 0.7552753333333334, 0.49795200000000006, 0.9894566666666668, 1.6635165000000003, 0.4815106666666667, 1.4862686666666667, 1.5376066666666668, 1.8534405000000003, 0.5272880000000001, 0.2637706666666667, 0.9959040000000001, 0.3654333333333334, 0.7915653333333335, 1.193542, 0.35746600000000006, 1.178874, 0.5272880000000001, 1.6521988333333337, 0.810217, 0.270218, 0.8974650000000001, 0.358986, 0.6227566666666668, 0.4907446666666667, 0.2413886666666667, 0.2637706666666667, 0.810217, 1.5940335000000003, 0.7403540000000001, 0.4947283333333334, 0.7552753333333334).finished(); typedef boost::mpl::list, mat_fixture<2, 3>, mat_fixture<3, 3>> test_types; BOOST_FIXTURE_TEST_CASE_TEMPLATE(material_hyper_elastic_test, F, test_types, F) { BOOST_CHECK_CLOSE(F::lambda, E*nu/((1+nu)*(1-2*nu)), tol); BOOST_CHECK_CLOSE(F::mu, E/2/(1+nu), tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_elastic_tensor, F, test_types, F) { Eigen::Matrix Cv; F::VoigtConv::fourth_to_voigt(F::C, Cv); auto & E = F::Young; auto & nu = F::Poisson; auto prefactor = E/(1+nu)/(1-2*nu); const auto diff{vsize(F::DimM) - F::DimM}; for (Dim_t i = 0; i < F::DimM; ++i) { BOOST_CHECK_CLOSE(Cv(i, i ), prefactor*(1-nu), 1e-12); for (Dim_t j = i + 1; j < F::DimM; ++j ) { BOOST_CHECK_CLOSE(Cv(i, j), prefactor*nu, 1e-12); BOOST_CHECK_CLOSE(Cv(j, i), prefactor*nu, 1e-12); } } for (Dim_t i = 0; i < diff; ++i) { BOOST_CHECK_CLOSE(Cv(i+F::DimM, i+F::DimM), prefactor*(1-2*nu)/2, 1e-12); for (Dim_t j = i + 1; j < diff; ++j ) { BOOST_CHECK_CLOSE(Cv(i+F::DimM, j+F::DimM), 0, 1e-12); } for (Dim_t j = 0; j < diff; ++j) { BOOST_CHECK_CLOSE(Cv(i+F::DimM, j ), 0, 1e-12); BOOST_CHECK_CLOSE(Cv(i, j+F::DimM), 0, 1e-12); } } } //----------------------------------------------------------------------------// BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constitutive_law, F, test_types, F) { //using SecondArray = typename MaterialHyperElastic::SecondArray; //using FourthArray = typename MaterialHyperElastic::FourthArray; using SecondTens = typename MaterialHyperElastic::SecondTens; using FourthTens = typename MaterialHyperElastic::FourthTens; //auto F_t = getMap(add_dim(F::F)); //auto K_t = Eigen::TensorMap(const_cast(F::K.data()), F::DimM, F::DimM, F::DimM, F::DimM, 1, 1); //std::cout << F_t; //std::cout << K_t; this->add_pixel(Index_t{0}); this->initialize(); std::array sizes; for (auto & s: sizes) { s=1; } auto F_t = F::get_second_array_map(F::F.data(), sizes); //std::cout << "hello: " << std::endl << F_t << std::endl; SecondTens P_; FourthTens K_; auto P_t = F::get_second_array_map(P_.data(), sizes); auto K_t = F::get_fourth_array_map(K_.data(), sizes); for (Dim_t i = 0; i < F::DimM; ++i) { for (Dim_t j = 0; j < F::DimM; ++j) { P_(i, j) = 10*(i+1)+j+1; } } F::compute_First_Piola_Kirchhoff_stress(F_t, P_t, K_t); auto K_v = VoigtConversion::template fourth_to_voigt(K_); Real P_error = (P_-F::P).norm(); if (P_error >= tol) { std::cout << "First Piola-Kirhoff stress wrong:" << std::endl << "P ref = " << std::endl << F::P << std::endl << "P comp = " << std::endl << P_ << std::endl << "P diff = " << std::endl << F::P-P_ << std::endl<< std::endl; } BOOST_CHECK_LT(P_error, tol); Real K_error = (K_v - F::K).norm(); if (K_error >= tol) { std::cout << "Tangent wrong: "<< std::endl << "K ref = " << std::endl << F::K < + struct sys_fixture:public SystemBase { + sys_fixture() + :SystemBase(std_size, + std_nb_pix){} + + const static std::array std_size; + const static Index_t std_nb_pix; + }; + + template<> const std::array sys_fixture<2, 2>::std_size{5.5, 6.5}; + template<> const std::array sys_fixture<2, 3>::std_size{5.5, 6.5}; + template<> const std::array sys_fixture<3, 3>::std_size{5.5, 6.5, 7.5}; + template<> const Index_t<2> sys_fixture<2, 2>::std_nb_pix{7, 9}; + template<> const Index_t<2> sys_fixture<2, 3>::std_nb_pix{7, 9}; + template<> const Index_t<3> sys_fixture<3, 3>::std_nb_pix{7, 9, 11}; + //----------------------------------------------------------------------------// + using test_systems = boost::mpl::list, + sys_fixture<2, 3>, + sys_fixture<3, 3>>; + + //----------------------------------------------------------------------------// + BOOST_AUTO_TEST_CASE(fftfreq_test) { + auto length = 15.2; + auto n1{8}, n2{9}; + + Eigen::RowVectorXd ref1(n1); + ref1 << 0.0, 0.06578947368421052, 0.13157894736842105, 0.19736842105263158, + -0.2631578947368421, -0.19736842105263158, -0.13157894736842105, + -0.06578947368421052; + Eigen::RowVectorXd ref2(n2); + ref2 << 0.0, 0.06578947368421052, 0.13157894736842105, 0.19736842105263158, + 0.2631578947368421, -0.2631578947368421, -0.19736842105263158, + -0.13157894736842105, -0.06578947368421052; + + auto challenge1(Projection::fft_freqs(n1, length)); + auto challenge2(Projection::fft_freqs(n2, length)); + Eigen::RowVectorXd er1 = ref1-challenge1; + Real error1 = (ref1 - challenge1).norm(); + Real error2 = (ref2 - challenge2).norm(); + BOOST_CHECK_LT(error1, tol); + BOOST_CHECK_LT(error2, tol); + + } + + //----------------------------------------------------------------------------// + BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, F, test_systems, F) { + } + +} // muSpectre diff --git a/tests/tests.hh b/tests/tests.hh new file mode 100644 index 0000000..235d0c5 --- /dev/null +++ b/tests/tests.hh @@ -0,0 +1,34 @@ +/** + * file tests.hh + * + * @author Till Junge + * + * @date 10 May 2017 + * + * @brief common defs for tests + * + * @section LICENCE + * + * Copyright (C) 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. + */ + +namespace muSpectre { + + const Real tol = 1e-14*100; //it's in percent + +} // muSpectre