epsolver.cpp

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

#include "epsolver.h"
#include "global.h"
#include "globmat.h"
#include "mechprint.h"
#include "loadcase.h"
#include "dloadcase.h"
#include "nssolver.h"
#include "springel.h"
#include "intpoints.h"
#include "auxdatout.h"
#include "vector.h"
#include "matrix.h"
#include "mathem.h"
#include "gmatrix.h"
#include "tensor.h"
#include "vecttens.h"

#include <math.h>
#include <string.h>

/**
  This function controls run of earth pressure problem computation.
*/
void solve_epressure ()
{
  long i;
  double *rhs, *lhs;


  for (i=0;i<Lsrs->nlc;i++){

    rhs=Lsrs->give_rhs (i);
    lhs=Lsrs->give_lhs (i);
    //  load vectors
    nullv (rhs+i*2*Ndofm, 2*Ndofm);
    if (Mp->tepsol < epplast_sol)
    {
      Mb->lc[i].assemble (i,rhs+i*Ndofm,1.0);
      Mb->dlc[i].assemble (rhs+i*Ndofm, lhs+i*Ndofm);
    }
    else
    {
      Mb->lc[i*2+0].assemble (i*2+0,rhs+(i*2+0)*Ndofm,1.0);
      Mb->lc[i*2+1].assemble (i*2+1,rhs+(i*2+1)*Ndofm,1.0);
    }
    //  solver
    epressure_solver (i);
  }
}

/**
  This function calls appropriate solver for the earth pressure problem

  @param lcid - load case id.
*/
void epressure_solver (long lcid)
{
  switch (Mp->tepsol){
  case gep_sol:{
    general_epressure (lcid);
    break;
  }
  case gepvarsup_sol:{
    general_epressure_varsup (lcid);
    break;
  }
  case epplast_sol:{
    femplast_epressure (lcid);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown solver of nonlinear equation system is required");
    fprintf (stderr,"\n in function nonlinear_solver (solver.cpp).\n");
  }
  }

}

/**
  This function solves erath pressure problem using GLPT. Iteration of load is used.

  @param lcid - load case id.
*/
void general_epressure (long lcid)
{
  long i, j, n, ni;
  double zero, err, norf, dl;
  double *r,*f,*fi,*fb, *fc, lambda, rlam;


  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->niep;
  //  computer zero
  zero = Mp->zero;
  // required norm of unbalanced forces
  err = Mp->errep;
  // required dynamic load step
  dl = Mp->stepep;

  fb = new double [n];
  fi = new double [n];
  memset (fb,0,n*sizeof(double));
  memset (fi,0,n*sizeof(double));

  //  initialization phase
  r  = Lsrs->give_lhs (lcid);
  f  = Lsrs->give_rhs (2*lcid);
  fc = Lsrs->give_rhs (2*lcid+1);

  //***************************
  //  main iteration loop  ****
  //***************************
  lambda = 0;
  for (i=0;i<ni;i++){

    print_init(i, "wt");
    //  vectors of loading
    nullv (f+lcid*n, n);
    Mb->lc [lcid].assemble (0,f+2*lcid*n,1.0);
    Mb->dlc[lcid].assemble (f+2*lcid*n, r+lcid*n);
    //  backup of load vector
    for (j=0;j<n;j++){
      f[j] = (f[j] - fb[j]); // difference between new and backupped load vector 
      if (lambda < 1.0) // in case that the load is not fully applied
        f[j] *= dl/(1.0-lambda);     // use only required increment
      fc[j] = fb[j]; // store previously reached load
      fb[j] += f[j]; // add new increment to previously reached load
    }
    lambda += dl;
    if (lambda > 1.0)
      lambda = 1.0;
    //    memset (r,0,n*sizeof(double));
    rlam = arclengthrv(lcid, 1.0, Mp->rlerrep);    
    if (fabs(rlam-1.0) > Mp->rlerrep)
    {
      fprintf(stderr, "\n\n Arclength could not reach required lambda\n");
      fprintf(stderr, "(file %s, line %d)\n", __FILE__, __LINE__);
      abort();
    }
    //  computation of internal forces, updating strains in the integration points
//    internal_forces (lcid,fi);
    print_step(lcid, i+1, lambda, fb);

    nullv (f+lcid*n, n);
    Mb->lc [lcid].assemble (0,f+2*lcid*n,1.0);
    Mb->dlc[lcid].assemble (f+2*lcid*n, r+lcid*n);
    for (j = 0; j < n; j++)
      fi[j] = f[j] - fb[j];
    norf = ss(fi, fi, n);
    norf /= ss(f, f, n);
    fprintf (stdout, "\n\n MAIN ITERATION LOOP : NORF=%e, lambda = %e\n", norf, lambda);
    fprintf (stdout, "********************************************\n\n");
    if ((norf < err) && (lambda == 1.0))
      break;
    print_close();
  }
  print_close();
  delete [] fi;  delete [] fb;
}


