grid2d.h

crack modeling with xfem

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

#ifndef MPM_GRID2D_H
#define MPM_GRID2D_H

// Local
#include "defs.h"

class Grid2D
{
public:
	/* Constructor. */
	Grid2D (double xMin, double yMin, size_t nRow, size_t nCol, double Lx, double Ly, ptIsFixed pIsFixed);

	// Methods
	size_t nNodes        () const { return _n.Size(); }   ///< Return the number of nodes = nRow*nCol
	size_t nCells        () const { return _c.Size(); }   ///< Return the number of cells = (nRow-1)*(nCol-1)
	bool   IsInside      (Vector3D const & P);            ///< Check if a given point (P) is inside the grid
	void   ClearState    ();                              ///< Reset all state values in the grid, such as mass, velocity, momentum, etc.
	void   SetFixed      (size_t i, FixType Type);        ///< Set fixed node # i
	void   SetLoading    (size_t i, Vector3D const & Fe); ///< Set external force (loading) on node # i
	void   FixNodesVel   ();                              ///< Fix nodes by clearing v (velocities)
	void   FixNodesDqDt  ();                              ///< Fix nodes by clearing dqdt (momenta rate)
	void   ApplyLoadings (double M);                      ///< Apply loading using data on loaded nodes and scaling by the multiplier M
	void   Refine        ();                              ///< Refine the grid (bi-section)

	// Access methods
	double                  xMin ()             const { return _xmin;          } ///< Returns Min x
	double                  yMin ()             const { return _ymin;          } ///< Returns Min y
	double                  xMax ()             const { return _xmax;          } ///< Returns Max x
	double                  yMax ()             const { return _ymax;          } ///< Returns Max y
	size_t                  nRow ()             const { return _nrow;          } ///< Returns Num of rows
	size_t                  nCol ()             const { return _ncol;          } ///< Returns Num of columns
	double                  L    (size_t iComp) const { return _L(iComp);      } ///< Returns Cell-length 0=>x, 1=>y
	Vector3D const &        L    ()             const { return _L;             } ///< Returns Cell-length 0=>x, 1=>y
	Array<Vector3D> const & N    ()             const { return _n;             } ///< Access all nodes (coordinates)
	Vector3D const &        N    (size_t i)     const { return _n[i];          } ///< Returns access to a node position (read)
	Connec2D const &        C    (size_t i)     const { return _c[i];          } ///< Returns access to a connectivity (read)
	size_t                  nFix ()             const { return _fixed.Size();  } ///< Returns the number of fixed nodes
	Vector3D const &        FixN (size_t i)     const { return _n[_fixed[i]];  } ///< Access a fixed node (coordinates)
	FixType                 FixT (size_t i)     const { return _fixtype[i];    } ///< Access the type of a fixed node
	size_t                  nLd  ()             const { return _loaded.Size(); } ///< Returns the number of loaded nodes
	Vector3D const &        LdN  (size_t i)     const { return _n[_loaded[i]]; } ///< Access a loaded node (coordinates)
	Vector3D const &        Ld   (size_t i)     const { return _loads[i];      } ///< Access a loading

	// Set methods
	double   & m    (size_t i) { return _m   [i]; } ///< Nodes masses
	Vector3D & q    (size_t i) { return _q   [i]; } ///< Nodes momenta (x,y comps)
	Vector3D & v    (size_t i) { return _v   [i]; } ///< Nodes velocities (x,y comps)
	Vector3D & fe   (size_t i) { return _fe  [i]; } ///< Nodes external forces (x,y comps)
	Vector3D & fi   (size_t i) { return _fi  [i]; } ///< Nodes internal forces (x,y comps)
	Vector3D & dqdt (size_t i) { return _dqdt[i]; } ///< Rate of nodes momenta (x,y comps)

