1 #include <petsc/private/glvisviewerimpl.h> 2 #include <petsc/private/petscimpl.h> 3 #include <petsc/private/dmpleximpl.h> 4 #include <petscbt.h> 5 #include <petscdmplex.h> 6 #include <petscsf.h> 7 #include <petscds.h> 8 9 typedef struct { 10 PetscInt nf; 11 VecScatter *scctx; 12 } GLVisViewerCtx; 13 14 static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx) 15 { 16 GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx; 17 PetscInt i; 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 for (i=0;i<ctx->nf;i++) { 22 ierr = VecScatterDestroy(&ctx->scctx[i]);CHKERRQ(ierr); 23 } 24 ierr = PetscFree(ctx->scctx);CHKERRQ(ierr); 25 ierr = PetscFree(vctx);CHKERRQ(ierr); 26 PetscFunctionReturn(0); 27 } 28 29 static PetscErrorCode DMPlexSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx) 30 { 31 GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx; 32 PetscInt f; 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 for (f=0;f<nf;f++) { 37 ierr = VecScatterBegin(ctx->scctx[f],(Vec)oX,(Vec)oXfield[f],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 38 ierr = VecScatterEnd(ctx->scctx[f],(Vec)oX,(Vec)oXfield[f],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 39 } 40 PetscFunctionReturn(0); 41 } 42 43 /* for FEM, it works for H1 fields only and extracts dofs at cell vertices, discarding any other dof */ 44 PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject odm, PetscViewer viewer) 45 { 46 DM dm = (DM)odm; 47 Vec xlocal,xfield; 48 PetscDS ds; 49 IS globalNum,isfield; 50 PetscBT vown; 51 char **fieldname = NULL,**fec_type = NULL; 52 const PetscInt *gNum; 53 PetscInt *nlocal,*bs,*idxs,*dims; 54 PetscInt f,maxfields,nfields,c,totc,totdofs,Nv,cum,i; 55 PetscInt dim,cStart,cEnd,cEndInterior,vStart,vEnd; 56 GLVisViewerCtx *ctx; 57 PetscSection s; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 62 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 63 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 64 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 65 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 66 ierr = DMPlexGetCellNumbering(dm,&globalNum);CHKERRQ(ierr); 67 ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr); 68 ierr = PetscBTCreate(vEnd-vStart,&vown);CHKERRQ(ierr); 69 for (c = cStart, totc = 0; c < cEnd; c++) { 70 if (gNum[c-cStart] >= 0) { 71 PetscInt i,numPoints,*points = NULL; 72 73 totc++; 74 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 75 for (i=0;i<numPoints*2;i+= 2) { 76 if ((points[i] >= vStart) && (points[i] < vEnd)) { 77 ierr = PetscBTSet(vown,points[i]-vStart);CHKERRQ(ierr); 78 } 79 } 80 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 81 } 82 } 83 for (f=0,Nv=0;f<vEnd-vStart;f++) if (PetscLikely(PetscBTLookup(vown,f))) Nv++; 84 85 ierr = DMCreateLocalVector(dm,&xlocal);CHKERRQ(ierr); 86 ierr = VecGetLocalSize(xlocal,&totdofs);CHKERRQ(ierr); 87 ierr = DMGetDefaultSection(dm,&s);CHKERRQ(ierr); 88 ierr = PetscSectionGetNumFields(s,&nfields);CHKERRQ(ierr); 89 for (f=0,maxfields=0;f<nfields;f++) { 90 PetscInt bs; 91 92 ierr = PetscSectionGetFieldComponents(s,f,&bs);CHKERRQ(ierr); 93 maxfields += bs; 94 } 95 ierr = PetscCalloc6(maxfields,&fieldname,maxfields,&nlocal,maxfields,&bs,maxfields,&dims,maxfields,&fec_type,totdofs,&idxs);CHKERRQ(ierr); 96 ierr = PetscNew(&ctx);CHKERRQ(ierr); 97 ierr = PetscCalloc1(maxfields,&ctx->scctx);CHKERRQ(ierr); 98 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 99 if (ds) { 100 for (f=0;f<nfields;f++) { 101 const char* fname; 102 char name[256]; 103 PetscObject disc; 104 size_t len; 105 106 ierr = PetscSectionGetFieldName(s,f,&fname);CHKERRQ(ierr); 107 ierr = PetscStrlen(fname,&len);CHKERRQ(ierr); 108 if (len) { 109 ierr = PetscStrcpy(name,fname);CHKERRQ(ierr); 110 } else { 111 ierr = PetscSNPrintf(name,256,"Field%d",f);CHKERRQ(ierr); 112 } 113 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 114 if (disc) { 115 PetscClassId id; 116 PetscInt Nc; 117 char fec[64]; 118 119 ierr = PetscObjectGetClassId(disc, &id);CHKERRQ(ierr); 120 if (id == PETSCFE_CLASSID) { 121 PetscFE fem = (PetscFE)disc; 122 PetscDualSpace sp; 123 PetscDualSpaceType spname; 124 PetscInt order; 125 PetscBool islag,continuous,H1 = PETSC_TRUE; 126 127 ierr = PetscFEGetNumComponents(fem,&Nc);CHKERRQ(ierr); 128 ierr = PetscFEGetDualSpace(fem,&sp);CHKERRQ(ierr); 129 ierr = PetscDualSpaceGetType(sp,&spname);CHKERRQ(ierr); 130 ierr = PetscStrcmp(spname,PETSCDUALSPACELAGRANGE,&islag);CHKERRQ(ierr); 131 if (!islag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dual space"); 132 ierr = PetscDualSpaceLagrangeGetContinuity(sp,&continuous);CHKERRQ(ierr); 133 if (!continuous) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Discontinuous space visualization currently unsupported"); 134 ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr); 135 if (continuous && order > 0) { 136 ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: H1_%dD_P1",dim);CHKERRQ(ierr); 137 } else { 138 H1 = PETSC_FALSE; 139 ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%dD_P%d",dim,order);CHKERRQ(ierr); 140 } 141 ierr = PetscStrallocpy(name,&fieldname[ctx->nf]);CHKERRQ(ierr); 142 bs[ctx->nf] = Nc; 143 dims[ctx->nf] = dim; 144 if (H1) { 145 nlocal[ctx->nf] = Nc * Nv; 146 ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr); 147 ierr = VecCreateSeq(PETSC_COMM_SELF,Nv*Nc,&xfield);CHKERRQ(ierr); 148 for (i=0,cum=0;i<vEnd-vStart;i++) { 149 PetscInt j,off; 150 151 if (PetscUnlikely(!PetscBTLookup(vown,i))) continue; 152 ierr = PetscSectionGetFieldOffset(s,i+vStart,f,&off);CHKERRQ(ierr); 153 for (j=0;j<Nc;j++) idxs[cum++] = off + j; 154 } 155 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),Nv*Nc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr); 156 } else { 157 nlocal[ctx->nf] = Nc * totc; 158 ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr); 159 ierr = VecCreateSeq(PETSC_COMM_SELF,Nc*totc,&xfield);CHKERRQ(ierr); 160 for (i=0,cum=0;i<cEnd-cStart;i++) { 161 PetscInt j,off; 162 163 if (PetscUnlikely(gNum[i] < 0)) continue; 164 ierr = PetscSectionGetFieldOffset(s,i+cStart,f,&off);CHKERRQ(ierr); 165 for (j=0;j<Nc;j++) idxs[cum++] = off + j; 166 } 167 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc*Nc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr); 168 } 169 ierr = VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);CHKERRQ(ierr); 170 ierr = VecDestroy(&xfield);CHKERRQ(ierr); 171 ierr = ISDestroy(&isfield);CHKERRQ(ierr); 172 ctx->nf++; 173 } else if (id == PETSCFV_CLASSID) { 174 PetscInt c; 175 176 ierr = PetscFVGetNumComponents((PetscFV)disc,&Nc);CHKERRQ(ierr); 177 ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%dD_P0",dim);CHKERRQ(ierr); 178 for (c = 0; c < Nc; c++) { 179 char comp[256]; 180 ierr = PetscSNPrintf(comp,256,"%s-Comp%d",name,c);CHKERRQ(ierr); 181 ierr = PetscStrallocpy(comp,&fieldname[ctx->nf]);CHKERRQ(ierr); 182 bs[ctx->nf] = 1; /* Does PetscFV support components with different block size? */ 183 nlocal[ctx->nf] = totc; 184 dims[ctx->nf] = dim; 185 ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr); 186 ierr = VecCreateSeq(PETSC_COMM_SELF,totc,&xfield);CHKERRQ(ierr); 187 for (i=0,cum=0;i<cEnd-cStart;i++) { 188 PetscInt off; 189 190 if (PetscUnlikely(gNum[i])<0) continue; 191 ierr = PetscSectionGetFieldOffset(s,i+cStart,f,&off);CHKERRQ(ierr); 192 idxs[cum++] = off + c; 193 } 194 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr); 195 ierr = VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);CHKERRQ(ierr); 196 ierr = VecDestroy(&xfield);CHKERRQ(ierr); 197 ierr = ISDestroy(&isfield);CHKERRQ(ierr); 198 ctx->nf++; 199 } 200 } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",f); 201 } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing discretization for field %D",f); 202 } 203 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Needs a DS attached to the DM"); 204 ierr = PetscBTDestroy(&vown);CHKERRQ(ierr); 205 ierr = VecDestroy(&xlocal);CHKERRQ(ierr); 206 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 207 208 /* customize the viewer */ 209 ierr = PetscViewerGLVisSetFields(viewer,ctx->nf,(const char**)fieldname,(const char**)fec_type,nlocal,bs,dims,DMPlexSampleGLVisFields_Private,ctx,DestroyGLVisViewerCtx_Private);CHKERRQ(ierr); 210 for (f=0;f<ctx->nf;f++) { 211 ierr = PetscFree(fieldname[f]);CHKERRQ(ierr); 212 ierr = PetscFree(fec_type[f]);CHKERRQ(ierr); 213 } 214 ierr = PetscFree6(fieldname,nlocal,bs,dims,fec_type,idxs);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 typedef enum {MFEM_POINT=0,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE,MFEM_TETRAHEDRON,MFEM_CUBE,MFEM_UNDEF} MFEM_cid; 219 220 MFEM_cid mfem_table_cid[4][7] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF, MFEM_UNDEF}, 221 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF, MFEM_UNDEF}, 222 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE ,MFEM_UNDEF, MFEM_UNDEF}, 223 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_TETRAHEDRON,MFEM_UNDEF, MFEM_CUBE } }; 224 225 static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt p, PetscInt *mid, PetscInt *cid) 226 { 227 DMLabel dlabel; 228 PetscInt depth,csize; 229 PetscErrorCode ierr; 230 231 PetscFunctionBegin; 232 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 233 ierr = DMLabelGetValue(dlabel,p,&depth);CHKERRQ(ierr); 234 ierr = DMPlexGetConeSize(dm,p,&csize);CHKERRQ(ierr); 235 if (label) { 236 ierr = DMLabelGetValue(label,p,mid);CHKERRQ(ierr); 237 } else *mid = 1; 238 *cid = mfem_table_cid[depth][csize]; 239 PetscFunctionReturn(0); 240 } 241 242 static PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscSection csec, PetscInt *nv, int vids[]) 243 { 244 PetscInt dim,sdim,dof = 0,off = 0,i,q,vStart,vEnd,numPoints,*points = NULL; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 249 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 250 sdim = dim; 251 if (csec) { 252 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 253 ierr = PetscSectionGetOffset(csec,vStart,&off);CHKERRQ(ierr); 254 off = off/sdim; 255 ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr); 256 } 257 if (!dof) { 258 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 259 for (i=0,q=0;i<numPoints*2;i+= 2) 260 if ((points[i] >= vStart) && (points[i] < vEnd)) 261 vids[q++] = points[i]-vStart+off; 262 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 263 } else { 264 ierr = PetscSectionGetOffset(csec,p,&off);CHKERRQ(ierr); 265 ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr); 266 for (q=0;q<dof/sdim;q++) vids[q] = off/sdim + q; 267 } 268 *nv = q; 269 PetscFunctionReturn(0); 270 } 271 272 /* 273 ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR 274 Higher order meshes are also supported 275 */ 276 static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer) 277 { 278 DMLabel label; 279 PetscSection coordSection,parentSection; 280 Vec coordinates,hovec; 281 IS globalNum = NULL; 282 const PetscScalar *array; 283 const PetscInt *gNum; 284 PetscInt bf,p,sdim,dim,depth,novl; 285 PetscInt cStart,cEnd,cEndInterior,vStart,vEnd,nvert; 286 PetscMPIInt commsize; 287 PetscBool localized,ovl,isascii,enable_boundary,enable_ncmesh; 288 PetscBT pown; 289 PetscErrorCode ierr; 290 PetscContainer glvis_container; 291 PetscBool periodic, enabled = PETSC_TRUE; 292 const char *fmt; 293 294 PetscFunctionBegin; 295 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 296 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 297 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 298 if (!isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERASCII"); 299 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&commsize);CHKERRQ(ierr); 300 if (commsize > 1) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Use single sequential viewers for parallel visualization"); 301 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 302 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 303 if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 304 305 /* get container: determines if a process visualizes is portion of the data or not */ 306 ierr = PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);CHKERRQ(ierr); 307 if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis container"); 308 { 309 PetscViewerGLVisInfo glvis_info; 310 ierr = PetscContainerGetPointer(glvis_container,(void**)&glvis_info);CHKERRQ(ierr); 311 enabled = glvis_info->enabled; 312 fmt = glvis_info->fmt; 313 } 314 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 315 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 316 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 317 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 318 ierr = DMGetPeriodicity(dm,&periodic,NULL,NULL,NULL);CHKERRQ(ierr); 319 ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr); 320 if (periodic && !localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Coordinates need to be localized"); 321 ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 322 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 323 ierr = DMGetCoordinatesLocal(dm,&coordinates);CHKERRQ(ierr); 324 325 /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh 326 DMPlex does not currently support HO meshes, so there's no API for this */ 327 ierr = PetscObjectQuery((PetscObject)dm,"_glvis_mesh_coords",(PetscObject*)&hovec);CHKERRQ(ierr); 328 329 /* 330 a couple of sections of the mesh specification are disabled 331 - boundary: unless we want to visualize boundary attributes or we have a periodic mesh 332 the boundary is not needed for proper mesh visualization 333 - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package 334 and be able to derefine the mesh 335 */ 336 enable_boundary = periodic; 337 enable_ncmesh = PETSC_FALSE; 338 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"GLVis PetscViewer DMPlex Options","PetscViewer");CHKERRQ(ierr); 339 ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary","Enable boundary section in mesh representation; useful for debugging purposes",NULL,enable_boundary,&enable_boundary,NULL);CHKERRQ(ierr); 340 ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_ncmesh","Enable vertex_parents section in mesh representation; useful for debugging purposes",NULL,enable_ncmesh,&enable_ncmesh,NULL);CHKERRQ(ierr); 341 ierr = PetscOptionsEnd();CHKERRQ(ierr); 342 343 /* Identify possible cells in the overlap */ 344 gNum = NULL; 345 novl = 0; 346 ovl = PETSC_FALSE; 347 pown = NULL; 348 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr); 349 if (commsize > 1) { 350 ierr = DMPlexGetCellNumbering(dm,&globalNum);CHKERRQ(ierr); 351 ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr); 352 for (p=cStart; p<cEnd; p++) { 353 if (gNum[p-cStart] < 0) { 354 ovl = PETSC_TRUE; 355 novl++; 356 } 357 } 358 if (ovl) { 359 /* it may happen that pown get not destroyed, if the user closes the window while this function is running. 360 TODO: garbage collector? attach pown to dm? */ 361 ierr = PetscBTCreate(PetscMax(cEnd-cStart,vEnd-vStart),&pown);CHKERRQ(ierr); 362 } 363 } 364 365 if (!enabled) { 366 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 367 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 368 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 369 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 370 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 371 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 372 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 373 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 374 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 375 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 376 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 377 if (globalNum) { 378 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 379 } 380 PetscFunctionReturn(0); 381 } 382 383 /* header */ 384 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 385 386 /* topological dimension */ 387 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 388 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 389 390 /* elements */ 391 /* TODO: is this the label we want to use for marking material IDs? 392 We should decide to have a single marker for these stuff 393 Proposal: DMSetMaterialLabel? 394 */ 395 ierr = DMGetLabel(dm,"Cell Sets",&label);CHKERRQ(ierr); 396 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 397 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);CHKERRQ(ierr); 398 for (p=cStart;p<cEnd;p++) { 399 int vids[8]; 400 PetscInt i,nv = 0,cid = -1,mid = 1; 401 402 if (ovl) { 403 if (gNum[p-cStart] < 0) continue; 404 else { 405 ierr = PetscBTSet(pown,p-cStart);CHKERRQ(ierr); 406 } 407 } 408 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);CHKERRQ(ierr); 409 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,p,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 410 ierr = DMPlexInvertCell(dim,nv,vids);CHKERRQ(ierr); 411 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 412 for (i=0;i<nv;i++) { 413 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[i]);CHKERRQ(ierr); 414 } 415 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 416 } 417 418 /* boundary */ 419 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 420 if (!enable_boundary) { 421 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 422 } else { 423 DMLabel perLabel; 424 PetscBT bfaces; 425 PetscInt fStart,fEnd,fEndInterior,*fcells; 426 PetscInt *faces = NULL,fpc = 0,vpf = 0, vpc = 0; 427 PetscInt fv1[] = {0,1}, 428 fv2tri[] = {0,1, 429 1,2, 430 2,0}, 431 fv2quad[] = {0,1, 432 1,2, 433 2,3, 434 3,0}, 435 fv3tet[] = {0,1,2, 436 0,3,1, 437 0,2,3, 438 2,1,3}, 439 fv3hex[] = {0,1,2,3, 440 4,5,6,7, 441 0,3,5,4, 442 2,1,7,6, 443 3,2,6,5, 444 0,4,7,1}; 445 446 /* determine orientation of boundary mesh */ 447 if (cEnd-cStart) { 448 ierr = DMPlexGetConeSize(dm,cStart,&fpc);CHKERRQ(ierr); 449 switch(dim) { 450 case 1: 451 if (fpc != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case faces per cell %D",fpc); 452 faces = fv1; 453 vpf = 1; 454 vpc = 2; 455 break; 456 case 2: 457 switch (fpc) { 458 case 3: 459 faces = fv2tri; 460 vpf = 2; 461 vpc = 3; 462 break; 463 case 4: 464 faces = fv2quad; 465 vpf = 2; 466 vpc = 4; 467 break; 468 default: 469 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 470 break; 471 } 472 break; 473 case 3: 474 switch (fpc) { 475 case 4: 476 faces = fv3tet; 477 vpf = 3; 478 vpc = 4; 479 break; 480 case 6: 481 faces = fv3hex; 482 vpf = 4; 483 vpc = 8; 484 break; 485 default: 486 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 487 break; 488 } 489 break; 490 default: 491 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim"); 492 break; 493 } 494 } 495 ierr = DMPlexGetHybridBounds(dm,NULL,&fEndInterior,NULL,NULL);CHKERRQ(ierr); 496 ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);CHKERRQ(ierr); 497 fEnd = fEndInterior < 0 ? fEnd : fEndInterior; 498 ierr = PetscBTCreate(fEnd-fStart,&bfaces);CHKERRQ(ierr); 499 ierr = DMPlexGetMaxSizes(dm,NULL,&p);CHKERRQ(ierr); 500 ierr = PetscMalloc1(p,&fcells);CHKERRQ(ierr); 501 ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr); 502 if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */ 503 ierr = DMCreateLabel(dm,"glvis_periodic_cut");CHKERRQ(ierr); 504 ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr); 505 ierr = DMLabelSetDefaultValue(perLabel,1);CHKERRQ(ierr); 506 for (p=cStart;p<cEnd;p++) { 507 PetscInt dof; 508 ierr = PetscSectionGetDof(coordSection,p,&dof);CHKERRQ(ierr); 509 if (dof) { 510 PetscInt v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL; 511 PetscScalar *vals = NULL; 512 513 if (dof%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Incompatible number of cell dofs %D and space dimension %D",dof,sdim); 514 if (dof/sdim != vpc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Incompatible number of cell dofs %D and vertices %D",dof/sdim,vpc); 515 ierr = DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 516 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 517 for (v=0;v<cellClosureSize;v++) 518 if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) { 519 vidxs = cellClosure + 2*v; 520 break; 521 } 522 if (!vidxs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing vertices"); 523 for (v=0;v<vpc;v++) { 524 PetscInt s; 525 for (s=0;s<sdim;s++) { 526 if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+vpc*sdim])>PETSC_MACHINE_EPSILON) { 527 ierr = DMLabelSetValue(perLabel,vidxs[2*v],2);CHKERRQ(ierr); 528 } 529 } 530 } 531 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 532 ierr = DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 533 } 534 } 535 if (dim > 1) { 536 PetscInt eEnd,eStart,eEndInterior; 537 ierr = DMPlexGetHybridBounds(dm,NULL,NULL,&eEndInterior,NULL);CHKERRQ(ierr); 538 ierr = DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);CHKERRQ(ierr); 539 eEnd = eEndInterior < 0 ? eEnd : eEndInterior; 540 for (p=eStart;p<eEnd;p++) { 541 const PetscInt *cone; 542 PetscInt coneSize,i; 543 PetscBool ispe = PETSC_TRUE; 544 545 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 546 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 547 for (i=0;i<coneSize;i++) { 548 PetscInt v; 549 550 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 551 ispe = (PetscBool)(ispe && (v==2)); 552 } 553 if (ispe && coneSize) { 554 ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr); 555 } 556 } 557 if (dim > 2) { 558 for (p=fStart;p<fEnd;p++) { 559 const PetscInt *cone; 560 PetscInt coneSize,i; 561 PetscBool ispe = PETSC_TRUE; 562 563 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 564 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 565 for (i=0;i<coneSize;i++) { 566 PetscInt v; 567 568 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 569 ispe = (PetscBool)(ispe && (v==2)); 570 } 571 if (ispe && coneSize) { 572 ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr); 573 } 574 } 575 } 576 } 577 } 578 for (p=fStart;p<fEnd;p++) { 579 const PetscInt *support; 580 PetscInt supportSize; 581 PetscBool isbf = PETSC_FALSE; 582 583 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 584 if (ovl) { 585 PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE; 586 PetscInt i; 587 588 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 589 for (i=0;i<supportSize;i++) { 590 if (PetscLikely(PetscBTLookup(pown,support[i]-cStart))) has_owned = PETSC_TRUE; 591 else has_ghost = PETSC_TRUE; 592 } 593 isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost)); 594 } else { 595 isbf = (PetscBool)(supportSize == 1); 596 } 597 if (!isbf && perLabel) { 598 const PetscInt *cone; 599 PetscInt coneSize,i; 600 601 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 602 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 603 isbf = PETSC_TRUE; 604 for (i=0;i<coneSize;i++) { 605 PetscInt v,d; 606 607 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 608 ierr = DMLabelGetDefaultValue(perLabel,&d);CHKERRQ(ierr); 609 isbf = (PetscBool)(isbf && v != d); 610 } 611 } 612 if (isbf) { 613 ierr = PetscBTSet(bfaces,p-fStart);CHKERRQ(ierr); 614 } 615 } 616 /* count boundary faces */ 617 for (p=fStart,bf=0;p<fEnd;p++) { 618 if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) { 619 const PetscInt *support; 620 PetscInt supportSize,c; 621 622 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 623 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 624 if (ovl) { 625 for (c=0;c<supportSize;c++) { 626 if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) { 627 bf++; 628 } 629 } 630 } else bf += supportSize; 631 } 632 } 633 /* TODO: is this the label we want to use for marking boundary facets? 634 We should decide to have a single marker for these stuff 635 Proposal: DMSetBoundaryLabel? 636 */ 637 ierr = DMGetLabel(dm,"marker",&label);CHKERRQ(ierr); 638 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",bf);CHKERRQ(ierr); 639 for (p=fStart;p<fEnd;p++) { 640 if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) { 641 const PetscInt *support; 642 PetscInt supportSize,c,nc = 0; 643 644 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 645 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 646 if (ovl) { 647 for (c=0;c<supportSize;c++) { 648 if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) { 649 fcells[nc++] = support[c]; 650 } 651 } 652 } else for (c=0;c<supportSize;c++) fcells[nc++] = support[c]; 653 for (c=0;c<nc;c++) { 654 const PetscInt *cone; 655 int vids[8]; 656 PetscInt i,cell,cl,nv,cid = -1,mid = -1; 657 658 cell = fcells[c]; 659 ierr = DMPlexGetCone(dm,cell,&cone);CHKERRQ(ierr); 660 for (cl=0;cl<fpc;cl++) 661 if (cone[cl] == p) 662 break; 663 664 /* face material id and type */ 665 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);CHKERRQ(ierr); 666 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 667 /* vertex ids */ 668 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,cell,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 669 for (i=0;i<vpf;i++) { 670 ierr = PetscViewerASCIIPrintf(viewer," %D",vids[faces[cl*vpf+i]]);CHKERRQ(ierr); 671 } 672 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 673 } 674 bf = bf-nc; 675 } 676 } 677 if (bf) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Remaining boundary faces %D",bf); 678 ierr = PetscBTDestroy(&bfaces);CHKERRQ(ierr); 679 ierr = PetscFree(fcells);CHKERRQ(ierr); 680 } 681 682 /* mark owned vertices */ 683 if (ovl) { 684 ierr = PetscBTMemzero(vEnd-vStart,pown);CHKERRQ(ierr); 685 for (p=cStart;p<cEnd;p++) { 686 PetscInt i,closureSize,*closure = NULL; 687 688 if (gNum[p-cStart] < 0) continue; 689 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 690 for (i=0;i<closureSize;i++) { 691 const PetscInt pp = closure[2*i]; 692 693 if (pp >= vStart && pp < vEnd) { 694 ierr = PetscBTSet(pown,pp-vStart);CHKERRQ(ierr); 695 } 696 } 697 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 698 } 699 } 700 if (globalNum) { 701 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 702 } 703 704 /* vertex_parents (Non-conforming meshes) */ 705 parentSection = NULL; 706 if (enable_ncmesh) { 707 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 708 } 709 if (parentSection) { 710 PetscInt vp,gvp; 711 712 for (vp=0,p=vStart;p<vEnd;p++) { 713 DMLabel dlabel; 714 PetscInt parent,depth; 715 716 if (PetscUnlikely(ovl && !PetscBTLookup(pown,p-vStart))) continue; 717 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 718 ierr = DMLabelGetValue(dlabel,p,&depth);CHKERRQ(ierr); 719 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 720 if (parent != p) vp++; 721 } 722 ierr = MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 723 if (gvp) { 724 PetscInt maxsupp; 725 PetscBool *skip = NULL; 726 727 ierr = PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");CHKERRQ(ierr); 728 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vp);CHKERRQ(ierr); 729 ierr = DMPlexGetMaxSizes(dm,NULL,&maxsupp);CHKERRQ(ierr); 730 ierr = PetscMalloc1(maxsupp,&skip);CHKERRQ(ierr); 731 for (p=vStart;p<vEnd;p++) { 732 DMLabel dlabel; 733 PetscInt parent; 734 735 if (PetscUnlikely(ovl && !PetscBTLookup(pown,p-vStart))) continue; 736 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 737 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 738 if (parent != p) { 739 int vids[8] = { -1, -1, -1, -1, -1, -1, -1, -1 }; /* silent overzealous clang static analyzer */ 740 PetscInt i,nv,size,n,numChildren,depth = -1; 741 const PetscInt *children; 742 ierr = DMPlexGetConeSize(dm,parent,&size);CHKERRQ(ierr); 743 switch (size) { 744 case 2: /* edge */ 745 nv = 0; 746 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 747 ierr = DMPlexInvertCell(dim,nv,vids);CHKERRQ(ierr); 748 ierr = PetscViewerASCIIPrintf(viewer,"%D",p-vStart);CHKERRQ(ierr); 749 for (i=0;i<nv;i++) { 750 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[i]);CHKERRQ(ierr); 751 } 752 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 753 vp--; 754 break; 755 case 4: /* face */ 756 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 757 for (n=0;n<numChildren;n++) { 758 ierr = DMLabelGetValue(dlabel,children[n],&depth);CHKERRQ(ierr); 759 if (!depth) { 760 const PetscInt *hvsupp,*hesupp,*cone; 761 PetscInt hvsuppSize,hesuppSize,coneSize; 762 PetscInt hv = children[n],he = -1,f; 763 764 ierr = PetscMemzero(skip,maxsupp*sizeof(PetscBool));CHKERRQ(ierr); 765 ierr = DMPlexGetSupportSize(dm,hv,&hvsuppSize);CHKERRQ(ierr); 766 ierr = DMPlexGetSupport(dm,hv,&hvsupp);CHKERRQ(ierr); 767 for (i=0;i<hvsuppSize;i++) { 768 PetscInt ep; 769 ierr = DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);CHKERRQ(ierr); 770 if (ep != hvsupp[i]) { 771 he = hvsupp[i]; 772 } else { 773 skip[i] = PETSC_TRUE; 774 } 775 } 776 if (he == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize); 777 ierr = DMPlexGetCone(dm,he,&cone);CHKERRQ(ierr); 778 vids[0] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 779 ierr = DMPlexGetSupportSize(dm,he,&hesuppSize);CHKERRQ(ierr); 780 ierr = DMPlexGetSupport(dm,he,&hesupp);CHKERRQ(ierr); 781 for (f=0;f<hesuppSize;f++) { 782 PetscInt j; 783 784 ierr = DMPlexGetCone(dm,hesupp[f],&cone);CHKERRQ(ierr); 785 ierr = DMPlexGetConeSize(dm,hesupp[f],&coneSize);CHKERRQ(ierr); 786 for (j=0;j<coneSize;j++) { 787 PetscInt k; 788 for (k=0;k<hvsuppSize;k++) { 789 if (hvsupp[k] == cone[j]) { 790 skip[k] = PETSC_TRUE; 791 break; 792 } 793 } 794 } 795 } 796 for (i=0;i<hvsuppSize;i++) { 797 if (!skip[i]) { 798 ierr = DMPlexGetCone(dm,hvsupp[i],&cone);CHKERRQ(ierr); 799 vids[1] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 800 } 801 } 802 ierr = PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);CHKERRQ(ierr); 803 for (i=0;i<2;i++) { 804 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[i]-vStart);CHKERRQ(ierr); 805 } 806 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 807 vp--; 808 } 809 } 810 break; 811 default: 812 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %d",size); 813 } 814 } 815 } 816 ierr = PetscFree(skip);CHKERRQ(ierr); 817 } 818 if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp); 819 } 820 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 821 822 /* vertices */ 823 if (hovec) { /* higher-order meshes */ 824 const char *fec; 825 PetscInt i,n; 826 827 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 828 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);CHKERRQ(ierr); 829 ierr = PetscViewerASCIIPrintf(viewer,"nodes\n");CHKERRQ(ierr); 830 ierr = PetscObjectGetName((PetscObject)hovec,&fec);CHKERRQ(ierr); 831 ierr = PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");CHKERRQ(ierr); 832 ierr = PetscViewerASCIIPrintf(viewer,"FiniteElementCollection: %s\n",fec);CHKERRQ(ierr); 833 ierr = PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);CHKERRQ(ierr); 834 ierr = PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n");CHKERRQ(ierr); /*Ordering::byVDIM*/ 835 ierr = VecGetArrayRead(hovec,&array);CHKERRQ(ierr); 836 ierr = VecGetLocalSize(hovec,&n);CHKERRQ(ierr); 837 if (n%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Size of local coordinate vector %D incompatible with space dimension %D",n,sdim); 838 for (i=0;i<n/sdim;i++) { 839 PetscInt s; 840 for (s=0;s<sdim;s++) { 841 ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[i*sdim+s]));CHKERRQ(ierr); 842 } 843 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 844 } 845 ierr = VecRestoreArrayRead(hovec,&array);CHKERRQ(ierr); 846 } else { 847 ierr = VecGetLocalSize(coordinates,&nvert);CHKERRQ(ierr); 848 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 849 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);CHKERRQ(ierr); 850 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 851 ierr = VecGetArrayRead(coordinates,&array);CHKERRQ(ierr); 852 for (p=0;p<nvert/sdim;p++) { 853 PetscInt s; 854 for (s=0;s<sdim;s++) { 855 ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[p*sdim+s]));CHKERRQ(ierr); 856 } 857 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 858 } 859 ierr = VecRestoreArrayRead(coordinates,&array);CHKERRQ(ierr); 860 } 861 PetscFunctionReturn(0); 862 } 863 864 /* dispatching, prints through the socket by prepending the mesh keyword to the usual ASCII dump */ 865 PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer) 866 { 867 PetscErrorCode ierr; 868 PetscBool isglvis,isascii; 869 870 PetscFunctionBegin; 871 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 872 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 873 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);CHKERRQ(ierr); 874 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 875 if (!isglvis && !isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERGLVIS or VIEWERASCII"); 876 if (isglvis) { 877 PetscViewer view; 878 PetscViewerGLVisType type; 879 880 ierr = PetscViewerGLVisGetType_Private(viewer,&type);CHKERRQ(ierr); 881 ierr = PetscViewerGLVisGetDMWindow_Private(viewer,&view);CHKERRQ(ierr); 882 if (view) { /* in the socket case, it may happen that the connection failed */ 883 if (type == PETSC_VIEWER_GLVIS_SOCKET) { 884 PetscMPIInt size,rank; 885 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr); 886 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 887 ierr = PetscViewerASCIIPrintf(view,"parallel %D %D\nmesh\n",size,rank);CHKERRQ(ierr); 888 } 889 ierr = DMPlexView_GLVis_ASCII(dm,view);CHKERRQ(ierr); 890 ierr = PetscViewerFlush(view);CHKERRQ(ierr); 891 if (type == PETSC_VIEWER_GLVIS_SOCKET) { 892 PetscInt dim; 893 const char* name; 894 895 ierr = PetscObjectGetName((PetscObject)dm,&name);CHKERRQ(ierr); 896 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 897 ierr = PetscViewerGLVisInitWindow_Private(view,PETSC_TRUE,dim,name);CHKERRQ(ierr); 898 ierr = PetscBarrier((PetscObject)dm);CHKERRQ(ierr); 899 } 900 } 901 } else { 902 ierr = DMPlexView_GLVis_ASCII(dm,viewer);CHKERRQ(ierr); 903 } 904 PetscFunctionReturn(0); 905 } 906