mpm1d.cpp

crack modeling with xfem

			// Plot kinetic energy
			fl_color (FL_RED);
			for (int ti=1; ti<_mw->_ti+1; ++ti)
			{
				int x0=_x(_mw->_dt*(ti-1));  int y0=_y(_mw->_Ekp[ti-1]);
				int x1=_x(_mw->_dt*ti    );  int y1=_y(_mw->_Ekp[ti  ]);
				fl_line  (x0,y0, x1,y1);
			}

			//fl_color (FL_BLACK);
			//fl_line (_x(_Xmin),_y(_Ymin), _x(_Xmax),_y(_Ymax));
		}

	private:
		// Data
		double             _Xmin; // Real coordinates
		double             _Ymin; // Real coordinates
		double             _Xmax; // Real coordinates
		double             _Ymax; // Real coordinates
		double             _sfx;  // Scale factor: x=xmin+(X-Xmin)*sf and y=ymin+ymax-(Y-Ymin)*sf, where x and y are screen coordinates
		double             _sfy;  // Scale factor: x=xmin+(X-Xmin)*sf and y=ymin+ymax-(Y-Ymin)*sf, where x and y are screen coordinates
		MainWindow const * _mw;   // Access to the MainWindow

		// Private methods
		inline int _x (double   X) const { return static_cast<int>((x()+1)        +_sfx*(X-_Xmin)); } // X in screen coordinates
		inline int _y (double   Y) const { return static_cast<int>((y()+1)+(h()-2)-_sfy*(Y-_Ymin)); } // Y in screen coordinates
		inline int _l (double   L) const { return static_cast<int>(0.5*(_sfx+_sfy)*L)             ; } // Length in screen coordinates

	}; // class PlotEnergy }}}

	friend class DrawArea;
	friend class PlotVeloc;
	friend class PlotEnergy;

	// Constructor
	MainWindow(int width=1024, int height=824)
		: Fl_Double_Window (width,height,"Material Point Method (2D)"), // {{{
		  _prob    (1),
		  _tf      (0.0),
		  _dt      (0.01),
		  _nt      (0),
		  _ti      (0),
		  _b       (0),
		  _zero    (sqrt(DBL_EPSILON)),
		  _npcell  (2),
		  _mpm     (true),
		  _usf     (true),
		  _defgrad (false),
		  _outvp   (0),
		  _ntstep  (1),

		  _E  (1000),

		  _np      (0),
		  _p       (NULL),
		  _mp      (NULL),
		  _vp      (NULL),
		  _sp      (NULL),
		  _ep      (NULL),
		  _pnodes  (NULL),
		  _S       (NULL),
		  _G       (NULL),
		  _fixedn  (NULL),

		  _nn    (0),
		  _l     (0),
		  _n     (NULL),
		  _mn    (NULL),
		  _qn    (NULL),
		  _vn    (NULL),
		  _fen   (NULL),
		  _fin   (NULL),
		  _dqndt (NULL),

		  _W   (NULL),
		  _Esp (NULL),
		  _Ekp (NULL),
		  _Ovp (NULL),
		  _Osp (NULL)
	{

		// Add widgets to window
		begin();
			// Widgets
			_wnpcell  = new Fl_Input       ( 70, 2,  30, 25, "NPcell =");
			_wprob    = new Fl_Input       (180, 2,  30, 25, "Prob # =");
			_wsetprob = new Fl_Button      (220, 2, 100, 25, "Set &Prob");
			_wmpm     = new Fl_Check_Button(340, 2,  30, 25, "MPM");
			_wusf     = new Fl_Check_Button(400, 2,  30, 25, "USF");
			_wdefgrad = new Fl_Check_Button(455, 2,  30, 25, "e=dx/dX");
			_wntstep  = new Fl_Input       (620, 2,  60, 25, "NtStep =");
			_wstep    = new Fl_Button      (680, 2,  80, 25, "&Step");
			_wrun     = new Fl_Button      (760, 2,  80, 25, "&Run" );
			_wquit    = new Fl_Button      (840, 2,  80, 25, "&Quit");
			// Initialize widgets
			char buf[256];
			snprintf(buf, 256, "%d", _npcell); _wnpcell->value(buf);
			snprintf(buf, 256, "%d", _prob  ); _wprob  ->value(buf);
			snprintf(buf, 256, "%d", _ntstep); _wntstep->value(buf);
			if (_mpm)     _wmpm    ->set();
			if (_usf)     _wusf    ->set();
			if (_defgrad) _wdefgrad->set();
			// Callbacks
			_wsetprob->callback(_setprob_cb, this);
			_wstep   ->callback(_step_cb   , this);
			_wrun    ->callback(_run_cb    , this);
			_wquit   ->callback(_quit_cb   , this);
			// Draw area
			_wda = new DrawArea  (/*xmin*/5    , /*ymin*/30        , /*width*/w()/2-10, /*height*/h()-35  );
			_wpv = new PlotVeloc (/*xmin*/w()/2, /*ymin*/30        , /*width*/w()/2-5,  /*height*/h()/3-5 );
			_wps = new PlotStress(/*xmin*/w()/2, /*ymin*/  h()/3+30, /*width*/w()/2-5,  /*height*/h()/3-10);
			_wpe = new PlotEnergy(/*xmin*/w()/2, /*ymin*/2*h()/3+25, /*width*/w()/2-5,  /*height*/h()/3-29);
		end();

		// Set focus
		_wrun->take_focus();

		// Set resizable and show window
		//resizable(this);
		resizable(_wda);
		show();

	} // }}}

	// Destructor
	~MainWindow()
	{ // {{{

		if (_p      !=NULL) delete [] _p;
		if (_mp     !=NULL) delete [] _mp;
		if (_vp     !=NULL) delete [] _vp;
		if (_sp     !=NULL) delete [] _sp;
		if (_ep     !=NULL) delete [] _ep;
		if (_pnodes !=NULL) delete [] _pnodes;
		if (_S      !=NULL) delete [] _S;
		if (_G      !=NULL) delete [] _G;
		if (_fixedn !=NULL) delete [] _fixedn;

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

		if (_W  !=NULL) delete [] _W; 
		if (_Esp!=NULL) delete [] _Esp; 
		if (_Ekp!=NULL) delete [] _Ekp; 
		if (_Ovp!=NULL) delete [] _Ovp; 
		if (_Osp!=NULL) delete [] _Osp; 

	} // }}}

