plelemlt.cpp

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

  //  transformation of nodal forces
  //  (in the case of nodal coordinate systems)
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    //globloctransf (nfor,v,tmat);
    glvectortransf (nfor,v,tmat);
    copyv (v,nfor);
  }
}

/**
   function computes nodal forces caused by temperature changes
   
   this function is used in plane stress/strain elements (function is called
   by function res_eigstrain_forces) and shell elements

   @param eid - element id
   @param ri,ci - row and column indices
   @param nfor - array containing nodal forces
   @param x,y - nodal coordinates
   
   30.11.2002, JK
*/
void planeelemlt::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y)
{
  integratedquant iq;
  iq=eigstress;
  
  //  computation of eigenstresses
  compute_eigstress (lcid,eid,ri,ci);
  
  //  integration of stresses over the element
  elem_integration (iq,lcid,eid,ri,ci,nfor,x,y);
}



/**
   function computes internal forces

   this function is used in plane stress/strain elements (function is called
   by function res_internal_forces) and shell elements

   @param lcid - number of load case
   @param eid - element id
   @param ri,ci - row and column indices
   @param ifor - vector of internal forces
   @param x,y - node coordinates

   JK, 27.11.2006
*/
void planeelemlt::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
{
  integratedquant iq;
  iq=locstress;
  
  //  computation of stresses
  compute_nlstress (lcid,eid,ri,ci);
  
  //  integration of stresses over the element
  elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);

}

/**
   function computes internal forces
   
   @param lcid - load case id
   @param eid - element id
   @param ifor - internal forces
   
   JK, 27.11.2006
*/
void planeelemlt::res_internal_forces (long lcid,long eid,vector &ifor)
{
  long transf;
  ivector nodes (nne);
  vector v(ndofe),x(nne),y(nne);
  
  Mt->give_node_coord2d (x,y,eid);
  
  internal_forces (lcid,eid,0,0,ifor,x,y);

  //  transformation of nodal forces
  //  (in the case of nodal coordinate systems)
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    //globloctransf (ifor,v,tmat);
    glvectortransf (ifor,v,tmat);
    copyv (v,ifor);
  }
}

/**
   function computes nonlocal internal forces

   this function is used in plane stress/strain elements (function is called
   by function res_nonloc_internal_forces) and shell elements

   @param lcid - number of load case
   @param eid - element id
   @param ri,ci - row and column indices
   @param ifor - %vector of internal forces
   @param x,y - node coordinates
   
   JK, 27.11.2006
*/
void planeelemlt::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
{
  integratedquant iq;
  iq=nonlocstress;
  
  //  computation of nonlocal stresses
  compute_nonloc_nlstress (lcid,eid,ri,ci);
  
  //  integration of stresses over the element
  elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
}

/**
   function computes nonlocal internal forces

   @param lcid - number of load case
   @param eid - element id
   @param ri,ci - row and column indices
   @param ifor - %vector of internal forces
   
   JK, 27.11.2006
*/
void planeelemlt::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
{
  long transf;
  ivector nodes (nne);
  vector v(ndofe),x(nne),y(nne);
  
  Mt->give_node_coord2d (x,y,eid);

  nonloc_internal_forces (lcid,eid,0,0,ifor,x,y);

  //  transformation of nodal forces
  //  (in the case of nodal coordinate systems)
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    //globloctransf (ifor,v,tmat);
    glvectortransf (ifor,v,tmat);
    copyv (v,ifor);
  }
}








/**
   function returns coordinates of integration points

   @param eid - element id
   @param ipp - integration point pointer
   @param ri - row index
   @param ci - column index
   @param coord - vector of coordinates

   19.1.2002
*/
void planeelemlt::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
{
  long i,ii;
  vector x(nne),y(nne),w(intordsm[ri][ci]),gp1(intordsm[ri][ci]),gp2(intordsm[ri][ci]);
  
  gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ri][ci]);
  Mt->give_node_coord2d (x,y,eid);
  ii=Mt->elements[eid].ipp[ri][ci];
  
  for (i=0;i<intordsm[ri][ci];i++){
    if (ii==ipp){
      coord[0]=approx_nat (gp1[i],gp2[i],x);
      coord[1]=approx_nat (gp1[i],gp2[i],y);
      coord[2]=0.0;
    }
    ii++;
  }

}

