mpoints2d.h

crack modeling with xfem

			}
			DyadDot (_p[p], _sg[p].G[k], _F[p],  newF); // deformation gradient
			k++;
		}

		// Update deformation gradient
		_F[p] = newF;

		// Update displacements
		_u[p] = _p[p] - _p_ini[p];

		// Update left-bottom node id
		double xlb = _g->N(_lbn[p])(0); // x-left-bottom
		double ylb = _g->N(_lbn[p])(1); // y-left-bottom
		if (_p[p](0)<xlb || _p[p](0)>xlb+3.0*_g->L(0) || _p[p](1)<ylb || _p[p](1)>ylb+3.0*_g->L(1))
		{
			std::cout << "Error: material point jumped over one cell" << std::endl;
			return false;
		}
		if (!_g->IsInside(_p[p]))
		{
			std::cout << "Error: material point moved outside the grid" << std::endl;
			return false;
		}
		if (_p[p](0)<xlb+    _g->L(0)) _lbn[p]--;           // moved to the left
		if (_p[p](0)>xlb+2.0*_g->L(0)) _lbn[p]++;           // moved to the right
		if (_p[p](1)<ylb+    _g->L(1)) _lbn[p]-=_g->nCol(); // moved to the bottom
		if (_p[p](1)>ylb+2.0*_g->L(1)) _lbn[p]+=_g->nCol(); // moved to the top
	}
	
	// Continue
	return true;
}

inline void MPoints2D::ResetDeform()
{
	for (int i=0; i<_p.Size(); ++i)
	{
		_F[i] = 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0;
		if (_mdl[i]!=NULL) _mdl[i]->Eps() = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0;
	}
}


/* private */

inline void MPoints2D::_initialize()
{
	// Check input
	if (!(_npcell==1 || _npcell==4)) throw new Fatal("MPoints2D::_initialize: Number of points per cell (NPcell=%d) must be 1 or 4.",_npcell);

	// Gauss points
	double gv     = (_npcell==1 ? 0.0 : 0.5); // Gauss point value (must be 0.5 because the definition of lp)
	double ksi[4] = {-1,1,1,-1};              // Natural coordinates
	double eta[4] = {-1,-1,1,1};              // Natural coordinates
	double kgp[4] = {-gv,gv,gv,-gv};          // ksi Gauss point
	double egp[4] = {-gv,-gv,gv,gv};          // eta Gauss point

	// Loop over cells
	_p    .Resize(0);
	_lbn  .Resize(0);
	_t    .Resize(0);
	_has_t.Resize(0);
	for (size_t c=0; c<_g->nCells(); ++c)
	{
		// For each cell c, generate _npcell mat points
		for (size_t j=0; j<_npcell; ++j)
		{
			// For each node j, sum the contribution of node c
			Vector3D p; p = 0.0;
			Array<Vector3D> nodes(4);
			for (int l=0; l<4; ++l)
			{
				nodes[l] = _g->N(_g->C(c)(l));
				double shape = (1.0+ksi[j]*kgp[l])*(1.0+eta[j]*egp[l])/4.0;
				p(0) += shape*nodes[l](0);
				p(1) += shape*nodes[l](1);
			}

			// Check if the mat point is inside the geometry
			if ((*_is_point_in_geom)(p))
			{
				// Coords
				_p  .Push (p);                        // coordinates
				_lbn.Push (_g->C(c)(0)-1-_g->nCol()); // left-bottom node in the left-bottom cell

				// Boundary conds
				_t    .Push(0.0);
				_has_t.Push((*_has_traction) (_p[_p.Size()-1], nodes, _g->L(), _npcell, _t[_p.Size()-1]));
			}
		}
	}

	// Initialize arrays
	_p_ini.Resize(_p.Size());
	_l    .Resize(_p.Size());
	_m    .Resize(_p.Size());
	_f    .Resize(_p.Size());
	_u    .Resize(_p.Size());
	_v    .Resize(_p.Size());
	_F    .Resize(_p.Size());
	_sg   .Resize(_p.Size());
	_mdl  .Resize(_p.Size());
	for (int i=0; i<_p.Size(); ++i)
	{
		_p_ini [i] = _p[i];
		if (_npcell==1) _l[i] = _g->L(0)/2.0, _g->L(1)/2.0, 0.0;
		else            _l[i] = _g->L(0)/4.0, _g->L(1)/4.0, 0.0;
		double Vp = 4.0*_l[i](0)*_l[i](1); // particle volume
		_m     [i] = Vp * (*_density) (_p[i]);
		_f     [i] = 0.0, 0.0, 0.0;
		_u     [i] = 0.0, 0.0, 0.0;
		(*_ini_velocity) (_p[i], _v[i]);
		_F     [i] = 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0;
		_mdl   [i] = (*_allocate_model) (_p[i]);
	}
}

