Page MenuHomec4science
Files F91168032

main.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Nov 8, 14:25

main.cpp
View Options

#include <Eigen/Eigen>
#include <GooseFEM/GooseFEM.h>
#include <GooseFEM/ParaView.h>
#include <GMatElastic/Cartesian3d.h>
#include <highfive/H5Easy.hpp>
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<double,2> coor = mesh.coor();
xt::xtensor<size_t,2> conn = mesh.conn();
xt::xtensor<size_t,2> dofs = mesh.dofs();
xt::xtensor<size_t,2> 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<size_t,2> control_dofs = control.controlDofs();
xt::xtensor<size_t,1> 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<size_t,1> 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<double,2> disp = xt::zeros<double>(coor.shape()); // nodal displacement
xt::xtensor<double,2> fint = xt::zeros<double>(coor.shape()); // internal force
xt::xtensor<double,2> fext = xt::zeros<double>(coor.shape()); // external force
xt::xtensor<double,2> fres = xt::zeros<double>(coor.shape()); // residual force
// element vectors / matrix
xt::xtensor<double,3> ue = xt::empty<double>({nelem, nne, ndim});
xt::xtensor<double,3> fe = xt::empty<double>({nelem, nne, ndim});
xt::xtensor<double,3> Ke = xt::empty<double>({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<size_t,1> ehard = xt::ravel(xt::view(elmat, xt::range(0,2*10), xt::range(0,2*10)));
xt::xtensor<size_t,2> Ihard = xt::zeros<size_t>({nelem, nip});
xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul;
xt::xtensor<size_t,2> Isoft = xt::ones<size_t>({nelem, nip}) - Ihard;
mat.setElastic(Isoft, 10.0, 1.0);
mat.setElastic(Ihard, 10.0, 10.0);
// solve
// -----
// allocate tensors
xt::xtensor<double,4> Eps = xt::empty<double>({nelem, nip, tdim, tdim});
xt::xtensor<double,4> Sig = xt::empty<double>({nelem, nip, tdim, tdim});
xt::xtensor<double,6> C = xt::empty<double>({nelem, nip, tdim, tdim, tdim, tdim});
// allocate system matrix
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<double,4> 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<double,3> Sigelem = xt::average(Sig, dV, {1});
xt::xtensor<double,3> Epselem = xt::average(Eps, dV, {1});
// - macroscopic strain and stress
xt::xtensor_fixed<double, xt::xshape<3,3>> Sigbar = xt::average(Sig, dV, {0,1});
xt::xtensor_fixed<double, xt::xshape<3,3>> 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;
}

Event Timeline

c4science · Help