diff --git a/docs/examples/statics/FixedDisplacements_LinearElastic/example.py b/docs/examples/statics/FixedDisplacements_LinearElastic/example.py index 160959a..1fcb97a 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/example.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/example.py @@ -1,181 +1,182 @@ 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 = K.Solve(fres, disp) +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/FixedDisplacements_LinearElastic/manual_partition.py b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py index b8c336f..46a47b5 100644 --- a/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py +++ b/docs/examples/statics/FixedDisplacements_LinearElastic/manual_partition.py @@ -1,197 +1,198 @@ 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 = K.Solve_u(fres_u, u_p) +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 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/example.py b/docs/examples/statics/MixedPeriodic_LinearElastic/example.py index d56432a..d4eaf42 100644 --- a/docs/examples/statics/MixedPeriodic_LinearElastic/example.py +++ b/docs/examples/statics/MixedPeriodic_LinearElastic/example.py @@ -1,188 +1,189 @@ 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 = K.Solve(fres, disp) +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/python/Matrix.hpp b/python/Matrix.hpp index 6d395e6..ad5256d 100644 --- a/python/Matrix.hpp +++ b/python/Matrix.hpp @@ -1,87 +1,93 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include #include "../include/GooseFEM/GooseFEM.h" -// ================================================================================================= - namespace py = pybind11; -// ================================================================================================= - void init_Matrix(py::module& m) { -py::class_>(m, "Matrix") +py::class_(m, "Matrix") .def(py::init< const xt::xtensor&, const xt::xtensor&>(), "Sparse matrix", py::arg("conn"), py::arg("dofs")) .def("nelem", - &GooseFEM::Matrix<>::nelem, + &GooseFEM::Matrix::nelem, "Number of element") .def("nne", - &GooseFEM::Matrix<>::nne, + &GooseFEM::Matrix::nne, "Number of nodes per element") .def("nnode", - &GooseFEM::Matrix<>::nnode, + &GooseFEM::Matrix::nnode, "Number of nodes") .def("ndim", - &GooseFEM::Matrix<>::ndim, + &GooseFEM::Matrix::ndim, "Number of dimensions") .def("ndof", - &GooseFEM::Matrix<>::ndof, + &GooseFEM::Matrix::ndof, "Number of degrees-of-freedom") .def("assemble", - &GooseFEM::Matrix<>::assemble, + &GooseFEM::Matrix::assemble, "Assemble matrix from 'elemmat", py::arg("elemmat")) .def("dofs", - &GooseFEM::Matrix<>::dofs, + &GooseFEM::Matrix::dofs, "Return degrees-of-freedom") .def("Dot", - py::overload_cast&>(&GooseFEM::Matrix<>::Dot, py::const_), + py::overload_cast&>(&GooseFEM::Matrix::Dot, py::const_), "Dot", py::arg("x")) .def("Dot", - py::overload_cast&>(&GooseFEM::Matrix<>::Dot, py::const_), + py::overload_cast&>(&GooseFEM::Matrix::Dot, py::const_), "Dot", py::arg("x")) + .def("__repr__", + [](const GooseFEM::Matrix&){ + return ""; }); + +py::class_>(m, "MatrixSolver") + + .def(py::init<>(), + "Sparse matrix solver") + .def("Solve", - py::overload_cast&>(&GooseFEM::Matrix<>::Solve), + py::overload_cast&>( + &GooseFEM::MatrixSolver<>::Solve), "Solve", + py::arg("matrix"), py::arg("b")) .def("Solve", - py::overload_cast&>(&GooseFEM::Matrix<>::Solve), + py::overload_cast&>( + &GooseFEM::MatrixSolver<>::Solve), "Solve", + py::arg("matrix"), py::arg("b")) .def("__repr__", - [](const GooseFEM::Matrix<>&){ - return ""; }); + [](const GooseFEM::MatrixSolver<>&){ + return ""; }); } - -// ================================================================================================= - diff --git a/python/MatrixPartitioned.hpp b/python/MatrixPartitioned.hpp index c9cfe94..9259e5c 100644 --- a/python/MatrixPartitioned.hpp +++ b/python/MatrixPartitioned.hpp @@ -1,127 +1,138 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include #include "../include/GooseFEM/GooseFEM.h" -// ================================================================================================= - namespace py = pybind11; -// ================================================================================================= - void init_MatrixPartitioned(py::module& m) { -py::class_>(m, "MatrixPartitioned") +py::class_(m, "MatrixPartitioned") .def(py::init< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>(), "Sparse, partitioned, matrix", py::arg("conn"), py::arg("dofs"), py::arg("iip")) .def("nelem", - &GooseFEM::MatrixPartitioned<>::nelem, + &GooseFEM::MatrixPartitioned::nelem, "Number of element") .def("nne", - &GooseFEM::MatrixPartitioned<>::nne, + &GooseFEM::MatrixPartitioned::nne, "Number of nodes per element") .def("nnode", - &GooseFEM::MatrixPartitioned<>::nnode, + &GooseFEM::MatrixPartitioned::nnode, "Number of nodes") .def("ndim", - &GooseFEM::MatrixPartitioned<>::ndim, + &GooseFEM::MatrixPartitioned::ndim, "Number of dimensions") .def("ndof", - &GooseFEM::MatrixPartitioned<>::ndof, + &GooseFEM::MatrixPartitioned::ndof, "Number of degrees-of-freedom") .def("nnu", - &GooseFEM::MatrixPartitioned<>::nnu, + &GooseFEM::MatrixPartitioned::nnu, "Number of unknown degrees-of-freedom") .def("nnp", - &GooseFEM::MatrixPartitioned<>::nnp, + &GooseFEM::MatrixPartitioned::nnp, "Number of prescribed degrees-of-freedom") .def("assemble", - &GooseFEM::MatrixPartitioned<>::assemble, + &GooseFEM::MatrixPartitioned::assemble, "Assemble matrix from 'elemmat", py::arg("elemmat")) .def("dofs", - &GooseFEM::MatrixPartitioned<>::dofs, + &GooseFEM::MatrixPartitioned::dofs, "Return degrees-of-freedom") .def("iiu", - &GooseFEM::MatrixPartitioned<>::iiu, + &GooseFEM::MatrixPartitioned::iiu, "Return unknown degrees-of-freedom") .def("iip", - &GooseFEM::MatrixPartitioned<>::iip, + &GooseFEM::MatrixPartitioned::iip, "Return prescribed degrees-of-freedom") - .def("Solve", - py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitioned<>::Solve), - "Solve", - py::arg("b"), - py::arg("x")) - - .def("Solve", - py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitioned<>::Solve), - "Solve", - py::arg("b"), - py::arg("x")) - - .def("Solve_u", - py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitioned<>::Solve_u), - "Solve_u", - py::arg("b_u"), - py::arg("x_p")) - .def("Reaction", py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitioned<>::Reaction, py::const_), + &GooseFEM::MatrixPartitioned::Reaction, py::const_), "Reaction", py::arg("x"), py::arg("b")) .def("Reaction", py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitioned<>::Reaction, py::const_), + &GooseFEM::MatrixPartitioned::Reaction, py::const_), "Reaction", py::arg("x"), py::arg("b")) .def("Reaction_p", py::overload_cast&, const xt::xtensor&>( - &GooseFEM::MatrixPartitioned<>::Reaction_p, py::const_), + &GooseFEM::MatrixPartitioned::Reaction_p, py::const_), "Reaction_p", py::arg("x_u"), py::arg("x_p")) .def("__repr__", - [](const GooseFEM::MatrixPartitioned<>&){ + [](const GooseFEM::MatrixPartitioned&){ return ""; }); -} +py::class_>(m, "MatrixPartitionedSolver") + + .def(py::init<>(), + "Sparse, partitioned, matrix solver") -// ================================================================================================= + .def("Solve", + py::overload_cast&, + const xt::xtensor&>( + &GooseFEM::MatrixPartitionedSolver<>::Solve), + "Solve", + py::arg("matrix"), + py::arg("b"), + py::arg("x")) + .def("Solve", + py::overload_cast&, + const xt::xtensor&>( + &GooseFEM::MatrixPartitionedSolver<>::Solve), + "Solve", + py::arg("matrix"), + py::arg("b"), + py::arg("x")) + + .def("Solve_u", + py::overload_cast&, + const xt::xtensor&>( + &GooseFEM::MatrixPartitionedSolver<>::Solve_u), + "Solve_u", + py::arg("matrix"), + py::arg("b_u"), + py::arg("x_p")) + + .def("__repr__", + [](const GooseFEM::MatrixPartitionedSolver<>&){ + return ""; }); + +}