problems.h

crack modeling with xfem

	size_t      Outp;    ///< point/particle number to output data
	Vector3D    B;       ///< current b-body mass
	double      M;       ///< Current multiplier for external forces

	// Output arrays
	Array<double>         Out_t;       // time array
	Array<double>         Out_Wx_t;    // Body weights (in time)
	Array<double>         Out_sE_t;    // Out strain Energy (in time)
	Array<double>         Out_kE_t;    // Out kinetic Energy (in time)
	Array<double>         Out_tE_t;    // Out total Energy (in time)
	Array<Array<double> > Out_p_vx_t;  // Out particle x-velocity (in time) Size=np
	Array<Array<double> > Out_p_vy_t;  // Out particle y-velocity (in time) Size=np
	Array<Array<double> > Out_p_cvx_t; // Out correct particle x-velocity (in time) Size=np
	Array<Array<double> > Out_p_sxx_t; // Out particle xx-stress (in time) Size=np
	Array<Array<double> > Out_p_syy_t; // Out particle yy-stress (in time) Size=np
	Array<Array<double> > Out_p_exx_t; // Out particle xx-strain (in time) Size=np
	Array<Array<double> > Out_p_eyy_t; // Out particle yy-strain (in time) Size=np
	Array<Array<double> > Out_p_fx_t;  // Out external force x (in time) Size=np
	Array<Array<double> > Out_p_fy_t;  // Out external force y (in time) Size=np
	Array<Array<double> > Out_p_ux_t;  // Out displacement x (in time) Size=np
	Array<Array<double> > Out_p_uy_t;  // Out displacement y (in time) Size=np
};

/** (Re)set output arrays. */
void ReSetOutputArrays (Problem2D & Prob)
{
	// Erase output arrays
	Prob.Out_t      .Resize(0);
	Prob.Out_Wx_t   .Resize(0);
	Prob.Out_sE_t   .Resize(0);
	Prob.Out_kE_t   .Resize(0);
	Prob.Out_tE_t   .Resize(0);
	Prob.Out_p_vx_t .Resize(Prob.P->nPoints());
	Prob.Out_p_vy_t .Resize(Prob.P->nPoints());
	Prob.Out_p_cvx_t.Resize(Prob.P->nPoints());
	Prob.Out_p_sxx_t.Resize(Prob.P->nPoints());
	Prob.Out_p_syy_t.Resize(Prob.P->nPoints());
	Prob.Out_p_exx_t.Resize(Prob.P->nPoints());
	Prob.Out_p_eyy_t.Resize(Prob.P->nPoints());
	Prob.Out_p_fx_t .Resize(Prob.P->nPoints());
	Prob.Out_p_fy_t .Resize(Prob.P->nPoints());
	Prob.Out_p_ux_t .Resize(Prob.P->nPoints());
	Prob.Out_p_uy_t .Resize(Prob.P->nPoints());
}

