meshtransfer.cpp

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

   @param ipcoord - array of coordinates of new mesh int. points, size is mm_new->tnip x dim
   @param passedel1 - empty array, size is gt_old->ne
   @param passedel2 - zero array, size is gt_new->ne
   @param nelheap = 1
   @param elheap - array of size gt_new->ne, elheap[0] = 0
   @param newelheap - empty array, size is gt_new->ne
   @param susel - empty array, size is gt_old->ne
   @param newsusel - answer = 1D array, size is gt_old->ne
   @param parentel - answer = 1D array, size is mm_new->tnip
   
   created  3.12.2002, Ladislav Svoboda, [e-mail]
 */
void locintpoints (gtopology *gt_old,gtopology *gt_new,mechtop *mt_new,double **ipcoord,long *passedel1,long *passedel2,long nelheap,long *elheap,long *newelheap,long *susel,long *newsusel,long *parentel)
{
  long i,j,k,l,el,nip,ipp,nnewelheap,nsusel;
  
  nnewelheap = 0;
  for (i=0;i<nelheap;i++){
    nsusel = mt_new->give_nip(elheap[i],0,0);
    ipp = mt_new->elements[elheap[i]].ipp[0][0];
    for (j=0;j<nsusel;j++)
      susel[j] = parentel[ipp++];
    
    for (j=0;j<gt_new->nadjelnod[elheap[i]];j++){
      el = gt_new->adjelel[elheap[i]][j];
      
      if (!passedel2[el]){
	nip = 0;
	ipp = mt_new->elements[el].nb;
	for (k=0;k<ipp;k++)
	  for (l=0;l<ipp;l++)
	    nip += mt_new->elements[el].nip[k][l];
	
	ipp = mt_new->elements[el].ipp[0][0];
	for (k=0;k<nip;k++){
	  for (l=0;l<gt_old->ne;l++)
	    passedel1[l] = 0;
	  
	  parentel[ipp+k] = whereispoint (gt_old,nsusel,susel,newsusel,passedel1,ipcoord[ipp+k][0],ipcoord[ipp+k][1],'f');
	}
	
	passedel2[el] = 1;
	newelheap[nnewelheap++] = el;
      }
    }
  }
  
  if (!nnewelheap)
    return ;
  
  locintpoints (gt_old,gt_new,mt_new,ipcoord,passedel1,passedel2,nnewelheap,newelheap,elheap,susel,newsusel,parentel);
  return ;
}

/**
   Function runs function 'locnodes'.
   
   @param gt_old - gtopology of old mesh
   @param gt_new - gtopology of new mesh
   @param parentel_nod - answer = array of parent elements, size is mt_new->nn
   
   created  6.12.2002, Ladislav Svoboda, [e-mail]
 */
void findout_parentel_nod (gtopology *gt_old,gtopology *gt_new,long *parentel_nod,long dim)
{
  locnodes (gt_old,gt_new,parentel_nod,dim);
}

/**
   For every node (from new mesh) function finds out "parent" element (from old mesh),
   in which the "new" node lays.
   
   @param gt_old - gtopology of old mesh
   @param gt_new - gtopology of new mesh
   @param parentel - answer = allocated array of parent elements, size is mt_new->nn
   
   created  4.12.2002, Ladislav Svoboda, [e-mail]
*/
void locnodes (gtopology *gt_old,gtopology *gt_new,long *parentel,long dim)
{
  long i,j,k,n,nod,endheap;
  ivector oldsusel(gt_old->ne),newsusel(gt_old->ne),passedel(gt_old->ne),nodheap(gt_new->nn),nsusel(gt_new->nn);
  imatrix susel;
  
  if (dim==2)  allocm (gt_new->nn,18,susel);
  if (dim==3)  allocm (gt_new->nn,120,susel);
  
  fillv (-1,parentel,gt_new->nn);
  
  endheap = 1;
  nsusel[0] = 1;
  susel[0][0] = 0;
  nodheap[0] = 0;
  
  for (i=0;i<gt_new->nn;i++){
    nod = nodheap[i];
    
    nullv (passedel);
    copyv (susel[nod],oldsusel.a,nsusel[nod]);
    
    parentel[nod] = whereispoint (gt_old,nsusel[nod],oldsusel.a,newsusel.a,passedel.a,gt_new->gnodes[nod].x,gt_new->gnodes[nod].y,'f');
    
    for (j=0;j<gt_new->nadjelnod[nod];j++)
      for (k=0;k<gt_new->gelements[gt_new->adjelnod[nod][j]].nne;k++){
	n = gt_new->gelements[gt_new->adjelnod[nod][j]].nodes[k];
	if (parentel[n] == -1){
	  parentel[n] = -2;
	  nodheap[endheap++] = n;
	}
	if (parentel[n] == -2)
	  susel[n][nsusel[n]++] = parentel[nod];
      }
    
  }
}

