plelemlq.cpp

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

  ivector nodes (nne);
  vector v(ndofe),x(nne),y(nne);
  
  Mt->give_node_coord2d (x,y,eid);
  
  gl_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);
  }
}




void planeelemlq::nodeforces (long eid,long *le,double *nv,vector &nf)
{
  long i;
  double ww,jac,xi,eta;
  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);
    eta = 1.0;
    for (i=0;i<intordb;i++){
      xi = gp[i];
      ww = w[i];
      
      bf_matrix (n,xi,eta);
      
      jac1d_2d (jac,x,y,xi,0);
      jac *= ww;

      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);
    xi = -1.0;
    for (i=0;i<intordb;i++){
      eta = gp[i];
      ww = w[i];

      bf_matrix (n,xi,eta);
      
      jac1d_2d (jac,x,y,eta,1);
      jac *= ww;
      
      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);
    eta = -1.0;
    for (i=0;i<intordb;i++){
      xi = gp[i];
      ww = w[i];
      
      bf_matrix (n,xi,eta);
      
      jac1d_2d (jac,x,y,xi,2);
      jac *= ww;
      
      nnj (am.a,n.a,jac,n.m,n.n);
    }
    fillv (0.0,av);
    av[4]=nv[8];  av[5]=nv[9];  av[6]=nv[10];  av[7]=nv[11];
    mxv (am,av,v);  addv (nf,v,nf);
  }
  if (le[3]==1){
    fillm (0.0,am);
    xi = 1.0;
    for (i=0;i<intordb;i++){
      eta = gp[i];
      ww = w[i];
      
      bf_matrix (n,xi,eta);
      
      jac1d_2d (jac,x,y,eta,3);
      jac *= ww;

      nnj (am.a,n.a,jac,n.m,n.n);
    }
    fillv (0.0,av);
    av[6]=nv[12];  av[7]=nv[13];  av[0]=nv[14];  av[1]=nv[15];
    mxv (am,av,v);  addv (nf,v,nf);
  }
}



void planeelemlq::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
{
  long i, j, k, l, ipp;
  long ii, jj, nv = nodval.n;
  long nstra, ncompstr, ncompeqother;
  double xi, eta, ipval;
  long idstra, idstre, idoth, idic;
  vector w, gp, anv(nne);

  nstra = idstra = idstre = idoth = idic = 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],gp);
        allocv (intordsm[ii][jj],w);
        gauss_points (gp.a,w.a,intordsm[ii][jj]);
        for (k = 0; k < intordsm[ii][jj]; k++)
        {
          xi=gp[k];
          for (l = 0; l < intordsm[ii][jj]; l++)
          {
            eta=gp[l];
            //  value in integration point
            ipval = approx (xi,eta,anv);
            ncompstr =  Mm->ip[ipp].ncompstr;
            ncompeqother = Mm->ip[ipp].ncompeqother;
            if ((ictn[0] & inistrain) && (j < ncompstr))
            {
              Mm->ip[ipp].strain[idstra] += ipval;
              ipp++;
              continue;
            }
            if ((ictn[0] & inistress) && (j < nstra + ncompstr))
            {
              Mm->ip[ipp].stress[idstre] += ipval;
              ipp++;
              continue;
            }
            if ((ictn[0] & iniother) && (j < nstra + ncompeqother))
            {
              Mm->ip[ipp].eqother[idoth] += ipval;
              ipp++;
              continue;
            }
            if ((ictn[0] & 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->ic[ipp][idic] += ipval;
              ipp++;
              continue;
            }
            ipp++;

          }
        }
        destrv (gp);  destrv (w);
      }
    }
    ipp=Mt->elements[eid].ipp[ri][ci];
    ncompstr =  Mm->ip[ipp].ncompstr;
    ncompeqother = Mm->ip[ipp].ncompeqother;
    if ((ictn[0] & inistrain) && (j < ncompstr))
    {
      nstra++;
      idstra++;
    }
    if ((ictn[0] & inistress) && (j < nstra + ncompstr))
    {      
      nstra++;
      idstre++;
    }  
    if ((ictn[0] & iniother)  && (j < nstra + ncompeqother))
    {
      nstra++;
      idoth++;
    }
    if ((ictn[0] & inicond) && (j < nv))
      idic++;
  }
}


