diff --git a/docs/Matrix.rst b/docs/Matrix.rst new file mode 100644 index 0000000..df108cf --- /dev/null +++ b/docs/Matrix.rst @@ -0,0 +1,78 @@ + +************* +Sparse Matrix +************* + +Linear solver +============= + +The classes ``GooseFEM:::MatrixPartitioned`` and ``GooseFEM:::MatrixPartitionedTyings`` make use of a library to solver the linear system (stored as a sparse matrix). In particular the Eigen library and its plug-ins are used. To use the library's default solver: + +.. code-block:: cpp + + #include + #include + + int main() + { + ... + + GooseFEM::MatrixPartitioned<> K(...); + + ... + + return 0; + } + +The default solver can, however, be quite slow. Therefore Eigen has quite some `plug-ins `_ for the solver. GooseFEM allows the use of Eigen's Sparse Solver Concept to use such plug-ins. For example, to use SuiteSparse: + +.. code-block:: cpp + + #include + #include + #include + + int main() + { + ... + + GooseFEM::MatrixPartitioned>> K(...); + + ... + + return 0; + } + +SuiteSparse +----------- + +Preparation +^^^^^^^^^^^ + +1. `Download SuiteSparse `_. + +2. Extract the downloaded ``SuiteSparse-X.X.X.tar.gz```. + +3. Compile the library by: + + .. code-block:: bash + + cd /path/to/SuiteSparse + make + +Compiling GooseFEM code manually +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: bash + + clang++ -I/path/to/include/eigen3 -I/path/to/lapack/include -L/path/to/lapack/lib -I/path/to/openblas/include -L/path/to/openblas/lib -lopenblas -I/path/to/SuiteSparse/include -L/path/to/SuiteSparse/lib -lumfpack -lamd -lcholmod -lsuitesparseconfig -lm -std=c++14 -Wall -Wextra -pedantic -march=native -O3 -o example example.cpp + +.. tip:: + + * Replace ``-I/path/to/include/eigen3`` with + + .. code-block:: bash + + `pkg-config --cflags eigen3` + + diff --git a/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp b/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp index e06f24b..190bca7 100644 --- a/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp +++ b/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp @@ -1,170 +1,164 @@ #include #include -#include +#include // ------------------------------------------------------------------------------------------------- +namespace GM = GMatElastoPlasticQPot::Cartesian2d; namespace GF = GooseFEM; namespace QD = GooseFEM::Element::Quad4; -namespace GM = GMatElastoPlasticQPot::Cartesian2d; +namespace H5 = H5Easy; // ------------------------------------------------------------------------------------------------- inline double sqdot(const xt::xtensor &M, const xt::xtensor &V) { double out = 0.; - for ( size_t i = 0 ; i < M.size() ; ++i ) + for (size_t i = 0; i < M.size(); ++i) out += M(i) * V(i) * V(i); return out; } // ------------------------------------------------------------------------------------------------- int main() { // 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 // get mesh & quadrature GF::Mesh::Quad4::Regular mesh(nx,nx,1.); xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofsPeriodic(); GF::Vector vector(conn, dofs); QD::Quadrature quad(vector.AsElement(coor)); - size_t nnode = mesh.nnode(); 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 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}); // material definition GM::Matrix material({nelem, nip}); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::xtensor Isoft = xt::ones ({nelem, nip}); - - for ( size_t e = 0 ; e < nx*nx/4 ; ++e ) - for ( size_t q = 0 ; q < nip ; ++q ) - Ihard(e,q) = 1; - + 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.check(); // mass matrix QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w()); xt::xtensor rho = 1.0 * xt::ones({nelem, nip}); GF::MatrixDiagonal M(conn, dofs); M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho)); xt::xtensor mass = M.AsDiagonal(); // update in macroscopic deformation gradient xt::xtensor dFbar = xt::zeros({2,2}); dFbar(0,1) = gamma; - for ( size_t n = 0 ; n < nnode ; ++n ) - for ( size_t j = 0 ; j < ndim ; ++j ) - for ( size_t k = 0 ; k < ndim ; ++k ) - u(n,j) += dFbar(j,k) * ( coor(n,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 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(); // loop over increments - for ( size_t inc = 0 ; inc < static_cast(Epot.size()) ; ++inc ) + for (size_t inc = 0; inc < static_cast(Epot.size()); ++inc) { // store history xt::noalias(v_n) = v; xt::noalias(a_n) = a; // new displacement xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt,2.) * a; // compute strain/strain, and corresponding internal force vector.asElement(u, ue); quad.symGradN_vector(ue, Eps); - material.Sig(Eps, Sig); + material.stress(Eps, Sig); quad.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // compute residual force & solve xt::noalias(fres) = fext - fint; M.solve(fres, a); // new velocity - xt::noalias(v) = v_n + .5 * dt * ( a_n + a ); + xt::noalias(v) = v_n + .5 * dt * (a_n + a); // store output variables - xt::xtensor E = material.energy(Eps); + 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 - - HighFive::File file("example.hdf5", HighFive::File::Overwrite); - - xt::dump(file, "/global/Epot",Epot ); - xt::dump(file, "/global/Ekin",Ekin ); - xt::dump(file, "/global/t" ,t ); - xt::dump(file, "/mesh/conn" ,conn ); - xt::dump(file, "/mesh/coor" ,coor ); - xt::dump(file, "/mesh/disp" ,u ); + 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; } diff --git a/docs/examples/Dynamics/laminateElastic/noDamping/main_velocityVerlet.cpp b/docs/examples/Dynamics/laminateElastic/noDamping/main_velocityVerlet.cpp index edff410..46c2df3 100644 --- a/docs/examples/Dynamics/laminateElastic/noDamping/main_velocityVerlet.cpp +++ b/docs/examples/Dynamics/laminateElastic/noDamping/main_velocityVerlet.cpp @@ -1,190 +1,184 @@ #include #include -#include +#include // ------------------------------------------------------------------------------------------------- +namespace GM = GMatElastoPlasticQPot::Cartesian2d; namespace GF = GooseFEM; namespace QD = GooseFEM::Element::Quad4; -namespace GM = GMatElastoPlasticQPot::Cartesian2d; +namespace H5 = H5Easy; // ------------------------------------------------------------------------------------------------- inline double sqdot(const xt::xtensor &M, const xt::xtensor &V) { double out = 0.; - for ( size_t i = 0 ; i < M.size() ; ++i ) + for (size_t i = 0; i < M.size(); ++i) out += M(i) * V(i) * V(i); return out; } // ------------------------------------------------------------------------------------------------- int main() { // 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 // get mesh & quadrature GF::Mesh::Quad4::Regular mesh(nx,nx,1.); xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofsPeriodic(); GF::Vector vector(conn, dofs); QD::Quadrature quad(vector.AsElement(coor)); - size_t nnode = mesh.nnode(); 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 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}); // material definition GM::Matrix material({nelem, nip}); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::xtensor Isoft = xt::ones ({nelem, nip}); - - for ( size_t e = 0 ; e < nx*nx/4 ; ++e ) - for ( size_t q = 0 ; q < nip ; ++q ) - Ihard(e,q) = 1; - + 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.check(); // mass matrix QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w()); xt::xtensor rho = 1.0 * xt::ones({nelem, nip}); GF::MatrixDiagonal M(conn, dofs); M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho)); xt::xtensor mass = M.AsDiagonal(); // update in macroscopic deformation gradient xt::xtensor dFbar = xt::zeros({2,2}); dFbar(0,1) = gamma; - for ( size_t n = 0 ; n < nnode ; ++n ) - for ( size_t j = 0 ; j < ndim ; ++j ) - for ( size_t k = 0 ; k < ndim ; ++k ) - u(n,j) += dFbar(j,k) * ( coor(n,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 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(); // loop over increments - for ( size_t inc = 0 ; inc < static_cast(Epot.size()) ; ++inc ) + for (size_t inc = 0; inc < static_cast(Epot.size()); ++inc) { // store history xt::noalias(v_n) = v; xt::noalias(a_n) = a; // new displacement xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt,2.) * a; // compute strain/strain, and corresponding internal force vector.asElement(u, ue); quad.symGradN_vector(ue, Eps); - material.Sig(Eps, Sig); + material.stress(Eps, Sig); quad.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // estimate new velocity xt::noalias(v) = v_n + dt * a_n; // compute residual force & solve xt::noalias(fres) = fext - fint; M.solve(fres, a); // re-estimate new velocity - xt::noalias(v) = v_n + .5 * dt * ( a_n + a ); + xt::noalias(v) = v_n + .5 * dt * (a_n + a); // compute residual force & solve xt::noalias(fres) = fext - fint; M.solve(fres, a); // new velocity - xt::noalias(v) = v_n + .5 * dt * ( a_n + a ); + xt::noalias(v) = v_n + .5 * dt * (a_n + a); // compute residual force & solve xt::noalias(fres) = fext - fint; M.solve(fres, a); // store output variables - xt::xtensor E = material.energy(Eps); + 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 - - HighFive::File file("example.hdf5", HighFive::File::Overwrite); - - xt::dump(file, "/global/Epot",Epot ); - xt::dump(file, "/global/Ekin",Ekin ); - xt::dump(file, "/global/t" ,t ); - xt::dump(file, "/mesh/conn" ,conn ); - xt::dump(file, "/mesh/coor" ,coor ); - xt::dump(file, "/mesh/disp" ,u ); + 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; } diff --git a/docs/examples/Mesh/MeshHex8-FineLayer-nodes.py b/docs/examples/Mesh/MeshHex8-FineLayer-nodes.py index 51b7f6e..0664db4 100644 --- a/docs/examples/Mesh/MeshHex8-FineLayer-nodes.py +++ b/docs/examples/Mesh/MeshHex8-FineLayer-nodes.py @@ -1,151 +1,130 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Hex8.FineLayer(9,17,27) # initialize all node sets nodesets = dict( nodesFront = np.zeros((mesh.nnode()),dtype='int'), nodesBack = np.zeros((mesh.nnode()),dtype='int'), nodesLeft = np.zeros((mesh.nnode()),dtype='int'), nodesRight = np.zeros((mesh.nnode()),dtype='int'), nodesBottom = np.zeros((mesh.nnode()),dtype='int'), nodesTop = np.zeros((mesh.nnode()),dtype='int'), nodesFrontFace = np.zeros((mesh.nnode()),dtype='int'), nodesBackFace = np.zeros((mesh.nnode()),dtype='int'), nodesLeftFace = np.zeros((mesh.nnode()),dtype='int'), nodesRightFace = np.zeros((mesh.nnode()),dtype='int'), nodesBottomFace = np.zeros((mesh.nnode()),dtype='int'), nodesTopFace = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopRightCorner = np.zeros((mesh.nnode()),dtype='int'), ) # define node-sets nodesets['nodesFront' ][mesh.nodesFront() ] = 1 nodesets['nodesBack' ][mesh.nodesBack() ] = 1 nodesets['nodesLeft' ][mesh.nodesLeft() ] = 1 nodesets['nodesRight' ][mesh.nodesRight() ] = 1 nodesets['nodesBottom' ][mesh.nodesBottom() ] = 1 nodesets['nodesTop' ][mesh.nodesTop() ] = 1 nodesets['nodesFrontFace' ][mesh.nodesFrontFace() ] = 1 nodesets['nodesBackFace' ][mesh.nodesBackFace() ] = 1 nodesets['nodesLeftFace' ][mesh.nodesLeftFace() ] = 1 nodesets['nodesRightFace' ][mesh.nodesRightFace() ] = 1 nodesets['nodesBottomFace' ][mesh.nodesBottomFace() ] = 1 nodesets['nodesTopFace' ][mesh.nodesTopFace() ] = 1 nodesets['nodesFrontBottomEdge' ][mesh.nodesFrontBottomEdge() ] = 1 nodesets['nodesFrontTopEdge' ][mesh.nodesFrontTopEdge() ] = 1 nodesets['nodesFrontLeftEdge' ][mesh.nodesFrontLeftEdge() ] = 1 nodesets['nodesFrontRightEdge' ][mesh.nodesFrontRightEdge() ] = 1 nodesets['nodesBackBottomEdge' ][mesh.nodesBackBottomEdge() ] = 1 nodesets['nodesBackTopEdge' ][mesh.nodesBackTopEdge() ] = 1 nodesets['nodesBackLeftEdge' ][mesh.nodesBackLeftEdge() ] = 1 nodesets['nodesBackRightEdge' ][mesh.nodesBackRightEdge() ] = 1 nodesets['nodesBottomLeftEdge' ][mesh.nodesBottomLeftEdge() ] = 1 nodesets['nodesBottomRightEdge' ][mesh.nodesBottomRightEdge() ] = 1 nodesets['nodesTopLeftEdge' ][mesh.nodesTopLeftEdge() ] = 1 nodesets['nodesTopRightEdge' ][mesh.nodesTopRightEdge() ] = 1 nodesets['nodesFrontBottomOpenEdge' ][mesh.nodesFrontBottomOpenEdge() ] = 1 nodesets['nodesFrontTopOpenEdge' ][mesh.nodesFrontTopOpenEdge() ] = 1 nodesets['nodesFrontLeftOpenEdge' ][mesh.nodesFrontLeftOpenEdge() ] = 1 nodesets['nodesFrontRightOpenEdge' ][mesh.nodesFrontRightOpenEdge() ] = 1 nodesets['nodesBackBottomOpenEdge' ][mesh.nodesBackBottomOpenEdge() ] = 1 nodesets['nodesBackTopOpenEdge' ][mesh.nodesBackTopOpenEdge() ] = 1 nodesets['nodesBackLeftOpenEdge' ][mesh.nodesBackLeftOpenEdge() ] = 1 nodesets['nodesBackRightOpenEdge' ][mesh.nodesBackRightOpenEdge() ] = 1 nodesets['nodesBottomLeftOpenEdge' ][mesh.nodesBottomLeftOpenEdge() ] = 1 nodesets['nodesBottomRightOpenEdge' ][mesh.nodesBottomRightOpenEdge() ] = 1 nodesets['nodesTopLeftOpenEdge' ][mesh.nodesTopLeftOpenEdge() ] = 1 nodesets['nodesTopRightOpenEdge' ][mesh.nodesTopRightOpenEdge() ] = 1 nodesets['nodesFrontBottomLeftCorner' ][mesh.nodesFrontBottomLeftCorner() ] = 1 nodesets['nodesFrontBottomRightCorner'][mesh.nodesFrontBottomRightCorner()] = 1 nodesets['nodesFrontTopLeftCorner' ][mesh.nodesFrontTopLeftCorner() ] = 1 nodesets['nodesFrontTopRightCorner' ][mesh.nodesFrontTopRightCorner() ] = 1 nodesets['nodesBackBottomLeftCorner' ][mesh.nodesBackBottomLeftCorner() ] = 1 nodesets['nodesBackBottomRightCorner' ][mesh.nodesBackBottomRightCorner() ] = 1 nodesets['nodesBackTopLeftCorner' ][mesh.nodesBackTopLeftCorner() ] = 1 nodesets['nodesBackTopRightCorner' ][mesh.nodesBackTopRightCorner() ] = 1 # add DOF-numbers after eliminating periodicity nodesets['dofsPeriodic'] = mesh.dofsPeriodic()[:,0] -# open data file -name = 'MeshHex8-FineLayer-nodes' -file = h5py.File(name+'.hdf5','w') - -# write nodal positions and a dummy connectivity -file['coor'] = mesh.coor() -file['conn'] = np.arange(mesh.nnode()) +# filename of the HDF5-file +fname = 'MeshHex8-FineLayer-nodes.hdf5' -# write node-sets -for key,value in nodesets.items(): - file[key] = value +# write HDF-file containing the data -# ======================================== write XDMF-file ========================================= +with h5py.File(fname, 'w') as data: -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), -) + data.file['coor'] = mesh.coor() + data.file['conn'] = mesh.conn() -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Nodes") + for key, value in nodesets.items(): + data[key] = value -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Polyvertex", NumberOfElements='{nnode:d}'.format(**dims), NodesPerElement="1") -data = etree.SubElement(conn, "DataItem", Dimensions='{nnode:d} 1'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) +# write XDMF-file with metadata -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XYZ") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) +xdmf = pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Hexahedron, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +) -# add attributes -for key in nodesets: - attr = etree.SubElement(grid, "Attribute", Name=key, AttributeType="Scalar", Center="Node") - data = etree.SubElement(attr, "DataItem", Dimensions='{nnode:d}'.format(**dims), Format="HDF") - data.text = "{name:s}.hdf5:/{key:s}".format(name=name,key=key) +for key, value in nodesets.items(): + xdmf.push_back(pv.Attribute(fname, "/"+key, key, pv.AttributeType.Node, value.shape)) -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +xdmf.write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshHex8-FineLayer.py b/docs/examples/Mesh/MeshHex8-FineLayer.py index c00f71b..80be274 100644 --- a/docs/examples/Mesh/MeshHex8-FineLayer.py +++ b/docs/examples/Mesh/MeshHex8-FineLayer.py @@ -1,56 +1,34 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Hex8.FineLayer(9,17,27) -# open data file -name = 'MeshHex8-FineLayer' -file = h5py.File(name+'.hdf5','w') +# filename of the HDF5-file +fname = 'MeshHex8-FineLayer.hdf5' -# element set +# define element set elementsMiddleLayer = np.zeros((mesh.nelem()),dtype='int') elementsMiddleLayer[mesh.elementsMiddleLayer()] = 1 -# write nodal coordinates and connectivity -file['coor'] = mesh.coor() -file['conn'] = mesh.conn() -file['elementsMiddleLayer'] = elementsMiddleLayer - -# ======================================== write XDMF-file ========================================= +# write HDF-file containing the data -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne(), -) +with h5py.File(fname, 'w') as data: -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Mesh") + data['coor'] = mesh.coor() + data['conn'] = mesh.conn() + data['elementsMiddleLayer'] = elementsMiddleLayer -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Hexahedron", NumberOfElements='{nelem:d}'.format(**dims)) -data = etree.SubElement(conn, "DataItem", Dimensions='{nelem:d} {nne:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) +# write XDMF-file with metadata -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XYZ") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) +xdmf = pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Hexahedron, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +) -# add attributes -attr = etree.SubElement(grid, "Attribute", Name="elementsMiddleLayer", AttributeType="Scalar", Center="Cell") -data = etree.SubElement(attr, "DataItem", Dimensions='{nelem:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/elementsMiddleLayer".format(name) +xdmf.push_back(pv.Attribute(fname, "/elementsMiddleLayer", "elementsMiddleLayer", pv.AttributeType.Cell, elementsMiddleLayer.shape)) -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +xdmf.write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshHex8-Regular-nodes.py b/docs/examples/Mesh/MeshHex8-Regular-nodes.py index 03cbdbd..81b78b2 100644 --- a/docs/examples/Mesh/MeshHex8-Regular-nodes.py +++ b/docs/examples/Mesh/MeshHex8-Regular-nodes.py @@ -1,151 +1,130 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Hex8.Regular(6,8,12) # initialize all node sets nodesets = dict( nodesFront = np.zeros((mesh.nnode()),dtype='int'), nodesBack = np.zeros((mesh.nnode()),dtype='int'), nodesLeft = np.zeros((mesh.nnode()),dtype='int'), nodesRight = np.zeros((mesh.nnode()),dtype='int'), nodesBottom = np.zeros((mesh.nnode()),dtype='int'), nodesTop = np.zeros((mesh.nnode()),dtype='int'), nodesFrontFace = np.zeros((mesh.nnode()),dtype='int'), nodesBackFace = np.zeros((mesh.nnode()),dtype='int'), nodesLeftFace = np.zeros((mesh.nnode()),dtype='int'), nodesRightFace = np.zeros((mesh.nnode()),dtype='int'), nodesBottomFace = np.zeros((mesh.nnode()),dtype='int'), nodesTopFace = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBackRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesFrontBottomRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesFrontTopRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackBottomRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBackTopRightCorner = np.zeros((mesh.nnode()),dtype='int'), ) # define node-sets nodesets['nodesFront' ][mesh.nodesFront() ] = 1 nodesets['nodesBack' ][mesh.nodesBack() ] = 1 nodesets['nodesLeft' ][mesh.nodesLeft() ] = 1 nodesets['nodesRight' ][mesh.nodesRight() ] = 1 nodesets['nodesBottom' ][mesh.nodesBottom() ] = 1 nodesets['nodesTop' ][mesh.nodesTop() ] = 1 nodesets['nodesFrontFace' ][mesh.nodesFrontFace() ] = 1 nodesets['nodesBackFace' ][mesh.nodesBackFace() ] = 1 nodesets['nodesLeftFace' ][mesh.nodesLeftFace() ] = 1 nodesets['nodesRightFace' ][mesh.nodesRightFace() ] = 1 nodesets['nodesBottomFace' ][mesh.nodesBottomFace() ] = 1 nodesets['nodesTopFace' ][mesh.nodesTopFace() ] = 1 nodesets['nodesFrontBottomEdge' ][mesh.nodesFrontBottomEdge() ] = 1 nodesets['nodesFrontTopEdge' ][mesh.nodesFrontTopEdge() ] = 1 nodesets['nodesFrontLeftEdge' ][mesh.nodesFrontLeftEdge() ] = 1 nodesets['nodesFrontRightEdge' ][mesh.nodesFrontRightEdge() ] = 1 nodesets['nodesBackBottomEdge' ][mesh.nodesBackBottomEdge() ] = 1 nodesets['nodesBackTopEdge' ][mesh.nodesBackTopEdge() ] = 1 nodesets['nodesBackLeftEdge' ][mesh.nodesBackLeftEdge() ] = 1 nodesets['nodesBackRightEdge' ][mesh.nodesBackRightEdge() ] = 1 nodesets['nodesBottomLeftEdge' ][mesh.nodesBottomLeftEdge() ] = 1 nodesets['nodesBottomRightEdge' ][mesh.nodesBottomRightEdge() ] = 1 nodesets['nodesTopLeftEdge' ][mesh.nodesTopLeftEdge() ] = 1 nodesets['nodesTopRightEdge' ][mesh.nodesTopRightEdge() ] = 1 nodesets['nodesFrontBottomOpenEdge' ][mesh.nodesFrontBottomOpenEdge() ] = 1 nodesets['nodesFrontTopOpenEdge' ][mesh.nodesFrontTopOpenEdge() ] = 1 nodesets['nodesFrontLeftOpenEdge' ][mesh.nodesFrontLeftOpenEdge() ] = 1 nodesets['nodesFrontRightOpenEdge' ][mesh.nodesFrontRightOpenEdge() ] = 1 nodesets['nodesBackBottomOpenEdge' ][mesh.nodesBackBottomOpenEdge() ] = 1 nodesets['nodesBackTopOpenEdge' ][mesh.nodesBackTopOpenEdge() ] = 1 nodesets['nodesBackLeftOpenEdge' ][mesh.nodesBackLeftOpenEdge() ] = 1 nodesets['nodesBackRightOpenEdge' ][mesh.nodesBackRightOpenEdge() ] = 1 nodesets['nodesBottomLeftOpenEdge' ][mesh.nodesBottomLeftOpenEdge() ] = 1 nodesets['nodesBottomRightOpenEdge' ][mesh.nodesBottomRightOpenEdge() ] = 1 nodesets['nodesTopLeftOpenEdge' ][mesh.nodesTopLeftOpenEdge() ] = 1 nodesets['nodesTopRightOpenEdge' ][mesh.nodesTopRightOpenEdge() ] = 1 nodesets['nodesFrontBottomLeftCorner' ][mesh.nodesFrontBottomLeftCorner() ] = 1 nodesets['nodesFrontBottomRightCorner'][mesh.nodesFrontBottomRightCorner()] = 1 nodesets['nodesFrontTopLeftCorner' ][mesh.nodesFrontTopLeftCorner() ] = 1 nodesets['nodesFrontTopRightCorner' ][mesh.nodesFrontTopRightCorner() ] = 1 nodesets['nodesBackBottomLeftCorner' ][mesh.nodesBackBottomLeftCorner() ] = 1 nodesets['nodesBackBottomRightCorner' ][mesh.nodesBackBottomRightCorner() ] = 1 nodesets['nodesBackTopLeftCorner' ][mesh.nodesBackTopLeftCorner() ] = 1 nodesets['nodesBackTopRightCorner' ][mesh.nodesBackTopRightCorner() ] = 1 # add DOF-numbers after eliminating periodicity nodesets['dofsPeriodic'] = mesh.dofsPeriodic()[:,0] -# open data file -name = 'MeshHex8-Regular-nodes' -file = h5py.File(name+'.hdf5','w') - -# write nodal positions and a dummy connectivity -file['coor'] = mesh.coor() -file['conn'] = np.arange(mesh.nnode()) +# filename of the HDF5-file +fname = 'MeshHex8-Regular-nodes.hdf5' -# write node-sets -for key,value in nodesets.items(): - file[key] = value +# write HDF-file containing the data -# ======================================== write XDMF-file ========================================= +with h5py.File(fname, 'w') as data: -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), -) + data.file['coor'] = mesh.coor() + data.file['conn'] = mesh.conn() -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Nodes") + for key, value in nodesets.items(): + data[key] = value -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Polyvertex", NumberOfElements='{nnode:d}'.format(**dims), NodesPerElement="1") -data = etree.SubElement(conn, "DataItem", Dimensions='{nnode:d} 1'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) +# write XDMF-file with metadata -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XYZ") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) +xdmf = pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Hexahedron, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +) -# add attributes -for key in nodesets: - attr = etree.SubElement(grid, "Attribute", Name=key, AttributeType="Scalar", Center="Node") - data = etree.SubElement(attr, "DataItem", Dimensions='{nnode:d}'.format(**dims), Format="HDF") - data.text = "{name:s}.hdf5:/{key:s}".format(name=name,key=key) +for key, value in nodesets.items(): + xdmf.push_back(pv.Attribute(fname, "/"+key, key, pv.AttributeType.Node, value.shape)) -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +xdmf.write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshHex8-Regular.py b/docs/examples/Mesh/MeshHex8-Regular.py index 032a523..9ecd659 100644 --- a/docs/examples/Mesh/MeshHex8-Regular.py +++ b/docs/examples/Mesh/MeshHex8-Regular.py @@ -1,46 +1,25 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Hex8.Regular(6,8,12) -# open data file -name = 'MeshHex8-Regular' -file = h5py.File(name+'.hdf5','w') - -# write nodal coordinates and connectivity -file['coor'] = mesh.coor() -file['conn'] = mesh.conn() - -# ======================================== write XDMF-file ========================================= - -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne(), -) - -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Mesh") - -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Hexahedron", NumberOfElements='{nelem:d}'.format(**dims)) -data = etree.SubElement(conn, "DataItem", Dimensions='{nelem:d} {nne:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) - -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XYZ") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) - -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +# filename of the HDF5-file +fname = 'MeshHex8-Regular.hdf5' + +# write HDF-file containing the data + +with h5py.File(fname, 'w') as data: + + data['coor'] = mesh.coor() + data['conn'] = mesh.conn() + +# write XDMF-file with metadata + +pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Hexahedron, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +).write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshQuad4-FineLayer-nodes.py b/docs/examples/Mesh/MeshQuad4-FineLayer-nodes.py index 034cbc9..0f18200 100644 --- a/docs/examples/Mesh/MeshQuad4-FineLayer-nodes.py +++ b/docs/examples/Mesh/MeshQuad4-FineLayer-nodes.py @@ -1,87 +1,66 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Quad4.FineLayer(9,17) # initialize all node sets nodesets = dict( nodesBottomEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopEdge = np.zeros((mesh.nnode()),dtype='int'), nodesLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBottomRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesTopLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesTopRightCorner = np.zeros((mesh.nnode()),dtype='int'), ) # define node-sets nodesets['nodesBottomEdge' ][mesh.nodesBottomEdge() ] = 1 nodesets['nodesTopEdge' ][mesh.nodesTopEdge() ] = 1 nodesets['nodesLeftEdge' ][mesh.nodesLeftEdge() ] = 1 nodesets['nodesRightEdge' ][mesh.nodesRightEdge() ] = 1 nodesets['nodesBottomOpenEdge' ][mesh.nodesBottomOpenEdge() ] = 1 nodesets['nodesTopOpenEdge' ][mesh.nodesTopOpenEdge() ] = 1 nodesets['nodesLeftOpenEdge' ][mesh.nodesLeftOpenEdge() ] = 1 nodesets['nodesRightOpenEdge' ][mesh.nodesRightOpenEdge() ] = 1 nodesets['nodesBottomLeftCorner' ][mesh.nodesBottomLeftCorner() ] = 1 nodesets['nodesBottomRightCorner'][mesh.nodesBottomRightCorner()] = 1 nodesets['nodesTopLeftCorner' ][mesh.nodesTopLeftCorner() ] = 1 nodesets['nodesTopRightCorner' ][mesh.nodesTopRightCorner() ] = 1 # add DOF-numbers after eliminating periodicity nodesets['dofsPeriodic'] = mesh.dofsPeriodic()[:,0] -# open data file -name = 'MeshQuad4-FineLayer-nodes' -file = h5py.File(name+'.hdf5','w') - -# write nodal positions and a dummy connectivity -file['coor'] = mesh.coor() -file['conn'] = np.arange(mesh.nnode()) +# filename of the HDF5-file +fname = 'MeshQuad4-FineLayer-nodes.hdf5' -# write node-sets -for key,value in nodesets.items(): - file[key] = value +# write HDF-file containing the data -# ======================================== write XDMF-file ========================================= +with h5py.File(fname, 'w') as data: -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), -) + data.file['coor'] = mesh.coor() + data.file['conn'] = mesh.conn() -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Nodes") + for key, value in nodesets.items(): + data[key] = value -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Polyvertex", NumberOfElements='{nnode:d}'.format(**dims), NodesPerElement="1") -data = etree.SubElement(conn, "DataItem", Dimensions='{nnode:d} 1'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) +# write XDMF-file with metadata -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XY") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) +xdmf = pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Quadrilateral, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +) -# add attributes -for key in nodesets: - attr = etree.SubElement(grid, "Attribute", Name=key, AttributeType="Scalar", Center="Node") - data = etree.SubElement(attr, "DataItem", Dimensions='{nnode:d}'.format(**dims), Format="HDF") - data.text = "{name:s}.hdf5:/{key:s}".format(name=name,key=key) +for key, value in nodesets.items(): + xdmf.push_back(pv.Attribute(fname, "/"+key, key, pv.AttributeType.Node, value.shape)) -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +xdmf.write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshQuad4-FineLayer.py b/docs/examples/Mesh/MeshQuad4-FineLayer.py index 3932dd9..a8a3dcc 100644 --- a/docs/examples/Mesh/MeshQuad4-FineLayer.py +++ b/docs/examples/Mesh/MeshQuad4-FineLayer.py @@ -1,56 +1,34 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Quad4.FineLayer(9,17) -# open data file -name = 'MeshQuad4-FineLayer' -file = h5py.File(name+'.hdf5','w') +# filename of the HDF5-file +fname = 'MeshQuad4-FineLayer.hdf5' -# element set +# define element set elementsMiddleLayer = np.zeros((mesh.nelem()),dtype='int') elementsMiddleLayer[mesh.elementsMiddleLayer()] = 1 -# write nodal coordinates and connectivity -file['coor'] = mesh.coor() -file['conn'] = mesh.conn() -file['elementsMiddleLayer'] = elementsMiddleLayer - -# ======================================== write XDMF-file ========================================= +# write HDF-file containing the data -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne(), -) +with h5py.File(fname, 'w') as data: -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Mesh") + data['coor'] = mesh.coor() + data['conn'] = mesh.conn() + data['elementsMiddleLayer'] = elementsMiddleLayer -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Quadrilateral", NumberOfElements='{nelem:d}'.format(**dims)) -data = etree.SubElement(conn, "DataItem", Dimensions='{nelem:d} {nne:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) +# write XDMF-file with metadata -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XY") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) +xdmf = pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Quadrilateral, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +) -# add attributes -attr = etree.SubElement(grid, "Attribute", Name="elementsMiddleLayer", AttributeType="Scalar", Center="Cell") -data = etree.SubElement(attr, "DataItem", Dimensions='{nelem:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/elementsMiddleLayer".format(name) +xdmf.push_back(pv.Attribute(fname, "/elementsMiddleLayer", "elementsMiddleLayer", pv.AttributeType.Cell, elementsMiddleLayer.shape)) -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +xdmf.write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshQuad4-Regular-nodes.py b/docs/examples/Mesh/MeshQuad4-Regular-nodes.py index 1eda6d0..746bb34 100644 --- a/docs/examples/Mesh/MeshQuad4-Regular-nodes.py +++ b/docs/examples/Mesh/MeshQuad4-Regular-nodes.py @@ -1,87 +1,66 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Quad4.Regular(9,11) # initialize all node sets nodesets = dict( nodesBottomEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopEdge = np.zeros((mesh.nnode()),dtype='int'), nodesLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBottomRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesTopLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesTopRightCorner = np.zeros((mesh.nnode()),dtype='int'), ) # define node-sets nodesets['nodesBottomEdge' ][mesh.nodesBottomEdge() ] = 1 nodesets['nodesTopEdge' ][mesh.nodesTopEdge() ] = 1 nodesets['nodesLeftEdge' ][mesh.nodesLeftEdge() ] = 1 nodesets['nodesRightEdge' ][mesh.nodesRightEdge() ] = 1 nodesets['nodesBottomOpenEdge' ][mesh.nodesBottomOpenEdge() ] = 1 nodesets['nodesTopOpenEdge' ][mesh.nodesTopOpenEdge() ] = 1 nodesets['nodesLeftOpenEdge' ][mesh.nodesLeftOpenEdge() ] = 1 nodesets['nodesRightOpenEdge' ][mesh.nodesRightOpenEdge() ] = 1 nodesets['nodesBottomLeftCorner' ][mesh.nodesBottomLeftCorner() ] = 1 nodesets['nodesBottomRightCorner'][mesh.nodesBottomRightCorner()] = 1 nodesets['nodesTopLeftCorner' ][mesh.nodesTopLeftCorner() ] = 1 nodesets['nodesTopRightCorner' ][mesh.nodesTopRightCorner() ] = 1 # add DOF-numbers after eliminating periodicity nodesets['dofsPeriodic'] = mesh.dofsPeriodic()[:,0] -# open data file -name = 'MeshQuad4-Regular-nodes' -file = h5py.File(name+'.hdf5','w') - -# write nodal positions and a dummy connectivity -file['coor'] = mesh.coor() -file['conn'] = np.arange(mesh.nnode()) +# filename of the HDF5-file +fname = 'MeshQuad4-Regular-nodes.hdf5' -# write node-sets -for key,value in nodesets.items(): - file[key] = value +# write HDF-file containing the data -# ======================================== write XDMF-file ========================================= +with h5py.File(fname, 'w') as data: -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), -) + data.file['coor'] = mesh.coor() + data.file['conn'] = mesh.conn() -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Nodes") + for key, value in nodesets.items(): + data[key] = value -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Polyvertex", NumberOfElements='{nnode:d}'.format(**dims), NodesPerElement="1") -data = etree.SubElement(conn, "DataItem", Dimensions='{nnode:d} 1'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) +# write XDMF-file with metadata -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XY") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) +xdmf = pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Quadrilateral, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +) -# add attributes -for key in nodesets: - attr = etree.SubElement(grid, "Attribute", Name=key, AttributeType="Scalar", Center="Node") - data = etree.SubElement(attr, "DataItem", Dimensions='{nnode:d}'.format(**dims), Format="HDF") - data.text = "{name:s}.hdf5:/{key:s}".format(name=name,key=key) +for key, value in nodesets.items(): + xdmf.push_back(pv.Attribute(fname, "/"+key, key, pv.AttributeType.Node, value.shape)) -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +xdmf.write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshQuad4-Regular.py b/docs/examples/Mesh/MeshQuad4-Regular.py index 6c998f1..013ca73 100644 --- a/docs/examples/Mesh/MeshQuad4-Regular.py +++ b/docs/examples/Mesh/MeshQuad4-Regular.py @@ -1,46 +1,25 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Quad4.Regular(9,11) -# open data file -name = 'MeshQuad4-Regular' -file = h5py.File(name+'.hdf5','w') - -# write nodal coordinates and connectivity -file['coor'] = mesh.coor() -file['conn'] = mesh.conn() - -# ======================================== write XDMF-file ========================================= - -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne(), -) - -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Mesh") - -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Quadrilateral", NumberOfElements='{nelem:d}'.format(**dims)) -data = etree.SubElement(conn, "DataItem", Dimensions='{nelem:d} {nne:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) - -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XY") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) - -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +# filename of the HDF5-file +fname = 'MeshQuad4-Regular.hdf5' + +# write HDF-file containing the data + +with h5py.File(fname, 'w') as data: + + data['coor'] = mesh.coor() + data['conn'] = mesh.conn() + +# write XDMF-file with metadata + +pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Quadrilateral, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +).write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshTri3-Regular-nodes.py b/docs/examples/Mesh/MeshTri3-Regular-nodes.py index 1b23429..5d2e99e 100644 --- a/docs/examples/Mesh/MeshTri3-Regular-nodes.py +++ b/docs/examples/Mesh/MeshTri3-Regular-nodes.py @@ -1,87 +1,66 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Tri3.Regular(9,11) # initialize all node sets nodesets = dict( nodesBottomEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopEdge = np.zeros((mesh.nnode()),dtype='int'), nodesLeftEdge = np.zeros((mesh.nnode()),dtype='int'), nodesRightEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesTopOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesLeftOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesRightOpenEdge = np.zeros((mesh.nnode()),dtype='int'), nodesBottomLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesBottomRightCorner = np.zeros((mesh.nnode()),dtype='int'), nodesTopLeftCorner = np.zeros((mesh.nnode()),dtype='int'), nodesTopRightCorner = np.zeros((mesh.nnode()),dtype='int'), ) # define node-sets nodesets['nodesBottomEdge' ][mesh.nodesBottomEdge() ] = 1 nodesets['nodesTopEdge' ][mesh.nodesTopEdge() ] = 1 nodesets['nodesLeftEdge' ][mesh.nodesLeftEdge() ] = 1 nodesets['nodesRightEdge' ][mesh.nodesRightEdge() ] = 1 nodesets['nodesBottomOpenEdge' ][mesh.nodesBottomOpenEdge() ] = 1 nodesets['nodesTopOpenEdge' ][mesh.nodesTopOpenEdge() ] = 1 nodesets['nodesLeftOpenEdge' ][mesh.nodesLeftOpenEdge() ] = 1 nodesets['nodesRightOpenEdge' ][mesh.nodesRightOpenEdge() ] = 1 nodesets['nodesBottomLeftCorner' ][mesh.nodesBottomLeftCorner() ] = 1 nodesets['nodesBottomRightCorner'][mesh.nodesBottomRightCorner()] = 1 nodesets['nodesTopLeftCorner' ][mesh.nodesTopLeftCorner() ] = 1 nodesets['nodesTopRightCorner' ][mesh.nodesTopRightCorner() ] = 1 # add DOF-numbers after eliminating periodicity nodesets['dofsPeriodic'] = mesh.dofsPeriodic()[:,0] -# open data file -name = 'MeshTri3-Regular-nodes' -file = h5py.File(name+'.hdf5','w') - -# write nodal positions and a dummy connectivity -file['coor'] = mesh.coor() -file['conn'] = np.arange(mesh.nnode()) +# filename of the HDF5-file +fname = 'MeshTri3-Regular-nodes.hdf5' -# write node-sets -for key,value in nodesets.items(): - file[key] = value +# write HDF-file containing the data -# ======================================== write XDMF-file ========================================= +with h5py.File(fname, 'w') as data: -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), -) + data.file['coor'] = mesh.coor() + data.file['conn'] = mesh.conn() -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Nodes") + for key, value in nodesets.items(): + data[key] = value -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Polyvertex", NumberOfElements='{nnode:d}'.format(**dims), NodesPerElement="1") -data = etree.SubElement(conn, "DataItem", Dimensions='{nnode:d} 1'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) +# write XDMF-file with metadata -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XY") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) +xdmf = pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Triangle, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +) -# add attributes -for key in nodesets: - attr = etree.SubElement(grid, "Attribute", Name=key, AttributeType="Scalar", Center="Node") - data = etree.SubElement(attr, "DataItem", Dimensions='{nnode:d}'.format(**dims), Format="HDF") - data.text = "{name:s}.hdf5:/{key:s}".format(name=name,key=key) +for key, value in nodesets.items(): + xdmf.push_back(pv.Attribute(fname, "/"+key, key, pv.AttributeType.Node, value.shape)) -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +xdmf.write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Mesh/MeshTri3-Regular.py b/docs/examples/Mesh/MeshTri3-Regular.py index 40e0ed1..bf95ea8 100644 --- a/docs/examples/Mesh/MeshTri3-Regular.py +++ b/docs/examples/Mesh/MeshTri3-Regular.py @@ -1,46 +1,25 @@ -import h5py -import numpy as np -import GooseFEM as gf -import lxml.etree as etree - -# ====================== create fictitious configuration + store to HDF5-file ====================== +import h5py, os +import numpy as np +import GooseFEM as gf +import GooseFEM.ParaView.HDF5 as pv # create mesh object mesh = gf.Mesh.Tri3.Regular(9,11) -# open data file -name = 'MeshTri3-Regular' -file = h5py.File(name+'.hdf5','w') - -# write nodal coordinates and connectivity -file['coor'] = mesh.coor() -file['conn'] = mesh.conn() - -# ======================================== write XDMF-file ========================================= - -# get mesh dimensions -dims = dict( - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne(), -) - -# initialize file -root = etree.fromstring('') -domain = etree.SubElement(root, "Domain") -grid = etree.SubElement(domain, "Grid", Name="Mesh") - -# add connectivity -conn = etree.SubElement(grid, "Topology", TopologyType="Triangle", NumberOfElements='{nelem:d}'.format(**dims)) -data = etree.SubElement(conn, "DataItem", Dimensions='{nelem:d} {nne:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/conn".format(name) - -# add coordinates -coor = etree.SubElement(grid, "Geometry", GeometryType="XY") -data = etree.SubElement(coor, "DataItem", Dimensions='{nnode:d} {ndim:d}'.format(**dims), Format="HDF") -data.text = "{0:s}.hdf5:/coor".format(name) - -# write to file -open(name+'.xdmf','wb').write(etree.tostring(root, pretty_print=True)) +# filename of the HDF5-file +fname = 'MeshTri3-Regular.hdf5' + +# write HDF-file containing the data + +with h5py.File(fname, 'w') as data: + + data['coor'] = mesh.coor() + data['conn'] = mesh.conn() + +# write XDMF-file with metadata + +pv.Mesh( + pv.Connectivity(fname, "/conn", pv.ElementType.Triangle, mesh.conn().shape), + pv.Coordinates(fname, "/coor", mesh.coor().shape), +).write(os.path.splitext(fname)[0] + '.xdmf') diff --git a/docs/examples/Statics/FixedDisplacements/LinearElasticity/main.cpp b/docs/examples/Statics/FixedDisplacements/LinearElasticity/main.cpp index c1a7df4..edefdc0 100644 --- a/docs/examples/Statics/FixedDisplacements/LinearElasticity/main.cpp +++ b/docs/examples/Statics/FixedDisplacements/LinearElasticity/main.cpp @@ -1,147 +1,145 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5,5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLeft = mesh.nodesLeftEdge(); xt::xtensor nodesRight = mesh.nodesRightEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBottom = mesh.nodesBottomEdge(); // fixed displacements DOFs // ------------------------ xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesRight ), 0), xt::view(dofs, xt::keep(nodesTop ), 1), xt::view(dofs, xt::keep(nodesLeft ), 0), xt::view(dofs, xt::keep(nodesBottom), 1) )); // simulation variables // -------------------- // vector definition GooseFEM::VectorPartitioned vector(conn, dofs, iip); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // element/material definition // --------------------------- GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); GMatElastic::Cartesian3d::Matrix mat(nelem, nip, 1., 1.); // solve // ----- // allocate tensors size_t d = 3; xt::xtensor Eps = xt::empty({nelem, nip, d, d }); xt::xtensor Sig = xt::empty({nelem, nip, d, d }); xt::xtensor C = xt::empty({nelem, nip, d, d, d, d}); // allocate system matrix - GooseFEM::MatrixPartitioned K(conn, dofs, iip); + GooseFEM::MatrixPartitioned<> K(conn, dofs, iip); // 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(nodesRight ), 0) = +0.1; xt::view(disp, xt::keep(nodesTop ), 1) = -0.1; xt::view(disp, xt::keep(nodesLeft ), 0) = 0.0; xt::view(disp, xt::keep(nodesBottom), 1) = 0.0; // residual xt::noalias(fres) = fext - fint; // solve K.solve(fres, disp); // post-process // ------------ // compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.copy_p(fint, fext); // residual xt::noalias(fres) = fext - fint; // print residual std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl; // average stress per node xt::xtensor dV = elem.DV(2); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output - HighFive::File file("main.h5", HighFive::File::Overwrite); - xt::dump(file, "/coor", coor); xt::dump(file, "/conn", conn); xt::dump(file, "/disp", disp); xt::dump(file, "/Sig" , SigAv); return 0; } diff --git a/docs/examples/Statics/FixedDisplacements/LinearElasticity/main_manualPartition.cpp b/docs/examples/Statics/FixedDisplacements/LinearElasticity/main_manualPartition.cpp index 45e995d..c6fbbbb 100644 --- a/docs/examples/Statics/FixedDisplacements/LinearElasticity/main_manualPartition.cpp +++ b/docs/examples/Statics/FixedDisplacements/LinearElasticity/main_manualPartition.cpp @@ -1,164 +1,162 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5,5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLeft = mesh.nodesLeftEdge(); xt::xtensor nodesRight = mesh.nodesRightEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBottom = mesh.nodesBottomEdge(); // fixed displacements DOFs // ------------------------ xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesRight ), 0), xt::view(dofs, xt::keep(nodesTop ), 1), xt::view(dofs, xt::keep(nodesLeft ), 0), xt::view(dofs, xt::keep(nodesBottom), 1) )); // simulation variables // -------------------- // vector definition GooseFEM::VectorPartitioned vector(conn, dofs, iip); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // DOF values xt::xtensor u_u = xt::zeros({vector.nnu()}); xt::xtensor fres_u = xt::zeros({vector.nnu()}); xt::xtensor fext_p = xt::zeros({vector.nnp()}); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // element/material definition // --------------------------- GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); GMatElastic::Cartesian3d::Matrix mat(nelem, nip, 1., 1.); // solve // ----- // allocate tensors size_t d = 3; xt::xtensor Eps = xt::empty({nelem, nip, d, d }); xt::xtensor Sig = xt::empty({nelem, nip, d, d }); xt::xtensor C = xt::empty({nelem, nip, d, d, d, d}); // allocate system matrix - GooseFEM::MatrixPartitioned K(conn, dofs, iip); + GooseFEM::MatrixPartitioned<> K(conn, dofs, iip); // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent - mat.Tangent(Eps, Sig, C); + mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // set fixed displacements xt::xtensor u_p = xt::concatenate(xt::xtuple( +0.1 * xt::ones({nodesRight .size()}), -0.1 * xt::ones({nodesTop .size()}), 0.0 * xt::ones({nodesLeft .size()}), 0.0 * xt::ones({nodesBottom.size()}) )); // residual xt::noalias(fres) = fext - fint; // partition vector.asDofs_u(fres, fres_u); // solve K.solve_u(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.Sig(Eps, Sig); + mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.asDofs_p(fint, fext_p); // residual xt::noalias(fres) = fext - fint; // partition vector.asDofs_u(fres, fres_u); // print residual std::cout << xt::sum(xt::abs(fres_u))[0] / xt::sum(xt::abs(fext_p))[0] << std::endl; // average stress per node xt::xtensor dV = elem.DV(2); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output - HighFive::File file("main.h5", HighFive::File::Overwrite); - xt::dump(file, "/coor", coor); xt::dump(file, "/conn", conn); xt::dump(file, "/disp", disp); xt::dump(file, "/Sig" , SigAv); return 0; } diff --git a/docs/examples/Statics/Periodic/ElastoPlasticFiniteStrainSimo/main.cpp b/docs/examples/Statics/Periodic/ElastoPlasticFiniteStrainSimo/main.cpp index e723697..91112b3 100644 --- a/docs/examples/Statics/Periodic/ElastoPlasticFiniteStrainSimo/main.cpp +++ b/docs/examples/Statics/Periodic/ElastoPlasticFiniteStrainSimo/main.cpp @@ -1,238 +1,238 @@ #include #include #include #include #include namespace GM = GMatElastoPlasticFiniteStrainSimo::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor du = xt::zeros(coor.shape()); // iterative displacement update xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // DOF values xt::xtensor Fext = xt::zeros({tyings.nni()}); xt::xtensor Fint = xt::zeros({tyings.nni()}); // element/material definition // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem0(vector.AsElement(coor)); GF::Element::Quad4::QuadraturePlanar elem (vector.AsElement(coor)); size_t nip = elem0.nip(); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setLinearHardening(Isoft, 1.0, 1.0, 0.05, 0.05); mat.setElastic(Ihard, 1.0, 1.0); // solve // ----- // allocate tensors xt::xtensor I2 = mat.I2(); xt::xtensor F = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); // allocate internal variables double res; // some shear strain history xt::xtensor dgamma = 0.001 * xt::ones({101}); dgamma(0) = 0.0; // initialise output // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // loop over increments for (size_t inc = 0; inc < dgamma.size(); ++inc) { // update history mat.increment(); // iterate for (size_t iter = 0; ; ++iter) { // deformation gradient tensor vector.asElement(disp, ue); elem0.gradN_vector_T(ue, F); F += I2; // stress & tangent mat.tangent(F, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // residual xt::noalias(fres) = fext - fint; // check for convergence (skip the zeroth iteration, as the residual still vanishes) if (iter > 0) { // - internal/external force as DOFs (account for periodicity) vector.asDofs_i(fext, Fext); vector.asDofs_i(fint, Fint); // - extract reaction force vector.copy_p(Fint, Fext); // - norm of the residual and the reaction force double nfres = xt::sum(xt::abs(Fext-Fint))[0]; double nfext = xt::sum(xt::abs(Fext))[0]; // - relative residual, for convergence check if (nfext) res = nfres / nfext; else res = nfres; // - print progress to screen std::cout << "inc = " << inc << ", iter = " << iter << ", res = " << res << std::endl; // - check for convergence if (res < 1.0e-5) break; // - safe-guard from infinite loop if (iter > 20) throw std::runtime_error("Maximal number of iterations exceeded"); } // initialise displacement update du.fill(0.0); // set fixed displacements if (iter == 0) du(control_nodes(0),1) = dgamma(inc); // solve K.solve(fres, du); // add displacement update disp += du; // update shape functions elem.update_x(vector.AsElement(coor+disp)); } // post-process // - compute strain and stress vector.asElement(disp, ue); elem0.gradN_vector_T(ue, F); F += I2; GM::strain(F, Eps); mat.stress(F, Sig); // - integration point volume xt::xtensor dV = elem.DV(2); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: increment numbers H5::dump(file, "/stored", inc, {inc}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar), {inc}); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar), {inc}); // - write to output-file: element quantities H5::dump(file, "/sigeq/" + std::to_string(inc), GM::Sigeq(Sigelem)); H5::dump(file, "/epseq/" + std::to_string(inc), GM::Epseq(Epselem)); H5::dump(file, "/disp/" + std::to_string(inc), PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp/" + std::to_string(inc), "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq/" + std::to_string(inc), "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq/" + std::to_string(inc), "Eq. strain" , PV::AttributeType::Cell), } )); } // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/docs/examples/Statics/Periodic/ElastoPlasticity/main.cpp b/docs/examples/Statics/Periodic/ElastoPlasticity/main.cpp index 981411f..a78b5a3 100644 --- a/docs/examples/Statics/Periodic/ElastoPlasticity/main.cpp +++ b/docs/examples/Statics/Periodic/ElastoPlasticity/main.cpp @@ -1,228 +1,228 @@ #include #include #include #include #include namespace GM = GMatElastoPlastic::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor du = xt::zeros(coor.shape()); // iterative displacement update xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // DOF values xt::xtensor Fext = xt::zeros({tyings.nni()}); xt::xtensor Fint = xt::zeros({tyings.nni()}); // element/material definition // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); xt::xtensor dV = elem.DV(2); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setLinearHardening(Isoft, 1.0, 1.0, 0.05, 0.05); mat.setElastic(Ihard, 1.0, 1.0); // solve // ----- // allocate tensors xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); // allocate internal variables double res; // some shear strain history xt::xtensor dgamma = 0.001 * xt::ones({101}); dgamma(0) = 0.0; // initialise output // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // loop over increments for (size_t inc = 0; inc < dgamma.size(); ++inc) { // update history mat.increment(); // iterate for (size_t iter = 0; ; ++iter) { // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // residual xt::noalias(fres) = fext - fint; // check for convergence (skip the zeroth iteration, as the residual still vanishes) if (iter > 0) { // - internal/external force as DOFs (account for periodicity) vector.asDofs_i(fext, Fext); vector.asDofs_i(fint, Fint); // - extract reaction force vector.copy_p(Fint, Fext); // - norm of the residual and the reaction force double nfres = xt::sum(xt::abs(Fext-Fint))[0]; double nfext = xt::sum(xt::abs(Fext))[0]; // - relative residual, for convergence check if (nfext) res = nfres / nfext; else res = nfres; // - print progress to screen std::cout << "inc = " << inc << ", iter = " << iter << ", res = " << res << std::endl; // - check for convergence if (res < 1.0e-5) break; // - safe-guard from infinite loop if (iter > 20) throw std::runtime_error("Maximal number of iterations exceeded"); } // initialise displacement update du.fill(0.0); // set fixed displacements if (iter == 0) du(control_nodes(0),1) = dgamma(inc); // solve K.solve(fres, du); // add displacement update disp += du; } // post-process // - compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: increment numbers H5::dump(file, "/stored", inc, {inc}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar), {inc}); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar), {inc}); // - write to output-file: element quantities H5::dump(file, "/sigeq/" + std::to_string(inc), GM::Sigeq(Sigelem)); H5::dump(file, "/epseq/" + std::to_string(inc), GM::Epseq(Epselem)); H5::dump(file, "/disp/" + std::to_string(inc), PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp/" + std::to_string(inc), "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq/" + std::to_string(inc), "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq/" + std::to_string(inc), "Eq. strain" , PV::AttributeType::Cell), } )); } // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/docs/examples/Statics/Periodic/LinearElasticity/main.cpp b/docs/examples/Statics/Periodic/LinearElasticity/main.cpp index efa2817..f9fbfe6 100644 --- a/docs/examples/Statics/Periodic/LinearElasticity/main.cpp +++ b/docs/examples/Statics/Periodic/LinearElasticity/main.cpp @@ -1,170 +1,170 @@ #include #include #include #include #include namespace GM = GMatElastic::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // element/material definition // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setElastic(Isoft, 10.0, 1.0); mat.setElastic(Ihard, 10.0, 10.0); // solve // ----- // allocate tensors xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // residual xt::noalias(fres) = fext - fint; // set fixed displacements disp(control_nodes(0),1) = 0.1; - + // solve K.solve(fres, disp); // post-process // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // - integration point volume xt::xtensor dV = elem.DV(2); // - compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar)); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar)); // - write to output-file: element quantities H5::dump(file, "/sigeq", GM::Sigeq(Sigelem)); H5::dump(file, "/epseq", GM::Epseq(Epselem)); H5::dump(file, "/disp" , PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp" , "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq", "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq", "Eq. strain" , PV::AttributeType::Cell), } )); // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/docs/examples/Statics/Periodic/NonLinearElasticity/main.cpp b/docs/examples/Statics/Periodic/NonLinearElasticity/main.cpp index 485c321..e5495a7 100644 --- a/docs/examples/Statics/Periodic/NonLinearElasticity/main.cpp +++ b/docs/examples/Statics/Periodic/NonLinearElasticity/main.cpp @@ -1,214 +1,214 @@ #include #include #include #include #include namespace GM = GMatNonLinearElastic::Cartesian3d; namespace GF = GooseFEM; namespace PV = GooseFEM::ParaView::HDF5; namespace H5 = H5Easy; int main() { // mesh // ---- // define mesh GF::Mesh::Quad4::Regular mesh(5*10, 5*10); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); xt::xtensor elmat = mesh.elementMatrix(); // periodicity and fixed displacements DOFs // ---------------------------------------- // add control nodes GF::Tyings::Control control(coor, dofs); coor = control.coor(); dofs = control.dofs(); xt::xtensor control_dofs = control.controlDofs(); xt::xtensor control_nodes = control.controlNodes(); // extract fixed DOFs: // - all control nodes: to prescribe the deformation gradient // - one node of the mesh: to remove rigid body modes xt::xtensor iip = xt::concatenate(xt::xtuple( xt::reshape_view(control_dofs, {ndim*ndim}), xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim}) )); // get DOF-tyings, reorganise system GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip); dofs = tyings.dofs(); // simulation variables // -------------------- // vector definition: // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi()); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); // nodal displacement xt::xtensor du = xt::zeros(coor.shape()); // iterative displacement update xt::xtensor fint = xt::zeros(coor.shape()); // internal force xt::xtensor fext = xt::zeros(coor.shape()); // external force xt::xtensor fres = xt::zeros(coor.shape()); // residual force // element vectors / matrix xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // DOF values xt::xtensor Fext = xt::zeros({tyings.nni()}); xt::xtensor Fint = xt::zeros({tyings.nni()}); // element/material definition // --------------------------- // FEM quadrature GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); // material model // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed GM::Matrix mat(nelem, nip); size_t tdim = mat.ndim(); // some artificial material definition xt::xtensor ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10))); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setNonLinearElastic(Isoft, 10.0, 0.1, 0.1, 2.0); mat.setNonLinearElastic(Ihard, 10.0, 1.0, 0.1, 2.0); // solve // ----- // allocate tensors xt::xtensor Eps = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor Sig = xt::empty({nelem, nip, tdim, tdim}); xt::xtensor C = xt::empty({nelem, nip, tdim, tdim, tdim, tdim}); // allocate system matrix - GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp()); + GF::MatrixPartitionedTyings<> K(conn, dofs, tyings.Cdu(), tyings.Cdp()); // allocate internal variables double res; // iterate for (size_t iter = 0; ; ++iter) { // strain vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); // stress & tangent mat.tangent(Eps, Sig, C); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // stiffness matrix elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke); K.assemble(Ke); // residual xt::noalias(fres) = fext - fint; // check for convergence (skip the zeroth iteration, as the residual still vanishes) if (iter > 0) { // - internal/external force as DOFs (account for periodicity) vector.asDofs_i(fext, Fext); vector.asDofs_i(fint, Fint); // - extract reaction force vector.copy_p(Fint, Fext); // - norm of the residual and the reaction force double nfres = xt::sum(xt::abs(Fext-Fint))[0]; double nfext = xt::sum(xt::abs(Fext))[0]; // - relative residual, for convergence check if (nfext) res = nfres / nfext; else res = nfres; // - print progress to screen std::cout << "iter = " << iter << ", res = " << res << std::endl; // - check for convergence if (res < 1.0e-5) break; // - safe-guard from infinite loop if (iter > 20) throw std::runtime_error("Maximal number of iterations exceeded"); } // initialise displacement update du.fill(0.0); // set fixed displacements if (iter == 0) du(control_nodes(0),1) = 0.1; // solve K.solve(fres, du); // add displacement update disp += du; } // post-process // - output-file containing data HighFive::File file("main.h5", HighFive::File::Overwrite); // - ParaView meta-data PV::TimeSeries xdmf; // - write mesh H5::dump(file, "/conn", conn); H5::dump(file, "/coor", coor); // - integration point volume xt::xtensor dV = elem.DV(2); // - compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // - element average stress xt::xtensor Sigelem = xt::average(Sig, dV, {1}); xt::xtensor Epselem = xt::average(Eps, dV, {1}); // - macroscopic strain and stress xt::xtensor_fixed> Sigbar = xt::average(Sig, dV, {0,1}); xt::xtensor_fixed> Epsbar = xt::average(Eps, dV, {0,1}); // - write to output-file: macroscopic response H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar)); H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar)); // - write to output-file: element quantities H5::dump(file, "/sigeq", GM::Sigeq(Sigelem)); H5::dump(file, "/epseq", GM::Epseq(Epselem)); H5::dump(file, "/disp" , PV::as3d(disp)); // - update ParaView meta-data xdmf.push_back(PV::Increment( PV::Connectivity(file, "/conn", mesh.getElementType()), PV::Coordinates (file, "/coor"), { PV::Attribute(file, "/disp" , "Displacement", PV::AttributeType::Node), PV::Attribute(file, "/sigeq", "Eq. stress" , PV::AttributeType::Cell), PV::Attribute(file, "/epseq", "Eq. strain" , PV::AttributeType::Cell), } )); // write ParaView meta-data xdmf.write("main.xdmf"); return 0; } diff --git a/docs/examples/Statics/PeriodicPartial/LinearElasticity/main.cpp b/docs/examples/Statics/PeriodicPartial/LinearElasticity/main.cpp index 86f1848..8f53f81 100644 --- a/docs/examples/Statics/PeriodicPartial/LinearElasticity/main.cpp +++ b/docs/examples/Statics/PeriodicPartial/LinearElasticity/main.cpp @@ -1,159 +1,159 @@ #include #include #include #include int main() { // mesh // ---- // define mesh GooseFEM::Mesh::Quad4::Regular mesh(5,5); // mesh dimensions size_t nelem = mesh.nelem(); size_t nne = mesh.nne(); size_t ndim = mesh.ndim(); // mesh definitions xt::xtensor coor = mesh.coor(); xt::xtensor conn = mesh.conn(); xt::xtensor dofs = mesh.dofs(); // node sets xt::xtensor nodesLeft = mesh.nodesLeftOpenEdge(); xt::xtensor nodesRight = mesh.nodesRightOpenEdge(); xt::xtensor nodesTop = mesh.nodesTopEdge(); xt::xtensor nodesBottom = mesh.nodesBottomEdge(); // periodicity and fixed displacements DOFs // ---------------------------------------- for ( size_t i = 0 ; i < nodesLeft.size() ; ++i ) for ( size_t j = 0 ; j < coor.shape()[1] ; ++j ) dofs(nodesRight(i),j) = dofs(nodesLeft(i),j); dofs = GooseFEM::Mesh::renumber(dofs); xt::xtensor iip = xt::concatenate(xt::xtuple( xt::view(dofs, xt::keep(nodesBottom), 0), xt::view(dofs, xt::keep(nodesBottom), 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); // nodal quantities xt::xtensor disp = xt::zeros(coor.shape()); xt::xtensor fint = xt::zeros(coor.shape()); xt::xtensor fext = xt::zeros(coor.shape()); xt::xtensor fres = xt::zeros(coor.shape()); // element vectors xt::xtensor ue = xt::empty({nelem, nne, ndim}); xt::xtensor fe = xt::empty({nelem, nne, ndim}); xt::xtensor Ke = xt::empty({nelem, nne*ndim, nne*ndim}); // element/material definition // --------------------------- GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor)); size_t nip = elem.nip(); GMatElastic::Cartesian3d::Matrix mat(nelem, nip); xt::xtensor Ihard = xt::zeros({nelem, nip}); xt::view(Ihard, xt::keep(0,1,5,6), xt::all()) = 1; xt::xtensor Isoft = xt::ones({nelem, nip}) - Ihard; mat.setElastic(Isoft,10.,1. ); mat.setElastic(Ihard,10.,10.); // solve // ----- // allocate tensors size_t d = 3; xt::xtensor Eps = xt::empty({nelem, nip, d, d }); xt::xtensor Sig = xt::empty({nelem, nip, d, d }); xt::xtensor C = xt::empty({nelem, nip, d, d, d, d}); // allocate system matrix - GooseFEM::MatrixPartitioned K(conn, dofs, iip); + GooseFEM::MatrixPartitioned<> K(conn, dofs, iip); // 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), xt::keep(0)) = +0.1; // residual xt::noalias(fres) = fext - fint; // solve K.solve(fres, disp); // post-process // ------------ // compute strain and stress vector.asElement(disp, ue); elem.symGradN_vector(ue, Eps); mat.stress(Eps, Sig); // internal force elem.int_gradN_dot_tensor2_dV(Sig, fe); vector.assembleNode(fe, fint); // apply reaction force vector.copy_p(fint, fext); // residual xt::noalias(fres) = fext - fint; // print residual std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl; // average stress per node xt::xtensor dV = elem.DV(2); xt::xtensor SigAv = xt::average(Sig, dV, {1}); // write output HighFive::File file("main.h5", HighFive::File::Overwrite); xt::dump(file, "/coor", coor); xt::dump(file, "/conn", conn); xt::dump(file, "/disp", disp); xt::dump(file, "/Sig" , SigAv); return 0; } diff --git a/include/GooseFEM/MatrixPartitioned.h b/include/GooseFEM/MatrixPartitioned.h index 381ba97..bd7b12b 100644 --- a/include/GooseFEM/MatrixPartitioned.h +++ b/include/GooseFEM/MatrixPartitioned.h @@ -1,172 +1,173 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MATRIXPARTITIONED_H #define GOOSEFEM_MATRIXPARTITIONED_H // ------------------------------------------------------------------------------------------------- #include "config.h" #include #include #include // ================================================================================================= namespace GooseFEM { // ------------------------------------------------------------------------------------------------- +template >> class MatrixPartitioned { public: // Constructors MatrixPartitioned() = default; MatrixPartitioned( const xt::xtensor &conn, const xt::xtensor &dofs, const xt::xtensor &iip); // Dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs size_t nnu() const; // number of unknown DOFs size_t nnp() const; // number of prescribed DOFs // DOF lists xt::xtensor dofs() const; // DOFs xt::xtensor iiu() const; // unknown DOFs xt::xtensor iip() const; // prescribed DOFs // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] void assemble(const xt::xtensor &elemmat); // Solve: // x_u = A_uu \ ( b_u - A_up * x_p ) void solve( const xt::xtensor &b, xt::xtensor &x); // modified with "x_u" void solve( const xt::xtensor &b, xt::xtensor &x); // modified with "x_u" void solve_u( const xt::xtensor &b_u, const xt::xtensor &x_p, xt::xtensor &x_u); // overwritten // Get right-hand-size for corresponding to the prescribed DOFs: // b_p = A_pu * x_u + A_pp * x_p = A_pp * x_p void reaction( const xt::xtensor &x, xt::xtensor &b) const; // modified with "b_p" void reaction( const xt::xtensor &x, xt::xtensor &b) const; // modified with "b_p" void reaction_p( const xt::xtensor &x_u, const xt::xtensor &x_p, xt::xtensor &b_p) const; // overwritten // Auto-allocation of the functions above xt::xtensor Solve( const xt::xtensor &b, const xt::xtensor &x); xt::xtensor Solve( const xt::xtensor &b, const xt::xtensor &x); xt::xtensor Solve_u( const xt::xtensor &b_u, const xt::xtensor &x_p); xt::xtensor Reaction( const xt::xtensor &x, const xt::xtensor &b) const; xt::xtensor Reaction( const xt::xtensor &x, const xt::xtensor &b) const; xt::xtensor Reaction_p( const xt::xtensor &x_u, const xt::xtensor &x_p) const; private: // The matrix Eigen::SparseMatrix m_Auu; Eigen::SparseMatrix m_Aup; Eigen::SparseMatrix m_Apu; Eigen::SparseMatrix m_App; // Matrix entries std::vector> m_Tuu; std::vector> m_Tup; std::vector> m_Tpu; std::vector> m_Tpp; // Solver (re-used to solve different RHS) - Eigen::SimplicialLDLT> m_solver; + Solver m_solver; // Signal changes to data compare to the last inverse bool m_factor=false; // Bookkeeping xt::xtensor m_conn; // connectivity [nelem, nne ] xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] xt::xtensor m_part; // DOF-numbers per node, renumbered [nnode, ndim] xt::xtensor m_iiu; // unknown DOFs [nnu] xt::xtensor m_iip; // prescribed DOFs [nnp] // Dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs size_t m_nnu; // number of unknown DOFs size_t m_nnp; // number of prescribed DOFs // Compute inverse (automatically evaluated by "solve") void factorize(); // Convert arrays (Eigen version of VectorPartitioned, which contains public functions) Eigen::VectorXd asDofs_u(const xt::xtensor &dofval ) const; Eigen::VectorXd asDofs_u(const xt::xtensor &nodevec) const; Eigen::VectorXd asDofs_p(const xt::xtensor &dofval ) const; Eigen::VectorXd asDofs_p(const xt::xtensor &nodevec) const; }; // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #include "MatrixPartitioned.hpp" // ================================================================================================= #endif diff --git a/include/GooseFEM/MatrixPartitioned.hpp b/include/GooseFEM/MatrixPartitioned.hpp index 770237c..c6f1ed7 100644 --- a/include/GooseFEM/MatrixPartitioned.hpp +++ b/include/GooseFEM/MatrixPartitioned.hpp @@ -1,416 +1,449 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MATRIXPARTITIONED_HPP #define GOOSEFEM_MATRIXPARTITIONED_HPP // ------------------------------------------------------------------------------------------------- #include "MatrixPartitioned.h" #include "Mesh.h" // ================================================================================================= namespace GooseFEM { // ------------------------------------------------------------------------------------------------- -inline MatrixPartitioned::MatrixPartitioned( +template +inline MatrixPartitioned::MatrixPartitioned( const xt::xtensor &conn, const xt::xtensor &dofs, const xt::xtensor &iip) : m_conn(conn), m_dofs(dofs), m_iip(iip) { m_nelem = m_conn.shape()[0]; m_nne = m_conn.shape()[1]; m_nnode = m_dofs.shape()[0]; m_ndim = m_dofs.shape()[1]; m_iiu = xt::setdiff1d(dofs, iip); m_ndof = xt::amax(m_dofs)[0] + 1; m_nnp = m_iip.size(); m_nnu = m_iiu.size(); m_part = Mesh::Reorder({m_iiu, m_iip}).get(m_dofs); m_Tuu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tup.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tpu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tpp.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Auu.resize(m_nnu,m_nnu); m_Aup.resize(m_nnu,m_nnp); m_Apu.resize(m_nnp,m_nnu); m_App.resize(m_nnp,m_nnp); GOOSEFEM_ASSERT(xt::amax(m_conn)[0] + 1 == m_nnode); GOOSEFEM_ASSERT(xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0]); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitioned::nelem() const +template +inline size_t MatrixPartitioned::nelem() const { return m_nelem; } -inline size_t MatrixPartitioned::nne() const +template +inline size_t MatrixPartitioned::nne() const { return m_nne; } -inline size_t MatrixPartitioned::nnode() const +template +inline size_t MatrixPartitioned::nnode() const { return m_nnode; } -inline size_t MatrixPartitioned::ndim() const +template +inline size_t MatrixPartitioned::ndim() const { return m_ndim; } -inline size_t MatrixPartitioned::ndof() const +template +inline size_t MatrixPartitioned::ndof() const { return m_ndof; } -inline size_t MatrixPartitioned::nnu() const +template +inline size_t MatrixPartitioned::nnu() const { return m_nnu; } -inline size_t MatrixPartitioned::nnp() const +template +inline size_t MatrixPartitioned::nnp() const { return m_nnp; } -inline xt::xtensor MatrixPartitioned::dofs() const +template +inline xt::xtensor MatrixPartitioned::dofs() const { return m_dofs; } -inline xt::xtensor MatrixPartitioned::iiu() const +template +inline xt::xtensor MatrixPartitioned::iiu() const { return m_iiu; } -inline xt::xtensor MatrixPartitioned::iip() const +template +inline xt::xtensor MatrixPartitioned::iip() const { return m_iip; } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::factorize() +template +inline void MatrixPartitioned::factorize() { if (!m_factor) return; m_solver.compute(m_Auu); m_factor = false; } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::assemble(const xt::xtensor &elemmat) +template +inline void MatrixPartitioned::assemble(const xt::xtensor &elemmat) { GOOSEFEM_ASSERT(elemmat.shape() ==\ std::decay_t::shape_type({m_nelem, m_nne*m_ndim, m_nne*m_ndim})); m_Tuu.clear(); m_Tup.clear(); m_Tpu.clear(); m_Tpp.clear(); for (size_t e = 0 ; e < m_nelem ; ++e) { for (size_t m = 0 ; m < m_nne ; ++m) { for (size_t i = 0 ; i < m_ndim ; ++i) { size_t di = m_part(m_conn(e,m),i); for (size_t n = 0 ; n < m_nne ; ++n) { for (size_t j = 0 ; j < m_ndim ; ++j) { size_t dj = m_part(m_conn(e,n),j); if (di < m_nnu and dj < m_nnu) m_Tuu.push_back(Eigen::Triplet(di ,dj ,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (di < m_nnu) m_Tup.push_back(Eigen::Triplet(di ,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (dj < m_nnu) m_Tpu.push_back(Eigen::Triplet(di-m_nnu,dj ,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else m_Tpp.push_back(Eigen::Triplet(di-m_nnu,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j))); } } } } } m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end()); m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end()); m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end()); m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end()); m_factor = true; } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::solve( +template +inline void MatrixPartitioned::solve( const xt::xtensor &b, xt::xtensor &x) { GOOSEFEM_ASSERT(b.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); GOOSEFEM_ASSERT(x.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); this->factorize(); Eigen::VectorXd B_u = this->asDofs_u(b); Eigen::VectorXd X_p = this->asDofs_p(x); Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_Aup * X_p)); #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) for (size_t i = 0 ; i < m_ndim ; ++i) if ( m_part(m,i) < m_nnu ) x(m,i) = X_u(m_part(m,i)); } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::solve( +template +inline void MatrixPartitioned::solve( const xt::xtensor &b, xt::xtensor &x) { GOOSEFEM_ASSERT(b.size() == m_ndof); GOOSEFEM_ASSERT(x.size() == m_ndof); this->factorize(); Eigen::VectorXd B_u = this->asDofs_u(b); Eigen::VectorXd X_p = this->asDofs_p(x); Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_Aup * X_p)); #pragma omp parallel for for (size_t d = 0 ; d < m_nnu ; ++d) x(m_iiu(d)) = X_u(d); } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::solve_u( +template +inline void MatrixPartitioned::solve_u( const xt::xtensor &b_u, const xt::xtensor &x_p, xt::xtensor &x_u) { GOOSEFEM_ASSERT(b_u.size() == m_nnu); GOOSEFEM_ASSERT(x_p.size() == m_nnp); GOOSEFEM_ASSERT(x_u.size() == m_nnu); this->factorize(); Eigen::VectorXd B_u(m_nnu,1); Eigen::VectorXd X_p(m_nnp,1); std::copy(b_u.begin(), b_u.end(), B_u.data()); std::copy(x_p.begin(), x_p.end(), X_p.data()); Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_Aup * X_p)); std::copy(X_u.data(), X_u.data()+m_nnu, x_u.begin()); } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::reaction( +template +inline void MatrixPartitioned::reaction( const xt::xtensor &x, xt::xtensor &b) const { GOOSEFEM_ASSERT(x.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); GOOSEFEM_ASSERT(b.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); Eigen::VectorXd X_u = this->asDofs_u(x); Eigen::VectorXd X_p = this->asDofs_p(x); Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p; #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) for (size_t i = 0 ; i < m_ndim ; ++i) if (m_part(m,i) >= m_nnu) b(m,i) = B_p(m_part(m,i)-m_nnu); } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::reaction( +template +inline void MatrixPartitioned::reaction( const xt::xtensor &x, xt::xtensor &b) const { GOOSEFEM_ASSERT(x.size() == m_ndof); GOOSEFEM_ASSERT(b.size() == m_ndof); Eigen::VectorXd X_u = this->asDofs_u(x); Eigen::VectorXd X_p = this->asDofs_p(x); Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p; #pragma omp parallel for for (size_t d = 0 ; d < m_nnp ; ++d) b(m_iip(d)) = B_p(d); } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitioned::reaction_p( +template +inline void MatrixPartitioned::reaction_p( const xt::xtensor &x_u, const xt::xtensor &x_p, xt::xtensor &b_p) const { GOOSEFEM_ASSERT(x_u.size() == m_nnu); GOOSEFEM_ASSERT(x_p.size() == m_nnp); GOOSEFEM_ASSERT(b_p.size() == m_nnp); Eigen::VectorXd X_u(m_nnu,1); Eigen::VectorXd X_p(m_nnp,1); std::copy(x_u.begin(), x_u.end(), X_u.data()); std::copy(x_p.begin(), x_p.end(), X_p.data()); Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p; std::copy(B_p.data(), B_p.data()+m_nnp, b_p.begin()); } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitioned::Solve( +template +inline xt::xtensor MatrixPartitioned::Solve( const xt::xtensor &b, const xt::xtensor &x) { xt::xtensor out = x; this->solve(b, out); return out; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitioned::Solve( +template +inline xt::xtensor MatrixPartitioned::Solve( const xt::xtensor &b, const xt::xtensor &x) { xt::xtensor out = x; this->solve(b, out); return out; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitioned::Solve_u( +template +inline xt::xtensor MatrixPartitioned::Solve_u( const xt::xtensor &b_u, const xt::xtensor &x_p) { xt::xtensor x_u = xt::empty({m_nnu}); this->solve_u(b_u, x_p, x_u); return x_u; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitioned::Reaction( +template +inline xt::xtensor MatrixPartitioned::Reaction( const xt::xtensor &x, const xt::xtensor &b) const { xt::xtensor out = b; this->reaction(x, out); return out; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitioned::Reaction( +template +inline xt::xtensor MatrixPartitioned::Reaction( const xt::xtensor &x, const xt::xtensor &b) const { xt::xtensor out = b; this->reaction(x, out); return out; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitioned::Reaction_p( +template +inline xt::xtensor MatrixPartitioned::Reaction_p( const xt::xtensor &x_u, const xt::xtensor &x_p) const { xt::xtensor b_p = xt::empty({m_nnp}); this->reaction_p(x_u, x_p, b_p); return b_p; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitioned::asDofs_u(const xt::xtensor &dofval) const +template +inline Eigen::VectorXd MatrixPartitioned::asDofs_u( + const xt::xtensor &dofval) const { assert(dofval.size() == m_ndof); Eigen::VectorXd dofval_u(m_nnu,1); #pragma omp parallel for for (size_t d = 0 ; d < m_nnu ; ++d) dofval_u(d) = dofval(m_iiu(d)); return dofval_u; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitioned::asDofs_u(const xt::xtensor &nodevec) const +template +inline Eigen::VectorXd MatrixPartitioned::asDofs_u( + const xt::xtensor &nodevec) const { assert(nodevec.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); Eigen::VectorXd dofval_u(m_nnu,1); #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) for (size_t i = 0 ; i < m_ndim ; ++i) if (m_part(m,i) < m_nnu) dofval_u(m_part(m,i)) = nodevec(m,i); return dofval_u; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitioned::asDofs_p(const xt::xtensor &dofval) const +template +inline Eigen::VectorXd MatrixPartitioned::asDofs_p( + const xt::xtensor &dofval) const { assert(dofval.size() == m_ndof); Eigen::VectorXd dofval_p(m_nnp,1); #pragma omp parallel for for (size_t d = 0 ; d < m_nnp ; ++d) dofval_p(d) = dofval(m_iip(d)); return dofval_p; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitioned::asDofs_p(const xt::xtensor &nodevec) const +template +inline Eigen::VectorXd MatrixPartitioned::asDofs_p( + const xt::xtensor &nodevec) const { assert(nodevec.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); Eigen::VectorXd dofval_p(m_nnp,1); #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) for (size_t i = 0 ; i < m_ndim ; ++i) if (m_part(m,i) >= m_nnu) dofval_p(m_part(m,i)-m_nnu) = nodevec(m,i); return dofval_p; } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/GooseFEM/MatrixPartitionedTyings.h b/include/GooseFEM/MatrixPartitionedTyings.h index 5e07396..930b3d0 100644 --- a/include/GooseFEM/MatrixPartitionedTyings.h +++ b/include/GooseFEM/MatrixPartitionedTyings.h @@ -1,166 +1,167 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_H #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_H // ------------------------------------------------------------------------------------------------- #include "config.h" #include #include #include // ================================================================================================= namespace GooseFEM { // ------------------------------------------------------------------------------------------------- +template >> class MatrixPartitionedTyings { public: // Constructors MatrixPartitionedTyings() = default; MatrixPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp); // Dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs size_t nnu() const; // number of independent, unknown DOFs size_t nnp() const; // number of independent, prescribed DOFs size_t nni() const; // number of independent DOFs size_t nnd() const; // number of dependent DOFs // DOF lists xt::xtensor dofs() const; // DOFs xt::xtensor iiu() const; // independent, unknown DOFs xt::xtensor iip() const; // independent, prescribed DOFs xt::xtensor iii() const; // independent DOFs xt::xtensor iid() const; // dependent DOFs // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] void assemble(const xt::xtensor &elemmat); // Solve: // A' = A_ii + K_id * C_di + C_di^T * K_di + C_di^T * K_dd * C_di // b' = b_i + C_di^T * b_d // x_u = A'_uu \ ( b'_u - A'_up * x_p ) // x_i = [x_u, x_p] // x_d = C_di * x_i void solve( const xt::xtensor &b, xt::xtensor &x); // modified with "x_u", "x_d" void solve( const xt::xtensor &b, xt::xtensor &x); // modified with "x_u", "x_d" void solve_u( const xt::xtensor &b_u, const xt::xtensor &b_d, const xt::xtensor &x_p, xt::xtensor &x_u); // overwritten private: // The matrix Eigen::SparseMatrix m_Auu; Eigen::SparseMatrix m_Aup; Eigen::SparseMatrix m_Apu; Eigen::SparseMatrix m_App; Eigen::SparseMatrix m_Aud; Eigen::SparseMatrix m_Apd; Eigen::SparseMatrix m_Adu; Eigen::SparseMatrix m_Adp; Eigen::SparseMatrix m_Add; // The matrix for which the tyings have been applied Eigen::SparseMatrix m_ACuu; Eigen::SparseMatrix m_ACup; Eigen::SparseMatrix m_ACpu; Eigen::SparseMatrix m_ACpp; // Matrix entries std::vector> m_Tuu; std::vector> m_Tup; std::vector> m_Tpu; std::vector> m_Tpp; std::vector> m_Tud; std::vector> m_Tpd; std::vector> m_Tdu; std::vector> m_Tdp; std::vector> m_Tdd; // Solver (re-used to solve different RHS) - Eigen::SimplicialLDLT> m_solver; + Solver m_solver; // Signal changes to data compare to the last inverse bool m_factor=false; // Bookkeeping xt::xtensor m_conn; // connectivity [nelem, nne ] xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] xt::xtensor m_iiu; // unknown DOFs [nnu] xt::xtensor m_iip; // prescribed DOFs [nnp] xt::xtensor m_iid; // dependent DOFs [nnd] // Dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs size_t m_nnu; // number of independent, unknown DOFs size_t m_nnp; // number of independent, prescribed DOFs size_t m_nni; // number of independent DOFs size_t m_nnd; // number of dependent DOFs // Tyings Eigen::SparseMatrix m_Cdu; Eigen::SparseMatrix m_Cdp; Eigen::SparseMatrix m_Cud; Eigen::SparseMatrix m_Cpd; // Compute inverse (automatically evaluated by "solve") void factorize(); // Convert arrays (Eigen version of VectorPartitioned, which contains public functions) Eigen::VectorXd asDofs_u(const xt::xtensor &dofval ) const; Eigen::VectorXd asDofs_u(const xt::xtensor &nodevec) const; Eigen::VectorXd asDofs_p(const xt::xtensor &dofval ) const; Eigen::VectorXd asDofs_p(const xt::xtensor &nodevec) const; Eigen::VectorXd asDofs_d(const xt::xtensor &dofval ) const; Eigen::VectorXd asDofs_d(const xt::xtensor &nodevec) const; }; // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #include "MatrixPartitionedTyings.hpp" // ================================================================================================= #endif diff --git a/include/GooseFEM/MatrixPartitionedTyings.hpp b/include/GooseFEM/MatrixPartitionedTyings.hpp index 6abf492..1590494 100644 --- a/include/GooseFEM/MatrixPartitionedTyings.hpp +++ b/include/GooseFEM/MatrixPartitionedTyings.hpp @@ -1,438 +1,462 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP // ------------------------------------------------------------------------------------------------- #include "MatrixPartitionedTyings.h" // ================================================================================================= namespace GooseFEM { // ------------------------------------------------------------------------------------------------- -inline MatrixPartitionedTyings::MatrixPartitionedTyings( +template +inline MatrixPartitionedTyings::MatrixPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp) : m_conn(conn), m_dofs(dofs), m_Cdu(Cdu), m_Cdp(Cdp) { GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows()); m_nnu = static_cast(m_Cdu.cols()); m_nnp = static_cast(m_Cdp.cols()); m_nnd = static_cast(m_Cdp.rows()); m_nni = m_nnu + m_nnp; m_ndof = m_nni + m_nnd; m_iiu = xt::arange(m_nnu); m_iip = xt::arange(m_nnu, m_nnu + m_nnp); m_iid = xt::arange(m_nni, m_nni + m_nnd); m_nelem = m_conn.shape()[0]; m_nne = m_conn.shape()[1]; m_nnode = m_dofs.shape()[0]; m_ndim = m_dofs.shape()[1]; m_Cud = m_Cdu.transpose(); m_Cpd = m_Cdp.transpose(); m_Tuu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tup.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tpu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tpp.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tud.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tpd.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tdu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tdp.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Tdd.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); m_Auu.resize(m_nnu,m_nnu); m_Aup.resize(m_nnu,m_nnp); m_Apu.resize(m_nnp,m_nnu); m_App.resize(m_nnp,m_nnp); m_Aud.resize(m_nnu,m_nnd); m_Apd.resize(m_nnp,m_nnd); m_Adu.resize(m_nnd,m_nnu); m_Adp.resize(m_nnd,m_nnp); m_Add.resize(m_nnd,m_nnd); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)[0] + 1); } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitionedTyings::nelem() const +template +inline size_t MatrixPartitionedTyings::nelem() const { return m_nelem; } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitionedTyings::nne() const +template +inline size_t MatrixPartitionedTyings::nne() const { return m_nne; } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitionedTyings::nnode() const +template +inline size_t MatrixPartitionedTyings::nnode() const { return m_nnode; } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitionedTyings::ndim() const +template +inline size_t MatrixPartitionedTyings::ndim() const { return m_ndim; } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitionedTyings::ndof() const +template +inline size_t MatrixPartitionedTyings::ndof() const { return m_ndof; } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitionedTyings::nnu() const +template +inline size_t MatrixPartitionedTyings::nnu() const { return m_nnu; } // ------------------------------------------------------------------------------------------------- -inline size_t MatrixPartitionedTyings::nnp() const +template +inline size_t MatrixPartitionedTyings::nnp() const { return m_nnp; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitionedTyings::dofs() const +template +inline xt::xtensor MatrixPartitionedTyings::dofs() const { return m_dofs; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitionedTyings::iiu() const +template +inline xt::xtensor MatrixPartitionedTyings::iiu() const { return m_iiu; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitionedTyings::iip() const +template +inline xt::xtensor MatrixPartitionedTyings::iip() const { return m_iip; } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitionedTyings::iii() const +template +inline xt::xtensor MatrixPartitionedTyings::iii() const { return xt::arange(m_nni); } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor MatrixPartitionedTyings::iid() const +template +inline xt::xtensor MatrixPartitionedTyings::iid() const { return m_iid; } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitionedTyings::factorize() +template +inline void MatrixPartitionedTyings::factorize() { if (!m_factor) return; m_ACuu = m_Auu + m_Aud * m_Cdu + m_Cud * m_Adu + m_Cud * m_Add * m_Cdu; m_ACup = m_Aup + m_Aud * m_Cdp + m_Cud * m_Adp + m_Cud * m_Add * m_Cdp; // m_ACpu = m_Apu + m_Apd * m_Cdu + m_Cpd * m_Adu + m_Cpd * m_Add * m_Cdu; // m_ACpp = m_App + m_Apd * m_Cdp + m_Cpd * m_Adp + m_Cpd * m_Add * m_Cdp; m_solver.compute(m_ACuu); m_factor = false; } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitionedTyings::assemble(const xt::xtensor& elemmat) +template +inline void MatrixPartitionedTyings::assemble(const xt::xtensor& elemmat) { GOOSEFEM_ASSERT(elemmat.shape() ==\ std::decay_t::shape_type({m_nelem, m_nne*m_ndim, m_nne*m_ndim})); m_Tuu.clear(); m_Tup.clear(); m_Tpu.clear(); m_Tpp.clear(); m_Tud.clear(); m_Tpd.clear(); m_Tdu.clear(); m_Tdp.clear(); m_Tdd.clear(); for (size_t e = 0 ; e < m_nelem ; ++e) { for (size_t m = 0 ; m < m_nne ; ++m) { for (size_t i = 0 ; i < m_ndim ; ++i) { size_t di = m_dofs(m_conn(e,m),i); for (size_t n = 0 ; n < m_nne ; ++n) { for (size_t j = 0 ; j < m_ndim ; ++j) { size_t dj = m_dofs(m_conn(e,n),j); if (di < m_nnu and dj < m_nnu) m_Tuu.push_back(Eigen::Triplet(di ,dj ,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (di < m_nnu and dj < m_nni) m_Tup.push_back(Eigen::Triplet(di ,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (di < m_nnu) m_Tud.push_back(Eigen::Triplet(di ,dj-m_nni,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (di < m_nni and dj < m_nnu) m_Tpu.push_back(Eigen::Triplet(di-m_nnu,dj ,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (di < m_nni and dj < m_nni) m_Tpp.push_back(Eigen::Triplet(di-m_nnu,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (di < m_nni) m_Tpd.push_back(Eigen::Triplet(di-m_nnu,dj-m_nni,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (dj < m_nnu) m_Tdu.push_back(Eigen::Triplet(di-m_nni,dj ,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else if (dj < m_nni) m_Tdp.push_back(Eigen::Triplet(di-m_nni,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j))); else m_Tdd.push_back(Eigen::Triplet(di-m_nni,dj-m_nni,elemmat(e,m*m_ndim+i,n*m_ndim+j))); } } } } } m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end()); m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end()); m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end()); m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end()); m_Aud.setFromTriplets(m_Tud.begin(), m_Tud.end()); m_Apd.setFromTriplets(m_Tpd.begin(), m_Tpd.end()); m_Adu.setFromTriplets(m_Tdu.begin(), m_Tdu.end()); m_Adp.setFromTriplets(m_Tdp.begin(), m_Tdp.end()); m_Add.setFromTriplets(m_Tdd.begin(), m_Tdd.end()); m_factor = true; } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitionedTyings::solve( +template +inline void MatrixPartitionedTyings::solve( const xt::xtensor& b, xt::xtensor& x) { GOOSEFEM_ASSERT(b.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); GOOSEFEM_ASSERT(x.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); this->factorize(); Eigen::VectorXd B_u = this->asDofs_u(b); Eigen::VectorXd B_d = this->asDofs_d(b); Eigen::VectorXd X_p = this->asDofs_p(x); B_u += m_Cud * B_d; Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p)); Eigen::VectorXd X_d = m_Cdu * X_u + m_Cdp * X_p; #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) { for (size_t i = 0 ; i < m_ndim ; ++i) { if (m_dofs(m,i) < m_nnu) x(m,i) = X_u(m_dofs(m,i)); else if (m_dofs(m,i) >= m_nni) x(m,i) = X_d(m_dofs(m,i)-m_nni); } } } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitionedTyings::solve( +template +inline void MatrixPartitionedTyings::solve( const xt::xtensor& b, xt::xtensor& x) { GOOSEFEM_ASSERT(b.size() == m_ndof); GOOSEFEM_ASSERT(x.size() == m_ndof); this->factorize(); Eigen::VectorXd B_u = this->asDofs_u(b); Eigen::VectorXd B_d = this->asDofs_d(b); Eigen::VectorXd X_p = this->asDofs_p(x); Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p)); Eigen::VectorXd X_d = m_Cdu * X_u + m_Cdp * X_p; #pragma omp parallel for for (size_t d = 0 ; d < m_nnu ; ++d) x(m_iiu(d)) = X_u(d); #pragma omp parallel for for (size_t d = 0 ; d < m_nnd ; ++d) x(m_iid(d)) = X_d(d); } // ------------------------------------------------------------------------------------------------- -inline void MatrixPartitionedTyings::solve_u( +template +inline void MatrixPartitionedTyings::solve_u( const xt::xtensor& b_u, const xt::xtensor& b_d, const xt::xtensor& x_p, xt::xtensor& x_u) { GOOSEFEM_ASSERT(b_u.size() == m_nnu); GOOSEFEM_ASSERT(b_d.size() == m_nnd); GOOSEFEM_ASSERT(x_p.size() == m_nnp); GOOSEFEM_ASSERT(x_u.size() == m_nnu); this->factorize(); Eigen::VectorXd B_u(m_nnu,1); Eigen::VectorXd B_d(m_nnd,1); Eigen::VectorXd X_p(m_nnp,1); std::copy(b_u.begin(), b_u.end(), B_u.data()); std::copy(b_d.begin(), b_d.end(), B_d.data()); std::copy(x_p.begin(), x_p.end(), X_p.data()); Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p)); std::copy(X_u.data(), X_u.data()+m_nnu, x_u.begin()); } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_u( +template +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_u( const xt::xtensor& dofval) const { assert(dofval.size() == m_ndof); Eigen::VectorXd dofval_u(m_nnu,1); #pragma omp parallel for for (size_t d = 0 ; d < m_nnu ; ++d) dofval_u(d) = dofval(m_iiu(d)); return dofval_u; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_u( +template +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_u( const xt::xtensor& nodevec) const { assert(nodevec.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); Eigen::VectorXd dofval_u(m_nnu,1); #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) for (size_t i = 0 ; i < m_ndim ; ++i) if (m_dofs(m,i) < m_nnu) dofval_u(m_dofs(m,i)) = nodevec(m,i); return dofval_u; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_p( +template +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_p( const xt::xtensor& dofval) const { assert(dofval.size() == m_ndof); Eigen::VectorXd dofval_p(m_nnp,1); #pragma omp parallel for for (size_t d = 0 ; d < m_nnp ; ++d) dofval_p(d) = dofval(m_iip(d)); return dofval_p; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_p( +template +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_p( const xt::xtensor& nodevec) const { assert(nodevec.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); Eigen::VectorXd dofval_p(m_nnp,1); #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) for (size_t i = 0 ; i < m_ndim ; ++i) if (m_dofs(m,i) >= m_nnu and m_dofs(m,i) < m_nni) dofval_p(m_dofs(m,i)-m_nnu) = nodevec(m,i); return dofval_p; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d( +template +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d( const xt::xtensor& dofval) const { assert(dofval.size() == m_ndof); Eigen::VectorXd dofval_d(m_nnd,1); #pragma omp parallel for for (size_t d = 0 ; d < m_nnd ; ++d) dofval_d(d) = dofval(m_iip(d)); return dofval_d; } // ------------------------------------------------------------------------------------------------- -inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d( +template +inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d( const xt::xtensor& nodevec) const { assert(nodevec.shape() ==\ std::decay_t::shape_type({m_nnode, m_ndim})); Eigen::VectorXd dofval_d(m_nnd,1); #pragma omp parallel for for (size_t m = 0 ; m < m_nnode ; ++m) for (size_t i = 0 ; i < m_ndim ; ++i) if (m_dofs(m,i) >= m_nni) dofval_d(m_dofs(m,i)-m_nni) = nodevec(m,i); return dofval_d; } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/GooseFEM/ParaView.h b/include/GooseFEM/ParaView.h index 142c1d6..aea4da4 100644 --- a/include/GooseFEM/ParaView.h +++ b/include/GooseFEM/ParaView.h @@ -1,413 +1,418 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_PARAVIEW_H #define GOOSEFEM_PARAVIEW_H // ------------------------------------------------------------------------------------------------- #include "config.h" #include #ifndef GOOSEFEM_NO_HIGHFIVE #include #endif // ================================================================================================= namespace GooseFEM { namespace ParaView { namespace HDF5 { // ------------------------------------------------------------------------------------------------- inline std::string join(const std::vector& lines, const std::string& sep="\n"); // ------------------------------------------------------------------------------------------------- inline std::string indent(size_t n); // ------------------------------------------------------------------------------------------------- xt::xtensor as3d(const xt::xtensor &data); // ------------------------------------------------------------------------------------------------- // http://xdmf.org/index.php/XDMF_Model_and_Format enum class ElementType { Triangle, Quadrilateral, Hexahedron }; // ------------------------------------------------------------------------------------------------- enum class AttributeType { Cell, Node }; // ------------------------------------------------------------------------------------------------- inline ElementType convert(GooseFEM::Mesh::ElementType type); // ------------------------------------------------------------------------------------------------- inline std::string to_string(ElementType type); inline std::string to_string(AttributeType type); // ------------------------------------------------------------------------------------------------- class Connectivity { public: // constructors Connectivity() = default; #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) Connectivity( const H5Easy::File& data, const std::string& dataset, GooseFEM::Mesh::ElementType type); #endif #ifndef GOOSEFEM_NO_HIGHFIVE // (requires opening & reading the HDF5-file) Connectivity( const std::string& fname, const std::string& dataset, GooseFEM::Mesh::ElementType type); #endif Connectivity( const std::string& fname, const std::string& dataset, GooseFEM::Mesh::ElementType type, const std::vector& shape); #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) Connectivity( const H5Easy::File& data, const std::string& dataset, ElementType type); #endif #ifndef GOOSEFEM_NO_HIGHFIVE // (requires opening & reading the HDF5-file) Connectivity( const std::string& fname, const std::string& dataset, ElementType type); #endif Connectivity( const std::string& fname, const std::string& dataset, ElementType type, const std::vector& shape); // info size_t nelem() const; size_t nne() const; std::vector shape() const; std::string fname() const; // check #ifndef GOOSEFEM_NO_HIGHFIVE // (requires opening & reading the HDF5-file) void checkShape(); #endif // return lines of XDMF file (indented, indentation starts at 0) std::vector xdmf(size_t indent=4) const; private: // get shape from opened HDF5-file #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) void readShape(const H5Easy::File& data); #endif // back-end for constructors (allows single implementation for several overloads) #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) void init( const H5Easy::File& data, const std::string& dataset, ElementType type); #endif #ifndef GOOSEFEM_NO_HIGHFIVE // (requires opening & reading the HDF5-file) void init( const std::string& fname, const std::string& dataset, ElementType type); #endif void init( const std::string& fname, const std::string& dataset, ElementType type, const std::vector& shape); // internal data ElementType m_type; std::string m_fname; std::string m_dataset; std::vector m_shape; #ifndef GOOSEFEM_NO_HIGHFIVE bool m_verified=false; // if true: shape read from file, not need to check #endif }; // ------------------------------------------------------------------------------------------------- class Coordinates { public: // constructors Coordinates() = default; #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) Coordinates( const H5Easy::File& data, const std::string& dataset); #endif #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) Coordinates( const std::string& fname, const std::string& dataset); // (requires reading the HDF5) #endif Coordinates( const std::string& fname, const std::string& dataset, const std::vector& shape); // info size_t nnode() const; size_t ndim() const; std::vector shape() const; std::string fname() const; // check #ifndef GOOSEFEM_NO_HIGHFIVE // (requires opening & reading the HDF5-file) void checkShape(); #endif // return lines of XDMF file (indented, indentation starts at 0) std::vector xdmf(size_t indent=4) const; private: // get shape from opened HDF5-file #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) void readShape(const H5Easy::File& data); #endif // internal data std::string m_fname; std::string m_dataset; std::vector m_shape; #ifndef GOOSEFEM_NO_HIGHFIVE bool m_verified=false; // if true: shape read from file, not need to check #endif }; // ------------------------------------------------------------------------------------------------- class Attribute { public: // constructors Attribute() = default; #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) Attribute( const H5Easy::File& data, const std::string& dataset, const std::string& name, AttributeType type); #endif #ifndef GOOSEFEM_NO_HIGHFIVE // (requires opening & reading the HDF5-file) Attribute( const std::string& fname, const std::string& dataset, const std::string& name, AttributeType type); #endif Attribute( const std::string& fname, const std::string& dataset, const std::string& name, AttributeType type, const std::vector& shape); // info std::vector shape() const; std::string fname() const; // check #ifndef GOOSEFEM_NO_HIGHFIVE // (requires opening & reading the HDF5-file) void checkShape(); #endif std::vector xdmf(size_t indent=4) const; private: // get shape from opened HDF5-file #ifndef GOOSEFEM_NO_HIGHFIVE // (requires reading the HDF5-file) void readShape(const H5Easy::File& data); #endif // internal data AttributeType m_type; std::string m_fname; std::string m_dataset; std::string m_name; std::vector m_shape; #ifndef GOOSEFEM_NO_HIGHFIVE bool m_verified=false; // if true: shape read from file, not need to check #endif }; // ------------------------------------------------------------------------------------------------- class Mesh { public: // constructors Mesh() = default; Mesh( const Connectivity& conn, const Coordinates& coor); + // add attribute to this increment + + void push_back(const Attribute& data); + // return lines of XDMF file (indented, indentation starts at 0) std::vector xdmf(size_t indent=4) const; // write readable XDMF file void write(const std::string& fname, size_t n_indent=4) const; private: // internal data Connectivity m_conn; - Coordinates m_coor; + Coordinates m_coor; + std::vector m_attr; }; // ------------------------------------------------------------------------------------------------- class Increment { public: // constructors Increment() = default; Increment( const Connectivity& conn, const Coordinates& coor); Increment( const Connectivity& conn, const Coordinates& coor, const std::vector& attr); // add attribute to this increment void push_back(const Connectivity& data); void push_back(const Coordinates& data); void push_back(const Attribute& data); // return lines of XDMF file (indented, indentation starts at 0) std::vector xdmf(size_t indent=4) const; private: // internal data std::vector m_conn; - std::vector m_coor; - std::vector m_attr; + std::vector m_coor; + std::vector m_attr; }; // ------------------------------------------------------------------------------------------------- class TimeSeries { public: // constructors TimeSeries() = default; TimeSeries( const std::vector& data); // add increment void push_back(const Increment& data); // return lines of XDMF file (indented, indentation starts at 0) std::vector xdmf(size_t indent=4) const; // write readable XDMF file void write(const std::string& fname, size_t n_indent=4) const; private: // internal data std::vector m_data; }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #include "ParaView.hpp" // ================================================================================================= #endif diff --git a/include/GooseFEM/ParaView.hpp b/include/GooseFEM/ParaView.hpp index b3699cb..6483e21 100644 --- a/include/GooseFEM/ParaView.hpp +++ b/include/GooseFEM/ParaView.hpp @@ -1,770 +1,804 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_PARAVIEW_HPP #define GOOSEFEM_PARAVIEW_HPP // ------------------------------------------------------------------------------------------------- #include "ParaView.h" // ================================================================================================= namespace GooseFEM { namespace ParaView { namespace HDF5 { // ------------------------------------------------------------------------------------------------- inline std::string join(const std::vector& lines, const std::string &sep) { if (lines.size() == 1) return lines[0]; std::string out = ""; for (auto line: lines) { if (out.size() == 0) { out += line; continue; } if (line[0] == sep[0]) out += line; else if (out[out.size()-1] == sep[0]) out += line; else out += sep + line; } return out; } // ------------------------------------------------------------------------------------------- ------ inline std::string indent(size_t n) { std::string out = ""; for (size_t i = 0; i < n; ++i) out += " "; return out; } // ------------------------------------------------------------------------------------------------- xt::xtensor as3d(const xt::xtensor &data) { GOOSEFEM_ASSERT(data.shape(1) > 0 && data.shape(1) < 4) if (data.shape(1) == 3ul) return data; xt::xtensor out = xt::zeros({data.shape(0), 3ul}); if (data.shape(1) == 2ul) xt::view(out, xt::all(), xt::keep(0,1)) = data; if (data.shape(1) == 1ul) xt::view(out, xt::all(), xt::keep(0)) = data; return out; } // ------------------------------------------------------------------------------------------------- inline ElementType convert(GooseFEM::Mesh::ElementType type) { if (type == GooseFEM::Mesh::ElementType::Tri3) return ElementType::Triangle; if (type == GooseFEM::Mesh::ElementType::Quad4) return ElementType::Quadrilateral; if (type == GooseFEM::Mesh::ElementType::Hex8) return ElementType::Hexahedron; throw std::runtime_error("Unknown GooseFEM::Mesh::ElementType"); } // ------------------------------------------------------------------------------------------------- inline std::string to_string(ElementType type) { if (type == ElementType::Triangle) return "Triangle"; if (type == ElementType::Quadrilateral) return "Quadrilateral"; if (type == ElementType::Hexahedron) return "Hexahedron"; throw std::runtime_error("Unknown GooseFEM::ParaView::HDF5::ElementType"); } // ------------------------------------------------------------------------------------------------- inline std::string to_string(AttributeType type) { if (type == AttributeType::Cell) return "Cell"; if (type == AttributeType::Node) return "Node"; throw std::runtime_error("Unknown GooseFEM::ParaView::HDF5::AttributeType"); } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Connectivity::Connectivity( const H5Easy::File& data, const std::string& dataset, GooseFEM::Mesh::ElementType type) { init(data, dataset, convert(type)); } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Connectivity::Connectivity( const std::string& fname, const std::string& dataset, GooseFEM::Mesh::ElementType type) { init(fname, dataset, convert(type)); } #endif // ------------------------------------------------------------------------------------------------- inline Connectivity::Connectivity( const std::string& fname, const std::string& dataset, GooseFEM::Mesh::ElementType type, const std::vector& shape) { init(fname, dataset, convert(type), shape); } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Connectivity::Connectivity( const H5Easy::File& data, const std::string& dataset, ElementType type) { init(data, dataset, type); } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Connectivity::Connectivity( const std::string& fname, const std::string& dataset, ElementType type) { init(fname, dataset, type); } #endif // ------------------------------------------------------------------------------------------------- inline Connectivity::Connectivity( const std::string& fname, const std::string& dataset, ElementType type, const std::vector& shape) { init(fname, dataset, type, shape); } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Connectivity::init( const H5Easy::File& data, const std::string& dataset, ElementType type) { m_type = type; m_dataset = dataset; this->readShape(data); } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Connectivity::init( const std::string& fname, const std::string& dataset, ElementType type) { m_type = type; m_fname = fname; m_dataset = dataset; H5Easy::File data(m_fname, H5Easy::File::ReadOnly); this->readShape(data); } #endif // ------------------------------------------------------------------------------------------------- inline void Connectivity::init( const std::string& fname, const std::string& dataset, ElementType type, const std::vector& shape) { m_type = type; m_fname = fname; m_dataset = dataset; m_shape = shape; GOOSEFEM_ASSERT(m_shape.size() == 2); #ifdef GOOSEFEM_ENABLE_ASSERT - if (m_type == ElementType::Quadrilateral) + if (m_type == ElementType::Triangle) { + GOOSEFEM_ASSERT(m_shape[1] == 3); + } + else if (m_type == ElementType::Quadrilateral) { GOOSEFEM_ASSERT(m_shape[1] == 4); + } + else if (m_type == ElementType::Hexahedron) { + GOOSEFEM_ASSERT(m_shape[1] == 8); + } #endif } // ------------------------------------------------------------------------------------------------- inline size_t Connectivity::nelem() const { return m_shape[0]; } // ------------------------------------------------------------------------------------------------- inline size_t Connectivity::nne() const { return m_shape[1]; } // ------------------------------------------------------------------------------------------------- inline std::vector Connectivity::shape() const { return m_shape; } // ------------------------------------------------------------------------------------------------- inline std::string Connectivity::fname() const { return m_fname; } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Connectivity::checkShape() { if (m_verified) return; H5Easy::File data(m_fname, H5Easy::File::ReadOnly); GOOSEFEM_CHECK(m_shape == H5Easy::getShape(data, m_dataset)); m_verified = true; } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Connectivity::readShape(const H5Easy::File& data) { if (m_fname.size() == 0) m_fname = data.getName(); m_shape = H5Easy::getShape(data, m_dataset); GOOSEFEM_ASSERT(m_shape.size() == 2); #ifdef GOOSEFEM_ENABLE_ASSERT - if (m_type == ElementType::Quadrilateral) + if (m_type == ElementType::Triangle) { + GOOSEFEM_ASSERT(m_shape[1] == 3); + } + else if (m_type == ElementType::Quadrilateral) { GOOSEFEM_ASSERT(m_shape[1] == 4); + } + else if (m_type == ElementType::Hexahedron) { + GOOSEFEM_ASSERT(m_shape[1] == 8); + } #endif m_verified = true; } #endif // ------------------------------------------------------------------------------------------------- inline std::vector Connectivity::xdmf(size_t n_indent) const { std::vector out; out.push_back( ""); out.push_back(indent(n_indent) + "" + m_fname + ":" + m_dataset + ""); out.push_back(""); return out; } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Coordinates::Coordinates( const H5Easy::File& data, const std::string& dataset) : m_dataset(dataset) { this->readShape(data); } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Coordinates::Coordinates( const std::string& fname, const std::string& dataset) : m_fname(fname), m_dataset(dataset) { H5Easy::File data(m_fname, H5Easy::File::ReadOnly); this->readShape(data); } #endif // ------------------------------------------------------------------------------------------------- inline Coordinates::Coordinates( const std::string& fname, const std::string& dataset, const std::vector& shape) : m_fname(fname), m_dataset(dataset), m_shape(shape) { GOOSEFEM_ASSERT(m_shape.size() == 2); } // ------------------------------------------------------------------------------------------------- inline size_t Coordinates::nnode() const { return m_shape[0]; } // ------------------------------------------------------------------------------------------------- inline size_t Coordinates::ndim() const { return m_shape[1]; } // ------------------------------------------------------------------------------------------------- inline std::vector Coordinates::shape() const { return m_shape; } // ------------------------------------------------------------------------------------------------- inline std::string Coordinates::fname() const { return m_fname; } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Coordinates::checkShape() { if (m_verified) return; H5Easy::File data(m_fname, H5Easy::File::ReadOnly); GOOSEFEM_CHECK(m_shape == H5Easy::getShape(data, m_dataset)); m_verified = true; } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Coordinates::readShape(const H5Easy::File& data) { if (m_fname.size() == 0) m_fname = data.getName(); m_shape = H5Easy::getShape(data, m_dataset); GOOSEFEM_ASSERT(m_shape.size() == 2); m_verified = true; } #endif // ------------------------------------------------------------------------------------------------- inline std::vector Coordinates::xdmf(size_t n_indent) const { std::vector out; - out.push_back(""); + if (m_shape[1] == 1) + out.push_back(""); + else if (m_shape[1] == 2) + out.push_back(""); + else if (m_shape[1] == 3) + out.push_back(""); out.push_back(indent(n_indent) + "" + m_fname + ":" + m_dataset + ""); out.push_back(")"); return out; } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Attribute::Attribute( const H5Easy::File& data, const std::string& dataset, const std::string& name, AttributeType type) : m_type(type), m_dataset(dataset), m_name(name) { this->readShape(data); } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline Attribute::Attribute( const std::string& fname, const std::string& dataset, const std::string& name, AttributeType type) : m_type(type), m_fname(fname), m_dataset(dataset), m_name(name) { H5Easy::File data(m_fname, H5Easy::File::ReadOnly); this->readShape(data); } #endif // ------------------------------------------------------------------------------------------------- inline Attribute::Attribute( const std::string& fname, const std::string& dataset, const std::string& name, AttributeType type, const std::vector& shape) : m_type(type), m_fname(fname), m_dataset(dataset), m_name(name), m_shape(shape) { GOOSEFEM_ASSERT(m_shape.size() > 0); } // ------------------------------------------------------------------------------------------------- inline std::vector Attribute::shape() const { return m_shape; } // ------------------------------------------------------------------------------------------------- inline std::string Attribute::fname() const { return m_fname; } // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Attribute::checkShape() { if (m_verified) return; H5Easy::File data(m_fname, H5Easy::File::ReadOnly); GOOSEFEM_CHECK(m_shape == H5Easy::getShape(data, m_dataset)); m_verified = true; } #endif // ------------------------------------------------------------------------------------------------- #ifndef GOOSEFEM_NO_HIGHFIVE inline void Attribute::readShape(const H5Easy::File& data) { if (m_fname.size() == 0) m_fname = data.getName(); m_shape = H5Easy::getShape(data, m_dataset); GOOSEFEM_ASSERT(m_shape.size() > 0); m_verified = true; } #endif // ------------------------------------------------------------------------------------------------- inline std::vector Attribute::xdmf(size_t n_indent) const { GOOSEFEM_ASSERT(m_shape.size() > 0); GOOSEFEM_ASSERT(m_shape.size() < 3); std::vector out; if (m_shape.size() == 1) { out.push_back( ""); out.push_back(indent(n_indent) + "" + m_fname + ":" + m_dataset + ""); } else if (m_shape.size() == 2) { out.push_back( ""); out.push_back(indent(n_indent) + "" + m_fname + ":" + m_dataset + ""); } out.push_back(")"); return out; } // ------------------------------------------------------------------------------------------------- inline Mesh::Mesh(const Connectivity& conn, const Coordinates& coor) : m_conn(conn), m_coor(coor) { } // ------------------------------------------------------------------------------------------------- +inline void Mesh::push_back(const Attribute& data) +{ + m_attr.push_back(data); +} + +// ------------------------------------------------------------------------------------------------- + inline std::vector Mesh::xdmf(size_t n_indent) const { std::vector out; { std::vector lines = m_conn.xdmf(n_indent); for (auto& line: lines) out.push_back(line); } { std::vector lines = m_coor.xdmf(n_indent); for (auto& line: lines) out.push_back(line); } + for (auto& i: m_attr) + { + std::vector lines = i.xdmf(n_indent); + + for (auto& line: lines) + out.push_back(line); + } + return out; } // ------------------------------------------------------------------------------------------------- inline void Mesh::write(const std::string& fname, size_t n_indent) const { - TimeSeries({Increment(m_conn, m_coor)}).write(fname, n_indent); + TimeSeries({Increment(m_conn, m_coor, m_attr)}).write(fname, n_indent); } // ------------------------------------------------------------------------------------------------- inline Increment::Increment( const Connectivity& conn, const Coordinates& coor) { m_conn.push_back(conn); m_coor.push_back(coor); } // ------------------------------------------------------------------------------------------------- inline Increment::Increment( const Connectivity& conn, const Coordinates& coor, - const std::vector& attr) : + const std::vector& attr) : m_attr(attr) { m_conn.push_back(conn); m_coor.push_back(coor); } // ------------------------------------------------------------------------------------------------- inline void Increment::push_back(const Connectivity& data) { m_conn.push_back(data); } // ------------------------------------------------------------------------------------------------- inline void Increment::push_back(const Coordinates& data) { m_coor.push_back(data); } // ------------------------------------------------------------------------------------------------- inline void Increment::push_back(const Attribute& data) { m_attr.push_back(data); } // ------------------------------------------------------------------------------------------------- inline std::vector Increment::xdmf(size_t n_indent) const { std::vector out; for (auto& i: m_conn) { std::vector lines = i.xdmf(n_indent); for (auto& line: lines) out.push_back(line); } for (auto& i: m_coor) { std::vector lines = i.xdmf(n_indent); for (auto& line: lines) out.push_back(line); } for (auto& i: m_attr) { std::vector lines = i.xdmf(n_indent); for (auto& line: lines) out.push_back(line); } return out; } // ------------------------------------------------------------------------------------------------- inline TimeSeries::TimeSeries(const std::vector& data) : m_data(data) { } // ------------------------------------------------------------------------------------------------- inline void TimeSeries::push_back(const Increment& data) { m_data.push_back(data); } // ------------------------------------------------------------------------------------------------- inline std::vector TimeSeries::xdmf(size_t n_indent) const { std::vector out; for (size_t inc = 0; inc < m_data.size(); ++inc) { out.push_back( ""); out.push_back(indent(n_indent) + ")"); } return out; } // ------------------------------------------------------------------------------------------------- inline void TimeSeries::write(const std::string& fname, size_t n_indent) const { std::ofstream myfile; myfile.open(fname); myfile << "" << std::endl; myfile << indent(n_indent) + "" << std::endl; myfile << indent(n_indent * 2) + "" << std::endl; std::vector lines = this->xdmf(n_indent); for (auto& line: lines) myfile << indent(n_indent * 3) + line << std::endl; myfile << indent(n_indent * 2) + "" << std::endl; myfile << indent(n_indent) + "" << std::endl; myfile << "" << std::endl; myfile.close(); } // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/python/ParaView.hpp b/python/ParaView.hpp index c79390e..0d0a1d0 100644 --- a/python/ParaView.hpp +++ b/python/ParaView.hpp @@ -1,186 +1,190 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include #define GOOSEFEM_NO_HIGHFIVE #include "../include/GooseFEM/GooseFEM.h" #include "../include/GooseFEM/ParaView.h" // ================================================================================================= namespace py = pybind11; // ================================================================================================= void init_ParaViewHDF5(py::module &m) { m.def("as3d", &GooseFEM::ParaView::HDF5::as3d); // ------------------------------------------------------------------------------------------------- py::enum_(m, "ElementType", "ElementType") .value("Triangle" , GooseFEM::ParaView::HDF5::ElementType::Triangle ) .value("Quadrilateral", GooseFEM::ParaView::HDF5::ElementType::Quadrilateral) .value("Hexahedron" , GooseFEM::ParaView::HDF5::ElementType::Hexahedron ) .export_values(); // ------------------------------------------------------------------------------------------------- py::enum_(m, "AttributeType", "AttributeType") .value("Cell", GooseFEM::ParaView::HDF5::AttributeType::Cell) .value("Node", GooseFEM::ParaView::HDF5::AttributeType::Node) .export_values(); // ------------------------------------------------------------------------------------------------- py::class_(m, "Connectivity") .def(py::init< const std::string&, const std::string&, GooseFEM::ParaView::HDF5::ElementType, const std::vector&>(), "Connectivity", py::arg("fname"), py::arg("dataset"), py::arg("ElementType"), py::arg("shape") ) .def("nelem", &GooseFEM::ParaView::HDF5::Connectivity::nelem) .def("nne", &GooseFEM::ParaView::HDF5::Connectivity::nne) .def("shape", &GooseFEM::ParaView::HDF5::Connectivity::shape) .def("fname", &GooseFEM::ParaView::HDF5::Connectivity::fname) .def("xdmf", &GooseFEM::ParaView::HDF5::Connectivity::xdmf, py::arg("indent")=4) .def("__repr__", [](const GooseFEM::ParaView::HDF5::Connectivity &){ return ""; }); // ------------------------------------------------------------------------------------------------- py::class_(m, "Coordinates") .def(py::init< const std::string&, const std::string&, const std::vector&>(), "Coordinates", py::arg("fname"), py::arg("dataset"), py::arg("shape") ) .def("nnode", &GooseFEM::ParaView::HDF5::Coordinates::nnode) .def("ndim", &GooseFEM::ParaView::HDF5::Coordinates::ndim) .def("shape", &GooseFEM::ParaView::HDF5::Coordinates::shape) .def("fname", &GooseFEM::ParaView::HDF5::Coordinates::fname) .def("xdmf", &GooseFEM::ParaView::HDF5::Coordinates::xdmf, py::arg("indent")=4) .def("__repr__", [](const GooseFEM::ParaView::HDF5::Coordinates &){ return ""; }); // ------------------------------------------------------------------------------------------------- py::class_(m, "Attribute") .def(py::init< const std::string&, const std::string&, const std::string&, GooseFEM::ParaView::HDF5::AttributeType, const std::vector&>(), "Attribute", py::arg("fname"), py::arg("dataset"), py::arg("name"), py::arg("AttributeType"), py::arg("shape") ) .def("shape", &GooseFEM::ParaView::HDF5::Attribute::shape) .def("fname", &GooseFEM::ParaView::HDF5::Attribute::fname) .def("xdmf", &GooseFEM::ParaView::HDF5::Attribute::xdmf, py::arg("indent")=4) .def("__repr__", [](const GooseFEM::ParaView::HDF5::Attribute &){ return ""; }); // ------------------------------------------------------------------------------------------------- py::class_(m, "Mesh") .def(py::init< const GooseFEM::ParaView::HDF5::Connectivity&, const GooseFEM::ParaView::HDF5::Coordinates&>(), "Mesh", py::arg("conn"), py::arg("coor") ) + .def("push_back", + py::overload_cast( + &GooseFEM::ParaView::HDF5::Mesh::push_back)) + .def("xdmf", &GooseFEM::ParaView::HDF5::Mesh::xdmf, py::arg("indent")=4) .def("write", &GooseFEM::ParaView::HDF5::Mesh::write, py::arg("fname"), py::arg("indent")=4) .def("__repr__", [](const GooseFEM::ParaView::HDF5::Mesh &){ return ""; }); // ------------------------------------------------------------------------------------------------- py::class_(m, "Increment") .def(py::init< const GooseFEM::ParaView::HDF5::Connectivity&, const GooseFEM::ParaView::HDF5::Coordinates&>(), "Increment", py::arg("conn"), py::arg("coor") ) .def("push_back", py::overload_cast( &GooseFEM::ParaView::HDF5::Increment::push_back)) .def("push_back", py::overload_cast( &GooseFEM::ParaView::HDF5::Increment::push_back)) .def("push_back", py::overload_cast( &GooseFEM::ParaView::HDF5::Increment::push_back)) .def("xdmf", &GooseFEM::ParaView::HDF5::Increment::xdmf, py::arg("indent")=4) .def("__repr__", [](const GooseFEM::ParaView::HDF5::Increment &){ return ""; }); // ------------------------------------------------------------------------------------------------- py::class_(m, "TimeSeries") .def(py::init<>(), "TimeSeries" ) .def("push_back", &GooseFEM::ParaView::HDF5::TimeSeries::push_back) .def("xdmf", &GooseFEM::ParaView::HDF5::TimeSeries::xdmf, py::arg("indent")=4) .def("write", &GooseFEM::ParaView::HDF5::TimeSeries::write, py::arg("fname"), py::arg("indent")=4) .def("__repr__", [](const GooseFEM::ParaView::HDF5::TimeSeries &){ return ""; }); }