lintet.cpp

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

    
    av[6] = nv[9+1*3+0];
    av[7] = nv[9+1*3+1];
    av[8] = nv[9+1*3+2];
    
    av[9]  = nv[9+2*3+0];
    av[10] = nv[9+2*3+1];
    av[11] = nv[9+2*3+2];
    
    mxv (am,av,v);
    addv (v,nf,nf);
  }
  if (is[2]==1){
    for (i=0;i<intordb;i++){
      xi=gp1[i]; eta=gp2[i];  zeta=0.0;
      
      jac2d_3d (jac,x,y,z,xi,eta,2);
      bf_matrix (n,xi,eta,zeta);
      jac = jac*w[i];
      nnj (am.a,n.a,jac,n.m,n.n);
    }
    
    av[3] = nv[18+0*3+0];
    av[4] = nv[18+0*3+1];
    av[5] = nv[18+0*3+2];
    
    av[0] = nv[18+1*3+0];
    av[1] = nv[18+1*3+1];
    av[2] = nv[18+1*3+2];
    
    av[9]  = nv[18+2*3+0];
    av[10] = nv[18+2*3+1];
    av[11] = nv[18+2*3+2];
    
    mxv (am,av,v);
    addv (v,nf,nf);
  }
  if (is[3]==1){
    for (i=0;i<intordb;i++){
      xi=gp1[i]; eta=gp2[i];  zeta=1.0-xi-eta;
      
      jac2d_3d (jac,x,y,z,xi,eta,3);
      bf_matrix (n,xi,eta,zeta);
      jac = jac*w[i];
      nnj (am.a,n.a,jac,n.m,n.n);
    }
    
    av[3] = nv[27+0*3+0];
    av[4] = nv[27+0*3+1];
    av[5] = nv[27+0*3+2];
    
    av[6] = nv[27+1*3+0];
    av[7] = nv[27+1*3+1];
    av[8] = nv[27+1*3+2];
    
    av[0] = nv[27+2*3+0];
    av[1] = nv[27+2*3+1];
    av[2] = nv[27+2*3+2];
    
    mxv (am,av,v);
    addv (v,nf,nf);
  }
}


