mpm1d.cpp

crack modeling with xfem

			{
				_Ovp[ti*_np+p] = 0.0;
				_Osp[ti*_np+p] = 0.0;
			}
		}

		// Initialize control variables
		_ti = 0;   // current time increment
		_b  = 0.0; // current gravity

		// Initialize output arrays
		_W  [0] = 0.0;
		_Esp[0] = 0.0;
		_Ekp[0] = 0.0;
		_Ovp[0] = _vp[_outvp];
		_Osp[0] = 0.0;
		for (int p=0; p<_np; ++p) _Ekp[0] += 0.5*_mp[p]*_vp[p]*_vp[p];

		// Set Draw area
		_wda->SetMainWindow   (this);
		_wda->SetRadiusNodes  (_l/8);
		_wda->SetRadiusPoints (_l/9);
		_wda->redraw          ();

		// Set Plot velocity
		_wpv->SetMainWindow (this);
		_wpv->redraw        ();

		// Set Plot stress
		_wps->SetMainWindow (this);
		_wps->redraw        ();

		// Set Plot energy
		_wpe->SetMainWindow (this);
		_wpe->redraw        ();

		// Wait
		Fl::wait(0);

	} // }}}

	// Run callback
	inline void _run()
	{ // {{{

		if (_p!=NULL && _n!=NULL)
		{
			_ti     = 0;
			_ntstep = _nt;
			char buf[256];
			snprintf(buf, 256, "%d", _ntstep); _wntstep->value(buf);
			_step();
		}

	} // }}}

	// Step callback
	inline double _step()
	{ // {{{

		if (_p!=NULL && _n!=NULL)
		{
			// Get input
			_ntstep  = atoi(_wntstep->value());
			_mpm     = (_wmpm    ->value()==1 ? true : false);
			_usf     = (_wusf    ->value()==1 ? true : false);
			_defgrad = (_wdefgrad->value()==1 ? true : false);

			// Initialize deformation gradient
			for (int p=0; p<_np; ++p) _ep[p] = 1.0;

			// Run
			for (int tstep=0; tstep<_ntstep && _ti<_nt; ++tstep)
			{
				// 1) Time increment and body forces
				_ti++;
				_b = (*_ptb)(_ti*_dt);

				// 2) Calculate interpolation functions
				for (int p=0; p<_np; ++p)
				{
					for (int i=0; i<4; ++i)
					{
						int n = _pnodes[p](i); // my node
						if (n>=0 && n<_nn)
						{
							_S[p](i) = _Snp (n,p); // interpolation function
							_G[p](i) = _Gnp (n,p); // deriv. interp. function
						}
					}
				}

				// 3) Clear grid
				for (int n=0; n<_nn; ++n)
				{
					_mn   [n] = 0.0; // nodes masses
					_qn   [n] = 0.0; // nodes momenta
					_vn   [n] = 0.0; // nodes velocities
					_fen  [n] = 0.0; // nodes external forces
					_fin  [n] = 0.0; // nodes internal forces
					_dqndt[n] = 0.0; // nodes momenta rate
				}

				// 4) Initialize grid state from points (mass and momentum)
				for (int p=0; p<_np; ++p)
				{
					for (int i=0; i<4; ++i)
					{
						int n = _pnodes[p](i); // my node
						if (n>=0 && n<_nn)
						{
							_mn[n] = _mn[n] + _S[p](i)*_mp[p];        // node mass
							_qn[n] = _qn[n] + _S[p](i)*_mp[p]*_vp[p]; // node momentum
						}
					}
				}

				// 5) Update strains and stress (USF - update stress first)
				double dEsp = 0.0; // Energy-strain-points increment
				if (_usf) _stress_update(dEsp);

				// 6) Calculate internal and external forces
				for (int p=0; p<_np; ++p)
				{
					for (int i=0; i<4; ++i)
					{
						int n = _pnodes[p](i); // my node
						if (n>=0 && n<_nn)
						{
							_fin[n] += _G[p](i)*_sp[p]*2*_lp; // node internal force (Vp=2lp)
							_fen[n] += _S[p](i)*_mp[p]*_b;    // node external force
						}
					}
				}

				// 7) Calculate nodes momenta rate
				for (int n=0; n<_nn; ++n) _dqndt[n] = _fen[n] - _fin[n];

				// 8) Set fixed nodes (again)
				for (int i=0; i<_nfixedn; ++i) _dqndt[_fixedn[i]] = 0.0;

				// 9) Update nodes momenta
				for (int n=0; n<_nn; ++n) _qn[n] += _dqndt[n]*_dt;

				// 10) Update strains and stress (USL - update stress last)
				if (!_usf) _stress_update(dEsp);

				// 11) Update material points (position and velocity)
				for (int p=0; p<_np; ++p)
				{
					// Update position and velocity
					for (int i=0; i<4; ++i)
					{
						int n = _pnodes[p](i); // my node
						if (n>=0 && n<_nn)
						{
							_p [p] += _dt*_S[p](i)*_qn[n]/_mn[n];    // point position
							_vp[p] += _dt*_S[p](i)*_dqndt[n]/_mn[n]; // point velocity
						}
					}
					// Update list of nodes to contribute to
					     if (_p[p]<_n[_pnodes[p](1)]) _pnodes[p] -= 1; // moved to the left
					else if (_p[p]>_n[_pnodes[p](2)]) _pnodes[p] += 1; // moved to the right
					if (_pnodes[p](1)<0)
					{
						cout << "mat point moved outside grid (to the left)" << endl;
						return -1;
					}
					if (_pnodes[p](2)>=_nn)
					{
						cout << "mat point moved outside grid (to the right)" << endl;
						return -1;
					}
				}

				// 12) Output arrays
				_W  [_ti] = _W  [_ti-1] + fabs(_b)*_len0;
				_Esp[_ti] = _Esp[_ti-1] + dEsp;
				_Ekp[_ti] = 0.0;
				for (int p=0; p<_np; ++p)
				{
					_Ekp[_ti] += 0.5*_mp[p]*_vp[p]*_vp[p];
					_Ovp[_ti*_np+p] = _vp[p];
					_Osp[_ti*_np+p] = _sp[p];
				}

				// 13) Plot stresses
				_wps->SetOutTimeInc(_ti);
				_wps->redraw       ();
				Fl::wait           (0);
			}

			// Plot results
			_wda->redraw ();
			_wpv->redraw ();
			_wpe->redraw ();
			Fl::wait     (0);
		}
		else return -1.0;

	} // }}}

	// Stress-update
	void _stress_update(double & dEsp) // returns dEsp (increment of Energy-strain-points)
	{ // {{{

		// Calculate nodes velocities (used only for stress update)
		for (int n=0; n<_nn; ++n) _vn[n] = _qn[n]/_mn[n];

		// Set fixed nodes
		for (int i=0; i<_nfixedn; ++i) _vn[_fixedn[i]] = 0.0;

		// Stress/strain update
		for (int p=0; p<_np; ++p)
		{
			double sp0 = _sp[p]; // initial point stress
			double de  = 0.0;    // mat point strain increment
			double ff  = 1.0;    // deformation gradient F multiplier
			for (int i=0; i<4; ++i)
			{
				int n = _pnodes[p](i); // my node
				if (n>=0 && n<_nn)
				{
					ff += _G[p](i)*_vn[n]*_dt;                       // def. gradient mult.
					de += 0.5*(_G[p](i)*_vn[n]+_vn[n]*_G[p](i))*_dt; // strain increment
				}
			}
			if (_defgrad)
			{
				_ep[p] = _ep[p]*ff;       // update def. grad.
				_sp[p] = _E*(_ep[p]-1.0); // update stress
			}
			else
			{
				_ep[p] += de;    // update strain
				_sp[p] += _E*de; // update stress
			}
			dEsp += 0.5*(sp0+_sp[p])*de*2*_lp; // Energy-strain-points (Vp=2lp)
		}

	} // }}}

	// Quit callback
	inline void _quit() { hide(); }

	// Shape functions (MPM)
	inline double _mpm_Snp(int n, int p)
	{ // {{{
		return (_p[p]-_n[n]<=-_l  ? 0.0                  :
		       (_p[p]-_n[n]>  _l  ? 0.0                  :
		       (_p[p]-_n[n]<= 0.0 ? 1.0+(_p[p]-_n[n])/_l :
		                            1.0-(_p[p]-_n[n])/_l )));
	} // }}}

	// Derivative of shape functions (MPM)
	inline double _mpm_Gnp(int n, int p)
	{ // {{{
		return (_p[p]-_n[n]<=-_l  ?  0.0    :
		       (_p[p]-_n[n]>  _l  ?  0.0    :
		       (_p[p]-_n[n]<= 0.0 ?  1.0/_l :
		                            -1.0/_l )));
	} // }}}

	// Shape functions (GMPM)
	inline double _Snp(int n, int p)
	{ // {{{
		if (_mpm) return _mpm_Snp(n,p);
		return (_p[p]-_n[n]<=-_l-_lp ? 0.0                                             :
		       (_p[p]-_n[n]>  _l+_lp ? 0.0                                             :
		       (_p[p]-_n[n]<=-_l+_lp ? (pow(_l+_lp+(_p[p]-_n[n]),2.0))/(4.0*_l*_lp)    :
		       (_p[p]-_n[n]<=   -_lp ? 1.0+(_p[p]-_n[n])/_l                            :
		       (_p[p]-_n[n]<=    _lp ? 1.0-(pow(_p[p]-_n[n],2.0)+_lp*_lp)/(2.0*_l*_lp) :
		       (_p[p]-_n[n]<= _l-_lp ? 1.0-(_p[p]-_n[n])/_l                            :
		                               (pow(_l+_lp-(_p[p]-_n[n]),2.0))/(4.0*_l*_lp)    ))))));
	} // }}}

	// Derivative of shape functions (GMPM)
	inline double _Gnp(int n, int p)
	{ // {{{
		if (_mpm) return _mpm_Gnp(n,p);
		return (_p[p]-_n[n]<=-_l-_lp ? 0.0                                :
		       (_p[p]-_n[n]>  _l+_lp ? 0.0                                :
		       (_p[p]-_n[n]<=-_l+_lp ? (_l+_p[p]-_n[n]+_lp)/(2.0*_l*_lp)  :
		       (_p[p]-_n[n]<=   -_lp ? 1.0/_l                             :
		       (_p[p]-_n[n]<=    _lp ? -(_p[p]-_n[n])/(_l*_lp)            :
		       (_p[p]-_n[n]<= _l-_lp ? -1.0/_l                            :
		                               (-_l+_p[p]-_n[n]-_lp)/(2.0*_l*_lp) ))))));
	} // }}}

	// Regenerate grid
	void _re_gen_grid(double Xmin, int nCells, double L)
	{ // {{{

		// Deallocate previous grid
		if (_n    !=NULL) delete [] _n;      _n     = NULL;
		if (_mn   !=NULL) delete [] _mn;     _mn    = NULL;
		if (_qn   !=NULL) delete [] _qn;     _qn    = NULL;
		if (_vn   !=NULL) delete [] _vn;     _vn    = NULL;
		if (_fen  !=NULL) delete [] _fen;    _fen   = NULL;
		if (_fin  !=NULL) delete [] _fin;    _fin   = NULL;
		if (_dqndt!=NULL) delete [] _dqndt;  _dqndt = NULL;

		// Cell length
		_l = L;

		// Set limits
		_xmin = Xmin;
		_xmax = _xmin+nCells*_l;

		// Allocate nodes
		_nn = nCells+1;
		_n  = new double [_nn];

		// Set nodes coordinates
		for (int n=0; n<_nn; ++n) _n[n] = n*_l;

		// Allocate and initialize (zero) nodes data
		_mn    = new double [_nn];
		_qn    = new double [_nn];
		_vn    = new double [_nn];
		_fen   = new double [_nn];
		_fin   = new double [_nn];
		_dqndt = new double [_nn];
		for (int i=0; i<_nn; ++i)
		{
			_mn   [i] = 0.0;
			_qn   [i] = 0.0;
			_vn   [i] = 0.0;
			_fen  [i] = 0.0;
			_fin  [i] = 0.0;
			_dqndt[i] = 0.0;
		}

	} // }}}

}; // class MainWindow

/////////////////////////////////////////////////////////////////////////////////////////////// MainWindow /////

int main(int argc, char **argv) try
{
	MainWindow win;
	Fl::run();
	return 0;
}
catch (char const * m) // {{{
{
	std::cout << "Fatal: " << m << std::endl;
	exit (1);
}
catch (...)
{
	std::cout << "Some exception (...) ocurred\n";
} //}}} 

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