scaldamcc.cpp

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

#include "scaldamcc.h"
#include "matrix.h"
#include "vector.h"
#include "tensor.h"
#include "elastisomat.h"
#include "global.h"
#include "intpoints.h"
#include "vecttens.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <iostream.h>

#define nijac 20
#define limit 1.0e-8


/**
  This constructor inializes attributes to zero values.
*/
scaldamcc::scaldamcc (void)
{
  k0 = at = bt = ac = bc = beta = k = 0.0;
}

/**
  This destructor is only for the formal purposes.
*/
scaldamcc::~scaldamcc (void)
{

}

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

  @param in - pointer to the opned text file

*/
void scaldamcc::read (XFILE *in)
{
  xfscanf (in,"%lf %lf %lf %lf %lf %lf %d",&k0,&at, &bt, &ac, &bc, &beta, (int *)(&ft));
  if (ft == vonmises)
    xfscanf(in, "%lf", &k);
}

/**
  This function computes parameters for the damage function
  Different types of the parameter computing are given by the
  attribute ft.

  @param  ipp - integration point number
  @param  eps - %vector of the strains
  @param  kappa - %vector where the computed parameters are stored

*/
void scaldamcc::damfuncpar(long ipp, vector &eps, vector &kappa)
{
  vector poseps;
  vector princeps;
  matrix pvect;
  matrix epst;
  matrix dev;
  double i1e, j2e, nu, e, tmp;
  long ncomp, idem;

  switch (ft)
  {
    case vonmises :
      ncomp=Mm->ip[ipp].ncompstr;
      allocm(3, 3, epst);
      allocm(3, 3, dev);
      vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
      i1e = first_invar(epst);
      deviator (epst,dev);
      j2e = second_invar (dev);
      idem = Mm->ip[ipp].gemid();
      if (Mm->ip[ipp].tm[idem] != elisomat)
      {
        fprintf(stderr, "\n\n Invalid type of elastic material is required");
        fprintf(stderr, "\n  in function scaldamcc::damfuncpar, (%s, line %d)\n", __FILE__, __LINE__);
      }
      e  = Mm->eliso[idem].e;
      nu = Mm->eliso[idem].nu;
      kappa[0]  = (k-1.0)/(2.0*k)/(1.0-2.0*nu)*i1e;
      tmp = (k-1.0)*(k-1.0)/(1.0-2.0*nu)/(1.0-2.0*nu)*i1e*i1e -12.0*k/(1.0+nu)/(1.0+nu)*j2e;
      kappa[0] += 1.0/(2.0*k)*sqrt(tmp);
      break;
    case normazar :
      allocm(3, 3, epst);
      allocm(3, 3, pvect);
      allocv(3, princeps);
      allocv(3,poseps);
      fillv(0.0, princeps);
      vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
      princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3);
      extractposv (princeps, poseps);
      scprd(poseps, poseps, tmp);
      kappa[0] = sqrt(tmp);
      break;
    default :
    {
      fprintf(stderr, "\n\n Unknown damage parameter function type");
      fprintf(stderr, "\n  in function scaldamcc::damfuncpar, (%s, line %d)\n", __FILE__, __LINE__);
    }
  }
  return;
}

