cpsolver.cpp

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

    */

    //if (nce!=0 || ncd!=0){
      //  assembling of displacement vector for new code numbers
      Mt->lhs_restore (r,n);
      //  assembling of vector of previous nodal forces for new code numbers
      Mt->rhs_restore (fp,n);
      //}

      Gtm->node_mark (n);

    Mm->updateipval();

    //  cleaning of the stiffness matrix
    //if (ncd>0){
    if (Smat != NULL){
      delete Smat;
      Smat=NULL;
    }    
    //}  


    /***************************************************/
    // TKO GRAPHICS
    /* if (nce || ncd){
       if (idcp > 0){
       print_close();
       }
       //compute_req_val (lcid);
       
       fprintf (stdout,"\n tisk do souboru");
       print_init(idcp, "wt");
       print_step(lcid, i, Mp->time, f);
       print_flush();
       idcp++;
       }
    // TKO GRAPHICS
    */
    /***************************************************/
    
    //  prescribed nodal forces
    mefel_right_hand_side (lcid,f);
    
    //  computation of sums of force components in particular directions
    //  it is used in special creep and consolidation models
    //Mb->comp_sum (f);

    //  forces caused by stress increments
    Mp->phase=1;
    internal_forces (lcid,fi);
    
    //fprintf (Out,"\n\n kontrola prave strany \n");
    
    //  right hand side assembling and storage of actual prescribed forces
    for (j=0;j<n;j++){
      //fprintf (Out,"\n f %4ld  %10.5le       fp %4ld  %10.5le    rozdil %10.5le",j,f[j],j,fp[j],f[j]-fp[j]);

      fb[j]=f[j]-fp[j]+fi[j];
      //fb[j]=f[j]-fp[j];
      fp[j]=f[j];
    }
    
    //  allocation and assembling of the stiffness matrix
    stiffness_matrix (lcid);

    //  solution of K(r).v=F
    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);
    //Smat->solve_system (Gtm,dr,fb);
    
    //fprintf (Out,"\n\n kontrola reseni \n");
    
    //  new displacement vector
    for (j=0;j<n;j++){
      r[j]+=dr[j];
      //fprintf (Out,"\n dr %4ld  %10.5le      r %4ld  %10.5le",j,dr[j],j,r[j]);
    }

    //  storage of displacements to nodes (possible changes in code numbers)
    //old:???????
    Mt->lhs_save (r,1);
    //Mt->lhs_save (dr,1);
    //  storage of prescribed forces to nodes (possible changes in code numbers)
    Mt->rhs_save (f);
    
    
    //  computation of internal forces
    Mp->phase=2;
    internal_forces (lcid,fi);


    /***************************************************/
    
    //old
    // computation of values required for output    
    compute_req_val (lcid);
    
    //  time increment
    prev_time=Mp->time;
    Mp->time = Mp->timecon.newtime ();
    i++;
    
    if (Mp->timecon.isitimptime ()==1){
      mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
    }
    
    
    //fprintf (stdout,"\n tisk do souboru");
    //print_init(i, "wt");
    print_step(lcid, i, Mp->time, f);
    print_flush();
    
    /***************************************************/
    /***************************************************/
    
    /*
    //new:
    norfa=ss(f,f,n);
    for (j=0;j<n;j++){
      fb[j] = f[j] - fi[j];
    }
    //  norm of vector of unbalanced forces
    norf=ss(fb,fb,n);
    if (norfa > 0.0)
      norf /= norfa;
    
    if (norf<err){
      //  backup of nodal values because code numbers may be changed
      //new:
      Mt->lhs_save (r,1);
      Mt->rhs_save (f);
      
      // computation of values required for output    
      compute_req_val (lcid);
      //  time increment
      prev_time=Mp->time;
      Mp->time = Mp->timecon.newtime ();
      Mp->timecon.take_values (&Mp->timecon);
      
      if (Mp->timecon.isitimptime ()==1){
	mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
      }
      
      print_step(lcid, i, Mp->time, f);
      print_flush();
      i++;	
    }
    
    // iteration of unbalanced forces caused by time independent models
    //  internal iteration loop
    else{
      for (j=0;j<ini;j++)
	{
	  //  solution of K(r).v=F
	  Smat->solve_system (dr,fb);
	  
	  for (k=0;k<n;k++)
	    r[k]+=dr[k];
	  
	  //  computation of internal forces
	  internal_forces (lcid,fi);
	  
	  //  vector of unbalanced forces
	  for (k=0;k<n;k++)
	    fb[k]=f[k]-fi[k];
	  
	  //  norm of vector of unbalanced forces
	  norf=ss(fb,fb,n);
	  if (norfa > 0.0)
	    norf /= norfa;
	  
	  if (Mespr==1) fprintf (stdout,"\n Time increment %lf   inner loop j %ld     norf=%e",Mp->time,j,norf);
	  if (norf<err)  break;
	}
      if (j==ini || norf>err)
	{
	  if (Mespr==1) fprintf (stdout,"\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
	  
	  //  backup of nodal values because code numbers may be changed
	  //new:
	  Mt->lhs_save (r,1);
	  Mt->rhs_save (f);
	  
	  // equilibrium state was reached not - added by TKr
	  compute_req_val (lcid);
	  //  time increment
	  prev_time = Mp->time;
	  Mp->time = Mp->timecon.newtime ();
	  Mp->timecon.take_values (&Mp->timecon);
	  
	  if (Mp->timecon.isitimptime ()==1)
	    mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
	  print_step(lcid, i, Mp->time, f);
	  print_flush();
	  i++;
	}
      else
	{
	  //  backup of nodal values because code numbers may be changed
	  //new:
	  Mt->lhs_save (r,1);
	  Mt->rhs_save (f);
	  
	  // equilibrium state was reached
	  compute_req_val (lcid);
	  //  time increment
	  prev_time = Mp->time;
	  Mp->time = Mp->timecon.newtime ();
	  Mp->timecon.take_values (&Mp->timecon);
	  
	  if (Mp->timecon.isitimptime ()==1)
	    mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
	  print_step(lcid, i, Mp->time, f);
	  print_flush();
	  i++;
	}
    }
    */

  }while(Mp->time<end_time);
  
  print_close ();
  
  //fclose (kon);
  
  
  if (Mp->timecon.nit > 0)
    fclose (aux);
  
  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f, delete [] fp;
}

