meshtransfer.cpp

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

  }
  
  
  delete [] nchilds;
  for (i=0;i<gt_old->ne;i++)  delete [] childip[i];
  delete [] childip;
  delete [] x;
  delete [] y;
  for (i=0;i<maxnchild;i++)  pointval1[i] = NULL;  delete [] pointval1;
  for (i=0;i<maxnchild;i++)  pointval2[i] = NULL;  delete [] pointval2;
  for (i=0;i<maxnchild;i++)  pointval3[i] = NULL;  delete [] pointval3;
}

/**
   Function transforms values from int. points of old mesh into int. points of new mesh.
   
   @param gt_old - gtopology of old mesh
   @param mm_new - mechmat of new mesh
   @param apcoord - array of coordinates of new mesh int. points, size is mm_new->tnip x dim
   @param parentel - array of parent elements of new mesh int. points, size is mt_new->tnip
   @param nipcomp - number of members (of arrays Mm->ip[].strain , Mm->ip[].stress and Mm->ip[].eqother), which are approximated;
                    nipcomp[0] <0;Mm->ip[].ncompstr>, nipcomp[1] <0;Mm->ip[].ncompstr>, nipcomp[2] <0;Mm->ip[].ncompother>
		    if nipcomp = -1 => all members
   @param mm_old - mechmat of old mesh
   @param dim - dimension of solved problem
   
   created  7.12.2002, Ladislav Svoboda, [e-mail]
 */
void transfvalues_ip_direct (gtopology *gt_old,mechmat *mm_new,const double **apcoord,const long *parentel,long nipcomp[],const mechmat *mm_old,long dim)
{
  long i,j,nvals;
  double **spvalue,**apvalue;
  
  if (nipcomp[0]==-1 || nipcomp[0] > Mm->ip[0].ncompstr  ) nipcomp[0] = Mm->ip[0].ncompstr;
  if (nipcomp[1]==-1 || nipcomp[1] > Mm->ip[0].ncompstr  ) nipcomp[1] = Mm->ip[0].ncompstr;
  if (nipcomp[2]==-1 || nipcomp[2] > Mm->ip[0].ncompother) nipcomp[2] = Mm->ip[0].ncompother;
  
  nvals = nipcomp[0] + nipcomp[1] + nipcomp[2];
  
  spvalue = new double* [gt_old->ne];
  
  for (i=0;i<gt_old->ne;i++){
    switch (gt_old->gelements[i].auxinf){
    case 312:{
      spvalue[i] = new double [1*nvals];
      
      for (j=0;j<nipcomp[0];j++)  spvalue[i][j                      ] = mm_old->ip[i].strain[j];
      for (j=0;j<nipcomp[1];j++)  spvalue[i][j+nipcomp[0]           ] = mm_old->ip[i].stress[j];
      for (j=0;j<nipcomp[2];j++)  spvalue[i][j+nipcomp[0]+nipcomp[1]] = mm_old->ip[i].eqother[j];
      
      break;
    }
    default:{
	fprintf (stderr,"\n\n wrong dimdegnne in function transfvalues_ip_direct (%s, line %d)",__FILE__,__LINE__);
	break;
    }}
  }
  
  //jede jen pro 312 !!!!!!!!!!!!!!!!!!!!
  apvalue = new double* [mm_new->tnip];
  for (i=0;i<mm_new->tnip;i++)
    apvalue[i] = mm_new->ip[i].eqother;    // !!! zde to nereflektuje styl nstra x nstre x nother
  
  least_square spr(dim,nvals,gt_old,((Mp->adaptflag & 16) ? 16 : 0));
  spr.L2_sp2sp (Out,spvalue,mm_new->tnip,parentel,apcoord,apvalue);
  
  
  destrm (spvalue,gt_old->ne);
  for (i=0;i<mm_new->tnip;i++) apvalue[i] = NULL;
  delete [] apvalue;
}

/**
   Function transforms displacement from nodes of old mesh into nodes of new mesh by direct interpolation.
   
   @param po - problem of old mesh
   @param pn - problem of new mesh
   @param lcid - id of load case
   @param ndofn - number of DOFs in node
   @param parentel - array of parent (old)elements of new mesh nodes, size is p_new->mt->nn
   
   created  7.12.2002, Ladislav Svoboda, [e-mail]
 */
//void transfvalues_nod (problem *po,problem *pn,long lcid,long ndofn,long *parentel)


/**
   Function transforms values from nodes of old mesh into nodes of new mesh.
   
   @param po - problem of old mesh
   @param pn - problem of new mesh
   @param lcid - id of load case
   @param dim - dimension of solved problem
   @param ndofn - number of DOFs in node
   @param parentel - array of parent (old)elements of new mesh nodes, size is p_new->mt->nn
   
   created  7.12.2002, Ladislav Svoboda, [e-mail]
 */
