plelemqt.cpp

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

      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 (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(gp1);  destrv (gp2);  destrv (w);
      }
    }
    if (ictn[i] & inistrain) nstra++;
  }
}


/**
   function computes volume appropriate to integration point
   
   2.3.2004, JK
*/
void planeelemqt::ipvolume (long eid,long ri,long ci)
{
  long i,ii,jj,ipp;
  double jac;
  vector x(nne),y(nne),gp1,gp2,w;

  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],gp1);
      allocv (intordsm[ii][jj],gp2);
      allocv (intordsm[ii][jj],w);

      gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
      ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
      
      for (i=0;i<intordsm[ii][jj];i++){
	jac_2d (jac,x,y,gp1[i],gp2[i]);
	jac*=w[i];
	
	Mm->storeipvol (ipp,jac);
	
	ipp++;
      }
      destrv (gp1);  destrv (gp2);  destrv (w);
    }
  }
}


////////////////////       /* 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 planeelemqt::ntdbr_vector (long eid,vector &ntdbr)
{
  long intord = 4;
  long i,k,l,ipp,ri,ci,lcid;
  double thick,jac,**stra;
  ivector nodes(nne);
  vector x(nne),y(nne),gp1(intord),gp2(intord),w(intord);
  vector t(nne),sig(tncomp),eps(tncomp),bf(nne);
  matrix d(tncomp,tncomp);

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

  stra = new double* [nne];
  for (i=0;i<nne;i++){
    stra[i] = new double [tncomp];
  }
  elem_strains (stra,lcid,eid,ri,ci);

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

  gauss_points_tr (gp1.a,gp2.a,w.a,intord);

  fillv(0.0,ntdbr);

  for (i=0;i<intord;i++){
    appval (gp1[i],gp2[i],0,tncomp,eps,stra);
    mxv (d,eps,sig);

    bf_quad_3_2d (bf.a,gp1[i],gp2[i]);

    thick = approx (gp1[i],gp2[i],t);
    jac_2d (jac,x,y,gp1[i],gp2[i]);
    jac *= w[i]*thick;

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

  }

  for (i=0;i<nne;i++)  delete [] stra[i];
                       delete [] stra;
}

/**
   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 planeelemqt::ntn_matrix (long eid,matrix &ntn)
{
  long intord = 6;
  long i,k,l;
  double thick,jac;
  ivector nodes(nne);
  vector x(nne),y(nne),t(nne),gp1(intord),gp2(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_tr (gp1.a,gp2.a,w.a,intord);
  
  fillm(0.0,ntn);
  
  for (i=0;i<intord;i++){
    bf_quad_3_2d (bf.a,gp1[i],gp2[i]);
    
    thick = approx (gp1[i],gp2[i],t);
    jac_2d (jac,x,y,gp1[i],gp2[i]);
    jac *= w[i]*thick;
    
    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 planeelemqt::compute_error (long eid,double &volume,double &e2,double &u2,double &sizel,double *rsigfull)
{
  long intorde = 7;
  long i,ipp,ri,ci,lcid;
  double area,thick,jac,contr,**stra;
  double zero=1.0e-20;
  ivector nodes(nne);
  vector x(nne),y(nne),t(nne),gp1(intorde),gp2(intorde),w(intorde),bf(nne);
  vector sig_star(tncomp),sig_roof(tncomp),sig_err(tncomp),eps(tncomp);
  matrix d(tncomp,tncomp),dinv(tncomp,tncomp);
  
  ri = ci = lcid = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  
  stra = new double* [nne];
  for (i=0;i<nne;i++){
    stra[i] = new double [tncomp];
  }
  elem_strains (stra,lcid,eid,ri,ci);
  
  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);
  
  gauss_points_tr (gp1.a,gp2.a,w.a,intorde);
  
  e2 = u2 = 0;
  for (i=0;i<intorde;i++){
    thick = approx (gp1[i],gp2[i],t);
    jac_2d (jac,x,y,gp1[i],gp2[i]);
    jac *= w[i]*thick;
    
    bf_quad_3_2d (bf.a,gp1[i],gp2[i]);
    give_der_star (bf,rsigfull,nodes,sig_star,Mt->nn);
    
    appval (gp1[i],gp2[i],0,tncomp,eps,stra);
    mxv (d,eps,sig_roof);
    
    subv (sig_star,sig_roof,sig_err);
    
    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]);
  volume += area/2.0;
  sizel = sqrt(1.1547005384*area);
  
  for (i=0;i<nne;i++)  delete [] stra[i];
                       delete [] stra;
}

/**
   Only for Termit
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
 */
