BeamSolver.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef BEAM_SOLVER_H
#define BEAM_SOLVER_H

#include <malloc.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include "Solver.h"

class BeamSolver : public Solver {
public:
  enum BoundaryCondition { Clamped, FixedZeroMoment, Free };
  
  BeamSolver() : E(0.), nu(0.), I(0.), h(0.), rho(0.), l(1.), M(12),
                 time(0.), steps(0), lower(Clamped), upper(Clamped), dtret(1.),
                 P(0), w(0), wdot(0), pcur(0), pold(0) {
    OutputName="-";
    CheckpointName="beam";
  }  
  virtual ~BeamSolver() {}
  
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"BeamSolver");
    RegisterAt(LocCtrl,"YoungMod",E);
    RegisterAt(LocCtrl,"PoissonRatio",nu);    
    RegisterAt(LocCtrl,"Inertia",I);
    RegisterAt(LocCtrl,"Thickness",h);
    RegisterAt(LocCtrl,"Density",rho);
    RegisterAt(LocCtrl,"Length",l);
    RegisterAt(LocCtrl,"NPoints",M);
    RegisterAt(LocCtrl,"LowerBoundary",lower);
    RegisterAt(LocCtrl,"UpperBoundary",upper);
    RegisterAt(LocCtrl,"dt",dtret);    
    RegisterAt(LocCtrl,"OutputName",OutputName);
    RegisterAt(LocCtrl,"CheckpointName",CheckpointName);
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual bool setup() {
    if (E<=0. || h<=0. || rho<=0. || l<=0. || M<=2) 
      return false;
    if (I==0.) I=h*h*h/(12.*(1.-nu*nu));

    N=M-2;
    dx=l/(M-1.);
    P=matrix(N,N);
    w=vector(M);
    wdot=vector(M);
    pcur=vector(M);
    pold=vector(M);
    
    for (register int i=0;i<N;i++)
      for (register int j=0;j<N;j++)
        P[i][j]=0.;
    for (register int i=2;i<N-2;i++) {
      P[i][i-2]=1.; P[i][i-1]=-4.; P[i][i]=6.; P[i][i+1]=-4.; P[i][i+2]=1.; 
    }
    
    if (lower == Free) {
      P[0][0]=1.;  P[0][1]=-2.; P[0][2]=1.;  
      P[1][0]=-2.; P[1][1]= 5.; P[1][2]=-4.; P[1][3]=1.;
    }
    else if (lower == FixedZeroMoment) {
      P[0][0]=5.;  P[0][1]=-4.; P[0][2]=1.;  
      P[1][0]=-4.; P[1][1]= 6.; P[1][2]=-4.; P[1][3]=1.;
    }
    else {
      P[0][0]=7.;  P[0][1]=-4.; P[0][2]=1.;  
      P[1][0]=-4.; P[1][1]= 6.; P[1][2]=-4.; P[1][3]=1.;
    }
    
    if (upper == Free) {
      P[N-2][N-4]=1.; P[N-2][N-3]=-4.; P[N-2][N-2]=5.;  P[N-2][N-1]=-2.; 
                      P[N-1][N-3]=1.;  P[N-1][N-2]=-2.; P[N-1][N-1]=1.;
    }
    else if (upper == FixedZeroMoment) {
                      P[N-2][N-4]=1.; P[N-2][N-3]=-4.; P[N-2][N-2]=6.;  P[N-2][N-1]=-4.; 
      P[N-1][N-3]=1.;  P[N-1][N-2]=-4.; P[N-1][N-1]=5.;
    }
    else {
      P[N-2][N-4]=1.; P[N-2][N-3]=-4.; P[N-2][N-2]=6.;  P[N-2][N-1]=-4.; 
                      P[N-1][N-3]=1.;  P[N-1][N-2]=-4.; P[N-1][N-1]=7.;
    }
    return true;
  }    
  
  virtual void finish() {
    free_matrix(P,N,N);
    free_vector(w); 
    free_vector(wdot); 
    free_vector(pcur); 
    free_vector(pold); 
  } 

  virtual void Initialize(double& t, double& dt) {
    if (P) 
      for (register int i=0;i<N;i++) {
        pold[i]=0.; pcur[i]=0.; w[i]=0.; wdot[i]=0.;
      }     
    t=0.;
    dt=dtret;
  }

  virtual int NSteps() { return steps; }

  virtual void Advance(double& t, double& dt) {
    if (P) {
      double tau=0.5*dt;
      double zeta=tau*(E*I)/(dx*dx*dx*dx);
      double** A=matrix(2*N,2*N);
      double* f=vector(2*N);
      
      for (register int i=0;i<2*N;i++)
        for (register int j=0;j<2*N;j++)
          A[i][j] = 0.;
      for (register int i=0;i<N;i++)
        for (register int j=0;j<N;j++)
          A[i][j] = zeta*P[i][j];
      for (register int i=0;i<N;i++) {
        A[i][N+i]=rho*h; A[N+i][i]=1.; A[N+i][N+i]=-tau; 
      }
      
      for (register int i=0;i<N;i++) {
        f[i]=rho*h*wdot[i+1]-tau*(pcur[i+1]+pold[i+1])-zeta*xyT(P[i],w+1,N);
        f[N+i]=w[i+1]+tau*wdot[i+1];
      }
      for (register int i=0;i<M;i++)
        pold[i]=pcur[i];
      
      int *perm;
      perm = lr_factor(A,2*N);
      lr_solve(A,f,perm,2*N);
      for (register int i=0;i<N;i++) {
        w[i+1]=f[i]; wdot[i+1]=f[N+i];
      }
      double w0=w[0];
      double wN=w[N+1];
      EvalBoundaryPoints();
      wdot[0]=2.*(w[0]-w0)/dt-wdot[0]; 
      wdot[N+1]=2.*(w[N+1]-wN)/dt-wdot[N+1];        
      free_matrix(A,2*N,2*N);
      free_vector(f); 
    }

    t+=dt;
    time=t;
    steps++;
    dt=dtret;
  }

  virtual void Output() {
    if (!P) return;
    if (OutputName.c_str()[0] == '-') 
      return;
    char fname[256];
    sprintf(fname,"%s_%d.txt",OutputName.c_str(),NSteps());
    std::ofstream outfile(fname, std::ios::out);
    for (register int i=0; i<NumNodes(); i++)
      outfile << i*dx << "   " << Deflection(i) << "   " << Velocity(i)  
              << "   " << CurLoad(i) << "   " << OldLoad(i) << std::endl;
    outfile.close();
  }    

  virtual void Restart(double& t, double& dt) {
    if (!P) return;
    char fname[256];
    sprintf(fname,"%s.cp",CheckpointName.c_str());
    std::ifstream ifs(fname, std::ios::in);
    int Mr;
    ifs.read((char*)&Mr,sizeof(int));
    assert(Mr==M);
    M=Mr;
    ifs.read((char*)&time,sizeof(double));
    ifs.read((char*)w,M*sizeof(double));
    ifs.read((char*)wdot,M*sizeof(double));
    ifs.read((char*)pcur,M*sizeof(double));
    ifs.read((char*)pold,M*sizeof(double));
    ifs.close();
    dt=dtret;
  }

  virtual void Checkpointing() {
    if (!P) return;
    char fname[256];
    sprintf(fname,"%s.cp",CheckpointName.c_str());
    std::ofstream ofs(fname, std::ios::out);
    ofs.write((char*)&M,sizeof(int));
    ofs.write((char*)&time,sizeof(double));
    ofs.write((char*)w,M*sizeof(double));
    ofs.write((char*)wdot,M*sizeof(double));
    ofs.write((char*)pcur,M*sizeof(double));
    ofs.write((char*)pold,M*sizeof(double));
    ofs.close();
  }

  void EvalBoundaryPoints() {
    if (!w) return;
    if (lower == Free) 
      w[0]=2.*w[1]-w[2];
    else 
      w[0]=0.;

    if (upper == Free)
      w[N+1]=2.*w[N]-w[N-1];
    else 
      w[N+1]=0;
  }
  
  void SetDeflection(double *v) {
    if (!w) return;
    for (register int i=0;i<M;i++)
      w[i]=v[i];
    EvalBoundaryPoints();
  }
  double *Deflection() const { return w; }
  inline const double& Deflection(const int& i) const { return w[i]; }
  double MaxDeflection(double& xmax, double& wmax) {
    wmax=0.;
    if (w) {
      int imax=0;
      for (register int i=0; i<NumNodes(); i++)
        if (fabs(Deflection(i))>fabs(wmax)) { 
          wmax=Deflection(i); imax=i; 
        }
      xmax=imax*dx;
    }
    return wmax;
  }
    
  void SetVelocity(double *v) {
    if (!wdot) return;
    for (register int i=0;i<M;i++)
      wdot[i]=v[i];
    if (lower == Free) 
      wdot[0]=0.;
    if (upper == Free)
      wdot[N+1]=0;
  }
  inline double *Velocity() const { return wdot; }
  inline const double& Velocity(const int& i) const { return wdot[i]; }

  inline double *OldLoad() { return pold; }
  inline double& OldLoad(const int& i) { return pold[i]; }
  inline double *CurLoad() { return pcur; }
  inline double& CurLoad(const int& i) { return pcur[i]; }
    
  inline int NumNodes() const { return M; }
  inline const double& StepSize() const { return dx; }
  inline double CurrentTime() const { return time; }
  
  void ComputeStatic() {
    if (!P) return;
    double zeta=(E*I)/(dx*dx*dx*dx);
    double** A=matrix(N,N);
    double* f=vector(N);

    for (register int i=0;i<N;i++)
      for (register int j=0;j<N;j++)
        A[i][j]=zeta*P[i][j];

    for (register int i=0;i<N;i++)
      f[i]=-pcur[i+1]; 
    for (register int i=0;i<M;i++)
      pold[i]=pcur[i];

    int *perm;
    perm = lr_factor(A,N); 
    lr_solve(A,f,perm,N); 
    
    for (register int i=0;i<N;i++) {
      w[i+1]=f[i]; wdot[i+1]=0.;
    }
    EvalBoundaryPoints();
    wdot[0]=0.; wdot[N+1]=0.;

    free_matrix(A,N,N);
    free_vector(f); 
  }
  
