mpm2d.cpp

crack modeling with xfem

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

// STL
#include <iostream>
#include <cmath>
#include <cfloat> // for DBL_EPSILON

// FLTK
#include <FL/Fl.H>
#include <FL/Fl_Group.H>
#include <FL/Fl_Box.H>
#include <FL/Fl_Button.H>
#include <FL/Fl_Round_Button.H>
#include <FL/Fl_Check_Button.H>
#include <FL/Fl_Input.H>
#include <FL/Fl_Double_Window.H>

// Local
#include "defs.h"
#include "array.h"
#include "fatal.h"
#include "plotxy.h"
#include "drawarea2d.h"
#include "problems.h"
#include "infobox.h"
#include "tiled.h"

using std::cout;
using std::endl;

/////////////////////////////////////////////////////////////////////////////////////////// TLframe /////

class TLframe : public Fl_Group
{
public:
	TLframe (int X, int Y, int W, int H, Problem2D const * Prob) : Fl_Group (X,Y,W,H,NULL), _pb(Prob)
	{
		int hp = (H-1)/2;
		_plt_f_u = new PlotXY (X+1         , Y    , (W-1)/2 , hp , "Force x Displacement"    , "Displacement" , "Force");
		_plt_s_e = new PlotXY (X+1+(W-1)/2 , Y    , (W-1)/2 , hp , "Stress x Strain"         , "Strain"       , "Stress");
		_plt_v_u = new PlotXY (X+1         , Y+hp , (W-1)/2 , hp , "Velocity x Displacement" , "Displacement" , "Velocity");
		_plt_f_t = new PlotXY (X+1+(W-1)/2 , Y+hp , (W-1)/2 , hp , "Force x Time"            , "Time"         , "Force");
		end();

		// Set Plot force x displacement
		_plt_f_u->AddCurve ("fx/ux");
		_plt_f_u->AddCurve ("fy/uy");
		_plt_f_u->SetCurve (0).Typ = CT_LINES;
		_plt_f_u->SetCurve (0).Clr = FL_BLUE;
		_plt_f_u->SetCurve (1).Typ = CT_LINES;
		_plt_f_u->SetCurve (1).Clr = FL_GREEN;

		// Set Plot stress x strain
		_plt_s_e->AddCurve ("exx/sxx");
		_plt_s_e->AddCurve ("eyy/syy");
		_plt_s_e->SetCurve (0).Typ = CT_LINES;
		_plt_s_e->SetCurve (0).Clr = FL_BLUE;
		_plt_s_e->SetCurve (1).Typ = CT_LINES;
		_plt_s_e->SetCurve (1).Clr = FL_GREEN;

		// Set Plot velocity x displacement
		_plt_v_u->AddCurve ("vx/ux");
		_plt_v_u->AddCurve ("vy/uy");
		_plt_v_u->SetCurve (0).Typ = CT_LINES;
		_plt_v_u->SetCurve (0).Clr = FL_BLUE;
		_plt_v_u->SetCurve (1).Typ = CT_LINES;
		_plt_v_u->SetCurve (1).Clr = FL_GREEN;

		// Set Plot displacement
		_plt_f_t->AddCurve ("fx");
		_plt_f_t->AddCurve ("fy");
		_plt_f_t->SetCurve (0).Typ = CT_LINES;
		_plt_f_t->SetCurve (0).Clr = FL_BLUE;
		_plt_f_t->SetCurve (1).Typ = CT_LINES;
		_plt_f_t->SetCurve (1).Clr = FL_GREEN;
	}
	void SetCurves ()
	{
		_plt_f_u->SetXY (0, &_pb->Out_p_ux_t [_pb->Outp], &_pb->Out_p_fx_t [_pb->Outp]);
		_plt_f_u->SetXY (1, &_pb->Out_p_uy_t [_pb->Outp], &_pb->Out_p_fy_t [_pb->Outp]);
		_plt_s_e->SetXY (0, &_pb->Out_p_exx_t[_pb->Outp], &_pb->Out_p_sxx_t[_pb->Outp]);
		_plt_s_e->SetXY (1, &_pb->Out_p_eyy_t[_pb->Outp], &_pb->Out_p_syy_t[_pb->Outp]);
		_plt_v_u->SetXY (0, &_pb->Out_p_ux_t[_pb->Outp], &_pb->Out_p_vx_t[_pb->Outp]);
		_plt_v_u->SetXY (1, &_pb->Out_p_uy_t[_pb->Outp], &_pb->Out_p_vy_t[_pb->Outp]);
		_plt_f_t->SetXY (0, &_pb->Out_t, &_pb->Out_p_fx_t[_pb->Outp]);
		_plt_f_t->SetXY (1, &_pb->Out_t, &_pb->Out_p_fy_t[_pb->Outp]);
	}
	void draw ()
	{
		_plt_f_u->draw();
		_plt_s_e->draw();
		_plt_v_u->draw();
		_plt_f_t->draw();
		fl_color (0xefebe700);
		fl_rectf (x()       , y()+h()-2 , w() , 3);
		fl_rectf (x()+w()-1 , y()       , 3   , h());
	}
private:
	PlotXY          * _plt_f_u; // plot force x displacement
	PlotXY          * _plt_s_e; // plot stress x strain
	PlotXY          * _plt_v_u; // plot velocity x displacement
	PlotXY          * _plt_f_t; // plot force x time
	Problem2D const * _pb;      // problem
}; // class TLframe