void planeelemqt::give_sig_roof (matrix &d,double *areacoord,double *x,double *y,long *etnodes,double **xyrr,vector &sig_roof,vector &sig_star)
{                                   //pozor na zmenu poradi ip
  long i,nnecm;
  double xx,yy,detcm,shear,xicm,etacm,auxd;
  vector areacoordcm(3),eps(3),xcm,ycm,rcm,nodsig_x,nodsig_y,nodsig_xy;
  matrix gm;

  nnecm = etnodes[1];

  allocv (nnecm,xcm);
  allocv (nnecm,ycm);
  allocv (nnecm,nodsig_x);
  allocv (nnecm,nodsig_y);
  allocv (nnecm,nodsig_xy);
  allocv (2*nnecm,rcm);
  allocm (3,2*nnecm,gm);

  for (i=0;i<nnecm;i++){
    xcm[i]       = xyrr[etnodes[i+2]][0];
    ycm[i]       = xyrr[etnodes[i+2]][1];
    rcm[2*i]     = xyrr[etnodes[i+2]][2];
    rcm[2*i+1]   = xyrr[etnodes[i+2]][3];
    nodsig_x[i]  = xyrr[etnodes[i+2]][4];
    nodsig_y[i]  = xyrr[etnodes[i+2]][5];
    nodsig_xy[i] = xyrr[etnodes[i+2]][6];
  }

  xx = yy = 0.0;
  for (i=0;i<3;i++){
    xx += areacoord[i] * x[i];
    yy += areacoord[i] * y[i];
  }

  switch (etnodes[0]){
  case planeelementlt:{
    detcm = (xcm[1]-xcm[0])*(ycm[2]-ycm[0])-(xcm[2]-xcm[0])*(ycm[1]-ycm[0]);
    areacoordcm[0] = ( (xcm[1]*ycm[2] - xcm[2]*ycm[1]) + (ycm[1] - ycm[2])*xx + (xcm[2] - xcm[1])*yy)/detcm;
    areacoordcm[1] = ( (xcm[2]*ycm[0] - xcm[0]*ycm[2]) + (ycm[2] - ycm[0])*xx + (xcm[0] - xcm[2])*yy)/detcm;
    areacoordcm[2] = ( (xcm[0]*ycm[1] - xcm[1]*ycm[0]) + (ycm[0] - ycm[1])*xx + (xcm[1] - xcm[0])*yy)/detcm;  //= 1-ar[0]-ar[1]

    Pelt->geom_matrix (gm,xcm,ycm);
    mxv (gm,rcm,eps);
    mxv (d,eps,sig_roof);

    sig_star[0] = Pelt->approx (areacoordcm,nodsig_x);
    sig_star[1] = Pelt->approx (areacoordcm,nodsig_y);
    sig_star[2] = Pelt->approx (areacoordcm,nodsig_xy);
    break;
  }
  case planeelementlq:{
    nc_lin_4_2d (1e-6,xx,yy,xcm.a,ycm.a,xicm,etacm);
    
    Pelq->geom_matrix (gm,xcm,ycm,0.0,0.0,auxd);
    mxv (gm,rcm,eps);
    mxv (d,eps,sig_roof);
    shear=sig_roof[2];
  
    Pelq->geom_matrix (gm,xcm,ycm,xicm,etacm,auxd);
    mxv (gm,rcm,eps);
    mxv (d,eps,sig_roof);
    sig_roof[2]=shear;

    sig_star[0] = Pelq->approx (xicm,etacm,nodsig_x);
    sig_star[1] = Pelq->approx (xicm,etacm,nodsig_y);
    sig_star[2] = Pelq->approx (xicm,etacm,nodsig_xy);
    break;
  }
  case planeelementqq:{
    nc_lin_4_2d (1e-6,xx,yy,xcm.a,ycm.a,xicm,etacm);   // funguje jen pro subparametriclou sit = rovne strany ctyruhelnika

    Peqq->geom_matrix (gm,xcm,ycm,xicm,etacm,auxd);
    mxv (gm,rcm,eps);
    mxv (d,eps,sig_roof);

    sig_star[0] = Peqq->approx (xicm,etacm,nodsig_x);
    sig_star[1] = Peqq->approx (xicm,etacm,nodsig_y);
    sig_star[2] = Peqq->approx (xicm,etacm,nodsig_xy);
    break;
  }
  case planeelementqt:{
    detcm = (xcm[1]-xcm[0])*(ycm[2]-ycm[0])-(xcm[2]-xcm[0])*(ycm[1]-ycm[0]);
    areacoordcm[0] = ( (xcm[1]*ycm[2] - xcm[2]*ycm[1]) + (ycm[1] - ycm[2])*xx + (xcm[2] - xcm[1])*yy)/detcm;
    areacoordcm[1] = ( (xcm[2]*ycm[0] - xcm[0]*ycm[2]) + (ycm[2] - ycm[0])*xx + (xcm[0] - xcm[2])*yy)/detcm;
    areacoordcm[2] = ( (xcm[0]*ycm[1] - xcm[1]*ycm[0]) + (ycm[0] - ycm[1])*xx + (xcm[1] - xcm[0])*yy)/detcm;

    //Pesqt->geom_matrix (gm,xcm,ycm,areacoordcm);
    mxv (gm,rcm,eps);
    mxv (d,eps,sig_roof);

    sig_star[0] = approx (xicm,etacm,nodsig_x);
    sig_star[1] = approx (xicm,etacm,nodsig_y);
    sig_star[2] = approx (xicm,etacm,nodsig_xy);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function");
    fprintf (stderr," planeelemqt::give_sig_roof (%s, line %d)",__FILE__,__LINE__);
  }
  }
}

