nonlocplast.cpp

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

#include <string.h>

#include "nonlocplast.h"
#include "global.h"
#include "mechmat.h"
#include "intpoints.h"
#include "matrix.h"
#include "vector.h"
#include "vecttens.h"

/**
  This constructor inializes attributes to zero values.
*/
nonlocplast::nonlocplast(void)
{
  r = 0;
  waf = wavrg(0);
}

/**
  This destructor is only for the formal purposes.
*/
nonlocplast::~nonlocplast(void)
{
  r = 0;
  waf = wavrg(0);
}

/**
  This function reads material parameters from the opened text file given
  by the parameter in.

  @param in - pointer to the opned text file

*/
long nonlocplast::read(XFILE *in)
{
  if (xfscanf(in, "%le %d", &r, (int*)&waf) == 2)
    return 0;

  fprintf(stderr, "\n Error reading parameters of nonlocal plasticity\n");
  fprintf(stderr, " in function nonlocplast::read, (file %s, line %d)\n", __FILE__, __LINE__);
  return 1;

}

/**
  This function allocates additional data arrays on the given integration point for the computed
  averaged internal parameters.

  @param ipp - integration point number

*/
void nonlocplast::allocation (long ipp)
{
  switch (waf){
  case avggamma:{
    Mm->ip[ipp].nonloc = new double [Mm->ip[ipp].ncompother];
    memset (Mm->ip[ipp].nonloc, 0, sizeof(*Mm->ip[ipp].nonloc)*Mm->ip[ipp].ncompother);
    break;
  }
  case avgepsp:{
    Mm->ip[ipp].nonloc = new double [Mm->ip[ipp].ncompstr];
    memset (Mm->ip[ipp].nonloc, 0, sizeof(*Mm->ip[ipp].nonloc)*Mm->ip[ipp].ncompstr);
    break;
  }
  default:{
    fprintf(stderr, "\n Unknown type of the averageing variable is required.\n");
    fprintf(stderr, " in function nonlocplast::allocation, (file %s, line %d)\n", __FILE__, __LINE__);
  }
  }
}


/**
  This function averages values of the gamma or plastic strains in the given integration point.
  Which values are averaged is defined by the data member waf.

  @param ipp - integration point number

*/
void nonlocplast::average (long ipp)
{
  double rr, alpha, sum = 0.0;
  long ncompo, ncomps = Mm->ip[ipp].ncompstr;
  matrix d(ncomps,ncomps), sigt(3,3), dgdst(3,3);
  double gamma, agamma = 0.0;
  vector dgds(ncomps),aepsp(ncomps), epsa(ncomps), epsn(ncomps), epsp(ncomps), sig(ncomps);
  vector q(1);
  long i, j, nad, aip;

  for (i=0;i<ncomps;i++){
    Mm->ip[ipp].nonloc[i]=0.0;
  }

  //nad .... number of adjacent points
  nad=Mt->nadjip[ipp];

  //loop over all adjacent points
  for (i = 0; i < nad; i++)
  {
    //number of adj. int. point
    aip=Mt->adjip[ipp][i];
    //compute distance of int. points ipp<->aip
    rr=Mt->dist[ipp][i];
    //compute nonlocal weight
    alpha = (rr >= r) ? 0 : (1-rr*rr/(r*r));
    // multiply by the ip volume
    alpha *= Mm->ipv[aip];
    sum += alpha;

    switch (waf)
    {
      case avggamma :
        ncomps = Mm->ip[aip].ncompstr;
        gamma = Mm->ip[aip].other[ncomps];
        agamma += alpha * gamma;
        break;
      case avgepsp :
        for (j = 0; j < ncomps; j++)
          Mm->ip[ipp].nonloc[j] += alpha*Mm->ip[aip].other[j];
        break;
      default :
        fprintf(stderr, "\n Unknown type of the averageing variable is required.\n");
        fprintf(stderr, " in function nonlocplast::nlstresses, (file %s, line %d)\n", __FILE__, __LINE__);
        break;
    }
  }//loop i
  ncomps = Mm->ip[ipp].ncompstr;
  ncompo = Mm->ip[ipp].ncompother;
  switch (waf)
  {
    case avggamma :
    {
      //  initial values
      for (i = 0; i < ncomps; i++)
      {
        epsn[i]=Mm->ip[ipp].strain[i];
        epsp[i]=Mm->ip[ipp].eqother[i];
      }
      //  elastic strain
      subv (epsn,epsp,epsa);
      //  elastic stiffness matrix
      Mm->elmatstiff (d,ipp);
      //  stress computation
      mxv (d,epsa,sig);
      // values of averaged plastic strains are computed
      vector_tensor (sig, sigt, Mm->ip[ipp].ssst, stress);
      // obtaining hardening parameters
      q[0] = 0.0;
      if (ncompo-ncomps-1 > 0)
        q[0] = Mm->ip[ipp].nonloc[ncomps+1];
      Mm->dgdsigma (ipp, 0, sigt, q, dgds);
      vector_tensor (dgds,dgdst,Mm->ip[ipp].ssst, strain);
      cmulm((agamma-Mm->ip[ipp].eqother[ncomps])/sum, dgdst);
      tensor_vector(aepsp, dgdst, Mm->ip[ipp].ssst, strain);
      // storing averaged values
      for (i = 0; i < ncomps; i++)
        Mm->ip[ipp].nonloc[i] = aepsp[i];
      Mm->ip[ipp].nonloc[ncomps] = agamma/sum;
      break;
    }
    case avgepsp :
      // storing averaged values
      for (i = 0; i < ncomps; i++)
        Mm->ip[ipp].nonloc[i] /= sum;
      break;
    default :
      fprintf(stderr, "\n Unknown type of the averageing variable is required.\n");
      fprintf(stderr, " in function nonlocplast::nlstresses, (file %s, line %d)\n", __FILE__, __LINE__);
      break;

  }
}