/**
   function computes volume appropriate to integration point
   
   2.3.2004, JK
*/
void planeelemlq::ipvolume (long eid,long ri,long ci)
{
  long i,j,ii,jj,ipp;
  double xi,eta,jac;
  vector x(nne),y(nne),w,gp;
  
  Mt->give_node_coord2d (x,y,eid);

  for (ii=0;ii<nb;ii++){
    for (jj=0;jj<nb;jj++){
      if (intordsm[ii][jj]==0)  continue;
      
      allocv (intordsm[ii][jj],w);
      allocv (intordsm[ii][jj],gp);

      gauss_points (gp.a,w.a,intordsm[ii][jj]);

      ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];

      for (i=0;i<intordsm[ii][jj];i++){
	xi=gp[i];
	for (j=0;j<intordsm[ii][jj];j++){
	  eta=gp[j];

	  jac_2d (jac,x,y,xi,eta);
	  jac*=w[i]*w[j];

	  Mm->storeipvol (ipp,jac);
	    
	  ipp++;
	}
      }
      destrv (gp);  destrv (w);
    }
  }
}

/**
   function interpolates the nodal values to the integration points on the element
   
   @param eid - element id
   
   21.6.2004, JK
*/
void planeelemlq::intpointval (long eid,vector &nodval,vector &ipval)
{
  long i,j,ii,jj,k;
  double xi,eta;
  vector w,gp;
  
  k=0;
  for (ii=0;ii<nb;ii++){
    for (jj=0;jj<nb;jj++){
      if (intordsm[ii][jj]==0)  continue;
      
      allocv (intordsm[ii][jj],w);
      allocv (intordsm[ii][jj],gp);

      gauss_points (gp.a,w.a,intordsm[ii][jj]);

      for (i=0;i<intordsm[ii][jj];i++){
	xi=gp[i];
	for (j=0;j<intordsm[ii][jj];j++){
	  eta=gp[j];
	  
	  ipval[k]=approx (xi,eta,nodval);
	  k++;
	}
      }
      
      destrv (w);
      destrv (gp);
    }
  }
}





/**
   function defines meaning of DOFs at nodes
   
   @param eid - element id
   
   JK, 21.8.2005
*/
void planeelemlq::define_meaning (long eid)
{
  ivector cn(ndofe),nod(nne);
  
  Mt->give_elemnodes (eid,nod);
  Mt->give_code_numbers (eid,cn.a);

  //  displacement in x direction
  if (cn[0]>0)  Mt->nodes[nod[0]].meaning[0]=1;
  //  displacement in y direction
  if (cn[1]>0)  Mt->nodes[nod[0]].meaning[1]=2;
  //  displacement in x direction
  if (cn[2]>0)  Mt->nodes[nod[1]].meaning[0]=1;
  //  displacement in y direction
  if (cn[3]>0)  Mt->nodes[nod[1]].meaning[1]=2;
  //  displacement in x direction
  if (cn[4]>0)  Mt->nodes[nod[2]].meaning[0]=1;
  //  displacement in y direction
  if (cn[5]>0)  Mt->nodes[nod[2]].meaning[1]=2;
  //  displacement in x direction
  if (cn[6]>0)  Mt->nodes[nod[3]].meaning[0]=1;
  //  displacement in y direction
  if (cn[7]>0)  Mt->nodes[nod[3]].meaning[1]=2;

}