inline void MPoints2D::_stress_update(double Dt, double & DsE)
{
	// Calculate nodes velocities (used only for stress update)
	for (size_t n=0; n<_g->nNodes(); ++n) { if (_g->m(n)>MINMASS) _g->v(n) = _g->q(n)/_g->m(n); }

	// Set fixed nodes (velocity)
	_g->FixNodesVel();

	// Stress/strain update
	DsE = 0.0;
	for (int p=0; p<_p.Size(); ++p)
	{
		if (_mdl[p]!=NULL)
		{
			// Variables
			STensor2 sp0; sp0 = _mdl[p]->Sig(); // initial point stress
			ATensor2 L;   L   = 0.0;            // mat point velocity gradient

			// Compute the velocity gradient
			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
				DyadUp (_g->v(n), _sg[p].G[k],  L); // L += v dyad G
				k++;
			}
			
			// Compute the strain increment == rate of deformation D*Dt
			STensor2 de; Sym (Dt, L,  de); // de = Dt*sym(L)

			// Update the stress state
			_mdl[p]->Update (de);

			// Strain energy
			double Vp = 4.0*_l[p](0)*_l[p](1); // particle volume
			DsE += 0.5*blitz::dot((sp0+_mdl[p]->Sig()),de)*Vp; // Strain Energy on points
		}
	}
}

inline double MPoints2D::_mpm_snp(int n, int p, int icomp) const
{
	double d = _p[p](icomp)-_g->N(n)(icomp);
	double L = _g->L(icomp);
	return (d<=-L  ? 0.0     :
		   (d>  L  ? 0.0     :
		   (d<=0.0 ? 1.0+d/L :
					 1.0-d/L )));
}

inline double MPoints2D::_mpm_gnp(int n, int p, int icomp) const
{
	double d = _p[p](icomp)-_g->N(n)(icomp);
	double L = _g->L(icomp);
	return (d<=-L   ?  0.0   :
		   (d>  L   ?  0.0   :
		   (d<= 0.0 ?  1.0/L :
					  -1.0/L )));
}

inline double MPoints2D::_snp(int n, int p, int icomp) const
{
	if (_mpm) return _mpm_snp(n,p,icomp);
	double d  = _p[p](icomp)-_g->N(n)(icomp);
	double L  = _g->L(icomp);
	double lp = _l[p](icomp);
	return (d<=-L-lp ? 0.0                        :
		   (d>  L+lp ? 0.0                        :
		   (d<=-L+lp ? pow(L+lp+d,2.0)/(4.0*L*lp) :
		   (d<=  -lp ? 1.0+d/L                    :
		   (d<=   lp ? 1.0-(d*d+lp*lp)/(2.0*L*lp) :
		   (d<= L-lp ? 1.0-d/L                    :
					   pow(L+lp-d,2.0)/(4.0*L*lp) ))))));
}

inline double MPoints2D::_gnp(int n, int p, int icomp) const
{
	if (_mpm) return _mpm_gnp(n,p,icomp);
	double d  = _p[p](icomp)-_g->N(n)(icomp);
	double L  = _g->L(icomp);
	double lp = _l[p](icomp);
	return (d<=-L-lp ? 0.0                  :
		   (d>  L+lp ? 0.0                  :
		   (d<=-L+lp ? (L+d+lp)/(2.0*L*lp)  :
		   (d<=  -lp ? 1.0/L                :
		   (d<=   lp ? -d/(L*lp)            :
		   (d<= L-lp ? -1.0/L               :
					   (-L+d-lp)/(2.0*L*lp) ))))));
}

inline void MPoints2D::_Snp(int n, int p, double & S) const
{
	S = _snp(n,p,0)*_snp(n,p,1);
}

inline void MPoints2D::_Gnp(int n, int p, Vector3D & G) const 
{
	G = _snp(n,p,1)*_gnp(n,p,0),
	    _snp(n,p,0)*_gnp(n,p,1),
	    0.0;
}


#endif // MPM_MPOINTS2D_H

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