fem_simulator.hpp

Go to the documentation of this file.

#pragma once
 
#include "tetmesh.hpp"
 
namespace netdem {
 
// use tetradron element with 4 nodes
class FEMSimulator {
public:
  double neo_k{6.94e5};
 
  double neo_mu{5.21e5};
 
  double density{500.0};
 
  double damp_coef{0.7};
 
  Vec3d gravity_coef{0.0, 0.0, 0.0};
 
  double timestep{1.0e-4};
 
  VecXT<Vec3d> nodes;
  VecXT<Vec4i> elements;
 
  VecXT<Vec3d> nodes_ref;
 
  VecXT<Vec3i> bound_facets;
  VecXT<int> bound_nodes;
 
  VecXT<double> elemental_vol;
  VecXT<VecNT<double, 6>> elemental_stress;
 
  VecXT<double> nodal_vols;
  VecXT<Vec3d> nodal_vels;
 
  VecXT<Vec3d> bc_facet_forces;
 
  VecXT<VecNT<double, 6>> bc_nodal_vels;
 
public:
  FEMSimulator();
 
  void SetMesh(TetMesh const &tetmesh);
 
  void Init();
 
  void SetBCNodalVelocity(int nid, double vx, double vy, double vz,
                          bool use_prescribed_vx, bool use_prescribed_vy,
                          bool use_prescribed_vz);
 
  void AddBCFacetForce(int bc_fid, double fx, double fy, double fz);
 
  void SetNodalVels(double v_x, double v_y, double v_z);
 
  void ClearBoundaryCondition();
 
  void Solve(double dt);
 
  VecXT<Vec3d> GetNodalPositions(VecXT<int> nids);
 
  VecXT<Vec3d> GetNodalDisps(VecXT<int> nids);
 
  VecXT<Vec3d> GetNodalVels(VecXT<int> nids);
 
  void SaveAsVTK(std::string const &file_name);
 
protected:
  VecXT<Vec3d> nodal_forces_int;
 
  VecXT<Vec3d> nodal_forces_ext;
 
  VecXT<Vec3d> nodal_vels_ave;
 
  void Advance(double dt);
 
  void InitInitialCondition();
 
  double GetElementVolume(Vec3d const &v0, Vec3d const &v1, Vec3d const &v2,
                          Vec3d const &v3);
 
  Mat3d GetDeformationGradient(Vec3d const &v0_new, Vec3d const &v1_new,
                               Vec3d const &v2_new, Vec3d const &v3_new,
                               Vec3d const &v0_ref, Vec3d const &v1_ref,
                               Vec3d const &v2_ref, Vec3d const &v3_ref);
 
  Mat3d GetCauchyStress(Mat3d const &def_grad);
 
  MatNd<4, 3> GetInternalForces(Mat3d const &cauchy_stress, Vec3d const &v0,
                                Vec3d const &v1, Vec3d const &v2,
                                Vec3d const &v3);
};
 
} // namespace netdem