mpoints2d.h

crack modeling with xfem

/* 2008 (c) Dorival M. Pedroso */

#ifndef MPM_MPOINTS2D_H
#define MPM_MPOINTS2D_H

// STL
#include <cmath>  // for sqrt, pow, etc.
#include <cfloat> // for DBL_EPSILON

// Local
#include "defs.h"
#include "array.h"
#include "grid2d.h"
#include "model2d.h"

class MPoints2D
{
public:
	// Constants
	static double MINMASS; ///< Mininum mass of the point/particle

	/** Alternative constructor. nPCell=number of points per cell */
	MPoints2D (size_t nPcell, Grid2D * G,
	           ptIsPointInGeom pIsPointInGeom,
	           ptDensity       pDensity,
	           ptAllocateModel pAllocateModel,
	           ptIniVelocity   pIniVelocity,
	           ptHasTraction   pHasTraction);

	// Methods
	size_t nPoints      () const { return _p.Size();     } ///< Return the number of material points
	size_t nPointsOnBry () const { return _ipbry.Size(); } ///< Return the number of material points on boundary
	void   SetGrid      (Grid2D * G) { _g = G; }           ///< Set access (read/write) to the grid
	bool   TimeUpdate   (double Dt, Vector3D const & b, double LdM, double & DsE); ///< Advance one step in time, considering this new b-body forces. Returns DsE, the increment of strain Energy.
	void   ResetDeform  ();                                ///< Clear strains and reset deformation gradient
	void   ReInit       () { _initialize(); }              ///< Reinitialize points (during initial refinement)

	// Access methods
	Array<Vector3D>       & Ps   ()               { return _p;             } ///< Returns access to all material points (coordinates) (write)
	Array<Vector3D> const & P    ()         const { return _p;             } ///< Returns access to all material points (coordinates) (read)
	Array<Vector3D> const & v    ()         const { return _v;             } ///< Returns access to all mat points velocities
	Vector3D        const & P    (size_t i) const { return _p[i];          } ///< Returns access to a material point (read)
	Vector3D        const & B    (size_t i) const { return _p[_ipbry[i]];  } ///< Returns access to a material point which is on boundary (read)
	Vector3D        const & l    (size_t i) const { return _l[i];          } ///< Length of the point/particle
	double          const & m    (size_t i) const { return _m[i];          } ///< Mass of the point/particle
	Vector3D        const & f    (size_t i) const { return _f[i];          } ///< Force on the point/particle
	Vector3D        const & t    (size_t i) const { return _t[i];          } ///< Tractions on point/particle
	Vector3D        const & u    (size_t i) const { return _u[i];          } ///< Displacement of the point/particle
	Vector3D        const & v    (size_t i) const { return _v[i];          } ///< Velocity of the point/particle
	STensor2        const & s    (size_t i) const { return _mdl[i]->Sig(); } ///< Stress on the point/particle
	STensor2        const & e    (size_t i) const { return _mdl[i]->Eps(); } ///< Strain on the point/particle
	bool                    HasT (size_t i) const { return _has_t[i];      } ///< Is mat point on boundary with applied tractions ?

private:
	// Data
	bool              _mpm;              ///< Use original MPM shape functions
	bool              _usf;              ///< Update stress first ?
	int               _npcell;           ///< Number of points per cell (initial only)
	Grid2D          * _g;                ///< Access to the grid
	ptIsPointInGeom   _is_point_in_geom; ///< Is point in geometry ?
	ptDensity         _density;          ///< Density function
	ptAllocateModel   _allocate_model;   ///< Allocate constitutive model
	ptIniVelocity     _ini_velocity;     ///< Initial velocity
	ptHasTraction     _has_traction;     ///< Has traction BC ?

	// Geometry
	Array<Vector3D> _p;     ///< Mat points current positions (x,y coords)
	Array<Vector3D> _p_ini; ///< Mat points initial positions (x,y coords)
	Array<int>      _ipbry; ///< Index of the mat points on boundary

	// Data (size == _np)
	Array<Vector3D>        _l;   ///< Mat points half-lengths (Vp=4*lx*ly)
	Array<double>          _m;   ///< Mat points masses
	Array<Vector3D>        _f;   ///< Mat points force (x,y comps)
	Array<Vector3D>        _u;   ///< Mat points displacements (x,y comps)
	Array<Vector3D>        _v;   ///< Mat points velocities (x,y comps)
	Array<ATensor2>        _F;   ///< Mat points def grad
	Array<size_t>          _lbn; ///< Left-bottom grid-nodes in the left-bottom cell of the area where each point contributes to
	Array<ShapeAndGrads2D> _sg;  ///< Shape and gradient functions
	Array<Model2D *>       _mdl; ///< Constitutive models for each material point

	// Boundary conditions (size == np)
	Array<bool>     _has_t; ///< Is mat point on boundary with applied traction ?
	Array<Vector3D> _t;     ///< Mat points tractions

