mefelinit.cpp

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "mefelinit.h"

#include "global.h"
#include "globmat.h"
#include "sequent.h"
#include "intpoints.h"
#include "gtopology.h"
#include "meshtransfer.h"
#include "adaptivity.h"
#include "mechprint.h"
#include "nssolver.h"
#include "element.h"
#include "elemswitch.h"
#include "iotools.h"
#include "xfile.h"

/**
   function reads data from input file and initiates all
   necessary variables and arrays
   
   @param argv -
   @param stochd - pointer to %stochdriver
   
*/
void mefel_init (char *argv,stochdriver *stochd)
{
  XFILE *in;
  long merr;

  fprintf (stdout,"\n\n ******************************************************\n");
  fprintf (stderr,"\n\n ******************************************************\n");
  
  Mp  = new probdesc;
  Gtm = new gtopology;
  Mt  = new mechtop;
  Mm  = new mechmat;
  Mc  = new mechcrsec;
  Mb  = new mechbclc;
  Outdm = new outdriverm;
  
  in = xfopen(argv,"r");
  if (in==NULL){
    fprintf (stderr,"\n Input file has not been specified.");
    fprintf (stderr,"\n Try it again!\n\n");
    abort ();
  }
  in->maxlnsize = 1000;
  in->warrning = 1;
  in->kwdmode = ignore;
  //in->kwdmode = sequent_mode;
  in->ignorecase = 1;
  
  // **************************
  // **  basic problem data  **
  // **************************
  Mp->read (in);
  
  // **************************************
  // **  input of node and element data  **
  // **************************************
  Mt->read (in);
  
  /*
  Out=fopen("dupnode","w");

  Gtm->edges (Out);
  Gtm->edgenode_sorting ();
  Gtm->prev_next_edges ();
  Gtm->edge_series ();
  Gtm->edge_dirvect ();
  Gtm->edge_normvect ();
  
  long i,j,a,n;
  
  
  //fprintf (Out,"\n\n kontrola prev a next \n");
  //for (i=0;i<Gtm->nged;i++){
    //fprintf (Out,"\n hrana %6ld   prev %6ld  next %6ld",i,Gtm->gedges[i].prev,Gtm->gedges[i].next);
  //}
  //fprintf (Out,"\n\n");

  
  fprintf (Out,"\n\n kontrola souslednosti\n");
  long *av;
  av = new long [Gtm->nged];
  for (i=0;i<Gtm->nged;i++){
    av[i]=-1;
  }
  for (i=0;i<Gtm->nged;i++){
    if (Gtm->gedges[i].prev==-1){
      a=i;
      n=Gtm->gedges[i].next;
      fprintf (Out,"\n %6ld -> %6ld   %lf %lf %lf      %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],
	       Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]);
      av[i]=1;
      break;
    }
  }
  long ui=Gtm->nged;
  
  for (i=1;i<ui;i++){
    a=n;
    n=Gtm->gedges[n].next;
    av[a]=1;
    fprintf (Out,"\n %6ld -> %6ld   %lf %lf %lf      %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],
	     Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]);
    if (n<0){
      for (j=0;j<Gtm->nged;j++){
	if (Gtm->gedges[j].prev==-1 && av[j]==-1){
	  a=j;
	  n=Gtm->gedges[j].next;
	  av[j]=1;
	  ui--;
	  fprintf (Out,"\n %6ld -> %6ld   %lf %lf %lf      %lf %lf %lf",a,n,Gtm->gedges[a].dv[0],Gtm->gedges[a].dv[1],Gtm->gedges[a].dv[2],
		   Gtm->gedges[a].nv[0],Gtm->gedges[a].nv[1],Gtm->gedges[a].nv[2]);
	  break;
	}
      }
    }
  }



  //a=0; n=Gtm->gedges[a].next;
  //fprintf (Out,"\n\n kontrola serie \n");
  //for (i=0;i<Gtm->nged;i++){
    //fprintf (Out,"\n%ld -> %ld",a,n);
    //a=n;
    //n=Gtm->gedges[a].next;
  //}


  for (i=0;i<Gtm->nged;i++){
    fprintf (Out,"\n\n hrana cislo  %ld",i);
    Gtm->gedges[i].print (Out);
  }

  
  fprintf (Out,"\n\n\n kontrola serii \n");
  for (i=0;i<Gtm->nged;i++){
    fprintf (Out,"\n edge %6ld je v serii  %6ld",i,Gtm->edgeser[i]);
  }

  fclose (Out);
  */


  // mesh control TKr
  if (Mp->tprob == growing_mech_structure){
    if (Mespr==1){
      merr=Mt->mesh_check();
      if (merr > 0){
	fprintf (stdout,"Mesh error = %ld\n\n\n",merr);
	exit(0);
      }
    }
  }

  // zacatek TKo
  if (Mp->tprob == earth_pressure)
    {
      Gtm->bckcn = new long*[Gtm->nn];
      memset(Gtm->bckcn, 0, sizeof(*Gtm->bckcn)*Gtm->nn);
      Gtm->backup_codnum();
    }
  //  konec TKo

  // *************************************
  //  generation of ordering of unknowns
  //  (generation of code numbers)
  //  number of degrees of freedom
  // *************************************
  //Ndofm = Gtm->gencodnum ();
  //Out=fopen("gtmzapis","w");
  Ndofm = Gtm->codenum_generation (Out);
  //fclose (Out);
  
  if (Mespr==1)
    fprintf (stdout,"\n number of degrees of freedom    %ld",Ndofm);
  //  zacatek TKr
  //Mt->store_code_num_elem(); // enable combination of plane elements, springs and bar elems with beams//tady??!!!
  //  konec TKr
  
  //  zacatek JK
  if (Mp->tprob==layered_linear_statics){
    //  definition of type of layered problem
    if (Mt->elements[0].te<20)  Mp->tlm=1;
    if (Mt->elements[0].te>20)  Mp->tlm=2;
    
    //  generation of ordering of Lagrange multipliers
    Mt->gencodnumlagrmult (Ndofm);
    Gtm->initiate_elemmult ();
    if (Mespr==1)
      fprintf (stdout,"\n total number of degrees of freedom    %ld",Ndofm);
  }


  if (Mp->tprob==lin_floating_subdomain){
    //  auxiliary informations about ordering of unknowns
    Gtm->flsub.ndof_subdom (Gtm->gnodes);
  }
  //  konec JK  
  
  // ***********************************************
  // **  input of material models and parameters  **
  // ***********************************************

  //  reading of material models and parameters
  //  allocation of integation points
  //  definition of the number of material models (nm) on integration points
  //  definition of material types (tm) on integration points
  //  definition of number of material model (idm) on integration points
  Mm->read (in);

  //  definition of variables ncompstr on intergation points
  //  definition of variables ssst on intergation points
  Mm->init_ip_1 ();
  
  // **********************************************
  // **  input of cross section characteristics  **
  // **********************************************
  Mc->read (in);
  
  // ************************
  // **  input of loading  **
  // ************************
  Mb->read (in);
  
  // ********************************************************
  // **  input of additional data necessary for problems   **
  // **  with changing number of nodes, elements and DOFs  **
  // ********************************************************
  if (Mp->tprob == growing_mech_structure){
    Gtm->read_gf (in);
    Mt->alloc_growstr ();
  }
  
  // ******************************************
  // **  definition of output from the code  **
  // ******************************************
  Outdm->read(in);
  
  // ************************************************
  // **  assembling of output files from the code  **
  // ************************************************
  char fname[1020];
  char *path;
  char *name;
  char *suffix;
  Mp->filename_decomposition (Outdm->outfn,path,name,suffix);
  Mp->path = path;
  Mp->filename=name;
  Mp->suffix=suffix;
  sprintf(fname, "%s%s.log", path, name);

  delete [] path;
  delete [] name;
  delete [] suffix;
  Mp->path = NULL;
  Mp->filename=NULL;
  Mp->suffix=NULL;

  Out = fopen(fname,"wt");
  
  
  //  zacatek JK
  //  toto jsou uzivatelsky definovane body
  //Mm->stra.init(strain);  
  //Mm->stre.init(stress);  
  //  konec JK
  
  //  zacatek TKo
  //  definition of nonlocal models and models for temperature
  Mm->init_ip_2 ();
  //  konec TKo
  
  
  //  allocation of arrays defined on nodes
  Mt->alloc_nodes_arrays ();
  
  //  number of components in arrays used in outdriver
  //  maximum number of components in array strain/stress on nodes and elements
  Mt->give_maxncompstr(Mm->max_ncompstrn, Mm->max_ncompstre);
  //  maximum number of components in array other on nodes and elements
  Mt->give_maxncompo(Mm->max_nncompo, Mm->max_encompo);
  //  konec TKo
  
  //  zacatek JK
  //if (Mp->straincomp==1){
  //Mm->stra.alloc(Mb->nlc);
  //}
  //if (Mp->stresscomp==1){
  //Mm->stre.alloc(Mb->nlc);
  //}
  //  konec JK
  
  
  //  prescribed displacements on elements investigation
  Mt->elemprescdisp ();
  //  prescribed temperatures on elements investigation
  if (Mp->tprob == linear_statics){
    Mt->elempresctemp (0);
  }
  //  reactions on elements investigation
  if (Mp->reactcomp==1){
    Mt->comreac ();
  }

  // *************************************************************
  //  vectors on left and right hand side of system of equations
  // *************************************************************
  Lsrs = new lhsrhs;
  //  allocation of arrays lhs, rhs (tdlhs, stdlhs if necessary)
  Lsrs->alloc ();
  //  function reads initial conditons
  Lsrs->initcond (in);


  //  computation of initial values at integration points from initial nodal values
  if (Mb->ico > 0)
    Mb->inicipval();
  
  
  
  if (Mp->matmodel == nonlocal){    
    //  list of adjacent elements
    Gtm->adjacelem (Out);
    //  list of adjacent integration points
    adjacip ();
    
    ipvolume ();
  }
  
  
  // adaptivity
  if (Mp->adaptivity)
    Ada = new adaptivity ();
  
  //  nondeterministic computation
  if (Mp->stochasticcalc){
    stochd->read (in);
    St = stochd;
  };
  
  //  initiation and computation of eigenstrains
  if (Mp->temperature==1){
    Mm->alloceigstrains ();
    Mm->alloctempstrains ();
  }
  
  
  Mt->nodedisplacements ();
  
  
  xfclose (in);

}