diff --git a/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp b/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp index 46c2df3..fe3ff87 100644 --- a/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp +++ b/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp @@ -1,184 +1,180 @@ -#include #include +#include #include -// ------------------------------------------------------------------------------------------------- - namespace GM = GMatElastoPlasticQPot::Cartesian2d; namespace GF = GooseFEM; namespace QD = GooseFEM::Element::Quad4; namespace H5 = H5Easy; -// ------------------------------------------------------------------------------------------------- - -inline double sqdot(const xt::xtensor &M, const xt::xtensor &V) +inline double sqdot(const xt::xtensor& M, const xt::xtensor& V) { - double out = 0.; + double ret = 0.; - for (size_t i = 0; i < M.size(); ++i) - out += M(i) * V(i) * V(i); + for (size_t i = 0; i < M.size(); ++i) { + ret += M(i) * V(i) * V(i); + } - return out; + return ret; } -// ------------------------------------------------------------------------------------------------- - int main() { - // simulation parameters + // simulation parameters - double T = 60. ; // total time - double dt = 1.e-2; // time increment - size_t nx = 60 ; // number of elements in both directions - double gamma = .05 ; // displacement step + double T = 60.0; // total time + double dt = 1.0e-2; // time increment + size_t nx = 60; // number of elements in both directions + double gamma = 0.05; // displacement step - // get mesh & quadrature + // get mesh & quadrature - GF::Mesh::Quad4::Regular mesh(nx,nx,1.); + GF::Mesh::Quad4::Regular mesh(nx, nx, 1.0); - xt::xtensor coor = mesh.coor(); - xt::xtensor conn = mesh.conn(); - xt::xtensor dofs = mesh.dofsPeriodic(); + xt::xtensor coor = mesh.coor(); + xt::xtensor conn = mesh.conn(); + xt::xtensor dofs = mesh.dofsPeriodic(); - GF::Vector vector(conn, dofs); + GF::Vector vector(conn, dofs); - QD::Quadrature quad(vector.AsElement(coor)); + QD::Quadrature quad(vector.AsElement(coor)); - size_t ndim = mesh.ndim(); - size_t nne = mesh.nne(); - size_t nelem = mesh.nelem(); - size_t nip = quad.nip(); + size_t ndim = mesh.ndim(); + size_t nne = mesh.nne(); + size_t nelem = mesh.nelem(); + size_t nip = quad.nip(); - xt::xtensor u = xt::zeros(coor.shape()); - xt::xtensor v = xt::zeros(coor.shape()); - xt::xtensor a = xt::zeros(coor.shape()); - xt::xtensor v_n = xt::zeros(coor.shape()); - xt::xtensor a_n = 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()); + xt::xtensor u = xt::zeros(coor.shape()); + xt::xtensor v = xt::zeros(coor.shape()); + xt::xtensor a = xt::zeros(coor.shape()); + xt::xtensor v_n = xt::zeros(coor.shape()); + xt::xtensor a_n = 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()); - xt::xtensor ue = xt::zeros({nelem, nne, ndim}); - xt::xtensor fe = xt::zeros({nelem, nne, ndim}); + xt::xtensor ue = xt::zeros({nelem, nne, ndim}); + xt::xtensor fe = xt::zeros({nelem, nne, ndim}); - xt::xtensor Eps = xt::zeros({nelem, nip, ndim, ndim}); - xt::xtensor Sig = xt::zeros({nelem, nip, ndim, ndim}); + xt::xtensor Eps = xt::zeros({nelem, nip, ndim, ndim}); + xt::xtensor Sig = xt::zeros({nelem, nip, ndim, ndim}); - // material definition + // material definition - GM::Matrix material({nelem, nip}); + GM::Matrix material({nelem, nip}); - xt::xtensor Ihard = xt::zeros({nelem, nip}); - xt::xtensor Isoft = xt::ones ({nelem, nip}); - xt::view(Ihard, xt::range(0, nx * nx / 4), xt::all()) = 1ul; - Isoft -= Ihard; + xt::xtensor Ihard = xt::zeros({nelem, nip}); + xt::xtensor Isoft = xt::ones({nelem, nip}); + xt::view(Ihard, xt::range(0, nx * nx / 4), xt::all()) = 1ul; + Isoft -= Ihard; - material.setElastic(Ihard, 100., 10.); - material.setElastic(Isoft, 100., 1.); + material.setElastic(Ihard, 100.0, 10.0); + material.setElastic(Isoft, 100.0, 1.0); - material.check(); + material.check(); - // mass matrix + // mass matrix - QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w()); + QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w()); - xt::xtensor rho = 1.0 * xt::ones({nelem, nip}); + xt::xtensor rho = 1.0 * xt::ones({nelem, nip}); - GF::MatrixDiagonal M(conn, dofs); + GF::MatrixDiagonal M(conn, dofs); - M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho)); + M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho)); - xt::xtensor mass = M.AsDiagonal(); + xt::xtensor mass = M.AsDiagonal(); - // update in macroscopic deformation gradient + // update in macroscopic deformation gradient - xt::xtensor dFbar = xt::zeros({2,2}); + xt::xtensor dFbar = xt::zeros({2, 2}); - dFbar(0,1) = gamma; + dFbar(0, 1) = gamma; - for (size_t j = 0; j < ndim; ++j ) - for (size_t k = 0; k < ndim; ++k ) - xt::view(u, xt::all(), j) += dFbar(j, k) * (xt::view(coor, xt::all(), k) - coor(0, k)); + for (size_t j = 0; j < ndim; ++j) { + for (size_t k = 0; k < ndim; ++k) { + xt::view(u, xt::all(), j) += dFbar(j, k) * (xt::view(coor, xt::all(), k) - coor(0, k)); + } + } - // output variables + // output variables - xt::xtensor Epot = xt::zeros({static_cast(T/dt)}); - xt::xtensor Ekin = xt::zeros({static_cast(T/dt)}); - xt::xtensor t = xt::zeros({static_cast(T/dt)}); + xt::xtensor Epot = xt::zeros({static_cast(T / dt)}); + xt::xtensor Ekin = xt::zeros({static_cast(T / dt)}); + xt::xtensor t = xt::zeros({static_cast(T / dt)}); - xt::xtensor dV = quad.DV(); + xt::xtensor dV = quad.DV(); - // loop over increments + // loop over increments - for (size_t inc = 0; inc < static_cast(Epot.size()); ++inc) - { - // store history + for (size_t inc = 0; inc < static_cast(Epot.size()); ++inc) { - xt::noalias(v_n) = v; - xt::noalias(a_n) = a; + // store history - // new displacement + xt::noalias(v_n) = v; + xt::noalias(a_n) = a; - xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt,2.) * a; + // new displacement - // compute strain/strain, and corresponding internal force + xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt, 2.0) * a; - vector.asElement(u, ue); - quad.symGradN_vector(ue, Eps); - material.stress(Eps, Sig); - quad.int_gradN_dot_tensor2_dV(Sig, fe); - vector.assembleNode(fe, fint); + // compute strain/strain, and corresponding internal force - // estimate new velocity + vector.asElement(u, ue); + quad.symGradN_vector(ue, Eps); + material.stress(Eps, Sig); + quad.int_gradN_dot_tensor2_dV(Sig, fe); + vector.assembleNode(fe, fint); - xt::noalias(v) = v_n + dt * a_n; + // estimate new velocity - // compute residual force & solve + xt::noalias(v) = v_n + dt * a_n; - xt::noalias(fres) = fext - fint; + // compute residual force & solve - M.solve(fres, a); + xt::noalias(fres) = fext - fint; - // re-estimate new velocity + M.solve(fres, a); - xt::noalias(v) = v_n + .5 * dt * (a_n + a); + // re-estimate new velocity - // compute residual force & solve + xt::noalias(v) = v_n + 0.5 * dt * (a_n + a); - xt::noalias(fres) = fext - fint; + // compute residual force & solve - M.solve(fres, a); + xt::noalias(fres) = fext - fint; - // new velocity + M.solve(fres, a); - xt::noalias(v) = v_n + .5 * dt * (a_n + a); + // new velocity - // compute residual force & solve + xt::noalias(v) = v_n + 0.5 * dt * (a_n + a); - xt::noalias(fres) = fext - fint; + // compute residual force & solve - M.solve(fres, a); + xt::noalias(fres) = fext - fint; - // store output variables + M.solve(fres, a); - xt::xtensor E = material.Energy(Eps); - xt::xtensor V = vector.AsDofs(v); + // store output variables - t (inc) = static_cast(inc) * dt; - Ekin(inc) = 0.5 * sqdot(mass,V); - Epot(inc) = xt::sum(E*dV)[0]; + xt::xtensor E = material.Energy(Eps); + xt::xtensor V = vector.AsDofs(v); - } + t(inc) = static_cast(inc) * dt; + Ekin(inc) = 0.5 * sqdot(mass, V); + Epot(inc) = xt::sum(E * dV)[0]; + } - // write output variables to file - H5::File file("example.hdf5", H5::File::Overwrite); - H5::dump(file, "/global/Epot", Epot); - H5::dump(file, "/global/Ekin", Ekin); - H5::dump(file, "/global/t" , t ); - H5::dump(file, "/mesh/conn" , conn); - H5::dump(file, "/mesh/coor" , coor); - H5::dump(file, "/mesh/disp" , u ); + // write output variables to file + H5::File file("example.hdf5", H5::File::Overwrite); + H5::dump(file, "/global/Epot", Epot); + H5::dump(file, "/global/Ekin", Ekin); + H5::dump(file, "/global/t", t); + H5::dump(file, "/mesh/conn", conn); + H5::dump(file, "/mesh/coor", coor); + H5::dump(file, "/mesh/disp", u); - return 0; + return 0; } diff --git a/docs/examples/dynamics/Elastic-VelocityVerlet/plot.py b/docs/examples/dynamics/Elastic-VelocityVerlet/plot.py index 5a365ea..5e1320f 100644 --- a/docs/examples/dynamics/Elastic-VelocityVerlet/plot.py +++ b/docs/examples/dynamics/Elastic-VelocityVerlet/plot.py @@ -1,49 +1,49 @@ import h5py import numpy as np import matplotlib.pyplot as plt import GooseMPL as gplt plt.style.use(['goose', 'goose-latex']) # -------------------------------------------------------------------------------------------------- with h5py.File('example.hdf5', 'r') as data: - t = data['/global/t'][...] + t = data['/global/t'][...] - Epot = data['/global/Epot'][...] - Ekin = data['/global/Ekin'][...] + Epot = data['/global/Epot'][...] + Ekin = data['/global/Ekin'][...] - coor = data['/mesh/coor'][...] - conn = data['/mesh/conn'][...] - disp = data['/mesh/disp'][...] + coor = data['/mesh/coor'][...] + conn = data['/mesh/conn'][...] + disp = data['/mesh/disp'][...] # -------------------------------------------------------------------------------------------------- fig, axes = gplt.subplots(ncols=2) ax = axes[0] ax.plot(t, Ekin + Epot, color='k', label=r'$E_\mathrm{pot} + E_\mathrm{kin}$') ax.plot(t, Epot, color='r', label=r'$E_\mathrm{pot}$') ax.plot(t, Ekin, color='b', label=r'$E_\mathrm{kin}$') ax.legend(loc='lower right') ax.set_xlabel(r'$t$') ax.set_ylabel(r'$E$') ax = axes[1] -im = gplt.patch(coor=coor+disp, conn=conn) +im = gplt.patch(coor=coor + disp, conn=conn) -Lx = np.max(coor[:,0])-np.min(coor[:,0]) -Ly = np.max(coor[:,1])-np.min(coor[:,1]) +Lx = np.max(coor[:, 0]) - np.min(coor[:, 0]) +Ly = np.max(coor[:, 1]) - np.min(coor[:, 1]) -ax.set_xlim([-0.1*Lx, 1.3*Lx]) -ax.set_ylim([-0.1*Ly, 1.3*Ly]) +ax.set_xlim([-0.1 * Lx, 1.3 * Lx]) +ax.set_ylim([-0.1 * Ly, 1.3 * Ly]) ax.set_aspect('equal') plt.show() diff --git a/docs/examples/dynamics/Elastic-Verlet/main.cpp b/docs/examples/dynamics/Elastic-Verlet/main.cpp index 190bca7..6f40094 100644 --- a/docs/examples/dynamics/Elastic-Verlet/main.cpp +++ b/docs/examples/dynamics/Elastic-Verlet/main.cpp @@ -1,164 +1,160 @@ -#include #include +#include #include -// ------------------------------------------------------------------------------------------------- - namespace GM = GMatElastoPlasticQPot::Cartesian2d; namespace GF = GooseFEM; namespace QD = GooseFEM::Element::Quad4; namespace H5 = H5Easy; -// ------------------------------------------------------------------------------------------------- - -inline double sqdot(const xt::xtensor &M, const xt::xtensor &V) +inline double sqdot(const xt::xtensor& M, const xt::xtensor& V) { - double out = 0.; + double ret = 0.0; - for (size_t i = 0; i < M.size(); ++i) - out += M(i) * V(i) * V(i); + for (size_t i = 0; i < M.size(); ++i) { + ret += M(i) * V(i) * V(i); + } - return out; + return ret; } -// ------------------------------------------------------------------------------------------------- - int main() { - // simulation parameters + // simulation parameters - double T = 60. ; // total time - double dt = 1.e-2; // time increment - size_t nx = 60 ; // number of elements in both directions - double gamma = .05 ; // displacement step + double T = 60.0; // total time + double dt = 1.0e-2; // time increment + size_t nx = 60; // number of elements in both directions + double gamma = 0.05; // displacement step - // get mesh & quadrature + // get mesh & quadrature - GF::Mesh::Quad4::Regular mesh(nx,nx,1.); + GF::Mesh::Quad4::Regular mesh(nx, nx, 1.0); - xt::xtensor coor = mesh.coor(); - xt::xtensor conn = mesh.conn(); - xt::xtensor dofs = mesh.dofsPeriodic(); + xt::xtensor coor = mesh.coor(); + xt::xtensor conn = mesh.conn(); + xt::xtensor dofs = mesh.dofsPeriodic(); - GF::Vector vector(conn, dofs); + GF::Vector vector(conn, dofs); - QD::Quadrature quad(vector.AsElement(coor)); + QD::Quadrature quad(vector.AsElement(coor)); - size_t ndim = mesh.ndim(); - size_t nne = mesh.nne(); - size_t nelem = mesh.nelem(); - size_t nip = quad.nip(); + size_t ndim = mesh.ndim(); + size_t nne = mesh.nne(); + size_t nelem = mesh.nelem(); + size_t nip = quad.nip(); - xt::xtensor u = xt::zeros(coor.shape()); - xt::xtensor v = xt::zeros(coor.shape()); - xt::xtensor a = xt::zeros(coor.shape()); - xt::xtensor v_n = xt::zeros(coor.shape()); - xt::xtensor a_n = 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()); + xt::xtensor u = xt::zeros(coor.shape()); + xt::xtensor v = xt::zeros(coor.shape()); + xt::xtensor a = xt::zeros(coor.shape()); + xt::xtensor v_n = xt::zeros(coor.shape()); + xt::xtensor a_n = 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()); - xt::xtensor ue = xt::zeros({nelem, nne, ndim}); - xt::xtensor fe = xt::zeros({nelem, nne, ndim}); + xt::xtensor ue = xt::zeros({nelem, nne, ndim}); + xt::xtensor fe = xt::zeros({nelem, nne, ndim}); - xt::xtensor Eps = xt::zeros({nelem, nip, ndim, ndim}); - xt::xtensor Sig = xt::zeros({nelem, nip, ndim, ndim}); + xt::xtensor Eps = xt::zeros({nelem, nip, ndim, ndim}); + xt::xtensor Sig = xt::zeros({nelem, nip, ndim, ndim}); - // material definition + // material definition - GM::Matrix material({nelem, nip}); + GM::Matrix material({nelem, nip}); - xt::xtensor Ihard = xt::zeros({nelem, nip}); - xt::xtensor Isoft = xt::ones ({nelem, nip}); - xt::view(Ihard, xt::range(0, nx * nx / 4), xt::all()) = 1ul; - Isoft -= Ihard; + xt::xtensor Ihard = xt::zeros({nelem, nip}); + xt::xtensor Isoft = xt::ones({nelem, nip}); + xt::view(Ihard, xt::range(0, nx * nx / 4), xt::all()) = 1ul; + Isoft -= Ihard; - material.setElastic(Ihard, 100., 10.); - material.setElastic(Isoft, 100., 1.); + material.setElastic(Ihard, 100.0, 10.0); + material.setElastic(Isoft, 100.0, 1.0); - material.check(); + material.check(); - // mass matrix + // mass matrix - QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w()); + QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w()); - xt::xtensor rho = 1.0 * xt::ones({nelem, nip}); + xt::xtensor rho = 1.0 * xt::ones({nelem, nip}); - GF::MatrixDiagonal M(conn, dofs); + GF::MatrixDiagonal M(conn, dofs); - M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho)); + M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho)); - xt::xtensor mass = M.AsDiagonal(); + xt::xtensor mass = M.AsDiagonal(); - // update in macroscopic deformation gradient + // update in macroscopic deformation gradient - xt::xtensor dFbar = xt::zeros({2,2}); + xt::xtensor dFbar = xt::zeros({2, 2}); - dFbar(0,1) = gamma; + dFbar(0, 1) = gamma; - for (size_t j = 0; j < ndim; ++j ) - for (size_t k = 0; k < ndim; ++k ) - xt::view(u, xt::all(), j) += dFbar(j, k) * (xt::view(coor, xt::all(), k) - coor(0, k)); + for (size_t j = 0; j < ndim; ++j) { + for (size_t k = 0; k < ndim; ++k) { + xt::view(u, xt::all(), j) += dFbar(j, k) * (xt::view(coor, xt::all(), k) - coor(0, k)); + } + } - // output variables + // output variables - xt::xtensor Epot = xt::zeros({static_cast(T/dt)}); - xt::xtensor Ekin = xt::zeros({static_cast(T/dt)}); - xt::xtensor t = xt::zeros({static_cast(T/dt)}); + xt::xtensor Epot = xt::zeros({static_cast(T / dt)}); + xt::xtensor Ekin = xt::zeros({static_cast(T / dt)}); + xt::xtensor t = xt::zeros({static_cast(T / dt)}); - xt::xtensor dV = quad.DV(); + xt::xtensor dV = quad.DV(); - // loop over increments + // loop over increments - for (size_t inc = 0; inc < static_cast(Epot.size()); ++inc) - { - // store history + for (size_t inc = 0; inc < static_cast(Epot.size()); ++inc) { - xt::noalias(v_n) = v; - xt::noalias(a_n) = a; + // store history - // new displacement + xt::noalias(v_n) = v; + xt::noalias(a_n) = a; - xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt,2.) * a; + // new displacement - // compute strain/strain, and corresponding internal force + xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt, 2.0) * a; - vector.asElement(u, ue); - quad.symGradN_vector(ue, Eps); - material.stress(Eps, Sig); - quad.int_gradN_dot_tensor2_dV(Sig, fe); - vector.assembleNode(fe, fint); + // compute strain/strain, and corresponding internal force - // compute residual force & solve + vector.asElement(u, ue); + quad.symGradN_vector(ue, Eps); + material.stress(Eps, Sig); + quad.int_gradN_dot_tensor2_dV(Sig, fe); + vector.assembleNode(fe, fint); - xt::noalias(fres) = fext - fint; + // compute residual force & solve - M.solve(fres, a); + xt::noalias(fres) = fext - fint; - // new velocity + M.solve(fres, a); - xt::noalias(v) = v_n + .5 * dt * (a_n + a); + // new velocity - // store output variables + xt::noalias(v) = v_n + 0.5 * dt * (a_n + a); - xt::xtensor E = material.Energy(Eps); - xt::xtensor V = vector.AsDofs(v); + // store output variables - t (inc) = static_cast(inc) * dt; - Ekin(inc) = 0.5 * sqdot(mass,V); - Epot(inc) = xt::sum(E*dV)[0]; + xt::xtensor E = material.Energy(Eps); + xt::xtensor V = vector.AsDofs(v); - } + t(inc) = static_cast(inc) * dt; + Ekin(inc) = 0.5 * sqdot(mass, V); + Epot(inc) = xt::sum(E * dV)[0]; + } - // write output variables to file - H5::File file("example.hdf5", H5::File::Overwrite); - H5::dump(file, "/global/Epot", Epot); - H5::dump(file, "/global/Ekin", Ekin); - H5::dump(file, "/global/t" , t ); - H5::dump(file, "/mesh/conn" , conn); - H5::dump(file, "/mesh/coor" , coor); - H5::dump(file, "/mesh/disp" , u ); + // write output variables to file + H5::File file("example.hdf5", H5::File::Overwrite); + H5::dump(file, "/global/Epot", Epot); + H5::dump(file, "/global/Ekin", Ekin); + H5::dump(file, "/global/t", t); + H5::dump(file, "/mesh/conn", conn); + H5::dump(file, "/mesh/coor", coor); + H5::dump(file, "/mesh/disp", u); - return 0; + return 0; } diff --git a/docs/examples/dynamics/Elastic-Verlet/plot.py b/docs/examples/dynamics/Elastic-Verlet/plot.py index 5a365ea..5e1320f 100644 --- a/docs/examples/dynamics/Elastic-Verlet/plot.py +++ b/docs/examples/dynamics/Elastic-Verlet/plot.py @@ -1,49 +1,49 @@ import h5py import numpy as np import matplotlib.pyplot as plt import GooseMPL as gplt plt.style.use(['goose', 'goose-latex']) # -------------------------------------------------------------------------------------------------- with h5py.File('example.hdf5', 'r') as data: - t = data['/global/t'][...] + t = data['/global/t'][...] - Epot = data['/global/Epot'][...] - Ekin = data['/global/Ekin'][...] + Epot = data['/global/Epot'][...] + Ekin = data['/global/Ekin'][...] - coor = data['/mesh/coor'][...] - conn = data['/mesh/conn'][...] - disp = data['/mesh/disp'][...] + coor = data['/mesh/coor'][...] + conn = data['/mesh/conn'][...] + disp = data['/mesh/disp'][...] # -------------------------------------------------------------------------------------------------- fig, axes = gplt.subplots(ncols=2) ax = axes[0] ax.plot(t, Ekin + Epot, color='k', label=r'$E_\mathrm{pot} + E_\mathrm{kin}$') ax.plot(t, Epot, color='r', label=r'$E_\mathrm{pot}$') ax.plot(t, Ekin, color='b', label=r'$E_\mathrm{kin}$') ax.legend(loc='lower right') ax.set_xlabel(r'$t$') ax.set_ylabel(r'$E$') ax = axes[1] -im = gplt.patch(coor=coor+disp, conn=conn) +im = gplt.patch(coor=coor + disp, conn=conn) -Lx = np.max(coor[:,0])-np.min(coor[:,0]) -Ly = np.max(coor[:,1])-np.min(coor[:,1]) +Lx = np.max(coor[:, 0]) - np.min(coor[:, 0]) +Ly = np.max(coor[:, 1]) - np.min(coor[:, 1]) -ax.set_xlim([-0.1*Lx, 1.3*Lx]) -ax.set_ylim([-0.1*Ly, 1.3*Ly]) +ax.set_xlim([-0.1 * Lx, 1.3 * Lx]) +ax.set_ylim([-0.1 * Ly, 1.3 * Ly]) ax.set_aspect('equal') plt.show() diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp b/docs/examples/statics/FixedDisplacements_LinearElastic/example.cpp index d430c18..ff164c4 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 #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(); + 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(); + 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::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()); + 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}); + 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}); + 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.AsTensor<2>(elem.dV()); - xt::xtensor SigAv = xt::average(Sig, dV, {1}); + 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/example.py b/docs/examples/statics/FixedDisplacements_LinearElastic/example.py index 1fcb97a..aacf7fa 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/example.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/example.py @@ -1,182 +1,180 @@ import sys import GooseFEM import GMatElastic import numpy as np # mesh # ---- # define mesh mesh = GooseFEM.Mesh.Quad4.Regular(5, 5) # mesh dimensions nelem = mesh.nelem() nne = mesh.nne() ndim = mesh.ndim() # mesh definitions coor = mesh.coor() conn = mesh.conn() dofs = mesh.dofs() # node sets nodesLft = mesh.nodesLeftEdge() nodesRgt = mesh.nodesRightEdge() nodesTop = mesh.nodesTopEdge() nodesBot = mesh.nodesBottomEdge() # fixed displacements DOFs # ------------------------ iip = np.concatenate(( dofs[nodesRgt, 0], dofs[nodesTop, 1], dofs[nodesLft, 0], dofs[nodesBot, 1] )) # simulation variables # -------------------- # vector definition vector = GooseFEM.VectorPartitioned(conn, dofs, iip) # allocate system matrix K = GooseFEM.MatrixPartitioned(conn, dofs, iip) Solver = GooseFEM.MatrixPartitionedSolver() # nodal quantities disp = np.zeros(coor.shape) fint = np.zeros(coor.shape) fext = np.zeros(coor.shape) fres = np.zeros(coor.shape) # element vectors ue = np.empty((nelem, nne, ndim)) fe = np.empty((nelem, nne, ndim)) Ke = np.empty((nelem, nne * ndim, nne * ndim)) # element/material definition # --------------------------- # element definition elem = GooseFEM.Element.Quad4.QuadraturePlanar(vector.AsElement(coor)) nip = elem.nip() # material definition mat = GMatElastic.Cartesian3d.Matrix(nelem, nip, 1.0, 1.0) # integration point tensors Eps = np.empty((nelem, nip, 3, 3)) Sig = np.empty((nelem, nip, 3, 3)) C = np.empty((nelem, nip, 3, 3, 3, 3)) # solve # ----- # strain ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) # stress & tangent (Sig, C) = mat.Tangent(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # stiffness matrix Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C) K.assemble(Ke) # set fixed displacements disp[nodesRgt, 0] = +0.1 disp[nodesTop, 1] = -0.1 disp[nodesLft, 0] = 0.0 disp[nodesBot, 1] = 0.0 # residual fres = fext - fint # solve disp = Solver.Solve(K, fres, disp) # post-process # ------------ # compute strain and stress ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) Sig = mat.Stress(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # apply reaction force fext = vector.Copy_p(fint, fext) # residual fres = fext - fint # print residual print(np.sum(np.abs(fres)) / np.sum(np.abs(fext))) # average stress per node dV = elem.DV(2) Sig = np.average(Sig, weights=dV, axis=1) # skip plot with "--no-plot" command line argument # ------------------------------------------------ if len(sys.argv) == 2: if sys.argv[1] == "--no-plot": sys.exit(0) # plot # ---- -import matplotlib.pyplot as plt import GooseMPL as gplt +import matplotlib.pyplot as plt plt.style.use(['goose', 'goose-latex']) # extract dimension nelem = conn.shape[0] # tensor products def ddot22(A2, B2): return np.einsum('eij, eji -> e', A2, B2) - def ddot42(A4, B2): return np.einsum('eijkl, elk -> eij', A4, B2) - def dyad22(A2, B2): return np.einsum('eij, ekl -> eijkl', A2, B2) # identity tensors i = np.eye(3) I = np.einsum('ij, e', i, np.ones([nelem])) I4 = np.einsum('ijkl, e -> eijkl', np.einsum('il, jk', i, i), np.ones([nelem])) I4rt = np.einsum('ijkl, e -> eijkl', np.einsum('ik,jl', i, i), np.ones([nelem])) I4s = 0.5 * (I4 + I4rt) II = dyad22(I, I) I4d = I4s - II / 3.0 # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3.0 / 2.0 * ddot22(Sigd, Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor + disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0, 0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) plt.savefig('plot.pdf') plt.close() diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.cpp index 130b44a..cbb4dfa 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 #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(); + 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(); + 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::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()); + 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()}); + 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}); + 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}); + 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( + 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.AsTensor<2>(elem.dV()); - xt::xtensor SigAv = xt::average(Sig, dV, {1}); + 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.py b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py index 46a47b5..cdf902e 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py @@ -1,198 +1,196 @@ import sys import GooseFEM import GMatElastic import numpy as np # mesh # ---- # define mesh mesh = GooseFEM.Mesh.Quad4.Regular(5, 5) # mesh dimensions nelem = mesh.nelem() nne = mesh.nne() ndim = mesh.ndim() # mesh definitions coor = mesh.coor() conn = mesh.conn() dofs = mesh.dofs() # node sets nodesLft = mesh.nodesLeftEdge() nodesRgt = mesh.nodesRightEdge() nodesTop = mesh.nodesTopEdge() nodesBot = mesh.nodesBottomEdge() # fixed displacements DOFs # ------------------------ iip = np.concatenate(( dofs[nodesRgt, 0], dofs[nodesTop, 1], dofs[nodesLft, 0], dofs[nodesBot, 1] )) # simulation variables # -------------------- # vector definition vector = GooseFEM.VectorPartitioned(conn, dofs, iip) # allocate system matrix K = GooseFEM.MatrixPartitioned(conn, dofs, iip) Solver = GooseFEM.MatrixPartitionedSolver() # nodal quantities disp = np.zeros(coor.shape) fint = np.zeros(coor.shape) fext = np.zeros(coor.shape) fres = np.zeros(coor.shape) # DOF values fext_p = np.zeros(vector.nnp()) fres_u = np.zeros(vector.nnu()) fext_p = np.zeros(vector.nnp()) # element vectors ue = np.empty((nelem, nne, ndim)) fe = np.empty((nelem, nne, ndim)) Ke = np.empty((nelem, nne * ndim, nne * ndim)) # element/material definition # --------------------------- # element definition elem = GooseFEM.Element.Quad4.QuadraturePlanar(vector.AsElement(coor)) nip = elem.nip() # material definition mat = GMatElastic.Cartesian3d.Matrix(nelem, nip, 1.0, 1.0) # integration point tensors Eps = np.empty((nelem, nip, 3, 3)) Sig = np.empty((nelem, nip, 3, 3)) C = np.empty((nelem, nip, 3, 3, 3, 3)) # solve # ----- # strain ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) # stress & tangent (Sig, C) = mat.Tangent(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # stiffness matrix Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C) K.assemble(Ke) # set fixed displacements u_p = np.concatenate(( +0.1 * np.ones(nodesRgt.size), -0.1 * np.ones(nodesTop.size), np.zeros(nodesLft.size), np.zeros(nodesBot.size) )) # residual fres = fext - fint # partition fres_u = vector.AsDofs_u(fres) # solve u_u = Solver.Solve_u(K, fres_u, u_p) # assemble to nodal vector disp = vector.AsNode(u_u, u_p) # post-process # ------------ # compute strain and stress ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) Sig = mat.Stress(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # apply reaction force fext_p = vector.AsDofs_p(fint) # residual fres = fext - fint # partition fres_u = vector.AsDofs_u(fres) # print residual print(np.sum(np.abs(fres_u)) / np.sum(np.abs(fext_p))) # average stress per node dV = elem.DV(2) Sig = np.average(Sig, weights=dV, axis=1) # skip plot with "--no-plot" command line argument # ------------------------------------------------ if len(sys.argv) == 2: if sys.argv[1] == "--no-plot": sys.exit(0) # plot # ---- -import GooseMPL as gplt import matplotlib.pyplot as plt +import GooseMPL as gplt plt.style.use(['goose', 'goose-latex']) # extract dimension nelem = conn.shape[0] # tensor products def ddot22(A2, B2): return np.einsum('eij, eji -> e', A2, B2) - def ddot42(A4, B2): return np.einsum('eijkl, elk -> eij', A4, B2) - def dyad22(A2, B2): return np.einsum('eij, ekl -> eijkl', A2, B2) # identity tensors i = np.eye(3) I = np.einsum('ij, e', i, np.ones([nelem])) I4 = np.einsum('ijkl, e -> eijkl', np.einsum('il, jk', i, i), np.ones([nelem])) I4rt = np.einsum('ijkl, e -> eijkl', np.einsum('ik,jl', i, i), np.ones([nelem])) I4s = 0.5 * (I4 + I4rt) II = dyad22(I, I) I4d = I4s - II / 3.0 # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3.0 / 2.0 * ddot22(Sigd, Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor + disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0, 0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) plt.savefig('plot.pdf') plt.close() diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/plot.py b/docs/examples/statics/FixedDisplacements_LinearElastic/plot.py index 8447c77..52fd916 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/plot.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/plot.py @@ -1,51 +1,49 @@ import h5py import matplotlib.pyplot as plt import GooseMPL as gplt import numpy as np plt.style.use(['goose', 'goose-latex']) # load data with h5py.File('output.h5', 'r') as data: coor = data['/coor'][...] conn = data['/conn'][...] disp = data['/disp'][...] Sig = data['/Sig'][...] nelem = conn.shape[0] # tensor products def ddot22(A2, B2): return np.einsum('eij, eji -> e', A2, B2) - def ddot42(A4, B2): return np.einsum('eijkl, elk -> eij', A4, B2) - def dyad22(A2, B2): return np.einsum('eij, ekl -> eijkl', A2, B2) # identity tensors i = np.eye(3) I = np.einsum('ij, e', i, np.ones([nelem])) I4 = np.einsum('ijkl, e -> eijkl', np.einsum('il, jk', i, i), np.ones([nelem])) I4rt = np.einsum('ijkl, e -> eijkl', np.einsum('ik,jl', i, i), np.ones([nelem])) I4s = 0.5 * (I4 + I4rt) II = dyad22(I, I) I4d = I4s - II / 3.0 # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3.0 / 2.0 * ddot22(Sigd, Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor + disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0, 0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) plt.savefig('plot.pdf') plt.close() diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp b/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp index 2d5356f..e8da101 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 #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(); + 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(); + 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::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()); + 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}); + 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::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; + 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}); + 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.AsTensor<2>(elem.dV()); - xt::xtensor SigAv = xt::average(Sig, dV, {1}); + 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.py b/docs/examples/statics/MixedPeriodic_LinearElastic/example.py index d4eaf42..fd10013 100644 --- a/docs/examples/statics/MixedPeriodic_LinearElastic/example.py +++ b/docs/examples/statics/MixedPeriodic_LinearElastic/example.py @@ -1,189 +1,187 @@ import sys import GooseFEM import GMatElastic import numpy as np # mesh # ---- # define mesh mesh = GooseFEM.Mesh.Quad4.Regular(5, 5) # mesh dimensions nelem = mesh.nelem() nne = mesh.nne() ndim = mesh.ndim() # mesh definitions coor = mesh.coor() conn = mesh.conn() dofs = mesh.dofs() # node sets nodesLft = mesh.nodesLeftOpenEdge() nodesRgt = mesh.nodesRightOpenEdge() nodesTop = mesh.nodesTopEdge() nodesBot = mesh.nodesBottomEdge() # periodicity and fixed displacements DOFs # ---------------------------------------- for j in range(coor.shape[1]): dofs[nodesRgt, j] = dofs[nodesLft, j] dofs = GooseFEM.Mesh.renumber(dofs) iip = np.concatenate(( dofs[nodesBot, 0], dofs[nodesBot, 1], dofs[nodesTop, 0], dofs[nodesTop, 1] )) # simulation variables # -------------------- # vector definition vector = GooseFEM.VectorPartitioned(conn, dofs, iip) # allocate system matrix K = GooseFEM.MatrixPartitioned(conn, dofs, iip) Solver = GooseFEM.MatrixPartitionedSolver() # nodal quantities disp = np.zeros(coor.shape) fint = np.zeros(coor.shape) fext = np.zeros(coor.shape) fres = np.zeros(coor.shape) # element vectors ue = np.empty((nelem, nne, ndim)) fe = np.empty((nelem, nne, ndim)) Ke = np.empty((nelem, nne * ndim, nne * ndim)) # element/material definition # --------------------------- # element definition elem = GooseFEM.Element.Quad4.QuadraturePlanar(vector.AsElement(coor)) nip = elem.nip() # material definition mat = GMatElastic.Cartesian3d.Matrix(nelem, nip) Ihard = np.zeros((nelem, nip), dtype='int') Ihard[[0, 1, 5, 6], :] = 1 Isoft = np.ones((nelem, nip), dtype='int') - Ihard mat.setElastic(Isoft, 10.0, 1.0) mat.setElastic(Ihard, 10.0, 10.0) # integration point tensors Eps = np.empty((nelem, nip, 3, 3)) Sig = np.empty((nelem, nip, 3, 3)) C = np.empty((nelem, nip, 3, 3, 3, 3)) # solve # ----- # strain ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) # stress & tangent (Sig, C) = mat.Tangent(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # stiffness matrix Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C) K.assemble(Ke) # set fixed displacements disp[nodesTop, 0] = +0.1 # residual fres = fext - fint # solve disp = Solver.Solve(K, fres, disp) # post-process # ------------ # compute strain and stress ue = vector.AsElement(disp) Eps = elem.SymGradN_vector(ue) Sig = mat.Stress(Eps) # internal force fe = elem.Int_gradN_dot_tensor2_dV(Sig) fint = vector.AssembleNode(fe) # apply reaction force fext = vector.Copy_p(fint, fext) # residual fres = fext - fint # print residual print(np.sum(np.abs(fres)) / np.sum(np.abs(fext))) # average stress per node dV = elem.DV(2) Sig = np.average(Sig, weights=dV, axis=1) # skip plot with "--no-plot" command line argument # ------------------------------------------------ if len(sys.argv) == 2: if sys.argv[1] == "--no-plot": sys.exit(0) # plot # ---- import matplotlib.pyplot as plt import GooseMPL as gplt plt.style.use(['goose', 'goose-latex']) # extract dimension nelem = conn.shape[0] # tensor products def ddot22(A2, B2): return np.einsum('eij, eji -> e', A2, B2) - def ddot42(A4, B2): return np.einsum('eijkl, elk -> eij', A4, B2) - def dyad22(A2, B2): return np.einsum('eij, ekl -> eijkl', A2, B2) # identity tensors i = np.eye(3) I = np.einsum('ij, e', i, np.ones([nelem])) I4 = np.einsum('ijkl, e -> eijkl', np.einsum('il, jk', i, i), np.ones([nelem])) I4rt = np.einsum('ijkl, e -> eijkl', np.einsum('ik,jl', i, i), np.ones([nelem])) I4s = 0.5 * (I4 + I4rt) II = dyad22(I, I) I4d = I4s - II / 3.0 # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3.0 / 2.0 * ddot22(Sigd, Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor + disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0, 0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) plt.savefig('plot.pdf') plt.close() diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py b/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py index 8447c77..52fd916 100644 --- a/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py +++ b/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py @@ -1,51 +1,49 @@ import h5py import matplotlib.pyplot as plt import GooseMPL as gplt import numpy as np plt.style.use(['goose', 'goose-latex']) # load data with h5py.File('output.h5', 'r') as data: coor = data['/coor'][...] conn = data['/conn'][...] disp = data['/disp'][...] Sig = data['/Sig'][...] nelem = conn.shape[0] # tensor products def ddot22(A2, B2): return np.einsum('eij, eji -> e', A2, B2) - def ddot42(A4, B2): return np.einsum('eijkl, elk -> eij', A4, B2) - def dyad22(A2, B2): return np.einsum('eij, ekl -> eijkl', A2, B2) # identity tensors i = np.eye(3) I = np.einsum('ij, e', i, np.ones([nelem])) I4 = np.einsum('ijkl, e -> eijkl', np.einsum('il, jk', i, i), np.ones([nelem])) I4rt = np.einsum('ijkl, e -> eijkl', np.einsum('ik,jl', i, i), np.ones([nelem])) I4s = 0.5 * (I4 + I4rt) II = dyad22(I, I) I4d = I4s - II / 3.0 # compute equivalent stress Sigd = ddot42(I4d, Sig) sigeq = np.sqrt(3.0 / 2.0 * ddot22(Sigd, Sigd)) # plot fig, ax = plt.subplots() gplt.patch(coor=coor + disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0, 0.1)) gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax) plt.savefig('plot.pdf') plt.close() diff --git a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp index 26ec70c..627a7c8 100644 --- a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp +++ b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp @@ -1,229 +1,232 @@ -#include +#include #include +#include +#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.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; + // 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.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), + })); } - // 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; + // write ParaView meta-data + xdmf.write("main.xdmf"); + + return 0; } diff --git a/docs/examples/statics/Periodic_ElastoPlastic/plot.py b/docs/examples/statics/Periodic_ElastoPlastic/plot.py index 2ddfaa0..fefa86f 100644 --- a/docs/examples/statics/Periodic_ElastoPlastic/plot.py +++ b/docs/examples/statics/Periodic_ElastoPlastic/plot.py @@ -1,33 +1,32 @@ import h5py import matplotlib.pyplot as plt -import GooseMPL as gplt -import numpy as np +import GooseMPL as gplt +import numpy as np plt.style.use(['goose', 'goose-latex']) # open file file = h5py.File('main.h5', 'r') # read stored increments incs = file['/stored'][...] # read fields -coor = file['/coor' ][...] -conn = file['/conn' ][...] -disp = file['/disp' ][str(np.max(incs))][...][:,:2] +coor = file['/coor'][...] +conn = file['/conn'][...] +disp = file['/disp'][str(np.max(incs))][...][:, :2] Sigeq = file['/sigeq'][str(np.max(incs))][...] sigeq = file['/macroscopic/sigeq'][...] epseq = file['/macroscopic/epseq'][...] # plot fig, axes = gplt.subplots(ncols=2) axes[0].plot(epseq, sigeq) -gplt.patch(coor=coor+disp, conn=conn, cindex=Sigeq, cmap='jet', axis=axes[1], clim=(0.0, 0.1)) +gplt.patch(coor=coor + disp, conn=conn, cindex=Sigeq, cmap='jet', axis=axes[1], clim=(0.0, 0.1)) plt.show() - diff --git a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp index 17c87df..9588c44 100644 --- a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp +++ b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp @@ -1,239 +1,240 @@ #include +#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)); + // 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.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), + })); } - // 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.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; + // write ParaView meta-data + xdmf.write("main.xdmf"); + + return 0; } diff --git a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/plot.py b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/plot.py index 2ddfaa0..fefa86f 100644 --- a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/plot.py +++ b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/plot.py @@ -1,33 +1,32 @@ import h5py import matplotlib.pyplot as plt -import GooseMPL as gplt -import numpy as np +import GooseMPL as gplt +import numpy as np plt.style.use(['goose', 'goose-latex']) # open file file = h5py.File('main.h5', 'r') # read stored increments incs = file['/stored'][...] # read fields -coor = file['/coor' ][...] -conn = file['/conn' ][...] -disp = file['/disp' ][str(np.max(incs))][...][:,:2] +coor = file['/coor'][...] +conn = file['/conn'][...] +disp = file['/disp'][str(np.max(incs))][...][:, :2] Sigeq = file['/sigeq'][str(np.max(incs))][...] sigeq = file['/macroscopic/sigeq'][...] epseq = file['/macroscopic/epseq'][...] # plot fig, axes = gplt.subplots(ncols=2) axes[0].plot(epseq, sigeq) -gplt.patch(coor=coor+disp, conn=conn, cindex=Sigeq, cmap='jet', axis=axes[1], clim=(0.0, 0.1)) +gplt.patch(coor=coor + disp, conn=conn, cindex=Sigeq, cmap='jet', axis=axes[1], clim=(0.0, 0.1)) plt.show() - diff --git a/docs/examples/statics/Periodic_LinearElastic/main.cpp b/docs/examples/statics/Periodic_LinearElastic/main.cpp index cb51d1d..d0e9891 100644 --- a/docs/examples/statics/Periodic_LinearElastic/main.cpp +++ b/docs/examples/statics/Periodic_LinearElastic/main.cpp @@ -1,171 +1,168 @@ #include +#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.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; + // 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.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_LinearElastic/plot.py b/docs/examples/statics/Periodic_LinearElastic/plot.py index 6864116..18bccbd 100644 --- a/docs/examples/statics/Periodic_LinearElastic/plot.py +++ b/docs/examples/statics/Periodic_LinearElastic/plot.py @@ -1,26 +1,25 @@ import h5py import matplotlib.pyplot as plt -import GooseMPL as gplt -import numpy as np +import GooseMPL as gplt +import numpy as np plt.style.use(['goose', 'goose-latex']) # open file file = h5py.File('main.h5', 'r') # read fields -coor = file['/coor' ][...] -conn = file['/conn' ][...] -disp = file['/disp' ][...][:,:2] +coor = file['/coor'][...] +conn = file['/conn'][...] +disp = file['/disp'][...][:, :2] Sigeq = file['/sigeq'][...] # plot fig, ax = plt.subplots() -gplt.patch(coor=coor+disp, conn=conn, cindex=Sigeq, cmap='jet', axis=ax, clim=(0, 0.5)) +gplt.patch(coor=coor + disp, conn=conn, cindex=Sigeq, cmap='jet', axis=ax, clim=(0, 0.5)) plt.show() - diff --git a/docs/examples/statics/Periodic_NonLinearElastic/main.cpp b/docs/examples/statics/Periodic_NonLinearElastic/main.cpp index 10d420f..f529ef3 100644 --- a/docs/examples/statics/Periodic_NonLinearElastic/main.cpp +++ b/docs/examples/statics/Periodic_NonLinearElastic/main.cpp @@ -1,215 +1,214 @@ #include +#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; + // 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 - 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.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), + // ----- + + // 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; } - )); - - // write ParaView meta-data - xdmf.write("main.xdmf"); - return 0; + // 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.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/plot.py b/docs/examples/statics/Periodic_NonLinearElastic/plot.py index 0242896..f93b8ce 100644 --- a/docs/examples/statics/Periodic_NonLinearElastic/plot.py +++ b/docs/examples/statics/Periodic_NonLinearElastic/plot.py @@ -1,26 +1,25 @@ import h5py import matplotlib.pyplot as plt -import GooseMPL as gplt -import numpy as np +import GooseMPL as gplt +import numpy as np plt.style.use(['goose', 'goose-latex']) # open file file = h5py.File('main.h5', 'r') # read fields -coor = file['/coor' ][...] -conn = file['/conn' ][...] -disp = file['/disp' ][...][:,:2] +coor = file['/coor'][...] +conn = file['/conn'][...] +disp = file['/disp'][...][:, :2] Sigeq = file['/sigeq'][...] # plot fig, ax = plt.subplots() -gplt.patch(coor=coor+disp, conn=conn, cindex=Sigeq, cmap='jet', axis=ax, clim=(0, 0.1)) +gplt.patch(coor=coor + disp, conn=conn, cindex=Sigeq, cmap='jet', axis=ax, clim=(0, 0.1)) plt.show() -