////////////////////       /* 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 planeelemlq::ntdbr_vector (long eid,vector &ntdbr)
{
  long intord = 2;
  long i,j,k,l,ipp,ippsh,ri,ci,lcid,body;
  double thick,xi,eta,jac,sh;
  ivector nodes(nne),cn;
  vector x(nne),y(nne),gp(intord),w(intord),t(nne),eps(tncomp),sig(tncomp),bf(nne),r;
  matrix d(tncomp,tncomp),gm;

  ri = ci = lcid = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];

  body = Ada->body[5];

  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord2d (x,y,eid);
  Mc->give_thickness (eid,nodes,t);

  Mm->matstiff (d,ipp);

  if (body){
    ippsh = ipp+4;
  }
  else{
    lcid = 0;
    allocv (ndofe,cn);
    allocv (ndofe,r);
    allocm (tncomp,ndofe,gm);
    
    Mt->give_code_numbers (eid,cn.a);
    eldispl (lcid,eid,r.a,cn.a,ndofe);
    
    geom_matrix (gm,x,y,0.0,0.0,jac);
    mxv (gm,r,eps);
    sh = eps[2];
  }
  
  gauss_points (gp.a,w.a,intord);
  
  fillv (0.0,ntdbr);
  
  for (i=0;i<intord;i++){
    xi=gp[i];
    for (j=0;j<intord;j++){
      eta=gp[j];
      
      if (body){
	eps[0] = Mm->ip[ipp  ].strain[0];
	eps[1] = Mm->ip[ipp  ].strain[1];
	eps[2] = Mm->ip[ippsh].strain[2];
	ipp++;
	jac_2d (jac,x,y,xi,eta);
      }
      else{
	geom_matrix (gm,x,y,xi,eta,jac);
	mxv (gm,r,eps);
	eps[2] = sh;
      }
      mxv (d,eps,sig);
      
      bf_lin_4_2d (bf.a,xi,eta);
      
      thick = approx (xi,eta,t);
      jac*=thick*w[i]*w[j];

      for (k=0;k<tncomp;k++)
        for (l=0;l<nne;l++)
          ntdbr[k*nne+l] += jac * bf[l] * sig[k];
      
    }
  }
}

/**
   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 planeelemlq::ntn_matrix (long eid,matrix &ntn)
{
  long intord = 2; //intordmm;
  long i,j,k,l;
  double thick,xi,eta,jac;
  ivector nodes(nne);
  vector x(nne),y(nne),t(nne),gp(intord),w(intord),bf(nne);
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord2d (x,y,eid);
  Mc->give_thickness (eid,nodes,t);
  
  gauss_points (gp.a,w.a,intord);
  
  fillm (0.0,ntn);
  
  for (i=0;i<intord;i++){
    xi=gp[i];
    for (j=0;j<intord;j++){
      eta=gp[j];
      
      thick = approx (xi,eta,t);
      jac_2d (jac,x,y,xi,eta);
      jac*=thick*w[i]*w[j];
      
      bf_lin_4_2d (bf.a,xi,eta);
      
      for (k=0;k<nne;k++)
        for (l=0;l<nne;l++)
          ntn[k][l] += jac * bf[k] * bf[l]; 
      
    }
  }
}

/**
   Function 1.
   1. computes "e2" - square of energy norm of error of solution on element
      e2 = integral{ (sig_star - sig_roof)T * D_inv * (sig_star - sig_roof) }
      sig_star = recovered stress, values in nodes are defined by "rsigfull" and over element are interpolated by base functions
      sig_roof = stress obtained by FEM
      D_inv = inverse stiffness matrix of material
   2. computes "u2" - square of energy norm of strain on element
      u2 = integral{ epsT * D * eps }
      eps = strain obtained by FEM
      D = stiffness matrix of material
   3. computes area of element and adds to "volume" (sum of areas of previous elements)
   4. computes "sizel" (characteristic size of element)
   
   @param eid - element id
   @param volume - sum of areas of previous elements
   @param e2 - empty(returned) value; 
   @param u2 - empty(returned) value; 
   @param sizel - empty(returned) value; 
   @param rsigfull - 1Darray of stress in nodes; dimmension is ncomp*gt->nn after this manner {(val[0]; ... ;val[gt->nn])[0] ; ... ; (val[0]; ... ;val[gt->nn])[ncomp]}
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
 */