/**
   Function finds out element including "point".
   First of all, the elements in "susel" are investigated.
   
   @param gt - gtopology, where the element is finding out
   @param nsusel - number suspected elements
   @param susel - array of suspected elements, size is gt->ne
   @param newsusel - empty array, size is gt->ne
   @param passedel - zero !!!!! array, size is gt->ne
   @param x,y - coordinates of point

   created  3.12.2002, Ladislav Svoboda, [e-mail]
 */
long whereispoint (gtopology *gt,long nsusel,long *susel,long *newsusel,long *passedel,double x,double y,char flag)
{
  long i,j,el,nnewsusel;
  
  if (flag=='f')
    for (i=0;i<nsusel;i++)
      if ( !passedel[susel[i]]++ && (ispointinel (gt->gnodes,gt->gelements[susel[i]].nodes,x,y,gt->gelements[susel[i]].auxinf)) )
	return (susel[i]);
  
  nnewsusel = 0;
  for (i=0;i<nsusel;i++){
    for (j=0;j<gt->nadjelnod[susel[i]];j++){
      el = gt->adjelel[susel[i]][j];
      if (!passedel[el]){
	if (ispointinel (gt->gnodes,gt->gelements[el].nodes,x,y,gt->gelements[el].auxinf))
	  return (el);
	
	passedel[el] = 1;
	newsusel[nnewsusel++] = el;
      }
    }
  }
  
  if (!nnewsusel){
    fprintf (stderr,"\n\n!!!!!!!!!  Point does not lay in domain (%s, line %d)\n\n",__FILE__,__LINE__);
    return (-1);
  }
  
  return ( whereispoint (gt,nnewsusel,newsusel,susel,passedel,x,y,'n') );
}

/**
   Function finds out, whether 'node' lays in 'element'.
   Node is defined by its coordinates - x,y.
   Element is defined by its node numbers - 'elnod' and dimdegnne - 'ndd'
   
   @param nod - array of all gnodes
   @param elnod - array of node numbers of element
   @param x,y - coordinates of node
   @param ndd - nnedimdeg of element
   
   created  3.12.2002, Ladislav Svoboda, [e-mail]
 */
long ispointinel (gnode *nod,long *elnod,double x,double y,long ndd)
{
  switch (ndd){
  case 312:{
    if ( isnodonlhsofline (nod[elnod[0]],nod[elnod[1]],x,y) &&
	 isnodonlhsofline (nod[elnod[1]],nod[elnod[2]],x,y) &&
	 isnodonlhsofline (nod[elnod[2]],nod[elnod[0]],x,y)     )
      return (1);
    break;
  }
  case 412:{
    if ( isnodonlhsofline (nod[elnod[0]],nod[elnod[1]],x,y) &&
	 isnodonlhsofline (nod[elnod[1]],nod[elnod[2]],x,y) &&
	 isnodonlhsofline (nod[elnod[2]],nod[elnod[3]],x,y) &&
	 isnodonlhsofline (nod[elnod[3]],nod[elnod[0]],x,y)     )
      return (1);
    break;
  }
  case 622:{
    if ( isnodonlhsof3pcurve (nod[elnod[0]],nod[elnod[3]],nod[elnod[1]],x,y) &&
	 isnodonlhsof3pcurve (nod[elnod[1]],nod[elnod[4]],nod[elnod[2]],x,y) &&
	 isnodonlhsof3pcurve (nod[elnod[2]],nod[elnod[5]],nod[elnod[0]],x,y)     )
      return (1);
    break;
  }
  case 822:{
    if ( isnodonlhsof3pcurve (nod[elnod[0]],nod[elnod[4]],nod[elnod[1]],x,y) &&
	 isnodonlhsof3pcurve (nod[elnod[1]],nod[elnod[5]],nod[elnod[2]],x,y) &&
	 isnodonlhsof3pcurve (nod[elnod[2]],nod[elnod[6]],nod[elnod[3]],x,y) &&
	 isnodonlhsof3pcurve (nod[elnod[3]],nod[elnod[7]],nod[elnod[0]],x,y)     )
      return (1);
    break;
  }
  //case 413:{ break; }
  default:{
    fprintf (stderr,"\n\n wrong dimdegnne in function ispointinel (%s, line %d)",__FILE__,__LINE__);
    break;
  }}
  
  return (0);
}

/**
   Function finds out, whether point (xx,yy) lays on the left hand side of 'line'.
   Line is defined by two nodes n1,n2.
   
   created  3.12.2002, Ladislav Svoboda, [e-mail]
 */