/////////////////////////////////////////////////////////////////////////////////////////// TRframe /////

class TRframe : public Fl_Group
{
public:
	TRframe (int X, int Y, int W, int H, Problem2D const * Prob) : Fl_Group (X,Y,W,H,NULL), _pb(Prob)
	{
		int hp = (H-1)/2;
		_plt_E_t = new PlotXY (X+1 , Y   , W-1 , hp , "Energy x Time"       , "Time" , "Energy");
		_plt_u_t = new PlotXY (X+1 , Y+hp, W-1 , hp , "Displacement x Time" , "Time" , "Displ.");
		end();

		// Set Plot energy
		_plt_E_t->AddCurve ("strain");
		_plt_E_t->AddCurve ("kinetic");
		_plt_E_t->AddCurve ("total");
		_plt_E_t->SetCurve (0).Typ = CT_LINES;
		_plt_E_t->SetCurve (0).Clr = FL_BLUE;
		_plt_E_t->SetCurve (1).Typ = CT_LINES;
		_plt_E_t->SetCurve (1).Clr = FL_GREEN;
		_plt_E_t->SetCurve (2).Typ = CT_LINES;
		_plt_E_t->SetCurve (2).Clr = FL_RED;

		// Set Plot displacement
		_plt_u_t->AddCurve ("ux");
		_plt_u_t->AddCurve ("uy");
		_plt_u_t->SetCurve (0).Typ = CT_LINES;
		_plt_u_t->SetCurve (0).Clr = FL_BLUE;
		_plt_u_t->SetCurve (1).Typ = CT_LINES;
		_plt_u_t->SetCurve (1).Clr = FL_GREEN;
	}
	void SetCurves ()
	{
		_plt_E_t->SetXY (0, &_pb->Out_t, &_pb->Out_sE_t);
		_plt_E_t->SetXY (1, &_pb->Out_t, &_pb->Out_kE_t);
		_plt_E_t->SetXY (2, &_pb->Out_t, &_pb->Out_tE_t);
		_plt_u_t->SetXY (0, &_pb->Out_t, &_pb->Out_p_ux_t [_pb->Outp]);
		_plt_u_t->SetXY (1, &_pb->Out_t, &_pb->Out_p_uy_t [_pb->Outp]);
	}
	void draw ()
	{
		_plt_E_t->draw();
		_plt_u_t->draw();
		fl_color (0xefebe700);
		fl_rectf (x(),y()+h()-2,w(),3);
	}
private:
	PlotXY          * _plt_E_t; // plot energy x time
	PlotXY          * _plt_u_t; // plot displacement x time
	Problem2D const * _pb;      // problem
}; // class TRframe

/////////////////////////////////////////////////////////////////////////////////////////// BRframe /////

