plelemqq.cpp

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

      coord[k][1]=approx (xi,eta,y);
      coord[k++][2]=0.0;
    }
  }
}

/**
   function computes nodal forces caused by edge load
   
   @param eid - element id
   @param le - list of loaded edges
   @param nv - nodal values of edge loads, it contains 24 components
   @param nf - %vector of nodal forces

   nv contains 2 values at each node on each edge, each edge contains 3 nodes,
   element contains 4 edges, therefore, nv contains 4*3*2 = 24 components

   JK, modification 23.11.2006
*/
void planeelemqq::res_nodeforces (long eid,long *le,double *nv,vector &nf)
{
  vector x(nne),y(nne);
  Mt->give_node_coord2d (x,y,eid);
  nodeforces (eid,le,nv,nf,x,y);
}

/**
   function computes nodal forces caused by edge load
   
   this function is used in plane stress/strain elements (function is called
   by function res_nodeforces) and shell elements

   @param eid - element id
   @param le - list of loaded edges
   @param nv - nodal values of edge loads, it contains 24 components
   @param nf - %vector of nodal forces
   @param x,y - node coordinates

   nv contains 2 values at each node on each edge, each edge contains 3 nodes,
   element contains 4 edges, therefore, nv contains 4*3*2 = 24 components

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

  gauss_points (gp.a,w.a,intordb);

  if (le[0]==1){
    //  first edge is loaded
    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[8]=nv[2];  av[9]=nv[3];  av[2]=nv[4];  av[3]=nv[5];
    mxv (am,av,v);  addv (nf,v,nf);
  }
  if (le[1]==1){
    //  second edge is loaded
    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[6];  av[3]=nv[7];  av[10]=nv[8];  av[11]=nv[9];  av[4]=nv[10];  av[5]=nv[11];
    mxv (am,av,v);  addv (nf,v,nf);
  }
  if (le[2]==1){
    //  third edge is loaded
    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[12];  av[5]=nv[13];  av[12]=nv[14];  av[13]=nv[15];  av[6]=nv[16];  av[7]=nv[17];
    mxv (am,av,v);  addv (nf,v,nf);
  }
  if (le[3]==1){
    //  fourth edge is loaded
    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[18];  av[7]=nv[19];  av[14]=nv[20];  av[15]=nv[21];  av[0]=nv[22];  av[1]=nv[23];
    mxv (am,av,v);  addv (nf,v,nf);
  }
}


/**
   function computes displacement, strain and stress in the centre=middle of element
   
   @param eid - number of element
   @paran dd  - decomposed deformation
   @param valel - array of displacements, strains and stresses in the middle of all elements
   
   created  20.11.2002, Ladislav Svoboda, [e-mail]
 */
void planeelemqq::midpoints (long eid,long dd,double *valel)
{
  long i,j,lcid;
  double r[2];
  vector nodval(nne),der(tncomp);
  ivector nodes(nne);
  
  lcid = 0;
  Mt->give_elemnodes (eid,nodes);  
  
  for (i=0;i<2;i++){
    for (j=0;j<nne;j++){
      noddispl (lcid,r,nodes[j]);
      nodval[j] = r[i];
    }
    valel[i] = approx (0.0,0.0,nodval);
    if (dd) valel[i+2] = valel[i];
  }

  //appstrain (lcid,eid,0.0,0.0,0,tncomp,der);
  for (i=0;i<tncomp;i++)
    valel[i+2*(1+dd)] = der[i];
  
  //appstress (lcid,eid,0.0,0.0,0,tncomp,der);
  for (i=0;i<tncomp;i++)
    valel[i+2*(1+dd)+tncomp] = der[i];
}


void planeelemqq::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;
  double xi, eta, ipval;
  vector w, gp, 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],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);
            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 < nv))
            {
              Mm->ip[ipp].other[j] += ipval;
              ipp++;
              continue;
            }
            ipp++;

          }
        }
        destrv (gp);  destrv (w);
      }
    }
    if (ictn[i] & inistrain) nstra++;
  }
}