/** Set problem. */
void SetProblem (Problem2D & Prob)
{
	// Deallocate previous grid and matpoints
	if (Prob.G!=NULL) delete Prob.G;
	if (Prob.P!=NULL) delete Prob.P;

	// Set problem data
	switch (Prob.Id)
	{
		case 1:
		{
			// Allocate grid
			Prob.G = new Grid2D (/*xMin*/0, /*yMin*/0, /*nRow*/4, /*nCol*/4, /*Lx*/1.0, /*Ly*/0.1, &CbIsFixed1);

			// Allocate material points
			Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom1, &CbDensity1, &CbAllocMdl1, &CbIniVeloc1, &CbHasTraction1);

			// Set callbacks
			Prob.pB      = &CbB1;
			Prob.pM      = &CbLdM1;
			Prob.pVeloc  = &CbVel1;
			Prob.pStress = &CbStress1;

			// Set constants
			Prob.t     = 0.0;          // current time
			Prob.tf    = 5.0;          // final time
			Prob.Dt    = 0.1/(2.0*PI); // time increment
			Prob.Dtout = Prob.Dt;      // time increment for output
			Prob.Outp  = 0;            // point/particle number to output data
			Prob.B     = 0.0;          // current b-body mass
			Prob.M     = 0.0;          // multiplier for external forces

			// Boundary conditions
			for (int i=0; i<16; ++i) Prob.G->SetFixed (i, FIX_Y);
			Prob.G->SetFixed( 0, FIX_XY);  Prob.G->SetFixed( 1, FIX_XY);
			Prob.G->SetFixed( 4, FIX_XY);  Prob.G->SetFixed( 5, FIX_XY);
			Prob.G->SetFixed( 8, FIX_XY);  Prob.G->SetFixed( 9, FIX_XY);
			Prob.G->SetFixed(12, FIX_XY);  Prob.G->SetFixed(13, FIX_XY);

			break;
		}
		case 2:
		{
			// Allocate grid
			Prob.G = new Grid2D (/*xMin*/0, /*yMin*/0, /*nRow*/9, /*nCol*/8, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed2);

			// Allocate material points
			Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom2, &CbDensity2, &CbAllocMdl2, &CbIniVeloc2, &CbHasTraction2);

			// Set callbacks
			Prob.pB      = &CbB2;
			Prob.pM      = &CbLdM2;
			Prob.pVeloc  = &CbVel2;
			Prob.pStress = &CbStress2;

			// Set constants
			Prob.t     = 0.0;   // current time
			Prob.tf    = 1.0;   // final time
			Prob.Dt    = 0.001; // time increment
			Prob.Dtout = 0.05;  // time increment for output
			Prob.Outp  = 0;     // point/particle number to output data
			Prob.B     = 0.0;   // current b-body mass
			Prob.M     = 0.0;   // multiplier for external forces

			// Boundary conditions
			Vector3D l; l = 0.0, 1.0, 0.0;
			for (int i=0; i<Prob.G->nNodes(); ++i)
			{
				if (Prob.G->N(i)(0)>0.9 && Prob.G->N(i)(0)<6.1
				 && Prob.G->N(i)(1)>0.9 && Prob.G->N(i)(1)<7.1)
				{
					if (Prob.G->N(i)(0)>0.9 && Prob.G->N(i)(0)<1.1) Prob.G->SetFixed   (i, FIX_X);
					if (Prob.G->N(i)(1)>0.9 && Prob.G->N(i)(1)<1.1) Prob.G->SetFixed   (i, FIX_Y);
					if (Prob.G->N(i)(0)>0.9 && Prob.G->N(i)(0)<1.1
					 && Prob.G->N(i)(1)>0.9 && Prob.G->N(i)(1)<1.1) Prob.G->SetFixed   (i, FIX_XY);
					if (Prob.G->N(i)(1)>6.9 && Prob.G->N(i)(1)<7.1) Prob.G->SetLoading (i, l);
				}
			}

			break;
		}
		case 3:
		{
			// Allocate grid
			Prob.G = new Grid2D (/*xMin*/-0.10, /*yMin*/-0.10, /*nRow*/24, /*nCol*/24, /*Lx*/0.05, /*Ly*/0.05, &CbIsFixed3);

			// Allocate material points
			Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom3, &CbDensity3, &CbAllocMdl3, &CbIniVeloc3, &CbHasTraction3);

			// Set callbacks
			Prob.pB      = &CbB3;
			Prob.pM      = &CbLdM3;
			Prob.pVeloc  = &CbVel3;
			Prob.pStress = &CbStress3;

			// Set constants
			Prob.t     = 0.0;   // current time
			Prob.tf    = 10.0;  // final time
			Prob.Dt    = 0.001; // time increment
			Prob.Dtout = 0.05;  // time increment for output
			Prob.Outp  = 51;    // point/particle number to output data
			Prob.B     = 0.0;   // current b-body mass
			Prob.M     = 0.0;   // multiplier for external forces

			// Boundary conditions
			for (int i=0; i<Prob.G->nNodes(); ++i)
			{
				if (Prob.G->N(i)(0)<0.10) Prob.G->SetFixed(i,FIX_XY); // left
				if (Prob.G->N(i)(0)>0.89) Prob.G->SetFixed(i,FIX_XY); // right
				if (Prob.G->N(i)(1)<0.10) Prob.G->SetFixed(i,FIX_XY); // bottom
				if (Prob.G->N(i)(1)>0.89) Prob.G->SetFixed(i,FIX_XY); // top
			}

			break;
		}
		case 4:
		{
			// Allocate grid
			Prob.G = new Grid2D (/*xMin*/-0.25, /*yMin*/-0.25, /*nRow*/40, /*nCol*/30, /*Lx*/0.05, /*Ly*/0.05, &CbIsFixed4);

			// Allocate material points
			Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom4, &CbDensity4, &CbAllocMdl4, &CbIniVeloc4, &CbHasTraction4);

			// Set callbacks
			Prob.pB      = &CbB4;
			Prob.pM      = &CbLdM4;
			Prob.pVeloc  = &CbVel4;
			Prob.pStress = &CbStress4;

			// Set constants
			Prob.t     = 0.0;   // current time
			Prob.tf    = 40.0;  // final time
			Prob.Dt    = 0.001; // time increment
			Prob.Dtout = 0.05;  // time increment for output
			Prob.Outp  = 0;     // point/particle number to output data
			Prob.B     = 0.0;   // current b-body mass
			Prob.M     = 0.0;   // multiplier for external forces

			break;
		}
		case 5:
		{
			// Allocate grid
			Prob.G = new Grid2D (/*xMin*/-0.25, /*yMin*/-0.25, /*nRow*/7, /*nCol*/7, /*Lx*/0.25, /*Ly*/0.25, &CbIsFixed5);

			// Allocate material points
			Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom5, &CbDensity5, &CbAllocMdl5, &CbIniVeloc5, &CbHasTraction5);

			// Set callbacks
			Prob.pB      = &CbB5;
			Prob.pM      = &CbLdM5;
			Prob.pVeloc  = &CbVel5;
			Prob.pStress = &CbStress5;

			// Set constants
			Prob.t     = 0.0;    // current time
			Prob.tf    = 4.0;    // final time
			Prob.Dt    = 0.0005; // time increment
			Prob.Dtout = 0.01;   // time increment for output
			Prob.Outp  = 0;      // point/particle number to output data
			Prob.B     = 0.0;    // current b-body mass
			Prob.M     = 0.0;    // multiplier for external forces

			break;
		}
		default: { throw new Fatal("Prob # %d is not available.",Prob.Id); }
	}
	ReSetOutputArrays (Prob);
}

#endif // MPM_PROBLEMS2D_H

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