meshtransfer.cpp

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

#include "meshtransfer.h"

#include "global.h"
#include "globmat.h"
#include "intpoints.h"
#include "loadcase.h"
#include "element.h"
#include "plelemlt.h"
#include "plelemqt.h"
#include "plelemlq.h"
#include "plelemqq.h"
#include "nssolver.h"
#include "stochdriver.h"
#include "mefelinit.h"

#include "gtopology.h"
#include "gelement.h"
#include "gnode.h"
#include "gmatrix.h"
#include "mathem.h"
#include "problem.h"
#include "least_square.h"

#include <time.h>
#include <math.h>
#include <stdlib.h>

////////////////////       /* termitovo */       ////////////////////////////////////

/*
vse je delano pro 2D planestress prvky s ndofn=2
pro jine to neni kontrolovano !!!!!!!!!!!!!!!!!!
prvek plelemqq je predpokladan s rovnymi hranami(subparametricky), neb v fci 'nc_quad_4_2d' je zatim odkaz jen na 'nc_lin_4_2d'
*/
/*
ze site na sit se prenaseji posuny a prvnich ncompstr hodnot z pole other , ne prvnich ncompother neb gamma a hardening nepotrebuju
*/
/*
fce nezavisle na jine globalni promenne:
give_ndofn             ok
give_dof               ok
give_nip               ok
give_nne               ok
gt->give_nodes         ok
gt->give_node_coord2d  ok
*/
/**
 */
void newmeshread (char *filename,long lcid)
{
  long ip_dir,nod_spr;
  ip_dir = 0;
  nod_spr = 1;
  
  if ((Mp->adaptflag & 16))
    nod_spr = 0;
  
  
  // **********************************************************************
  // begin of step 1  -  aproximation `values` into nodes
  // **********************************************************************
  
  long ndofn,dim,nipcomp[3];
  double *nodvals_ip;
  
  dim = 2;
  ndofn = 2;
  
  nipcomp[0] = 0;   // number of transferd components of array ip[]->strain
  nipcomp[1] = 0;   // number of transferd components of array ip[]->stress
  nipcomp[2] = 5;   // number of transferd components of array ip[]->eqother        !!!  bez hardeningu
  
  if (!ip_dir)
    give_nodvals_ip (dim,nipcomp,nodvals_ip);
  else
    nodvals_ip = NULL;
  
  // **********************************************************************
  // begin of step 2  -  loading of new mesh
  // **********************************************************************
  
  problem *p_old = new problem;
  problem *p_new = new problem;
  
  p_old->globinic ();
  
  Mp = NULL; Gtm = NULL; Mt = NULL; Mm = NULL; Mc = NULL; Mb = NULL; Lsrs = NULL; Smat = NULL;
  fclose (Out);
  
  stochdriver *stochd;
  mefel_init (filename,stochd);
  
  p_new->globinic ();
  
  // **********************************************************************
  // begin of step 3  -  konverse
  // **********************************************************************
  
  long *parentel_ip,*parentel_nod;
  double **ipcoord_new;
  
  parentel_ip = new long [p_new->mm->tnip];
  parentel_nod = new long [p_new->mt->nn];
  allocm (p_new->mm->tnip,3,ipcoord_new);
  
  
  findout_parentel_ip (p_old->gt,Gtm,Mm,Mt,ipcoord_new,parentel_ip);
  findout_parentel_nod (p_old->gt,Gtm,parentel_nod,dim);
  
  //if (Mespr==1)   fprintf (stdout,"\n ***  Transformation of meshes  ***\n");
  
  if (!ip_dir)  transfvalues_ip_indirect (p_old->gt,Mm,                 ipcoord_new,parentel_ip,nipcomp,nodvals_ip );
  else          transfvalues_ip_direct   (p_old->gt,Mm,(const double **)ipcoord_new,parentel_ip,nipcomp,p_old->mm,dim);
  
  if (!nod_spr) transfvalues_nod (p_old,p_new,lcid,dim,ndofn,parentel_nod,'n');
  else          transfvalues_nod (p_old,p_new,lcid,dim,ndofn,parentel_nod,'y');
  
  
  delete [] parentel_ip;
  delete [] parentel_nod;
  destrm (ipcoord_new,Mm->tnip);
  
  // **********************************************************************
  // end of step 3  -  konverse
  // **********************************************************************
  
  p_old->dealoc ();
  p_new->deinic ();
  
  // **********************************************************************
  // end of step 2  -  loading of new mesh
  // **********************************************************************
  
  if (nodvals_ip!=NULL)  delete [] nodvals_ip;
  
  // **********************************************************************
  // end of step 1  -  aproximation `values` into nodes
  // **********************************************************************
  
  
  mefel_right_hand_side (lcid,Lsrs->rhs);
  
}

