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 2681 2009-07-31 04:30:53Z 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 assert(a1>0);
00221
00222 if (flg==PETSC_FALSE) {
00223 gdata->n_prop=PetscMax(gdata->n_prop,a1);
00224
00225 gdata->eleprop[id]=--a1;
00226
00227 gdata->elevert[id*NV+0]=--n1;
00228 gdata->elevert[id*NV+1]=--n2;
00229 gdata->elevert[id*NV+2]=--n3;
00230 gdata->elevert[id*NV+3]=--n4;
00231 }
00232 else {
00233
00234 assert(gdata->eleprop[id]==--a1);
00235
00236 assert(gdata->elevert[id*NV+0]==--n1);
00237 assert(gdata->elevert[id*NV+1]==--n2);
00238 assert(gdata->elevert[id*NV+2]==--n3);
00239 assert(gdata->elevert[id*NV+3]==--n4);
00240 }
00241 }
00242
00243 PetscPrintf(PETSC_COMM_WORLD,"Number of element prop.: %i\n",gdata->n_prop);
00244
00245 if (col<=0) {
00246 PetscFClose(PETSC_COMM_WORLD,fd);
00247 MagparFunctionLogReturn(0);
00248 }
00249
00250
00251
00252
00253
00254
00255 fscanf(fd,"%i\n",&a1);
00256 assert(a1==t_dnv);
00257
00258 if (a1<col+2) {
00259 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00260 "Data in %s corrupt! Cannot read data from column %i. Only %i data columns found! Exiting!\n",
00261 fmesh,col,a1
00262 );
00263 }
00264
00265
00266 fgets(msg,256,fd);
00267
00268
00269 for (i=0; i<t_dnv; i++) {
00270 fgets(msg,256,fd);
00271 if ((col-1)<=i && i<=(col+2-1)) {
00272 PetscPrintf(PETSC_COMM_WORLD,"reading: %s",msg);
00273 }
00274 }
00275
00276 for (i=0; i<gdata->n_vert; i++) {
00277
00278 fscanf(fd,"%i",&id);
00279
00280
00281 for (j=0; j<col-1; j++) {
00282 #ifdef PETSC_USE_SINGLE
00283 if (fscanf(fd,"%f",&d1)!=1) {
00284 #else
00285 if (fscanf(fd,"%lf",&d1)!=1) {
00286 #endif
00287 SETERRQ3(PETSC_ERR_FILE_UNEXPECTED,
00288 "Data in %s corrupt! Cannot read data from column %i. Only %i data columns found! Exiting!\n",
00289 fmesh,col,j
00290 );
00291 }
00292 }
00293
00294
00295 #ifdef PETSC_USE_SINGLE
00296 fscanf(fd,"%f %f %f",&d1,&d2,&d3);
00297 #else
00298 fscanf(fd,"%lf %lf %lf",&d1,&d2,&d3);
00299 #endif
00300
00301 fgets(msg,256,fd);
00302
00303
00304 id--;
00305
00306 ierr = VecSetValue(Vdest,ND*id+0,d1,INSERT_VALUES);CHKERRQ(ierr);
00307 ierr = VecSetValue(Vdest,ND*id+1,d2,INSERT_VALUES);CHKERRQ(ierr);
00308 ierr = VecSetValue(Vdest,ND*id+2,d3,INSERT_VALUES);CHKERRQ(ierr);
00309 }
00310
00311 ierr = VecAssemblyBegin(Vdest);CHKERRQ(ierr);
00312 ierr = VecAssemblyEnd(Vdest);CHKERRQ(ierr);
00313
00314 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00315
00316
00317 MagparFunctionLogReturn(0);
00318 }
00319