/**
  This function solves erath pressure problem using GLPT. Temporary supports are used for
  better algorithm convergence.

  @param lcid - load case id.
*/
void general_epressure_varsup (long lcid)
{
  long i, j, n, ni;
  double zero, err, norf;
  double *r,*rb,*f,*fi,*fb, lambda;
  double **br;


  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->niep;
  //  computer zero
  zero = Mp->zero;
  // required norm of unbalanced forces
  err = Mp->errep;

  rb = new double [n];
  fb = new double [n];
  fi = new double [n];
  br = new double*[Mt->nn];
  memset (rb,0,n*sizeof(double));
  memset (fb,0,n*sizeof(double));
  memset (fi,0,n*sizeof(double));
  memset (br,0,Mt->nn*sizeof(*br));
  for (i = 0; i < Mt->nn; i++)
    br[i] = new double[Gtm->give_ndofn(i)];

  //  initialization phase
  r  = Lsrs->give_lhs (lcid);
  f  = Lsrs->give_rhs (lcid);

  //***************************
  //  main iteration loop  ****
  //***************************
  for (i=0;i<ni;i++) 
  {
    print_init(i, "wt");
    //  vectors of loading
    nullv (f+lcid*n, n);
    Mb->lc [lcid].assemble (0,f+lcid*n,1.0);
    Mb->dlc[lcid].assemble (f+lcid*n, r+lcid*n);
    //  backup of load vector
    for (j=0;j<n;j++){
      fb[j] = f[j];
    }

    memset (r,0,n*sizeof(double));
    lambda = arclengthrv(lcid, 1.0, Mp->rlerrep);

    //  computation of internal forces, updating strains in the integration points
    internal_forces (lcid,fi);
    for (j = 0; j < n; j++)
      f[j] = fb[j] - fi[j];
    print_step(lcid, i+1, lambda, fb);
    norf = ss(f, f, n);
    fprintf (stdout, "\n MAIN ITERATION LOOP : NORF=%e", norf);
    if ((norf < err) && (i >= Mp->nselnodep))
      break;
    if (i < Mp->nselnodep)
    {
      // restore original support state
      Gtm->restore_codnum();
      fprintf(stdout, "\n removing support in the node number %ld", Mp->selnodep[i]+1);
      for (j = 0; j <= i; j++)
      {
        // remove support in the x direction from the selected nodes depending current iteration step
        Gtm->gnodes[Mp->selnodep[j]].cn[0] = 1;
      }
      // generation new code numbers
      n = Ndofm = Gtm->gencodnum();
      // backup of displacements in the nodes
      for (j = 0; j < Mt->nn; j++)
        noddispl(lcid, br[j], j);
      //  initialization phase for next step
      //  left and right hand sides should be deleted
      delete Lsrs;
      delete Smat;  // size of global matrix will grow -> we can delete it
      Smat = NULL;
      delete [] fi;
      delete [] fb;
      delete [] rb;
      Lsrs = new lhsrhs;
      Lsrs->alloc ();
      rb = new double [n];
      fb = new double [n];
      fi = new double [n];
      memset (rb,0,n*sizeof(double));
      memset (fb,0,n*sizeof(double));
      memset (fi,0,n*sizeof(double));

      r  = Lsrs->give_lhs (lcid);
      f  = Lsrs->give_rhs (lcid);
      restore_displ(lcid, br);
    }
    print_close();
  }
  print_close();
  delete [] fi;  delete [] fb; delete [] rb;
  for (i = 0; i < Gtm->nn; i++)
    delete [] br;
  delete [] br;
}