/*




    if (nce!=0 || ncd!=0)
    
    //  toto umisteni je jen pokusne
    Mm->updateipval();
    
    fprintf (stdout,"\n\n cas %le",Mp->time);
    fprintf (stdout,"\n pocet zmen na prvcich   %ld",nce);

    Gtm->print_time_functions (kon,Mp->time);



    if (nce>0){
      //  nce>0 means changed elements
      
      //  attained displacements have to be saved at new elements
      //Mt->initial_displ ();
      
      
      fprintf (stdout,"\n\n pocitame korekci posunu na zapnute casti v case %le",Mp->time);
      
      
      //  searching for new nodes
      ncd = Gtm->search_newdofs (Mp->time,prev_time);
      
      
      
      
      
      
      //  definition of DOFs
      Ndofm = Gtm->gencodnum ();
      n=Ndofm;
      if (Smat != NULL){
	delete Smat;
	Smat=NULL;
      }
      

      rhs_graphic (fb,lcid=0,n);

      //  assembling of tangent stiffness matrix
      stiffness_matrix (lcid);
      
      //  solution of K(r).v=F
      Smat->solve_system (dr,fb);
      
      
      //  backup of nodal values because code numbers may be changed
      Mt->lhs_save (dr,2);
      
      
      //  array leso is restored
      //  it was changed in function Gtm->search_newelem (time)
      //Gtm->restore_leso ();
      //Gtm->cn_restore ();
      
      
      Gtm->update_dofs (Mp->time);
    }








    ncd = Gtm->search_newdofs (Mp->time,prev_time);
    Gtm->update_dofs (Mp->time);
    Ndofm = Gtm->gencodnum ();
    n=Ndofm;
    
    //  update of nodes, DOFs and elements
    //Gtm->update_mesh (Mp->time,prev_time,ncn,nce);
    
    //  set up of correct values in array leso
    Gtm->update_elemstat (Mp->time);
    
    
    Gtm->print_time_functions (kon,Mp->time);
    
    fprintf (stdout,"\n cas %le  zmeny v DOF %ld    zmeny v prvcich %ld",Mp->time,ncd,nce);
    
    
    
    //  definition of DOFs
    if (ncd>0){
      //Ndofm = Gtm->gencodnum ();
      //n=Ndofm;

      if (Smat != NULL){
	delete Smat;
	Smat=NULL;
      }

      //  assembling of lhs array for new code numbers
      Mt->lhs_restore (r,n);
      Mt->rhs_restore (fp,n);
    }
    


    long k;
    fprintf (kon,"\n\n\n\n novy zapis (cas %le)\n\n\n\n",Mp->time);
    for (j=0;j<Mt->nn;j++){
      if (Gtm->lnso[j]==1){
	fprintf (kon,"%5ld     ",j+1);
	for (k=0;k<Mt->give_ndofn (j);k++){
	  fprintf (kon,"%5ld   ",Gtm->gnodes[j].cn[k]);
	}
	fprintf (kon,"  %le   %le   %le\n",Gtm->gnodes[j].x,Gtm->gnodes[j].y,Gtm->gnodes[j].z);
      }
    }

    //fclose (kon);
    //abort ();

    
    fprintf (stdout,"\n pocet neznamych  %ld",n);
    
    if (nce || ncd){
      if (idcp > 0){
	print_close();
      }
      compute_req_val (lcid);
      
      fprintf (stdout,"\n tisk do souboru");
      print_init(idcp, "wt");
      print_step(lcid, i, Mp->time, f);
      print_flush();
      idcp++;
    }
    
    //fprintf (stdout,"\n time %e,      NDOF %ld\n",Mp->time,Ndofm);
    
    //for (j=0;j<29;j++){
    //fprintf (stdout,"%ld %ld\n",j,Gtm->gf[j].getval_long (Mp->time));
    //}
    


    //  assignment of initial displacements
    //  it works for problem with growing number of elements
    //  it does not work generally for problems with increasing and decreasing number of elements
    if (nce>0){
      Mt->initial_displ ();
    }


    fprintf (kon,"\n\n cas  %lf \n\n",Mp->time);
    fprintf (kon,"\n\n kontrola kodovych cisel \n\n");
    long k,ndofn;
    for (j=0;j<Mt->nn;j++){
      fprintf (kon,"\n uzel %ld",j+1);
      ndofn=Gtm->gnodes[j].ndofn;
      for (k=0;k<ndofn;k++){
	fprintf (kon,"  %ld",Gtm->gnodes[j].cn[k]);
      }
    }
    fprintf (kon,"\n\n");
    
    fprintf (kon,"\n\n\n kontrola pocatecnich posunu na prvcich \n");
    for (j=0;j<Gtm->ne;j++){
      fprintf (kon,"\n element %ld\n",j+1);
      long ndofe=Gtm->give_ndofe (j);
      for (k=0;k<ndofe;k++){
	fprintf (kon,"  %lf",Mt->elements[j].initdispl[k]);
      }
    }
    fprintf (kon,"\n\n");
    

    
    Mm->updateipval();
    //fprintf (stdout,"\n time %e,      NDOF %ld",Mp->time,Ndofm);
    //    fprintf (Out,"\n\n\n\n time %e",Mp->time);
    
    //  prescribed nodal forces
    mefel_right_hand_side (lcid,f);

    // initialization of material models, initial temperatures are
    // computed during right hand side assembling -> right hand side
    // have to be reassembled due to correct temprstrains
    if (i==0)
    {
      initmaterialmodels();
      mefel_right_hand_side (lcid,f);
    }
    
    //Mb->comp_sum (f);
    

    //  unbalanced forces
    Mp->phase=1;
    //Mm->est=eigstrain;
    internal_forces (lcid,fi);
    
    fprintf (kon,"\n\n\n kontrola velkych vektoru \n");
    for (j=0;j<n;j++){
      fprintf (kon,"\n %4ld  f %12.8lf   fp %12.8lf   fi %12.8lf   fb %12.8lf",j,f[j],fp[j],fi[j],f[j]-fp[j]);
    }
    fprintf (kon,"\n\n");

    //  right hand side assembling
    for (j=0;j<n;j++){
      //fb[j]=f[j]-fp[j]+fi[j];
      fb[j]=f[j]-fp[j];
      fp[j]=f[j];
    }
    
    //fprintf (Out,"\n\n kontrola prave strany \n");
    //for (i=0;i<n;i++){
    //fprintf (Out,"\n fb %5ld   %le",i,fb[i]);
    //}
    
    
    //  assembling of tangent stiffness matrix
    stiffness_matrix (lcid);
    
    //  solution of K(r).v=F
    Smat->solve_system (dr,fb);
    
    //  new displacement vector
    for (j=0;j<n;j++){
      r[j]+=dr[j];
    }
    //  backup of nodal values because code numbers may be changed
    Mt->lhs_save (r,1);
    Mt->rhs_save (f);
    
    fprintf (kon,"\n\n kontrola reseni \n\n");
    for (j=0;j<Mt->nn;j++){
      fprintf (kon,"\n uzel %ld",j+1);
      ndofn=Mt->give_ndofn (j);
      long cn;
      for (k=0;k<ndofn;k++){
	cn=Mt->give_dof (j,k);
	if (cn==0)
	  fprintf (kon,"  %lf",0.0);
	else
	  fprintf (kon,"  %lf",dr[cn-1]);
      }
      fprintf (kon,"   ");
      for (k=0;k<ndofn;k++){
	cn=Mt->give_dof (j,k);
	if (cn==0)
	  fprintf (kon,"  %lf",0.0);
	else
	  fprintf (kon,"  %lf",r[cn-1]);
      }

    }
    
    
    
    
    

    long kj,ndofn;
    double *jk;
    fprintf (Out,"\n\n\n cas %le\n",Mp->time);
    for (j=0;j<Mt->nn;j++){
      ndofn = Mt->give_ndofn (j);
      jk = new double [ndofn];
      noddispl (0,jk,j);
      fprintf (Out,"\n\n uzel  %5ld",j+1);
      for (kj=0;kj<ndofn;kj++){
	fprintf (Out,"\n slozka %5ld   d %le   id %le   d+id %le",kj,jk[kj],Mt->nodes[j].nvgraph[kj],jk[kj]+Mt->nodes[j].nvgraph[kj]);
      }
      delete jk;
    }

    
    
    Mm->tempstr_eigstr();
    //  computation of internal forces
    Mp->phase=2;
    internal_forces (lcid,fi);

    // computation of values required for output    
    compute_req_val (lcid);

    //  time increment
    prev_time=Mp->time;
    Mp->time = Mp->timecon.newtime ();
    i++;
    
    
    if (Mp->timecon.isitimptime ()==1){
      mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
    }

    
    fprintf (stdout,"\n tisk do souboru");
    print_init(i, "wt");
    print_step(lcid, i, Mp->time, f);
    print_flush();
    //print_close();
    
  }while(Mp->time<end_time);
  
  print_close ();
  
  fclose (kon);


  if (Mp->timecon.nit > 0)
    fclose (aux);

  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f, delete [] fp;
}

*/