/**
   function computes nodal forces caused by presure on surface
   
   @param lcid - number of load case
   @param eid - element id
   @param is - identification of surfaces
   @param nv - nodal values
   @param nf - nodal forces
   
   12.12.2005, JK
*/
void lintet::node_forces_surf (long lcid,long eid,long *is,double *nv,vector &nf)
{
  long i;
  double xi,eta,zeta,jac;
  double *tnv;
  vector x(nne),y(nne),z(nne),v(ndofe),av(ndofe),gp1(intordb),gp2(intordb),w(intordb);
  matrix am(ndofe,ndofe),n(napfun,ndofe);
  
  tnv = new double [9];
  
  //  coordinates of element nodes
  Mt->give_node_coord3d (x,y,z,eid);
  
  gauss_points_tr (gp1.a,gp2.a,w.a,intordb);
  
  if (is[0]>0 ){
    for (i=0;i<intordb;i++){
      xi=0.0; eta=gp1[i];  zeta=gp2[i];
      
      jac2d_3d (jac,x,y,z,eta,zeta,0);
      bf_matrix (n,xi,eta,zeta);
      jac = jac*w[i];
      nnj (am.a,n.a,jac,n.m,n.n);
    }
    
    if (is[0]==1){
      av[6] = nv[0*3+0];
      av[7] = nv[0*3+1];
      av[8] = nv[0*3+2];
      
      av[3] = nv[1*3+0];
      av[4] = nv[1*3+1];
      av[5] = nv[1*3+2];
      
      av[9] = nv[2*3+0];
      av[10] = nv[2*3+1];
      av[11] = nv[2*3+2];
    }


    if (is[0]==2){
      
      av[0] = nv[0*3+0];
      av[1] = nv[0*3+1];
      av[2] = nv[0*3+2];
      
      av[3] = nv[1*3+0];
      av[4] = nv[1*3+1];
      av[5] = nv[1*3+2];
      
      av[6] = nv[2*3+0];
      av[7] = nv[2*3+1];
      av[8] = nv[2*3+2];
      
      locglob_nodeval (0,av,tnv,x,y,z);
      
      av[6] = tnv[0];
      av[7] = tnv[1];
      av[8] = tnv[2];
      
      av[3] = tnv[3];
      av[4] = tnv[4];
      av[5] = tnv[5];
      
      av[9] = tnv[6];
      av[10] = tnv[7];
      av[11] = tnv[8];
    }


    mxv (am,av,v);
    addv (v,nf,nf);
  }
  if (is[1]>0){
    for (i=0;i<intordb;i++){
      xi=gp1[i]; eta=0.0;  zeta=gp2[i];
      
      jac2d_3d (jac,x,y,z,xi,zeta,1);
      bf_matrix (n,xi,eta,zeta);
      jac = jac*w[i];
      nnj (am.a,n.a,jac,n.m,n.n);
    }
    
    if (is[1]==1){
      av[0] = nv[9+0*3+0];
      av[1] = nv[9+0*3+1];
      av[2] = nv[9+0*3+2];
      
      av[6] = nv[9+1*3+0];
      av[7] = nv[9+1*3+1];
      av[8] = nv[9+1*3+2];
      
      av[9] = nv[9+2*3+0];
      av[10] = nv[9+2*3+1];
      av[11] = nv[9+2*3+2];
    }
    
    if (is[1]==2){
      
      av[0] = nv[9+0*3+0];
      av[1] = nv[9+0*3+1];
      av[2] = nv[9+0*3+2];
      
      av[3] = nv[9+1*3+0];
      av[4] = nv[9+1*3+1];
      av[5] = nv[9+1*3+2];
      
      av[6] = nv[9+2*3+0];
      av[7] = nv[9+2*3+1];
      av[8] = nv[9+2*3+2];
      
      locglob_nodeval (1,av,tnv,x,y,z);
      
      av[0] = tnv[0];
      av[1] = tnv[1];
      av[2] = tnv[2];
      
      av[6] = tnv[3];
      av[7] = tnv[4];
      av[8] = tnv[5];
      
      av[9] = tnv[6];
      av[10] = tnv[7];
      av[11] = tnv[8];
    }
    
    mxv (am,av,v);
    addv (v,nf,nf);
  }
  if (is[2]>0){
    for (i=0;i<intordb;i++){
      xi=gp1[i]; eta=gp2[i];  zeta=0.0;
      
      jac2d_3d (jac,x,y,z,xi,eta,2);
      bf_matrix (n,xi,eta,zeta);
      jac = jac*w[i];
      nnj (am.a,n.a,jac,n.m,n.n);
    }
    
    if (is[2]==1){
      av[3] = nv[18+0*3+0];
      av[4] = nv[18+0*3+1];
      av[5] = nv[18+0*3+2];
      
      av[0] = nv[18+1*3+0];
      av[1] = nv[18+1*3+1];
      av[2] = nv[18+1*3+2];
      
      av[9] = nv[18+2*3+0];
      av[10] = nv[18+2*3+1];
      av[11] = nv[18+2*3+2];

    }

    if (is[2]==2){
      
      av[0] = nv[18+0*3+0];
      av[1] = nv[18+0*3+1];
      av[2] = nv[18+0*3+2];
      
      av[3] = nv[18+1*3+0];
      av[4] = nv[18+1*3+1];
      av[5] = nv[18+1*3+2];
      
      av[6] = nv[18+2*3+0];
      av[7] = nv[18+2*3+1];
      av[8] = nv[18+2*3+2];
      
      locglob_nodeval (2,av,tnv,x,y,z);
      
      av[3] = tnv[0];
      av[4] = tnv[1];
      av[5] = tnv[2];
      
      av[0] = tnv[3];
      av[1] = tnv[4];
      av[2] = tnv[5];
      
      av[9] = tnv[6];
      av[10] = tnv[7];
      av[11] = tnv[8];
    }

    mxv (am,av,v);
    addv (v,nf,nf);
  }
  if (is[3]>0){
    for (i=0;i<intordb;i++){
      xi=gp1[i]; eta=gp2[i];  zeta=1.0-xi-eta;
      
      jac2d_3d (jac,x,y,z,xi,eta,3);
      bf_matrix (n,xi,eta,zeta);
      jac = jac*w[i];
      nnj (am.a,n.a,jac,n.m,n.n);
    }
    
    if (is[3]==1){
      av[3] = nv[27+0*3+0];
      av[4] = nv[27+0*3+1];
      av[5] = nv[27+0*3+2];
      
      av[6] = nv[27+1*3+0];
      av[7] = nv[27+1*3+1];
      av[8] = nv[27+1*3+2];
      
      av[0] = nv[27+2*3+0];
      av[1] = nv[27+2*3+1];
      av[2] = nv[27+2*3+2];

    }

    if (is[3]==2){
      
      av[0] = nv[27+0*3+0];
      av[1] = nv[27+0*3+1];
      av[2] = nv[27+0*3+2];
      
      av[3] = nv[27+1*3+0];
      av[4] = nv[27+1*3+1];
      av[5] = nv[27+1*3+2];
      
      av[6] = nv[27+2*3+0];
      av[7] = nv[27+2*3+1];
      av[8] = nv[27+2*3+2];
      
      locglob_nodeval (3,av,tnv,x,y,z);
      
      av[3] = tnv[0];
      av[4] = tnv[1];
      av[5] = tnv[2];
      
      av[6] = tnv[3];
      av[7] = tnv[4];
      av[8] = tnv[5];
      
      av[0] = tnv[6];
      av[1] = tnv[7];
      av[2] = tnv[8];
    }
    
    mxv (am,av,v);
    addv (v,nf,nf);
  }
  
  delete [] tnv;
}