/**
  This function solves erath pressure problem using FEM plasticity models. 
  Iteration of rigid spring support is used. Spring supports are incrementally 
  removed during each itreation step. Removed supports are expected to be ordered
  as last elements, each step is removed the last element and number of elements is 
  decremented. So the user should order removed spring elements so the first removed
  spring element is the last element (element with the greatest number) and so on.
  Newton-Raphson method is used for given plasticity problem solution. This solver supposes
  the static load equal to dead weight of excaved soil is applied to the bottom of 
  excaved part of the soil body and this static load is decremented in the each step when 
  the spring support is removed.

  @param lcid - load case id.
*/
void femplast_epressure (long lcid)
{
  long i,j,idxp, idxs,ni, *nod, ndofn, *cnum;
  double *rhs, *rb, *react, *fpb, *fsb, *forig, lambda;
  vector ifor;
  ivector cn;

  vector sig;
  matrix sigt(3,3);
  double *p = new double[Mm->tnip];
  double *porig = new double[Mm->tnip];
  char fn[1001];
  char fn2[1001];
  FILE *femcad;

  ni = Mp->niep+1;
  //***************************
  //  main iteration loop  ****
  //***************************
  rb    = new double[Ndofm];
  fpb   = new double[Ndofm];
  fsb   = new double[Ndofm];
  forig = new double[Ndofm];
  nod   = new long  [Mp->niep];
  react = new double[Mp->niep];
  cnum  = new long  [Mp->niep]; 
  for (i=0;i<ni;i++)
  {
    print_init(i+1, "wt");
    if (i < 1)
      fprintf (stdout, "\n MAIN ITERATION LOOP : step %ld\n", i+1);
    else
      fprintf (stdout, "\n MAIN ITERATION LOOP : step %ld, elem. %ld was removed\n", i+1, Mt->ne+1);
    fprintf(stdout, " ---------------------------------------\n\n");
    rhs=Lsrs->give_rhs (lcid);
    idxs = (lcid*2+1)*Ndofm;
    idxp = (lcid*2+0)*Ndofm;
    if (i == 0)
    {
      copyv (rhs+idxp, fpb, Ndofm);
      copyv (rhs+idxs, fsb, Ndofm);
      addv(rhs+idxp, fsb, Ndofm);
      nullv (rhs+idxs, Ndofm);
      nullv (forig, Ndofm);
    }
    if (i > 0)
    {      
      nullv (rhs+idxp, Ndofm);
      for (j = 0; j < Ndofm; j++)
      {
        if (i == 1)
          rhs[idxp+j] = (Mp->loadcoef[i-1]-1.0)*fsb[j];
        else
          rhs[idxp+j] = (Mp->loadcoef[i-1]-Mp->loadcoef[i-2])*fsb[j];
      }
/*      if (cnum[i-1] > 0)
	rhs[idxp+cnum[i-1]-1] += react[i-1];*/
    }
    lambda = newton_raphsonep(lcid, forig);
//    nonlinear_solver (lcid);
    for (j=0; j < Ndofm; j++)   
      forig[j] += lambda * (rhs+idxp)[j];
    
    Mm->computenlstresses(0);
    for (j=0; j<Mm->tnip; j++)
    {
      allocv(Mm->ip[j].ncompstr, sig);
      Mm->givestress (0,j,sig);
      vector_tensor (sig, sigt, Mm->ip[j].ssst, stress);
      p[j] = first_invar (sigt)/3;
      destrv(sig);
    }
    if (i == 0)
      copyv(p, porig,Ndofm);
    else
    {
      for (j=0; j<Mm->tnip; j++)
      {
        p[j] -= porig[j];      
      }
      sprintf(fn,"%s%s.jk0", Mp->path, Mp->filename);
      sprintf(fn2,"%s%s.jk0.%ld", Mp->path, Mp->filename,i);
      femcad = fopen(fn, "at");
      char dlcid[50];
      sprintf(dlcid, "%ld", i+1); 
      write_elemscalar(femcad, p, "p", dlcid);
      fclose(femcad);
      rename (fn, fn2);
    }
	     

    if (i == Mp->niep)
      break;
    nod[i] = Gtm->gelements[Gtm->ne-1].nodes[0];
    ndofn = Mt->give_ndofn(nod[i]);
    allocv(ndofn, ifor);
    allocv(ndofn, cn);
    Spring->res_internal_forces (lcid, Mt->ne-1, ifor);
    react[i] = ifor[0];
    Mt->give_node_code_numbers (nod[i],cn.a);
    cnum[i] = cn[0];
    // removing currently last spring support element.
    Mt->ne--;
    Gtm->ne--;
    destrv(ifor);
    destrv(cn);
    print_close();
  }
  print_close();
  delete [] rb;
  delete [] fpb;
  delete [] fsb;
  delete [] nod;
  delete [] react;
  delete [] cnum;
}




/**
  This function extracts nodal displacements bckr to the displacement %vector of Lsrs.
  @param lcid  - number of load case
  @param bckr  - array with displacement in the each node
  //  3.12.2002
*/
void restore_displ(long lcid, double **bckr)
{
  long i, j, ii;

  for (i = 0; i < Gtm->nn; i++)
  {
    for (j = 0; j < Gtm->gnodes[i].ndofn; j++)
    {
      ii=Gtm->gnodes[i].cn[j];
      if (ii >  0)   Lsrs->lhs[lcid*Ndofm+ii-1] = bckr[i][j];
    }
  }
}

/**
  This function solves system of nonlinear algebraic equations
  by the arc-length method. This method was modified to stop iteration
  on the given value rlambda of the lambda parameter. Rlambda value is reached with
  required error rlerr. Only one right hand side vector is supported with respect
  to nonlinearity and absence of superposition.


  @param lcid    - load case id
  @param rlambda - required value of lambda parameter
  @param rlerr   - required error between lambda and rlambda

  @retval The function returns reached lambda parameter.

*/
double arclengthrv (long lcid, double rlambda, double rlerr)
{
  long i,j,k,n,ni,ini,stop,modif,li,newmesh, numr;
  double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr;
  double lambdao,ss1,ss2,ss3,ss4,ss5;
  double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
  long cn, ecn[6];
  char file[255];
  char *mode;
//  FILE *gr;
  FILE *grout;
  
  
  
  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->nlman.nial;