problems.h

crack modeling with xfem

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

#ifndef MPM_PROBLEMS2D_H
#define MPM_PROBLEMS2D_H

// Local
#include "defs.h"
#include "array.h"
#include "grid2d.h"
#include "model2d.h"
#include "mpoints2d.h"

// Fixed nodes?
bool CbIsFixed1 (Vector3D const & N, FixType & FType) { return false; }
bool CbIsFixed2 (Vector3D const & N, FixType & FType) { return false; }
bool CbIsFixed3 (Vector3D const & N, FixType & FType) { return false; }
bool CbIsFixed4 (Vector3D const & N, FixType & FType) { return false; }
bool CbIsFixed5 (Vector3D const & N, FixType & FType)
{
	if (N(0)>-0.01 && N(0)< 0.01 && N(1)>-0.01 && N(1)<0.01) { FType=FIX_XY; return true; }
	if (N(0)>-0.01 && N(0)< 0.01 && N(1)> 0.01 && N(1)<1.01) { FType=FIX_X;  return true; }
	if (N(0)< 1.01 && N(0)>-0.01 && N(1)>-0.01 && N(1)<0.01) { FType=FIX_Y;  return true; }
	return false;
}

// Is point inside geometry callbacks
bool CbIsPointInGeom1 (Vector3D const & P)
{

	if (P(0)>=1.0 && P(0)<=2.0 && P(1)>=0.10 && P(1)<=0.2) return true;
	else                                                   return false;

}
bool CbIsPointInGeom2 (Vector3D const & P)
{
	if (P(0)>0.9 && P(0)<6.1 && P(1)>0.9 && P(1)<7.1) return true;
	else                                              return false;
}
bool CbIsPointInGeom3 (Vector3D const & P)
{
	// Cylinders
	double xc1=0.25; double yc1=0.25; double r1=0.14;
	double xc2=0.65; double yc2=0.65; double r2=0.14;

	// Bottom-left cylinder
	double d1 = sqrt(pow(P(0)-xc1,2.0)+pow(P(1)-yc1,2.0));

	// Upper-right cylinder
	double d2 = sqrt(pow(P(0)-xc2,2.0)+pow(P(1)-yc2,2.0));

	if (d1<=r1 || d2<=r2) return true;
	else
	{
		if (P(0)<0.05 && P(0)>0    && P(1)>0    && P(1)<0.95) return true; // left
		if (P(0)>0.9  && P(0)<0.95 && P(1)>0    && P(1)<0.95) return true; // right
		if (P(0)>0    && P(0)<0.95 && P(1)<0.05 && P(1)>0   ) return true; // bottom
		if (P(0)>0    && P(0)<0.95 && P(1)>0.9  && P(1)<0.95) return true; // top
		return false;
	}
}
bool CbIsPointInGeom4 (Vector3D const & P)
{
	if (P(0)>0.0 && P(0)<1.0 && P(1)>0.0 && P(1)<1.0) return true;
	else                                              return false;
}
bool CbIsPointInGeom5 (Vector3D const & P)
{
	if (P(0)>0.0 && P(0)<1.0 && P(1)>0.0 && P(1)<1.0) return true;
	else                                              return false;
}

// Density callbacks
double CbDensity1 (Vector3D const & P) { return 1.0; }
double CbDensity2 (Vector3D const & P) { return 1.0; }
double CbDensity3 (Vector3D const & P) { return 1.0; }
double CbDensity4 (Vector3D const & P) { return 1.0; }
double CbDensity5 (Vector3D const & P) { return 1.0; }

// Allocate model callbacks
Model2D * CbAllocMdl1 (Vector3D const & P)
{
	Array<double> prms(2);
	Array<double> inis(6);
	prms[0] = 4.0*PI*PI; // E
	prms[1] = 0.0;       // nu
	for (int i=0; i<6; ++i) inis[i] = 0.0; // Sig_ini
	return new LinEls2D (prms, inis);
}
Model2D * CbAllocMdl2 (Vector3D const & P)
{
	Array<double> prms(2);
	Array<double> inis(6);
	prms[0] = 10000.0;    // E
	prms[1] = 0.25;       // nu
	for (int i=0; i<6; ++i) inis[i] = 0.0; // Sig_ini
	return new LinEls2D (prms, inis);
}
Model2D * CbAllocMdl3 (Vector3D const & P)
{
	Array<double> prms(2);
	Array<double> inis(6);
	prms[0] = 1.0;    // E
	prms[1] = 0.2;    // nu
	for (int i=0; i<6; ++i) inis[i] = 0.0; // Sig_ini
	return new LinEls2D (prms, inis);
}
Model2D * CbAllocMdl4 (Vector3D const & P)
{
	Array<double> prms(2);
	Array<double> inis(6);
	prms[0] = 2.1;  // E
	prms[1] = 0.3;  // nu
	for (int i=0; i<6; ++i) inis[i] = 0.0; // Sig_ini
	return new LinEls2D (prms, inis);
}
Model2D * CbAllocMdl5 (Vector3D const & P)
{
	Array<double> prms(2);
	Array<double> inis(6);
	prms[0] = 210.0; // E
	prms[1] = 0.3;   // nu
	for (int i=0; i<6; ++i) inis[i] = 0.0; // Sig_ini
	return new LinEls2D (prms, inis);
}