protected:
  void error_exit(char *msg,char *file, int line) {
    printf("%s: file %s, line %d\n",msg,file,line);
    exit(1);
  }
  
  double *vector(int n) {
    double *vec;
    vec=(double *)calloc(n,sizeof(double));
    return vec;
  }
  
  void free_vector(double *vec) {
    if(vec) free(vec);
  }
  
  double **matrix(int n,int m) {
    double **mat;
    int i,j;
    mat=(double **)calloc(n,sizeof(double *));
    if (mat) 
      for(i=0; i<n; i++)
        if (! (mat[i]=(double *)calloc(m,sizeof(double))) ) {
          for(j=0;j<i;j++) 
            free(mat[i]);
          free(mat);
          mat=(double **)NULL;
          break;
        }
    return mat;
  }
  
  void free_matrix(double **mat,int n,int m) {
    int i;
    for(i=0;i<n;i++) free(mat[i]);
    free(mat);
  }
  
  double xyT(double *x,double *y,int n) {
    int i;
    double ret; 
    ret=0.0;
    for(i=0; i<n; i++)
      ret += x[i]*y[i]; 
    return ret;
  }
  
  double **xTy(double *x,double *y,int n) {
    double **ret;
    int i,j;  
    if ( ! (ret=matrix(n,n)) )
      error_exit("Out  of memory",__FILE__,__LINE__);
    for(i=0; i<n; i++)
      for(j=0; j<n; j++)
        ret[i][j] = x[i]*y[j];
    return ret;
  }
  
  void scale(double *x,double a, int n) {
    int i;
    for(i=0; i<n; i++) 
      x[i] *= a;
    return;
  }  

  void xpay(double *x,double a,double *y,int n) {
    int i;
    for(i=0; i<n; i++)
      x[i] += a*y[i];
    return;
  }
  
  void pivot(double **a,int *perm,int n,int k) {
    int i,mi;
    double *tmp,max;
    max=fabs(a[k][k]);
    mi=k;
    for (i=k+1; i<n; i++)
      if (fabs(a[i][k]) > max ) { 
        max=fabs(a[i][k]);
        mi=i;
      }
    if (max == 0.0) 
      error_exit("Matrix singular",__FILE__,__LINE__);
    if (mi != k) {
      tmp=a[k];a[k]=a[mi];a[mi]=tmp;
      perm[k]=mi;
    }
    return;
  }
  
  int *lr_factor(double **a,int n) {
    int i,k;
    int *perm;
    perm = (int *)calloc(n,sizeof(int));
    if (!perm)
      error_exit("Out of memory",__FILE__,__LINE__);
    for(k=0; k<n; k++) 
      perm[k]=k;
    for(k=0; k<n-1; k++) {
      pivot(a,perm,n,k);
      scale(a[k]+k+1,1.0/a[k][k],n-(k+1));
      for(i=k+1; i<n; i++) 
        xpay(a[i]+k+1,-a[i][k],a[k]+k+1,n-(k+1));
    }
    return perm;
  }

  double *lr_solve(double **a,double *b,int * perm,int n) {
    int i,j;
    double l;
    for(i=0; i<n; i++) {
      j=perm[i];
      if (i!=j) {
        l=b[i];b[i]=b[j];b[j]=l;
      }
    }
    for(i=0;i<n;i++)
      b[i] = (b[i] - xyT(a[i],b,i))/a[i][i];
    for(i=n-1; i>=0; i--) {
      b[i] = b[i] - xyT(a[i]+i+1,b+i+1,n-(i+1));
    }
    return b;
  }
  
  void print_vec(double *x,int n) {
    int i;
    for(i=0; i<n; i++) 
      printf("%6.2g%c",x[i],i<n-1?' ':'\n');
  }
  
  void print_mat(double **a,int n) {
    int i,j;
    for(i=0; i<n; i++) 
      for(j=0; j<n; j++) 
        printf("%6.2g%c",a[i][j],j<n-1?' ':'\n');
  }

protected:
  double E;   // Young modulus
  double nu;  // Poisson ratio (for plate strips)
  double I;   // Moment of inertia
  double h;   // Beam thickness
  double rho; // Material density
  double l;   // Beam length
  int M, N;
  double dx, time;
  int steps;
  int lower, upper;
  double dtret;
  std::string OutputName;
  std::string CheckpointName;
  double** P; // Stiffness matrix
  double*  w; // Deflection
  double*  wdot; // Normal velocity
  double*  pcur; // Current loading vector
  double*  pold; // Previous loading vector
  ControlDevice LocCtrl;
};


#endif

---