diff --git a/.travis.yml b/.travis.yml index 0184050..ec1f419 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,98 +1,93 @@ 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 cmake - - conda install -c conda-forge xsimd - - conda install -c conda-forge xtensor - - conda install -c conda-forge python - - conda install -c conda-forge numpy - - conda install -c conda-forge pyxtensor - - conda install -c conda-forge catch2 - - conda install -c conda-forge eigen - - conda install -c conda-forge gmatelastic - - conda install -c conda-forge python-gmatelastic - - conda install -c conda-forge highfive + - 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 # 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 - ./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 + # - ./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/docs/examples/CMakeLists.txt b/docs/examples/CMakeLists.txt index 9525068..67c6041 100644 --- a/docs/examples/CMakeLists.txt +++ b/docs/examples/CMakeLists.txt @@ -1,77 +1,81 @@ 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) +option(DEBUG "Enable all assertions" OFF) 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) - # create executable set(exec "statics_MixedPeriodic_LinearElastic_example") set(source "statics/MixedPeriodic_LinearElastic/example.cpp") +add_executable(${exec} ${source}) +target_link_libraries(${exec} PRIVATE libraries) + +# create executable +set(exec "statics_Periodic_ElastoPlastic_main") +set(source "statics/Periodic_ElastoPlastic/main.cpp") add_executable(${exec} ${source}) +target_link_libraries(${exec} PRIVATE libraries) +# create executable + +set(exec "statics_Periodic_ElastoPlasticFiniteStrainSimo_main") +set(source "statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp") +add_executable(${exec} ${source}) +target_link_libraries(${exec} PRIVATE libraries) + +# create executable + +set(exec "statics_Periodic_LinearElastic_main") +set(source "statics/Periodic_LinearElastic/main.cpp") +add_executable(${exec} ${source}) +target_link_libraries(${exec} PRIVATE libraries) + +# create executable + +set(exec "statics_Periodic_NonLinearElastic_main") +set(source "statics/Periodic_NonLinearElastic/main.cpp") +add_executable(${exec} ${source}) target_link_libraries(${exec} PRIVATE libraries) -configure_file( - "statics/MixedPeriodic_LinearElastic/example.py" - "statics_MixedPeriodic_LinearElastic_example.py" COPYONLY) -configure_file( - "statics/MixedPeriodic_LinearElastic/plot.py" - "statics_MixedPeriodic_LinearElastic_plot.py" COPYONLY) diff --git a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp index a78b5a3..bef11f3 100644 --- a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp +++ b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp @@ -1,228 +1,229 @@ #include #include #include #include #include namespace GM = GMatElastoPlastic::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // 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(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor du = xt::zeros(coor.shape()); // iterative displacement update xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix 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}); // DOF values xt::xtensor Fext = xt::zeros({tyings.nni()}); xt::xtensor Fint = xt::zeros({tyings.nni()}); // element/material definition // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); xt::xtensor dV = elem.DV(2); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setLinearHardening(Isoft, 1.0, 1.0, 0.05, 0.05); mat.setElastic(Ihard, 1.0, 1.0); // solve // ----- // allocate tensors xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyingsSolver<> Solver; // allocate internal variables double res; // some shear strain history xt::xtensor dgamma = 0.001 * xt::ones({101}); dgamma(0) = 0.0; // initialise output // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // loop over increments for (size_t inc = 0; inc < dgamma.size(); ++inc) { // update history mat.increment(); // iterate for (size_t iter = 0; ; ++iter) { // 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); // residual xt::noalias(fres) = fext - fint; // check for convergence (skip the zeroth iteration, as the residual still vanishes) if (iter > 0) { // - internal/external force as DOFs (account for periodicity) vector.asDofs_i(fext, Fext); vector.asDofs_i(fint, Fint); // - extract reaction force vector.copy_p(Fint, Fext); // - norm of the residual and the reaction force double nfres = xt::sum(xt::abs(Fext-Fint))[0]; double nfext = xt::sum(xt::abs(Fext))[0]; // - relative residual, for convergence check if (nfext) res = nfres / nfext; else res = nfres; // - print progress to screen std::cout << "inc = " << inc << ", iter = " << iter << ", res = " << res << std::endl; // - check for convergence if (res < 1.0e-5) break; // - safe-guard from infinite loop if (iter > 20) throw std::runtime_error("Maximal number of iterations exceeded"); } // initialise displacement update du.fill(0.0); // set fixed displacements if (iter == 0) du(control_nodes(0),1) = dgamma(inc); // solve - K.solve(fres, du); + Solver.solve(K, fres, du); // add displacement update disp += du; } // post-process // - compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: increment numbers H5::dump(file, "/stored", inc, {inc}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar), {inc}); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar), {inc}); // - write to output-file: element quantities H5::dump(file, "/sigeq/" + std::to_string(inc), GM::Sigeq(Sigelem)); H5::dump(file, "/epseq/" + std::to_string(inc), GM::Epseq(Epselem)); H5::dump(file, "/disp/" + std::to_string(inc), PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp/" + std::to_string(inc), "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq/" + std::to_string(inc), "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq/" + std::to_string(inc), "Eq. strain" , PV::AttributeType::Cell), } )); } // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp index 91112b3..2588625 100644 --- a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp +++ b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp @@ -1,238 +1,239 @@ #include #include #include #include #include namespace GM = GMatElastoPlasticFiniteStrainSimo::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // 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(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor du = xt::zeros(coor.shape()); // iterative displacement update xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix 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}); // DOF values xt::xtensor Fext = xt::zeros({tyings.nni()}); xt::xtensor Fint = xt::zeros({tyings.nni()}); // element/material definition // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem0(vector.AsElement(coor)); GF::Element::Quad4::QuadraturePlanar elem (vector.AsElement(coor)); size_t nip = elem0.nip(); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setLinearHardening(Isoft, 1.0, 1.0, 0.05, 0.05); mat.setElastic(Ihard, 1.0, 1.0); // solve // ----- // allocate tensors xt::xtensor I2 = mat.I2(); xt::xtensor F = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyingsSolver<> Solver; // allocate internal variables double res; // some shear strain history xt::xtensor dgamma = 0.001 * xt::ones({101}); dgamma(0) = 0.0; // initialise output // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // loop over increments for (size_t inc = 0; inc < dgamma.size(); ++inc) { // update history mat.increment(); // iterate for (size_t iter = 0; ; ++iter) { // deformation gradient tensor vector.asElement(disp, ue); elem0.gradN_vector_T(ue, F); F += I2; // stress & tangent mat.tangent(F, 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); // residual xt::noalias(fres) = fext - fint; // check for convergence (skip the zeroth iteration, as the residual still vanishes) if (iter > 0) { // - internal/external force as DOFs (account for periodicity) vector.asDofs_i(fext, Fext); vector.asDofs_i(fint, Fint); // - extract reaction force vector.copy_p(Fint, Fext); // - norm of the residual and the reaction force double nfres = xt::sum(xt::abs(Fext-Fint))[0]; double nfext = xt::sum(xt::abs(Fext))[0]; // - relative residual, for convergence check if (nfext) res = nfres / nfext; else res = nfres; // - print progress to screen std::cout << "inc = " << inc << ", iter = " << iter << ", res = " << res << std::endl; // - check for convergence if (res < 1.0e-5) break; // - safe-guard from infinite loop if (iter > 20) throw std::runtime_error("Maximal number of iterations exceeded"); } // initialise displacement update du.fill(0.0); // set fixed displacements if (iter == 0) du(control_nodes(0),1) = dgamma(inc); // solve - K.solve(fres, du); + Solver.solve(K, fres, du); // add displacement update disp += du; // update shape functions elem.update_x(vector.AsElement(coor+disp)); } // post-process // - compute strain and stress vector.asElement(disp, ue); elem0.gradN_vector_T(ue, F); F += I2; GM::strain(F, Eps); mat.stress(F, Sig); // - integration point volume xt::xtensor dV = elem.DV(2); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: increment numbers H5::dump(file, "/stored", inc, {inc}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar), {inc}); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar), {inc}); // - write to output-file: element quantities H5::dump(file, "/sigeq/" + std::to_string(inc), GM::Sigeq(Sigelem)); H5::dump(file, "/epseq/" + std::to_string(inc), GM::Epseq(Epselem)); H5::dump(file, "/disp/" + std::to_string(inc), PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp/" + std::to_string(inc), "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq/" + std::to_string(inc), "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq/" + std::to_string(inc), "Eq. strain" , PV::AttributeType::Cell), } )); } // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/docs/examples/statics/Periodic_LinearElastic/main.cpp b/docs/examples/statics/Periodic_LinearElastic/main.cpp index f9fbfe6..702e3a0 100644 --- a/docs/examples/statics/Periodic_LinearElastic/main.cpp +++ b/docs/examples/statics/Periodic_LinearElastic/main.cpp @@ -1,170 +1,171 @@ #include #include #include #include #include namespace GM = GMatElastic::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // 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(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix 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 // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setElastic(Isoft, 10.0, 1.0); mat.setElastic(Ihard, 10.0, 10.0); // solve // ----- // allocate tensors xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyingsSolver<> Solver; // 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); // residual xt::noalias(fres) = fext - fint; // set fixed displacements disp(control_nodes(0),1) = 0.1; // solve - K.solve(fres, disp); + Solver.solve(K, fres, disp); // post-process // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // - integration point volume xt::xtensor dV = elem.DV(2); // - compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar)); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar)); // - write to output-file: element quantities H5::dump(file, "/sigeq", GM::Sigeq(Sigelem)); H5::dump(file, "/epseq", GM::Epseq(Epselem)); H5::dump(file, "/disp" , PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp" , "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq", "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq", "Eq. strain" , PV::AttributeType::Cell), } )); // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/docs/examples/statics/Periodic_NonLinearElastic/main.cpp b/docs/examples/statics/Periodic_NonLinearElastic/main.cpp index e5495a7..46667d8 100644 --- a/docs/examples/statics/Periodic_NonLinearElastic/main.cpp +++ b/docs/examples/statics/Periodic_NonLinearElastic/main.cpp @@ -1,214 +1,215 @@ #include #include #include #include #include namespace GM = GMatNonLinearElastic::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // 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(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor du = xt::zeros(coor.shape()); // iterative displacement update xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix 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}); // DOF values xt::xtensor Fext = xt::zeros({tyings.nni()}); xt::xtensor Fint = xt::zeros({tyings.nni()}); // element/material definition // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setNonLinearElastic(Isoft, 10.0, 0.1, 0.1, 2.0); mat.setNonLinearElastic(Ihard, 10.0, 1.0, 0.1, 2.0); // solve // ----- // allocate tensors xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyingsSolver<> Solver; // allocate internal variables double res; // iterate for (size_t iter = 0; ; ++iter) { // 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); // residual xt::noalias(fres) = fext - fint; // check for convergence (skip the zeroth iteration, as the residual still vanishes) if (iter > 0) { // - internal/external force as DOFs (account for periodicity) vector.asDofs_i(fext, Fext); vector.asDofs_i(fint, Fint); // - extract reaction force vector.copy_p(Fint, Fext); // - norm of the residual and the reaction force double nfres = xt::sum(xt::abs(Fext-Fint))[0]; double nfext = xt::sum(xt::abs(Fext))[0]; // - relative residual, for convergence check if (nfext) res = nfres / nfext; else res = nfres; // - print progress to screen std::cout << "iter = " << iter << ", res = " << res << std::endl; // - check for convergence if (res < 1.0e-5) break; // - safe-guard from infinite loop if (iter > 20) throw std::runtime_error("Maximal number of iterations exceeded"); } // initialise displacement update du.fill(0.0); // set fixed displacements if (iter == 0) du(control_nodes(0),1) = 0.1; // solve - K.solve(fres, du); + Solver.solve(K, fres, du); // add displacement update disp += du; } // post-process // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // - integration point volume xt::xtensor dV = elem.DV(2); // - compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar)); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar)); // - write to output-file: element quantities H5::dump(file, "/sigeq", GM::Sigeq(Sigelem)); H5::dump(file, "/epseq", GM::Epseq(Epselem)); H5::dump(file, "/disp" , PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp" , "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq", "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq", "Eq. strain" , PV::AttributeType::Cell), } )); // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/include/GooseFEM/MatrixPartitionedTyings.h b/include/GooseFEM/MatrixPartitionedTyings.h index 8b3fa0c..effbdec 100644 --- a/include/GooseFEM/MatrixPartitionedTyings.h +++ b/include/GooseFEM/MatrixPartitionedTyings.h @@ -1,150 +1,177 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_H #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_H #include "config.h" #include #include #include namespace GooseFEM { -template >> +// forward declaration +template class MatrixPartitionedTyingsSolver; + class MatrixPartitionedTyings { public: // Constructors MatrixPartitionedTyings() = default; MatrixPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp); // Dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs size_t nnu() const; // number of independent, unknown DOFs size_t nnp() const; // number of independent, prescribed DOFs size_t nni() const; // number of independent DOFs size_t nnd() const; // number of dependent DOFs // DOF lists xt::xtensor dofs() const; // DOFs xt::xtensor iiu() const; // independent, unknown DOFs xt::xtensor iip() const; // independent, prescribed DOFs xt::xtensor iii() const; // independent DOFs xt::xtensor iid() const; // dependent DOFs // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] void assemble(const xt::xtensor& elemmat); - // Solve: - // A' = A_ii + K_id * C_di + C_di^T * K_di + C_di^T * K_dd * C_di - // b' = b_i + C_di^T * b_d - // x_u = A'_uu \ ( b'_u - A'_up * x_p ) - // x_i = [x_u, x_p] - // x_d = C_di * x_i - void solve(const xt::xtensor& b, xt::xtensor& x); // updates x_u and x_d - void solve(const xt::xtensor& b, xt::xtensor& x); // updates x_u and x_d - - void solve_u( - const xt::xtensor& b_u, - const xt::xtensor& b_d, - const xt::xtensor& x_p, - xt::xtensor& x_u); - - // Auto-allocation of the functions above - xt::xtensor Solve(const xt::xtensor& b, const xt::xtensor& x); - xt::xtensor Solve(const xt::xtensor& b, const xt::xtensor& x); - - xt::xtensor Solve_u( - const xt::xtensor& b_u, - const xt::xtensor& b_d, - const xt::xtensor& x_p); - private: // The matrix Eigen::SparseMatrix m_Auu; Eigen::SparseMatrix m_Aup; Eigen::SparseMatrix m_Apu; Eigen::SparseMatrix m_App; Eigen::SparseMatrix m_Aud; Eigen::SparseMatrix m_Apd; Eigen::SparseMatrix m_Adu; Eigen::SparseMatrix m_Adp; Eigen::SparseMatrix m_Add; // The matrix for which the tyings have been applied Eigen::SparseMatrix m_ACuu; Eigen::SparseMatrix m_ACup; Eigen::SparseMatrix m_ACpu; Eigen::SparseMatrix m_ACpp; // Matrix entries std::vector> m_Tuu; std::vector> m_Tup; std::vector> m_Tpu; std::vector> m_Tpp; std::vector> m_Tud; std::vector> m_Tpd; std::vector> m_Tdu; std::vector> m_Tdp; std::vector> m_Tdd; - // Solver (re-used to solve different RHS) - Solver m_solver; - - // Signal changes to data compare to the last inverse - bool m_factor = false; + // Signal changes to data + bool m_changed = true; // Bookkeeping xt::xtensor m_conn; // connectivity [nelem, nne ] xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] xt::xtensor m_iiu; // unknown DOFs [nnu] xt::xtensor m_iip; // prescribed DOFs [nnp] xt::xtensor m_iid; // dependent DOFs [nnd] // Dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs size_t m_nnu; // number of independent, unknown DOFs size_t m_nnp; // number of independent, prescribed DOFs size_t m_nni; // number of independent DOFs size_t m_nnd; // number of dependent DOFs // Tyings Eigen::SparseMatrix m_Cdu; Eigen::SparseMatrix m_Cdp; Eigen::SparseMatrix m_Cud; Eigen::SparseMatrix m_Cpd; - // Compute inverse (automatically evaluated by "solve") - void factorize(); + // grant access to solver class + template friend class MatrixPartitionedTyingsSolver; // Convert arrays (Eigen version of VectorPartitioned, which contains public functions) Eigen::VectorXd asDofs_u(const xt::xtensor& dofval) const; Eigen::VectorXd asDofs_u(const xt::xtensor& nodevec) const; Eigen::VectorXd asDofs_p(const xt::xtensor& dofval) const; Eigen::VectorXd asDofs_p(const xt::xtensor& nodevec) const; Eigen::VectorXd asDofs_d(const xt::xtensor& dofval) const; Eigen::VectorXd asDofs_d(const xt::xtensor& nodevec) const; }; +template >> +class MatrixPartitionedTyingsSolver { +public: + // Constructors + MatrixPartitionedTyingsSolver() = default; + + // Solve: + // A' = A_ii + K_id * C_di + C_di^T * K_di + C_di^T * K_dd * C_di + // b' = b_i + C_di^T * b_d + // x_u = A'_uu \ ( b'_u - A'_up * x_p ) + // x_i = [x_u, x_p] + // x_d = C_di * x_i + void solve( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b, + xt::xtensor& x); // updates x_u and x_d + + void solve( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b, + xt::xtensor& x); // updates x_u and x_d + + void solve_u( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b_u, + const xt::xtensor& b_d, + const xt::xtensor& x_p, + xt::xtensor& x_u); + + // Auto-allocation of the functions above + xt::xtensor Solve( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b, + const xt::xtensor& x); + + xt::xtensor Solve( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b, + const xt::xtensor& x); + + xt::xtensor Solve_u( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b_u, + const xt::xtensor& b_d, + const xt::xtensor& x_p); + +private: + Solver m_solver; // solver + bool m_factor = false; // signal to force factorization + void factorize(MatrixPartitionedTyings& matrix); // compute inverse (evaluated by "solve") +}; + } // namespace GooseFEM #include "MatrixPartitionedTyings.hpp" #endif diff --git a/include/GooseFEM/MatrixPartitionedTyings.hpp b/include/GooseFEM/MatrixPartitionedTyings.hpp index 676e5de..7f4d1ed 100644 --- a/include/GooseFEM/MatrixPartitionedTyings.hpp +++ b/include/GooseFEM/MatrixPartitionedTyings.hpp @@ -1,475 +1,459 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP #include "MatrixPartitionedTyings.h" namespace GooseFEM { -template -inline MatrixPartitionedTyings::MatrixPartitionedTyings( +inline MatrixPartitionedTyings::MatrixPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp) : m_conn(conn), m_dofs(dofs), m_Cdu(Cdu), m_Cdp(Cdp) { GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows()); m_nnu = static_cast(m_Cdu.cols()); m_nnp = static_cast(m_Cdp.cols()); m_nnd = static_cast(m_Cdp.rows()); m_nni = m_nnu + m_nnp; m_ndof = m_nni + m_nnd; m_iiu = xt::arange(m_nnu); m_iip = xt::arange(m_nnu, m_nnu + m_nnp); m_iid = xt::arange(m_nni, m_nni + m_nnd); m_nelem = m_conn.shape(0); m_nne = m_conn.shape(1); m_nnode = m_dofs.shape(0); m_ndim = m_dofs.shape(1); m_Cud = m_Cdu.transpose(); m_Cpd = m_Cdp.transpose(); m_Tuu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tup.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tpu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tpp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tud.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tpd.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tdu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tdp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Tdd.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim); m_Auu.resize(m_nnu, m_nnu); m_Aup.resize(m_nnu, m_nnp); m_Apu.resize(m_nnp, m_nnu); m_App.resize(m_nnp, m_nnp); m_Aud.resize(m_nnu, m_nnd); m_Apd.resize(m_nnp, m_nnd); m_Adu.resize(m_nnd, m_nnu); m_Adp.resize(m_nnd, m_nnp); m_Add.resize(m_nnd, m_nnd); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)[0] + 1); } -template -inline size_t MatrixPartitionedTyings::nelem() const +inline size_t MatrixPartitionedTyings::nelem() const { return m_nelem; } -template -inline size_t MatrixPartitionedTyings::nne() const +inline size_t MatrixPartitionedTyings::nne() const { return m_nne; } -template -inline size_t MatrixPartitionedTyings::nnode() const +inline size_t MatrixPartitionedTyings::nnode() const { return m_nnode; } -template -inline size_t MatrixPartitionedTyings::ndim() const +inline size_t MatrixPartitionedTyings::ndim() const { return m_ndim; } -template -inline size_t MatrixPartitionedTyings::ndof() const +inline size_t MatrixPartitionedTyings::ndof() const { return m_ndof; } -template -inline size_t MatrixPartitionedTyings::nnu() const +inline size_t MatrixPartitionedTyings::nnu() const { return m_nnu; } -template -inline size_t MatrixPartitionedTyings::nnp() const +inline size_t MatrixPartitionedTyings::nnp() const { return m_nnp; } -template -inline size_t MatrixPartitionedTyings::nni() const +inline size_t MatrixPartitionedTyings::nni() const { return m_nni; } -template -inline size_t MatrixPartitionedTyings::nnd() const +inline size_t MatrixPartitionedTyings::nnd() const { return m_nnd; } -template -inline xt::xtensor MatrixPartitionedTyings::dofs() const +inline xt::xtensor MatrixPartitionedTyings::dofs() const { return m_dofs; } -template -inline xt::xtensor MatrixPartitionedTyings::iiu() const +inline xt::xtensor MatrixPartitionedTyings::iiu() const { return m_iiu; } -template -inline xt::xtensor MatrixPartitionedTyings::iip() const +inline xt::xtensor MatrixPartitionedTyings::iip() const { return m_iip; } -template -inline xt::xtensor MatrixPartitionedTyings::iii() const +inline xt::xtensor MatrixPartitionedTyings::iii() const { return xt::arange(m_nni); } -template -inline xt::xtensor MatrixPartitionedTyings::iid() const +inline xt::xtensor MatrixPartitionedTyings::iid() const { return m_iid; } -template -inline void MatrixPartitionedTyings::factorize() -{ - if (!m_factor) { - return; - } - - m_ACuu = m_Auu + m_Aud * m_Cdu + m_Cud * m_Adu + m_Cud * m_Add * m_Cdu; - m_ACup = m_Aup + m_Aud * m_Cdp + m_Cud * m_Adp + m_Cud * m_Add * m_Cdp; - // m_ACpu = m_Apu + m_Apd * m_Cdu + m_Cpd * m_Adu + m_Cpd * m_Add * m_Cdu; - // m_ACpp = m_App + m_Apd * m_Cdp + m_Cpd * m_Adp + m_Cpd * m_Add * m_Cdp; - - m_solver.compute(m_ACuu); - - m_factor = false; -} - -template -inline void MatrixPartitionedTyings::assemble(const xt::xtensor& elemmat) +inline void MatrixPartitionedTyings::assemble(const xt::xtensor& elemmat) { GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); m_Tuu.clear(); m_Tup.clear(); m_Tpu.clear(); m_Tpp.clear(); m_Tud.clear(); m_Tpd.clear(); m_Tdu.clear(); m_Tdp.clear(); m_Tdd.clear(); for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { size_t di = m_dofs(m_conn(e, m), i); for (size_t n = 0; n < m_nne; ++n) { for (size_t j = 0; j < m_ndim; ++j) { size_t dj = m_dofs(m_conn(e, n), j); if (di < m_nnu && dj < m_nnu) { m_Tuu.push_back(Eigen::Triplet( di, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else if (di < m_nnu && dj < m_nni) { m_Tup.push_back(Eigen::Triplet( di, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else if (di < m_nnu) { m_Tud.push_back(Eigen::Triplet( di, dj - m_nni, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else if (di < m_nni && dj < m_nnu) { m_Tpu.push_back(Eigen::Triplet( di - m_nnu, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else if (di < m_nni && dj < m_nni) { m_Tpp.push_back(Eigen::Triplet( di - m_nnu, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else if (di < m_nni) { m_Tpd.push_back(Eigen::Triplet( di - m_nnu, dj - m_nni, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else if (dj < m_nnu) { m_Tdu.push_back(Eigen::Triplet( di - m_nni, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else if (dj < m_nni) { m_Tdp.push_back(Eigen::Triplet( di - m_nni, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else { m_Tdd.push_back(Eigen::Triplet( di - m_nni, dj - m_nni, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } } } } } } m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end()); m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end()); m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end()); m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end()); m_Aud.setFromTriplets(m_Tud.begin(), m_Tud.end()); m_Apd.setFromTriplets(m_Tpd.begin(), m_Tpd.end()); m_Adu.setFromTriplets(m_Tdu.begin(), m_Tdu.end()); m_Adp.setFromTriplets(m_Tdp.begin(), m_Tdp.end()); m_Add.setFromTriplets(m_Tdd.begin(), m_Tdd.end()); - - m_factor = true; -} - -template -inline void -MatrixPartitionedTyings::solve(const xt::xtensor& b, xt::xtensor& x) -{ - GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim})); - GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim})); - - this->factorize(); - - Eigen::VectorXd B_u = this->asDofs_u(b); - Eigen::VectorXd B_d = this->asDofs_d(b); - Eigen::VectorXd X_p = this->asDofs_p(x); - - B_u += m_Cud * B_d; - - Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p)); - Eigen::VectorXd X_d = m_Cdu * X_u + m_Cdp * X_p; - - #pragma omp parallel for - for (size_t m = 0; m < m_nnode; ++m) { - for (size_t i = 0; i < m_ndim; ++i) { - if (m_dofs(m, i) < m_nnu) { - x(m, i) = X_u(m_dofs(m, i)); - } - else if (m_dofs(m, i) >= m_nni) { - x(m, i) = X_d(m_dofs(m, i) - m_nni); - } - } - } -} - -template -inline void -MatrixPartitionedTyings::solve(const xt::xtensor& b, xt::xtensor& x) -{ - GOOSEFEM_ASSERT(b.size() == m_ndof); - GOOSEFEM_ASSERT(x.size() == m_ndof); - - this->factorize(); - - Eigen::VectorXd B_u = this->asDofs_u(b); - Eigen::VectorXd B_d = this->asDofs_d(b); - Eigen::VectorXd X_p = this->asDofs_p(x); - - Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p)); - Eigen::VectorXd X_d = m_Cdu * X_u + m_Cdp * X_p; - - #pragma omp parallel for - for (size_t d = 0; d < m_nnu; ++d) { - x(m_iiu(d)) = X_u(d); - } - - #pragma omp parallel for - for (size_t d = 0; d < m_nnd; ++d) { - x(m_iid(d)) = X_d(d); - } -} - -template -inline void MatrixPartitionedTyings::solve_u( - const xt::xtensor& b_u, - const xt::xtensor& b_d, - const xt::xtensor& x_p, - xt::xtensor& x_u) -{ - GOOSEFEM_ASSERT(b_u.size() == m_nnu); - GOOSEFEM_ASSERT(b_d.size() == m_nnd); - GOOSEFEM_ASSERT(x_p.size() == m_nnp); - GOOSEFEM_ASSERT(x_u.size() == m_nnu); - - this->factorize(); - - Eigen::VectorXd B_u(m_nnu, 1); - Eigen::VectorXd B_d(m_nnd, 1); - Eigen::VectorXd X_p(m_nnp, 1); - - std::copy(b_u.begin(), b_u.end(), B_u.data()); - std::copy(b_d.begin(), b_d.end(), B_d.data()); - std::copy(x_p.begin(), x_p.end(), X_p.data()); - - Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p)); - - std::copy(X_u.data(), X_u.data() + m_nnu, x_u.begin()); -} - -template -inline xt::xtensor MatrixPartitionedTyings::Solve( - const xt::xtensor& b, const xt::xtensor& x) -{ - xt::xtensor out = x; - this->solve(b, out); - return out; -} - -template -inline xt::xtensor MatrixPartitionedTyings::Solve( - const xt::xtensor& b, const xt::xtensor& x) -{ - xt::xtensor out = x; - this->solve(b, out); - return out; -} - -template -inline xt::xtensor MatrixPartitionedTyings::Solve_u( - const xt::xtensor& b_u, - const xt::xtensor& b_d, - const xt::xtensor& x_p) -{ - xt::xtensor x_u = xt::empty({m_nnu}); - this->solve_u(b_u, b_d, x_p, x_u); - return x_u; + m_changed = true; } -template -inline Eigen::VectorXd -MatrixPartitionedTyings::asDofs_u(const xt::xtensor& dofval) const +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_u(const xt::xtensor& dofval) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); Eigen::VectorXd dofval_u(m_nnu, 1); #pragma omp parallel for for (size_t d = 0; d < m_nnu; ++d) { dofval_u(d) = dofval(m_iiu(d)); } return dofval_u; } -template -inline Eigen::VectorXd -MatrixPartitionedTyings::asDofs_u(const xt::xtensor& nodevec) const +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_u(const xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); Eigen::VectorXd dofval_u(m_nnu, 1); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_dofs(m, i) < m_nnu) { dofval_u(m_dofs(m, i)) = nodevec(m, i); } } } return dofval_u; } -template -inline Eigen::VectorXd -MatrixPartitionedTyings::asDofs_p(const xt::xtensor& dofval) const +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_p(const xt::xtensor& dofval) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); Eigen::VectorXd dofval_p(m_nnp, 1); #pragma omp parallel for for (size_t d = 0; d < m_nnp; ++d) { dofval_p(d) = dofval(m_iip(d)); } return dofval_p; } -template -inline Eigen::VectorXd -MatrixPartitionedTyings::asDofs_p(const xt::xtensor& nodevec) const +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_p(const xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); Eigen::VectorXd dofval_p(m_nnp, 1); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_dofs(m, i) >= m_nnu && m_dofs(m, i) < m_nni) { dofval_p(m_dofs(m, i) - m_nnu) = nodevec(m, i); } } } return dofval_p; } -template -inline Eigen::VectorXd -MatrixPartitionedTyings::asDofs_d(const xt::xtensor& dofval) const +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d(const xt::xtensor& dofval) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); Eigen::VectorXd dofval_d(m_nnd, 1); #pragma omp parallel for for (size_t d = 0; d < m_nnd; ++d) { dofval_d(d) = dofval(m_iip(d)); } return dofval_d; } -template -inline Eigen::VectorXd -MatrixPartitionedTyings::asDofs_d(const xt::xtensor& nodevec) const +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d(const xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); Eigen::VectorXd dofval_d(m_nnd, 1); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_dofs(m, i) >= m_nni) { dofval_d(m_dofs(m, i) - m_nni) = nodevec(m, i); } } } return dofval_d; } +template +inline void MatrixPartitionedTyingsSolver::factorize(MatrixPartitionedTyings& matrix) +{ + if (!matrix.m_changed && !m_factor) { + return; + } + + matrix.m_ACuu = matrix.m_Auu + matrix.m_Aud * matrix.m_Cdu + matrix.m_Cud * matrix.m_Adu + + matrix.m_Cud * matrix.m_Add * matrix.m_Cdu; + + matrix.m_ACup = matrix.m_Aup + matrix.m_Aud * matrix.m_Cdp + matrix.m_Cud * matrix.m_Adp + + matrix.m_Cud * matrix.m_Add * matrix.m_Cdp; + + // matrix.m_ACpu = matrix.m_Apu + matrix.m_Apd * matrix.m_Cdu + matrix.m_Cpd * matrix.m_Adu + // + matrix.m_Cpd * matrix.m_Add * matrix.m_Cdu; + + // matrix.m_ACpp = matrix.m_App + matrix.m_Apd * matrix.m_Cdp + matrix.m_Cpd * matrix.m_Adp + // + matrix.m_Cpd * matrix.m_Add * matrix.m_Cdp; + + m_solver.compute(matrix.m_ACuu); + m_factor = false; + matrix.m_changed = false; +} + +template +inline void MatrixPartitionedTyingsSolver::solve( + MatrixPartitionedTyings& matrix, const xt::xtensor& b, xt::xtensor& x) +{ + GOOSEFEM_ASSERT(xt::has_shape(b, {matrix.m_nnode, matrix.m_ndim})); + GOOSEFEM_ASSERT(xt::has_shape(x, {matrix.m_nnode, matrix.m_ndim})); + + this->factorize(matrix); + + Eigen::VectorXd B_u = matrix.asDofs_u(b); + Eigen::VectorXd B_d = matrix.asDofs_d(b); + Eigen::VectorXd X_p = matrix.asDofs_p(x); + + B_u += matrix.m_Cud * B_d; + + Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_ACup * X_p)); + Eigen::VectorXd X_d = matrix.m_Cdu * X_u + matrix.m_Cdp * X_p; + + #pragma omp parallel for + for (size_t m = 0; m < matrix.m_nnode; ++m) { + for (size_t i = 0; i < matrix.m_ndim; ++i) { + if (matrix.m_dofs(m, i) < matrix.m_nnu) { + x(m, i) = X_u(matrix.m_dofs(m, i)); + } + else if (matrix.m_dofs(m, i) >= matrix.m_nni) { + x(m, i) = X_d(matrix.m_dofs(m, i) - matrix.m_nni); + } + } + } +} + +template +inline void MatrixPartitionedTyingsSolver::solve( + MatrixPartitionedTyings& matrix, const xt::xtensor& b, xt::xtensor& x) +{ + GOOSEFEM_ASSERT(b.size() == matrix.m_ndof); + GOOSEFEM_ASSERT(x.size() == matrix.m_ndof); + + this->factorize(matrix); + + Eigen::VectorXd B_u = matrix.asDofs_u(b); + Eigen::VectorXd B_d = matrix.asDofs_d(b); + Eigen::VectorXd X_p = matrix.asDofs_p(x); + + Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_ACup * X_p)); + Eigen::VectorXd X_d = matrix.m_Cdu * X_u + matrix.m_Cdp * X_p; + + #pragma omp parallel for + for (size_t d = 0; d < matrix.m_nnu; ++d) { + x(matrix.m_iiu(d)) = X_u(d); + } + + #pragma omp parallel for + for (size_t d = 0; d < matrix.m_nnd; ++d) { + x(matrix.m_iid(d)) = X_d(d); + } +} + +template +inline void MatrixPartitionedTyingsSolver::solve_u( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b_u, + const xt::xtensor& b_d, + const xt::xtensor& x_p, + xt::xtensor& x_u) +{ + GOOSEFEM_ASSERT(b_u.size() == matrix.m_nnu); + GOOSEFEM_ASSERT(b_d.size() == matrix.m_nnd); + GOOSEFEM_ASSERT(x_p.size() == matrix.m_nnp); + GOOSEFEM_ASSERT(x_u.size() == matrix.m_nnu); + + this->factorize(matrix); + + Eigen::VectorXd B_u(matrix.m_nnu, 1); + Eigen::VectorXd B_d(matrix.m_nnd, 1); + Eigen::VectorXd X_p(matrix.m_nnp, 1); + + std::copy(b_u.begin(), b_u.end(), B_u.data()); + std::copy(b_d.begin(), b_d.end(), B_d.data()); + std::copy(x_p.begin(), x_p.end(), X_p.data()); + + Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_ACup * X_p)); + + std::copy(X_u.data(), X_u.data() + matrix.m_nnu, x_u.begin()); +} + +template +inline xt::xtensor MatrixPartitionedTyingsSolver::Solve( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b, + const xt::xtensor& x) +{ + xt::xtensor out = x; + this->solve(matrix, b, out); + return out; +} + +template +inline xt::xtensor MatrixPartitionedTyingsSolver::Solve( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b, + const xt::xtensor& x) +{ + xt::xtensor out = x; + this->solve(matrix, b, out); + return out; +} + +template +inline xt::xtensor MatrixPartitionedTyingsSolver::Solve_u( + MatrixPartitionedTyings& matrix, + const xt::xtensor& b_u, + const xt::xtensor& b_d, + const xt::xtensor& x_p) +{ + xt::xtensor x_u = xt::empty({matrix.m_nnu}); + this->solve_u(matrix, b_u, b_d, x_p, x_u); + return x_u; +} + } // namespace GooseFEM #endif diff --git a/python/MatrixPartitionedTyings.hpp b/python/MatrixPartitionedTyings.hpp index ed15013..105f4a3 100644 --- a/python/MatrixPartitionedTyings.hpp +++ b/python/MatrixPartitionedTyings.hpp @@ -1,128 +1,148 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include #include "../include/GooseFEM/GooseFEM.h" // ================================================================================================= namespace py = pybind11; // ================================================================================================= void init_MatrixPartitionedTyings(py::module& m) { -py::class_>(m, "MatrixPartitionedTyings") +py::class_(m, "MatrixPartitionedTyings") .def(py::init< const xt::xtensor&, const xt::xtensor&, const Eigen::SparseMatrix&, const Eigen::SparseMatrix&>(), "Sparse, partitioned, matrix", py::arg("conn"), py::arg("dofs"), py::arg("Cdu"), py::arg("Cdp")) .def("nelem", - &GooseFEM::MatrixPartitionedTyings<>::nelem, + &GooseFEM::MatrixPartitionedTyings::nelem, "Number of element") .def("nne", - &GooseFEM::MatrixPartitionedTyings<>::nne, + &GooseFEM::MatrixPartitionedTyings::nne, "Number of nodes per element") .def("nnode", - &GooseFEM::MatrixPartitionedTyings<>::nnode, + &GooseFEM::MatrixPartitionedTyings::nnode, "Number of nodes") .def("ndim", - &GooseFEM::MatrixPartitionedTyings<>::ndim, + &GooseFEM::MatrixPartitionedTyings::ndim, "Number of dimensions") .def("ndof", - &GooseFEM::MatrixPartitionedTyings<>::ndof, + &GooseFEM::MatrixPartitionedTyings::ndof, "Number of degrees-of-freedom") .def("nnu", - &GooseFEM::MatrixPartitionedTyings<>::nnu, + &GooseFEM::MatrixPartitionedTyings::nnu, "Number of unknown degrees-of-freedom") .def("nnp", - &GooseFEM::MatrixPartitionedTyings<>::nnp, + &GooseFEM::MatrixPartitionedTyings::nnp, "Number of prescribed degrees-of-freedom") .def("nni", - &GooseFEM::MatrixPartitionedTyings<>::nni, + &GooseFEM::MatrixPartitionedTyings::nni, "Number of independent degrees-of-freedom") .def("nnd", - &GooseFEM::MatrixPartitionedTyings<>::nnd, + &GooseFEM::MatrixPartitionedTyings::nnd, "Number of dependent degrees-of-freedom") .def("assemble", - &GooseFEM::MatrixPartitionedTyings<>::assemble, + &GooseFEM::MatrixPartitionedTyings::assemble, "Assemble matrix from 'elemmat", py::arg("elemmat")) .def("dofs", - &GooseFEM::MatrixPartitionedTyings<>::dofs, + &GooseFEM::MatrixPartitionedTyings::dofs, "Return degrees-of-freedom") .def("iiu", - &GooseFEM::MatrixPartitionedTyings<>::iiu, + &GooseFEM::MatrixPartitionedTyings::iiu, "Return unknown degrees-of-freedom") .def("iip", - &GooseFEM::MatrixPartitionedTyings<>::iip, + &GooseFEM::MatrixPartitionedTyings::iip, "Return prescribed degrees-of-freedom") .def("iii", - &GooseFEM::MatrixPartitionedTyings<>::iii, + &GooseFEM::MatrixPartitionedTyings::iii, "Return independent degrees-of-freedom") .def("iid", - &GooseFEM::MatrixPartitionedTyings<>::iid, + &GooseFEM::MatrixPartitionedTyings::iid, "Return dependent degrees-of-freedom") + .def("__repr__", + [](const GooseFEM::MatrixPartitionedTyings&){ + return ""; }); + + +py::class_>(m, "MatrixPartitionedTyingsSolver") + + .def(py::init<>(), + "Sparse, partitioned, matrix solver") + .def("Solve", - py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitionedTyings<>::Solve), + py::overload_cast< + GooseFEM::MatrixPartitionedTyings&, + const xt::xtensor&, + const xt::xtensor&>( + &GooseFEM::MatrixPartitionedTyingsSolver<>::Solve), "Solve", + py::arg("matrix"), py::arg("b"), py::arg("x")) .def("Solve", - py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitionedTyings<>::Solve), + py::overload_cast< + GooseFEM::MatrixPartitionedTyings&, + const xt::xtensor&, + const xt::xtensor&>( + &GooseFEM::MatrixPartitionedTyingsSolver<>::Solve), "Solve", + py::arg("matrix"), py::arg("b"), py::arg("x")) .def("Solve_u", py::overload_cast< + GooseFEM::MatrixPartitionedTyings&, const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>( - &GooseFEM::MatrixPartitionedTyings<>::Solve_u), + &GooseFEM::MatrixPartitionedTyingsSolver<>::Solve_u), "Solve_u", + py::arg("matrix"), py::arg("b_u"), py::arg("b_d"), py::arg("x_p")) .def("__repr__", - [](const GooseFEM::MatrixPartitionedTyings<>&){ - return ""; }); + [](const GooseFEM::MatrixPartitionedTyingsSolver<>&){ + return ""; }); } // =================================================================================================