soilplateq.cpp

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

    allocv (ndofe,aux);
    allocm (ndofe,ndofe,tmat);
    transfmat (nodes,tmat);
    //locglobtransf (aux,r,tmat);
    lgvectortransf (aux,r,tmat);
    copyv (aux,r);
    destrv (aux);
    destrm (tmat);
  }

  mainip_strains (lcid,eid,0,0,x,y,r);
}


void soilplateq::mainip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
{
  long i,j,ii,ipp;
  double a;
  vector gp,w,eps,natcoord(2),l(3);
  matrix gm,atd(16,12);

  // translation to center it is necesary for w which in local coord. and no in relation
  a=(x[0]+x[1]+x[2]+x[3])/4.;
  x[0]=x[0]-a;
  x[1]=x[1]-a;
  x[2]=x[2]-a;
  x[3]=x[3]-a;
  a=(y[0]+y[1]+y[2]+y[3])/4.;
  y[0]=y[0]-a;
  y[1]=y[1]-a;
  y[2]=y[2]-a;
  y[3]=y[3]-a;
  atd_matrix (atd,x,y);
  
//  fprintf (Out,"\n\n deformace");  fprintf (Out,"\n\n prvek cislo %ld",eid);
 for (ii=0;ii<nb;ii++){
    ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
    allocv (intordsm[ii][ii],gp);  allocv (intordsm[ii][ii],w);
    allocv (ncomp[ii],eps);
    allocm (ncomp[ii],ndofe,gm);
    gauss_points (gp.a,w.a,intordsm[ii][ii]);

    for (i=0;i<intordsm[ii][ii];i++){
        l[0]=gp[i];
		for (j=0;j<intordsm[ii][ii];j++){
  		  l[1]=gp[j];
		  if(ii==0) geom_matrix (gm,atd,x,y,l);
		  else if(ii==1) geom_matrix_shear (gm,atd,x,y,l);
          mxv (gm,r,eps);
          Mm->storestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
          ipp++;
//          fprintf (Out,"\n"); for (ij=0;ij<ncomp[ii];ij++){fprintf (Out,"%20.10e",eps[ij]);}
		}
	}
    destrm (gm);  destrv (eps);  destrv (gp);
  }
  
}

void soilplateq::elem_strains (double **stra,long lcid,long eid,long ri,long ci)
{
  long i,j,ii,ipp;
  double *lsm,*lhs,*rhs;
  vector nxi(nne),neta(nne),gp,w,eps,natcoord(2),l(3);
  ivector cn,nodes(nne);

  lsm = new double [9];
  nullv (lsm,9);

  nodecoord (nxi,neta);
  Mt->give_elemnodes (eid,nodes);
  

//  fprintf (Out,"\n\n deformace prvek cislo %ld",eid);
  for (ii=0;ii<nb;ii++){
    ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
    allocv (ncomp[ii],eps);
    allocv (intordsm[ii][ii],gp);    allocv (intordsm[ii][ii],w);
    lhs = new double [ncomp[ii]*3];
    rhs = new double [ncomp[ii]*3];
    gauss_points (gp.a,w.a,intordsm[ii][ii]);
    nullv (rhs,ncomp[ii]*3);

    for (i=0;i<intordsm[ii][ii];i++){
        l[0]=gp[i];
		for (j=0;j<intordsm[ii][ii];j++){
  		  l[1]=gp[j];
	      Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
          natcoord[0]=l[0];  natcoord[1]=l[1];
          matassem_lsm (lsm,natcoord);
          rhsassem_lsm (rhs,natcoord,eps);
//          fprintf (Out,"\n");
//		  for (long ij=0;ij<ncomp[ii];ij++){
//		    fprintf (Out,"%20.10e",eps[ij]);
//			}
        ipp++;
		}
	}
    solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
    nodal_values (stra,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii]);
    delete [] lhs;  delete [] rhs;
    destrv (eps);  destrv (w);  destrv (gp);
  }
  delete [] lsm;
}

void soilplateq::nod_strains (long lcid,long eid,long ri,long ci)
{
  long i,j,ii,ipp;
  double *lsm,*lhs,*rhs;
  vector nxi(nne),neta(nne),gp,w,eps,natcoord(2),l(3);
  ivector cn,nodes(nne);

  lsm = new double [9];

  nodecoord (nxi,neta);
  Mt->give_elemnodes (eid,nodes);
  

//  fprintf (Out,"\n\n deformace");
//  fprintf (Out,"\n\n prvek cislo %ld",eid);
  for (ii=0;ii<nb;ii++){
    ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
    allocv (ncomp[ii],eps);
    allocv (intordsm[ii][ii],gp);    allocv (intordsm[ii][ii],w);
    lhs = new double [ncomp[ii]*3];
    rhs = new double [ncomp[ii]*3];
    gauss_points (gp.a,w.a,intordsm[ii][ii]);
    nullv (lsm,9);
    nullv (rhs,ncomp[ii]*3);

    for (i=0;i<intordsm[ii][ii];i++){
        l[0]=gp[i];
		for (j=0;j<intordsm[ii][ii];j++){
  		  l[1]=gp[j];
	      Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
          natcoord[0]=l[0];  natcoord[1]=l[1];
          matassem_lsm (lsm,natcoord);
          rhsassem_lsm (rhs,natcoord,eps);
//          fprintf (Out,"\n");
//		  for (ij=0;ij<ncomp[ii];ij++){
//		    fprintf (Out,"%20.10e",eps[ij]);
//			}
        ipp++;
		}
	}
    solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
    Mt->strain_nodal_values (nodes,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii],lcid);
    delete [] lhs;  delete [] rhs;
    destrv (eps);  destrv (w);  destrv (gp);
  }
  
  delete [] lsm;
}