/**
   function computes volume appropriate to integration point
   
   2.3.2004, JK
*/
void planeelemqq::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
   quadratic approximation functions are used

   @param eid - element id
   @param nodval - nodal values
   @param ipval - value at integration points
   
   21.6.2004, JK
*/
void planeelemqq::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 interpolates the nodal values to the integration points on the element
   linear approximation functions are used
   
   @param eid - element id
   @param nodval - nodal values
   @param ipval - value at integration points
   
   21.6.2004, JK
*/
void planeelemqq::intpointval2 (long eid,vector &nodval,vector &ipval)
{
  long i,j,ii,jj,k;
  double xi,eta;
  vector w,gp;
  vector modnodval(Pelq->nne);
  
  for (i=0;i<Pelq->nne;i++){
    modnodval[i]=nodval[i];
  }
  
  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]=Pelq->approx (xi,eta,modnodval);
	  k++;
	}
      }
      destrv (w);
      destrv (gp);
    }
  }
}



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

  ri = ci = lcid = 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);
  Mt->give_code_numbers (eid,cn.a);
  eldispl (lcid,eid,r.a,cn.a,ndofe);

  Mm->matstiff (d,ipp);

  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];

      geom_matrix (gm,x,y,xi,eta,jac);
      mxv (gm,r,eps);
      mxv (d,eps,sig);

      bf_quad_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 planeelemqq::ntn_matrix (long eid,matrix &ntn)
{
  long i,j,k,l,intord = 3;
  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_quad_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 planeelemqq::compute_error (long eid,double &volume,double &e2,double &u2,double &sizel,double *rsigfull)
{
  long intord = 3;
  long i,j,ipp,ri,ci,lcid;
  double area,thick,xi,eta,jac,contr;
  double zero=1.0e-20;
  ivector cn(ndofe),nodes(nne);
  vector x(nne),y(nne),t(nne),gp(intord),w(intord),bf(nne),r(ndofe);
  vector sig_star(tncomp),sig_roof(tncomp),sig_err(tncomp),eps(tncomp);
  matrix gm(tncomp,ndofe),d(tncomp,tncomp),dinv(tncomp,tncomp);

  ri = ci = lcid = 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);
  Mt->give_code_numbers (eid,cn.a);
  eldispl (lcid,eid,r.a,cn.a,ndofe);

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

  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];

      bf_quad_4_2d (bf.a,xi,eta);
      give_der_star (bf,rsigfull,nodes,sig_star,Mt->nn);
      
      geom_matrix (gm,x,y,xi,eta,jac);
      mxv (gm,r,eps);
      mxv (d,eps,sig_roof);

      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 planeelemqq::elchar (long eid,double *&spsig)
{
  long intord = 2;
  long i,j,ipp,ri,ci,body;
  double *sig;
  vector eps(tncomp);
  matrix d(tncomp,tncomp);
  
  spsig = new double[tncomp*4];
  
  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  Mm->matstiff (d,ipp);

  body = Ada->body[7];

  if (body){
    for (i=0;i<intord;i++)
      for (j=0;j<intord;j++){
	eps[0] = Mm->ip[ipp  ].strain[0];
	eps[1] = Mm->ip[ipp  ].strain[1];
	eps[2] = Mm->ip[ipp+4].strain[2];
	ipp++;
	sig = spsig + (i*intord+j)*tncomp;
	mxv (d.a,eps.a,sig,d.m,d.n);
      }
  }
  else{
    long lcid = 0;
    double jac;
    ivector cn(ndofe);
    vector gp(intord),w(intord),x(nne),y(nne),r(ndofe);
    matrix gm(tncomp,ndofe);
    
    Mt->give_node_coord2d (x,y,eid);
    Mt->give_code_numbers (eid,cn.a);
    eldispl (lcid,eid,r.a,cn.a,ndofe);
    gauss_points (gp.a,w.a,intord);
    
    for (i=0;i<intord;i++)
      for (j=0;j<intord;j++){
	geom_matrix (gm,x,y,gp[i],gp[j],jac);
	mxv (gm,r,eps);
	sig = spsig + (i*intord+j)*tncomp;
	mxv (d.a,eps.a,sig,d.m,d.n);
      }
  }
}
////////////////////       /* termitovo */       ////////////////////////////////////