diff --git a/.appveyor.yml b/.appveyor.yml index 9746a4c..6950942 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -1,56 +1,63 @@ build: false branches: only: - master platform: - x64 image: - Visual Studio 2017 - Visual Studio 2015 environment: matrix: - MINICONDA: C:\myname-conda init: - "ECHO %MINICONDA%" - if "%APPVEYOR_BUILD_WORKER_IMAGE%" == "Visual Studio 2015" set VCVARPATH="C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat" - if "%APPVEYOR_BUILD_WORKER_IMAGE%" == "Visual Studio 2015" set VCARGUMENT=%PLATFORM% - if "%APPVEYOR_BUILD_WORKER_IMAGE%" == "Visual Studio 2017" set VCVARPATH="C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvars64.bat" - echo "%VCVARPATH% %VCARGUMENT%" - "%VCVARPATH% %VCARGUMENT%" - ps: if($env:Platform -eq "x64"){Start-FileDownload 'http://repo.continuum.io/miniconda/Miniconda3-latest-Windows-x86_64.exe' C:\Miniconda.exe; echo "Done"} - ps: if($env:Platform -eq "x86"){Start-FileDownload 'http://repo.continuum.io/miniconda/Miniconda3-latest-Windows-x86.exe' C:\Miniconda.exe; echo "Done"} - cmd: C:\Miniconda.exe /S /D=C:\myname-conda - "set PATH=%MINICONDA%;%MINICONDA%\\Scripts;%MINICONDA%\\Library\\bin;%PATH%" install: # Set environment using Conda - conda config --set always_yes yes --set changeps1 no - conda update -q conda - conda info -a - conda install -c conda-forge mamba - - mamba install -c conda-forge cmake xtensor xsimd eigen python numpy pyxtensor catch2 h5py highfive gmatelastic gmatnonlinearelastic gmatelastoplastic gmatelastoplasticfinitestrainsimo - # - mamba install -c conda-forge cmake xtensor xsimd python numpy pyxtensor catch2 h5py highfive gmatelastic gmatnonlinearelastic gmatelastoplastic gmatelastoplasticfinitestrainsimo python-gmatelastic + - > + mamba install -c conda-forge + cmake + xtensor + xsimd + eigen + python + numpy + pyxtensor + catch2 + h5py + highfive + gmatelastic + gmatnonlinearelastic + gmatelastoplastic + gmatelastoplasticfinitestrainsimo + gmatelastoplasticqpot # Build/install the library - cmake -G "NMake Makefiles" -DCMAKE_INSTALL_PREFIX=%MINICONDA%\\LIBRARY -DCMAKE_BUILD_TYPE=RELEASE -DBUILD_TESTS=ON -DBUILD_EXAMPLES=OFF . - nmake - nmake install - python setup.py build - python setup.py install build_script: # Run tests - - .\test\main - # Run examples - # - .\docs\examples\statics_FixedDisplacements_LinearElastic_example - # - .\docs\examples\statics_FixedDisplacements_LinearElastic_manual_partition - # - .\docs\examples\statics_Periodic_ElastoPlasticFiniteStrainSimo_main - # - .\docs\examples\statics_Periodic_ElastoPlastic_main - # - .\docs\examples\statics_Periodic_LinearElastic_main - # - .\docs\examples\statics_Periodic_NonLinearElastic_main - # - python .\docs\examples\statics\FixedDisplacements_LinearElastic\example.py --no-plot - # - python .\docs\examples\statics\FixedDisplacements_LinearElastic\manual_partition.py --no-plot + - .\test\basic\test-basic + - .\test\gmat\test-gmat diff --git a/.travis.yml b/.travis.yml index ec1f419..8178ecb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,93 +1,111 @@ language: cpp dist: xenial env: matrix: fast_finish: true include: - os: linux addons: apt: sources: - ubuntu-toolchain-r-test packages: - g++-5 env: COMPILER=gcc GCC=5 - os: linux addons: apt: sources: - ubuntu-toolchain-r-test packages: - g++-6 env: COMPILER=gcc GCC=6 - os: linux addons: apt: sources: - ubuntu-toolchain-r-test packages: - g++-7 env: COMPILER=gcc GCC=7 - os: linux addons: apt: sources: - ubuntu-toolchain-r-test - llvm-toolchain-xenial-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: # 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" - hash -r - conda config --set always_yes yes --set changeps1 no - conda update -q conda - conda install -c conda-forge mamba - - mamba install -c conda-forge cmake xtensor xsimd eigen python numpy pyxtensor catch2 h5py highfive gmatelastic gmatnonlinearelastic gmatelastoplastic gmatelastoplasticfinitestrainsimo python-gmatelastic + - > + mamba install -c conda-forge + cmake + xtensor + xsimd + eigen + python + numpy + pyxtensor + catch2 + h5py + highfive + gmatelastic + gmatnonlinearelastic + gmatelastoplastic + gmatelastoplasticfinitestrainsimo + gmatelastoplasticqpot + python-gmatelastic # 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 - - ./test/main - # run examples + # Run tests + - ./test/basic/test-basic + - ./test/gmat/test-gmat + # Run examples - ./docs/examples/statics_FixedDisplacements_LinearElastic_example - ./docs/examples/statics_FixedDisplacements_LinearElastic_manual_partition # - ./docs/examples/statics_Periodic_LinearElastic_main # - ./docs/examples/statics_Periodic_NonLinearElastic_main # - ./docs/examples/statics_Periodic_ElastoPlastic_main # - ./docs/examples/statics_Periodic_ElastoPlasticFiniteStrainSimo_main - python ./docs/examples/statics/FixedDisplacements_LinearElastic/example.py --no-plot - python ./docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py --no-plot diff --git a/CMakeLists.txt b/CMakeLists.txt index ab47dcd..9ab6a39 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,91 +1,92 @@ # # (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM # cmake_minimum_required(VERSION 3.0) # Basic settings # ============== project(GooseFEM) option(BUILD_TESTS "Build tests" OFF) option(BUILD_EXAMPLES "Build examples" OFF) # Version # ======= file(STRINGS "${CMAKE_CURRENT_SOURCE_DIR}/include/GooseFEM/config.h" GooseFEM_version_defines REGEX "#define GOOSEFEM_VERSION_(MAJOR|MINOR|PATCH)") 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() set(GOOSEFEM_VERSION ${GOOSEFEM_VERSION_MAJOR}.${GOOSEFEM_VERSION_MINOR}.${GOOSEFEM_VERSION_PATCH}) message(STATUS "Building GooseFEM v${GOOSEFEM_VERSION}") # Set target # ========== find_package(xtensor REQUIRED) 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") if(BUILD_TESTS) - add_subdirectory(test) + add_subdirectory(test/basic) + add_subdirectory(test/gmat) endif() if(BUILD_EXAMPLES) add_subdirectory(docs/examples) endif() diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt deleted file mode 100644 index c3d3e57..0000000 --- a/test/CMakeLists.txt +++ /dev/null @@ -1,66 +0,0 @@ - -cmake_minimum_required(VERSION 3.0) - -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(ASSERT "Enable assertions" ON) -option(DEBUG "Enable all assertions" ON) - -find_package(Catch2 REQUIRED) -if(!WIN32) - find_package(GMatElastic REQUIRED) -endif() - -if(!WIN32) - add_executable(main - main.cpp - ElementHex8.cpp - ElementQuad4.cpp - Iterate.cpp - Matrix.cpp - MatrixDiagonal.cpp - Mesh.cpp - MeshQuad4.cpp - Vector.cpp - VectorPartitioned.cpp - Example_hybrid-material.cpp) -else() - add_executable(main - main.cpp - ElementHex8.cpp - ElementQuad4.cpp - Iterate.cpp - Matrix.cpp - MatrixDiagonal.cpp - Mesh.cpp - MeshQuad4.cpp - Vector.cpp - VectorPartitioned.cpp) -endif() - -if(!WIN32) - target_link_libraries(main PRIVATE Catch2::Catch2 GooseFEM GMatElastic) -else() - target_link_libraries(main PRIVATE Catch2::Catch2 GooseFEM) -endif() - -if (SIMD) - target_link_libraries(main PRIVATE xtensor::optimize xtensor::use_xsimd) -endif() - -if (WARNINGS) - target_link_libraries(main PRIVATE GooseFEM::compiler_warnings) -endif() - -if (ASSERT) - target_link_libraries(main PRIVATE GooseFEM::assert) -endif() - -if (DEBUG) - target_link_libraries(main PRIVATE GooseFEM::debug) -endif() diff --git a/test/basic/CMakeLists.txt b/test/basic/CMakeLists.txt new file mode 100644 index 0000000..798cb53 --- /dev/null +++ b/test/basic/CMakeLists.txt @@ -0,0 +1,44 @@ + +cmake_minimum_required(VERSION 3.0) + +if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR) + project(GooseFEM-test-basic) + find_package(GooseFEM REQUIRED CONFIG) +endif() + +option(WARNINGS "Show build warnings" ON) +option(SIMD "Enable xsimd" ON) +option(ASSERT "Enable assertions" ON) +option(DEBUG "Enable all assertions" ON) + +find_package(Catch2 REQUIRED) + +add_executable(test-basic + main.cpp + ElementHex8.cpp + ElementQuad4.cpp + Iterate.cpp + Matrix.cpp + MatrixDiagonal.cpp + Mesh.cpp + MeshQuad4.cpp + Vector.cpp + VectorPartitioned.cpp) + +target_link_libraries(test-basic PRIVATE Catch2::Catch2 GooseFEM) + +if (SIMD) + target_link_libraries(test-basic PRIVATE xtensor::optimize xtensor::use_xsimd) +endif() + +if (WARNINGS) + target_link_libraries(test-basic PRIVATE GooseFEM::compiler_warnings) +endif() + +if (ASSERT) + target_link_libraries(test-basic PRIVATE GooseFEM::assert) +endif() + +if (DEBUG) + target_link_libraries(test-basic PRIVATE GooseFEM::debug) +endif() diff --git a/test/ElementHex8.cpp b/test/basic/ElementHex8.cpp similarity index 96% rename from test/ElementHex8.cpp rename to test/basic/ElementHex8.cpp index b3bfa0a..2018efe 100644 --- a/test/ElementHex8.cpp +++ b/test/basic/ElementHex8.cpp @@ -1,177 +1,182 @@ -#include "support.h" +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::ElementHex8", "ElementHex8.h") { using T2 = xt::xtensor_fixed>; SECTION("dV - Gauss") { GooseFEM::Mesh::Hex8::Regular mesh(2, 2, 2); xt::xtensor coor = mesh.coor(); xt::xtensor top = mesh.nodesTop(); xt::xtensor right = mesh.nodesRight(); xt::xtensor back = mesh.nodesBack(); xt::view(coor, xt::keep(back), xt::keep(2)) += 1.; xt::view(coor, xt::keep(top), xt::keep(1)) += 1.; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Hex8::Quadrature quad(vec.AsElement(coor)); xt::xtensor dV = quad.DV(); xt::xtensor dV_tensor = quad.DV(2); REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(4)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(5)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(6)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(7)), 8. / 8.)); for (size_t e = 0; e < mesh.nelem(); ++e) for (size_t q = 0; q < quad.nip(); ++q) REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } SECTION("dV - Nodal") { GooseFEM::Mesh::Hex8::Regular mesh(2, 2, 2); xt::xtensor coor = mesh.coor(); xt::xtensor top = mesh.nodesTop(); xt::xtensor right = mesh.nodesRight(); xt::xtensor back = mesh.nodesBack(); xt::view(coor, xt::keep(back), xt::keep(2)) += 1.; xt::view(coor, xt::keep(top), xt::keep(1)) += 1.; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Hex8::Quadrature quad( vec.AsElement(coor), GooseFEM::Element::Hex8::Nodal::xi(), GooseFEM::Element::Hex8::Nodal::w()); xt::xtensor dV = quad.DV(); xt::xtensor dV_tensor = quad.DV(2); REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(4)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(5)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(6)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(7)), 8. / 8.)); for (size_t e = 0; e < mesh.nelem(); ++e) for (size_t q = 0; q < quad.nip(); ++q) REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } SECTION("int_N_scalar_NT_dV") { GooseFEM::Mesh::Hex8::Regular mesh(3, 3, 3); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Hex8::Quadrature quad( vec.AsElement(mesh.coor()), GooseFEM::Element::Hex8::Nodal::xi(), GooseFEM::Element::Hex8::Nodal::w()); xt::xtensor rho = xt::ones({mesh.nelem(), quad.nip()}); mat.assemble(quad.Int_N_scalar_NT_dV(rho)); xt::xtensor M = mat.AsDiagonal(); REQUIRE(M.size() == vec.ndof()); REQUIRE(xt::allclose(M, 1.)); } SECTION("symGradN_vector") { GooseFEM::Mesh::Hex8::FineLayer mesh(27, 27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Hex8::Quadrature quad(vec.AsElement(mesh.coor())); T2 F = xt::zeros({3, 3}); T2 EPS = xt::zeros({3, 3}); F(0, 1) = 0.1; EPS(0, 1) = 0.05; EPS(1, 0) = 0.05; xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) for (size_t i = 0; i < F.shape()[0]; ++i) for (size_t j = 0; j < F.shape()[1]; ++j) disp(n, i) += F(i, j) * coor(n, j); xt::xtensor eps = quad.SymGradN_vector(vec.AsElement(disp)); xt::xtensor dV = quad.DV(2); auto epsbar = xt::average(eps, dV, {0, 1}); REQUIRE(eps.shape()[0] == mesh.nelem()); REQUIRE(eps.shape()[1] == quad.nip()); REQUIRE(eps.shape()[2] == mesh.ndim()); REQUIRE(eps.shape()[3] == mesh.ndim()); for (size_t e = 0; e < mesh.nelem(); ++e) for (size_t q = 0; q < quad.nip(); ++q) REQUIRE(xt::allclose(xt::view(eps, e, q), EPS)); REQUIRE(xt::allclose(epsbar, EPS)); } SECTION("symGradN_vector, int_gradN_dot_tensor2s_dV") { GooseFEM::Mesh::Hex8::FineLayer mesh(27, 27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Hex8::Quadrature quad(vec.AsElement(mesh.coor())); T2 F = xt::zeros({3, 3}); F(0, 1) = 0.1; xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) for (size_t i = 0; i < F.shape()[0]; ++i) for (size_t j = 0; j < F.shape()[1]; ++j) disp(n, i) += F(i, j) * coor(n, j); xt::xtensor eps = quad.SymGradN_vector(vec.AsElement(disp)); xt::xtensor Fi = vec.AssembleDofs(quad.Int_gradN_dot_tensor2_dV(eps)); REQUIRE(Fi.size() == vec.ndof()); REQUIRE(xt::allclose(Fi, 0.)); } } diff --git a/test/ElementQuad4.cpp b/test/basic/ElementQuad4.cpp similarity index 96% rename from test/ElementQuad4.cpp rename to test/basic/ElementQuad4.cpp index 6f1fef7..01b2780 100644 --- a/test/ElementQuad4.cpp +++ b/test/basic/ElementQuad4.cpp @@ -1,165 +1,170 @@ -#include "support.h" +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::ElementQuad4", "ElementQuad4.h") { using T2 = xt::xtensor_fixed>; SECTION("dV - Gauss") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2); xt::xtensor coor = mesh.coor(); xt::xtensor top = mesh.nodesTopEdge(); xt::xtensor right = mesh.nodesRightEdge(); xt::view(coor, xt::keep(top), xt::keep(1)) += 1.; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(coor)); xt::xtensor dV = quad.DV(); xt::xtensor dV_tensor = quad.DV(2); REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4. / 4.)); for (size_t e = 0; e < mesh.nelem(); ++e) for (size_t q = 0; q < quad.nip(); ++q) REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } SECTION("dV - Nodal") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2); xt::xtensor coor = mesh.coor(); xt::xtensor top = mesh.nodesTopEdge(); xt::xtensor right = mesh.nodesRightEdge(); xt::view(coor, xt::keep(top), xt::keep(1)) += 1.; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Quad4::Quadrature quad( vec.AsElement(coor), GooseFEM::Element::Quad4::Nodal::xi(), GooseFEM::Element::Quad4::Nodal::w()); xt::xtensor dV = quad.DV(); xt::xtensor dV_tensor = quad.DV(2); REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4. / 4.)); for (size_t e = 0; e < mesh.nelem(); ++e) for (size_t q = 0; q < quad.nip(); ++q) REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } SECTION("int_N_scalar_NT_dV") { GooseFEM::Mesh::Quad4::Regular mesh(3, 3); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Quad4::Quadrature quad( vec.AsElement(mesh.coor()), GooseFEM::Element::Quad4::Nodal::xi(), GooseFEM::Element::Quad4::Nodal::w()); xt::xtensor rho = xt::ones({mesh.nelem(), quad.nip()}); mat.assemble(quad.Int_N_scalar_NT_dV(rho)); xt::xtensor M = mat.AsDiagonal(); REQUIRE(M.size() == vec.ndof()); REQUIRE(xt::allclose(M, 1.)); } SECTION("symGradN_vector") { GooseFEM::Mesh::Quad4::FineLayer mesh(27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor())); T2 F = xt::zeros({2, 2}); T2 EPS = xt::zeros({2, 2}); F(0, 1) = 0.1; EPS(0, 1) = 0.05; EPS(1, 0) = 0.05; xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) for (size_t i = 0; i < F.shape()[0]; ++i) for (size_t j = 0; j < F.shape()[1]; ++j) disp(n, i) += F(i, j) * coor(n, j); xt::xtensor eps = quad.SymGradN_vector(vec.AsElement(disp)); xt::xtensor dV = quad.DV(2); auto epsbar = xt::average(eps, dV, {0, 1}); REQUIRE(eps.shape()[0] == mesh.nelem()); REQUIRE(eps.shape()[1] == quad.nip()); REQUIRE(eps.shape()[2] == mesh.ndim()); REQUIRE(eps.shape()[3] == mesh.ndim()); for (size_t e = 0; e < mesh.nelem(); ++e) for (size_t q = 0; q < quad.nip(); ++q) REQUIRE(xt::allclose(xt::view(eps, e, q), EPS)); REQUIRE(xt::allclose(epsbar, EPS)); } SECTION("symGradN_vector, int_gradN_dot_tensor2s_dV") { GooseFEM::Mesh::Quad4::FineLayer mesh(27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor())); T2 F = xt::zeros({2, 2}); F(0, 1) = 0.1; xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) for (size_t i = 0; i < F.shape()[0]; ++i) for (size_t j = 0; j < F.shape()[1]; ++j) disp(n, i) += F(i, j) * coor(n, j); xt::xtensor eps = quad.SymGradN_vector(vec.AsElement(disp)); xt::xtensor Fi = vec.AssembleDofs(quad.Int_gradN_dot_tensor2_dV(eps)); REQUIRE(Fi.size() == vec.ndof()); REQUIRE(xt::allclose(Fi, 0.)); } } diff --git a/test/Iterate.cpp b/test/basic/Iterate.cpp similarity index 75% rename from test/Iterate.cpp rename to test/basic/Iterate.cpp index 3e21aa2..aa57071 100644 --- a/test/Iterate.cpp +++ b/test/basic/Iterate.cpp @@ -1,22 +1,28 @@ -#include "support.h" +#include +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::Iterate", "Iterate.h") { SECTION("StopList") { GooseFEM::Iterate::StopList stop(5); REQUIRE(stop.stop(5.e+0, 1.e-3) == false); REQUIRE(stop.stop(5.e+1, 1.e-3) == false); REQUIRE(stop.stop(5.e-1, 1.e-3) == false); REQUIRE(stop.stop(5.e-2, 1.e-3) == false); REQUIRE(stop.stop(5.e-3, 1.e-3) == false); REQUIRE(stop.stop(5.e-4, 1.e-3) == false); REQUIRE(stop.stop(5.e-4, 1.e-3) == false); REQUIRE(stop.stop(5.e-4, 1.e-3) == false); REQUIRE(stop.stop(5.e-4, 1.e-3) == false); REQUIRE(stop.stop(5.e-4, 1.e-3) == true); } } diff --git a/test/Matrix.cpp b/test/basic/Matrix.cpp similarity index 82% rename from test/Matrix.cpp rename to test/basic/Matrix.cpp index 9eeb74e..d1c8a29 100644 --- a/test/Matrix.cpp +++ b/test/basic/Matrix.cpp @@ -1,34 +1,40 @@ -#include "support.h" +#include +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::Matrix", "Matrix.h") { SECTION("solve") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); size_t nelem = mesh.nelem(); size_t nnode = mesh.nnode(); xt::xtensor a = xt::empty({nelem, nne * ndim, nne * ndim}); xt::xtensor b = xt::random::rand({nnode * ndim}); for (size_t e = 0; e < nelem; ++e) { xt::xtensor ae = xt::random::rand({nne * ndim, nne * ndim}); ae = (ae + xt::transpose(ae)) / 2.0; xt::view(a, e, xt::all(), xt::all()) = ae; } GooseFEM::Matrix A(mesh.conn(), mesh.dofs()); GooseFEM::MatrixSolver<> Solver; A.assemble(a); xt::xtensor C = A.Dot(b); xt::xtensor B = Solver.Solve(A, C); REQUIRE(B.size() == b.size()); REQUIRE(xt::allclose(B, b)); } } diff --git a/test/MatrixDiagonal.cpp b/test/basic/MatrixDiagonal.cpp similarity index 85% rename from test/MatrixDiagonal.cpp rename to test/basic/MatrixDiagonal.cpp index 68302eb..df81af1 100644 --- a/test/MatrixDiagonal.cpp +++ b/test/basic/MatrixDiagonal.cpp @@ -1,39 +1,44 @@ -#include "support.h" +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::MatrixDiagonal", "MatrixDiagonal.h") { SECTION("dot") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2); xt::xtensor a = xt::random::rand({mesh.nnode() * mesh.ndim()}); xt::xtensor b = xt::random::rand({mesh.nnode() * mesh.ndim()}); xt::xtensor c = a * b; GooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs()); A.set(a); xt::xtensor C = A.Dot(b); REQUIRE(C.size() == c.size()); REQUIRE(xt::allclose(C, c)); } SECTION("solve") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2); xt::xtensor a = xt::random::rand({mesh.nnode() * mesh.ndim()}); xt::xtensor b = xt::random::rand({mesh.nnode() * mesh.ndim()}); xt::xtensor c = a * b; GooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs()); A.set(a); xt::xtensor C = A.Dot(b); xt::xtensor B = A.Solve(C); REQUIRE(B.size() == b.size()); REQUIRE(xt::allclose(B, b)); } } diff --git a/test/Mesh.cpp b/test/basic/Mesh.cpp similarity index 66% rename from test/Mesh.cpp rename to test/basic/Mesh.cpp index 0e59706..21cdbed 100644 --- a/test/Mesh.cpp +++ b/test/basic/Mesh.cpp @@ -1,14 +1,20 @@ -#include "support.h" + +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::Mesh", "Mesh.h") { SECTION("edgesize") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2, 10.0); auto s = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn()); auto t = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn(), mesh.getElementType()); REQUIRE(xt::allclose(s, 10.0)); REQUIRE(xt::allclose(t, 10.0)); } } diff --git a/test/MeshQuad4.cpp b/test/basic/MeshQuad4.cpp similarity index 99% rename from test/MeshQuad4.cpp rename to test/basic/MeshQuad4.cpp index a089209..7de0cf1 100644 --- a/test/MeshQuad4.cpp +++ b/test/basic/MeshQuad4.cpp @@ -1,366 +1,372 @@ -#include "support.h" + +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::MeshQuad4", "MeshQuad4.h") { SECTION("Regular") { xt::xtensor coor_ = {{0., 0.}, {1., 0.}, {2., 0.}, {3., 0.}, {4., 0.}, {5., 0.}, {0., 1.}, {1., 1.}, {2., 1.}, {3., 1.}, {4., 1.}, {5., 1.}, {0., 2.}, {1., 2.}, {2., 2.}, {3., 2.}, {4., 2.}, {5., 2.}, {0., 3.}, {1., 3.}, {2., 3.}, {3., 3.}, {4., 3.}, {5., 3.}}; xt::xtensor conn_ = {{0, 1, 7, 6}, {1, 2, 8, 7}, {2, 3, 9, 8}, {3, 4, 10, 9}, {4, 5, 11, 10}, {6, 7, 13, 12}, {7, 8, 14, 13}, {8, 9, 15, 14}, {9, 10, 16, 15}, {10, 11, 17, 16}, {12, 13, 19, 18}, {13, 14, 20, 19}, {14, 15, 21, 20}, {15, 16, 22, 21}, {16, 17, 23, 22}}; size_t nelem_ = 15; size_t nnode_ = 24; size_t nne_ = 4; size_t ndim_ = 2; size_t nnodePeriodic_ = 15; xt::xtensor nodesBottomEdge_ = {0, 1, 2, 3, 4, 5}; xt::xtensor nodesTopEdge_ = {18, 19, 20, 21, 22, 23}; xt::xtensor nodesLeftEdge_ = {0, 6, 12, 18}; xt::xtensor nodesRightEdge_ = {5, 11, 17, 23}; xt::xtensor nodesBottomOpenEdge_ = {1, 2, 3, 4}; xt::xtensor nodesTopOpenEdge_ = {19, 20, 21, 22}; xt::xtensor nodesLeftOpenEdge_ = {6, 12}; xt::xtensor nodesRightOpenEdge_ = {11, 17}; size_t nodesBottomLeftCorner_ = 0; size_t nodesBottomRightCorner_ = 5; size_t nodesTopLeftCorner_ = 18; size_t nodesTopRightCorner_ = 23; size_t nodesLeftBottomCorner_ = 0; size_t nodesLeftTopCorner_ = 18; size_t nodesRightBottomCorner_ = 5; size_t nodesRightTopCorner_ = 23; xt::xtensor dofs_ = {{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {10, 11}, {12, 13}, {14, 15}, {16, 17}, {18, 19}, {20, 21}, {22, 23}, {24, 25}, {26, 27}, {28, 29}, {30, 31}, {32, 33}, {34, 35}, {36, 37}, {38, 39}, {40, 41}, {42, 43}, {44, 45}, {46, 47}}; xt::xtensor nodesPeriodic_ = { {0, 5}, {0, 23}, {0, 18}, {1, 19}, {2, 20}, {3, 21}, {4, 22}, {6, 11}, {12, 17}}; size_t nodesOrigin_ = 0; xt::xtensor dofsPeriodic_ = { {0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}, {10, 11}, {12, 13}, {14, 15}, {16, 17}, {18, 19}, {10, 11}, {20, 21}, {22, 23}, {24, 25}, {26, 27}, {28, 29}, {20, 21}, {0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}}; GooseFEM::Mesh::Quad4::Regular mesh(5, 3); xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); size_t nelem = mesh.nelem(); size_t nnode = mesh.nnode(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); size_t nnodePeriodic = mesh.nnode() - mesh.nodesPeriodic().shape(0); xt::xtensor nodesBottomEdge = mesh.nodesBottomEdge(); xt::xtensor nodesTopEdge = mesh.nodesTopEdge(); xt::xtensor nodesLeftEdge = mesh.nodesLeftEdge(); xt::xtensor nodesRightEdge = mesh.nodesRightEdge(); xt::xtensor nodesBottomOpenEdge = mesh.nodesBottomOpenEdge(); xt::xtensor nodesTopOpenEdge = mesh.nodesTopOpenEdge(); xt::xtensor nodesLeftOpenEdge = mesh.nodesLeftOpenEdge(); xt::xtensor nodesRightOpenEdge = mesh.nodesRightOpenEdge(); size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner(); size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner(); size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner(); size_t nodesTopRightCorner = mesh.nodesTopRightCorner(); size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner(); size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner(); size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner(); size_t nodesRightTopCorner = mesh.nodesRightTopCorner(); xt::xtensor dofs = mesh.dofs(); xt::xtensor nodesPeriodic = mesh.nodesPeriodic(); size_t nodesOrigin = mesh.nodesOrigin(); xt::xtensor dofsPeriodic = mesh.dofsPeriodic(); REQUIRE(nelem == nelem_); REQUIRE(nnode == nnode_); REQUIRE(nne == nne_); REQUIRE(ndim == ndim_); REQUIRE(nnodePeriodic == nnodePeriodic_); REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_); REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_); REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_); REQUIRE(nodesTopRightCorner == nodesTopRightCorner_); REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_); REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_); REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_); REQUIRE(nodesRightTopCorner == nodesRightTopCorner_); REQUIRE(nodesOrigin == nodesOrigin_); REQUIRE(xt::allclose(coor, coor_)); REQUIRE(xt::all(xt::equal(conn, conn_))); REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_))); REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_))); REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_))); REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_))); REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_))); REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_))); } SECTION("FineLayer") { xt::xtensor coor_ = { {0., 0.}, {3., 0.}, {6., 0.}, {9., 0.}, {0., 3.}, {3., 3.}, {6., 3.}, {9., 3.}, {0., 6.}, {3., 6.}, {6., 6.}, {9., 6.}, {1., 7.}, {2., 7.}, {4., 7.}, {5., 7.}, {7., 7.}, {8., 7.}, {0., 8.}, {1., 8.}, {2., 8.}, {3., 8.}, {4., 8.}, {5., 8.}, {6., 8.}, {7., 8.}, {8., 8.}, {9., 8.}, {0., 9.}, {1., 9.}, {2., 9.}, {3., 9.}, {4., 9.}, {5., 9.}, {6., 9.}, {7., 9.}, {8., 9.}, {9., 9.}, {0., 10.}, {1., 10.}, {2., 10.}, {3., 10.}, {4., 10.}, {5., 10.}, {6., 10.}, {7., 10.}, {8., 10.}, {9., 10.}, {0., 11.}, {1., 11.}, {2., 11.}, {3., 11.}, {4., 11.}, {5., 11.}, {6., 11.}, {7., 11.}, {8., 11.}, {9., 11.}, {0., 12.}, {1., 12.}, {2., 12.}, {3., 12.}, {4., 12.}, {5., 12.}, {6., 12.}, {7., 12.}, {8., 12.}, {9., 12.}, {0., 13.}, {1., 13.}, {2., 13.}, {3., 13.}, {4., 13.}, {5., 13.}, {6., 13.}, {7., 13.}, {8., 13.}, {9., 13.}, {1., 14.}, {2., 14.}, {4., 14.}, {5., 14.}, {7., 14.}, {8., 14.}, {0., 15.}, {3., 15.}, {6., 15.}, {9., 15.}, {0., 18.}, {3., 18.}, {6., 18.}, {9., 18.}, {0., 21.}, {3., 21.}, {6., 21.}, {9., 21.}}; xt::xtensor conn_ = { {0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {4, 5, 9, 8}, {5, 6, 10, 9}, {6, 7, 11, 10}, {8, 9, 13, 12}, {9, 21, 20, 13}, {12, 13, 20, 19}, {8, 12, 19, 18}, {9, 10, 15, 14}, {10, 24, 23, 15}, {14, 15, 23, 22}, {9, 14, 22, 21}, {10, 11, 17, 16}, {11, 27, 26, 17}, {16, 17, 26, 25}, {10, 16, 25, 24}, {18, 19, 29, 28}, {19, 20, 30, 29}, {20, 21, 31, 30}, {21, 22, 32, 31}, {22, 23, 33, 32}, {23, 24, 34, 33}, {24, 25, 35, 34}, {25, 26, 36, 35}, {26, 27, 37, 36}, {28, 29, 39, 38}, {29, 30, 40, 39}, {30, 31, 41, 40}, {31, 32, 42, 41}, {32, 33, 43, 42}, {33, 34, 44, 43}, {34, 35, 45, 44}, {35, 36, 46, 45}, {36, 37, 47, 46}, {38, 39, 49, 48}, {39, 40, 50, 49}, {40, 41, 51, 50}, {41, 42, 52, 51}, {42, 43, 53, 52}, {43, 44, 54, 53}, {44, 45, 55, 54}, {45, 46, 56, 55}, {46, 47, 57, 56}, {48, 49, 59, 58}, {49, 50, 60, 59}, {50, 51, 61, 60}, {51, 52, 62, 61}, {52, 53, 63, 62}, {53, 54, 64, 63}, {54, 55, 65, 64}, {55, 56, 66, 65}, {56, 57, 67, 66}, {58, 59, 69, 68}, {59, 60, 70, 69}, {60, 61, 71, 70}, {61, 62, 72, 71}, {62, 63, 73, 72}, {63, 64, 74, 73}, {64, 65, 75, 74}, {65, 66, 76, 75}, {66, 67, 77, 76}, {68, 69, 78, 84}, {69, 70, 79, 78}, {70, 71, 85, 79}, {78, 79, 85, 84}, {71, 72, 80, 85}, {72, 73, 81, 80}, {73, 74, 86, 81}, {80, 81, 86, 85}, {74, 75, 82, 86}, {75, 76, 83, 82}, {76, 77, 87, 83}, {82, 83, 87, 86}, {84, 85, 89, 88}, {85, 86, 90, 89}, {86, 87, 91, 90}, {88, 89, 93, 92}, {89, 90, 94, 93}, {90, 91, 95, 94}}; size_t nelem_ = 81; size_t nnode_ = 96; size_t nne_ = 4; size_t ndim_ = 2; size_t shape_x_ = 9; size_t shape_y_ = 21; xt::xtensor elementsMiddleLayer_ = {36, 37, 38, 39, 40, 41, 42, 43, 44}; xt::xtensor nodesBottomEdge_ = {0, 1, 2, 3}; xt::xtensor nodesTopEdge_ = {92, 93, 94, 95}; xt::xtensor nodesLeftEdge_ = {0, 4, 8, 18, 28, 38, 48, 58, 68, 84, 88, 92}; xt::xtensor nodesRightEdge_ = {3, 7, 11, 27, 37, 47, 57, 67, 77, 87, 91, 95}; xt::xtensor nodesBottomOpenEdge_ = {1, 2}; xt::xtensor nodesTopOpenEdge_ = {93, 94}; xt::xtensor nodesLeftOpenEdge_ = {4, 8, 18, 28, 38, 48, 58, 68, 84, 88}; xt::xtensor nodesRightOpenEdge_ = {7, 11, 27, 37, 47, 57, 67, 77, 87, 91}; size_t nodesBottomLeftCorner_ = 0; size_t nodesBottomRightCorner_ = 3; size_t nodesTopLeftCorner_ = 92; size_t nodesTopRightCorner_ = 95; size_t nodesLeftBottomCorner_ = 0; size_t nodesLeftTopCorner_ = 92; size_t nodesRightBottomCorner_ = 3; size_t nodesRightTopCorner_ = 95; xt::xtensor dofs_ = { {0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {10, 11}, {12, 13}, {14, 15}, {16, 17}, {18, 19}, {20, 21}, {22, 23}, {24, 25}, {26, 27}, {28, 29}, {30, 31}, {32, 33}, {34, 35}, {36, 37}, {38, 39}, {40, 41}, {42, 43}, {44, 45}, {46, 47}, {48, 49}, {50, 51}, {52, 53}, {54, 55}, {56, 57}, {58, 59}, {60, 61}, {62, 63}, {64, 65}, {66, 67}, {68, 69}, {70, 71}, {72, 73}, {74, 75}, {76, 77}, {78, 79}, {80, 81}, {82, 83}, {84, 85}, {86, 87}, {88, 89}, {90, 91}, {92, 93}, {94, 95}, {96, 97}, {98, 99}, {100, 101}, {102, 103}, {104, 105}, {106, 107}, {108, 109}, {110, 111}, {112, 113}, {114, 115}, {116, 117}, {118, 119}, {120, 121}, {122, 123}, {124, 125}, {126, 127}, {128, 129}, {130, 131}, {132, 133}, {134, 135}, {136, 137}, {138, 139}, {140, 141}, {142, 143}, {144, 145}, {146, 147}, {148, 149}, {150, 151}, {152, 153}, {154, 155}, {156, 157}, {158, 159}, {160, 161}, {162, 163}, {164, 165}, {166, 167}, {168, 169}, {170, 171}, {172, 173}, {174, 175}, {176, 177}, {178, 179}, {180, 181}, {182, 183}, {184, 185}, {186, 187}, {188, 189}, {190, 191}}; xt::xtensor nodesPeriodic_ = {{0, 3}, {0, 95}, {0, 92}, {1, 93}, {2, 94}, {4, 7}, {8, 11}, {18, 27}, {28, 37}, {38, 47}, {48, 57}, {58, 67}, {68, 77}, {84, 87}, {88, 91}}; size_t nodesOrigin_ = 0; xt::xtensor dofsPeriodic_ = { {0, 1}, {2, 3}, {4, 5}, {0, 1}, {6, 7}, {8, 9}, {10, 11}, {6, 7}, {12, 13}, {14, 15}, {16, 17}, {12, 13}, {18, 19}, {20, 21}, {22, 23}, {24, 25}, {26, 27}, {28, 29}, {30, 31}, {32, 33}, {34, 35}, {36, 37}, {38, 39}, {40, 41}, {42, 43}, {44, 45}, {46, 47}, {30, 31}, {48, 49}, {50, 51}, {52, 53}, {54, 55}, {56, 57}, {58, 59}, {60, 61}, {62, 63}, {64, 65}, {48, 49}, {66, 67}, {68, 69}, {70, 71}, {72, 73}, {74, 75}, {76, 77}, {78, 79}, {80, 81}, {82, 83}, {66, 67}, {84, 85}, {86, 87}, {88, 89}, {90, 91}, {92, 93}, {94, 95}, {96, 97}, {98, 99}, {100, 101}, {84, 85}, {102, 103}, {104, 105}, {106, 107}, {108, 109}, {110, 111}, {112, 113}, {114, 115}, {116, 117}, {118, 119}, {102, 103}, {120, 121}, {122, 123}, {124, 125}, {126, 127}, {128, 129}, {130, 131}, {132, 133}, {134, 135}, {136, 137}, {120, 121}, {138, 139}, {140, 141}, {142, 143}, {144, 145}, {146, 147}, {148, 149}, {150, 151}, {152, 153}, {154, 155}, {150, 151}, {156, 157}, {158, 159}, {160, 161}, {156, 157}, {0, 1}, {2, 3}, {4, 5}, {0, 1}}; GooseFEM::Mesh::Quad4::FineLayer mesh(9, 18); xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); size_t nelem = mesh.nelem(); size_t nnode = mesh.nnode(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); size_t shape_x = mesh.nelx(); size_t shape_y = mesh.nely(); xt::xtensor elementsMiddleLayer = mesh.elementsMiddleLayer(); xt::xtensor nodesBottomEdge = mesh.nodesBottomEdge(); xt::xtensor nodesTopEdge = mesh.nodesTopEdge(); xt::xtensor nodesLeftEdge = mesh.nodesLeftEdge(); xt::xtensor nodesRightEdge = mesh.nodesRightEdge(); xt::xtensor nodesBottomOpenEdge = mesh.nodesBottomOpenEdge(); xt::xtensor nodesTopOpenEdge = mesh.nodesTopOpenEdge(); xt::xtensor nodesLeftOpenEdge = mesh.nodesLeftOpenEdge(); xt::xtensor nodesRightOpenEdge = mesh.nodesRightOpenEdge(); size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner(); size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner(); size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner(); size_t nodesTopRightCorner = mesh.nodesTopRightCorner(); size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner(); size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner(); size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner(); size_t nodesRightTopCorner = mesh.nodesRightTopCorner(); xt::xtensor dofs = mesh.dofs(); xt::xtensor nodesPeriodic = mesh.nodesPeriodic(); size_t nodesOrigin = mesh.nodesOrigin(); xt::xtensor dofsPeriodic = mesh.dofsPeriodic(); REQUIRE(nelem == nelem_); REQUIRE(nnode == nnode_); REQUIRE(nne == nne_); REQUIRE(ndim == ndim_); REQUIRE(shape_x == shape_x_); REQUIRE(shape_y == shape_y_); REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_); REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_); REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_); REQUIRE(nodesTopRightCorner == nodesTopRightCorner_); REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_); REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_); REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_); REQUIRE(nodesRightTopCorner == nodesRightTopCorner_); REQUIRE(nodesOrigin == nodesOrigin_); REQUIRE(xt::allclose(coor, coor_)); REQUIRE(xt::all(xt::equal(elementsMiddleLayer, elementsMiddleLayer_))); REQUIRE(xt::all(xt::equal(conn, conn_))); REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_))); REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_))); REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_))); REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_))); REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_))); REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_))); REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_))); } SECTION("Map::RefineRegular") { GooseFEM::Mesh::Quad4::Regular mesh(5, 4); GooseFEM::Mesh::Quad4::Map::RefineRegular refine(mesh, 5, 3); xt::xtensor a = xt::random::rand({mesh.nelem()}); auto a_ = refine.mapToCoarse(refine.mapToFine(a)); REQUIRE(xt::allclose(a, xt::mean(a_, {1}))); xt::xtensor b = xt::random::rand(std::array{mesh.nelem(), 4ul}); auto b_ = refine.mapToCoarse(refine.mapToFine(b)); REQUIRE(xt::allclose(xt::mean(b, {1}), xt::mean(b_, {1}))); xt::xtensor c = xt::random::rand(std::array{mesh.nelem(), 4ul, 3ul, 3ul}); auto c_ = refine.mapToCoarse(refine.mapToFine(c)); REQUIRE(xt::allclose(xt::mean(c, {1}), xt::mean(c_, {1}))); } SECTION("Map::FineLayer2Regular") { GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5); GooseFEM::Mesh::Quad4::Map::FineLayer2Regular map(mesh); xt::xtensor a = xt::random::rand({mesh.nelem()}); auto a_ = map.mapToRegular(a); REQUIRE(xt::allclose(a, a_)); xt::xtensor b = xt::random::rand(std::array{mesh.nelem(), 4ul}); auto b_ = map.mapToRegular(b); REQUIRE(xt::allclose(b, b_)); xt::xtensor c = xt::random::rand(std::array{mesh.nelem(), 4ul, 3ul, 3ul}); auto c_ = map.mapToRegular(c); REQUIRE(xt::allclose(c, c_)); } } diff --git a/test/Vector.cpp b/test/basic/Vector.cpp similarity index 96% rename from test/Vector.cpp rename to test/basic/Vector.cpp index d40328b..9d5d7e7 100644 --- a/test/Vector.cpp +++ b/test/basic/Vector.cpp @@ -1,198 +1,203 @@ -#include "support.h" +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::Vector", "Vector.h") { SECTION("asDofs - nodevec") { // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); // velocity field // - allocate xt::xtensor v = xt::empty({mesh.nnode(), std::size_t(2)}); // - set periodic v(0, 0) = 1.0; v(0, 1) = 0.0; v(1, 0) = 1.0; v(1, 1) = 0.0; v(2, 0) = 1.0; v(2, 1) = 0.0; v(3, 0) = 1.5; v(3, 1) = 0.0; v(4, 0) = 1.5; v(4, 1) = 0.0; v(5, 0) = 1.5; v(5, 1) = 0.0; v(6, 0) = 1.0; v(6, 1) = 0.0; v(7, 0) = 1.0; v(7, 1) = 0.0; v(8, 0) = 1.0; v(8, 1) = 0.0; // convert to DOFs xt::xtensor V = vector.AsDofs(v); // check // - size REQUIRE(V.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); // - individual entries ISCLOSE(V(0), v(0, 0)); ISCLOSE(V(1), v(0, 1)); ISCLOSE(V(2), v(1, 0)); ISCLOSE(V(3), v(1, 1)); ISCLOSE(V(4), v(3, 0)); ISCLOSE(V(5), v(3, 1)); ISCLOSE(V(6), v(4, 0)); ISCLOSE(V(7), v(4, 1)); } SECTION("asDofs - elemvec") { // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); // velocity field // - allocate xt::xtensor v = xt::empty({mesh.nnode(), std::size_t(2)}); // - set periodic v(0, 0) = 1.0; v(0, 1) = 0.0; v(1, 0) = 1.0; v(1, 1) = 0.0; v(2, 0) = 1.0; v(2, 1) = 0.0; v(3, 0) = 1.5; v(3, 1) = 0.0; v(4, 0) = 1.5; v(4, 1) = 0.0; v(5, 0) = 1.5; v(5, 1) = 0.0; v(6, 0) = 1.0; v(6, 1) = 0.0; v(7, 0) = 1.0; v(7, 1) = 0.0; v(8, 0) = 1.0; v(8, 1) = 0.0; // convert to DOFs - element - DOFs xt::xtensor V = vector.AsDofs(vector.AsElement(vector.AsDofs(v))); // check // - size REQUIRE(V.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); // - individual entries ISCLOSE(V(0), v(0, 0)); ISCLOSE(V(1), v(0, 1)); ISCLOSE(V(2), v(1, 0)); ISCLOSE(V(3), v(1, 1)); ISCLOSE(V(4), v(3, 0)); ISCLOSE(V(5), v(3, 1)); ISCLOSE(V(6), v(4, 0)); ISCLOSE(V(7), v(4, 1)); } SECTION("asDofs - assembleDofs") { // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); // force field // - allocate xt::xtensor f = xt::empty({mesh.nnode(), std::size_t(2)}); // - set periodic f(0, 0) = -1.0; f(0, 1) = -1.0; f(1, 0) = 0.0; f(1, 1) = -1.0; f(2, 0) = 1.0; f(2, 1) = -1.0; f(3, 0) = -1.0; f(3, 1) = 0.0; f(4, 0) = 0.0; f(4, 1) = 0.0; f(5, 0) = 1.0; f(5, 1) = 0.0; f(6, 0) = -1.0; f(6, 1) = 1.0; f(7, 0) = 0.0; f(7, 1) = 1.0; f(8, 0) = 1.0; f(8, 1) = 1.0; // assemble as DOFs xt::xtensor F = vector.AssembleDofs(f); // check // - size REQUIRE(F.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); // - 'analytical' result ISCLOSE(F(0), 0); ISCLOSE(F(1), 0); ISCLOSE(F(2), 0); ISCLOSE(F(3), 0); ISCLOSE(F(4), 0); ISCLOSE(F(5), 0); ISCLOSE(F(6), 0); ISCLOSE(F(7), 0); } SECTION("asDofs - assembleNode") { // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); // force field // - allocate xt::xtensor f = xt::empty({mesh.nnode(), std::size_t(2)}); // - set periodic f(0, 0) = -1.0; f(0, 1) = -1.0; f(1, 0) = 0.0; f(1, 1) = -1.0; f(2, 0) = 1.0; f(2, 1) = -1.0; f(3, 0) = -1.0; f(3, 1) = 0.0; f(4, 0) = 0.0; f(4, 1) = 0.0; f(5, 0) = 1.0; f(5, 1) = 0.0; f(6, 0) = -1.0; f(6, 1) = 1.0; f(7, 0) = 0.0; f(7, 1) = 1.0; f(8, 0) = 1.0; f(8, 1) = 1.0; // convert to element, assemble as DOFs xt::xtensor F = vector.AssembleDofs(vector.AsElement(f)); // check // - size REQUIRE(F.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); // - 'analytical' result ISCLOSE(F(0), 0); ISCLOSE(F(1), 0); ISCLOSE(F(2), 0); ISCLOSE(F(3), 0); ISCLOSE(F(4), 0); ISCLOSE(F(5), 0); ISCLOSE(F(6), 0); ISCLOSE(F(7), 0); } } diff --git a/test/VectorPartitioned.cpp b/test/basic/VectorPartitioned.cpp similarity index 94% rename from test/VectorPartitioned.cpp rename to test/basic/VectorPartitioned.cpp index 7cad8b5..a6faa47 100644 --- a/test/VectorPartitioned.cpp +++ b/test/basic/VectorPartitioned.cpp @@ -1,126 +1,131 @@ -#include "support.h" +#include +#include +#include +#include + +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::VectorPartitioned", "VectorPartitioned.h") { SECTION("asDofs") { // mesh GooseFEM::Mesh::Quad4::Regular mesh(9, 9); // prescribed DOFs xt::xtensor dofs = mesh.dofs(); xt::xtensor nodesLeft = mesh.nodesLeftOpenEdge(); xt::xtensor nodesRight = mesh.nodesRightOpenEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBottom = mesh.nodesBottomEdge(); for (size_t i = 0; i < nodesLeft.size(); ++i) for (size_t j = 0; j < mesh.ndim(); ++j) dofs(nodesRight(i), j) = dofs(nodesLeft(i), j); xt::xtensor iip = xt::empty({4 * nodesBottom.size()}); size_t i = 0; for (auto& n : nodesTop) { iip(i) = dofs(n, 0); ++i; } for (auto& n : nodesTop) { iip(i) = dofs(n, 1); ++i; } for (auto& n : nodesBottom) { iip(i) = dofs(n, 0); ++i; } for (auto& n : nodesBottom) { iip(i) = dofs(n, 1); ++i; } // random displacement xt::xtensor u = xt::random::randn({mesh.nnode(), mesh.ndim()}); for (size_t i = 0; i < nodesLeft.size(); ++i) for (size_t j = 0; j < mesh.ndim(); ++j) u(nodesRight(i), j) = u(nodesLeft(i), j); // vector-definition GooseFEM::VectorPartitioned vector(mesh.conn(), dofs, iip); // convert xt::xtensor u_u = vector.AsDofs_u(u); xt::xtensor u_p = vector.AsDofs_p(u); REQUIRE(xt::allclose(u, vector.AsNode(u_u, u_p))); } SECTION("copy_u, copy_p") { // mesh GooseFEM::Mesh::Quad4::Regular mesh(9, 9); // prescribed DOFs xt::xtensor dofs = mesh.dofs(); xt::xtensor nodesLeft = mesh.nodesLeftOpenEdge(); xt::xtensor nodesRight = mesh.nodesRightOpenEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBottom = mesh.nodesBottomEdge(); for (size_t i = 0; i < nodesLeft.size(); ++i) for (size_t j = 0; j < mesh.ndim(); ++j) dofs(nodesRight(i), j) = dofs(nodesLeft(i), j); xt::xtensor iip = xt::empty({4 * nodesBottom.size()}); size_t i = 0; for (auto& n : nodesTop) { iip(i) = dofs(n, 0); ++i; } for (auto& n : nodesTop) { iip(i) = dofs(n, 1); ++i; } for (auto& n : nodesBottom) { iip(i) = dofs(n, 0); ++i; } for (auto& n : nodesBottom) { iip(i) = dofs(n, 1); ++i; } // random displacement xt::xtensor u = xt::random::randn({mesh.nnode(), mesh.ndim()}); for (size_t i = 0; i < nodesLeft.size(); ++i) for (size_t j = 0; j < mesh.ndim(); ++j) u(nodesRight(i), j) = u(nodesLeft(i), j); // vector-definition GooseFEM::VectorPartitioned vector(mesh.conn(), dofs, iip); // convert xt::xtensor v = xt::empty({mesh.nnode(), mesh.ndim()}); vector.copy_u(u, v); vector.copy_p(u, v); REQUIRE(xt::allclose(u, v)); } } diff --git a/test/main.cpp b/test/basic/main.cpp similarity index 100% copy from test/main.cpp copy to test/basic/main.cpp diff --git a/test/gmat/CMakeLists.txt b/test/gmat/CMakeLists.txt new file mode 100644 index 0000000..c7b460c --- /dev/null +++ b/test/gmat/CMakeLists.txt @@ -0,0 +1,45 @@ + +cmake_minimum_required(VERSION 3.0) + +if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR) + project(GooseFEM-test-gmat) + find_package(GooseFEM REQUIRED CONFIG) +endif() + +option(WARNINGS "Show build warnings" ON) +option(SIMD "Enable xsimd" ON) +option(ASSERT "Enable assertions" ON) +option(DEBUG "Enable all assertions" ON) + +find_package(Catch2 REQUIRED) + +if(NOT WIN32) + find_package(GMatElastic REQUIRED) + + add_executable(test-gmat + main.cpp + hybrid-elastic.cpp) + + target_link_libraries(test-gmat PRIVATE Catch2::Catch2 GooseFEM GMatElastic) +else() + add_executable(test-gmat + main.cpp) + + target_link_libraries(test-gmat PRIVATE Catch2::Catch2 GooseFEM) +endif() + +if (SIMD) + target_link_libraries(test-gmat PRIVATE xtensor::optimize xtensor::use_xsimd) +endif() + +if (WARNINGS) + target_link_libraries(test-gmat PRIVATE GooseFEM::compiler_warnings) +endif() + +if (ASSERT) + target_link_libraries(test-gmat PRIVATE GooseFEM::assert) +endif() + +if (DEBUG) + target_link_libraries(test-gmat PRIVATE GooseFEM::debug) +endif() diff --git a/test/Example_hybrid-material.cpp b/test/gmat/hybrid-elastic.cpp similarity index 98% rename from test/Example_hybrid-material.cpp rename to test/gmat/hybrid-elastic.cpp index f138a26..2b3f11d 100644 --- a/test/Example_hybrid-material.cpp +++ b/test/gmat/hybrid-elastic.cpp @@ -1,256 +1,261 @@ -#include "support.h" +#include #include +#include +#include +#include #include -TEST_CASE("Example_hybrid-material", "GooseFEM.h") +#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.0e-12)); + +TEST_CASE("hybrid-elastic", "GooseFEM.h") { SECTION("Vector/Matrix - GMatElastic") { size_t N = 5; GooseFEM::Mesh::Quad4::Regular mesh(N, N); size_t nelem = mesh.nelem(); xt::xtensor elem_a = xt::arange(N); xt::xtensor elem_b = xt::arange(N, nelem); size_t nelem_a = elem_a.size(); size_t nelem_b = elem_b.size(); xt::xtensor conn = mesh.conn(); xt::xtensor conn_a = xt::view(conn, xt::keep(elem_a), xt::all()); xt::xtensor conn_b = xt::view(conn, xt::keep(elem_b), xt::all()); GooseFEM::Vector vector(conn, mesh.dofs()); GooseFEM::Vector vector_a(conn_a, mesh.dofs()); GooseFEM::Vector vector_b(conn_b, mesh.dofs()); GooseFEM::Matrix K(conn, mesh.dofs()); GooseFEM::Matrix K_a(conn_a, mesh.dofs()); GooseFEM::Matrix K_b(conn_b, mesh.dofs()); xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::random::rand(coor.shape()); GooseFEM::Element::Quad4::QuadraturePlanar quad(vector.AsElement(coor)); GooseFEM::Element::Quad4::QuadraturePlanar quad_a(vector_a.AsElement(coor)); GooseFEM::Element::Quad4::QuadraturePlanar quad_b(vector_b.AsElement(coor)); size_t nip = quad.nip(); GMatElastic::Cartesian3d::Matrix mat(nelem, nip); GMatElastic::Cartesian3d::Matrix mat_a(nelem_a, nip, 3.0, 4.0); GMatElastic::Cartesian3d::Matrix mat_b(nelem_b, nip, 5.0, 6.0); xt::xtensor I = xt::empty({nelem, nip}); I.fill(0); xt::view(I, xt::keep(elem_a), xt::all()) = 1; mat.setElastic(I, 3.0, 4.0); I.fill(0); xt::view(I, xt::keep(elem_b), xt::all()) = 1; mat.setElastic(I, 5.0, 6.0); mat.check(); mat_a.check(); mat_b.check(); xt::xtensor Eps = quad.SymGradN_vector(vector.AsElement(disp)); xt::xtensor Eps_a = quad_a.SymGradN_vector(vector_a.AsElement(disp)); xt::xtensor Eps_b = quad_b.SymGradN_vector(vector_b.AsElement(disp)); xt::xtensor Sig = mat.Stress(Eps); xt::xtensor Sig_a = mat_a.Stress(Eps_a); xt::xtensor Sig_b = mat_b.Stress(Eps_b); xt::xtensor C; xt::xtensor C_a; xt::xtensor C_b; std::tie(Sig, C) = mat.Tangent(Eps); std::tie(Sig_a, C_a) = mat_a.Tangent(Eps_a); std::tie(Sig_b, C_b) = mat_b.Tangent(Eps_b); K.assemble(quad.Int_gradN_dot_tensor4_dot_gradNT_dV(C)); K_a.assemble(quad_a.Int_gradN_dot_tensor4_dot_gradNT_dV(C_a)); K_b.assemble(quad_b.Int_gradN_dot_tensor4_dot_gradNT_dV(C_b)); auto f = vector.AssembleNode(quad.Int_gradN_dot_tensor2_dV(Sig)); auto f_a = vector_a.AssembleNode(quad_a.Int_gradN_dot_tensor2_dV(Sig_a)); auto f_b = vector_b.AssembleNode(quad_b.Int_gradN_dot_tensor2_dV(Sig_b)); REQUIRE(xt::allclose(f, f_a + f_b)); REQUIRE(xt::allclose(f, K.Dot(disp))); REQUIRE(xt::allclose(f_a, K_a.Dot(disp))); REQUIRE(xt::allclose(f_b, K_b.Dot(disp))); auto fd = vector.AssembleDofs(quad.Int_gradN_dot_tensor2_dV(Sig)); auto fd_a = vector_a.AssembleDofs(quad_a.Int_gradN_dot_tensor2_dV(Sig_a)); auto fd_b = vector_b.AssembleDofs(quad_b.Int_gradN_dot_tensor2_dV(Sig_b)); REQUIRE(xt::allclose(fd, fd_a + fd_b)); REQUIRE(xt::allclose(fd, K.Dot(vector.AsDofs(disp)))); REQUIRE(xt::allclose(fd_a, K_a.Dot(vector.AsDofs(disp)))); REQUIRE(xt::allclose(fd_b, K_b.Dot(vector.AsDofs(disp)))); } SECTION("Vector/Matrix - GMatElastic - Periodic") { size_t N = 5; GooseFEM::Mesh::Quad4::Regular mesh(N, N); size_t nelem = mesh.nelem(); xt::xtensor elem_a = xt::arange(N); xt::xtensor elem_b = xt::arange(N, nelem); size_t nelem_a = elem_a.size(); size_t nelem_b = elem_b.size(); xt::xtensor conn = mesh.conn(); xt::xtensor conn_a = xt::view(conn, xt::keep(elem_a), xt::all()); xt::xtensor conn_b = xt::view(conn, xt::keep(elem_b), xt::all()); GooseFEM::Vector vector(conn, mesh.dofsPeriodic()); GooseFEM::Vector vector_a(conn_a, mesh.dofsPeriodic()); GooseFEM::Vector vector_b(conn_b, mesh.dofsPeriodic()); GooseFEM::Matrix K(conn, mesh.dofsPeriodic()); GooseFEM::Matrix K_a(conn_a, mesh.dofsPeriodic()); GooseFEM::Matrix K_b(conn_b, mesh.dofsPeriodic()); xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::random::rand(coor.shape()); disp = vector.AsNode(vector.AsDofs(disp)); GooseFEM::Element::Quad4::QuadraturePlanar quad(vector.AsElement(coor)); GooseFEM::Element::Quad4::QuadraturePlanar quad_a(vector_a.AsElement(coor)); GooseFEM::Element::Quad4::QuadraturePlanar quad_b(vector_b.AsElement(coor)); size_t nip = quad.nip(); GMatElastic::Cartesian3d::Matrix mat(nelem, nip); GMatElastic::Cartesian3d::Matrix mat_a(nelem_a, nip, 3.0, 4.0); GMatElastic::Cartesian3d::Matrix mat_b(nelem_b, nip, 5.0, 6.0); xt::xtensor I = xt::empty({nelem, nip}); I.fill(0); xt::view(I, xt::keep(elem_a), xt::all()) = 1; mat.setElastic(I, 3.0, 4.0); I.fill(0); xt::view(I, xt::keep(elem_b), xt::all()) = 1; mat.setElastic(I, 5.0, 6.0); mat.check(); mat_a.check(); mat_b.check(); xt::xtensor Eps = quad.SymGradN_vector(vector.AsElement(disp)); xt::xtensor Eps_a = quad_a.SymGradN_vector(vector_a.AsElement(disp)); xt::xtensor Eps_b = quad_b.SymGradN_vector(vector_b.AsElement(disp)); xt::xtensor Sig = mat.Stress(Eps); xt::xtensor Sig_a = mat_a.Stress(Eps_a); xt::xtensor Sig_b = mat_b.Stress(Eps_b); xt::xtensor C; xt::xtensor C_a; xt::xtensor C_b; std::tie(Sig, C) = mat.Tangent(Eps); std::tie(Sig_a, C_a) = mat_a.Tangent(Eps_a); std::tie(Sig_b, C_b) = mat_b.Tangent(Eps_b); K.assemble(quad.Int_gradN_dot_tensor4_dot_gradNT_dV(C)); K_a.assemble(quad_a.Int_gradN_dot_tensor4_dot_gradNT_dV(C_a)); K_b.assemble(quad_b.Int_gradN_dot_tensor4_dot_gradNT_dV(C_b)); auto f = vector.AssembleNode(quad.Int_gradN_dot_tensor2_dV(Sig)); auto f_a = vector_a.AssembleNode(quad_a.Int_gradN_dot_tensor2_dV(Sig_a)); auto f_b = vector_b.AssembleNode(quad_b.Int_gradN_dot_tensor2_dV(Sig_b)); REQUIRE(xt::allclose(f, f_a + f_b)); REQUIRE(xt::allclose(f, K.Dot(disp))); REQUIRE(xt::allclose(f_a, K_a.Dot(disp))); REQUIRE(xt::allclose(f_b, K_b.Dot(disp))); } SECTION("VectorPartitioned/Matrix - GMatElastic - Periodic") { size_t N = 5; GooseFEM::Mesh::Quad4::Regular mesh(N, N); size_t nelem = mesh.nelem(); size_t ndim = mesh.ndim(); xt::xtensor elem_a = xt::arange(N); xt::xtensor elem_b = xt::arange(N, nelem); size_t nelem_a = elem_a.size(); size_t nelem_b = elem_b.size(); auto dofs = mesh.dofs(); // periodicity in horizontal direction : eliminate 'dependent' DOFs auto left = mesh.nodesLeftOpenEdge(); auto right = mesh.nodesRightOpenEdge(); xt::view(dofs, xt::keep(right), 0) = xt::view(dofs, xt::keep(left), 0); xt::view(dofs, xt::keep(right), 1) = xt::view(dofs, xt::keep(left), 1); // fixed top and bottom auto top = mesh.nodesTopEdge(); auto bottom = mesh.nodesBottomEdge(); size_t nfix = top.size(); xt::xtensor iip = xt::empty({2 * ndim * nfix}); xt::view(iip, xt::range(0 * nfix, 1 * nfix)) = xt::view(dofs, xt::keep(bottom), 0); xt::view(iip, xt::range(1 * nfix, 2 * nfix)) = xt::view(dofs, xt::keep(bottom), 1); xt::view(iip, xt::range(2 * nfix, 3 * nfix)) = xt::view(dofs, xt::keep(top), 0); xt::view(iip, xt::range(3 * nfix, 4 * nfix)) = xt::view(dofs, xt::keep(top), 1); xt::xtensor conn = mesh.conn(); xt::xtensor conn_a = xt::view(conn, xt::keep(elem_a), xt::all()); xt::xtensor conn_b = xt::view(conn, xt::keep(elem_b), xt::all()); GooseFEM::VectorPartitioned vector(conn, dofs, iip); GooseFEM::VectorPartitioned vector_a(conn_a, dofs, iip); GooseFEM::VectorPartitioned vector_b(conn_b, dofs, iip); GooseFEM::Matrix K(conn, dofs); GooseFEM::Matrix K_a(conn_a, dofs); GooseFEM::Matrix K_b(conn_b, dofs); xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::random::rand(coor.shape()); disp = vector.AsNode(vector.AsDofs(disp)); GooseFEM::Element::Quad4::QuadraturePlanar quad(vector.AsElement(coor)); GooseFEM::Element::Quad4::QuadraturePlanar quad_a(vector_a.AsElement(coor)); GooseFEM::Element::Quad4::QuadraturePlanar quad_b(vector_b.AsElement(coor)); size_t nip = quad.nip(); GMatElastic::Cartesian3d::Matrix mat(nelem, nip); GMatElastic::Cartesian3d::Matrix mat_a(nelem_a, nip, 3.0, 4.0); GMatElastic::Cartesian3d::Matrix mat_b(nelem_b, nip, 5.0, 6.0); xt::xtensor I = xt::empty({nelem, nip}); I.fill(0); xt::view(I, xt::keep(elem_a), xt::all()) = 1; mat.setElastic(I, 3.0, 4.0); I.fill(0); xt::view(I, xt::keep(elem_b), xt::all()) = 1; mat.setElastic(I, 5.0, 6.0); mat.check(); mat_a.check(); mat_b.check(); xt::xtensor Eps = quad.SymGradN_vector(vector.AsElement(disp)); xt::xtensor Eps_a = quad_a.SymGradN_vector(vector_a.AsElement(disp)); xt::xtensor Eps_b = quad_b.SymGradN_vector(vector_b.AsElement(disp)); xt::xtensor Sig = mat.Stress(Eps); xt::xtensor Sig_a = mat_a.Stress(Eps_a); xt::xtensor Sig_b = mat_b.Stress(Eps_b); xt::xtensor C; xt::xtensor C_a; xt::xtensor C_b; std::tie(Sig, C) = mat.Tangent(Eps); std::tie(Sig_a, C_a) = mat_a.Tangent(Eps_a); std::tie(Sig_b, C_b) = mat_b.Tangent(Eps_b); K.assemble(quad.Int_gradN_dot_tensor4_dot_gradNT_dV(C)); K_a.assemble(quad_a.Int_gradN_dot_tensor4_dot_gradNT_dV(C_a)); K_b.assemble(quad_b.Int_gradN_dot_tensor4_dot_gradNT_dV(C_b)); auto f = vector.AssembleNode(quad.Int_gradN_dot_tensor2_dV(Sig)); auto f_a = vector_a.AssembleNode(quad_a.Int_gradN_dot_tensor2_dV(Sig_a)); auto f_b = vector_b.AssembleNode(quad_b.Int_gradN_dot_tensor2_dV(Sig_b)); REQUIRE(xt::allclose(f, f_a + f_b)); REQUIRE(xt::allclose(f, K.Dot(disp))); REQUIRE(xt::allclose(f_a, K_a.Dot(disp))); REQUIRE(xt::allclose(f_b, K_b.Dot(disp))); } } diff --git a/test/main.cpp b/test/gmat/main.cpp similarity index 100% rename from test/main.cpp rename to test/gmat/main.cpp diff --git a/test/support.h b/test/support.h deleted file mode 100644 index 22c98e5..0000000 --- a/test/support.h +++ /dev/null @@ -1,16 +0,0 @@ - -#ifndef SUPPORT_H -#define SUPPORT_H - -#include - -#include -#include - -#include - -#include - -#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); - -#endif