	// Get methods
	double   const & m    (size_t i) const { return _m   [i]; } ///< Nodes masses
	Vector3D const & q    (size_t i) const { return _q   [i]; } ///< Nodes momenta (x,y comps)
	Vector3D const & v    (size_t i) const { return _v   [i]; } ///< Nodes velocities (x,y comps)
	Vector3D const & fe   (size_t i) const { return _fe  [i]; } ///< Nodes external forces (x,y comps)
	Vector3D const & fi   (size_t i) const { return _fi  [i]; } ///< Nodes internal forces (x,y comps)
	Vector3D const & dqdt (size_t i) const { return _dqdt[i]; } ///< Rate of nodes momenta (x,y comps)

private:
	// Input
	double     _xmin;    ///< Min x
	double     _ymin;    ///< Min y
	size_t     _nrow;    ///< Num of rows
	size_t     _ncol;    ///< Num of columns
	Vector3D  _L;        ///< Cell-length 0=>x, 1=>y
	ptIsFixed _is_fixed; ///< Point to function which detects if a node is fixed

	// Geometry
	double          _xmax; ///< Max x
	double          _ymax; ///< Max y
	Array<Vector3D> _n;    ///< Nodes (x,y coords) (size=nRow*nCol)
	Array<Connec2D> _c;    ///< Connectivities (size=(nRow-1)*(nCol-1))

	// State values
	Array<double>   _m;    ///< Nodes masses
	Array<Vector3D> _q;    ///< Nodes momenta (x,y comps)
	Array<Vector3D> _v;    ///< Nodes velocities (x,y comps)
	Array<Vector3D> _fe;   ///< Nodes external forces (x,y comps)
	Array<Vector3D> _fi;   ///< Nodes internal forces (x,y comps)
	Array<Vector3D> _dqdt; ///< Rate of nodes momenta (x,y comps)

	// Fixed nodes
	Array<size_t>  _fixed;   ///< Indexes of the fixed nodes
	Array<FixType> _fixtype; ///< Type of the fixing

	// Loaded nodes
	Array<size_t>   _loaded; ///< Indexes of the loaded nodes
	Array<Vector3D> _loads;  ///< Loading values (applied/external forces)

	// Private methods
	void _refine(); ///< Initialize/refine the grid

}; // class Grid2D


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


/* public */

inline Grid2D::Grid2D(double xMin, double yMin, size_t nRow, size_t nCol, double Lx, double Ly, ptIsFixed pIsFixed)
	: _xmin     (xMin),
	  _ymin     (yMin),
	  _nrow     (nRow),
	  _ncol     (nCol),
	  _is_fixed (pIsFixed)
{
	// Save lengths
	_L = Lx, Ly, 0.0;

	// Do generate grid
	_refine ();
}

inline bool Grid2D::IsInside(Vector3D const & P)
{
	if (P(0)>=_xmin && P(0)<=_xmax && P(1)>=_ymin && P(1)<=_ymax) return true;
	else                                                          return false;
	
}

inline void Grid2D::ClearState()
{
	for (int i=0; i<_n.Size(); ++i)
	{
		_m   [i] = 0.0; // nodes masses
		_q   [i] = 0.0; // nodes momenta
		_v   [i] = 0.0; // nodes velocities
		_fe  [i] = 0.0; // nodes external forces
		_fi  [i] = 0.0; // nodes internal forces
		_dqdt[i] = 0.0; // nodes momenta rate
	}
}

inline void Grid2D::SetFixed(size_t i, FixType Type)
{
	long k = _fixed.Find(i);
	if (k<0) // not added yet
	{
		_fixed  .Push(i);
		_fixtype.Push(Type);
	}
	else
	{
		_fixed  [k] = i;
		_fixtype[k] = Type;
	}
}

inline void Grid2D::SetLoading(size_t i, Vector3D const & Fe)
{
	long k = _loaded.Find(i);
	if (k<0) // not added yet
	{
		_loaded.Push(i);
		_loads .Push(Fe);
	}
	else
	{
		_loaded[k] = i;
		_loads [k] = Fe;
	}
}

