readinp.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2009 Werner Scholz
00005 
00006     www:   http://www.magpar.net/
00007     email: magpar(at)magpar.net
00008 
00009     magpar is free software; you can redistribute it and/or modify
00010     it under the terms of the GNU General Public License as published by
00011     the Free Software Foundation; either version 2 of the License, or
00012     (at your option) any later version.
00013 
00014     magpar is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017     GNU General Public License for more details.
00018 
00019     You should have received a copy of the GNU General Public License
00020     along with magpar; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
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 /* Reads grid and/or data from AVS INP file */
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   /* check if we are supposed to read vertex data */
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   /* read
00082      number of vertices,
00083      number of elements,
00084      number of vertex data,
00085      number of element data,
00086      number of model data
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      Read in all vertices
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   /* read data of all vertices */
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      Read in all elements
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     /* our numbering starts with 0, not 1 */
00218     id--;
00219 
00220     assert(a1>0);
00221 
00222     if (flg==PETSC_FALSE) {
00223       gdata->n_prop=PetscMax(gdata->n_prop,a1);
00224       /* internal property ids start with 0 !!! */
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       /* internal property ids start with 0 !!! */
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      Read vertex data
00252   */
00253 
00254   /* read number of components of vertex data */
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   /* read remainder of this line */
00266   fgets(msg,256,fd);
00267 
00268   /* we assume that all components have size 1 ! */
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     /* read id */
00278     fscanf(fd,"%i",&id);
00279 
00280     /* skip useless columns */
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     /* read data */
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     /* read remainder of this line */
00301     fgets(msg,256,fd);
00302 
00303     /* our numbering starts with 0, not 1 */
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 

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2009 Werner Scholz