void lintet::locglob_nodeval (long is,vector &nv,double *tnv,vector &x,vector &y,vector &z)
{
  double norm;
  vector g1(3),g2(3),g3(3),lv(3),gv(3);
  matrix t(3,3);
  
  if (is==0){
    g1[0]=x[2]-x[3];
    g1[1]=y[2]-y[3];
    g1[2]=z[2]-z[3];
    
    g2[0]=x[1]-x[3];
    g2[1]=y[1]-y[3];
    g2[2]=z[1]-z[3];
  }
  if (is==1){
    g1[0]=x[0]-x[3];
    g1[1]=y[0]-y[3];
    g1[2]=z[0]-z[3];
    
    g2[0]=x[2]-x[3];
    g2[1]=y[2]-y[3];
    g2[2]=z[2]-z[3];
  }
  if (is==2){
    g1[0]=x[1]-x[3];
    g1[1]=y[1]-y[3];
    g1[2]=z[1]-z[3];
    
    g2[0]=x[0]-x[3];
    g2[1]=y[0]-y[3];
    g2[2]=z[0]-z[3];
  }
  if (is==3){
    g1[0]=x[1]-x[0];
    g1[1]=y[1]-y[0];
    g1[2]=z[1]-z[0];
    
    g2[0]=x[2]-x[0];
    g2[1]=y[2]-y[0];
    g2[2]=z[2]-z[0];
  }
  
  
  scprd (g1,g1,norm);
  norm=1.0/sqrt(norm);
  cmulv (norm,g1,g1);
  
  scprd (g1,g2,norm);
  g2[0]=g2[0]-norm*g1[0];
  g2[1]=g2[1]-norm*g1[1];
  g2[2]=g2[2]-norm*g1[2];
  
  scprd (g2,g2,norm);
  norm=1.0/sqrt(norm);
  cmulv (norm,g2,g2);
  
  g3[0]=g1[1]*g2[2]-g1[2]*g2[1];
  g3[1]=g1[2]*g2[0]-g1[0]*g2[2];
  g3[2]=g1[0]*g2[1]-g1[1]*g2[0];
  
  t[0][0]=g1[0];
  t[1][0]=g1[1];
  t[2][0]=g1[2];

  t[0][1]=g2[0];
  t[1][1]=g2[1];
  t[2][1]=g2[2];

  t[0][2]=g3[0];
  t[1][2]=g3[1];
  t[2][2]=g3[2];
  
  mxv (t.a,nv.a,tnv,3,3);
  mxv (t.a,nv.a+3,tnv+3,3,3);
  mxv (t.a,nv.a+6,tnv+6,3,3);
  
}

/**
   function interpolates the nodal values to the integration points on the element
   
   @param eid - element id
   
   JK, 22.4.2005
*/
void lintet::intpointval (long eid,vector &nodval,vector &ipval)
{
  vector volcoord(4);
  
  volcoord[0]=0.25;
  volcoord[1]=0.25;
  volcoord[2]=0.25;
  volcoord[3]=0.25;
  
  ipval[0]=approx (volcoord,nodval);
}


/**
   function computes contributions from eigenstrains to the right hand side
   
   @param eid - element id
   @param ifor - vector of internal forces

   JK, 22.4.2005
*/
void lintet::res_eigstrain_forces (long eid,vector &nfor)
{
  eigstrain_forces (eid,nfor);
}