long isnodonlhsofline (gnode &n1,gnode &n2,double xx,double yy)
{
  if (n1.x != n2.x)
    if ( (n2.x>n1.x ? 1.0:-1.0)*((xx-n1.x)*(n1.y-n2.y)/(n1.x-n2.x)+(n1.y-yy)) <= 1.0e-13 )
      return (1);
    else
      return (0);
  else
    if ( (n1.y-n2.y)*(xx-n1.x) > -1.0e-13 )
      return (1);
    else
      return (0);
}

/**
   Function finds out, whether point (xx,yy) lays on the left hand side of 'curve'.
   Curve is defined by three nodes n1,n2,n3.
   
   created  3.12.2002, Ladislav Svoboda, [e-mail]
 */
long isnodonlhsof3pcurve (gnode &n1,gnode &n2,gnode &n3,double xx,double yy)
{
  double c;

  c = (n1.x*(n3.y-n2.y)+n2.x*(n1.y-n3.y)+n3.x*(n2.y-n1.y))/(fabs(n1.x-n3.x)+fabs(n1.y-n3.y));
  if (-1.0e-13<c && c<1.0e-13)
    return (isnodonlhsofline (n1,n3,xx,yy));
  else{
    // docasne
    //fprintf (stdout,"\n\n\n hrana je krivocara \n\n\n");
    // docasne
    double s,x[3],y[2];
    
    c = sqrt((n3.x-n1.x)*(n3.x-n1.x)+(n3.y-n1.y)*(n3.y-n1.y));
    s = (n3.y-n1.y)/c;
    c = (n3.x-n1.x)/c;
    
    x[0] =  c*(  xx-n1.x) + s*(  yy-n1.y);  y[0] = -s*(  xx-n1.x) + c*(  yy-n1.y);
    x[1] =  c*(n2.x-n1.x) + s*(n2.y-n1.y);  y[1] = -s*(n2.x-n1.x) + c*(n2.y-n1.y);
    x[2] =  c*(n3.x-n1.x) + s*(n3.y-n1.y);
  
    if ( y[1]*x[0]*x[0]/x[1]/(x[1]-x[2]) + y[1]*x[2]*x[0]/x[1]/(x[2]-x[1]) -y[0] <= 1.0e-13 )
      return (1);
    
    return (0);
  }
}

/**
   Function transforms values from nodes of old mesh into int. points of new mesh.
   
   @param gt_old - gtopology of old mesh
   @param mm_new - mechmat of new mesh
   @param ipcoord - 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 nodvalue - 1D array of nodal values of old mesh, size is gt_old->nn * nvals
                     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])
   
   created  7.12.2002, Ladislav Svoboda, [e-mail]
 */
void transfvalues_ip_indirect (gtopology *gt_old,mechmat *mm_new,double **ipcoord,long *parentel,long nipcomp[],double *nodvalue)
{
  long i,j,maxnchild,nvals;
  long *nchilds,**childip;
  double *x,*y,**pointval1,**pointval2,**pointval3;
  
  
  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];
  
  nchilds = new long [gt_old->ne];
  memset (nchilds,0,gt_old->ne*sizeof(long));
  
  for (i=0;i<mm_new->tnip;i++)
    nchilds[parentel[i]]++;
  
  maxnchild = 0;
  childip = new long* [gt_old->ne];
  for (i=0;i<gt_old->ne;i++){
    childip[i] = new long [nchilds[i]];
    if (maxnchild<nchilds[i]) maxnchild = nchilds[i];
    nchilds[i] = 0;
  }
  
  for (i=0;i<mm_new->tnip;i++)
    childip[parentel[i]][nchilds[parentel[i]]++] = i;
  
  x = new double [maxnchild];
  y = new double [maxnchild];
  pointval1 = new double* [maxnchild];
  pointval2 = new double* [maxnchild];
  pointval3 = new double* [maxnchild];
  
  for (i=0;i<gt_old->ne;i++){
    for (j=0;j<nchilds[i];j++){
      x[j] = ipcoord[childip[i][j]][0];
      y[j] = ipcoord[childip[i][j]][1];
      
      if (nipcomp[0]) pointval1[j] = mm_new->ip[childip[i][j]].strain;
      if (nipcomp[1]) pointval2[j] = mm_new->ip[childip[i][j]].stress;
      if (nipcomp[2]) pointval3[j] = mm_new->ip[childip[i][j]].eqother;
    }
    
    if (nipcomp[0]) give_valuesinpoints (gt_old,i,nchilds[i],x,y,nipcomp[0],nodvalue                                     ,pointval1,'2');
    if (nipcomp[1]) give_valuesinpoints (gt_old,i,nchilds[i],x,y,nipcomp[1],nodvalue + gt_old->nn*(nipcomp[0]           ),pointval2,'2');
    if (nipcomp[2]) give_valuesinpoints (gt_old,i,nchilds[i],x,y,nipcomp[2],nodvalue + gt_old->nn*(nipcomp[0]+nipcomp[1]),pointval3,'2');