diff --git a/.travis.yml b/.travis.yml index c9e1716..26e9309 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,92 +1,98 @@ language: cpp dist: trusty env: matrix: fast_finish: true include: - os: linux addons: apt: sources: - ubuntu-toolchain-r-test packages: - - g++-4.9 - - libeigen3-dev - env: COMPILER=gcc GCC=4.9 + - g++-5 + env: COMPILER=gcc GCC=5 - os: linux addons: apt: sources: - ubuntu-toolchain-r-test packages: - - g++-5 - - libeigen3-dev - env: COMPILER=gcc GCC=5 + - g++-6 + env: COMPILER=gcc GCC=6 - os: linux addons: apt: sources: - ubuntu-toolchain-r-test packages: - - g++-6 - - libeigen3-dev - env: COMPILER=gcc GCC=6 + - g++-7 + env: COMPILER=gcc GCC=7 + - os: linux + addons: + apt: + sources: + - ubuntu-toolchain-r-test + - llvm-toolchain-trusty-6.0 + packages: + - clang-6.0 + env: COMPILER=clang CLANG=6.0 + - os: osx + osx_image: xcode8 + compiler: clang env: global: - MINCONDA_VERSION="latest" - MINCONDA_LINUX="Linux-x86_64" - MINCONDA_OSX="MacOSX-x86_64" before_install: - | # Configure build variables if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then if [[ "$COMPILER" == "gcc" ]]; then export CXX=g++-$GCC CC=gcc-$GCC; fi if [[ "$COMPILER" == "clang" ]]; then export CXX=clang++-$CLANG CC=clang-$CLANG; fi elif [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export CXX=clang++ CC=clang; fi install: - # Define the version of miniconda to download + # Set environment using Conda - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then MINCONDA_OS=$MINCONDA_LINUX; elif [[ "$TRAVIS_OS_NAME" == "osx" ]]; then MINCONDA_OS=$MINCONDA_OSX; fi - wget "http://repo.continuum.io/miniconda/Miniconda3-$MINCONDA_VERSION-$MINCONDA_OS.sh" -O miniconda.sh; - bash miniconda.sh -b -p $HOME/miniconda - - export PATH="$HOME/miniconda/bin:${PATH}" - - export INCLUDE_PATH="$HOME/miniconda/include:${INCLUDE_PATH}" + - export PATH="$HOME/miniconda/bin:$PATH" - hash -r - conda config --set always_yes yes --set changeps1 no - conda update -q conda - conda install cmake -c conda-forge - - conda install xtl -c conda-forge - conda install xsimd -c conda-forge - conda install xtensor -c conda-forge - # store relevant - - build_dir="${PWD}/build" - - example_dir="${PWD}/docs/examples" - # build path - - mkdir "${build_dir}" - - cd "${build_dir}" - # install catch - - git clone https://github.com/catchorg/Catch2.git - - cd Catch2 - - mkdir build - - cd build - - mkdir opt - - cmake .. -DCATCH_BUILD_TESTING=OFF -DCMAKE_INSTALL_PREFIX:PATH="${PWD}"/opt - - make install - - export INCLUDE_PATH="${PWD}"/opt/include:"${INCLUDE_PATH}" - # make test-cases - - cd "${build_dir}" - - cmake .. -DBUILD_TESTS=ON + - conda install python -c conda-forge + - conda install numpy -c conda-forge + - conda install pyxtensor -c conda-forge + - conda install catch2 -c conda-forge + - conda install eigen -c conda-forge + - conda install gmatelastic -c conda-forge + - conda install python-gmatelastic -c conda-forge + - conda install highfive -c conda-forge + # Build/install the library + - cmake . -DBUILD_TESTS=ON -DBUILD_EXAMPLES=ON - make + - sudo make install + - python setup.py build + - python setup.py install script: # run tests - - cd "${build_dir}/test" - - ./test + - ./test/main + # run examples + - ./docs/examples/statics_FixedDisplacements_LinearElastic_example + - python docs/examples/statics_FixedDisplacements_LinearElastic_example.py --no-plot + - ./docs/examples/statics_FixedDisplacements_LinearElastic_manual_partition + - python docs/examples/statics_FixedDisplacements_LinearElastic_manual_partition.py --no-plot diff --git a/CMakeLists.txt b/CMakeLists.txt index 9b3847a..c3635e7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,156 +1,118 @@ #==================================================================================================# # # # (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM # # # #==================================================================================================# -# initialization -# -------------- - -# required to specify the c++ standard cmake_minimum_required(VERSION 3.0) -# required for install -include(CMakePackageConfigHelpers) -include(GNUInstallDirs) +# Basic settings +# ============== -# options -option(PKGCONFIG "Build pkg-config" ON) -option(BUILD_TESTS "${PROJECT_NAME} test suite" OFF) +project(GooseFEM) -# project settings -# ---------------- +option(BUILD_TESTS "Build tests" OFF) +option(BUILD_EXAMPLES "Build examples" OFF) -# name -project(GooseFEM) +# Version +# ======= -# file that contains the version information -set(parse_version include/GooseFEM/config.h) - -# header files -set(headers - include/GooseFEM/config.h - include/GooseFEM/GooseFEM.h - include/GooseFEM/Mesh.hpp - include/GooseFEM/Mesh.h - include/GooseFEM/MeshTri3.hpp - include/GooseFEM/MeshTri3.h - include/GooseFEM/MeshQuad4.hpp - include/GooseFEM/MeshQuad4.h - include/GooseFEM/MeshHex8.hpp - include/GooseFEM/MeshHex8.h - include/GooseFEM/Element.hpp - include/GooseFEM/Element.h - include/GooseFEM/ElementQuad4.hpp - include/GooseFEM/ElementQuad4.h - include/GooseFEM/ElementQuad4Planar.hpp - include/GooseFEM/ElementQuad4Planar.h - include/GooseFEM/ElementQuad4Axisymmetric.hpp - include/GooseFEM/ElementQuad4Axisymmetric.h - include/GooseFEM/ElementHex8.hpp - include/GooseFEM/ElementHex8.h - include/GooseFEM/Vector.hpp - include/GooseFEM/Vector.h - include/GooseFEM/VectorPartitioned.hpp - include/GooseFEM/VectorPartitioned.h - include/GooseFEM/VectorPartitionedTyings.hpp - include/GooseFEM/VectorPartitionedTyings.h - include/GooseFEM/Matrix.hpp - include/GooseFEM/Matrix.h - include/GooseFEM/MatrixPartitioned.hpp - include/GooseFEM/MatrixPartitioned.h - include/GooseFEM/MatrixPartitionedTyings.hpp - include/GooseFEM/MatrixPartitionedTyings.h - include/GooseFEM/MatrixDiagonal.hpp - include/GooseFEM/MatrixDiagonal.h - include/GooseFEM/MatrixDiagonalPartitioned.hpp - include/GooseFEM/MatrixDiagonalPartitioned.h - include/GooseFEM/Iterate.hpp - include/GooseFEM/Iterate.h - include/GooseFEM/TyingsPeriodic.hpp - include/GooseFEM/TyingsPeriodic.h - include/GooseFEM/ParaView.hpp - include/GooseFEM/ParaView.h -) - -# automatically parse the version number -file(READ "${parse_version}" version) -string(REGEX MATCH "define[ \t]+GOOSEFEM_WORLD_VERSION[ \t]+([0-9]+)" _ "${version}") -set(world_version "${CMAKE_MATCH_1}") -string(REGEX MATCH "define[ \t]+GOOSEFEM_MAJOR_VERSION[ \t]+([0-9]+)" _ "${version}") -set(major_version "${CMAKE_MATCH_1}") -string(REGEX MATCH "define[ \t]+GOOSEFEM_MINOR_VERSION[ \t]+([0-9]+)" _ "${version}") -set(minor_version "${CMAKE_MATCH_1}") -set(GOOSEFEM_VERSION ${world_version}.${major_version}.${minor_version}) - -# print information to screen -message(STATUS "Building ${PROJECT_NAME} v${GOOSEFEM_VERSION}") - -# paths -# ----- - -set(GOOSEFEM_ROOT_DIR "${CMAKE_INSTALL_PREFIX}") -set(GOOSEFEM_INCLUDE_DIR "${CMAKE_INSTALL_PREFIX}/${INCLUDE_INSTALL_DIR}") -set(INCLUDE_INSTALL_DIR "${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}") -set(CMAKEPACKAGE_INSTALL_DIR "${CMAKE_INSTALL_DATADIR}/${PROJECT_NAME}") -set(PKGCONFIG_INSTALL_DIR "${CMAKE_INSTALL_DATADIR}/pkgconfig") -set(fcmake "GooseFEMConfig.cmake") -set(fpkg "GooseFEM.pc") - -# configure CMake -# --------------- - -configure_package_config_file( - ${CMAKE_CURRENT_SOURCE_DIR}/GooseFEMConfig.cmake.in - ${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfig.cmake - PATH_VARS GOOSEFEM_INCLUDE_DIR GOOSEFEM_ROOT_DIR - INSTALL_DESTINATION ${CMAKEPACKAGE_INSTALL_DIR} - NO_CHECK_REQUIRED_COMPONENTS_MACRO -) - -# install/build -# ------------- - -# build test cases -if(BUILD_TESTS) - add_subdirectory(test) -endif() +file( + STRINGS + "${CMAKE_CURRENT_SOURCE_DIR}/include/GooseFEM/config.h" + GooseFEM_version_defines + REGEX + "#define GOOSEFEM_VERSION_(MAJOR|MINOR|PATCH)") -# disable pkg-config for native Windows builds -if(WIN32 OR CMAKE_HOST_SYSTEM_NAME MATCHES Windows) - option(PKGCONFIG "Build pkg-config ${fpkg} file" OFF) -endif() +foreach(ver ${GooseFEM_version_defines}) + if(ver MATCHES "#define GOOSEFEM_VERSION_(MAJOR|MINOR|PATCH) +([^ ]+)$") + set(GOOSEFEM_VERSION_${CMAKE_MATCH_1} "${CMAKE_MATCH_2}" CACHE INTERNAL "") + endif() +endforeach() -# pkg-config -if(PKGCONFIG) - configure_file(${fpkg}.in ${fpkg} @ONLY) - install(FILES ${CMAKE_CURRENT_BINARY_DIR}/${fpkg} DESTINATION ${PKGCONFIG_INSTALL_DIR}) -endif() +set(GOOSEFEM_VERSION + ${GOOSEFEM_VERSION_MAJOR}.${GOOSEFEM_VERSION_MINOR}.${GOOSEFEM_VERSION_PATCH}) + +message(STATUS "Building GooseFEM v${GOOSEFEM_VERSION}") -# CMake -install(FILES ${CMAKE_CURRENT_BINARY_DIR}/${fcmake} DESTINATION ${CMAKEPACKAGE_INSTALL_DIR}) +# Set target +# ========== -# header files -install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/${headers} DESTINATION ${INCLUDE_INSTALL_DIR}) +find_package(xtensor REQUIRED) -# print information to screen -# --------------------------- +add_library(GooseFEM INTERFACE) + +target_include_directories(GooseFEM INTERFACE + $ + $) + +target_link_libraries(GooseFEM INTERFACE xtensor) + +# Installation +# ============ + +include(CMakePackageConfigHelpers) +include(GNUInstallDirs) + +install( + DIRECTORY + "${CMAKE_CURRENT_SOURCE_DIR}/include/" + DESTINATION + include) + +install( + TARGETS + GooseFEM + EXPORT + GooseFEM-targets) + +install( + EXPORT + GooseFEM-targets + FILE + GooseFEMTargets.cmake + DESTINATION + "${CMAKE_INSTALL_LIBDIR}/cmake/GooseFEM") + +set(_GOOSEFEM ${CMAKE_SIZEOF_VOID_P}) +unset(CMAKE_SIZEOF_VOID_P) + +write_basic_package_version_file( + "${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfigVersion.cmake" + VERSION + ${GOOSEFEM_VERSION} + COMPATIBILITY + AnyNewerVersion) + +set(CMAKE_SIZEOF_VOID_P ${_GOOSEFEM}) + +install( + FILES + "${CMAKE_CURRENT_SOURCE_DIR}/GooseFEMConfig.cmake" + "${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfigVersion.cmake" + DESTINATION + "${CMAKE_INSTALL_LIBDIR}/cmake/GooseFEM") + +configure_file( + "${CMAKE_CURRENT_SOURCE_DIR}/GooseFEM.pc.in" + "${CMAKE_CURRENT_BINARY_DIR}/GooseFEM.pc" + @ONLY) + +install( + FILES + "${CMAKE_CURRENT_BINARY_DIR}/GooseFEM.pc" + DESTINATION + "${CMAKE_INSTALL_LIBDIR}/pkgconfig/") + +# Add builds +# ========== + +include("GooseFEMConfig.cmake") -message(STATUS "") -message(STATUS "+---------------------------------------------------------------------------------") -message(STATUS "|") -message(STATUS "| Use 'make install' to install in") -message(STATUS "| ${CMAKE_INSTALL_PREFIX}") if(BUILD_TESTS) - message(STATUS "|") - message(STATUS "| Use 'make' and './test/test' to run the tests") + add_subdirectory(test) +endif() + +if(BUILD_EXAMPLES) + add_subdirectory(docs/examples) endif() -message(STATUS "|") -message(STATUS "| To specify a custom directory call") -message(STATUS "| cmake /path/to/${PROJECT_NAME} -DCMAKE_INSTALL_PREFIX=yourprefix") -message(STATUS "|") -message(STATUS "| For custom paths, add the following line to your '~/.bashrc'") -message(STATUS "| export PKG_CONFIG_PATH=${CMAKE_INSTALL_PREFIX}/share/pkgconfig:$PKG_CONFIG_PATH") -message(STATUS "|") -message(STATUS "+---------------------------------------------------------------------------------") -message(STATUS "") diff --git a/GooseFEM.pc.in b/GooseFEM.pc.in index 0c9ae5f..4ba0cf4 100644 --- a/GooseFEM.pc.in +++ b/GooseFEM.pc.in @@ -1,9 +1,7 @@ prefix=@CMAKE_INSTALL_PREFIX@ -exec_prefix=${prefix} +includedir=${prefix}/include -Name: GooseFEM -Description: Several Finite Element simulation, and some simple meshes and operations (C++/Python) -Requires: +Name: @PROJECT_NAME@ +Description: Finite Element simulation, some simple meshes and operations. Version: @GOOSEFEM_VERSION@ -Libs: -Cflags: -I${prefix}/@GOOSEFEM_INCLUDE_DIR@ +Cflags: -I${includedir} diff --git a/GooseFEMConfig.cmake b/GooseFEMConfig.cmake new file mode 100644 index 0000000..abe8508 --- /dev/null +++ b/GooseFEMConfig.cmake @@ -0,0 +1,76 @@ +# GooseFEM cmake module +# +# This module sets the target: +# +# GooseFEM +# +# In addition, it sets the following variables: +# +# GooseFEM_FOUND - true if GooseFEM found +# GooseFEM_VERSION - GooseFEM's version +# GooseFEM_INCLUDE_DIRS - the directory containing GooseFEM headers +# +# The following support targets are defined to simplify things: +# +# GooseFEM::compiler_warnings - enable compiler warnings +# GooseFEM::assert - enable GooseFEM assertions +# GooseFEM::debug - enable all assertions (slow) + +include(CMakeFindDependencyMacro) + +# Define target "GooseFEM" +# ======================== + +if(NOT TARGET GooseFEM) + include("${CMAKE_CURRENT_LIST_DIR}/GooseFEMTargets.cmake") + get_target_property(GooseFEM_INCLUDE_DIRS GooseFEM INTERFACE_INCLUDE_DIRECTORIES) +endif() + +# Find dependencies +# ================= + +find_dependency(xtensor) + +find_package(Eigen3 QUIET) + +if(NOT Eigen3_FOUND) + find_package(PkgConfig) + pkg_check_modules(EIGEN3 QUIET eigen3) +endif() + +if(Eigen3_FOUND) + target_include_directories(GooseFEM INTERFACE ${EIGEN3_INCLUDE_DIRS}) +endif() + +# Define support target "GooseFEM::compiler_warnings" +# =================================================== + +if(NOT TARGET GooseFEM::compiler_warnings) + add_library(GooseFEM::compiler_warnings INTERFACE IMPORTED) + if(MSVC) + target_compile_options(GooseFEM::compiler_warnings INTERFACE + /W4) + else() + target_compile_options(GooseFEM::compiler_warnings INTERFACE + -Wall + -Wextra + -pedantic + -Wno-unknown-pragmas) + endif() +endif() + +# Define support target "GooseFEM::assert" +# ======================================== + +if(NOT TARGET GooseFEM::assert) + add_library(GooseFEM::assert INTERFACE IMPORTED) + target_compile_definitions(GooseFEM::assert INTERFACE GOOSEFEM_ENABLE_ASSERT) +endif() + +# Define support target "GooseEYE::debug" +# ======================================= + +if(NOT TARGET GooseFEM::debug) + add_library(GooseFEM::debug INTERFACE IMPORTED) + target_compile_definitions(GooseFEM::debug INTERFACE GOOSEFEM_DEBUG) +endif() diff --git a/GooseFEMConfig.cmake.in b/GooseFEMConfig.cmake.in deleted file mode 100644 index 6bc7e68..0000000 --- a/GooseFEMConfig.cmake.in +++ /dev/null @@ -1,4 +0,0 @@ -set(GOOSEFEM_FOUND 1) - -set(GOOSEFEM_INCLUDE_DIR ${CMAKE_CURRENT_SOURCE_DIR}) -set(GOOSEFEM_INCLUDE_DIRS ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/docs/examples/CMakeLists.txt b/docs/examples/CMakeLists.txt new file mode 100644 index 0000000..cb06ae1 --- /dev/null +++ b/docs/examples/CMakeLists.txt @@ -0,0 +1,60 @@ + +cmake_minimum_required(VERSION 3.0) + +if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR) + project(GooseFEM-examples) + find_package(GooseFEM REQUIRED CONFIG) +endif() + +option(WARNINGS "Show build warnings" ON) +option(SIMD "Enable xsimd" ON) +option(DEBUG "Enable all assertions" ON) + +find_package(HDF5 REQUIRED) + +add_library(libraries INTERFACE IMPORTED) + +target_link_libraries(libraries INTERFACE GooseFEM ${HDF5_LIBRARIES}) + +if (SIMD) + target_link_libraries(libraries INTERFACE xtensor::optimize xtensor::use_xsimd) +endif() + +if (WARNINGS) + target_link_libraries(libraries INTERFACE GooseFEM::compiler_warnings) +endif() + +if (DEBUG) + target_link_libraries(libraries INTERFACE GooseFEM::debug) +endif() + + +# create executable + +set(exec "statics_FixedDisplacements_LinearElastic_example") +set(source "statics/FixedDisplacements_LinearElastic/example.cpp") + +add_executable(${exec} ${source}) + +target_link_libraries(${exec} PRIVATE libraries) + +configure_file( + "statics/FixedDisplacements_LinearElastic/example.py" + "statics_FixedDisplacements_LinearElastic_example.py" COPYONLY) + +configure_file( + "statics/FixedDisplacements_LinearElastic/plot.py" + "statics_FixedDisplacements_LinearElastic_plot.py" COPYONLY) + +# create executable + +set(exec "statics_FixedDisplacements_LinearElastic_manual_partition") +set(source "statics/FixedDisplacements_LinearElastic/manual_partition.cpp") + +add_executable(${exec} ${source}) + +target_link_libraries(${exec} PRIVATE libraries) + +configure_file( + "statics/FixedDisplacements_LinearElastic/manual_partition.py" + "statics_FixedDisplacements_LinearElastic_manual_partition.py" COPYONLY) diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/example/main.cpp b/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp similarity index 98% rename from docs/examples/statics/FixedDisplacements_LinearElastic/example/main.cpp rename to docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp index 8f8a101..d2360b4 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/example/main.cpp +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp @@ -1,139 +1,139 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5,5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLeft = mesh.nodesLeftEdge(); xt::xtensor nodesRight = mesh.nodesRightEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBottom = mesh.nodesBottomEdge(); // fixed displacements DOFs // ------------------------ xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesRight ), 0), xt::view(dofs, xt::keep(nodesTop ), 1), xt::view(dofs, xt::keep(nodesLeft ), 0), xt::view(dofs, xt::keep(nodesBottom), 1) )); // simulation variables // -------------------- // vector definition GooseFEM::VectorPartitioned vector(conn, dofs, iip); // allocate system matrix GooseFEM::MatrixPartitioned<> K(conn, dofs, iip); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // element/material definition // --------------------------- // element definition GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material definition GMatElastic::Cartesian3d::Matrix mat(nelem, nip, 1., 1.); // integration point tensors size_t d = 3; xt::xtensor Eps = xt::empty({nelem, nip, d, d }); xt::xtensor Sig = xt::empty({nelem, nip, d, d }); xt::xtensor C = xt::empty({nelem, nip, d, d, d, d}); // solve // ----- // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // set fixed displacements xt::view(disp, xt::keep(nodesRight ), 0) = +0.1; xt::view(disp, xt::keep(nodesTop ), 1) = -0.1; xt::view(disp, xt::keep(nodesLeft ), 0) = 0.0; xt::view(disp, xt::keep(nodesBottom), 1) = 0.0; // residual xt::noalias(fres) = fext - fint; // solve K.solve(fres, disp); // post-process // ------------ // compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.copy_p(fint, fext); // residual xt::noalias(fres) = fext - fint; // print residual std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl; // average stress per node xt::xtensor dV = elem.DV(2); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output - H5Easy::File file("main.h5", H5Easy::File::Overwrite); + H5Easy::File file("output.h5", H5Easy::File::Overwrite); H5Easy::dump(file, "/coor", coor); H5Easy::dump(file, "/conn", conn); H5Easy::dump(file, "/disp", disp); H5Easy::dump(file, "/Sig" , SigAv); return 0; } diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/example/main.py b/docs/examples/statics/FixedDisplacements_LinearElastic/example.py similarity index 91% rename from docs/examples/statics/FixedDisplacements_LinearElastic/example/main.py rename to docs/examples/statics/FixedDisplacements_LinearElastic/example.py index de064fa..8d98de8 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/example/main.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/example.py @@ -1,164 +1,173 @@ import GooseFEM import GMatElastic import numpy as np # mesh # ---- # define mesh mesh = GooseFEM.Mesh.Quad4.Regular(5,5) # mesh dimensions nelem = mesh.nelem() nne = mesh.nne() ndim = mesh.ndim() # mesh definitions coor = mesh.coor() conn = mesh.conn() dofs = mesh.dofs() # node sets nodesLeft = mesh.nodesLeftEdge() nodesRight = mesh.nodesRightEdge() nodesTop = mesh.nodesTopEdge() nodesBottom = mesh.nodesBottomEdge() # fixed displacements DOFs # ------------------------ iip = np.concatenate(( - dofs[nodesRight , 0], - dofs[nodesTop , 1], - dofs[nodesLeft , 0], - dofs[nodesBottom, 1] + dofs[nodesRight , 0], + dofs[nodesTop , 1], + dofs[nodesLeft , 0], + dofs[nodesBottom, 1] )) # simulation variables # -------------------- # vector definition vector = GooseFEM.VectorPartitioned(conn, dofs, iip) # allocate system matrix K = GooseFEM.MatrixPartitioned(conn, dofs, iip) # nodal quantities disp = np.zeros(coor.shape) fint = np.zeros(coor.shape) fext = np.zeros(coor.shape) fres = np.zeros(coor.shape) # element vectors ue = np.empty((nelem, nne, ndim)) fe = np.empty((nelem, nne, ndim)) Ke = np.empty((nelem, nne*ndim, nne*ndim)) # element/material definition # --------------------------- # element definition elem = GooseFEM.Element.Quad4.QuadraturePlanar(vector.AsElement(coor)) nip = elem.nip() # material definition mat = GMatElastic.Cartesian3d.Matrix(nelem, nip, 1., 1.) # integration point tensors d = 3 Eps = np.empty((nelem, nip, d, d )) Sig = np.empty((nelem, nip, d, d )) C = np.empty((nelem, nip, d, d, d, d)) # solve # ----- # strain ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) # stress & tangent (Sig, C) = mat.Tangent(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # stiffness matrix Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C) K.assemble(Ke) # set fixed displacements disp[nodesRight , 0] = +0.1 disp[nodesTop , 1] = -0.1 disp[nodesLeft , 0] = 0.0 disp[nodesBottom, 1] = 0.0 # residual fres = fext - fint # solve disp = K.Solve(fres, disp) # post-process # ------------ # compute strain and stress ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) Sig = mat.Stress(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # apply reaction force fext = vector.Copy_p(fint, fext) # residual fres = fext - fint # print residual print(np.sum(np.abs(fres)) / np.sum(np.abs(fext))) # average stress per node dV = elem.DV(2) Sig = np.average(Sig, weights=dV, axis=1) +# skip plot with "--no-plot" command line argument +# ------------------------------------------------ + +import sys + +if len(sys.argv) == 2: + if sys.argv[1] == "--no-plot": + sys.exit(0) + # plot # ---- import matplotlib.pyplot as plt import GooseMPL as gplt plt.style.use(['goose', 'goose-latex']) # extract dimension nelem = conn.shape[0] # tensor products ddot22 = lambda A2,B2: np.einsum('eij ,eji->e ',A2,B2) ddot42 = lambda A4,B2: np.einsum('eijkl,elk->eij ',A4,B2) dyad22 = lambda A2,B2: np.einsum('eij ,ekl->eijkl',A2,B2) # identity tensor (single tensor) i = np.eye(3) # identity tensors (grid) I = np.einsum('ij ,e' , i ,np.ones([nelem])) I4 = np.einsum('ijkl,e->eijkl',np.einsum('il,jk',i,i),np.ones([nelem])) I4rt = np.einsum('ijkl,e->eijkl',np.einsum('ik,jl',i,i),np.ones([nelem])) I4s = (I4+I4rt)/2. II = dyad22(I,I) I4d = I4s-II/3. # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3./2.*ddot22(Sigd,Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor+disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0,0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) -plt.savefig('main.pdf') +plt.savefig('plot.pdf') plt.close() diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/example/CMakeLists.txt b/docs/examples/statics/FixedDisplacements_LinearElastic/example/CMakeLists.txt deleted file mode 100644 index 0cf69f2..0000000 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/example/CMakeLists.txt +++ /dev/null @@ -1,64 +0,0 @@ -cmake_minimum_required(VERSION 3.1) - -# basic info -# ---------- - -# define a project name -project(main) - -# define empty list of libraries to link -set(PROJECT_LIBS "") - -# copy files to build folder -configure_file("plot.py" "plot.py" COPYONLY) - -# basic compiler options -# ---------------------- - -# set optimization level -set(CMAKE_BUILD_TYPE Release) - -# set C++ standard -set(CMAKE_CXX_STANDARD 14) -set(CMAKE_CXX_STANDARD_REQUIRED ON) - -# set warnings on -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic -g") - -# enable optimisation -# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native") -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DXTENSOR_USE_XSIMD=ON") - -# assertions -if(NOT ASSERT) - add_definitions(-DNDEBUG) -else() - add_definitions(-DXTENSOR_ENABLE_ASSERT=ON) -endif() - -# load packages -# ------------- - -# add custom include paths -if(NOT "$ENV{INCLUDE_PATH}" STREQUAL "") - string(REPLACE ":" ";" INCLUDE_LIST "$ENV{INCLUDE_PATH}") - include_directories(${INCLUDE_LIST}) -endif() - -# load HDF5 -find_package(HDF5 COMPONENTS CXX REQUIRED) -include_directories(${HDF5_INCLUDE_DIRS}) -set(PROJECT_LIBS ${PROJECT_LIBS} ${HDF5_LIBS} ${HDF5_LIBRARIES}) - -# load eigen -find_package(PkgConfig) -pkg_check_modules(EIGEN3 REQUIRED eigen3) -include_directories(SYSTEM ${EIGEN3_INCLUDE_DIRS}) - -# create executable -# ----------------- - -add_executable(${PROJECT_NAME} main.cpp) - -target_link_libraries(${PROJECT_NAME} ${PROJECT_LIBS}) - diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/main.cpp b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp similarity index 98% rename from docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/main.cpp rename to docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp index d27743a..012d907 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/main.cpp +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp @@ -1,155 +1,155 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5,5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLeft = mesh.nodesLeftEdge(); xt::xtensor nodesRight = mesh.nodesRightEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBottom = mesh.nodesBottomEdge(); // fixed displacements DOFs // ------------------------ xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesRight ), 0), xt::view(dofs, xt::keep(nodesTop ), 1), xt::view(dofs, xt::keep(nodesLeft ), 0), xt::view(dofs, xt::keep(nodesBottom), 1) )); // simulation variables // -------------------- // vector definition GooseFEM::VectorPartitioned vector(conn, dofs, iip); // allocate system matrix GooseFEM::MatrixPartitioned<> K(conn, dofs, iip); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // DOF values xt::xtensor u_u = xt::zeros({vector.nnu()}); xt::xtensor fres_u = xt::zeros({vector.nnu()}); xt::xtensor fext_p = xt::zeros({vector.nnp()}); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // element/material definition // --------------------------- // element definition GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material definition GMatElastic::Cartesian3d::Matrix mat(nelem, nip, 1., 1.); // integration point tensors size_t d = 3; xt::xtensor Eps = xt::empty({nelem, nip, d, d }); xt::xtensor Sig = xt::empty({nelem, nip, d, d }); xt::xtensor C = xt::empty({nelem, nip, d, d, d, d}); // solve // ----- // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // set fixed displacements xt::xtensor u_p = xt::concatenate(xt::xtuple( +0.1 * xt::ones({nodesRight .size()}), -0.1 * xt::ones({nodesTop .size()}), 0.0 * xt::ones({nodesLeft .size()}), 0.0 * xt::ones({nodesBottom.size()}) )); // residual xt::noalias(fres) = fext - fint; // partition vector.asDofs_u(fres, fres_u); // solve K.solve_u(fres_u, u_p, u_u); // assemble to nodal vector vector.asNode(u_u, u_p, disp); // post-process // ------------ // compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.asDofs_p(fint, fext_p); // residual xt::noalias(fres) = fext - fint; // partition vector.asDofs_u(fres, fres_u); // print residual std::cout << xt::sum(xt::abs(fres_u))[0] / xt::sum(xt::abs(fext_p))[0] << std::endl; // average stress per node xt::xtensor dV = elem.DV(2); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output - H5Easy::File file("main.h5", H5Easy::File::Overwrite); + H5Easy::File file("output.h5", H5Easy::File::Overwrite); H5Easy::dump(file, "/coor", coor); H5Easy::dump(file, "/conn", conn); H5Easy::dump(file, "/disp", disp); H5Easy::dump(file, "/Sig" , SigAv); return 0; } diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/main.py b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py similarity index 87% rename from docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/main.py rename to docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py index 2610e98..154e6ef 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/main.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py @@ -1,179 +1,189 @@ import GooseFEM import GMatElastic import numpy as np # mesh # ---- # define mesh mesh = GooseFEM.Mesh.Quad4.Regular(5,5) # mesh dimensions nelem = mesh.nelem() nne = mesh.nne() ndim = mesh.ndim() # mesh definitions coor = mesh.coor() conn = mesh.conn() dofs = mesh.dofs() # node sets nodesLeft = mesh.nodesLeftEdge() nodesRight = mesh.nodesRightEdge() nodesTop = mesh.nodesTopEdge() nodesBottom = mesh.nodesBottomEdge() # fixed displacements DOFs # ------------------------ iip = np.concatenate(( - dofs[nodesRight , 0], - dofs[nodesTop , 1], - dofs[nodesLeft , 0], - dofs[nodesBottom, 1] + dofs[nodesRight , 0], + dofs[nodesTop , 1], + dofs[nodesLeft , 0], + dofs[nodesBottom, 1] )) # simulation variables # -------------------- # vector definition vector = GooseFEM.VectorPartitioned(conn, dofs, iip) # allocate system matrix K = GooseFEM.MatrixPartitioned(conn, dofs, iip) # nodal quantities disp = np.zeros(coor.shape) fint = np.zeros(coor.shape) fext = np.zeros(coor.shape) fres = np.zeros(coor.shape) # DOF values fext_p = np.zeros(vector.nnp()) fres_u = np.zeros(vector.nnu()) fext_p = np.zeros(vector.nnp()) # element vectors ue = np.empty((nelem, nne, ndim)) fe = np.empty((nelem, nne, ndim)) Ke = np.empty((nelem, nne*ndim, nne*ndim)) # element/material definition # --------------------------- # element definition elem = GooseFEM.Element.Quad4.QuadraturePlanar(vector.AsElement(coor)) nip = elem.nip() # material definition mat = GMatElastic.Cartesian3d.Matrix(nelem, nip, 1., 1.) # integration point tensors d = 3 Eps = np.empty((nelem, nip, d, d )) Sig = np.empty((nelem, nip, d, d )) C = np.empty((nelem, nip, d, d, d, d)) # solve # ----- # strain ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) # stress & tangent (Sig, C) = mat.Tangent(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # stiffness matrix Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C) K.assemble(Ke) # set fixed displacements u_p = np.concatenate(( - +0.1 * np.ones(nodesRight .size), - -0.1 * np.ones(nodesTop .size), - 0.0 * np.ones(nodesLeft .size), - 0.0 * np.ones(nodesBottom.size) + +0.1 * np.ones(nodesRight .size), + -0.1 * np.ones(nodesTop .size), + 0.0 * np.ones(nodesLeft .size), + 0.0 * np.ones(nodesBottom.size) )) # residual fres = fext - fint # partition fres_u = vector.AsDofs_u(fres) # solve u_u = K.Solve_u(fres_u, u_p) # assemble to nodal vector disp = vector.AsNode(u_u, u_p) # post-process # ------------ # compute strain and stress ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) Sig = mat.Stress(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # apply reaction force fext_p = vector.AsDofs_p(fint) # residual fres = fext - fint # partition fres_u = vector.AsDofs_u(fres); # print residual print(np.sum(np.abs(fres_u)) / np.sum(np.abs(fext_p))) # average stress per node dV = elem.DV(2) Sig = np.average(Sig, weights=dV, axis=1) +# skip plot with "--no-plot" command line argument +# ------------------------------------------------ + +import sys + +if len(sys.argv) == 2: + if sys.argv[1] == "--no-plot": + sys.exit(0) + # plot # ---- import matplotlib.pyplot as plt import GooseMPL as gplt plt.style.use(['goose', 'goose-latex']) # extract dimension nelem = conn.shape[0] # tensor products ddot22 = lambda A2,B2: np.einsum('eij ,eji->e ',A2,B2) ddot42 = lambda A4,B2: np.einsum('eijkl,elk->eij ',A4,B2) dyad22 = lambda A2,B2: np.einsum('eij ,ekl->eijkl',A2,B2) # identity tensor (single tensor) i = np.eye(3) # identity tensors (grid) I = np.einsum('ij ,e' , i ,np.ones([nelem])) I4 = np.einsum('ijkl,e->eijkl',np.einsum('il,jk',i,i),np.ones([nelem])) I4rt = np.einsum('ijkl,e->eijkl',np.einsum('ik,jl',i,i),np.ones([nelem])) I4s = (I4+I4rt)/2. II = dyad22(I,I) I4d = I4s-II/3. # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3./2.*ddot22(Sigd,Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor+disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0,0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) -plt.show() +plt.savefig('plot.pdf') +plt.close() diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/CMakeLists.txt b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/CMakeLists.txt deleted file mode 100644 index 0cf69f2..0000000 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/CMakeLists.txt +++ /dev/null @@ -1,64 +0,0 @@ -cmake_minimum_required(VERSION 3.1) - -# basic info -# ---------- - -# define a project name -project(main) - -# define empty list of libraries to link -set(PROJECT_LIBS "") - -# copy files to build folder -configure_file("plot.py" "plot.py" COPYONLY) - -# basic compiler options -# ---------------------- - -# set optimization level -set(CMAKE_BUILD_TYPE Release) - -# set C++ standard -set(CMAKE_CXX_STANDARD 14) -set(CMAKE_CXX_STANDARD_REQUIRED ON) - -# set warnings on -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic -g") - -# enable optimisation -# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native") -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DXTENSOR_USE_XSIMD=ON") - -# assertions -if(NOT ASSERT) - add_definitions(-DNDEBUG) -else() - add_definitions(-DXTENSOR_ENABLE_ASSERT=ON) -endif() - -# load packages -# ------------- - -# add custom include paths -if(NOT "$ENV{INCLUDE_PATH}" STREQUAL "") - string(REPLACE ":" ";" INCLUDE_LIST "$ENV{INCLUDE_PATH}") - include_directories(${INCLUDE_LIST}) -endif() - -# load HDF5 -find_package(HDF5 COMPONENTS CXX REQUIRED) -include_directories(${HDF5_INCLUDE_DIRS}) -set(PROJECT_LIBS ${PROJECT_LIBS} ${HDF5_LIBS} ${HDF5_LIBRARIES}) - -# load eigen -find_package(PkgConfig) -pkg_check_modules(EIGEN3 REQUIRED eigen3) -include_directories(SYSTEM ${EIGEN3_INCLUDE_DIRS}) - -# create executable -# ----------------- - -add_executable(${PROJECT_NAME} main.cpp) - -target_link_libraries(${PROJECT_NAME} ${PROJECT_LIBS}) - diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/plot.py b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/plot.py deleted file mode 100644 index 0bf1209..0000000 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition/plot.py +++ /dev/null @@ -1,51 +0,0 @@ - -import h5py - -import matplotlib.pyplot as plt -import GooseMPL as gplt -import numpy as np - -plt.style.use(['goose', 'goose-latex']) - -# open file -file = h5py.File('main.h5', 'r') - -# read fields -coor = file['/coor'][...] -conn = file['/conn'][...] -disp = file['/disp'][...] -Sig = file['/Sig' ][...] - -# extract dimension -nelem = conn.shape[0] - -# tensor products -ddot22 = lambda A2,B2: np.einsum('eij ,eji->e ',A2,B2) -ddot42 = lambda A4,B2: np.einsum('eijkl,elk->eij ',A4,B2) -dyad22 = lambda A2,B2: np.einsum('eij ,ekl->eijkl',A2,B2) - -# identity tensor (single tensor) -i = np.eye(3) - -# identity tensors (grid) -I = np.einsum('ij ,e' , i ,np.ones([nelem])) -I4 = np.einsum('ijkl,e->eijkl',np.einsum('il,jk',i,i),np.ones([nelem])) -I4rt = np.einsum('ijkl,e->eijkl',np.einsum('ik,jl',i,i),np.ones([nelem])) -I4s = (I4+I4rt)/2. -II = dyad22(I,I) -I4d = I4s-II/3. - -# compute equivalent stress -Sigd = ddot42(I4d, Sig) -sigeq = np.sqrt(3./2.*ddot22(Sigd,Sigd)) - -# plot - -fig, ax = plt.subplots() - -gplt.patch(coor=coor+disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0,0.1)) - -gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) - -plt.show() - diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/example/plot.py b/docs/examples/statics/FixedDisplacements_LinearElastic/plot.py similarity index 87% rename from docs/examples/statics/FixedDisplacements_LinearElastic/example/plot.py rename to docs/examples/statics/FixedDisplacements_LinearElastic/plot.py index a5b511a..672914c 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/example/plot.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/plot.py @@ -1,49 +1,47 @@ import h5py import matplotlib.pyplot as plt import GooseMPL as gplt import numpy as np plt.style.use(['goose', 'goose-latex']) -# open file -file = h5py.File('main.h5', 'r') - # read fields -coor = file['/coor'][...] -conn = file['/conn'][...] -disp = file['/disp'][...] -Sig = file['/Sig' ][...] +with h5py.File('output.h5', 'r') as data: + coor = data['/coor'][...] + conn = data['/conn'][...] + disp = data['/disp'][...] + Sig = data['/Sig' ][...] # extract dimension nelem = conn.shape[0] # tensor products ddot22 = lambda A2,B2: np.einsum('eij ,eji->e ',A2,B2) ddot42 = lambda A4,B2: np.einsum('eijkl,elk->eij ',A4,B2) dyad22 = lambda A2,B2: np.einsum('eij ,ekl->eijkl',A2,B2) # identity tensor (single tensor) i = np.eye(3) # identity tensors (grid) I = np.einsum('ij ,e' , i ,np.ones([nelem])) I4 = np.einsum('ijkl,e->eijkl',np.einsum('il,jk',i,i),np.ones([nelem])) I4rt = np.einsum('ijkl,e->eijkl',np.einsum('ik,jl',i,i),np.ones([nelem])) I4s = (I4+I4rt)/2. II = dyad22(I,I) I4d = I4s-II/3. # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3./2.*ddot22(Sigd,Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor+disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0,0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) plt.savefig('plot.pdf') plt.close() diff --git a/docs/examples/statics_linear.rst b/docs/examples/statics_linear.rst index cc3170a..ad9a626 100644 --- a/docs/examples/statics_linear.rst +++ b/docs/examples/statics_linear.rst @@ -1,311 +1,307 @@ ************** Linear elastic ************** Consider a uniform linear elastic bar that is extended by a uniform fixed displacement on both sides. This problem can be modelled and discretised using symmetry as show below. In this example we will furthermore assume that the bar is sufficiently thick in the out-of-plane direction to be modelled using two-dimensional plane strain. .. image:: statics/FixedDisplacements_LinearElastic/problem-sketch.svg :width: 300px :align: center | Below an example is described line-by-line. The full example can be downloaded: -| :download:`main.cpp ` -| :download:`CMakeLists.txt ` -| :download:`plot.py ` +| :download:`example.cpp ` +| :download:`plot.py ` .. todo:: Compile and run instructions. .. note:: - This example is also available using the Python interface (:download:`main.py `). Compared to the C++ API, the Python API requires more data-allocation, in particular for the functions "AsElement" and "AssembleNode". See: :ref:`conventions_allocation`. + This example is also available using the Python interface (:download:`example.py `). Compared to the C++ API, the Python API requires more data-allocation, in particular for the functions "AsElement" and "AssembleNode". See: :ref:`conventions_allocation`. Include library =============== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 1-4 :emphasize-lines: 1-2 The first step is to include the header-only library. Note that for this example we also make use of a material model (`GMatElastic `_) and a library to write (and read) HDF5 files (`HighFive `_). Define mesh =========== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 11-28 :emphasize-lines: 2 A mesh is defined using GooseFEM. As observed the "mesh" is a class that has methods to extract the relevant information such as the nodal coordinates ("coor"), the connectivity ("conn"), the degrees-of-freedom per node ("dofs") and several node-sets that will be used to impose the sketched boundary conditions ("nodesLeft", "nodesRight", "nodesTop", "nodesBottom"). Note that: * The connectivity ("conn") contains information of which nodes, in which order, belong to which element. * The degrees-of-freedom per node ("dofs") contains information of how a nodal vector (a vector stored per node) can be transformed to a list of degrees-of-freedom as used in the linear system (although this can be mostly done automatically as we will see below). .. seealso:: * :ref:`conventions_terminology` * Details: :ref:`MeshQuad4` Define partitioning =================== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 33-38 We will reorder such that degrees-of-freedom are ordered such that .. math:: \texttt{u} = \begin{bmatrix} \texttt{u}_u \\ \texttt{u}_p \end{bmatrix} where the subscript :math:`u` and :math:`p` respectively denote *Unknown* and *Prescribed* degrees-of-freedom. To achieve this we start by collecting all prescribed degrees-of-freedom in "iip". (Avoid) Book-keeping ==================== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 44 :emphasize-lines: 1 To switch between the three of GooseFEM's data-representations, an instance of the "Vector" class is used. This instance, "vector", will enable us to switch between a vector field (e.g. the displacement) 1. collected per node, 2. collected per degree-of-freedom, or 3. collected per element. .. note:: The "Vector" class collects most, if not all, the burden of book-keeping. It is thus here that "conn", "dofs", and "iip" are used. In particular, * 'nodevec' :math:`\leftrightarrow` 'dofval' using "dofs" and "iip", * 'nodevec' :math:`\leftrightarrow` 'elemvec' using "conn". By contrast, most of GooseFEM's other methods receive the relevant representation, and consequently require no problem specific knowledge. They thus do not have to supplied with "conn", "dofs", or "iip". .. seealso:: * :ref:`conventions_vector` * :ref:`conventions_storage` * Details: :ref:`Vector` System matrix ============= -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 47 :emphasize-lines: 1 We now also allocate the system/stiffness system (stored as sparse matrix). Like vector, it can accept and return different vector representations, in addition to the ability to assemble from element system matrices. .. note:: Here, the default solver is used (which is the default template, hence the "<>"). To use other solvers see: :ref:`linear_solver`. .. seealso:: * :ref:`conventions_matrix` * Details: :ref:`Matrix` Allocate nodal vectors ====================== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 50-53 * "disp": nodal displacements * "fint": nodal internal forces * "fext": nodal external forces * "fres": nodal residual forces Allocate element vectors ======================== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 56-58 * "ue": displacement * "fe": force * "Ke": tangent matrix .. warning:: Upsizing (e.g. "disp" :math:`\rightarrow` "ue") can be done uniquely, but downsizing (e.g. "fe" :math:`\rightarrow` "fint") can be done in more than one way, see :ref:`conventions_vector_conversion`. We will get back to this point below. Element definition ================== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 64-65 At this moment the interpolation and quadrature is allocated. The shape functions and integration points (that can be customised) are stored in this class. As observed, no further information is needed than the number of elements and the nodal coordinates per element. Both are contained in the output of "vector.AsElement(coor)", which is an 'elemvec' of shape "[nelem, nne, ndim]". This illustrates that problem specific book-keeping is isolated to the main program, using "Vector" as tool. .. note:: The shape functions are computed when constructing this class, they are not recomputed when evaluating them. One can recompute them if the nodal coordinates change using ".update_x(...)", however, this is only relevant in a large deformation setting. .. seealso:: * :ref:`conventions_vector` * :ref:`conventions_storage` * Details: :ref:`Vector` * Details: :ref:`ElementQuad4` Material definition =================== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 68 We now define a uniform linear elastic material, using an external library that is tuned to GooseFEM. This material library will translate a strain tensor per integration point to a stress tensor per integration point and a stiffness tensor per integration point. .. seealso:: Material libraries tuned to GooseFEM include: * `GMatElastic `__ * `GMatElastoPlastic `__ * `GMatElastoPlasticFiniteStrainSimo `__ * `GMatElastoPlasticQPot `__ * `GMatElastoPlasticQPot3d `__ * `GMatNonLinearElastic `__ But other libraries can also be easily used with (simple) wrappers. Integration point tensors ========================= -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 71-74 These strain, stress, and stiffness tensors per integration point are now allocated. Note that these tensors are 3-d while our problem was 2-d. This is thanks to the plane strain assumption, and the element definition that ignores all third-dimension components. Compute strain ============== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 80-81 The strain per integration point is now computed using the current nodal displacements (stored as 'elemvec' in "ue") and the gradient of the shape functions. .. note:: "ue" is the output of "vector.asElement(disp, ue)". Using this syntax re-allocation of "ue" is avoided. If this optimisation is irrelevant for you problem (or if you are using the Python interface), please use the same function, but starting with a capital: .. code-block:: cpp ue = vector.AsElement(disp); Note that this allows the one-liner .. code-block:: cpp Eps = elem.SymGradN_vector(vector.AElement(disp)); Compute stress and tangent ========================== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 84 The stress and stiffness tensors are now computed for each integration point (completely independently) using the external material model. .. note:: "Sig" and "C" are the output variables that were preallocated in the main. Assemble system =============== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 86-92 The stress stored per integration point ("Sig") is now converted to nodal internal force vectors stored per element ("fe"). Using "vector" this 'elemvec' representation is then converted of a 'nodevec' representation in "fint". Likewise, the stiffness tensor stored for integration point ("C") are converted to system matrices stored per element ('elemmat') and finally assembled to the global stiffness matrix. .. warning:: Please note that downsizing ("fe" :math:`\rightarrow` "fint" and "Ke" :math:`\rightarrow` "K") can be done in two ways, and that "assemble..." is the right function here as it adds entries that occur more than once. In contrast "as..." would not result in what we want here. .. note:: Once more, "fe", "fint", and "Ke" are output variables. Less efficient, but shorter, is: .. code-block:: cpp // internal force fint = vector.AssembleNode(elem.Int_gradN_dot_tensor2_dV(Sig)); // stiffness matrix K.assemble(elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C)); Solve ===== -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 94-104 We now prescribe the displacement of the Prescribed degrees-of-freedom directly in the nodal displacements "disp" and compute the residual force. This is follows by partitioning and solving, all done internally in the "MatrixPartitioned" class. Post-process ============ Strain and stress ----------------- -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 110-112 The strain and stress per integration point are recomputed for post-processing. Residual force -------------- -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 114-125 We convince ourselves that the solution is indeed in mechanical equilibrium. Store & plot ------------ -.. literalinclude:: statics/FixedDisplacements_LinearElastic/example/main.cpp +.. literalinclude:: statics/FixedDisplacements_LinearElastic/example.cpp :language: cpp :lines: 127-136 -Finally we store some fields for plotting using :download:`plot.py `. +Finally we store some fields for plotting using :download:`plot.py `. Manual partitioning =================== To verify how partitioning and solving is done internally using the "MatrixPartitioned" class, the same example is provided where partitioning is done manually: -| :download:`main.cpp ` -| :download:`CMakeLists.txt ` -| :download:`plot.py ` -| :download:`main.py ` +| :download:`manual_partition.cpp ` diff --git a/include/GooseFEM/.gitignore b/include/GooseFEM/.gitignore deleted file mode 100644 index ab787b0..0000000 --- a/include/GooseFEM/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -test* -*_old.h -todo* diff --git a/include/GooseFEM/config.h b/include/GooseFEM/config.h index 2f848d6..4d64fb4 100644 --- a/include/GooseFEM/config.h +++ b/include/GooseFEM/config.h @@ -1,93 +1,112 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_CONFIG_H #define GOOSEFEM_CONFIG_H // ------------------------------------------------------------------------------------------------- +#ifndef GOOSEFEM_DEBUG + #ifndef NDEBUG + #define NDEBUG + #endif +#else + #ifndef GOOSEFEM_ENABLE_ASSERT + #define GOOSEFEM_ENABLE_ASSERT + #endif + #ifndef XTENSOR_ENABLE_ASSERT + #define XTENSOR_ENABLE_ASSERT + #endif +#endif + +// ------------------------------------------------------------------------------------------------- + #define _USE_MATH_DEFINES // to use "M_PI" from "math.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace xt::placeholders; // ------------------------------------------------------------------------------------------------- // dummy operation that can be use to suppress the "unused parameter" warnings #define UNUSED(p) ( (void)(p) ) // ------------------------------------------------------------------------------------------------- #ifndef NDEBUG #ifndef GOOSEFEM_ENABLE_ASSERT #define GOOSEFEM_ENABLE_ASSERT #endif #endif #ifdef GOOSEFEM_ENABLE_ASSERT #define GOOSEFEM_ASSERT(expr) GOOSEFEM_ASSERT_IMPL(expr, __FILE__, __LINE__) -#define GOOSEFEM_ASSERT_IMPL(expr, file, line) \ - if (!(expr)) \ - { \ - throw std::runtime_error(std::string(file) + ':' + std::to_string(line) + ": assertion failed (" #expr ") \n\t"); \ +#define GOOSEFEM_ASSERT_IMPL(expr, file, line) \ + if (!(expr)) \ + { \ + throw std::runtime_error( \ + std::string(file) + ':' + std::to_string(line) + \ + ": assertion failed (" #expr ") \n\t"); \ } #else #define GOOSEFEM_ASSERT(expr) #endif // ------------------------------------------------------------------------------------------------- #define GOOSEFEM_CHECK(expr) GOOSEFEM_CHECK_IMPL(expr, __FILE__, __LINE__) -#define GOOSEFEM_CHECK_IMPL(expr, file, line) \ - if (!(expr)) \ - { \ - throw std::runtime_error(std::string(file) + ':' + std::to_string(line) + ": assertion failed (" #expr ") \n\t"); \ +#define GOOSEFEM_CHECK_IMPL(expr, file, line) \ + if (!(expr)) \ + { \ + throw std::runtime_error( \ + std::string(file) + ':' + std::to_string(line) + \ + ": assertion failed (" #expr ") \n\t"); \ } // ------------------------------------------------------------------------------------------------- -#define GOOSEFEM_WORLD_VERSION 0 -#define GOOSEFEM_MAJOR_VERSION 2 -#define GOOSEFEM_MINOR_VERSION 3 +#define GOOSEFEM_VERSION_MAJOR 0 +#define GOOSEFEM_VERSION_MINOR 3 +#define GOOSEFEM_VERSION_PATCH 0 #define GOOSEFEM_VERSION_AT_LEAST(x,y,z) \ - (GOOSEFEM_WORLD_VERSION>x || (GOOSEFEM_WORLD_VERSION>=x && \ - (GOOSEFEM_MAJOR_VERSION>y || (GOOSEFEM_MAJOR_VERSION>=y && \ - GOOSEFEM_MINOR_VERSION>=z)))) + (GOOSEFEM_VERSION_MAJOR > x || (GOOSEFEM_VERSION_MAJOR >= x && \ + (GOOSEFEM_VERSION_MINOR > y || (GOOSEFEM_VERSION_MINOR >= y && \ + GOOSEFEM_VERSION_PATCH >= z)))) #define GOOSEFEM_VERSION(x,y,z) \ - (GOOSEFEM_WORLD_VERSION==x && \ - GOOSEFEM_MAJOR_VERSION==y && \ - GOOSEFEM_MINOR_VERSION==z) + (GOOSEFEM_VERSION_MAJOR == x && \ + GOOSEFEM_VERSION_MINOR == y && \ + GOOSEFEM_VERSION_PATCH == z) // ------------------------------------------------------------------------------------------------- #endif diff --git a/setup.py b/setup.py index f251676..d66173c 100644 --- a/setup.py +++ b/setup.py @@ -1,51 +1,57 @@ desc = ''' GooseFEM is a C++ module, wrapped in Python, that provides several predefined finite element meshes. The original C++ module also includes element definitions and several standard finite element simulations. ''' from setuptools import setup, Extension -import sys, re -import setuptools +import re +import os import pybind11 import pyxtensor header = open('include/GooseFEM/config.h','r').read() -world = re.split(r'(.*)(\#define GOOSEFEM_WORLD_VERSION\ )([0-9]+)(.*)',header)[3] -major = re.split(r'(.*)(\#define GOOSEFEM_MAJOR_VERSION\ )([0-9]+)(.*)',header)[3] -minor = re.split(r'(.*)(\#define GOOSEFEM_MINOR_VERSION\ )([0-9]+)(.*)',header)[3] +major = re.split(r'(.*)(\#define GOOSEFEM_VERSION_MAJOR\ )([0-9]+)(.*)', header)[3] +minor = re.split(r'(.*)(\#define GOOSEFEM_VERSION_MINOR\ )([0-9]+)(.*)', header)[3] +patch = re.split(r'(.*)(\#define GOOSEFEM_VERSION_PATCH\ )([0-9]+)(.*)', header)[3] -__version__ = '.'.join([world,major,minor]) +__version__ = '.'.join([major, minor, patch]) -ext_modules = [ - Extension( +include_dirs = [ + os.path.abspath('include/'), + pyxtensor.find_pyxtensor(), + pyxtensor.find_pybind11(), + pyxtensor.find_xtensor(), + pyxtensor.find_xtl(), + pyxtensor.find_eigen()] + +build = pyxtensor.BuildExt + +xsimd = pyxtensor.find_xsimd() + +if xsimd: + if len(xsimd) > 0: + include_dirs += [xsimd] + build.c_opts['unix'] += ['-march=native', '-DXTENSOR_USE_XSIMD'] + build.c_opts['msvc'] += ['/DXTENSOR_USE_XSIMD'] + +ext_modules = [Extension( 'GooseFEM', ['python/main.cpp'], - include_dirs=[ - pybind11 .get_include(False), - pybind11 .get_include(True ), - pyxtensor.get_include(False), - pyxtensor.get_include(True ), - pyxtensor.find_xtensor(), - pyxtensor.find_xtl(), - pyxtensor.find_eigen(), - ], - language='c++' - ), -] + include_dirs = include_dirs, + language = 'c++')] setup( - name = 'GooseFEM', - description = 'Finite element meshes, quadrature, and assembly tools', - long_description = desc, - version = __version__, - license = 'GPLv3', - author = 'Tom de Geus', - author_email = 'tom@geus.me', - url = 'https://github.com/tdegeus/GooseFEM', - ext_modules = ext_modules, - install_requires = ['pybind11>=2.2.0','pyxtensor>=0.0.1'], - cmdclass = {'build_ext': pyxtensor.BuildExt}, - zip_safe = False, -) + name = 'GooseFEM', + description = 'Finite element meshes, quadrature, and assembly tools', + long_description = desc, + version = __version__, + license = 'GPLv3', + author = 'Tom de Geus', + author_email = 'tom@geus.me', + url = 'https://github.com/tdegeus/GooseFEM', + ext_modules = ext_modules, + install_requires = ['pybind11>=2.2.0', 'pyxtensor>=0.1.1'], + cmdclass = {'build_ext': build}, + zip_safe = False) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8853fdc..852d7c1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,59 +1,29 @@ -cmake_minimum_required(VERSION 2.8.12) +cmake_minimum_required(VERSION 3.0) -project(test) +if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR) + project(GooseFEM-test) + find_package(GooseFEM REQUIRED CONFIG) +endif() option(WARNINGS "Show build warnings" ON) +option(SIMD "Enable xsimd" ON) +option(DEBUG "Enable all assertions" ON) -# basic compiler options -# ---------------------- - -# set optimization level -set(CMAKE_BUILD_TYPE Release) - -# set C++ standard -set(CMAKE_CXX_STANDARD 14) -set(CMAKE_CXX_STANDARD_REQUIRED ON) - -# switch on warnings -if(WARNINGS) - if(MSVC) - if(CMAKE_CXX_FLAGS MATCHES "/W[0-4]") - string(REGEX REPLACE "/W[0-4]" "/W4" CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}") - else() - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4") - endif() - else() - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic -Wno-unknown-pragmas") - endif() -endif() +find_package(Catch2 REQUIRED) + +add_executable(main main.cpp) -# load packages -# ------------- +target_link_libraries(main PRIVATE Catch2::Catch2 GooseFEM) -# add custom paths -if(NOT "$ENV{INCLUDE_PATH}" STREQUAL "") - string(REPLACE ":" ";" INCLUDE_LIST "$ENV{INCLUDE_PATH}") - include_directories(SYSTEM ${INCLUDE_LIST}) +if (SIMD) + target_link_libraries(main PRIVATE xtensor::optimize xtensor::use_xsimd) endif() -# load pkg-config -find_package(PkgConfig) - -# eigen3 -pkg_check_modules(EIGEN3 REQUIRED eigen3) -include_directories(SYSTEM ${EIGEN3_INCLUDE_DIRS}) - -# compile -# ------- - -add_executable(${PROJECT_NAME} - main.cpp - Vector.cpp - VectorPartitioned.cpp - MatrixDiagonal.cpp - ElementQuad4.cpp - ElementHex8.cpp - Iterate.cpp - MeshQuad4.cpp -) +if (WARNINGS) + target_link_libraries(main PRIVATE GooseFEM::compiler_warnings) +endif() + +if (DEBUG) + target_link_libraries(main PRIVATE GooseFEM::debug) +endif()