mpm1d.cpp

crack modeling with xfem

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

// STL
#include <iostream>
#include <cmath>
#include <vector>
#include <cfloat>    // for DBL_EPSILON
#include <algorithm> // for min_element and max_element

// FLTK
#include <FL/Fl.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_draw.H>
#include <FL/Fl_Double_Window.H>

// Blitz++
#include <blitz/tinyvec-et.h>
#include <blitz/tinymat.h>

// Constants
const double SQ2 = sqrt(2.0);
const double PI  = 3.141592653589793238462643;

// Structures
typedef blitz::TinyVector<int   ,4> MapsAry; // for points/nodes maps
typedef blitz::TinyVector<double,4> ShapAry; // for shape/derivs. S and G

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

// Gravity callbacks
static inline double Prob1B (double t) { return 0.0;      }
static inline double Prob2B (double t) { return 0.0;      }
static inline double Prob3B (double t) { return -400.0*t; }

// Correct velocities at the center of mass
static inline double Prob1CorrectVelCM (double t) { return 0.1*cos(2.0*PI*t);    }
static inline double Prob2CorrectVelCM (double t) { return 0.2*cos(PI*t/5.0)/PI; }
static inline double Prob3CorrectVelCM (double t) { return 0.0;                  }

// Correct stresses
static inline double Prob1CorrectStress (double t, double x) { return 0.0; }
static inline double Prob2CorrectStress (double t, double x) { return 0.0; }
static inline double Prob3CorrectStress (double t, double x)
{
	double b   = Prob3B(t);
	double len = 50.0+b*(1250.0/1.0e+6);
	return 1.0e+6*(sqrt(2.0*b*(len-x)/1.0e+6+1.0)-1.0);
}

// Typedefs for callback functions
typedef double (*ptB        )(double t);
typedef double (*ptCorVelCM )(double t);
typedef double (*ptCorStress)(double t, double x);

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

// MainWindow
class MainWindow : public Fl_Double_Window
{
public:
	// DrawArea structure
	class DrawArea : public Fl_Widget
	{ // {{{
	public:
		// Constructor
		DrawArea (int xmin, int ymin, int width, int height) // Screen coordinates
			: Fl_Widget (xmin,ymin,width,height,0),
			  _mw(NULL), _rn(0), _rp(0), _rr(0)
		{}

		// Get access to the MainWindow data
		void SetMainWindow   (MainWindow const * Mw) { _mw=Mw;                        }
		void SetRadiusPoints (double Rp)             { _rp=Rp; _rr=(_rp>_rn?_rp:_rn); }
		void SetRadiusNodes  (double Rn)             { _rn=Rn; _rr=(_rp>_rn?_rp:_rn); }

		// Internal methods
		void draw()
		{
			if (_mw==NULL)                      return;
			if (_mw->_p==NULL || _mw->_n==NULL) return;
			if (_mw->_nn<2)                     return;

			// Bounding box
			_Xmin = _mw->_n[0];
			_Ymin = -1.0;
			_Xmax = _mw->_n[1];
			_Ymax =  1.0;
			for (int n=2; n<_mw->_nn; ++n)
			{
				if (_mw->_n[n]<_Xmin) _Xmin = _mw->_n[n];
				if (_mw->_n[n]>_Xmax) _Xmax = _mw->_n[n];
			}
			_Xmin-=_rr; _Ymin-=_rr; _Xmax+=_rr; _Ymax+=_rr;

			// Scale factors
			double sfx = static_cast<double>((w()-2)/(_Xmax-_Xmin));
			double sfy = static_cast<double>((h()-2)/(_Ymax-_Ymin));
			_sf = (sfx>sfy ? sfy : sfx);

			// Clear the background
			fl_color (FL_WHITE);
			fl_rectf (x(),y(),w(),h());

			// Draw nodes
			fl_color (FL_BLACK);
			int r = _l(_rn);
			for (int n=0; n<_mw->_nn; ++n)
				fl_circle (_x(_mw->_n[n]), _y(0.0), r);

			// Draw material points
			r = _l(_rp);
			for (int p=0; p<_mw->_np; ++p)
			{
				fl_color (FL_RED);
				fl_pie   (_x(_mw->_p[p])-r, _y(0.0)-r, 2*r+1, 2*r+1, 0,360);
			}
		}

