mpm2d.cpp

crack modeling with xfem

		snprintf(buf, 256, "%g", _prob.Dt    ); _w_dt     ->value(buf);
		snprintf(buf, 256, "%g", _prob.Dtout ); _w_dtout  ->value(buf);
		if (_prob.MPM)                          _w_mpm    ->set();
		if (_prob.USF)                          _w_usf    ->set();
		if (_prob.SmallSt)                      _w_smallst->set();

		// Callbacks
		_w_setprob->callback (_setprob_cb, this);
		_w_refine ->callback (_refine_cb , this);
		_w_step   ->callback (_step_cb   , this);
		_w_run    ->callback (_run_cb    , this);

		// Set focus
		_w_run->take_focus();
	}
	
	// Methods
	void SetFrames (Fl_Group * TL, Fl_Group * TR, Fl_Group * BL, Fl_Group * BR)
	{
		_tl = static_cast<TLframe*>(TL);
		_tr = static_cast<TRframe*>(TR);
		_bl = static_cast<BLframe*>(BL);
		_br = static_cast<BRframe*>(BR);
	}

	// Destructor
	~Control()
	{
		if (_prob.G!=NULL) delete _prob.G;
		if (_prob.P!=NULL) delete _prob.P;
	}

protected:
	// No resizable
	void resize(int X, int Y, int W, int H) {}

private:
	// Data
	Problem2D       & _prob;
	Fl_Input        * _w_prob;
	Fl_Input        * _w_npcell;
	Fl_Check_Button * _w_mpm;
	Fl_Check_Button * _w_usf;
	Fl_Check_Button * _w_smallst;
	Fl_Input        * _w_tf;
	Fl_Input        * _w_dt;
	Fl_Input        * _w_dtout;
	Fl_Button       * _w_setprob;
	Fl_Button       * _w_refine;
	Fl_Button       * _w_step;
	Fl_Button       * _w_run;
	bool              _running;
	bool              _do_stop;

	// Frames
	TLframe * _tl;
	TRframe * _tr;
	BLframe * _bl;
	BRframe * _br;

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

	// Set problem
	void _setprob()
	{
		// Input
		_prob.Id     = atoi(_w_prob  ->value());
		_prob.NPcell = atoi(_w_npcell->value());
		if (!(_prob.NPcell==1 || _prob.NPcell==4))
		{
			InfoBox ("The number of points per cell (NPcell) must be 1 or 4");
			return;
		}

		// Set problem
		SetProblem (_prob);
		
		// Set GUI back
		char buf[256];
		snprintf(buf, 256, "%g", _prob.tf    ); _w_tf    ->value(buf);
		snprintf(buf, 256, "%g", _prob.Dt    ); _w_dt    ->value(buf);
		snprintf(buf, 256, "%g", _prob.Dtout ); _w_dtout ->value(buf);

		// Initial output
		_do_output (0.0);

		// Set draw area
		_bl->ReSetGridPointsAndOutp ();

		// Set plots
		_tl->SetCurves ();
		_tr->SetCurves ();
		_br->SetCurves ();

		// Redraw
		_do_redraw ();
	}

	void _refine()
	{
		if (_prob.G!=NULL)
		{
			_prob.t = 0.0;
			_prob.G->Refine   ();
			_prob.P->ReInit   ();
			ReSetOutputArrays (_prob);
			_tl->SetCurves    ();
			_tr->SetCurves    ();
			_br->SetCurves    ();
			_do_redraw        ();
		}
	}

	// Do one step
	bool _step()
	{
		if (_prob.G!=NULL && _prob.P!=NULL)
		{
			// Input
			_prob.tf      = atof(_w_tf   ->value());
			_prob.Dt      = atof(_w_dt   ->value());
			_prob.Dtout   = atof(_w_dtout->value());
			_prob.MPM     = (_w_mpm    ->value()==1 ? true : false);
			_prob.USF     = (_w_usf    ->value()==1 ? true : false);
			_prob.SmallSt = (_w_smallst->value()==1 ? true : false);
			if (fabs(_prob.t-_prob.tf)<1.0e-4 || _prob.t>_prob.tf) return false;

			// First output
			if (_prob.t==0.0) _do_output (0);

		   	// Next t for output
			double tout = _prob.t + _prob.Dtout;

			// Run
			double dse = 0.0;
			while (_prob.t<tout && _prob.t<_prob.tf)
			{
				// Time increment
				double dt = (_prob.t+_prob.Dt>tout ? tout-_prob.t : _prob.Dt);
				_prob.t += dt;

				// Body forces
				(*_prob.pB) (_prob.t, _prob.B);
				
				// Multiplier for external forces
				(*_prob.pM) (_prob.t, _prob.M);

				// Update
				double dset;
				if (!_prob.P->TimeUpdate (dt, _prob.B, _prob.M, dset)) return false;
				dse += dset;
			}

			// Do output
			_do_output (dse);

			// Redraw
			_do_redraw ();

			// Next step
			return true;
		}
		else return false;
	}

	// Run simulation
	void _run ()
	{
		if (_running) _do_stop = true;
		else
		{
			_running = true;
			_w_run->label("&Stop");
			while (!_do_stop && _step());
			_running = false;
			_do_stop = false;
			_w_run->label("&Run");
		}
	}

	// Do output (fill arrays)
	void _do_output (double DsE)
	{
		// Check
		if (_prob.P->nPoints()<1) return;

		// Time
		_prob.Out_t.Push (_prob.t);

		// Strain energy
		size_t n = _prob.Out_sE_t.Size(); // current number of elements
		double se = 0.0;                  // strain energy
		if (n>0) se = _prob.Out_sE_t[n-1]+DsE;
		_prob.Out_sE_t.Push (se);

		// Kinetic energy
		double ke = 0.0; // kinetic energy
		for (int i=0; i<_prob.P->nPoints(); ++i)
			ke += 0.5*_prob.P->m(i)*blitz::dot(_prob.P->v(i),_prob.P->v(i));
		_prob.Out_kE_t.Push (ke);

		// Total energy
		_prob.Out_tE_t.Push (se+ke);

		// Mat point state values
		for (size_t p=0; p<_prob.P->nPoints(); ++p)
		{
			// Mat point external forces
			_prob.Out_p_fx_t[p].Push (_prob.P->f(p)(0));
			_prob.Out_p_fy_t[p].Push (_prob.P->f(p)(1));

			// Mat point displacement
			_prob.Out_p_ux_t[p].Push (_prob.P->u(p)(0));
			_prob.Out_p_uy_t[p].Push (_prob.P->u(p)(1));

			// Mat point velocity
			_prob.Out_p_vx_t[p].Push (_prob.P->v(p)(0));
			_prob.Out_p_vy_t[p].Push (_prob.P->v(p)(1));
			
			// Mat point correct velocity
			if (_prob.pVeloc!=NULL)
			{
				Vector3D cvp; // correct velocity
				(*_prob.pVeloc) (_prob.t, _prob.P->P(p), cvp);
				_prob.Out_p_cvx_t[p].Push (cvp(0));
			}
			
			// Mat point stress
			_prob.Out_p_sxx_t[p].Push (_prob.P->s(p)(0));
			_prob.Out_p_syy_t[p].Push (_prob.P->s(p)(1));

			// Mat point strain
			_prob.Out_p_exx_t[p].Push (_prob.P->e(p)(0));
			_prob.Out_p_eyy_t[p].Push (_prob.P->e(p)(1));
		}
	}

	// Do redraw
	void _do_redraw ()
	{
		_tl->redraw();
		_tr->redraw();
		_bl->redraw();
		_br->redraw();
		Fl::wait(0);
	}

}; // class Control