/**
   function returns coordinates of integration points

   @param eid - element id
   @param ri - row index
   @param ci - column index
   @param coord - array of coordinates

   19.1.2002
*/
void planeelemlt::ipcoordblock (long eid,long ri,long ci,double **coord)
{
  vector x(nne),y(nne),w(intordsm[ri][ci]),gp1(intordsm[ri][ci]),gp2(intordsm[ri][ci]);
  
  gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ri][ci]);
  Mt->give_node_coord2d (x,y,eid);
  
  for (long i=0;i<intordsm[ri][ci];i++){
    coord[i][0]=approx_nat (gp1[i],gp2[i],x);
    coord[i][1]=approx_nat (gp1[i],gp2[i],y);
    coord[i][2]=0.0;
  }
}



void planeelemlt::nodeforces (long eid,long *le,double *nv,vector &nf)
{
  long i;
  double xi,eta,jac;
  vector x(nne),y(nne),gp(intordb),w(intordb),av(ndofe),v(ndofe);
  matrix n(napfun,ndofe),am(ndofe,ndofe);

  Mt->give_node_coord2d (x,y,eid);
  gauss_points (gp.a,w.a,intordb);

  if (le[0]==1){
    fillm (0.0,am);
    for (i=0;i<intordb;i++){
      xi=(1.0-gp[i])/2.0;  eta=(1.0+gp[i])/2.0;
      bf_matrix (n,xi,eta);

      jac1d_2d (jac,x,y,gp[i],0);
      jac*=w[i];

      nnj (am.a,n.a,jac,n.m,n.n);
    }
    fillv (0.0,av);
    av[0]=nv[0];  av[1]=nv[1];  av[2]=nv[2];  av[3]=nv[3];
    mxv (am,av,v);  addv (nf,v,nf);
  }
  if (le[1]==1){
    fillm (0.0,am);
    for (i=0;i<intordb;i++){
      xi=0.0;  eta=(1.0-gp[i])/2.0;
      bf_matrix (n,xi,eta);

      jac1d_2d (jac,x,y,gp[i],1);
      jac*=w[i];

      nnj (am.a,n.a,jac,n.m,n.n);
    }
    fillv (0.0,av);
    av[2]=nv[4];  av[3]=nv[5];  av[4]=nv[6];  av[5]=nv[7];
    mxv (am,av,v);  addv (nf,v,nf);
  }
  if (le[2]==1){
    fillm (0.0,am);
    for (i=0;i<intordb;i++){
      xi=(1.0+gp[i])/2.0;  eta=0.0;
      bf_matrix (n,xi,eta);

      jac1d_2d (jac,x,y,gp[i],2);
      jac*=w[i];

      nnj (am.a,n.a,jac,n.m,n.n);
    }
    fillv (0.0,av);
    av[4]=nv[8];  av[5]=nv[9];  av[0]=nv[10];  av[1]=nv[11];
    mxv (am,av,v);  addv (nf,v,nf);
  }
}



void planeelemlt::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
{
  long i, j, k, ipp;
  long ii, jj, nv = nodval.n;
  long nstra;
  double xi, eta, ipval;
  vector w, gp1, gp2, anv(nne);

  nstra = 0;
  for (j = 0; j < nv; j++) // for all initial values
  {
    for(i = 0; i < nne; i++) // for all nodes on element
      anv[i] = nodval[i][j];
    for (ii = 0; ii < nb; ii++)
    {
      for (jj = 0; jj < nb; jj++)
      {
        ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
        if (intordsm[ii][jj] == 0)
          continue;
        allocv (intordsm[ii][jj],gp1);
        allocv (intordsm[ii][jj],gp2);
        allocv (intordsm[ii][jj],w);
        gauss_points_tr (gp1.a, gp2.a, w.a, intordsm[ii][jj]);
        for (k = 0; k < intordsm[ii][jj]; k++)
        {
          xi=gp1[k];
          eta=gp2[k];
          //  value in integration point
          ipval = approx_nat (xi,eta,anv);
          if ((ictn[i] & inistrain) && (j < Mm->ip[ipp].ncompstr))
          {
            Mm->ip[ipp].strain[j] += ipval;
            ipp++;
            continue;
          }
          if ((ictn[i] & inistress) && (j < nstra + Mm->ip[ipp].ncompstr))
          {
            Mm->ip[ipp].stress[j] += ipval;
            ipp++;
            continue;
          }
          if ((ictn[i] & iniother) && (j < nstra + Mm->ip[ipp].ncompother))
          {
            Mm->ip[ipp].other[j] += ipval;
            ipp++;
            continue;
          }
          if ((ictn[i] & inicond) && (j < nv))
          {
            if (Mm->ic[ipp] == NULL)
	    {
              Mm->ic[ipp] = new double[nv-j];
              memset(Mm->ic[ipp], 0, sizeof(*Mm->ic)*(nv-j)); 
	    }
            Mm->ip[ipp].other[j] += ipval;
            ipp++;
            continue;
          }
          ipp++;
        }
        destrv(gp1);  destrv (gp2);  destrv (w);
      }
    }
    if (ictn[i] & inistrain) nstra++;
    if (ictn[i] & inistress) nstra++;
    if (ictn[i] & iniother)  nstra++;
  }
}


