Commit 01a0eb0d authored by François Trigaux's avatar François Trigaux
Browse files

Enabled dynamic alaysis

parent ee0f8280
[Beam]
analysis_flag = 2
nkp = 11
niter = 1
nstep = 10
ncond_pt = 2
nframe = 0
ncurv = 0
initFromFile = 0
L=60.0
dt=0.5
\ No newline at end of file
5.49178763e-10 -1.67119587e-11 -1.31517283e-11 8.30682986e-10 4.17793986e-10 5.76023831e-10
-1.67119587e-11 2.37592756e-09 2.86397280e-10 -4.44996430e-11 -1.87665262e-09 3.46902150e-11
-1.31517283e-11 2.86397280e-10 3.65937383e-08 7.59886140e-10 8.43481577e-10 3.53846314e-10
8.30682986e-10 -4.44996430e-11 7.59886140e-10 4.74790111e-08 2.26164544e-09 6.84960573e-10
4.17793986e-10 -1.87665262e-09 8.43481577e-10 2.26164544e-09 5.34979199e-08 -1.01551813e-09
5.76023831e-10 3.46902150e-11 3.53846314e-10 6.84960573e-10 -1.01551813e-09 2.91821980e-09
\ No newline at end of file
# NOTE ON THIS MAKEFILE
# As the library is not installed in a default repo, i.e. /usr/lib or /usr/local/lib, you need to give it the path to the library at runtime. This is the rpath argument.
DIR=$(shell dirname `pwd`)
default:
cp -r ../../lib ./
gcc -g -fPIC -c -o main.o main.c -I../../include
gcc -g -o gebt main.o -lgebt -L./lib
clean:
rm -f main.o
rm -f gebt
258.053 0.0 0.0 0.0 7.07839 -71.6871
0.0 258.053 0.0 -7.07839 0.0 0.0
0.0 0.0 258.053 71.6871 0.0 0.0
0.0 -7.07839 71.6871 48.59 0.0 0.0
7.07839 0.0 0.0 0.0 2.172 0.0
-71.6871 0.0 0.0 0.0 0.0 46.418
\ No newline at end of file
#include <stdlib.h>
#include <stdio.h>
#include "beam.h"
int main(int argc, char* argv[])
{
int i;
// Check the input argument
if(argc<2)
{
printf("Please specify the uniform load value as arguments\n");
printf("Example: ./gebt 100\n");
exit(EXIT_FAILURE);
}
// Initializing the Beam
printf("Initializing the beam...\n");
Beam *bm = malloc(sizeof(Beam));
initBeam(bm);
printf("...Done!\n");
// Creating an array with the uniform load for each member
printf("Creating the loads...\n");
double *loads = malloc(sizeof(double)*bm->nmemb);
for(i=0;i<bm->nmemb;i++)
{
loads[i] = atof(argv[1])*0.0;
}
loads[bm->nmemb-1] = atof(argv[1]);
beam_setLoads(bm,loads,2);
printf("...Done!\n");
char fname[32];
// Performing the analysis
printf("Performing the analysis...\n");
for (i=0;i<10;i++)
{
sprintf(fname,"output%d_",i);
printf("Global Step = %d...\n",i);
beam_analysis(bm);
beam_writeSolToFile(bm,fname);
printf("... Ok! \n");
}
return EXIT_SUCCESS;
}
# Launch the Euler-Bernoulli comparison for the beam
import numpy as np
import matplotlib.pyplot as plt
from readIni import IniReader
import subprocess as sp
ini = IniReader("Beam.ini");
nkp = ini.getValue("Beam","nkp",typeCast=int);
nmemb = nkp-1
display = 1
x = np.linspace(0,1);
xcheck = np.linspace(0,1,nkp);
##a = np.array([2.389e9, 1.524e6, 6.734e6, -3.382e7, -2.627e7, -4.736e8,
## 1.524e6, 4.334e8, -3.741e6, -2.935e5, 1.527e7, 3.835e5,
## 6.734e6, -3.741e6, 2.743e7, -4.592e5, -6.869e5, -4.742e6,
##-3.382e7, -2.935e5, -4.592e5, 2.167e7, -6.279e5, 1.430e6,
##-2.627e7, 1.527e7, -6.869e5, -6.279e5, 1.970e7, 1.209e7,
##-4.736e8, 3.835e5, -4.742e6, 1.430e6, 1.209e7, 4.406e8]);
##a = a.reshape((6,6));
##print(a);
##print(np.linalg.inv(a));
sp.run(["./gebt", "100.0"]);
for k in range(10):
for m in range(10):
data_pts = np.loadtxt("output%d_%03d.dat"%(k,m),max_rows=nkp);
data_mem = np.loadtxt("output%d_%03d.dat"%(k,m),skiprows=nkp);
pts = data_pts;
mem = data_mem;
if(display):
#plt.plot(mem[:,0]+mem[:,3],mem[:,4]+mem[:,1]);
plt.plot(pts[:,0]+pts[:,3],pts[:,5]+pts[:,2],'.-k');
#plt.plot(mem[:,0]+mem[:,3],mem[:,5]+mem[:,2],'.r');
plt.xlabel('x');
plt.ylabel('w(x)');
plt.title('%d %d'%(k,m));
plt.ylim([0,5]);
plt.pause(0.05);
plt.clf();
print("Done!");
if(display):
plt.show();
\ No newline at end of file
......@@ -17,7 +17,7 @@ SUBROUTINE Analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndistrfun,nc
& member,pt_condition,material,niter,nstep,sol_pt,sol_mb,error,ncond_mb, &
& ntimefun,frame,mb_condition,distr_fun,curvature,omega_a0,omega_a_tf, &
& v_root_a0,v_root_a_tf,simu_time,time_function,analysis_flag,init_cond, &
& nev,eigen_val,eigen_vec_pt,eigen_vec_mb)
& nev,eigen_val,eigen_vec_pt,eigen_vec_mb,initialize)
!DEC$ ATTRIBUTES DLLEXPORT :: Analysis
......@@ -67,6 +67,8 @@ INTEGER,INTENT(IN) :: analysis_flag
INTEGER,INTENT(INOUT):: nev ! upon return, nev is the number of converged eigen values might be the original nev given by the user in the inputs or nev+1 due to complex conjugate pairs
REAL(DBL),INTENT(OUT)::eigen_val(2,nev+1),eigen_vec_pt(nev+1,nkp,NDIM+NDOF_ND),eigen_vec_mb(nev+1,nelem,NDIM+ndof_el)
INTEGER,INTENT(IN) :: initialize
! Variables used in this subroutine
!---------------------------------------------------------------------------------------
INTEGER:: dof_con(nkp) ! the connecting condition for key point.
......@@ -147,37 +149,39 @@ IF(analysis_flag==2) THEN
two_divide_dt=2.0D0/time_incr !2/dt
!Obtaining the complete set of initial conditions needed for
!time marching later
!----------------------------------------
init_flag=1
!WRITE(*,*)"THE INITIAL STEP"
IF(DEBUG) WRITE(IOUT,*)"THE INITIAL STEP"
IF(niter==1) THEN ! solve the linear system for initial conditions
CALL LinearSolution(ndof_el,memb_info,v_root_a,omega_a,member,error,&
& ncond_mb,mb_condition,distr_fun,(dof_con),x,init_cond)
ELSE ! Solve the nonlinear system for initial conditions
!---------------------------------------------------
CALL NewtonRaphson(ndof_el,memb_info,v_root_a,omega_a,member,niter,error,&
& ncond_mb,mb_condition,distr_fun,(dof_con),x,init_cond)
ENDIF
IF(error/='') GOTO 9999
IF(initialize==1) THEN
!Obtaining the complete set of initial conditions needed for
!time marching later
!----------------------------------------
init_flag=1
!WRITE(*,*)"THE INITIAL STEP"
IF(DEBUG) WRITE(IOUT,*)"THE INITIAL STEP"
IF(niter==1) THEN ! solve the linear system for initial conditions
CALL LinearSolution(ndof_el,memb_info,v_root_a,omega_a,member,error,&
& ncond_mb,mb_condition,distr_fun,(dof_con),x,init_cond)
ELSE ! Solve the nonlinear system for initial conditions
!---------------------------------------------------
CALL NewtonRaphson(ndof_el,memb_info,v_root_a,omega_a,member,niter,error,&
& ncond_mb,mb_condition,distr_fun,(dof_con),x,init_cond)
ENDIF
IF(error/='') GOTO 9999
CALL ExtractElementValues(ndof_el,member,x,sol_mb(1,:,:))
PHDot=sol_mb(1,:,4:9) ! CTCabPDot, CTCabHDot
sol_mb(1,:,4:9)=init_cond(:,1:6) ! insert the initial coniditions for ui, thetai into the solution.
CALL ExtractElementValues(ndof_el,member,x,sol_mb(1,:,:))
PHDot=sol_mb(1,:,4:9) ! CTCabPDot, CTCabHDot
sol_mb(1,:,4:9)=init_cond(:,1:6) ! insert the initial coniditions for ui, thetai into the solution.
! set initial values as initial condition of x for time marching, replace Pdot, Hdot with u, theta
!----------------------------------------------------------------------------------------------
CALL InsertElementValues(ndof_el,member,x,init_cond)
! set initial values as initial condition of x for time marching, replace Pdot, Hdot with u, theta
!----------------------------------------------------------------------------------------------
CALL InsertElementValues(ndof_el,member,x,init_cond)
! obtain: 2/dt*value+\dot{value} at time t0
!-------------------------------------------------------
init_cond(:,1:6) =two_divide_dt*sol_mb(1,:,4:9)+init_cond(:,7:12) ! u, \theta
init_cond(:,7:12)=two_divide_dt*CTCabPH(niter,member,memb_info,sol_mb(1,:,4:))+PHDot ! P, H
! obtain: 2/dt*value+\dot{value} at time t0
!-------------------------------------------------------
init_cond(:,1:6) =two_divide_dt*sol_mb(1,:,4:9)+init_cond(:,7:12) ! u, \theta
init_cond(:,7:12)=two_divide_dt*CTCabPH(niter,member,memb_info,sol_mb(1,:,4:))+PHDot ! P, H
ENDIF
init_flag=2
......
......@@ -6,7 +6,7 @@ SUBROUTINE fortran_analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndis
& member,pt_condition,material,niter,nstep,sol_pt,sol_mb,error,ncond_mb, &
& ntimefun,frame,mb_condition,distr_fun,curvature,omega_a0,omega_a_tf, &
& v_root_a0,v_root_a_tf,simu_time,time_function,analysis_flag,init_cond, &
& nev,eigen_val,eigen_vec_pt,eigen_vec_mb) bind(c,name="fortran_analysis")
& nev,eigen_val,eigen_vec_pt,eigen_vec_mb,initialize) bind(c,name="fortran_analysis")
use,intrinsic :: iso_c_binding
USE InternalData
......@@ -54,6 +54,7 @@ SUBROUTINE fortran_analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndis
TYPE(C_PTR),INTENT(OUT)::eigen_val,eigen_vec_pt,eigen_vec_mb
CHARACTER(*),INTENT(OUT) ::error(:) ! error message
INTEGER(KIND=C_INT),INTENT(IN), VALUE :: initialize
CHARACTER(300) :: f_error=''
INTEGER :: i
......@@ -114,7 +115,7 @@ SUBROUTINE fortran_analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndis
& f_member,f_pt_condition,f_material,niter,nstep,f_sol_pt,f_sol_mb,f_error,ncond_mb, &
& ntimefun,f_frame,f_mb_condition,f_distr_fun,f_curvature,omega_a0,omega_a_tf, &
& v_root_a0,v_root_a_tf,simu_time,f_time_function,analysis_flag,f_init_cond, &
& nev,f_eigen_val,f_eigen_vec_pt,f_eigen_vec_mb)
& nev,f_eigen_val,f_eigen_vec_pt,f_eigen_vec_mb,initialize)
! --------------------------------------------------
! Transfer back to C pointers
......
......@@ -41,6 +41,10 @@ void initBeam(Beam *bm)
SET_OR_DEFAULT(&(bm->niter),%d,fileBeam,section,"niter",100);
SET_OR_DEFAULT(&(bm->nstep),%d,fileBeam,section,"nstep",1);
if(bm->analysis_flag>0){
SET(&(bm->dt),%lf,fileBeam,section,"dt");
}
bm->nmemb = bm->nkp-1; // specific for beams.
bm->nmate = bm->nkp; // 1 material per point
bm->ncond_mb = bm->nmemb; // each member has one condition
......@@ -49,6 +53,8 @@ void initBeam(Beam *bm)
bm->ntimefun = 0;
bm->initialize = 1;
if(bm->analysis_flag==0){
bm->ndof_el=NDOF_ND;
for (i = 0; i < NDIM; i++)
......@@ -225,7 +231,7 @@ void initBeam(Beam *bm)
// This must be initialize but has no consequence as no time-functions are used
bm->simu_time[0] = 0.0;
bm->simu_time[1] = 1.0;
bm->simu_time[1] = bm->nstep * bm->dt;
ALLOCATE1(bm->time_function,struct TimeFunction,bm->ntimefun);
if (bm->ntimefun>0)
......@@ -546,7 +552,7 @@ void initBeamFromGebtFile(Beam *bm)
}
void beam_writeSolToFile(Beam *bm)
void beam_writeSolToFile(Beam *bm,char* fileName)
{
int i,j,k;
char fname[32];
......@@ -554,7 +560,7 @@ void beam_writeSolToFile(Beam *bm)
for (k=0;k<bm->nstep;k++){
sprintf(fname,"output%03d.dat",k);
sprintf(fname,"%s%03d.dat",fileName,k);
fid = fopen(fname,"w");
if (fid==NULL)
{
......@@ -605,6 +611,7 @@ void beam_setLoads(Beam *bm, double *loads, int load_no)
void beam_analysis(Beam *bm)
{
int i;
char *error;
int nev;
......@@ -629,12 +636,19 @@ void beam_analysis(Beam *bm)
double *rav_eigen_vec_pt= malloc(sizeof(double));
double *rav_eigen_vec_mb= malloc(sizeof(double));
// The material tensor is modified by fortran_analysis. Thus, we will give it a copy
double *material_tmp = (double*) malloc(sizeof(double) * bm->nmate * (bm->ndof_el - NSTRN) * NSTRN);
for(i=0;i<bm->nmate * (bm->ndof_el - NSTRN) * NSTRN;i++)
{
material_tmp[i] = rav_material[i];
}
// Call the fortran function
fortran_analysis(bm->nkp,bm->nelem,bm->ndof_el,bm->nmemb,bm->ncond_pt,bm->nmate, bm->nframe,bm->ndistrfun,
bm->ncurv,&rav_coord,&rav_member,&(bm->pt_condition),&rav_material,bm->niter,bm->nstep,
bm->ncurv,&rav_coord,&rav_member,&(bm->pt_condition),&material_tmp,bm->niter,bm->nstep,
&rav_sol_pt,&rav_sol_mb,error,bm->ncond_mb,bm->ntimefun,&rav_frame,&(bm->mb_condition),&rav_distr_fun,
&rav_curvature,bm->omega_a0,bm->omega_a_tf,bm->v_root_a0,bm->v_root_a_tf,bm->simu_time,
&bm->time_function,bm->analysis_flag,&rav_init_cond,&nev,&rav_eigen_val,&rav_eigen_vec_pt,&rav_eigen_vec_mb);
&bm->time_function,bm->analysis_flag,&rav_init_cond,&nev,&rav_eigen_val,&rav_eigen_vec_pt,&rav_eigen_vec_mb,bm->initialize);
// Transform 1D pointers to 2D pointers if necessary
UNRAVEL2(bm,coord,bm->nkp,NDIM,double);
......@@ -650,14 +664,82 @@ void beam_analysis(Beam *bm)
//UNRAVEL3(eigen_vec_pt,nev+1,nkp,NDIM+NDOF_ND,double);
//UNRAVEL3(eigen_vec_mb,nev+1,nelem,NDIM+ndof_el,double);
// Set initialize to zero to perform dynamic analysis without redoing the initialization step
bm->initialize = 0;
free(rav_eigen_val);
free(rav_eigen_vec_mb);
free(rav_eigen_vec_pt);
free(material_tmp);
}
/*
This updates the initial condition, but when the analysis does its initial step, it does not work
void beam_updateInitialConditions(Beam *bm)
{
int i,j;
for (i=0;i<bm->nelem;i++)
{
// copy u_dot, theta_dot from the initial conditions
for (j=0;j<NSTRN;j++)
{
bm->init_cond[i][j+NSTRN] = bm->init_cond[i][j];
}
// copy u, theta from the solution member
for (j=0;j<NSTRN;j++)
{
bm->init_cond[i][j] = bm->sol_mb[bm->nstep-1][i][j+NDIM];
}
}
}
*/
#define WRITEF3(ptr,fid,type,n1,n2,n3) for (i_=0; i_<(n1); i_++){\
for (j_=0; j_<(n2); j_++){\
for (k_=0; k_<(n3); k_++){\
fprintf(fid,#type " ",ptr[i_][j_][k_]);}\
fprintf(fid,"\n");}\
fprintf(fid,"\n");}
#define WRITEF2(ptr,fid,type,n1,n2) for (i_=0; i_<(n1); i_++){\
for (j_=0; j_<(n2); j_++){\
fprintf(fid,#type " ",ptr[i_][j_]);}\
fprintf(fid,"\n");}
#define WRITEF1(ptr,fid,type,n1) for (i_=0; i_<(n1); i_++){\
fprintf(fid,#type " ",ptr[i_][j_]);\
fprintf(fid,"\n");}
void beam_writeBeamToFile(Beam *bm,char *fname)
{
int i,j,k;
FILE *fid = fopen(fname,"w");
fprintf(fid,"nkp = %d",bm->nkp);
fprintf(fid,"%d",bm->nelem);
fprintf(fid,"%d",bm->ndof_el);
fprintf(fid,"%d",bm->nmemb);
fprintf(fid,"%d",bm->ncond_pt);
fprintf(fid,"%d",bm->nmate);
fprintf(fid,"%d",bm->nframe);
fprintf(fid,"%d",bm->ndistrfun);
fprintf(fid,"%d",bm->ncurv);
WRITEF2(bm->coord, fid, %f, bm->nkp, NDIM);
WRITEF2(bm->member, fid, %d, bm->nmemb, MEMB_CONST);
WRITEF3(bm->material, fid, %f, bm->nmate, bm->ndof_el - NSTRN, NSTRN);
WRITEF3(bm->frame, fid, %f, bm->nframe, NDIM, NDIM);
WRITEF2(bm->distr_fun,fid, %f, bm->ndistrfun, NSTRN);
WRITEF2(bm->curvature,fid, %f, bm->ncurv, NDIM);
fclose(fid);
}
void freeBeam(Beam *bm)
{
......
......@@ -75,7 +75,7 @@ obj->VEC[i_][j_][k_] = rav_##VEC[k_*(nx)*(ny)+j_*(nx)+i_];\
}}}\
free(rav_##VEC)
extern int i_,j_,k_;
struct PrescriInf{
......@@ -144,14 +144,17 @@ typedef struct
// Personal add
double L; // Beam Length
int initFromFile;
double dt;
int initialize; // flag indicating if an initialisation step is required to find the compatible state variables
} Beam;
void initBeam(Beam *bm);
void beam_setLoads(Beam *bm, double *loads, int load_no);
void beam_writeSolToFile(Beam *bm);
void beam_writeSolToFile(Beam *bm, char *fileName);
void beam_analysis(Beam *bm);
//void beam_updateInitialConditions(Beam *bm);
void freeBeam(Beam *bm);
......@@ -160,7 +163,7 @@ void fortran_analysis(int nkp,int nelem,int ndof_el,int nmemb,int ncond_pt,int n
int **member,struct PrescriInf **pt_condition, double **material, int niter, int nstep, double **sol_pt, double **sol_mb,char *error, int ncond_mb,
int ntimefun,double **frame, struct PrescriInf **mb_condition,double **distr_fun, double **curvature, double omega_a0[], int omega_a_tf[],
double v_root_a0[], int v_root_a_tf[], double simu_time[], struct TimeFunction **time_function, int analysis_flag, double **init_cond,
int *nev, double **eigen_val, double **eigen_vec_pt,double **eigen_vec_mb);
int *nev, double **eigen_val, double **eigen_vec_pt,double **eigen_vec_mb, int initialize);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment