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