/**
   function computes volume appropriate to integration point
   
   2.3.2004, JK
*/
void planeelemlt::ipvolume (long eid,long ri,long ci)
{
  long ipp;
  double xi,eta,jac;
  vector x(nne),y(nne);
  
  Mt->give_node_coord2d (x,y,eid);
  ipp=Mt->elements[eid].ipp[ri][ci];
  xi=1.0/3.0;  eta=1.0/3.0;
  
  jac_2d (jac,x,y,xi,eta);
  
  Mm->storeipvol (ipp,jac);
}


////////////////////       /* termitovo */       ////////////////////////////////////

/**
   Function integrates function "NT*D*B*r" (= NT*Stress) over whole element.
   N is matrix of basic functions.
   !!! Values in 'ntdbr' are stored unusually. Not after this manner {(val[0]; ... ;val[tncomp])[0] ; ... ; (val[0]; ... ;val[tncomp])[nne]}
   but after this manner {(val[0]; ... ;val[nne])[0] ; ... ; (val[0]; ... ;val[nne])[ntncomp]}
   
   @param eid   - element id
   @param ntdbr - empty(returned) array, dimension is tncomp*nne
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
 */
void planeelemlt::ntdbr_vector (long eid,vector &ntdbr)
{
  long i,j,ipp,ri,ci;
  double det,thick,jac,w;
  ivector nodes(nne);
  vector x(nne),y(nne),areacoord(3),t(nne),eps(tncomp),sig(tncomp);
  matrix d(tncomp,tncomp);
  
  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord2d (x,y,eid);
  Mc->give_thickness (eid,nodes,t);
  
  fillv (1.0/3.0,areacoord);
  w = 0.5;
  
  det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
  thick = approx (areacoord,t);
  jac = w*thick*det;
  
  Mm->matstiff (d,ipp);
  Mm->givestrain (0,ipp,eps);
  mxv (d,eps,sig);
  
  for (i=0;i<tncomp;i++)
    for (j=0;j<nne;j++)
      ntdbr[j+i*nne] = jac * sig[i] * areacoord[j];
}

/**
   Function integrates function "NT*N" over whole element.
   N is matrix of basic functions.
   !!! Matrix N is composed of three same blocks, it is computed only one third.
   
   @param eid - element id
   @param ntn - empty(returned) 2Darray, dimension is nne x nne
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
 */
void planeelemlt::ntn_matrix (long eid,matrix &ntn)
{
  long i;
  double det,thick,jac;
  ivector nodes(nne);
  vector x(nne),y(nne),areacoord(3),t(nne);
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord2d (x,y,eid);
  Mc->give_thickness (eid,nodes,t);
  
  fillv (1.0/3.0,areacoord);

  det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
  thick = approx (areacoord,t);
  jac = thick*det/24.0;

  fillm (jac,ntn);

  for (i=0;i<nne;i++)
    ntn[i][i] *= 2;
}

/**
   prvne vykosti ty flagy a pak tam teprv soupni popis 
 */