	private:
		// Data
		double             _Xmin; // Real coordinates
		double             _Ymin; // Real coordinates
		double             _Xmax; // Real coordinates
		double             _Ymax; // Real coordinates
		double             _rn;   // Radius to draw nodes
		double             _rp;   // Radius to draw material points
		double             _rr;   // Bigger radius between _rn and _rp
		double             _sf;   // 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)        +_sf*(X-_Xmin)); } // X in screen coordinates
		inline int _y (double   Y) const { return static_cast<int>((y()+1)+(h()-2)-_sf*(Y-_Ymin)); } // Y in screen coordinates
		inline int _l (double   L) const { return static_cast<int>(_sf*L)                        ; } // Length in screen coordinates

	}; // class DrawArea }}}

	// PlotVeloc structure
	class PlotVeloc : public Fl_Widget
	{ // {{{
	public:
		// Constructor
		PlotVeloc (int xmin, int ymin, int width, int height) // Screen coordinates
			: Fl_Widget (xmin,ymin,width,height,0),
			  _mw(NULL)
		{}

		// Get access to the MainWindow data
		void SetMainWindow (MainWindow const * Mw   ) { _mw   = Mw;   }
		void SetPtrVelCM   (ptCorVelCM         pVcm ) { _pvcm = pVcm; }

		// Internal methods
		void draw()
		{
			if (_mw==NULL)       return;
			if (_mw->_nt<1)      return;
			if (_mw->_Ovp==NULL) return;

			// Bounding box
			_Xmin = 0.0;
			_Xmax = _mw->_dt*_mw->_ti; //_nt;
			_Ymin = _mw->_Ovp[0];
			_Ymax = _mw->_Ovp[0];
			for (int ti=0; ti<_mw->_ti+1; ++ti)
			{
				int i = ti*_mw->_np+_mw->_outvp;
				if (_mw->_Ovp[i]<_Ymin) _Ymin = _mw->_Ovp[i];
				if (_mw->_Ovp[i]>_Ymax) _Ymax = _mw->_Ovp[i];
			}

			// Scale factors
			_sfx = static_cast<double>((w()-2)/(_Xmax-_Xmin));
			_sfy = static_cast<double>((h()-2)/(_Ymax-_Ymin));

			// Clear the background
			fl_color (FL_WHITE);
			fl_rectf (x(),y(),w(),h());

			// Plot velocity (numerical)
			fl_color (FL_BLUE);
			for (int ti=0; ti<_mw->_ti; ++ti)
			{
				int i =  ti   *_mw->_np+_mw->_outvp;
				int j = (ti+1)*_mw->_np+_mw->_outvp;
				int x0=_x(_mw->_dt* ti   );  int y0=_y(_mw->_Ovp[i]);
				int x1=_x(_mw->_dt*(ti+1));  int y1=_y(_mw->_Ovp[j]);
				fl_line  (x0,y0, x1,y1);
			}

			// Plot velocity (analytical)
			fl_color (FL_RED);
			for (int ti=0; ti<_mw->_ti; ++ti)
			{
				double t0 = _mw->_dt* ti;
				double t1 = _mw->_dt*(ti+1);
				double v0 = (*_pvcm)(t0);
				double v1 = (*_pvcm)(t1);
				fl_line (_x(t0),_y(v0), _x(t1),_y(v1));
			}
		}

	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
		ptCorVelCM         _pvcm; // Pointer to the analytical function that calculates the velocity at CM

		// 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 PlotVeloc }}}

	// PlotStress structure
	class PlotStress : public Fl_Widget
	{ // {{{
	public:
		// Constructor
		PlotStress (int xmin, int ymin, int width, int height) // Screen coordinates
			: Fl_Widget (xmin,ymin,width,height,0),
			  _mw(NULL), _ti(0)
		{}

