GooseFEM::Dynamics::Diagonal

Note

Source:

#include <GooseFEM/DynamicsDiagonalPeriodic.h>
#include <GooseFEM/DynamicsDiagonalSemiPeriodic.h>
#include <GooseFEM/DynamicsDiagonalLinearStrainQuad4.h>

Overview

Principle

The philosophy is to provide some structure to efficiently run a finite element simulation which remains customizable. Even more customization can be obtained by copy/pasting the source and modifying it to your need. The idea that is followed involves a hierarchy of three classes, whereby a class that is higher in hierarchy writes to some field of the class that is lower in hierarchy, runs a function, and reads from some field. In general:

  • Discretized system (GooseFEM::Dynamics::Diagonal::Periodic, GooseFEM::Dynamics::Diagonal::SemiPeriodic).

    Defines the discretized system, and assembles the (inverse) mass matrix, the displacement dependent forces, and the velocity dependent forces from element arrays that are provided by the element definition.

  • Element definition (GooseFEM::Dynamics::Diagonal::LinearStrain::Qaud4)

    Provides the element arrays by performing numerical quadrature. At the integration point the constitutive response is probed from the quadrature point definition.

  • Quadrature point definition

    This class is not provided, and should be provided by the user.

Example

A simple example is:

[source: examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/main.cpp]

#include <Eigen/Eigen>
#include <cppmat/cppmat.h>
#include <GooseFEM/GooseFEM.h>
#include <GooseMaterial/AmorphousSolid/LinearStrain/Elastic/Cartesian2d.h>

// -------------------------------------------------------------------

using     MatS = GooseFEM::MatS;
using     MatD = GooseFEM::MatD;
using     ColD = GooseFEM::ColD;

using     T2   = cppmat::cartesian2d::tensor2 <double>;
using     T2s  = cppmat::cartesian2d::tensor2s<double>;

namespace GM   = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d;

// ===================================================================

class Quadrature
{
public:
  T2s eps, epsdot, sig;

  size_t nhard;
  GM::Material hard, soft;

  double Ebar, Vbar;

  Quadrature(size_t nhard);

  double density             (size_t elem, size_t k, double V);
  void   stressStrain        (size_t elem, size_t k, double V);
  void   stressStrainRate    (size_t elem, size_t k, double V);
  void   stressStrainPost    (size_t elem, size_t k, double V);
  void   stressStrainRatePost(size_t elem, size_t k, double V);
};

// -------------------------------------------------------------------

Quadrature::Quadrature(size_t _nhard)
{
  nhard    = _nhard;
  hard     = GM::Material(100.,10.);
  soft     = GM::Material(100., 1.);
}

// -------------------------------------------------------------------

double Quadrature::density(size_t elem, size_t k, double V)
{
  return 1.0;
}
// -------------------------------------------------------------------

void Quadrature::stressStrain(size_t elem, size_t k, double V)
{
  if ( elem < nhard ) sig = hard.stress(eps);
  else                sig = soft.stress(eps);
}
// -------------------------------------------------------------------

void Quadrature::stressStrainRate(size_t elem, size_t k, double V)
{
}
// -------------------------------------------------------------------

void Quadrature::stressStrainPost(size_t elem, size_t k, double V)
{
  Vbar += V;

  if ( elem < nhard ) Ebar += hard.energy(eps) * V;
  else                Ebar += soft.energy(eps) * V;
}

// -------------------------------------------------------------------

void Quadrature::stressStrainRatePost(size_t elem, size_t k, double V)
{
}

// ===================================================================

int main()
{
  // class which provides the mesh
  GooseFEM::Mesh::Quad4::Regular mesh(40,40,1.);

  // class which provides the constitutive response at each quadrature point
  auto  quadrature = std::make_shared<Quadrature>(40*40/4);

  // class which provides the response of each element
  using Elem = GooseFEM::Dynamics::Diagonal::LinearStrain::Quad4<Quadrature>;
  auto  elem = std::make_shared<Elem>(quadrature);

  // class which provides the system and an increment
  GooseFEM::Dynamics::Diagonal::Periodic<Elem> sim(
    elem,
    mesh.coor(),
    mesh.conn(),
    mesh.dofsPeriodic(),
    1.e-2,
    0.0
  );

  // loop over increments
  for ( ... )
  {
    // - set displacement of fixed DOFs
    ...

    // - compute time increment
    sim.Verlet();

    // - post-process
    quadrature->Ebar = 0.0;
    quadrature->Vbar = 0.0;

    sim.post();

    ...
  }

  return 0;
}

Pseudo-code

What is happening inside Verlet is evaluating the forces (and the mass matrix), and updating the displacements by solving the system. In pseudo-code:

  • Mass matrix:

    sim.computeMinv():
    {
      for e in elements:
    
        sim->elem->xe(i,j) = ...
        sim->elem->ue(i,j) = ...
    
        sim->elem->computeM(e):
        {
          for k in integration-points:
    
            sim->elem->M(...,...) += ... * sim->elem->quad->density(e,k,V)
        }
    
        M(...) += sim->elem->M(i,i)
    }
    
  • Displacement dependent force:

    sim.computeFu():
    {
      for e in elements:
    
        sim->elem->xe(i,j) = ...
        sim->elem->ue(i,j) = ...
    
        sim->elem->computeFu(e):
        {
          for k in integration-points:
    
            sim->elem->quad->eps(i,j) = ...
    
            sim->elem->quad->stressStrain(e,k,V)
    
            sim->elem->fu(...) += ... * sim->elem->quad->sig(i,j)
        }
    
        Fu(...) += sim->elem->fu(i)
    }
    
  • Velocity dependent force:

    sim.computeFv():
    {
      for e in elements:
    
        sim->elem->xe(i,j) = ...
        sim->elem->ue(i,j) = ...
        sim->elem->ve(i,j) = ...
    
        sim->elem->computeFv(e):
        {
          for k in integration-points:
    
            sim->elem->quad->epsdot(i,j) = ...
    
            sim->elem->quad->stressStrainRate(e,k,V)
    
            sim->elem->fv(...) += ... * sim->elem->quad->sig(i,j)
        }
    
        Fv(...) += sim->elem->fu(i)
    }
    

Signature

From this it is clear that:

  • GooseFEM::Dynamics::Diagonal::Periodic requires the following minimal signature from GooseFEM::Dynamics::Diagonal::LinearStrain::Qaud4:

    class Element
    {
    public:
      matrix M;                    // should have operator(i,j)
      column fu, fv;               // should have operator(i)
      matrix xe, ue, ve;           // should have operator(i,j)
    
      void computeM (size_t elem); // mass matrix                     <- quad->density
      void computeFu(size_t elem); // displacement dependent forces   <- quad->stressStrain
      void computeFv(size_t elem); // displacement dependent forces   <- quad->stressStrainRate
      void post     (size_t elem); // post-process                    <- quad->stressStrain(Rate)
    }
    
  • GooseFEM::Dynamics::Diagonal::LinearStrain::Qaud4 requires the minimal signature from Quadrature

    class Quadrature
    {
    public:
      tensor eps, epsdot, sig;     // should have operator(i,j)
    
      double density             (size_t elem, size_t k, double V);
      void   stressStrain        (size_t elem, size_t k, double V);
      void   stressStrainRate    (size_t elem, size_t k, double V);
      void   stressStrainPost    (size_t elem, size_t k, double V);
      void   stressStrainRatePost(size_t elem, size_t k, double V);
    }