	// Private methods
	void   _initialize    ();                                 ///< Initialize/generate mat points
	void   _stress_update (double Dt, double & DsE);          ///< Update stresses, returning DsE (increment of strain Energy)
	double _mpm_snp       (int n, int p, int icomp)    const; ///< 1D Shape functions (MPM)
	double _mpm_gnp       (int n, int p, int icomp)    const; ///< 1D Derivative of shape functions (MPM)
	double _snp           (int n, int p, int icomp)    const; ///< 1D Shape functions (GMPM)
	double _gnp           (int n, int p, int icomp)    const; ///< 1D Derivative of shape functions (GMPM)
	void   _Snp           (int n, int p, double & S)   const; ///< 2D shape functions
	void   _Gnp           (int n, int p, Vector3D & G) const; ///< 2D Derivative of shape functions

}; // class MPoints2D

double MPoints2D::MINMASS = sqrt(DBL_EPSILON);


/////////////////////////////////////////////////////////////////////////////////////////// Implementation /////


/* public */

inline MPoints2D::MPoints2D(size_t nPcell, Grid2D * G, ptIsPointInGeom pIsPointInGeom, ptDensity pDensity, ptAllocateModel pAllocateModel, ptIniVelocity pIniVelocity, ptHasTraction pHasTraction)
	: _mpm              (true),
	  _usf              (true),
	  _npcell           (nPcell),
	  _g                (G),
	  _is_point_in_geom (pIsPointInGeom),
	  _density          (pDensity),
	  _allocate_model   (pAllocateModel),
	  _ini_velocity     (pIniVelocity),
	  _has_traction     (pHasTraction)
{
	_initialize();
}

inline bool MPoints2D::TimeUpdate(double Dt, Vector3D const & b, double LdM, double & DsE)
{
	// Check if access to the grid was set
	if (_g==NULL) return false; // failed

	// 1) Calculate interpolation functions
	for (int p=0; p<_p.Size(); ++p)
	{
		int k = 0;
		for (int i=0; i<4; ++i)
		for (int j=0; j<4; ++j)
		{
			int n = _lbn[p]+((i*_g->nCol())+j); // grid node number
			if (n>=_g->nNodes())
			{
				std::cout << "MPoints2D::TimeUpdate: Material points can not be at the border of the grid." << std::endl;
				return false;
			}
			_Snp(n,p, _sg[p].S[k]); // calc interpolation function
			_Gnp(n,p, _sg[p].G[k]); // calc deriv. interp. function
			k++;
		}
	}

	// 2) Clear grid
	_g->ClearState();

	// 3) Initialize grid state from points (mass and momentum)
	for (int p=0; p<_p.Size(); ++p)
	{
		int k = 0;
		for (int i=0; i<4; ++i)
		for (int j=0; j<4; ++j)
		{
			int n = _lbn[p]+((i*_g->nCol())+j);  // grid node number
			_g->m(n) += _sg[p].S[k]*_m[p];       // node mass
			_g->q(n) += _sg[p].S[k]*_m[p]*_v[p]; // node momentum
			k++;
		}
	}

	// 4) Update strains and stress (USF - update stress first)
	if (_usf) _stress_update (Dt, DsE);
	
	// 5) Calculate internal and external forces
	for (int p=0; p<_p.Size(); ++p)
	{
		double Vp = 4.0*_l[p](0)*_l[p](1); // particle volume
		int k = 0;
		for (int i=0; i<4; ++i)
		for (int j=0; j<4; ++j)
		{
		   	// Grid node number
			int n = _lbn[p]+((i*_g->nCol())+j);

			// Internal forces
			if (_mdl[p]!=NULL)
				ScDotUp (Vp, _sg[p].G[k], _mdl[p]->Sig(),  _g->fi(n)); // fi += Vp * G . Sig

			// Body forces
			_g->fe(n) += _sg[p].S[k]*_m[p]*b;

			// Traction forces
			if (_has_t[p])
				_g->fe(n) += LdM*_sg[p].S[k]*_t[p];
			/*
			{
				Vector3D tmp;
				Vector3D nu; nu = _t[p]/Norm(_t[p]);
				Dot (_mdl[p]->Sig(), nu, tmp);
				_g->fe(n) += LdM*(_sg[p].S[k]*tmp - alp*(_u[p]-_t[p]));
			}
			*/

			// Next node contribution
			k++;
		}
	}

	// 6) Calculate nodes momenta rate
	for (int n=0; n<_g->nNodes(); ++n) _g->dqdt(n) = _g->fe(n) - _g->fi(n);

	// 7) Set fixed nodes (momentum)
	_g->FixNodesDqDt();

	// 8) Update nodes momenta
	for (int n=0; n<_g->nNodes(); ++n) _g->q(n) += _g->dqdt(n)*Dt;

	// 9) Update strains and stress (USL - update stress last)
	if (!_usf) _stress_update (Dt, DsE);

	// 10) Update material points (position, velocity, def. grad. and volume)
	for (int p=0; p<_p.Size(); ++p)
	{
		// Update p, v, f, and F
		ATensor2 newF; newF = 0.0;
		int k = 0;
		_f[p] = 0.0;
		for (int i=0; i<4; ++i)
		for (int j=0; j<4; ++j)
		{
			int n = _lbn[p]+((i*_g->nCol())+j); // grid node number
			if (_g->m(n)>MINMASS)
			{
				_p[p] += (Dt*_sg[p].S[k]/_g->m(n))*_g->q(n);    // point position
				_v[p] += (Dt*_sg[p].S[k]/_g->m(n))*_g->dqdt(n); // point velocity
				_f[p] += _sg[p].S[k]*_g->fe(n);                 // point force