void soilplateq::appstrain (long lcid,long eid,double xi,double eta,long fi,long ncomp,vector &eps)
{
  long i,j,k;
  ivector nodes;
  vector l(3),nodval;
  if (ncomp != eps.n){
    fprintf (stderr,"\n\n wrong interval of indices in function strain (%s, line %d).\n",__FILE__,__LINE__);
    abort ();
  }


  allocv (nne,nodes);
  allocv (nne,nodval);
  Mt->give_elemnodes (eid,nodes);
  k=0;
  for (i=fi;i<fi+ncomp;i++){
    for (j=0;j<nne;j++){
      nodval[j]=Mt->nodes[nodes[j]].strain[lcid*tncomp+i];
    }
    eps[k]=approx (xi,eta,nodval);
    k++;
  }
  
  destrv (nodes);  destrv (nodval);
}

void soilplateq::allip_strains (double **stra,long lcid,long eid,long ri,long ci)
{
  long i,j,ii,jj,ipp;
  vector eps(tncomp),gp,w,areacoord(3);
  
  for (ii=0;ii<nb;ii++){
    for (jj=0;jj<nb;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]);
      ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
      
      for (i=0;i<intordsm[ii][jj];i++){
    	for (j=0;j<intordsm[ii][jj];j++){
    	if (Mp->strainaver==0)
		  appval (gp[i],gp[j],0,tncomp,eps,stra);
	    if (Mp->strainaver==1)
	      appstrain (lcid,eid,gp[i],gp[j],0,tncomp,eps);
          Mm->storestrain (lcid,ipp,eps);
	      ipp++;
		}
      }
      destrv (w);  destrv (gp);
    }
  }
}