class BRframe : public Fl_Group
{
public:
	BRframe (int X, int Y, int W, int H, Problem2D const * Prob) : Fl_Group (X,Y,W,H,NULL), _pb(Prob)
	{
		int hp = (H-1)/3;
		_plt_v_t = new PlotXY (X+1, Y+1     , W-1 , hp, "Velocity x Time", "Time", "Velocity");
		_plt_s_t = new PlotXY (X+1, Y+1+  hp, W-1 , hp, "Stress x Time",   "Time", "Stress"  );
		_plt_e_t = new PlotXY (X+1, Y+1+2*hp, W-1 , hp, "Strain x Time",   "Time", "Strain"  );
		end();

		// Set Plot velocity
		_plt_v_t->AddCurve ("vx");
		_plt_v_t->AddCurve ("vy");
		_plt_v_t->AddCurve ("cvx");
		_plt_v_t->SetCurve (0).Typ = CT_LINES;
		_plt_v_t->SetCurve (0).Clr = FL_BLUE;
		_plt_v_t->SetCurve (1).Typ = CT_LINES;
		_plt_v_t->SetCurve (1).Clr = FL_GREEN;
		_plt_v_t->SetCurve (2).Typ = CT_LINES;
		_plt_v_t->SetCurve (2).Clr = FL_RED;

		// Set Plot stress
		_plt_s_t->AddCurve ("sxx");
		_plt_s_t->AddCurve ("syy");
		_plt_s_t->SetCurve (0).Typ = CT_LINES;
		_plt_s_t->SetCurve (0).Clr = FL_BLUE;
		_plt_s_t->SetCurve (1).Typ = CT_LINES;
		_plt_s_t->SetCurve (1).Clr = FL_GREEN;

		// Set Plot strain
		_plt_e_t->AddCurve ("exx");
		_plt_e_t->AddCurve ("eyy");
		_plt_e_t->SetCurve (0).Typ = CT_LINES;
		_plt_e_t->SetCurve (0).Clr = FL_BLUE;
		_plt_e_t->SetCurve (1).Typ = CT_LINES;
		_plt_e_t->SetCurve (1).Clr = FL_GREEN;
	}
	void SetCurves ()
	{
		_plt_v_t->SetXY (0, &_pb->Out_t, &_pb->Out_p_vx_t [_pb->Outp]);
		_plt_v_t->SetXY (1, &_pb->Out_t, &_pb->Out_p_vy_t [_pb->Outp]);
		_plt_v_t->SetXY (2, &_pb->Out_t, &_pb->Out_p_cvx_t[_pb->Outp]);
		_plt_s_t->SetXY (0, &_pb->Out_t, &_pb->Out_p_sxx_t[_pb->Outp]);
		_plt_s_t->SetXY (1, &_pb->Out_t, &_pb->Out_p_syy_t[_pb->Outp]);
		_plt_e_t->SetXY (0, &_pb->Out_t, &_pb->Out_p_exx_t[_pb->Outp]);
		_plt_e_t->SetXY (1, &_pb->Out_t, &_pb->Out_p_eyy_t[_pb->Outp]);
	}
	void draw ()
	{
		_plt_v_t->draw();
		_plt_s_t->draw();
		_plt_e_t->draw();
	}
private:
	PlotXY          * _plt_v_t; // plot velocity x time
	PlotXY          * _plt_s_t; // plot stress x time
	PlotXY          * _plt_e_t; // plot strain x time
	Problem2D const * _pb;      // problem

}; // class BRframe

/////////////////////////////////////////////////////////////////////////////////////////// BLframe /////

