00001
00002
00003
00004
00005
00006 #ifndef BEAM_SOLVER_H
00007 #define BEAM_SOLVER_H
00008
00016 #include <malloc.h>
00017 #include <stdio.h>
00018 #include <math.h>
00019 #include <stdlib.h>
00020 #include <assert.h>
00021 #include "Solver.h"
00022
00043 class BeamSolver : public Solver {
00044 public:
00045 enum BoundaryCondition { Clamped, FixedZeroMoment, Free };
00046
00047 BeamSolver() : E(0.), nu(0.), I(0.), h(0.), rho(0.), l(1.), M(12),
00048 time(0.), steps(0), lower(Clamped), upper(Clamped), dtret(1.),
00049 P(0), w(0), wdot(0), pcur(0), pold(0) {
00050 OutputName="-";
00051 CheckpointName="beam";
00052 }
00053 virtual ~BeamSolver() {}
00054
00055 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00056 LocCtrl = Ctrl.getSubDevice(prefix+"BeamSolver");
00057 RegisterAt(LocCtrl,"YoungMod",E);
00058 RegisterAt(LocCtrl,"PoissonRatio",nu);
00059 RegisterAt(LocCtrl,"Inertia",I);
00060 RegisterAt(LocCtrl,"Thickness",h);
00061 RegisterAt(LocCtrl,"Density",rho);
00062 RegisterAt(LocCtrl,"Length",l);
00063 RegisterAt(LocCtrl,"NPoints",M);
00064 RegisterAt(LocCtrl,"LowerBoundary",lower);
00065 RegisterAt(LocCtrl,"UpperBoundary",upper);
00066 RegisterAt(LocCtrl,"dt",dtret);
00067 RegisterAt(LocCtrl,"OutputName",OutputName);
00068 RegisterAt(LocCtrl,"CheckpointName",CheckpointName);
00069 }
00070 virtual void register_at(ControlDevice& Ctrl) {
00071 register_at(Ctrl, "");
00072 }
00073
00074 virtual bool setup() {
00075 if (E<=0. || h<=0. || rho<=0. || l<=0. || M<=2)
00076 return false;
00077 if (I==0.) I=h*h*h/(12.*(1.-nu*nu));
00078
00079 N=M-2;
00080 dx=l/(M-1.);
00081 P=matrix(N,N);
00082 w=vector(M);
00083 wdot=vector(M);
00084 pcur=vector(M);
00085 pold=vector(M);
00086
00087 for (register int i=0;i<N;i++)
00088 for (register int j=0;j<N;j++)
00089 P[i][j]=0.;
00090 for (register int i=2;i<N-2;i++) {
00091 P[i][i-2]=1.; P[i][i-1]=-4.; P[i][i]=6.; P[i][i+1]=-4.; P[i][i+2]=1.;
00092 }
00093
00094 if (lower == Free) {
00095 P[0][0]=1.; P[0][1]=-2.; P[0][2]=1.;
00096 P[1][0]=-2.; P[1][1]= 5.; P[1][2]=-4.; P[1][3]=1.;
00097 }
00098 else if (lower == FixedZeroMoment) {
00099 P[0][0]=5.; P[0][1]=-4.; P[0][2]=1.;
00100 P[1][0]=-4.; P[1][1]= 6.; P[1][2]=-4.; P[1][3]=1.;
00101 }
00102 else {
00103 P[0][0]=7.; P[0][1]=-4.; P[0][2]=1.;
00104 P[1][0]=-4.; P[1][1]= 6.; P[1][2]=-4.; P[1][3]=1.;
00105 }
00106
00107 if (upper == Free) {
00108 P[N-2][N-4]=1.; P[N-2][N-3]=-4.; P[N-2][N-2]=5.; P[N-2][N-1]=-2.;
00109 P[N-1][N-3]=1.; P[N-1][N-2]=-2.; P[N-1][N-1]=1.;
00110 }
00111 else if (upper == FixedZeroMoment) {
00112 P[N-2][N-4]=1.; P[N-2][N-3]=-4.; P[N-2][N-2]=6.; P[N-2][N-1]=-4.;
00113 P[N-1][N-3]=1.; P[N-1][N-2]=-4.; P[N-1][N-1]=5.;
00114 }
00115 else {
00116 P[N-2][N-4]=1.; P[N-2][N-3]=-4.; P[N-2][N-2]=6.; P[N-2][N-1]=-4.;
00117 P[N-1][N-3]=1.; P[N-1][N-2]=-4.; P[N-1][N-1]=7.;
00118 }
00119 return true;
00120 }
00121
00122 virtual void finish() {
00123 free_matrix(P,N,N);
00124 free_vector(w);
00125 free_vector(wdot);
00126 free_vector(pcur);
00127 free_vector(pold);
00128 }
00129
00130 virtual void Initialize(double& t, double& dt) {
00131 if (P)
00132 for (register int i=0;i<N;i++) {
00133 pold[i]=0.; pcur[i]=0.; w[i]=0.; wdot[i]=0.;
00134 }
00135 t=0.;
00136 dt=dtret;
00137 }
00138
00139 virtual int NSteps() { return steps; }
00140
00141 virtual void Advance(double& t, double& dt) {
00142 if (P) {
00143 double tau=0.5*dt;
00144 double zeta=tau*(E*I)/(dx*dx*dx*dx);
00145 double** A=matrix(2*N,2*N);
00146 double* f=vector(2*N);
00147
00148 for (register int i=0;i<2*N;i++)
00149 for (register int j=0;j<2*N;j++)
00150 A[i][j] = 0.;
00151 for (register int i=0;i<N;i++)
00152 for (register int j=0;j<N;j++)
00153 A[i][j] = zeta*P[i][j];
00154 for (register int i=0;i<N;i++) {
00155 A[i][N+i]=rho*h; A[N+i][i]=1.; A[N+i][N+i]=-tau;
00156 }
00157
00158 for (register int i=0;i<N;i++) {
00159 f[i]=rho*h*wdot[i+1]-tau*(pcur[i+1]+pold[i+1])-zeta*xyT(P[i],w+1,N);
00160 f[N+i]=w[i+1]+tau*wdot[i+1];
00161 }
00162 for (register int i=0;i<M;i++)
00163 pold[i]=pcur[i];
00164
00165 int *perm;
00166 perm = lr_factor(A,2*N);
00167 lr_solve(A,f,perm,2*N);
00168 for (register int i=0;i<N;i++) {
00169 w[i+1]=f[i]; wdot[i+1]=f[N+i];
00170 }
00171 double w0=w[0];
00172 double wN=w[N+1];
00173 EvalBoundaryPoints();
00174 wdot[0]=2.*(w[0]-w0)/dt-wdot[0];
00175 wdot[N+1]=2.*(w[N+1]-wN)/dt-wdot[N+1];
00176 free_matrix(A,2*N,2*N);
00177 free_vector(f);
00178 }
00179
00180 t+=dt;
00181 time=t;
00182 steps++;
00183 dt=dtret;
00184 }
00185
00186 virtual void Output() {
00187 if (!P) return;
00188 if (OutputName.c_str()[0] == '-')
00189 return;
00190 char fname[256];
00191 sprintf(fname,"%s_%d.txt",OutputName.c_str(),NSteps());
00192 std::ofstream outfile(fname, std::ios::out);
00193 for (register int i=0; i<NumNodes(); i++)
00194 outfile << i*dx << " " << Deflection(i) << " " << Velocity(i)
00195 << " " << CurLoad(i) << " " << OldLoad(i) << std::endl;
00196 outfile.close();
00197 }
00198
00199 virtual void Restart(double& t, double& dt) {
00200 if (!P) return;
00201 char fname[256];
00202 sprintf(fname,"%s.cp",CheckpointName.c_str());
00203 std::ifstream ifs(fname, std::ios::in);
00204 int Mr;
00205 ifs.read((char*)&Mr,sizeof(int));
00206 assert(Mr==M);
00207 M=Mr;
00208 ifs.read((char*)&time,sizeof(double));
00209 ifs.read((char*)w,M*sizeof(double));
00210 ifs.read((char*)wdot,M*sizeof(double));
00211 ifs.read((char*)pcur,M*sizeof(double));
00212 ifs.read((char*)pold,M*sizeof(double));
00213 ifs.close();
00214 dt=dtret;
00215 }
00216
00217 virtual void Checkpointing() {
00218 if (!P) return;
00219 char fname[256];
00220 sprintf(fname,"%s.cp",CheckpointName.c_str());
00221 std::ofstream ofs(fname, std::ios::out);
00222 ofs.write((char*)&M,sizeof(int));
00223 ofs.write((char*)&time,sizeof(double));
00224 ofs.write((char*)w,M*sizeof(double));
00225 ofs.write((char*)wdot,M*sizeof(double));
00226 ofs.write((char*)pcur,M*sizeof(double));
00227 ofs.write((char*)pold,M*sizeof(double));
00228 ofs.close();
00229 }
00230
00231 void EvalBoundaryPoints() {
00232 if (!w) return;
00233 if (lower == Free)
00234 w[0]=2.*w[1]-w[2];
00235 else
00236 w[0]=0.;
00237
00238 if (upper == Free)
00239 w[N+1]=2.*w[N]-w[N-1];
00240 else
00241 w[N+1]=0;
00242 }
00243
00244 void SetDeflection(double *v) {
00245 if (!w) return;
00246 for (register int i=0;i<M;i++)
00247 w[i]=v[i];
00248 EvalBoundaryPoints();
00249 }
00250 double *Deflection() const { return w; }
00251 inline const double& Deflection(const int& i) const { return w[i]; }
00252 double MaxDeflection(double& xmax, double& wmax) {
00253 wmax=0.;
00254 if (w) {
00255 int imax=0;
00256 for (register int i=0; i<NumNodes(); i++)
00257 if (fabs(Deflection(i))>fabs(wmax)) {
00258 wmax=Deflection(i); imax=i;
00259 }
00260 xmax=imax*dx;
00261 }
00262 return wmax;
00263 }
00264
00265 void SetVelocity(double *v) {
00266 if (!wdot) return;
00267 for (register int i=0;i<M;i++)
00268 wdot[i]=v[i];
00269 if (lower == Free)
00270 wdot[0]=0.;
00271 if (upper == Free)
00272 wdot[N+1]=0;
00273 }
00274 inline double *Velocity() const { return wdot; }
00275 inline const double& Velocity(const int& i) const { return wdot[i]; }
00276
00277 inline double *OldLoad() { return pold; }
00278 inline double& OldLoad(const int& i) { return pold[i]; }
00279 inline double *CurLoad() { return pcur; }
00280 inline double& CurLoad(const int& i) { return pcur[i]; }
00281
00282 inline int NumNodes() const { return M; }
00283 inline const double& StepSize() const { return dx; }
00284 inline double CurrentTime() const { return time; }
00285
00286 void ComputeStatic() {
00287 if (!P) return;
00288 double zeta=(E*I)/(dx*dx*dx*dx);
00289 double** A=matrix(N,N);
00290 double* f=vector(N);
00291
00292 for (register int i=0;i<N;i++)
00293 for (register int j=0;j<N;j++)
00294 A[i][j]=zeta*P[i][j];
00295
00296 for (register int i=0;i<N;i++)
00297 f[i]=-pcur[i+1];
00298 for (register int i=0;i<M;i++)
00299 pold[i]=pcur[i];
00300
00301 int *perm;
00302 perm = lr_factor(A,N);
00303 lr_solve(A,f,perm,N);
00304
00305 for (register int i=0;i<N;i++) {
00306 w[i+1]=f[i]; wdot[i+1]=0.;
00307 }
00308 EvalBoundaryPoints();
00309 wdot[0]=0.; wdot[N+1]=0.;
00310
00311 free_matrix(A,N,N);
00312 free_vector(f);
00313 }
00314
00315 protected:
00316 void error_exit(const char *msg,const char *file, int line) {
00317 printf("%s: file %s, line %d\n",msg,file,line);
00318 exit(1);
00319 }
00320
00321 double *vector(int n) {
00322 double *vec;
00323 vec=(double *)calloc(n,sizeof(double));
00324 return vec;
00325 }
00326
00327 void free_vector(double *vec) {
00328 if(vec) free(vec);
00329 }
00330
00331 double **matrix(int n,int m) {
00332 double **mat;
00333 int i,j;
00334 mat=(double **)calloc(n,sizeof(double *));
00335 if (mat)
00336 for(i=0; i<n; i++)
00337 if (! (mat[i]=(double *)calloc(m,sizeof(double))) ) {
00338 for(j=0;j<i;j++)
00339 free(mat[i]);
00340 free(mat);
00341 mat=(double **)NULL;
00342 break;
00343 }
00344 return mat;
00345 }
00346
00347 void free_matrix(double **mat,int n,int m) {
00348 int i;
00349 for(i=0;i<n;i++) free(mat[i]);
00350 free(mat);
00351 }
00352
00353 double xyT(double *x,double *y,int n) {
00354 int i;
00355 double ret;
00356 ret=0.0;
00357 for(i=0; i<n; i++)
00358 ret += x[i]*y[i];
00359 return ret;
00360 }
00361
00362 double **xTy(double *x,double *y,int n) {
00363 double **ret;
00364 int i,j;
00365 if ( ! (ret=matrix(n,n)) )
00366 error_exit("Out of memory",__FILE__,__LINE__);
00367 for(i=0; i<n; i++)
00368 for(j=0; j<n; j++)
00369 ret[i][j] = x[i]*y[j];
00370 return ret;
00371 }
00372
00373 void scale(double *x,double a, int n) {
00374 int i;
00375 for(i=0; i<n; i++)
00376 x[i] *= a;
00377 return;
00378 }
00379
00380 void xpay(double *x,double a,double *y,int n) {
00381 int i;
00382 for(i=0; i<n; i++)
00383 x[i] += a*y[i];
00384 return;
00385 }
00386
00387 void pivot(double **a,int *perm,int n,int k) {
00388 int i,mi;
00389 double *tmp,max;
00390 max=fabs(a[k][k]);
00391 mi=k;
00392 for (i=k+1; i<n; i++)
00393 if (fabs(a[i][k]) > max ) {
00394 max=fabs(a[i][k]);
00395 mi=i;
00396 }
00397 if (max == 0.0)
00398 error_exit("Matrix singular",__FILE__,__LINE__);
00399 if (mi != k) {
00400 tmp=a[k];a[k]=a[mi];a[mi]=tmp;
00401 perm[k]=mi;
00402 }
00403 return;
00404 }
00405
00406 int *lr_factor(double **a,int n) {
00407 int i,k;
00408 int *perm;
00409 perm = (int *)calloc(n,sizeof(int));
00410 if (!perm)
00411 error_exit("Out of memory",__FILE__,__LINE__);
00412 for(k=0; k<n; k++)
00413 perm[k]=k;
00414 for(k=0; k<n-1; k++) {
00415 pivot(a,perm,n,k);
00416 scale(a[k]+k+1,1.0/a[k][k],n-(k+1));
00417 for(i=k+1; i<n; i++)
00418 xpay(a[i]+k+1,-a[i][k],a[k]+k+1,n-(k+1));
00419 }
00420 return perm;
00421 }
00422
00423 double *lr_solve(double **a,double *b,int * perm,int n) {
00424 int i,j;
00425 double l;
00426 for(i=0; i<n; i++) {
00427 j=perm[i];
00428 if (i!=j) {
00429 l=b[i];b[i]=b[j];b[j]=l;
00430 }
00431 }
00432 for(i=0;i<n;i++)
00433 b[i] = (b[i] - xyT(a[i],b,i))/a[i][i];
00434 for(i=n-1; i>=0; i--) {
00435 b[i] = b[i] - xyT(a[i]+i+1,b+i+1,n-(i+1));
00436 }
00437 return b;
00438 }
00439
00440 void print_vec(double *x,int n) {
00441 int i;
00442 for(i=0; i<n; i++)
00443 printf("%6.2g%c",x[i],i<n-1?' ':'\n');
00444 }
00445
00446 void print_mat(double **a,int n) {
00447 int i,j;
00448 for(i=0; i<n; i++)
00449 for(j=0; j<n; j++)
00450 printf("%6.2g%c",a[i][j],j<n-1?' ':'\n');
00451 }
00452
00453 protected:
00454 double E;
00455 double nu;
00456 double I;
00457 double h;
00458 double rho;
00459 double l;
00460 int M, N;
00461 double dx, time;
00462 int steps;
00463 int lower, upper;
00464 double dtret;
00465 std::string OutputName;
00466 std::string CheckpointName;
00467 double** P;
00468 double* w;
00469 double* wdot;
00470 double* pcur;
00471 double* pold;
00472 ControlDevice LocCtrl;
00473 };
00474
00475
00476 #endif