1 2 #include <petsc.h> 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 #include <petscdmswarm.h> 6 #include <petsc/private/dmswarmimpl.h> 7 #include "data_bucket.h" 8 9 PetscErrorCode private_PetscViewerCreate_XDMF(MPI_Comm comm,const char filename[],PetscViewer *v) 10 { 11 long int *bytes; 12 PetscContainer container; 13 PetscViewer viewer; 14 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 ierr = PetscViewerCreate(comm,&viewer);CHKERRQ(ierr); 18 ierr = PetscViewerSetType(viewer,PETSCVIEWERASCII);CHKERRQ(ierr); 19 ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 20 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 21 22 ierr = PetscMalloc1(1,&bytes);CHKERRQ(ierr); 23 bytes[0] = 0; 24 ierr = PetscContainerCreate(comm,&container);CHKERRQ(ierr); 25 ierr = PetscContainerSetPointer(container,(void*)bytes);CHKERRQ(ierr); 26 ierr = PetscObjectCompose((PetscObject)viewer,"XDMFViewerContext",(PetscObject)container);CHKERRQ(ierr); 27 28 /* write xdmf header */ 29 PetscViewerASCIIPrintf(viewer,"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n"); 30 PetscViewerASCIIPrintf(viewer,"<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\" Version=\"2.99\">\n"); 31 /* write xdmf domain */ 32 PetscViewerASCIIPushTab(viewer); 33 PetscViewerASCIIPrintf(viewer,"<Domain>\n"); 34 *v = viewer; 35 PetscFunctionReturn(0); 36 } 37 38 PetscErrorCode private_PetscViewerDestroy_XDMF(PetscViewer *v) 39 { 40 PetscViewer viewer; 41 DM dm = NULL; 42 long int *bytes; 43 PetscContainer container = NULL; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 if (!v) PetscFunctionReturn(0); 48 viewer = *v; 49 50 ierr = PetscObjectQuery((PetscObject)viewer,"DMSwarm",(PetscObject*)&dm);CHKERRQ(ierr); 51 if (dm) { 52 PetscViewerASCIIPrintf(viewer,"</Grid>\n"); 53 PetscViewerASCIIPopTab(viewer); 54 } 55 56 /* close xdmf header */ 57 PetscViewerASCIIPrintf(viewer,"</Domain>\n"); 58 PetscViewerASCIIPopTab(viewer); 59 PetscViewerASCIIPrintf(viewer,"</Xdmf>\n"); 60 61 ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr); 62 if (container) { 63 ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr); 64 ierr = PetscFree(bytes);CHKERRQ(ierr); 65 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 66 } 67 68 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 69 *v = NULL; 70 PetscFunctionReturn(0); 71 } 72 73 PetscErrorCode private_CreateDataFileNameXDMF(const char filename[],char dfilename[]) 74 { 75 char *ext; 76 PetscBool flg; 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 ierr = PetscStrrchr(filename,'.',&ext);CHKERRQ(ierr); 81 ierr = PetscStrcmp("xmf",ext,&flg);CHKERRQ(ierr); 82 if (flg) { 83 size_t len; 84 char viewername_minus_ext[PETSC_MAX_PATH_LEN]; 85 86 ierr = PetscStrlen(filename,&len);CHKERRQ(ierr); 87 ierr = PetscSNPrintf(viewername_minus_ext,len-3,"%s",filename);CHKERRQ(ierr); 88 ierr = PetscSNPrintf(dfilename,PETSC_MAX_PATH_LEN-1,"%s_swarm_fields.pbin",viewername_minus_ext);CHKERRQ(ierr); 89 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"File extension must by .xmf"); 90 91 PetscFunctionReturn(0); 92 } 93 94 PetscErrorCode private_DMSwarmView_XDMF(DM dm,PetscViewer viewer) 95 { 96 PetscBool isswarm = PETSC_FALSE; 97 const char *viewername; 98 char datafile[PETSC_MAX_PATH_LEN]; 99 PetscViewer fviewer; 100 PetscInt k,ng,dim; 101 Vec dvec; 102 long int *bytes = NULL; 103 PetscContainer container = NULL; 104 const char *dmname; 105 PetscErrorCode ierr; 106 107 PetscFunctionBegin; 108 ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr); 109 if (container) { 110 ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr); 111 } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext"); 112 113 ierr = PetscObjectTypeCompare((PetscObject)dm,DMSWARM,&isswarm);CHKERRQ(ierr); 114 if (!isswarm) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Only valid for DMSwarm"); 115 116 ierr = PetscObjectCompose((PetscObject)viewer,"DMSwarm",(PetscObject)dm);CHKERRQ(ierr); 117 118 PetscViewerASCIIPushTab(viewer); 119 ierr = PetscObjectGetName((PetscObject)dm,&dmname);CHKERRQ(ierr); 120 if (!dmname) { 121 ierr = DMGetOptionsPrefix(dm,&dmname);CHKERRQ(ierr); 122 } 123 if (!dmname) { 124 PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n"); 125 } else { 126 PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n",dmname); 127 } 128 129 /* create a sub-viewer for topology, geometry and all data fields */ 130 /* name is viewer.name + "_swarm_fields.pbin" */ 131 ierr = PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);CHKERRQ(ierr); 132 ierr = PetscViewerSetType(fviewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 133 ierr = PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);CHKERRQ(ierr); 134 ierr = PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);CHKERRQ(ierr); 135 ierr = PetscViewerFileSetMode(fviewer,FILE_MODE_WRITE);CHKERRQ(ierr); 136 137 ierr = PetscViewerFileGetName(viewer,&viewername);CHKERRQ(ierr); 138 ierr = private_CreateDataFileNameXDMF(viewername,datafile);CHKERRQ(ierr); 139 ierr = PetscViewerFileSetName(fviewer,datafile);CHKERRQ(ierr); 140 141 ierr = DMSwarmGetSize(dm,&ng);CHKERRQ(ierr); 142 143 /* write topology header */ 144 PetscViewerASCIIPushTab(viewer); 145 PetscViewerASCIIPrintf(viewer,"<Topology Dimensions=\"%D\" TopologyType=\"Mixed\">\n",ng); 146 PetscViewerASCIIPushTab(viewer); 147 PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%D\" Seek=\"%D\">\n",ng*3,bytes[0]); 148 PetscViewerASCIIPushTab(viewer); 149 PetscViewerASCIIPrintf(viewer,"%s\n",datafile); 150 PetscViewerASCIIPopTab(viewer); 151 PetscViewerASCIIPrintf(viewer,"</DataItem>\n"); 152 PetscViewerASCIIPopTab(viewer); 153 PetscViewerASCIIPrintf(viewer,"</Topology>\n"); 154 PetscViewerASCIIPopTab(viewer); 155 156 /* write topology data */ 157 for (k=0; k<ng; k++) { 158 PetscInt pvertex[3]; 159 160 pvertex[0] = 1; 161 pvertex[1] = 1; 162 pvertex[2] = k; 163 ierr = PetscViewerBinaryWrite(fviewer,pvertex,3,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 164 } 165 bytes[0] += sizeof(PetscInt) * ng * 3; 166 167 /* write geometry header */ 168 PetscViewerASCIIPushTab(viewer); 169 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 170 switch (dim) { 171 case 1: 172 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for 1D"); 173 break; 174 case 2: 175 PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XY\">\n"); 176 break; 177 case 3: 178 PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XYZ\">\n"); 179 break; 180 } 181 PetscViewerASCIIPushTab(viewer); 182 PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",ng,dim,bytes[0]); 183 PetscViewerASCIIPushTab(viewer); 184 PetscViewerASCIIPrintf(viewer,"%s\n",datafile); 185 PetscViewerASCIIPopTab(viewer); 186 PetscViewerASCIIPrintf(viewer,"</DataItem>\n"); 187 PetscViewerASCIIPopTab(viewer); 188 PetscViewerASCIIPrintf(viewer,"</Geometry>\n"); 189 PetscViewerASCIIPopTab(viewer); 190 191 /* write geometry data */ 192 ierr = DMSwarmCreateGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);CHKERRQ(ierr); 193 ierr = VecView(dvec,fviewer);CHKERRQ(ierr); 194 ierr = DMSwarmDestroyGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);CHKERRQ(ierr); 195 bytes[0] += sizeof(PetscReal) * ng * dim; 196 197 ierr = PetscViewerDestroy(&fviewer);CHKERRQ(ierr); 198 199 PetscFunctionReturn(0); 200 } 201 202 PetscErrorCode private_VecView_Swarm_XDMF(Vec x,PetscViewer viewer) 203 { 204 long int *bytes = NULL; 205 PetscContainer container = NULL; 206 const char *viewername; 207 char datafile[PETSC_MAX_PATH_LEN]; 208 PetscViewer fviewer; 209 PetscInt N,bs; 210 const char *vecname; 211 char fieldname[PETSC_MAX_PATH_LEN]; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr); 216 if (container) { 217 ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr); 218 } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext"); 219 220 ierr = PetscViewerFileGetName(viewer,&viewername);CHKERRQ(ierr); 221 ierr = private_CreateDataFileNameXDMF(viewername,datafile);CHKERRQ(ierr); 222 223 /* re-open a sub-viewer for all data fields */ 224 /* name is viewer.name + "_swarm_fields.pbin" */ 225 ierr = PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);CHKERRQ(ierr); 226 ierr = PetscViewerSetType(fviewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 227 ierr = PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);CHKERRQ(ierr); 228 ierr = PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);CHKERRQ(ierr); 229 ierr = PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);CHKERRQ(ierr); 230 ierr = PetscViewerFileSetName(fviewer,datafile);CHKERRQ(ierr); 231 232 ierr = VecGetSize(x,&N);CHKERRQ(ierr); 233 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 234 N = N/bs; 235 ierr = PetscObjectGetName((PetscObject)x,&vecname);CHKERRQ(ierr); 236 if (!vecname) { 237 ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)x)->tag);CHKERRQ(ierr); 238 } else { 239 ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);CHKERRQ(ierr); 240 } 241 242 /* write data header */ 243 PetscViewerASCIIPushTab(viewer); 244 PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname); 245 PetscViewerASCIIPushTab(viewer); 246 if (bs == 1) { 247 PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]); 248 } else { 249 PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]); 250 } 251 PetscViewerASCIIPushTab(viewer); 252 PetscViewerASCIIPrintf(viewer,"%s\n",datafile); 253 PetscViewerASCIIPopTab(viewer); 254 PetscViewerASCIIPrintf(viewer,"</DataItem>\n"); 255 PetscViewerASCIIPopTab(viewer); 256 PetscViewerASCIIPrintf(viewer,"</Attribute>\n"); 257 PetscViewerASCIIPopTab(viewer); 258 259 /* write data */ 260 ierr = VecView(x,fviewer);CHKERRQ(ierr); 261 bytes[0] += sizeof(PetscReal) * N * bs; 262 263 ierr = PetscViewerDestroy(&fviewer);CHKERRQ(ierr); 264 265 PetscFunctionReturn(0); 266 } 267 268 /*@C 269 DMSwarmViewFieldsXDMF - Write a selection of DMSwarm fields to an XDMF3 file 270 271 Collective on DM 272 273 Input parameters: 274 + dm - the DMSwarm 275 . filename - the file name of the XDMF file (must have the extension .xmf) 276 . nfields - the number of fields to write into the XDMF file 277 - field_name_list - array of length nfields containing the textual name of fields to write 278 279 Level: beginner 280 281 Notes: 282 Only fields registered with data type PETSC_REAL can be written into the file 283 284 .seealso: DMSwarmViewXDMF() 285 @*/ 286 PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char *field_name_list[]) 287 { 288 PetscErrorCode ierr; 289 Vec dvec; 290 PetscInt f; 291 PetscViewer viewer; 292 293 PetscFunctionBegin; 294 ierr = private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);CHKERRQ(ierr); 295 ierr = private_DMSwarmView_XDMF(dm,viewer);CHKERRQ(ierr); 296 for (f=0; f<nfields; f++) { 297 ierr = DMSwarmCreateGlobalVectorFromField(dm,field_name_list[f],&dvec);CHKERRQ(ierr); 298 ierr = PetscObjectSetName((PetscObject)dvec,field_name_list[f]);CHKERRQ(ierr); 299 ierr = private_VecView_Swarm_XDMF(dvec,viewer);CHKERRQ(ierr); 300 ierr = DMSwarmDestroyGlobalVectorFromField(dm,field_name_list[f],&dvec);CHKERRQ(ierr); 301 } 302 ierr = private_PetscViewerDestroy_XDMF(&viewer);CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 /*@C 307 DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file 308 309 Collective on DM 310 311 Input parameters: 312 + dm - the DMSwarm 313 - filename - the file name of the XDMF file (must have the extension .xmf) 314 315 Level: beginner 316 317 Notes: 318 Only fields registered with data type PETSC_REAL will be written into the file 319 320 .seealso: DMSwarmViewFieldsXDMF() 321 @*/ 322 PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm,const char filename[]) 323 { 324 DM_Swarm *swarm = (DM_Swarm*)dm->data; 325 PetscErrorCode ierr; 326 Vec dvec; 327 PetscInt f; 328 PetscViewer viewer; 329 330 PetscFunctionBegin; 331 ierr = private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);CHKERRQ(ierr); 332 ierr = private_DMSwarmView_XDMF(dm,viewer);CHKERRQ(ierr); 333 for (f=4; f<swarm->db->nfields; f++) { /* only examine user defined fields - the first 4 are internally created by DMSwarmPIC */ 334 DataField field; 335 336 /* query field type - accept all those of type PETSC_DOUBLE */ 337 field = swarm->db->field[f]; 338 if (field->petsc_type != PETSC_DOUBLE) continue; 339 340 ierr = DMSwarmCreateGlobalVectorFromField(dm,field->name,&dvec);CHKERRQ(ierr); 341 ierr = PetscObjectSetName((PetscObject)dvec,field->name);CHKERRQ(ierr); 342 ierr = private_VecView_Swarm_XDMF(dvec,viewer);CHKERRQ(ierr); 343 ierr = DMSwarmDestroyGlobalVectorFromField(dm,field->name,&dvec);CHKERRQ(ierr); 344 } 345 ierr = private_PetscViewerDestroy_XDMF(&viewer);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348