00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: readinp.c 2988 2010-03-01 17:46:55Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "magpario.h"
00028
00029
00030
00031 int ReadINP(GridData *gdata, char *fmesh, Vec Vdest, int col)
00032 {
00033 int i,j;
00034
00035 char msg[256];
00036 FILE *fd;
00037
00038 int a1;
00039 int id, n1, n2, n3, n4, n5;
00040 int n_ele;
00041 int t_dnv;
00042 PetscReal d1,d2,d3;
00043
00044 PetscTruth flg;
00045
00046
00047 MagparFunctionLogBegin;
00048
00049 int rank,size;
00050 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00051 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00052
00053
00054 if (Vdest!=PETSC_NULL && col>0) {
00055 ierr = VecValid(Vdest,&flg);CHKERRQ(ierr);
00056 if (flg!=PETSC_TRUE) col=0;
00057 }
00058 else {
00059 col=0;
00060 }
00061
00062 if (rank) {
00063 gdata->ln_vert=0;
00064 gdata->ln_ele=0;
00065
00066 if (col>0) {
00067 VecAssemblyBegin(Vdest);
00068 VecAssemblyEnd(Vdest);
00069 }
00070
00071 MagparFunctionLogReturn(0);
00072 }
00073
00074 ierr = PetscFOpen(PETSC_COMM_WORLD,fmesh,"r",&fd);CHKERRQ(ierr);
00075 if (!rank && !fd) {
00076 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fmesh);
00077 }
00078 PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",fmesh);
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088 fscanf(fd,"%i %i %i %i %i\n",
00089 &n1,&n2,&n3,&n4,&n5
00090 );
00091 n_ele=n2;
00092
00093 PetscPrintf(PETSC_COMM_WORLD,"Number of grid vertices: %i\n",n1);
00094 PetscPrintf(PETSC_COMM_WORLD,"Number of grid elements: %i\n",n_ele);
00095
00096 t_dnv=n3;
00097
00098 PetscPrintf(PETSC_COMM_WORLD,"Number of vertex data: %i\n",n3);
00099 PetscPrintf(PETSC_COMM_WORLD,"Number of element data: %i\n",n4);
00100 PetscPrintf(PETSC_COMM_WORLD,"Number of element data2: %i\n",n5);
00101
00102
00103
00104
00105
00106 if (gdata->vertxyz==PETSC_NULL) {
00107 flg=PETSC_FALSE;
00108 PetscPrintf(PETSC_COMM_WORLD,
00109 "Allocating array for vertices and reading their coordinates\n");
00110
00111 gdata->n_vert=n1;
00112 gdata->ln_vert=gdata->n_vert;
00113
00114 ierr = PetscMalloc(ND*gdata->n_vert*sizeof(PetscReal),&gdata->vertxyz);CHKERRQ(ierr);
00115 }
00116 else {
00117 flg=PETSC_TRUE;
00118 if (gdata->n_vert!=n1) {
00119 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00120 "Data in %s corrupt! Number of vertices does not agree with previously read mesh (%i != %i).\n",
00121 fmesh,gdata->n_vert,n1
00122 );
00123 }
00124 if (gdata->ln_vert!=gdata->n_vert) {
00125 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00126 "Data in %s corrupt! Number of elements does not agree with previously read mesh (%i != %i).\n",
00127 fmesh,gdata->ln_vert,gdata->n_vert
00128 );
00129 }
00130 }
00131
00132
00133 for (i=0; i<gdata->n_vert; i++) {
00134 #ifdef PETSC_USE_SINGLE
00135 fscanf(fd,"%i %f %f %f\n",
00136 #else
00137 fscanf(fd,"%i %lf %lf %lf\n",
00138 #endif
00139 &id,&d1,&d2,&d3
00140 );
00141
00142 id--;
00143
00144 if (flg==PETSC_FALSE) {
00145 gdata->vertxyz[ND*id+0]=d1;
00146 gdata->vertxyz[ND*id+1]=d2;
00147 gdata->vertxyz[ND*id+2]=d3;
00148 }
00149 else {
00150 if (PetscAbsReal(d1)>D_EPS) {
00151 if (PetscAbsReal((gdata->vertxyz[ND*id+0]-d1)/d1) > 1e-4)
00152 {
00153 PetscPrintf(PETSC_COMM_WORLD,"Data in %s corrupt!\n",fmesh);
00154 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00155 "x-coordinates of vertex %i does not agree with previously read mesh (%g != %g).\n",
00156 i,
00157 gdata->vertxyz[ND*id+0],
00158 d1
00159 );
00160 }
00161 }
00162 if (PetscAbsReal(d2)>D_EPS) {
00163 if (PetscAbsReal((gdata->vertxyz[ND*id+1]-d2)/d2) > 1e-4)
00164 {
00165 PetscPrintf(PETSC_COMM_WORLD,"Data in %s corrupt!\n",fmesh);
00166 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00167 "y-coordinates of vertex %i does not agree with previously read mesh (%g != %g).\n",
00168 i,
00169 gdata->vertxyz[ND*id+1],
00170 d2
00171 );
00172 }
00173 }
00174 if (PetscAbsReal(d3)>D_EPS) {
00175 if (PetscAbsReal((gdata->vertxyz[ND*id+2]-d3)/d3) > 1e-4)
00176 {
00177 PetscPrintf(PETSC_COMM_WORLD,"Data in %s corrupt!\n",fmesh);
00178 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00179 "z-coordinates of vertex %i does not agree with previously read mesh (%g != %g).\n",
00180 i,
00181 gdata->vertxyz[ND*id+2],
00182 d2
00183 );
00184 }
00185 }
00186 }
00187 }
00188
00189
00190
00191
00192
00193 if (gdata->elevert==PETSC_NULL) {
00194 flg=PETSC_FALSE;
00195 PetscPrintf(PETSC_COMM_WORLD,"Allocating array for element connectivity and reading it\n");
00196
00197 gdata->n_ele=n_ele;
00198 gdata->ln_ele=gdata->n_ele;
00199
00200 gdata->n_prop=0;
00201
00202 ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&gdata->elevert);CHKERRQ(ierr);
00203 ierr = PetscMalloc(gdata->n_ele*sizeof(int),&gdata->eleprop);CHKERRQ(ierr);
00204 }
00205 else {
00206 flg=PETSC_TRUE;
00207 assert(gdata->n_ele==n_ele || n_ele==0);
00208 assert(gdata->ln_ele==gdata->n_ele);
00209 }
00210
00211
00212 for (i=0; i<n_ele; i++) {
00213 fscanf(fd,"%i %i tet %i %i %i %i\n",
00214 &id,&a1,&n1,&n2,&n3,&n4
00215 );
00216
00217
00218 id--;
00219
00220 if (a1<1) {
00221 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00222 "Data in %s corrupt! Property of element %i is %i, but has to be >0! Exiting!\n",
00223 fmesh,i,a1
00224 );
00225 }
00226
00227 if (flg==PETSC_FALSE) {
00228 gdata->n_prop=PetscMax(gdata->n_prop,a1);
00229
00230 gdata->eleprop[id]=--a1;
00231
00232 gdata->elevert[id*NV+0]=--n1;
00233 gdata->elevert[id*NV+1]=--n2;
00234 gdata->elevert[id*NV+2]=--n3;
00235 gdata->elevert[id*NV+3]=--n4;
00236 }
00237 else {
00238
00239 assert(gdata->eleprop[id]==--a1);
00240
00241 assert(gdata->elevert[id*NV+0]==--n1);
00242 assert(gdata->elevert[id*NV+1]==--n2);
00243 assert(gdata->elevert[id*NV+2]==--n3);
00244 assert(gdata->elevert[id*NV+3]==--n4);
00245 }
00246 }
00247
00248 PetscPrintf(PETSC_COMM_WORLD,"Number of element prop.: %i\n",gdata->n_prop);
00249
00250 if (col<=0) {
00251 PetscFClose(PETSC_COMM_WORLD,fd);
00252 MagparFunctionLogReturn(0);
00253 }
00254
00255
00256
00257
00258
00259
00260 fscanf(fd,"%i\n",&a1);
00261 assert(a1==t_dnv);
00262
00263 if (a1<col+2) {
00264 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00265 "Data in %s corrupt! Cannot read data from column %i. Only %i data columns found! Exiting!\n",
00266 fmesh,col,a1
00267 );
00268 }
00269
00270
00271 fgets(msg,256,fd);
00272
00273
00274 for (i=0; i<t_dnv; i++) {
00275 fgets(msg,256,fd);
00276 if ((col-1)<=i && i<=(col+2-1)) {
00277 PetscPrintf(PETSC_COMM_WORLD,"reading: %s",msg);
00278 }
00279 }
00280
00281 for (i=0; i<gdata->n_vert; i++) {
00282
00283 fscanf(fd,"%i",&id);
00284
00285
00286 for (j=0; j<col-1; j++) {
00287 #ifdef PETSC_USE_SINGLE
00288 if (fscanf(fd,"%f",&d1)!=1) {
00289 #else
00290 if (fscanf(fd,"%lf",&d1)!=1) {
00291 #endif
00292 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00293 "Data in %s corrupt! Cannot read data from column %i. Only %i data columns found! Exiting!\n",
00294 fmesh,col,j
00295 );
00296 }
00297 }
00298
00299
00300 #ifdef PETSC_USE_SINGLE
00301 fscanf(fd,"%f %f %f",&d1,&d2,&d3);
00302 #else
00303 fscanf(fd,"%lf %lf %lf",&d1,&d2,&d3);
00304 #endif
00305
00306 fgets(msg,256,fd);
00307
00308
00309 id--;
00310
00311 ierr = VecSetValue(Vdest,ND*id+0,d1,INSERT_VALUES);CHKERRQ(ierr);
00312 ierr = VecSetValue(Vdest,ND*id+1,d2,INSERT_VALUES);CHKERRQ(ierr);
00313 ierr = VecSetValue(Vdest,ND*id+2,d3,INSERT_VALUES);CHKERRQ(ierr);
00314 }
00315
00316 ierr = VecAssemblyBegin(Vdest);CHKERRQ(ierr);
00317 ierr = VecAssemblyEnd(Vdest);CHKERRQ(ierr);
00318
00319 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00320
00321
00322 MagparFunctionLogReturn(0);
00323 }
00324