mechmat.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

   @param d - stiffness %matrix
   @param ipp - number of integration point
   
   JK, 17.7.2001
*/
void mechmat::elmatstiff (matrix &d,long ipp)
  //  27.10.2001
{
  long i,ncomp,idem;

  idem = ip[ipp].gemid();
  switch (ip[ipp].tm[idem]){
  case elisomat:{
    i=ip[ipp].idm[idem];
    eliso[i].elmatstiff (d,ip[ipp].ssst);
    break;
  }
  case elgmat3d:{
    i=ip[ipp].idm[idem];
    ncomp=ip[ipp].ncompstr;
    elgm3d[i].elmatstiff (d);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown elastic material type is required");
    fprintf (stderr,"\n in function elmatstiff (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}

/**
   function computes elastic compliance %matrix of the material
   in the required integration point
   
   @param c - compliance %matrix
   @param ipp - number of integration point
   
   JK, 3.2.2003
*/
void mechmat::elmatcompl (matrix &c,long ipp)
  //  27.10.2001
{
  long i, idem;

  idem = ip[ipp].gemid();
  switch (ip[ipp].tm[idem]){
  case elisomat:{
    i=ip[ipp].idm[idem];
    eliso[i].matcompl (c,ip[ipp].ssst);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown elastic material type is required");
    fprintf (stderr,"\n in function elmatcompl (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}


// *****************************************************************
// *****************************************************************
//   PART CONTAINING FUNCTIONS DEALING WITH PLASTICITY
// *****************************************************************
// *****************************************************************

/**
   function evaluates yield function for stresses and internal variables
   
   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param q - internal variables (hardening parameters)
   
   28.10.2001
*/
double mechmat::yieldfunction (long ipp, long idpm, matrix &sig,vector &q)
{
  double f;
  
  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    f = spl1d[ip[ipp].idm[idpm]].yieldfunction (sig,q);
    break;
  }
  case jflow:{
    f = j2f[ip[ipp].idm[idpm]].yieldfunction (sig,q);
    break;
  }
  case jflow2:{
    f = j2f2[ip[ipp].idm[idpm]].yieldfunction (sig,q);
    break;
  }
  case mohrcoulparab:{
    f = mcpar[ip[ipp].idm[idpm]].yieldfunction (sig);
    break;
  }
  case boermaterial:{
    f = boerm[ip[ipp].idm[idpm]].yieldfunction (sig);
    break;
  }
  case druckerprager:{
    f = drprm[ip[ipp].idm[idpm]].yieldfunction (sig,q);
    break;
  }
  case chenplast:{
    f = chplast[ip[ipp].idm[idpm]].yieldfunction (sig,q);
    break;
  }
  case modcamclaymat:{
    f = cclay[ip[ipp].idm[idpm]].yieldfunction (sig,q);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown plasticity model is required in function");
    fprintf (stderr,"\n yieldfunction (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  return f;
}

/**
   function evaluates derivates of yield function with respect to stresses

   stresses must be assembled in tensor notation and stored in 3x3 %matrix
   function overwrites array sig, some material models use the array sig for auxiliary computations
   function stores derivatives in %vector dfds

   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param q - internal parameters (hardening parameters)
   @param dfds - derivatives of yield function with respect to stresses
   
   28.10.2001
*/
void mechmat::dfdsigma (long ipp, long idpm, matrix &sig, vector &q, vector &dfds)
{
  matrix dfdst(3,3);
  
  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    spl1d[ip[ipp].idm[idpm]].dfdsigma (sig,dfdst);
    break;
  }
  case jflow:{
    j2f[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
    break;
  }
  case jflow2:{
    j2f2[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
    break;
  }
  case mohrcoulparab:{
    mcpar[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
    break;
  }
  case boermaterial:{
    boerm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
    break;
  }
  case druckerprager:{
    drprm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
    break;
  }
  case chenplast:{
    chplast[ip[ipp].idm[idpm]].deryieldfdsigma (sig,q,dfdst);
    break;
  }
  case modcamclaymat:{
    cclay[ip[ipp].idm[idpm]].deryieldfsigma (sig, q, dfdst);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown plasticity model is required in function");
    fprintf (stderr,"\n dfdsigma (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
  //  conversion of tensor notation to vector notation
  tensor_vector (dfds,dfdst,ip[ipp].ssst,stress);
  
}

/**
   function evaluates derivates of plastic potential with respect to stresses

   stresses must be assembled in tensor notation and stored in 3x3 %matrix
   function overwrites array sig, some material models use the array sig for auxiliary computations
   function stores derivatives in %vector dgds

   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param q - internal parameters (hardening parameters)
   @param dfds - derivatives of plastic potential with respect to stresses
   
   28.10.2001
*/
void mechmat::dgdsigma (long ipp, long idpm, matrix &sig, vector &q, vector &dgds)
{
  matrix dgdst(3,3);
  
  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    spl1d[ip[ipp].idm[idpm]].dfdsigma (sig,dgdst);
    break;
  }
  case jflow:{
    j2f[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst);
    break;
  }
  case jflow2:{
    j2f2[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst);
    break;
  }
  case mohrcoulparab:{
    mcpar[ip[ipp].idm[idpm]].derplaspotsigma (sig,dgdst);
    break;
  }
  case boermaterial:{
    boerm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst);
    break;
  }
  case druckerprager:{
    drprm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst);
    break;
  }
  case chenplast:{
    chplast[ip[ipp].idm[idpm]].deryieldfdsigma (sig,q,dgdst);
    break;
  }
  case modcamclaymat:{
    cclay[ip[ipp].idm[idpm]].deryieldfsigma (sig, q, dgdst);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown plasticity model is required in function");
    fprintf (stderr,"\n dgdsigma (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
  //  conversion of tensor notation to vector notation
  tensor_vector (dgds,dgdst,ip[ipp].ssst,stress);
  
}

/**
   function evaluates second derivates of yield function with respect to stresses
   
   stresses must be assembled in tensor notation and stored in 3x3 %matrix
   function stores derivatives in %matrix dfdsds
   
   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param dfdsds - second derivatives of yield function with respect to stresses
   
   28.10.2001
*/
void mechmat::dfdsigmadsigma (long ipp, long idpm, matrix &sig, matrix &dfdsds)
{
  matrix dfdsdst(6,6);
  
  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    spl1d[ip[ipp].idm[idpm]].dfdsigmadsigma (dfdsdst);
    break;
  }
  case jflow2:{
    j2f2[ip[ipp].idm[idpm]].deryielddsigmadsigma (sig,dfdsdst);
    break;
  }
  case modcamclaymat:{
    cclay[ip[ipp].idm[idpm]].dderyieldfsigma (dfdsdst);
    break;
  }
  case chenplast:{
    chplast[ip[ipp].idm[idpm]].deryieldfdsigmadsigma (sig,dfdsdst);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown plasticity model is required in function");
    fprintf (stderr,"\n dfdsigmadsigma (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
  //  conversion from tensor notation to matrix notation
  tensor4_matrix (dfdsds,dfdsdst,ip[ipp].ssst);
  
}

/**
   function evaluates second derivates of plastic potential with respect to stresses
   
   stresses must be assembled in tensor notation and stored in 3x3 %matrix
   function stores derivatives in %matrix dgdsds
   
   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param dfdsds - second derivatives of plastic potential with respect to stresses
   
   28.10.2001
*/
void mechmat::dgdsigmadsigma (long ipp, long idpm, matrix &sig, matrix &dgdsds)
{
  matrix dgdsdst(6,6);
  
  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    spl1d[ip[ipp].idm[idpm]].dfdsigmadsigma (dgdsdst);
    break;
  }
  case jflow2:{
    j2f2[ip[ipp].idm[idpm]].deryielddsigmadsigma (sig,dgdsdst);
    break;
  }
  case modcamclaymat:{
    cclay[ip[ipp].idm[idpm]].dderyieldfsigma (dgdsdst);
    break;
  }
  case chenplast:{
    chplast[ip[ipp].idm[idpm]].deryieldfdsigmadsigma (sig,dgdsdst);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown plasticity model is required in function");
    fprintf (stderr,"\n dgdsigmadsigma (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
  //  conversion from tensor notation to matrix notation
  tensor4_matrix (dgdsds,dgdsdst,ip[ipp].ssst);
  
}


/**
   function evaluates derivates of yield function with respect to internal variables (hardening parameters)
   
   stresses must be assembled in tensor notation and stored in 3x3 %matrix
   function stores derivatives in %vector dq
   
   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param q - internal variables (hardening parameters)
   @param dq - derivatives of yield function with respect to internal variables
   
   28.10.2001
*/
void mechmat::dfdqpar (long ipp, long idpm,matrix &sig,vector &q,vector &dq)
{
  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    spl1d[ip[ipp].idm[idpm]].dfdqpar (dq);
    break;
  }
  case jflow:{
    j2f[ip[ipp].idm[idpm]].deryieldfq (dq);
    break;
  }
  case jflow2:{
    j2f2[ip[ipp].idm[idpm]].deryieldfq (dq);
    break;
  }
  case mohrcoulparab:{
    break;
  }
  case boermaterial:{
    break;
  }
  case druckerprager:{
    break;
  }
  case modcamclaymat:{
    cclay[ip[ipp].idm[idpm]].deryieldfq (sig, dq);
    break;
  }
  case chenplast:{
    chplast[ip[ipp].idm[idpm]].deryieldfdq (sig,q,dq);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown plasticity model is required in function dfdqpar (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   function evaluates derivates of plastic potential with respect to internal variables (hardening parameters)
   
   stresses must be assembled in tensor notation and stored in 3x3 %matrix
   function stores derivatives in %vector dq
   
   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param q - internal variables (hardening parameters)
   @param dq - derivatives of plastic potential with respect to internal variables
   
   28.10.2001
*/
void mechmat::dgdqpar (long ipp, long idpm,matrix &sig,vector &q,vector &dq)
{
  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    spl1d[ip[ipp].idm[idpm]].dfdqpar (dq);
    break;
  }
  case jflow:{
    j2f[ip[ipp].idm[idpm]].deryieldfq (dq);
    break;
  }
  case jflow2:{
    j2f2[ip[ipp].idm[idpm]].deryieldfq (dq);
    break;
  }
  case mohrcoulparab:{
    break;
  }
  case boermaterial:{
    break;
  }
  case druckerprager:{
    break;
  }
  case modcamclaymat:{
    cclay[ip[ipp].idm[idpm]].deryieldfq (sig, dq);
    break;
  }
  case chenplast:{
    chplast[ip[ipp].idm[idpm]].deryieldfdq (sig,q,dq);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown plasticity model is required in function dgdqpar (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   function evaluates derivates of yield function with respect to stresses and
   internal variables (hardening parameters)

   stresses must be assembled in tensor notation and stored in 3x3 %matrix
   function stores derivatives in %matrix dfdsdq
   
   @param ipp - integration point pointer
   @param idpm - id of plasticity model
   @param sig - stress components
   @param q - internal variables (hardening parameters)
   @param dfdsdq - derivatives of yield function with respect to stresses and internal variables
   
   28.10.2001
*/
void mechmat::dfdsigmadq (long ipp, long idpm,matrix &sig,vector &q,matrix &dfdsdq)
{
  //  number of hardening parameters
  long nh=dfdsdq.n;
  
  matrix dfdsdqt(6,nh);

  switch (ip[ipp].tm[idpm]){
  case simplas1d:{
    spl1d[ip[ipp].idm[idpm]].dfdsigmadq (dfdsdqt);
    break;
  }
  case jflow2:{
    j2f2[ip[ipp].idm[idpm]].deryieldfdsigmadq (dfdsdqt);
    break;
  }
  case modcamclaymat:{
    cclay[ip[ipp].idm[idpm]].deryieldfdsdq (dfdsdqt);