/**
  This function computes damage parameter omega which is the result of the
  damage function.

  @param ipp - integration point number
  @param kappa - %vector of the parameters of damage function
  @param eps   - strain %vector
*/
double scaldamcc::damfunction(long ipp, vector &kappa, vector &eps)
{
  long ncomp=Mm->ip[ipp].ncompstr;
  matrix epst, sigt, pvect, d(ncomp, ncomp),tmp3(ncomp, ncomp),tmp;
  vector princeps;
  vector poseps;
  vector tenseps;
  vector compeps;
  vector sigm(ncomp), princsig, tmp1, tmp2;
  vector negeps;
  double omega = 0.0, dt, dc, alphat, alphac;

  if (kappa[0] < k0)
  // elastic loading
    return 0.0;

  allocm(3, 3, epst);
  allocm(3, 3, sigt);
  allocm(3, 3, tmp);
  allocv(3, poseps);
  allocv(3, negeps);
  allocv(3, tenseps);
  allocv(3, compeps);
  allocv(3, princeps);
  allocv(3, princsig);
  allocv(3, tmp1);
  allocv(3, tmp2);
  allocm(3, 3, pvect);

  dt = 1.0 - k0 * (1.0 - at) / kappa[0] - at/exp(bt*(kappa[0]-k0));
  dc = 1.0 - k0 * (1.0 - ac) / kappa[0] - ac/exp(bc*(kappa[0]-k0));
  vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
  princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3);



  Mm->elmatstiff (d,ipp);
  mxv(d,eps,sigm);
  vector_tensor(sigm,sigt, Mm->ip[ipp].ssst, stress);
  princ_val(sigt,princsig,pvect,nijac,limit,Mp->zero,3);
  extractposv(princsig,tmp1);
  extractnegv(princsig,tmp2);
  invm(d,tmp3,Mp->zero);

  extractm(tmp,tmp3,0,3);
  mxv(tmp,tmp1,tenseps);
  mxv(tmp,tmp2,compeps);

  extractposv (princeps, poseps);
  subv(princeps, poseps, negeps);
  scprd(tenseps, poseps, alphat);
  alphat = alphat / kappa[0] / kappa[0];
  alphat = pow(alphat, beta);
  scprd(compeps, poseps, alphac);
  alphac = alphac / kappa[0] / kappa[0];
  alphac = pow(alphac, beta);

  omega = alphat*dt + alphac*dc;

  if ((omega > 1.0) || (omega < 0.0)) {
    fprintf(stderr, "\n\nError - iteration of omega doesn't convergate to range <0,1> %g ",omega);
    fprintf(stderr, "\n in function damfunction(), (%s, line %d)\n", __FILE__, __LINE__);
  }
  if(omega>1.0)
    omega = 1.0 ;
  if(omega<0.0)
    omega=0.0;

  return omega;
}

/**
  This function computes material stiffnes %matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number
  @param ido - index of internal variables for given material in the ipp other array

*/
void scaldamcc::matstiff (matrix &d,long ipp,long ido)
{
  double dp;

  switch (Mp->stmat)
  {
    case initial_stiff :
      Mm->elmatstiff (d,ipp);
      break;
    case tangent_stiff :
      Mm->elmatstiff (d,ipp);
      dp=Mm->ip[ipp].other[ido+1];
      cmulm (1.0-dp,d);
      break;
    default :
      fprintf(stderr, "\n\nError - unknown type of stifness matrix");
      fprintf(stderr, "\n in function scaldamcc::matstiff (%s, line %d)\n", __FILE__, __LINE__);
  }
}


/**
  This function computes correct stresses in the integration point and stores
  them into ip stress array.

  @param ipp - integration point pointer
  @param im  - index of material type for given ip
  @param ido - index of internal variables for given material in the ipp other array
  7.10.2001
*/
void scaldamcc::nlstresses (long ipp,long im,long ido)
{
  long i,ncomp=Mm->ip[ipp].ncompstr;
  double dp;
  vector epsn(ncomp),sigma(ncomp), kappa(1);

  //  initial values
  for (i=0;i<ncomp;i++){
    epsn[i] = Mm->ip[ipp].strain[i];
  }
  kappa[0] = Mm->ip[ipp].other[ido+0];

  // damage stress solver
  dp = Mm->scal_dam_sol (ipp, im, ido, epsn, kappa, sigma);

  //  new data storage
  for (i=0;i<ncomp;i++){
    Mm->ip[ipp].stress[i]=sigma[i];
  }
  Mm->ip[ipp].other[ido+2]=kappa[0];
  Mm->ip[ipp].other[ido+3]=dp;
}

void scaldamcc::updateval (long ipp,long im,long ido)
{
  long i,n = Mm->givencompeqother(ipp, im);

  for (i=0;i<n;i++){
    Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
  }

}

/**
  This function returns the value of the limit elastic strain .

*/
double scaldamcc::epsefunction ()
{
  return k0 ;
}