strut_xalm.cpp

Zobrazit dokumentaci tohoto souboru.

#include "strut_xalm.h"
#include <math.h>
#include <stdio.h>


namespace xalm {


//  funkce pro vypocet deformaci a vnitrnich sil na prvku
//  In:
//    x - vektor souradnic uzlu prvku  (x1,y1,x2,y2)
//    E-  modul pruznosti
//    r - vektor uzlovych posunu prvku (u1,v1,u2,v2)
// 
//  Out:
//    eps deformace
//    s   napeti
// 
void truss2d_postpro (const double *x, double E, const double *r, double &eps, double &s)
{
  double l2 = ( (x[2]-x[0])*(x[2]-x[0]) + (x[3]-x[1])*(x[3]-x[1]) );
  
  double x12 = x[0] - x[2];
  double y12 = x[1] - x[3];
  double x21 = x[2] - x[0];
  double y21 = x[3] - x[1];
  
  double u12 = r[0] - r[2];
  double v12 = r[1] - r[3];
  double u21 = r[2] - r[0];
  double v21 = r[3] - r[1];

  eps = (1/l2)*(x12*r[0]+y12*r[1]+x21*r[2]+y21*r[3]) + (1/(l2*2))*(u12*r[0]+v12*r[1]+u21*r[2]+v21*r[3]);
  s   = E * eps;
}


// funkce pro vypocet vektoru vnitrnich uzlovych sil fint na prvku
// In:
//   x - vektor souradnic uzlu prvku  (x1,y1,x2,y2)
//   r - vektor uzlovych posunu prvku (u1,v1,u2,v2)
//   s - vnitrni sily = napeti
//   A - plocha
//
// Out:
//   fint
//
void truss2d_fint (const double *x, const double *r, double s, double A, double *fint)
{
  double l = sqrt( (x[2]-x[0])*(x[2]-x[0]) + (x[3]-x[1])*(x[3]-x[1]) );
  
  double xu[4];
  for (int i = 0; i<4; i++)
    xu[i] = x[i] + r[i];
  
  fint[0] = s * (A/l) * (xu[0]-xu[2]);
  fint[1] = s * (A/l) * (xu[1]-xu[3]);
  fint[2] = s * (A/l) * (xu[2]-xu[0]);
  fint[3] = s * (A/l) * (xu[3]-xu[1]);
}

void truss2d (const double *x, double EA, double ke[4][4])
{
  double length = sqrt( (x[2]-x[0])*(x[2]-x[0]) + (x[3]-x[1])*(x[3]-x[1]) );
  
  double c = (x[2]-x[0]) / length;  double c2 = c*c;
  double s = (x[3]-x[1]) / length;  double s2 = s*s;
  
  double k = (EA/length);
  
  ke[0][0] = k * c2;    ke[0][1] = k * c*s;   ke[0][2] = k * -c2;   ke[0][3] = k * -c*s;
  ke[1][0] = k * c*s;   ke[1][1] = k * s2;    ke[1][2] = k * -c*s;  ke[1][3] = k * -s2;
  ke[2][0] = k * -c2;   ke[2][1] = k * -c*s;  ke[2][2] = k * c2;    ke[2][3] = k * c*s;
  ke[3][0] = k * -c*s;  ke[3][1] = k * -s2;   ke[3][2] = k * c*s;   ke[3][3] = k * s2;
}

double norm (double val)
{
  if (val < 0) val *= -1;
  return val;
}


void XALM_interface :: initialize (void)
{
  // ***   nastaveni promennych rodicovske tridy XALM   ***
  nsteps = 34;
  
  minStepLength = 0.01;
  maxStepLength = 0.2;
  initialStepLength = 0.2;
  
  Psi    = 0.0;       // displacement control on
  
  rtolf = rtold = 1.e-3;
  
  neq = 1;
  
  
  // lokalni pozadavky, zavisle na FEM baliku ktery inkluduje XALM knihovnu
  // sestavi matice tuhosti, vektor zatizeni atd.
  this->initialize_local();
}

void XALM_interface :: update_step (void)
{
  int step = give_step();
  
  tisk1[step+1] =  this->give_totalDisplacement()[0];
  tisk2[step+1] = -this->give_loadLevel();
  
  if (step == this->nsteps-1){
    FILE *file = fopen("graf1.xalm.dat", "w");
    for (long i = 0; i <=nsteps; i++)
      fprintf (file, "%g %g\n", tisk1[i], tisk2[i]);
    
    fclose(file);
  }
}


void XALM_interface :: initialize_local (void)
{
  // sestavit incrementalLoadVector
  incrementalLoadVector.resize_ignore_vals(neq);
  incrementalLoadVector[0] = -1;
  
  //
  initialLoadVector.resize_ignore_vals(neq);
  initialLoadVector[0] = 0;
  
  
  // ***   zde se nastavi popis vzperadla   ***
  
  // souradnice uzlu  1. uzel  2. uzel
  double xz[2][2] = { {0,0},   {4,3} };
  
  // vektor souradnic uzlu prvku (x1, y1, x2, y2)
  xze[0] = xz[0][0];
  xze[1] = xz[0][1];
  xze[2] = xz[1][0];
  xze[3] = xz[1][1];
  
  // modul pruznosti
  E = 1.0;
  
  // plocha
  A = 1;
  
  // typ matice tuhosti (false - elasticka, 1-tecna)
  tangent_stiff_mtrx = false;
  
  // matice tuhosti, pro vzperadlo ma rozmer 1x1
  k = 0;
  
  tisk1 = new double[nsteps+1];  memset (tisk1, 0, (nsteps+1)*sizeof(double));
  tisk2 = new double[nsteps+1];  memset (tisk2, 0, (nsteps+1)*sizeof(double));
  
  tisk1[0] = 0.0;
  tisk2[0] = 0.0;
}

void XALM_interface :: update_stiffness_matrix (const Dvctr *X)
{
  // vypocet tecne matice tuhosti, pokazde kdyz je vznesen pozadavek
  if (tangent_stiff_mtrx) {
    errol;
    //sestaveni linearizovane podminky rovnovahy
    //         //       ks = truss2d_initialstress(xze, A, s(1));
    //         //       kg = truss2d_initialdeformation (xze, r(lm(1,:)), E*A);
    //         //       if mtype == 1
    //         //         k = klin +kg+ks;
    //         //       else
    //         //         k = klin;
    //         //       end
    // tady dat asi vypocet    k = klin +kg+ks;
  }
  // vypocet elasticke matice tuhosti, jen v prvnim kroku
  else {
    if (this->give_step() == 0) {
      double klin[4][4];
      truss2d (xze, E*A, klin);
      k = klin[3][3];
    }
  }
}

void XALM_interface :: update_internal_forces (Dvctr *internalForces, const Dvctr *X)
{
  double eps, s;
  double fint[4];
  
  // vektor vsech posunu na prutu
  double r[4];
  r[0] = r[1] = r[2] = 0.0;
  r[3] = (*X)[0];
  
  // dopocteni deformaci a napeti
  truss2d_postpro (xze, E, r, eps, s);
  
  // dopocteni fint
  truss2d_fint (xze, r, s, A, fint);
  (*internalForces)[0] = fint[3];
}

void XALM_interface :: lineq_solve (Dvctr *X, const Dvctr *R)
{
  (*X)[0] = (*R)[0] / k;
}

void XALM_interface :: mtlb (void)
{
  neq = 1;
  
  // souradnice uzlu
  double xz[2][2] = {{0,0},{4,3}};
  
  // kodova cisla prvku
  //lm = [1,2,3,4];
  
  // element nodes
  int en[] {0,1};
  
  // typ matice tuhosti (0-elasticka, 1-tecna)
  //int mtype = 0;
  
  // mode (0...load displacement plot, 1...convergence plot)
  int mode = 0;
  
  // modul pruznosti
  double E = 1.0;
     
  // plocha
  double A = 1;
  
  
  
  // pocet predepsanych
  int pneq = 3;
  
  // vektor zatizeni
  double f[neq] = {0};
  
  // prirustek zatizeni
  double df[neq] = {0};
  
  // matice tuhosti
  //double k[neq][neq] = {0};
  
  // vektor neznamych posunuti
  double ru[neq] = {0};
  
  // vektor vsech posunu
  double r[neq+pneq] = {0};
  
  // residual vector
  double res[neq] = {0};
  
  // vektor vnitrnich sil = napeti
  //double s[1] = {0};
  
  //vektor pro ukladani historie chyby (graf konvergence)
  //double rh;
  
  
  f[0] = 0.00;
  df[0]=-0.01;
  double limit = 1.e-6;
  double eps, s;
  double fint[4];
  
  FILE *file = fopen("graf1.mumech.dat", "w");
  // kontrola posunem
  fprintf (file, "%g  %g\n", 0.0, 0.0);
  
  double xze[] = { xz[en[0]][0], xz[en[0]][1], xz[en[1]][0], xz[en[1]][1] };
  for (int i = 1; i<=35; i++) {
    r[0] = r[1] = r[2] = 0.0;
    r[3] = -i*0.2;
    truss2d_postpro (xze, E, r, eps, s);
    
    // dopocteni fint
    truss2d_fint (xze, r, s, A, fint);
    if (mode == 0) // load -displacement plot
      fprintf (file, "%g %g\n", r[3], fint[3]);
  }
  fclose(file);
  
  
  r[3] = 0;
  
  file = fopen("graf2.mumech.dat", "w");
  // kontrola posunem
  fprintf (file, "%lf   %lf\n", 0.0, 0.0);
  
  //cyklus prez prirustky zatizeni
  for (int i = 0; i<5; i++) {
    
    res[0] = df[0];
    double klin[4][4];
    truss2d (xze, E*A, klin);
    //                                   //inicializace citace iteraci
    double nite = 0;
    //rh=[];
    // iterace rovnovahy
    while (norm(res[0]) > limit) {
      nite = nite+1;
      //sestaveni linearizovane podminky rovnovahy
      //       ks = truss2d_initialstress(xze, A, s(1));
      //       kg = truss2d_initialdeformation (xze, r(lm(1,:)), E*A);
      //       if mtype == 1
      //         k = klin +kg+ks;
      //       else
      //         k = klin;
      //       end
      
      // vyresime soustavu klin * ru = res kde neznama je ru
      // ru[0] = klin[3][3] \ res[0];
      
      ru[0] = res[0] / klin[3][3];
      
      r[3] += ru[0];
      // dopocteni deformaci a napeti
      truss2d_postpro (xze, E, r, eps, s);
      
      // dopocteni fint
      truss2d_fint (xze, r, s, A, fint);
      
      res[0] = -1.0*(fint[3] - f[0] - df[0]);
      //rh(nite)=abs(res(1)); %remember for convergence plot
    }
    
    if (mode == 0) // load -displacement plot
      fprintf (file, "%g  %g\n", r[3], fint[3]);
    else
      ; //semilogy (rh)
    
    //     [nite; fint(4)]   ?????
    f[0] = f[0] + df[0];
    
  }
  
  fclose(file);
  
}

} // namespace xalm