/**
   Function approximates values stored on integration points (in arrays `strain` `stress` `eqother`) to nodes - for "global" mesh.
   Used interpolation method is 'spr_smoothing'.
   Values in nodes are returned by array nodvalues.
   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])
   Necessary precondition is Mm->ip[i].ncompstr and Mm->ip[i].ncompother are constant all over the domain !!!!
   
   @param dim - dimension of solved problem
   @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 nodvalues - answer ; nonalocated
   
   created  9.12.2002, Ladislav Svoboda, [e-mail]
 */
long give_nodvals_ip (long dim,long nipcomp[],double *&nodvalues)
{
  long i,j,k,nip,ipp,nvals,*nsp;
  double **spcoord,**spvalue,**ipcoord;
  
  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];
  
  if (nvals==0){
    nodvalues = NULL;
    return 0;
  }
  else
    nodvalues = new double [nvals*Mt->nn];
  
  nsp = new long [Mt->ne];
  spcoord = new double* [Mt->ne];
  spvalue = new double* [Mt->ne];
  ipcoord = new double* [Mm->tnip];
  for (i=0;i<Mm->tnip;i++)
    ipcoord[i] = new double [3];
  
  allipcoord (ipcoord);
  
  for (i=0;i<Mt->ne;i++){
    nip = 0;
    for (j=0;j<Mt->elements[i].nb;j++)
      for (k=0;k<Mt->elements[i].nb;k++)
	nip += Mt->elements[i].nip[j][k];
    
    nsp[i] = nip;
    
    ipp = Mt->elements[i].ipp[0][0];
    spcoord[i] = new double [nip*dim];
    spvalue[i] = new double [nip*nvals];
    for (j=0;j<nip;j++){
      for (k=0;k<dim;k++)
	spcoord[i][j*dim+k] = ipcoord[ipp+j][k];
      
      for (k=0;k<nipcomp[0];k++)
	spvalue[i][j*nvals + k                          ] = Mm->ip[ipp+j].strain[k];
      
      for (k=0;k<nipcomp[1];k++)
	spvalue[i][j*nvals + k + nipcomp[0]             ] = Mm->ip[ipp+j].stress[k];
      
      for (k=0;k<nipcomp[2];k++)
	spvalue[i][j*nvals + k + nipcomp[0] + nipcomp[1]] = Mm->ip[ipp+j].eqother[k];
      
    }
  }
  
  least_square spr(dim,nvals,Gtm,0+((Mp->adaptflag & 16) ? 16 : 0));   // co to je 1 ??
  
  spr.spr_optional (Out,nsp,spcoord,spvalue,nodvalues,1);
  
  delete [] nsp;
  for (i=0;i<Mt->ne;i++){ delete [] spcoord[i]; delete [] spvalue[i]; }
  delete [] spcoord;
  delete [] spvalue;
  for (i=0;i<Mm->tnip;i++) delete [] ipcoord[i];
  delete [] ipcoord;
  
  return nvals;
}