/**
   Only for Termit
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
 */
void planeelemqt::compute_error_fine (long eid,long *etnodes,double **xyrr,double *rsigfull,double &e2)
{
  long intord = 7;
  long i,ipp,ri,ci;
  double thick,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 sig_fine(tncomp),sig_star(tncomp),sig_roof(tncomp),sig_err(tncomp),bf(nne);
  matrix dinv(tncomp,tncomp),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);
  
  gauss_points_tr (gp1.a,gp2.a,w.a,intord);

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

  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 (gp1[i],gp2[i],t);
    jac_2d (jac,x,y,gp1[i],gp2[i]);
    jac *= w[i]*thick;

    //Mm->givestress (ipp,sig_fine);
    //ipp++;


    //bf_quad_3_2d (bf.a,gp1[i],gp2[i]);
    ac_bf_quad_3_2d (bf.a,areacoord.a);
    give_der_star (bf,rsigfull,nodes,sig_fine,Mt->nn);
    give_sig_roof (d,areacoord.a,x.a,y.a,etnodes,xyrr,sig_roof,sig_star);

    subv (sig_fine,sig_roof,sig_err);

    vxmxv (sig_err,dinv,contr);
    e2 += jac*contr;
  }
}

/**
   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 planeelemqt::elchar (long eid,double *&spsig)
{
  long intord = 3;
  long i,ipp,ri,ci,body;
  double *sig;
  vector eps(tncomp);
  matrix d(tncomp,tncomp);

  spsig = new double[tncomp*3];

  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  Mm->matstiff (d,ipp);

  body = Ada->body[2];
  
  if (body){
    for (i=0;i<intord;i++){
      eps[0] = Mm->ip[ipp  ].strain[0];
      eps[1] = Mm->ip[ipp  ].strain[1];
      eps[2] = Mm->ip[ipp+3].strain[2];
      ipp++;
      sig = spsig + i*tncomp;
      mxv (d.a,eps.a,sig,d.m,d.n);
    }
  }
  else{
    long lcid = 0;
    double **stra;
    vector gp1(intord),gp2(intord),w(intord);
    
    stra = new double* [nne];
    for (i=0;i<nne;i++){
      stra[i] = new double [tncomp];
    }
    elem_strains (stra,lcid,eid,ri,ci);
    
    gauss_points_tr (gp1.a,gp2.a,w.a,intord);
    
    for (i=0;i<intord;i++){
      appval (gp1[i],gp2[i],0,tncomp,eps,stra);
      sig = spsig + i*tncomp;
      mxv (d.a,eps.a,sig,d.m,d.n);
    }

    for (i=0;i<nne;i++)  delete [] stra[i];
                         delete [] stra;
  }
}
////////////////////       /* termitovo */       ////////////////////////////////////