void planeelemlt::compute_error (long eid,double &volume,double &e2,double &u2,double &sizel,double *rderfull,long flags)
{
  long intord = 3;
  long i,ipp,ri,ci;
  double thick,det,jac,contr;
  double zero=1.0e-20;
  ivector nodes(nne);
  vector x(nne),y(nne),areacoord(3),t(nne),gp1(intord),gp2(intord),w(intord);
  vector der1(tncomp),der2(tncomp),der_star(tncomp),der_err(tncomp),der_err2(tncomp);
  matrix d,dinv;
  
  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];

  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord2d (x,y,eid);
  Mc->give_thickness (eid,nodes,t);
  
  det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
  
  if ( !(flags & 8) || ((flags & 8) && !(flags & 4)) ){
    allocm (tncomp,tncomp,d);
    Mm->matstiff (d,ipp);
  }
  if ( !(flags & 8) && (flags & 2) ){
    allocm (tncomp,tncomp,dinv);
    invm (d,dinv,zero);
  }
  
  // compute u2
  fillv (1.0/3.0,areacoord);
  
  thick = approx (areacoord,t);
  jac = 0.5*thick*det;
  
  if (!(flags & 8)){         // linear
    if (!(flags & 2)){          // strain - der1 == strain
      Mm->givestrain (0,ipp,der1);
      mxv (d,der1,der2);
    }
    else                        // stress - der1 == stress
      if (!(flags & 4)){          // own
	Mm->givestrain (0,ipp,der2);
	mxv (d,der2,der1);
      }
      else{                       // other
	Mm->givestress (0,ipp,der1);
	mxv (dinv,der1,der2);
      }
  }
  else{                      // nonlinear - der1 == strain
    Mm->givestrain (0,ipp,der1);  // strain
    if (flags & 4)              // stress own
      Mm->givestress (0,ipp,der2);
    else                        // stress other
      mxv (d,der1,der2);
  }
  
  scprd (der1,der2,contr);
  u2 = jac*contr;
  
  // compute e2
  gauss_points_tr (gp1.a,gp2.a,w.a,intord);
  
  e2 = 0;
  for (i=0;i<intord;i++){
    areacoord[0]=gp1[i];
    areacoord[1]=gp2[i];
    areacoord[2]=1.0-gp1[i]-gp2[i];
    
    thick = approx (areacoord,t);
    jac=w[i]*thick*det;
    
    if (!(flags & 8)){         // linear
      give_der_star (areacoord,rderfull,nodes,der_star,Mt->nn);
      subv (der_star,der1,der_err);
      
      if (!(flags & 2))          // strain - der1 == strain
	vxmxv (der_err,d,contr);
      else                       // stress - der1 == stress
	vxmxv (der_err,dinv,contr);
      
    }
    else{                      // nonlinear
      give_der_star (areacoord,rderfull                ,nodes,der_star,Mt->nn);
      subv (der_star,der1,der_err);
      give_der_star (areacoord,rderfull + tncomp*Mt->nn,nodes,der_star,Mt->nn);
      subv (der_star,der2,der_err2);
      
      scprd (der_err,der_err2,contr);
    }
    
    e2 += jac*contr;
  }
  
  volume += det/2.0;
  sizel = sqrt(1.1547005384*det);
}

/**
   This function prepare array 'spder' for using in function `least_square::spr_default`.
   It allocates and fills up array 'spder' by derivatives(strain ro stress) in sampling points.
   For planeelemlt sampling points == main integration point == Gauss's point at the centroid of element.
   
   @param eid   - element id
   @param spder - allocated array 
   @param flags - determines by which derivate is spder filled
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
 */
void planeelemlt::elchar (long eid,double *&spder,long flags)
{
  long ipp,ri,ci;
  vector eps;
  matrix d;
  
  spder = new double[tncomp*1];
  
  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  
  
  if (!(flags & 2)){
    eps.n = tncomp;
    eps.a = spder;
    Mm->givestrain (0,ipp,eps);
    eps.a = NULL;
  }
  else if (!(flags & 4)){
    allocv (tncomp,eps);
    allocm (tncomp,tncomp,d);
    Mm->matstiff (d,ipp);
    Mm->givestrain (0,ipp,eps);
    mxv (d.a,eps.a,spder,d.m,d.n);
  }
  else{
    eps.n = tncomp;
    eps.a = spder;
    Mm->givestress (0,ipp,eps);
    eps.a = NULL;
  }
}

/**
   created  3.5.2002, Ladislav Svoboda, [e-mail]
 */
double planeelemlt::error (long eid,vector &n,double &a)
{
  long i;
  double answer;
  ivector nodes(nne);
  vector x(nne),y(nne),areacoord(3),nodval(nne);
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord2d (x,y,eid);
  
  fillv (1.0/3.0,areacoord);
  
  a = 0.5 * (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
  
  for (i=0;i<nne;i++)
    nodval[i] = n[nodes[i]];
  
  answer = approx (areacoord,nodval);
  answer = a*answer*answer;
  
  return answer;
}
////////////////////       /* termitovo */       ////////////////////////////////////