void transfvalues_nod (problem *po,problem *pn,long lcid,long dim,long ndofn,long *parentel,char spr)
{
  long i,j,k,nid,maxnchild;
  long *nchilds,**childnod;
  double **spcoord, *r_old, **r_new;
  vector x,y;
  
  // array of deformations of old mesh
  r_old = new double [ndofn * po->gt->nn];
  give_rfull (lcid,po,r_old);
  
  // allocation
  nchilds = new long [po->gt->ne];
  memset (nchilds,0,po->gt->ne*sizeof(long));
  
  for (i=0;i<pn->gt->nn;i++)
    nchilds[parentel[i]]++;
  
  childnod = new long* [po->gt->ne];
  r_new = new double* [po->gt->ne];
  for (i=0;i<po->gt->ne;i++){
    childnod[i] = new long [nchilds[i]];
    r_new[i] = new double [nchilds[i]*ndofn];
  }
  
  if (spr=='y'){
    spcoord = new double* [po->gt->ne];
    for (i=0;i<po->gt->ne;i++)  spcoord[i] = new double [nchilds[i]*dim];
  }
  else if (spr=='n'){
    maxnchild = 0;
    for (i=0;i<po->gt->ne;i++)  if (maxnchild<nchilds[i]) maxnchild = nchilds[i];
    allocv (maxnchild,x);  allocv (maxnchild,y);
  }
  else fprintf (stderr,"\n\n unknown spr flag is required in function transfvalues_nod (%s, line %d).\n",__FILE__,__LINE__);
  
  // transformation from r_old to r_new
  memset (nchilds,0,po->gt->ne*sizeof(long));
  for (i=0;i<pn->gt->nn;i++)
    childnod[parentel[i]][nchilds[parentel[i]]++] = i;
  
  // ********************
  /*
  double X,Y,f;
  
  for (i=0;i<po->mt->nn;i++){
    X = po->gt->gnodes[i].x;
    Y = po->gt->gnodes[i].y;
    f = -X*X*X*X/1000.0 - X*X*X/10.0 + X*X + 10*X -Y*Y*Y*Y/3.0 - Y*Y*Y/2.0 + 20*Y*Y + 30*Y - 212;
    
    r_old[i           ] = f;
    r_old[i+po->mt->nn] = f;
  }
  
  long hod,min;
  double sec = clock();
  */
  // ********************
  
  if (spr=='y'){
    for (i=0;i<po->gt->ne;i++){
      for (j=0;j<nchilds[i];j++){
	spcoord[i][j*dim  ] = pn->gt->gnodes[childnod[i][j]].x;
	spcoord[i][j*dim+1] = pn->gt->gnodes[childnod[i][j]].y;
	if (dim==3) spcoord[i][j*dim+2] = pn->gt->gnodes[childnod[i][j]].z;
      }
    }
    
    least_square spr(dim,ndofn,po->gt,((Mp->adaptflag & 16) ? 16 : 0));
    spr.L2_nod2sp (Out,nchilds,spcoord,r_old,r_new,'n');
  }
  else if (spr=='n'){
    for (i=0;i<po->gt->ne;i++){
      for (j=0;j<nchilds[i];j++){
	x[j] = pn->gt->gnodes[childnod[i][j]].x;
	y[j] = pn->gt->gnodes[childnod[i][j]].y;
      }      
      
      give_valuesinpoints (po->gt,i,nchilds[i],x.a,y.a,ndofn,r_old,&r_new[i],'1');
    }
  }
  else fprintf (stderr,"\n\n unknown spr flag is required in function transfvalues_nod (%s, line %d).\n",__FILE__,__LINE__);
  
  // ********************
  /*
  sec = (clock() - sec) / (double)CLOCKS_PER_SEC;
  hod = (long)sec/3600;  sec -= hod*3600;
  min = (long)sec/60;    sec -= min*60;
  fprintf (stdout,"\n -----------------------------------");
  fprintf (stdout,"\n Consumed time by TRANSF %2ld:%02ld:%05.2f",hod,min,sec);
  fprintf (stdout,"\n -----------------------------------\n");
  
  double a,ai,nordr,norr=0.0;
  vector dr(2*pn->gt->nn),r(pn->gt->nn);
  
  for (i=0;i<po->gt->ne;i++)
    for (j=0;j<nchilds[i];j++){
      nid = childnod[i][j];
      X = pn->gt->gnodes[nid].x;
      Y = pn->gt->gnodes[nid].y;
      r[nid] = -X*X*X*X/1000.0 - X*X*X/10.0 + X*X + 10*X -Y*Y*Y*Y/3.0 - Y*Y*Y/2.0 + 20*Y*Y + 30*Y - 212;
      
      dr[nid+pn->gt->nn] = r[nid] - r_new[i][j*2+1];
      dr[nid           ] = r[nid] - r_new[i][j*2  ];
    }
  
  sizev (dr,nordr);
  fprintf (stdout,"\n\n nordr     =  %25.15f \n\n",nordr);
  fprintf (stdout,"\n\n nordr/nn1 =  %25.15f \n\n",nordr/sqrt(769));
  fprintf (stdout,"\n\n nordr/nn2 =  %25.15f \n\n",nordr/sqrt(24048));
  
  nordr = a = 0.0;
  for (i=0;i<pn->gt->ne;i++){
    nordr += Pelt->error(i,dr,ai);
    a += ai;
  }
  
  nordr = sqrt(nordr/a);
  fprintf (stdout,"\n\n nordr(A)  =  %25.15f \n\n",nordr);
  
  exit (1);
  
  
  for (i=0;i<pn->gt->nn;i++)  norr += r[i];
  sizev (r,norr); norr *= 2;
  */
  // ********************
  
  // saving deformations from r_old to p_new->lsrs->lhs
  for (i=0;i<po->gt->ne;i++)
    for (j=0;j<nchilds[i];j++){
      nid = childnod[i][j];
      for (k=0;k<ndofn;k++)
	if (pn->gt->gnodes[nid].cn[k] > 0)
	  pn->lsrs->lhs[lcid*pn->lsrs->ndof + pn->gt->gnodes[nid].cn[k] - 1] = r_new[i][j*ndofn+k];
    }
  
  delete [] r_old;
  delete [] nchilds;
  destrm (childnod,po->gt->ne);
  destrm (r_new,po->gt->ne);
  if (spr=='y')  destrm (spcoord,po->gt->ne);
}