		// Get access to the MainWindow data
		void SetMainWindow (MainWindow const * Mw) { _mw = Mw; }
		void SetPtrStress  (ptCorStress        pS) { _ps = pS; }
		void SetOutTimeInc (int                Ti) { _ti = Ti; }

		// Internal methods
		void draw()
		{
			if (_mw==NULL)       return;
			if (_mw->_nt<1)      return;
			if (_mw->_Osp==NULL) return;

			// Bounding box
			_Xmin = 0.0;
			_Xmax = _mw->_len0;
			_Ymin = _mw->_Osp[0];
			_Ymax = _mw->_Osp[0];
			for (int p=0; p<_mw->_np; ++p)
			{
				int i = _ti*_mw->_np+p;
				if (_mw->_Osp[i]<_Ymin) _Ymin = _mw->_Osp[i];
				if (_mw->_Osp[i]>_Ymax) _Ymax = _mw->_Osp[i];
			}

			// Scale factors
			_sfx = static_cast<double>((w()-2)/(_Xmax-_Xmin));
			_sfy = static_cast<double>((h()-2)/(_Ymax-_Ymin));

			// Clear the background
			fl_color (FL_WHITE);
			fl_rectf (x(),y(),w(),h());

			// Plot stress (numerical)
			fl_color (FL_BLUE);
			for (int p=0; p<_mw->_np-1; ++p)
			{
				int i = _ti*_mw->_np+p;
				int j = _ti*_mw->_np+p+1;
				int x0=_x(_mw->_p[p  ]);  int y0=_y(_mw->_Osp[i]);
				int x1=_x(_mw->_p[p+1]);  int y1=_y(_mw->_Osp[j]);
				fl_line  (x0,y0, x1,y1);
			}

			// Plot stress (analytical)
			fl_color (FL_RED);
			for (int p=0; p<_mw->_np-1; ++p)
			{
				double t  = _ti*_mw->_dt;
				double x0 = _mw->_p[p  ];
				double x1 = _mw->_p[p+1];
				double s0 = (*_ps)(t,x0);
				double s1 = (*_ps)(t,x1);
				fl_line (_x(x0),_y(s0), _x(x1),_y(s1));
			}
		}

	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
		ptCorStress        _ps;   // Pointer to the analytical function that calculates the stress
		int                _ti;   // Output time increment

		// 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 PlotStress }}}

	// PlotEnergy structure
	class PlotEnergy : public Fl_Widget
	{ // {{{
	public:
		// Constructor
		PlotEnergy (int xmin, int ymin, int width, int height) // Screen coordinates
			: Fl_Widget (xmin,ymin,width,height,0),
			  _mw(NULL)
		{}

		// Get access to the MainWindow data
		void SetMainWindow (MainWindow const * Mw ) { _mw  = Mw;  }

		// Internal methods
		void draw()
		{
			if (_mw==NULL)                          return;
			if (_mw->_nt<1)                         return;
			if (_mw->_Esp==NULL || _mw->_Ekp==NULL) return;

			// Bounding box
			_Xmin = 0.0;
			_Xmax = _mw->_dt*_mw->_ti;
			_Ymin = _mw->_Esp[0];
			_Ymax = _mw->_Esp[0];
			for (int ti=0; ti<_mw->_ti+1; ++ti)
			{
				if (_mw->_Esp[ti]<_Ymin) _Ymin = _mw->_Esp[ti];
				if (_mw->_Ekp[ti]<_Ymin) _Ymin = _mw->_Ekp[ti];
				if (_mw->_Esp[ti]>_Ymax) _Ymax = _mw->_Esp[ti];
				if (_mw->_Ekp[ti]>_Ymax) _Ymax = _mw->_Ekp[ti];
			}

			// Scale factors
			_sfx = static_cast<double>((w()-2)/(_Xmax-_Xmin));
			_sfy = static_cast<double>((h()-2)/(_Ymax-_Ymin));

			// Clear the background
			fl_color (FL_WHITE);
			fl_rectf (x(),y(),w(),h());

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