void soilplateq::strains (long lcid,long eid,long ri,long ci)
{
  long i,naep,ncp,sid;
  double **stra;
  vector coord,eps;
  
  if (Mp->strainaver==0){
    stra = new double* [nne];
    for (i=0;i<nne;i++){
      stra[i] = new double [tncomp];
    }
    elem_strains (stra,lcid,eid,ri,ci);
  }
  
  switch (Mm->stra.tape[eid]){
  case nowhere:{
    break;
  }
  case intpts:{
    allip_strains (stra,lcid,eid,ri,ci);
    break;
  }
  case enodes:{
    break;
  }
  case userdefined:{
    //  number of auxiliary element points
    naep = Mm->stra.give_naep (eid);
    ncp = Mm->stra.give_ncomp (eid);
    sid = Mm->stra.give_sid (eid);
    allocv (ncp,eps);
    allocv (2,coord);
    for (i=0;i<naep;i++){
      Mm->stra.give_aepcoord (sid,i,coord);
      if (Mp->strainaver==0)
		appval (coord[0],coord[1],0,ncp,eps,stra);
      if (Mp->strainaver==1)
        appstrain (lcid,eid,coord[0],coord[1],0,ncp,eps);
      Mm->stra.storevalues(lcid,eid,i,eps);
    }
    destrv (eps);
    destrv (coord);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown strain point is required in function dstelem::strains (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
  if (Mp->strainaver==0){
    for (i=0;i<nne;i++){
      delete [] stra[i];
    }
    delete [] stra;
  }

}


void soilplateq::res_allip_stresses (long lcid,long eid)
{
   allip_stresses (lcid, eid, 0, 0);
}

/**
   function computes stresess
   
   @param lcid - load case id
   @param eid - element id
   @param ri,ci - row and column indices
   
   20.9.2006
*/
void soilplateq::allip_stresses (long lcid,long eid,long ri,long ci)
{
  long ipp,i,j;
  ivector nodes(nne);
  vector eps, sig;
  matrix cc;
  
  Mt->give_elemnodes (eid,nodes);
    
  // translation to center it is necesary for w which in local coord. and no in relation
//  elem_strains (stra,lcid,eid,ri,ci); 
  
  allocv (ncomp[0],eps);  allocv (ncomp[0],sig);
  allocm (ncomp[0],ncomp[0],cc);
        
  ipp=Mt->elements[eid].ipp[ri][ci];
 
//  fprintf (Out,"\n\n Eps,sig soilprvek cislo %ld\n",eid);
  for (i=0;i<intordsm[0][0];i++){
	for (j=0;j<intordsm[0][0];j++){
      Mm->matstiff (cc,ipp);
	  Mm->givestrain (lcid,ipp,cncomp[0],ncomp[0],eps);
	  mxv (cc,eps,sig);	
	  Mm->storestress (lcid,ipp,cncomp[0],sig);      
//  for (k=0;k<3;k++){
//    printf (Out," %15.5e",eps[k]); fprintf (Out," %15.5e",sig[k]);  }
//    fprintf (Out,"\n");
	  ipp++;  
	}
  }
  destrv (eps);  destrv (sig); destrm (cc); 

//  fprintf (Out,"\n\n Fint soil prvek cislo %ld\n",eid);
//  for (k=0;k<ndofe;k++){fprintf (Out," %15.5e",ifor[k]);}
  
}
void soilplateq::res_internal_forces (long lcid,long eid,vector &ifor)
{
   internal_forces (lcid, eid, 0, 0, ifor);
}

/**
   function computes internal forces
   
   @param lcid - load case id
   @param eid - element id
   @param ri,ci - row and column indices
   @param ifor - vector of internal forces
   
   20.12.2002
*/
void soilplateq::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
{
  long ipp,ii,i,j,k;
  double **stra;
  double a,a0,a1,jac;
  ivector nodes(nne);
  vector x(nne),y(nne),z(nne),f(ndofe),eps(3),sig(3);
  vector w,gp,l(2);
  matrix gm,cc,tran(3,3),atd(16,12);
  
  stra = new double* [nne];
  for (i=0;i<nne;i++){
    stra[i] = new double [tncomp];
  }
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);
    
  // translation to center it is necesary for w which in local coord. and no in relation
  a=(x[0]+x[1]+x[2]+x[3])/4.;
  x[0]=x[0]-a;
  x[1]=x[1]-a;
  x[2]=x[2]-a;
  x[3]=x[3]-a;
  a=(y[0]+y[1]+y[2]+y[3])/4.;
  y[0]=y[0]-a;
  y[1]=y[1]-a;
  y[2]=y[2]-a;
  y[3]=y[3]-a;
  //     jakobian
  a = (x[2]-x[0])*(y[3]-y[1])-(y[2]-y[0])*(x[3]-x[1]);
  a0=-(x[3]-x[2])*(y[1]-y[0])+(y[3]-y[2])*(x[1]-x[0]);
  a1=-(x[2]-x[1])*(y[3]-y[0])+(y[2]-y[1])*(x[3]-x[0]);
  atd_matrix (atd,x,y);
  fillv (0.0,ifor);
//  elem_strains (stra,lcid,eid,ri,ci); 
 

// fprintf (Out,"\n\n Eps,sig soilprvek cislo %ld\n",eid);
 for (ii=0;ii<nb;ii++){
  ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
  allocv (intordsm[ii][ii],gp);  allocv (intordsm[ii][ii],w);
  allocv (ncomp[ii],eps);
  allocm (ncomp[ii],ndofe,gm);
  gauss_points (gp.a,w.a,intordsm[ii][ii]);
  allocm (ncomp[ii],ncomp[ii],cc); 
  for (i=0;i<intordsm[0][0];i++){
	for (j=0;j<intordsm[0][0];j++){
	  l[0]=gp[i];  l[1]=gp[j];
	  jac=(a+a0*l[0]+a1*l[1])/8.*w[i]*w[j];
//	  appval (l[0],l[1],0,tncomp,eps,stra);
//	  Mm->storestrain (lcid,ipp,eps);	
//	  Mm->givestrain (lcid,ipp,cncomp[0],ncomp[0],eps);
      Mm->computenlstresses (ipp);
      Mm->givestress (lcid,ipp,eps);  //strain for consolidation
      Mm->matstiff (cc,ipp);
	  mtxv (cc,eps,sig);	
	  geom_matrix (gm,atd,x,y,l);
	  mtxv (gm,sig,f);	
	  for (k=0;k<f.n;k++){ ifor[k]+=f[k]*jac;}	
	  ipp++;  
//   for (k=0;k<3;k++){fprintf (Out," %15.5e",eps[k]); fprintf (Out," %15.5e",sig[k]);  }   fprintf (Out,"\n");
	}
  }
  destrm (cc); destrm (gm);  destrv (gp);  destrv (w);
 }
//  fprintf (Out,"\n\n Fint soil prvek cislo %ld\n",eid);
//  for (k=0;k<ndofe;k++){fprintf (Out," %15.5e",ifor[k]);}
  for (i=0;i<nne;i++){ delete [] stra[i];  }
  delete [] stra;    
}