/**
   Function returns array of deformations in ALL nodes of problem 'p'.
   Deformation position in answer array is (r[node_1][x] ... r[node_nn][x] , r[node_1][y] ... r[node_nn][y])
   
   @param lcid  - id of load case
   @param p     - pointr on problem
   @param rfull - answer - allocated 1D array, size is p->mt->nn * p->mm->ip[0].ndofn
   
   created  9.12.2002, Ladislav Svoboda, [e-mail]
 */
void give_rfull (long lcid,problem *p,double *rfull)
{
  long i,j,ndofn;
  double *r;

  ndofn = p->mt->give_ndofn (0);
  r = new double [ndofn];

  for (i=0;i<p->mt->nn;i++){
    noddispl (lcid,i,p,r);
    for (j=0;j<ndofn;j++)
      rfull[i+j*p->mt->nn] = r[j];
  }
  
  delete [] r;
}

/**
   function computes coordinates of all integration points in "global" mesh
   
   @param ipcoord - empty 1D array, size is Mt->tnip x 3
   
   created  5.12.2002, Ladislav Svoboda, [e-mail]
 */
void allipcoord (double **ipcoord)
{
  long i,ri,ci,nb,ipp;
  
  for (i=0;i<Mt->ne;i++){
    nb = Mt->elements[i].nb;
    for (ri=0;ri<nb;ri++){
      for (ci=0;ci<nb;ci++){
	ipp = Mt->elements[i].ipp[ri][ci];
	
	if (Mt->elements[i].nip[ri][ci]){
	  switch (Mt->give_elem_type(i)){
	  case planeelementlt:{ Pelt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; }
	  case planeelementqt:{ Peqt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; }
	  case planeelementlq:{ Pelq->ipcoordblock (i,ri,ci,ipcoord+ipp); break; }
	  case planeelementqq:{ Peqq->ipcoordblock (i,ri,ci,ipcoord+ipp); break; }
	  default:{
	    fprintf (stderr,"\n\n unknown element type is required in function  allipcoord (%s, line %d).\n",__FILE__,__LINE__);
	  }
	  }
	}
      }
    }
  }
}

/**
   Function alocates required arrays for function 'locintpoints',
   which finds out (for every integration point from new mesh) "parent" element (from old mesh),
   in which the integration point lays.
   Parent elements are returned in array "parentel_ip".
   
   @param gt_old - gtopology of old mesh
   @param gt_new - gtopology of new mesh
   @param mm_new - mechmat of new mesh
   @param mt_new - mechtop of new mesh
   @param ipcoord - array of coordinates of new mesh int. points, size is mm_new->tnip x dim
   @param parentel_ip - answer = 1D array, size is mt_new->tnip
   
   created  3.12.2002, Ladislav Svoboda, [e-mail]
 */
void findout_parentel_ip (gtopology *gt_old,gtopology *gt_new,mechmat *mm_new,mechtop *mt_new,double **ipcoord,long *parentel_ip)
{
  long *passedel1;  passedel1 = new long [gt_old->ne];
  long *passedel2;  passedel2 = new long [gt_new->ne];   memset (passedel2,0,gt_new->ne*sizeof(long));
  long *elheap;     elheap    = new long [gt_new->ne];   elheap[0] = 0;
  long *newelheap;  newelheap = new long [gt_new->ne];
  long *susel;      susel     = new long [gt_old->ne];
  long *newsusel;   newsusel  = new long [gt_old->ne];
  
  memset (parentel_ip,0,mm_new->tnip*sizeof(long));
  allipcoord (ipcoord);
  
  locintpoints (gt_old,gt_new,mt_new,ipcoord,passedel1,passedel2,1,elheap,newelheap,susel,newsusel,parentel_ip);
  
  delete [] passedel1;
  delete [] passedel2;
  delete [] elheap;
  delete [] newelheap;
  delete [] susel;
  delete [] newsusel;
}

/**
   For every integration point (from new mesh) function finds out "parent" element (from old mesh),
   in which the int. point lays.
   Parent elements are returned in array "parentel".

   @param gt_old - gtopology of old mesh
   @param gt_new - gtopology of new mesh
   @param mt_new - mechtop of new mesh