/**
   function computes contributions from eigenstrains to the right hand side
   
   @param eid - element id
   @param ifor - vector of internal forces

   JK, 22.4.2005
*/
void lintet::eigstrain_forces (long eid,vector &nfor)
{
  long i,ipp;
  double jac,det;
  vector x(nne),y(nne),z(nne),eigstr(tncomp),sig(tncomp),contr(ndofe);
  matrix d(tncomp,tncomp),gm(tncomp,ndofe);
  
  Mt->give_node_coord3d (x,y,z,eid);
  det = det3d (x.a,y.a,z.a);
  
  fillv (0.0,nfor);
  
  ipp=Mt->elements[eid].ipp[0][0];
  
  Mm->giveeigstrain (ipp,cncomp[0],ncomp[0],eigstr);
  
  //  matrix of stiffness of the material
  Mm->matstiff (d,ipp);
  
  mxv (d,eigstr,sig);
  
  geom_matrix (gm,x,y,z);
  
  mtxv (gm,sig,contr);
  
  jac = fabs(det)/6.0;

  cmulv (jac,contr);
  
  for (i=0;i<contr.n;i++){
    nfor[i]+=contr[i];
  }
  
}




////////////////////       /* 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 lintet::ntdbr_vector (long eid,vector &ntdbr)
{
  long i,j,ipp,ri,ci;
  double det,jac,w;
  vector x(nne),y(nne),z(nne),volcoord(4),eps(tncomp),sig(tncomp);
  matrix d(tncomp,tncomp);

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

  Mt->give_node_coord3d (x,y,z,eid);

  fillv (0.25,volcoord);
  w = 1.0/6.0;

  det = det3d (x.a,y.a,z.a);
  jac = w*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] * volcoord[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 lintet::ntn_matrix (long eid,matrix &ntn)
{
  long i;
  double det,jac;
  vector x(nne),y(nne),z(nne);
  
  Mt->give_node_coord3d (x,y,z,eid);
  det = det3d (x.a,y.a,z.a);
  jac = det/120.0;

  fillm (jac,ntn);

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

/**
   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 lintet::compute_error (long eid,double &volume,double &e2,double &u2,double &sizel,double *rsigfull)
{
  long intord = 4;
  long i,ipp,ri,ci;
  double det,jac,contr,zero=1.0e-20;
  ivector nodes (nne);
  vector x(nne),y(nne),z(nne),volcoord(4),gp1(intord),gp2(intord),gp3(intord),w(intord);
  vector sig_star(tncomp),sig_roof(tncomp),sig_err(tncomp),eps(tncomp);
  matrix d(tncomp,tncomp),dinv(tncomp,tncomp);
  
  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);

  det = det3d (x.a,y.a,z.a);

  Mm->matstiff (d,ipp);
  invm (d,dinv,zero);
  
  // compute u2
  fillv (0.25,volcoord);

  jac = det/6.0;

  Mm->givestrain (0,ipp,eps);

  mxv (d,eps,sig_roof);
  scprd (eps,sig_roof,contr);

  u2 = jac*contr;

  // compute e2
  gauss_points_tet (gp1.a,gp2.a,gp3.a,w.a,intord);

  e2 = 0;  
  for (i=0;i<intord;i++){
    volcoord[0]=gp1[i];
    volcoord[1]=gp2[i];
    volcoord[2]=gp3[i];
    volcoord[3]=1.0-gp1[i]-gp2[i]-gp3[i];

    jac=w[i]*det;

    give_der_star (volcoord,rsigfull,nodes,sig_star,Mt->nn);
    subv (sig_star,sig_roof,sig_err);
    
    vxmxv (sig_err,dinv,contr);
    e2 += jac*contr;
  }

  volume += det/6.0;
  sizel = pow (1.4142135624*det,1.0/3.0);
}

/**
   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 lintet::elchar (long eid,double *&spsig)
{
  long ipp,ri,ci;
  vector eps(tncomp);
  matrix d(tncomp,tncomp);
  
  spsig = new double[tncomp*1];
  
  ri = ci = 0;
  ipp = Mt->elements[eid].ipp[ri][ci];
  
  Mm->matstiff (d,ipp);
  Mm->givestrain (0,ipp,eps);
  mxv (d.a,eps.a,spsig,d.m,d.n);
}
////////////////////       /* termitovo */       ////////////////////////////////////