diff --git a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp index 6f27648..3736232 100644 --- a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp +++ b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp @@ -1,232 +1,232 @@ #include #include #include #include #include #include #include namespace GM = GMatElastoPlastic::Cartesian3d; namespace GF = GooseFEM; namespace PV = XDMFWrite_HighFive; 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.elementgrid(); // 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.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::Array<2> mat({nelem, nip}); size_t tdim = 3; // 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.setStrain(Eps); mat.stress(Sig); mat.tangent(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.setStrain(Eps); mat.stress(Sig); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress - auto Sigbar = xt::average(Sig, dV, {0, 1}); - auto Epsbar = xt::average(Eps, dV, {0, 1}); + xt::xtensor Sigbar = xt::average(Sig, dV, {0, 1}); + xt::xtensor 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), GF::as3d(disp)); // - update ParaView meta-data // xdmf.push_back({ // PV::Topology(file, "/conn", mesh.getElementType()), // PV::Geometry(file, "/coor"), // PV::Attribute(file, "/disp/" + std::to_string(inc), PV::AttributeCenter::Node, "Displacement"), // PV::Attribute(file, "/sigeq/" + std::to_string(inc), PV::AttributeCenter::Cell, "Eq. stress"), // PV::Attribute(file, "/epseq/" + std::to_string(inc), PV::AttributeCenter::Cell, "Eq. strain")}); } // write ParaView meta-data PV::write("main.xdmf", xdmf.get()); return 0; }