class BLframe : public Fl_Group
{
public:
	BLframe (int X, int Y, int W, int H, Problem2D * Prob) : Fl_Group (X,Y,W,H,NULL), _pb(Prob),
	 _tl(NULL), _tr(NULL), _br(NULL)
	{
		_da = new DrawArea2D (X+1,Y+1,W-2,H-1-3);
		end();

		// Set Draw area
		_da->RecalcSF (true);
		_da->SetTime  (&(Prob->t)); // time
		_da->SetM     (&(Prob->M)); // multiplier for external forces
		_da->callback (_execute_selection_cb, this); // method to be called when a point/node is selected in the DrawArea2D
	}
	void ReSetGridPointsAndOutp ()
	{
		_da->SetGrid   (_pb->G);
		_da->SetPoints (_pb->P);
		_da->ResetSelP (_pb->Outp);
	}
	void SetFrames (Fl_Group * TL, Fl_Group * TR, Fl_Group * BR)
	{
		_tl = static_cast<TLframe*>(TL);
		_tr = static_cast<TRframe*>(TR);
		_br = static_cast<BRframe*>(BR);
	}
	void draw ()
	{
		_da->damage (FL_DAMAGE_EXPOSE);
		_da->draw   ();
		fl_color    (0xefebe700);
		fl_rectf    (x()+w()-1 , y(), 3, h());
	}
private:
	// Data
	DrawArea2D * _da;      // draw area
	Problem2D  * _pb;      // problem

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

	// Callbacks (must be static)
	static void _execute_selection_cb (Fl_Widget * o, void * v) { ((BLframe*)v)->_execute_selection(); }

	// Private methods
	void _execute_selection ()
	{
		_pb->Outp = _da->SelP(); // selected point
		if (_tl!=NULL) { _tl->SetCurves(); _tl->redraw(); }
		if (_tr!=NULL) { _tr->SetCurves(); _tr->redraw(); }
		if (_br!=NULL) { _br->SetCurves(); _br->redraw(); }
		Fl::wait(0);
	}

}; // class BLframe

///////////////////////////////////////////////////////////////////////////////////////// CallBacks /////

inline Fl_Group * AllocTL (int X, int Y, int W, int H, void * Prob) { return new TLframe (X,Y,W,H,static_cast<Problem2D*>(Prob)); }
inline Fl_Group * AllocTR (int X, int Y, int W, int H, void * Prob) { return new TRframe (X,Y,W,H,static_cast<Problem2D*>(Prob)); }
inline Fl_Group * AllocBL (int X, int Y, int W, int H, void * Prob) { return new BLframe (X,Y,W,H,static_cast<Problem2D*>(Prob)); }
inline Fl_Group * AllocBR (int X, int Y, int W, int H, void * Prob) { return new BRframe (X,Y,W,H,static_cast<Problem2D*>(Prob)); }

/////////////////////////////////////////////////////////////////////////////////////////// Control /////

class Control : public Fl_Group
{
public:
	// Constructor
	Control (int X, int Y, int W, int H, Problem2D & Prob) : Fl_Group (X,Y,W,H,NULL), _prob(Prob),
	   _running(false), _do_stop(false)
	{
		// Allocate widgets
		_w_prob    = new Fl_Input       ( 50 , 2 , 30  , 25 , "Prob#");
		_w_npcell  = new Fl_Input       (135 , 2 , 30  , 25 , "NPcell");
		_w_mpm     = new Fl_Check_Button(185 , 2 , 30  , 25 , "MPM");
		_w_usf     = new Fl_Check_Button(245 , 2 , 30  , 25 , "USF");
		_w_smallst = new Fl_Check_Button(300 , 2 , 30  , 25 , "SmallST");
		_w_tf      = new Fl_Input       (410 , 2 , 30  , 25 , "tf");
		_w_dt      = new Fl_Input       (470 , 2 , 60  , 25 , "Dt");
		_w_dtout   = new Fl_Input       (585 , 2 , 60  , 25 , "Out Dt");
		_w_setprob = new Fl_Button      (671 , 2 , 110 , 25 , "(Re)set &Prob");
		_w_refine  = new Fl_Button      (781 , 2 , 80  , 25 , "Re&fine");
		_w_step    = new Fl_Button      (861 , 2 , 80  , 25 , "&Step");
		_w_run     = new Fl_Button      (941 , 2 , 80  , 25 , "&Run" );
		end();

		// Initialize widgets
		char buf[256];
		snprintf(buf, 256, "%d", _prob.Id    ); _w_prob   ->value(buf);
		snprintf(buf, 256, "%d", _prob.NPcell); _w_npcell ->value(buf);
		snprintf(buf, 256, "%g", _prob.tf    ); _w_tf     ->value(buf);