private:
	// Widgets
	Fl_Input        * _wnpcell;
	Fl_Input        * _wprob;
	Fl_Check_Button * _wmpm;
	Fl_Check_Button * _wusf;
	Fl_Check_Button * _wdefgrad;
	Fl_Input        * _wntstep;
	Fl_Button       * _wsetprob;
	Fl_Button       * _wstep;
	Fl_Button       * _wrun;
	Fl_Button       * _wquit;
	DrawArea        * _wda; 
	PlotVeloc       * _wpv; 
	PlotStress      * _wps; 
	PlotEnergy      * _wpe; 

	// Variables
	int      _prob;    // problem number
	double   _tf;      // final time
	double   _dt;      // time-step
	int      _nt;      // number of time-steps
	int      _ti;      // current time increment (current time: t=ti*dt)
	double   _b;       // current b
	double   _zero;    // machine epsilon (smaller positive)
	int      _npcell;  // initial number of particles per cell
	bool     _mpm;     // use original MPM shape functions ?
	bool     _usf;     // update stress first ?
	bool     _defgrad; // use deformation gradient instead of Green small strain ?
	int      _outvp;   // index of a point for output
	int      _ntstep;  // number of time-steps for method _step
	ptB      _ptb;     // pointer to a function that calculates b (gravity) for the current time t

	// Parameters
	double _E;    // Young
	double _len0; // initial bar length

	// Material points
	int       _np;      // number of mat points
	double    _lp;      // mat points length == _l/2
	double  * _p;       // (size=np) mat points (x coords)
	double  * _mp;      // (size=np) mat points masses
	double  * _vp;      // (size=np) mat points velocities (x comps)
	double  * _sp;      // (size=np) stresses (tensors) on each point
	double  * _ep;      // (size=np) strains (tensors) on each point
	MapsAry * _pnodes;  // (size=np) nodes that each point contributes to
	ShapAry * _S;       // (size=np) interpolation functions
	ShapAry * _G;       // (size=np) derivative of interp. functions

	// Cell/grid nodes
	double    _xmin;  // min x value
	double    _xmax;  // max x value
	int       _nn;    // number of nodes
	double    _l;     // cell-length
	double  * _n;     // (size=nn) nodes (x coords)
	double  * _mn;    // (size=nn) nodes masses
	double  * _qn;    // (size=nn) nodes momenta (x comps)
	double  * _vn;    // (size=nn) nodes velocities (x comps)
	double  * _fen;   // (size=nn) nodes external forces (x comps)
	double  * _fin;   // (size=nn) nodes internal forces (x comps)
	double  * _dqndt; // (size=nn) rate of nodes momenta

	// Output arrays
	double * _W;   // (size=nt) output weights
	double * _Esp; // (size=nt) output Energy-strain-point
	double * _Ekp; // (size=nt) output Energy-kinetic-point
	double * _Ovp; // (size=nt*np) output points velocities
	double * _Osp; // (size=nt*np) output points stresses

	// Boundary conditions
	int   _nfixedn; // number of fixed nodes
	int * _fixedn;  // fixed nodes

	// Callbacks (must be static)
	static void _setprob_cb (Fl_Widget * o, void * v) { ((MainWindow*)v)->_setprob  (); }
	static void _step_cb    (Fl_Widget * o, void * v) { ((MainWindow*)v)->_step     (); }
	static void _run_cb     (Fl_Widget * o, void * v) { ((MainWindow*)v)->_run      (); }
	static void _quit_cb    (Fl_Widget * o, void * v) { ((MainWindow*)v)->_quit     (); }

	// SetProb callback
	inline void _setprob()
	{ // {{{

		// Deallocate previous material points
		if (_p      !=NULL) delete [] _p;       _p      =NULL;
		if (_mp     !=NULL) delete [] _mp;      _mp     =NULL;
		if (_vp     !=NULL) delete [] _vp;      _vp     =NULL;
		if (_sp     !=NULL) delete [] _sp;      _sp     =NULL;
		if (_ep     !=NULL) delete [] _ep;      _ep     =NULL;
		if (_pnodes !=NULL) delete [] _pnodes;  _pnodes =NULL;
		if (_S      !=NULL) delete [] _S;       _S      =NULL;
		if (_G      !=NULL) delete [] _G;       _G      =NULL;
		if (_fixedn !=NULL) delete [] _fixedn;  _fixedn =NULL;
		if (_W      !=NULL) delete [] _W;       _W      =NULL;
		if (_Esp    !=NULL) delete [] _Esp;     _Esp    =NULL;
		if (_Ekp    !=NULL) delete [] _Ekp;     _Ekp    =NULL;
		if (_Ovp    !=NULL) delete [] _Ovp;     _Ovp    =NULL;
		if (_Osp    !=NULL) delete [] _Osp;     _Osp    =NULL;

		// Input
		_npcell = atoi(_wnpcell->value());
		_prob   = atoi(_wprob  ->value());

		// Set problem
		if (_prob==1 || _prob==2 || _prob==3)
		{
			// Constants
			int nc;
			if (_prob==1)
			{
				nc      = 1;
				_dt     = 0.1/(2.0*PI);
				_tf     = 10.0;
				_E      = 4.0*PI*PI;
				_outvp  = 0;
				_ptb    = &Prob1B;
				_wpv->SetPtrVelCM (&Prob1CorrectVelCM );
				_wps->SetPtrStress(&Prob1CorrectStress);
			}
			else if (_prob==2)
			{
				nc      = 25;
				_dt     = 0.01;
				_tf     = 100.0;
				_E      = 100.0;
				_outvp  = 24;
				_ptb    = &Prob2B;
				_wpv->SetPtrVelCM (&Prob2CorrectVelCM );
				_wps->SetPtrStress(&Prob2CorrectStress);
			}
			else if (_prob==3)
			{
				nc      = 50;
				_dt     = 1.0e-4;
				_tf     = 2.0;
				_E      = 10.0e+6;
				_outvp  = 49;
				_ptb    = &Prob3B;
				_wpv->SetPtrVelCM (&Prob3CorrectVelCM );
				_wps->SetPtrStress(&Prob3CorrectStress);
			}
			_len0 = nc*1.0;
			_np   = nc*_npcell;
			_lp   = 1.0/(2.0*_npcell); // 1.0==L

			// Allocate material points data
			_p      = new double  [_np];
			_mp     = new double  [_np];
			_vp     = new double  [_np];
			_sp     = new double  [_np];
			_ep     = new double  [_np];
			_pnodes = new MapsAry [_np];
			_S      = new ShapAry [_np];
			_G      = new ShapAry [_np];

			// Set points position
			_p[0] = _lp;
			for (int p=1; p<_np; ++p) _p[p]=_p[p-1]+2*_lp;

			// Set points mass, velocities, stress, and strains
			for (int p=0; p<_np; ++p)
			{
				_mp[p] = 2*_lp;
				     if (_prob==1) _vp[p] = 0.1;
				else if (_prob==2) _vp[p] = 0.1*sin(PI*_p[p]/50.0);
				else if (_prob==3) _vp[p] = 0.0;
				_sp[p] = 0.0;
				_ep[p] = 0.0;
			}

			// Boundary conditions
			_nfixedn   = 1;
			_fixedn    = new int [_nfixedn];
			_fixedn[0] = 0;

			// Regenerate the grid
			_re_gen_grid(0.0,nc,1.0);

			// Set points/nodes maps
			int p = 0;
			for (int n=0; n<_nn-1; ++n)
			{
				for (int i=0; i<_npcell; ++i)
				{
					_pnodes[p] = n-1,n,n+1,n+2;
					p++;
				}
			}
		}
		else return;

		// Allocate output arrays
		_nt  = static_cast<int>(_tf/_dt+1.0);
		_W   = new double [_nt+1      ]; // Bar weight
		_Esp = new double [_nt+1      ]; // Energy-strain-points
		_Ekp = new double [_nt+1      ]; // Energy-kinetic-points
		_Ovp = new double [(_nt+1)*_np]; // Output mat points velocity
		_Osp = new double [(_nt+1)*_np]; // Output mat points stress
		for (int ti=0; ti<_nt+1; ++ti)
		{
			_W  [ti] = 0.0;
			_Esp[ti] = 0.0;
			_Ekp[ti] = 0.0;
			for (int p=0; p<_np; ++p)