// Initial velocity
void CbIniVeloc1 (Vector3D const & P, Vector3D & v)
{
	if (P(0)>=1.0 && P(0)<=2.0 && P(1)>=0.10 && P(1)<=0.2) v = 0.1, 0.0, 0.0;
	else                                                   v = 0.0, 0.0, 0.0;
}
void CbIniVeloc2 (Vector3D const & P, Vector3D & v) { v = 0.0; }
void CbIniVeloc3 (Vector3D const & P, Vector3D & v)
{
	/*
	     if (P(0)<0.5 && P(1)<0.5 && P(0)>0.1 && P(1)>0.1) v =  0.15,  0.15, 0.0;
	else if (P(0)>0.5 && P(1)>0.5 && P(0)<0.9 && P(1)<0.9) v = -0.15, -0.1, 0.0;
	else                                                   v =  0.0,  0.0, 0.0;
	*/
	     if (P(0)<0.5 && P(1)<0.5 && P(0)>0.1 && P(1)>0.1) v =  0.1,  0.1, 0.0;
	else if (P(0)>0.5 && P(1)>0.5 && P(0)<0.9 && P(1)<0.9) v = -0.1, -0.1, 0.0;
	else                                                   v =  0.0,  0.0, 0.0;
}
void CbIniVeloc4 (Vector3D const & P, Vector3D & v) { v = 0.0; }
void CbIniVeloc5 (Vector3D const & P, Vector3D & v) { v = 0.0; }

bool CbHasTraction1 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; }
bool CbHasTraction2 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; }
bool CbHasTraction3 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; }
bool CbHasTraction4 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t)
{
	if (P(1)>0.87)
	{
		if (nPCell==1) t = 0.0, 0.25 , 0.0;
		else           t = 0.0, 0.125, 0.0;
		return true;
	}
	else return false;
}
bool CbHasTraction5 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t)
{
	if (N[0](1)>0.99-L(1))
	{
		if (P(1)>N[0](1)+L(1)/2.01)
		{
			if (nPCell==1) t = 0.0, 0.25 , 0.0;
			else           t = 0.0, 0.125, 0.0;
			return true;
		}
	}
	return false;
}

// Body mass callbacks
static inline void CbB1 (double t, Vector3D & B) { B=0.0; }
static inline void CbB2 (double t, Vector3D & B) { B=0.0; }
static inline void CbB3 (double t, Vector3D & B) { B=0.0; }
static inline void CbB4 (double t, Vector3D & B) { B=0.0; }
static inline void CbB5 (double t, Vector3D & B) { B=0.0; }

// Loading multiplier callbacks
static inline void CbLdM1 (double t, double & M) { M=0.0;    }
static inline void CbLdM2 (double t, double & M) { M=(t<5.0 ? t/5.0 : 1.0); }
static inline void CbLdM3 (double t, double & M) { M=0.0;    }
static inline void CbLdM4 (double t, double & M) { M=(t<5.0 ? t/5.0 : 1.0); }
static inline void CbLdM5 (double t, double & M) { M=(t<0.5 ? t/0.5 : 1.0); }

// Correct velocities callbacks
static inline void CbVel1 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.1*cos(2.0*PI*t),0,0; }
static inline void CbVel2 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; }
static inline void CbVel3 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; }
static inline void CbVel4 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; }
static inline void CbVel5 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; }

// Correct stresses callbacks
static inline void CbStress1 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; }
static inline void CbStress2 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; }
static inline void CbStress3 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; }
static inline void CbStress4 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; }
static inline void CbStress5 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; }

/** Problem related data. */
struct Problem2D
{
	// Data
	int         Id;      ///< Problem number
	int         NPcell;  ///< Number of points per cell
	Grid2D    * G;       ///< The grid
	MPoints2D * P;       ///< The material points
	ptB         pB;      ///< Pointer to a function which returns the body mass B
	ptLdM       pM;      ///< Pointer to a function which returns the loadings multiplier M
	ptVeloc     pVeloc;  ///< Pointer to a function which returns the correct velocity
	ptStress    pStress; ///< Pointer to a function which returns the correct stress
	bool        MPM;     ///< use original MPM
	bool        USF;     ///< update stress first ?
	bool        SmallSt; ///< small strains ?
	double      t;       ///< current time
	double      tf;      ///< final time
	double      Dt;      ///< time increment
	double      Dtout;   ///< time increment for output