inline void Grid2D::FixNodesVel()
{
	for (int i=0; i<_fixed.Size(); ++i)
	{
		switch (_fixtype[i])
		{
			case FIX_X   : { _v[_fixed[i]](0)=0.0;                                             break; }
			case FIX_Y   : { _v[_fixed[i]](1)=0.0;                                             break; }
			case FIX_Z   : { _v[_fixed[i]](2)=0.0;                                             break; }
			case FIX_XY  : { _v[_fixed[i]](0)=0.0; _v[_fixed[i]](1)=0.0;                       break; }
			case FIX_YZ  : { _v[_fixed[i]](1)=0.0; _v[_fixed[i]](2)=0.0;                       break; }
			case FIX_ZX  : { _v[_fixed[i]](2)=0.0; _v[_fixed[i]](0)=0.0;                       break; }
			case FIX_XYZ : { _v[_fixed[i]](0)=0.0; _v[_fixed[i]](1)=0.0; _v[_fixed[i]](2)=0.0; break; }
		}
	}
}

inline void Grid2D::FixNodesDqDt()
{
	for (size_t i=0; i<_fixed.Size(); ++i)
	{
		switch (_fixtype[i])
		{
			case FIX_X   : { _dqdt[_fixed[i]](0)=0.0;                                                   break; }
			case FIX_Y   : { _dqdt[_fixed[i]](1)=0.0;                                                   break; }
			case FIX_Z   : { _dqdt[_fixed[i]](2)=0.0;                                                   break; }
			case FIX_XY  : { _dqdt[_fixed[i]](0)=0.0; _dqdt[_fixed[i]](1)=0.0;                          break; }
			case FIX_YZ  : { _dqdt[_fixed[i]](1)=0.0; _dqdt[_fixed[i]](2)=0.0;                          break; }
			case FIX_ZX  : { _dqdt[_fixed[i]](2)=0.0; _dqdt[_fixed[i]](0)=0.0;                          break; }
			case FIX_XYZ : { _dqdt[_fixed[i]](0)=0.0; _dqdt[_fixed[i]](1)=0.0; _dqdt[_fixed[i]](2)=0.0; break; }
		}
	}
}

inline void Grid2D::ApplyLoadings(double M)
{
	for (size_t i=0; i<_loaded.Size(); ++i)
		_fe[_loaded[i]] = M*_loads[i];
}

inline void Grid2D::Refine()
{
	// New geometry
	_nrow = 2*_nrow-1;
	_ncol = 2*_ncol-1;
	_L    = _L(0)/2.0, _L(1)/2.0, 0.0;

	// Do refine
	_refine();
}


/* private */

inline void Grid2D::_refine()
{
	// Allocate nodes
	_n.Resize(_nrow*_ncol);

	// Set nodes coordinates
	size_t k = 0;
	for (int i=0; i<_nrow; ++i)
	for (int j=0; j<_ncol; ++j)
	{
		_n[k] = _xmin+j*_L(0),
		        _ymin+i*_L(1),
				0.0;
		k++;
	}
	_xmax = _xmin+(_ncol-1)*_L(0);
	_ymax = _ymin+(_nrow-1)*_L(1);

	// Allocate and initialize (zero) nodes data
	_m   .Resize(_n.Size()); // nodes masses
	_q   .Resize(_n.Size()); // nodes momenta
	_v   .Resize(_n.Size()); // nodes velocities
	_fe  .Resize(_n.Size()); // nodes external forces
	_fi  .Resize(_n.Size()); // nodes internal forces
	_dqdt.Resize(_n.Size()); // nodes momenta rate
	ClearState();

	// Allocate cells connectivities
	_c.Resize((_nrow-1)*(_ncol-1));

	// Set connectivities
	for (int i=0; i<_nrow-1; ++i)
	{
		for (int j=0; j<_ncol-1; ++j)
		{
			_c[i*(_ncol-1)+j] = i   *_ncol+j  ,
			                    i   *_ncol+j+1,
			                   (i+1)*_ncol+j+1,
			                   (i+1)*_ncol+j  ;
		}
	}

	// Set boundary conditions
	_fixed  .Resize(0);
	_fixtype.Resize(0);
	FixType type;
	for (size_t i=0; i<_n.Size(); ++i)
	{
		if ((*_is_fixed) (_n[i], type))
		{
			long k = _fixed.Find(i);
			if (k<0) // not added yet
			{
				_fixed  .Push(i);
				_fixtype.Push(type);
			}
			else
			{
				_fixed  [k] = i;
				_fixtype[k] = type;
			}
		}
	}
}


#endif // MPM_GRID2D_H

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