/**
   Function approximates 'values' from nodes into 'points'.
   All the point lays in 'element', which is determined by 'gt' and 'eid'.
   Values in points are returned by array pointval.
   
   @param npoints - number of points
   @param xx - array of x-coordinates of points, size is npoints
   @param yy - array of y-coordinates of points, size is npoints
   @param nval - number of values
   @param nodvalues - 1D array of values in nodes, size is nval*gt->nn, value position in array is
                      (val[0] on nod[0] ... val[0] on nod[gt->nn] ... val[nval] on nod[0] ... val[nval] on nod[gt->nn])
   @param pointval - answer = 2D array, size is npoints*nval
   
   created  7.12.2002, Ladislav Svoboda, [e-mail]
 */
void give_valuesinpoints (gtopology *gt,long eid,long npoints,double *xx,double *yy,long nval,double *nodvalues,double **pointval,char flag)
{
  long i,j,k,nne;
  double xi,eta,value;
  vector x,y,nodval;
  ivector nodes;
  
  nne = gt->give_nne (eid);
  
  allocv (nne,x);
  allocv (nne,y);
  allocv (nne,nodval);
  allocv (nne,nodes);
  
  gt->give_nodes (eid,nodes);
  gt->give_node_coord2d (x,y,eid);
  
  for (i=0;i<npoints;i++){
    
    switch (gt->gelements[eid].auxinf){
    case 312:{ nc_lin_3_2d          (xx[i],yy[i],x.a,y.a,xi,eta); break; }
    case 622:{ nc_quad_3_2d (0.00000001,xx[i],yy[i],x.a,y.a,xi,eta); break; }
    case 412:{ nc_lin_4_2d  (0.00001,xx[i],yy[i],x.a,y.a,xi,eta); break; }
    case 822:{ nc_quad_4_2d (0.00001,xx[i],yy[i],x.a,y.a,xi,eta); break; }
    default:{
      fprintf (stderr,"\n\n unknown nnedegdim is required in function give_valuesinpoints (%s, line %d).\n",__FILE__,__LINE__);
    }}
    
    for (j=0;j<nval;j++){
      for (k=0;k<nne;k++)
	nodval[k] = nodvalues[nodes[k]+j*gt->nn];
      
      switch (gt->gelements[eid].auxinf){
      case 312:{ value = Pelt->approx_nat (xi,eta,nodval);  break; }
      case 622:{ value = Peqt->approx     (xi,eta,nodval);  break; }
      case 412:{ value = Pelq->approx     (xi,eta,nodval);  break; }
      case 822:{ value = Peqq->approx     (xi,eta,nodval);  break; }
      default:{
	fprintf (stderr,"\n\n unknown nnedegdim is required in function give_valuesinpoints (%s, line %d).\n",__FILE__,__LINE__);
      }}
      
      if (flag=='2')
	pointval[i][j] = value;
      else if (flag=='1')
	pointval[0][i*nval+j] = value;
      else ;
    }
  }
}

/**
   function finds out node(from gt), which is the closest by required point
   
   @param x,y,z - coordinates of point
   
   created  30.1.2003, Ladislav Svoboda, [e-mail]
 */
long adjacnode (gtopology *gt,double x,double y,double z)
{
  long i,a = 0;
  double d,min = 1e50;
  
  for(i=0;i<gt->ne;i++){
    d = (gt->gnodes[i].x - x)*(gt->gnodes[i].x - x) + (gt->gnodes[i].y - y)*(gt->gnodes[i].y - y) + (gt->gnodes[i].z - z)*(gt->gnodes[i].z - z);
    if (d<min){
      min=d;
      a=i;
    }
  }
  return (a);
}
////////////////////       /* termitovo */       ////////////////////////////////////