void planeelemlq::compute_error (long eid,double &volume,double &e2,double &u2,double &sizel,double *rsigfull)
{
  long intord = 2;
  long i,j,ipp,ippsh,ri,ci,body,lcid;
  double thick,area,xi,eta,jac,contr;
  double zero=1.0e-20;
  ivector nodes(nne),cn;
  vector x(nne),y(nne),t(nne),gp(intord),w(intord),bf(nne),r;
  vector sig_star(tncomp),sig_roof(tncomp),sig_err(tncomp),eps(tncomp);
  matrix d(tncomp,tncomp),dinv(tncomp,tncomp),gm;
  
  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  
  body = Ada->body[5];

  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord2d (x,y,eid);
  Mc->give_thickness (eid,nodes,t);

  Mm->matstiff (d,ipp);
  invm (d,dinv,zero);

  if (body){
    ippsh = ipp+4;
  }
  else{
    lcid = 0;
    allocv (ndofe,cn);
    allocv (ndofe,r);
    allocm (tncomp,ndofe,gm);
    
    Mt->give_code_numbers (eid,cn.a);
    eldispl (lcid,eid,r.a,cn.a,ndofe);
    
    geom_matrix (gm,x,y,0.0,0.0,jac);
    mxv (gm,r,eps);
    area = eps[2];
  }
  
  gauss_points (gp.a,w.a,intord);
  
  e2 = u2 = 0;
  
  for (i=0;i<intord;i++){
    xi=gp[i];
    for (j=0;j<intord;j++){
      eta=gp[j];
      
      if (body){
	eps[0] = Mm->ip[ipp  ].strain[0];
	eps[1] = Mm->ip[ipp  ].strain[1];
	eps[2] = Mm->ip[ippsh].strain[2];
	ipp++;
	jac_2d (jac,x,y,xi,eta);
      }
      else{
	geom_matrix (gm,x,y,xi,eta,jac);
	mxv (gm,r,eps);
	eps[2] = area;
      }
      mxv (d,eps,sig_roof);
      
      bf_lin_4_2d (bf.a,xi,eta);
      give_der_star (bf,rsigfull,nodes,sig_star,Mt->nn);
      
      subv (sig_star,sig_roof,sig_err);
      
      thick = approx (xi,eta,t);
      jac*=thick*w[i]*w[j];
      
      vxmxv (sig_err,dinv,contr);
      e2 += jac*contr;
      
      vxmxv (eps,d,contr);
      u2 += jac*contr;
    }
  }
  
  area = ( (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]) + (x[2]-x[0])*(y[3]-y[0])-(x[3]-x[0])*(y[2]-y[0]) ) / 2.0;
  volume += area;
  sizel = sqrt(area);
}

/**
   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 planeelemlq::elchar (long eid,double *&spsig)
{
  long ipp,ri,ci,lcid;
  double jac;
  ivector cn(ndofe);
  vector x(nne),y(nne),r(ndofe),eps(tncomp);
  matrix gm(tncomp,ndofe),d(tncomp,tncomp);
  
  spsig = new double[tncomp];
  
  lcid = ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  
  Mm->matstiff (d,ipp);
  
  Mt->give_node_coord2d (x,y,eid);
  Mt->give_code_numbers (eid,cn.a);
  eldispl (lcid,eid,r.a,cn.a,ndofe);
  
  geom_matrix (gm,x,y,0.0,0.0,jac);
  mxv (gm,r,eps);
  mxv (d.a,eps.a,spsig,d.m,d.n);
}
////////////////////       /* termitovo */       ////////////////////////////////////


/**
   function extracts variables on element
   it is used in stochastic or fuzzy computations
   
   @param at - array of attributes (it describes which variables will be extracted)
   @param val - %vector of extracted values
   
   JK, 17.4.2007
*/
void planeelemlq::extract (atsel &at,vector &val)
{
  long i;
  
  for (i=0;i<at.num;i++){
    switch (at.atrib[i]){
    case 0:{
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function extract (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  
}