diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp b/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp index 2aae832..d430c18 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp @@ -1,138 +1,138 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5, 5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLft = mesh.nodesLeftEdge(); xt::xtensor nodesRgt = mesh.nodesRightEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBot = mesh.nodesBottomEdge(); // fixed displacements DOFs // ------------------------ xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesRgt), 0), xt::view(dofs, xt::keep(nodesTop), 1), xt::view(dofs, xt::keep(nodesLft), 0), xt::view(dofs, xt::keep(nodesBot), 1))); // simulation variables // -------------------- // vector definition GooseFEM::VectorPartitioned vector(conn, dofs, iip); // allocate system matrix GooseFEM::MatrixPartitioned K(conn, dofs, iip); GooseFEM::MatrixPartitionedSolver<> Solver; // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne * ndim, nne * ndim}); // element/material definition // --------------------------- // element definition GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material definition GMatElastic::Cartesian3d::Matrix mat(nelem, nip, 1.0, 1.0); // integration point tensors xt::xtensor Eps = xt::empty({nelem, nip, 3ul, 3ul}); xt::xtensor Sig = xt::empty({nelem, nip, 3ul, 3ul}); xt::xtensor C = xt::empty({nelem, nip, 3ul, 3ul, 3ul, 3ul}); // solve // ----- // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // set fixed displacements xt::view(disp, xt::keep(nodesRgt), 0) = +0.1; xt::view(disp, xt::keep(nodesTop), 1) = -0.1; xt::view(disp, xt::keep(nodesLft), 0) = 0.0; xt::view(disp, xt::keep(nodesBot), 1) = 0.0; // residual xt::noalias(fres) = fext - fint; // solve Solver.solve(K, fres, disp); // post-process // ------------ // compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.copy_p(fint, fext); // residual xt::noalias(fres) = fext - fint; // print residual std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl; // average stress per node - xt::xtensor dV = elem.DV(2); + xt::xtensor dV = elem.AsTensor<2>(elem.dV()); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output H5Easy::File file("output.h5", H5Easy::File::Overwrite); H5Easy::dump(file, "/coor", coor); H5Easy::dump(file, "/conn", conn); H5Easy::dump(file, "/disp", disp); H5Easy::dump(file, "/Sig", SigAv); return 0; } diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp index b3b619b..130b44a 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp @@ -1,153 +1,153 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5, 5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLft = mesh.nodesLeftEdge(); xt::xtensor nodesRgt = mesh.nodesRightEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBot = mesh.nodesBottomEdge(); // fixed displacements DOFs // ------------------------ xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesRgt), 0), xt::view(dofs, xt::keep(nodesTop), 1), xt::view(dofs, xt::keep(nodesLft), 0), xt::view(dofs, xt::keep(nodesBot), 1))); // simulation variables // -------------------- // vector definition GooseFEM::VectorPartitioned vector(conn, dofs, iip); // allocate system matrix GooseFEM::MatrixPartitioned K(conn, dofs, iip); GooseFEM::MatrixPartitionedSolver<> Solver; // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // DOF values xt::xtensor u_u = xt::zeros({vector.nnu()}); xt::xtensor fres_u = xt::zeros({vector.nnu()}); xt::xtensor fext_p = xt::zeros({vector.nnp()}); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne * ndim, nne * ndim}); // element/material definition // --------------------------- // element definition GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material definition GMatElastic::Cartesian3d::Matrix mat(nelem, nip, 1.0, 1.0); // integration point tensors xt::xtensor Eps = xt::empty({nelem, nip, 3ul, 3ul}); xt::xtensor Sig = xt::empty({nelem, nip, 3ul, 3ul}); xt::xtensor C = xt::empty({nelem, nip, 3ul, 3ul, 3ul, 3ul}); // solve // ----- // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // set fixed displacements xt::xtensor u_p = xt::concatenate(xt::xtuple( +0.1 * xt::ones({nodesRgt.size()}), -0.1 * xt::ones({nodesTop.size()}), xt::zeros({nodesLft.size()}), xt::zeros({nodesBot.size()}))); // residual xt::noalias(fres) = fext - fint; // partition vector.asDofs_u(fres, fres_u); // solve Solver.solve_u(K, fres_u, u_p, u_u); // assemble to nodal vector vector.asNode(u_u, u_p, disp); // post-process // ------------ // compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.asDofs_p(fint, fext_p); // residual xt::noalias(fres) = fext - fint; // partition vector.asDofs_u(fres, fres_u); // print residual std::cout << xt::sum(xt::abs(fres_u))[0] / xt::sum(xt::abs(fext_p))[0] << std::endl; // average stress per node - xt::xtensor dV = elem.DV(2); + xt::xtensor dV = elem.AsTensor<2>(elem.dV()); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output H5Easy::File file("output.h5", H5Easy::File::Overwrite); H5Easy::dump(file, "/coor", coor); H5Easy::dump(file, "/conn", conn); H5Easy::dump(file, "/disp", disp); H5Easy::dump(file, "/Sig", SigAv); return 0; } diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp b/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp index 77e79e1..2d5356f 100644 --- a/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp +++ b/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp @@ -1,146 +1,146 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5, 5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLft = mesh.nodesLeftOpenEdge(); xt::xtensor nodesRgt = mesh.nodesRightOpenEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBot = mesh.nodesBottomEdge(); // periodicity and fixed displacements DOFs // ---------------------------------------- for (size_t j = 0; j < coor.shape(1); ++j) { xt::view(dofs, xt::keep(nodesRgt), j) = xt::view(dofs, xt::keep(nodesLft), j); } dofs = GooseFEM::Mesh::renumber(dofs); xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesBot), 0), xt::view(dofs, xt::keep(nodesBot), 1), xt::view(dofs, xt::keep(nodesTop), 0), xt::view(dofs, xt::keep(nodesTop), 1))); // simulation variables // -------------------- // vector definition GooseFEM::VectorPartitioned vector(conn, dofs, iip); // allocate system matrix GooseFEM::MatrixPartitioned K(conn, dofs, iip); GooseFEM::MatrixPartitionedSolver<> Solver; // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne * ndim, nne * ndim}); // element/material definition // --------------------------- // element definition GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material definition GMatElastic::Cartesian3d::Matrix mat(nelem, nip); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(0, 1, 5, 6), xt::all()) = 1; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setElastic(Isoft, 10.0, 1.0); mat.setElastic(Ihard, 10.0, 10.0); // integration point tensors xt::xtensor Eps = xt::empty({nelem, nip, 3ul, 3ul}); xt::xtensor Sig = xt::empty({nelem, nip, 3ul, 3ul}); xt::xtensor C = xt::empty({nelem, nip, 3ul, 3ul, 3ul, 3ul}); // solve // ----- // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // set fixed displacements xt::view(disp, xt::keep(nodesTop), 0) = +0.1; // residual xt::noalias(fres) = fext - fint; // solve Solver.solve(K, fres, disp); // post-process // ------------ // compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.copy_p(fint, fext); // residual xt::noalias(fres) = fext - fint; // print residual std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl; // average stress per node - xt::xtensor dV = elem.DV(2); + xt::xtensor dV = elem.AsTensor<2>(elem.dV()); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output H5Easy::File file("output.h5", H5Easy::File::Overwrite); H5Easy::dump(file, "/coor", coor); H5Easy::dump(file, "/conn", conn); H5Easy::dump(file, "/disp", disp); H5Easy::dump(file, "/Sig", SigAv); return 0; } diff --git a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp index 2662f9a..26ec70c 100644 --- a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp +++ b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp @@ -1,229 +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, 5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); 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); + xt::xtensor dV = elem.AsTensor<2>(elem.dV()); // 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), xt::range(0,2))); 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::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 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 932b706..17c87df 100644 --- a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp +++ b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp @@ -1,239 +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, 5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); 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), xt::range(0,2))); 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::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 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); + xt::xtensor dV = elem.AsTensor<2>(elem.dV()); // - 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 702e3a0..cb51d1d 100644 --- a/docs/examples/statics/Periodic_LinearElastic/main.cpp +++ b/docs/examples/statics/Periodic_LinearElastic/main.cpp @@ -1,171 +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::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 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); + xt::xtensor dV = elem.AsTensor<2>(elem.dV()); // - 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 4d84eb8..10d420f 100644 --- a/docs/examples/statics/Periodic_NonLinearElastic/main.cpp +++ b/docs/examples/statics/Periodic_NonLinearElastic/main.cpp @@ -1,215 +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, 5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); 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), xt::range(0,2))); 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::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 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); + xt::xtensor dV = elem.AsTensor<2>(elem.dV()); // - 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; }