lssolver.cpp

来自「Finite element program for mechanical pr」· C++ 代码 · 共 222 行

#include "lssolver.h"
#include "global.h"
#include "globmat.h"
#include "loadcase.h"
#include "gmatrix.h"
#include "mechprint.h"
#include "node.h"
#include "elemswitch.h"
#include <string.h>

/**
   function solves linear statics problems
   
   function contains the following steps:
   1. stiffness matrix assembling
   2. right hand side assembling (load vectors assembling)
   3. solution of equation systems by choosen method
   4. strain computation (eligible)
   5. stress computation (eligible)
   6. reaction computation (eligible)
   7. results printing
   
   JK
*/
void solve_linear_statics ()
{
  long i;
  double *lhs,*rhs;
  
  //  stiffness matrix assembling
  stiffness_matrix (0);
  
  //  array is allocated in Lsrs
  lhs = Lsrs->give_lhs (0);
  rhs = Lsrs->give_rhs (0);
  
  //  loop over load cases
  for (i=0;i<Lsrs->nlc;i++){
    //  right hand side of system
    //Mb->lc[i].assemble (i,rhs+i*Ndof);
    mefel_right_hand_side (i,rhs);
  }
  
  Mm->tempstr_eigstr();
  
  
  
  /*
  Smat->printmat(Out);
  FILE *p;
  p = fopen ("vektor.txt","w");
  fprintf (p,"%ld\n",Ndofm);
  for (i=0;i<Ndofm;i++){
    fprintf (p,"%20.16le\n",rhs[i]);
  }
  fclose (p);
  */
  
  
  
  
  //Smat->printdiag(Out);
  
  //  solution of equation system
  for (i=0;i<Lsrs->nlc;i++){
    Mp->ssle.solve_system (Gtm,Smat,Lsrs->give_lhs(i),Lsrs->give_rhs(i),Out);
    //Smat->solve_system (Lsrs->give_lhs(i),Lsrs->give_rhs(i));
  }
  
  //Smat->printdiag(Out);
  
  
  /*
  p = fopen ("reseni.txt","w");
  fprintf (p,"%ld\n",Ndofm);
  for (i=0;i<Ndofm;i++){
    fprintf (p,"%20.16le\n",lhs[i]);
  }
  fclose (p);
  */

  //  zacatek Risoviny
  /*
    long j;
  long n=6;
  double *condmat,*condvect,*condsolv;
  condsolv = new double [n];
  condvect = new double [n];
  condmat = new double [n*n];
  nullvr (condmat,n*n);
  nullvr (condvect,n);
  nullvr (condsolv,n);

  //Smat->condense (Gtm,condmat,lhs,rhs,3,1);
  Smat->condense (Gtm,condmat,lhs,rhs,n,1);

  fprintf (Out,"\n\n kondenzovany vektor");
  j=0;
  for (i=Ndofm-n;i<Ndofm;i++){
    condvect[j]=rhs[i];
    fprintf (Out,"\n %e",condvect[j]);
    j++;
  }
  
  fprintf (Out,"\n\n kondenzovana matice");
  for (i=0;i<n;i++){
    fprintf (Out,"\n");
    for (j=0;j<n;j++){
      fprintf (Out," %e",condmat[i*n+j]);
    }
  }

  gemp (condmat,condsolv,condvect,n,1,1.0e-15,1);
  
  fprintf (Out,"\n\n vyresena cast vysledku");
  j=0;
  for (i=Ndofm-n;i<Ndofm;i++){
    lhs[i]=condsolv[j];
    fprintf (Out,"\n %e",condsolv[j]);
    j++;
  }
  
  //  back substitution on subdomains
  //Smat->condense (Gtm,condmat,lhs,rhs,3,2);
  Smat->condense (Gtm,condmat,lhs,rhs,n,2);
  
  delete [] condmat;
  delete [] condvect;
  delete [] condsolv;
  */
  //  konec Risoviny
  


  /*
  //  zacatek JK testu
  long j,nse=6,dim;
  long *se;
  double *rbm;
  
  se = new long [nse];
  rbm = new double [nse*Ndofm];
  
  Smat->kernel (rbm,dim,se,nse,1.0e-3,3);
  
  fprintf (Out,"\n\n\n RIGID BODY MOTIONS \n");
  fprintf (Out,"\n pocet RBM  %ld\n",dim);
  for (i=0;i<dim;i++){
    fprintf (Out," %ld",se[i]);
  }
  for (i=0;i<dim;i++){
    fprintf (Out,"\n\n RBM number %ld",i);
    for (j=0;j<Ndofm;j++){
      fprintf (Out,"\n %4ld %4ld  %e",i,j,rbm[i*Ndofm+j]);
    }
  }
  //  konec JK testu
  */


/*
  //  loop over load cases
  for (i=0;i<Lsrs->nlc;i++){
    
    //  strains
    if (Mp->straincomp==1){
      Mm->computestrains (i);
      Mm->stra.transformvalues (1);
    }
    
    //  stresses
    if (Mp->stresscomp==1){
      Mm->computestresses (i);
      Mm->stre.transformvalues(0);
    }
    
    //  reactions
    if (Mp->reactcomp==1){
      Mb->lc[i].compute_reactions (i);
    }
    

  }
*/
  
  //  reactions
  /*
  if (Mp->reactcomp==1){
    Mb->lc[0].compute_reactions (0);
  }
  */
  
  //  pondeli 30.1.2006
  //local_global_displ_transf (0);
  
  //for (i=0;i<Mt->nn;i++){
  //fprintf (Out,"\n uzel %ld   %ld",i,Mt->nodes[i].transf);
  //}
  
  print_init(-1, "wt");    
  for (i=0;i<Lsrs->nlc;i++){
    //  computes and prints required quantities
    //compute_ipstrains (i);
    //compute_ipstresses (i);
    compute_req_val (i);
    print_step(i, 0, 0.0, NULL);    
  }
  print_close();

  
  if (Mp->adaptivity)
    Ada->run (2,0,0); // 2 - for testing
  
  /*
  for (i=0;i<Ndofm;i++){
    printf ("\n posun %20.10lf",Lsrs->lhs[i]);
  }
  */
  
}