////////////////////////////////////////////////////////////////////////////////////////////// main /////

int main(int argc, char **argv) try
{
	// Set the backgroud color of all widgets
	Fl::set_color (FL_BACKGROUND_COLOR, 0xefebe700);

	// Initialize problem
	Problem2D prob;
	prob.Id      = 5;     // default problem
	prob.NPcell  = 4;     // number of points in each cell
	prob.G       = NULL;  // grid
	prob.P       = NULL;  // mat points
	prob.pB      = NULL;  // body mass
	prob.pVeloc  = NULL;  // correct velocity
	prob.pStress = NULL;  // correct stress
	prob.MPM     = true;  // use original MPM
	prob.USF     = true;  // update stress first ?
	prob.SmallSt = false; // small strains ?
	prob.t       = 0.0;   // current time
	prob.tf      = 1.0;   // final time
	prob.Dt      = 0.01;  // time increment
	prob.Dtout   = 0.1;   // time increment for output
	prob.Outp    = 0;     // point/particle number to output data
	prob.B       = 0.0;   // body mass

	// Allocate window, controll and layout
	Fl_Double_Window win (1024, 768);
	Control          ctr (0, 0, win.w(), 28, prob);
	FourTiled        lay (0,28, win.w(), win.h()-28, &prob, &AllocTL, &AllocTR, &AllocBL, &AllocBR, 0.67, 2./5.);

	// Set frames
	ctr.SetFrames (lay.TL(), lay.TR(), lay.BL(), lay.BR());
	static_cast<BLframe*>(lay.BL())->SetFrames (lay.TL(), lay.TR(), lay.BR());

	// Set and show window
	win.resizable (&lay);
	win.end       ();
	win.show      ();

	// Run
	Fl::run();

	// End
	return 0;
}
catch (Fatal * f)
{
	f->Cout();
	delete f;
	exit (1);
}
catch (char const * m)
{
	std::cout << "Fatal: " << m << std::endl;
	exit (1);
}
catch (...)
{
	std::cout << "Some exception (...) ocurred\n";
}

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