#include /* vector.c support routines for defining vector arrays it is a little funky due to the fact that I want to be able to interdigitate things as well EXAMPLE SYNTAX: vector {u,v} 100 x 1: u'=f(u(1),v(1))+du*(u(2)-u(1)) v'=g(u(1),v(1))+dv*(v(2)-v(1)) ; 2--nvec-1: u'=f(u(x),v(x))+du*(u(x+1)-2*u(x)+u(x-1)) v'=g(u(x),v(x))+dv*(v(x+1)-2*v(x)+v(x-1)) ; nvec: u'=f(u(nvec),v(nvec))+du*(u(nvec-1)-u(nvec)) v'=g(u(nvec),v(nvec))+dv*(v(nvec-1)-v(nvec)) ; init u=exp(-(x-nvec/2)^2) init v=1 endvector this is stuck in a structure for this particular vector group This is the structure for a component of the vector over a particular range of subscripts. If index<=0 then it refers to nvec+index typedef struct { int ncomp; // number of components int inlo,inhi; // indices 1 is first 0 mean to the end, -1 end-1 etc char **rhs; // strings for each component right-hand side int **comprhs; // compiled right-hand sides } VECRHS; Each of these lies in the bigger structure typedef struct { int length,ncomp,nrhs; // length of vectors , number of components, // number of subsets - dimension of VECRHS double *v,*vold; //actual values of the vector // ordered as u0,v0,u1,v1,....,un,vn uv componenst char **names,x[10]; // names of components and index name VECRHS *r; // holds the component wide rhs's } XPPVEC; so for the present example here is the structure XPPVEC xppvec[0]: 100 2 3 v[200],vold[200] u v r[0] : 2 1 1 rhs[0]=f(u(1),... rhs[1]=g(u(1),... comprhs[0]= ... comprhs[1]= ... r[1] : 2 2 -1 rhs[0]= rhs[1]=... ... r[2] : 2 0 0 rhs[0]=... rhs[1]=... */ #include #include #ifndef WCTYPE #include #else #include #endif typedef struct { int ncomp; // number of components int inlo,inhi,flag; // indices 1 is first 0 mean to the end, -1 end-1 etc /* flag determines which indices are there */ char **rhs; // strings for each component right-hand side int **comprhs; // compiled right-hand sides int type; /* fixed or ode */ } VECRHS; typedef struct { int length,ncomp,nrhs; // length of vectors , number of components, // number of subsets - dimension of VECRHS double *v,*vold; //actual values of the vector // ordered as u0,v0,u1,v1,....,un,vn uv componenst char **names,x[10]; // names of components and index name char **ics; // holds the initial data strings VECRHS r[50]; // holds the component wide rhs's } XPPVEC; #define VEC_ERROR -1 #define VEC_OK 1 #define END_VECTOR 2 #define END_VRHS 3 #define VEC_RHS 4 #define VEC_DOMAIN 5 #define VEC_INIT 6 #define VEC_NOOP 50 #define MAXVEC 100 #define VRHS xppvec[n].r[ir] XPPVEC xppvec[MAXVEC]; int nxppvec=0; main() { FILE *fp; char bob[256]; int err,i; fp=fopen("vectst.ode","r"); if(fp==NULL){ printf(" Not found \n"); exit(0); } while(!feof(fp)){ fgets(bob,255,fp); if(bob[0]=='v'&&bob[1]=='e'&&bob[2]=='c'&&bob[3]=='t'){ err=read_vec(fp,bob); if(err==VEC_ERROR){ printf("VECTOR ERROR \n"); fclose(fp); exit(0); } } else { printf("%s",bob); } } fclose(fp); printf("----------------------- FILE CLOSED ------------\n\n"); for(i=0;i %s \n",xppvec[n].names[j],xppvec[n].r[i].rhs[j]); } } read_vec(FILE *fp,char *line1) { char names[MAXVEC][10]; char index[10]; char rhs[256],lhs[50]; int type,lo,hi; int error=0; int n=nxppvec; char bob[256]; int nlen,ncomp; int command; int i,cnt,ir; vec_first_line(line1,names,&nlen,index,&ncomp); xppvec[n].ncomp=ncomp; strcpy(xppvec[n].x,index); xppvec[n].length=nlen; xppvec[n].names=(char **)malloc(ncomp*sizeof(char *)); xppvec[n].ics=(char **)malloc(ncomp*sizeof(char *)); for(i=0;i=0){ printf("vector domain defined more than once \n"); return VEC_ERROR; } VRHS.inlo=lo; VRHS.inhi=hi; VRHS.flag=type; break; case VEC_RHS: /* figure out which one */ for(i=0;i0)break; } if(isalpha(c)){ printf("Illegal dimension of vector \n %s \n",line1); return(0); } j++; if(j==nl){ printf("Illegal syntax in vector declaration \n %s \n",line1); return(0); } } num[i]=0; printf("Dimension=%d \n",atoi(num)); i=0; while(1){ c=line1[j]; if(isspace(c)){ if(i>0)break; } index[i]=c; i++; j++; if(j==nl)break; } if(i==0){ printf("Need to define an index in vector \n %s \n",line1); return(0); } index[i]=0; printf("index=%s \n",index); *nlen=atoi(num); *nc=ncomp; } grab_vec_ics(char *s,char *lhs,char *rhs) { int i=0,n=strlen(s); int j=0; char c; while(1){ c=s[i]; if(c=='='){ lhs[j]=0; i++; for(j=i;j=n)break; } return 0; } get_vect_domain(char *s,int *i1,int *i2) { int i=0,j=0; int dash=0; int flag=0; char c; char *z; char lo[100],hi[100]; int n=strlen(s); int nm1=n-1; lo[0]=0; hi[0]=0; while(1){ c=s[i]; if((c=='-')&&(i=n){ if(dash==0) lo[j]=0; else hi[j]=0; break; } } /* check to see if there are any */ if(strlen(lo)>0)flag=1; if(strlen(hi)>0)flag+=2; /* okay, now we have split them */ if(flag==1||flag==3){ z=strstr(lo,"NVEC"); /* check for the symbol nvec */ if(z==NULL) *i1=atoi(lo); else *i1=atoi(z+4); } if(dash==1){ z=strstr(hi,"NVEC"); if(z==NULL) *i2=atoi(hi); else *i2=atoi(z+4); } return(flag); /* 1 for lo 2 for hi 3 for both */ } get_rhs_lhs(char *s,char *lhs,char *rhs,int *type) { int n=strlen(s); int i=0,flag=0,j=0,ic,error=1; char c; rhs[0]=0; lhs[0]=0; de_space(s); if(strlen(s)==0)return 1; if(s[0]=='#')return 1; while(1){ c=s[i]; ic=c; if(ic==39){ *type=1; lhs[j]=0; strcpy(rhs,&s[i+2]); break; } if(c=='='){ *type=2; lhs[j]=0; strcpy(rhs,&s[i+1]); break; } lhs[j]=c; j++; i++; if(i>=n) break; } return 0; } /* support routines - in the full xpp */ de_space(s) char *s; { int n=strlen(s); int i,j=0; char ch; for(i=0;i