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,*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 (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 #undef __FUNCT__ 226 #define __FUNCT__ "DMPlexGetPointMFEMCellID_Internal" 227 PETSC_STATIC_INLINE PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt p, PetscInt *mid, PetscInt *cid) 228 { 229 DMLabel dlabel; 230 PetscInt depth,csize; 231 PetscErrorCode ierr; 232 233 PetscFunctionBegin; 234 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 235 ierr = DMLabelGetValue(dlabel,p,&depth);CHKERRQ(ierr); 236 ierr = DMPlexGetConeSize(dm,p,&csize);CHKERRQ(ierr); 237 if (label) { 238 ierr = DMLabelGetValue(label,p,mid);CHKERRQ(ierr); 239 } else { 240 *mid = 1; 241 } 242 *cid = mfem_table_cid[depth][csize]; 243 PetscFunctionReturn(0); 244 } 245 246 #undef __FUNCT__ 247 #define __FUNCT__ "DMPlexGetPointMFEMVertexIDs_Internal" 248 PETSC_STATIC_INLINE PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscInt *nv, int vids[]) 249 { 250 PetscInt dim,i,q,vStart,vEnd,numPoints,*points = NULL; 251 PetscErrorCode ierr; 252 253 PetscFunctionBegin; 254 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 255 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 256 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 257 for (i=0,q=0;i<numPoints*2;i+= 2) 258 if ((points[i] >= vStart) && (points[i] < vEnd)) 259 vids[q++] = points[i]-vStart; 260 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 261 ierr = DMPlexInvertCell(dim,q,vids);CHKERRQ(ierr); 262 *nv = q; 263 PetscFunctionReturn(0); 264 } 265 266 /* ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR */ 267 static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer) 268 { 269 DM cdm; 270 DMLabel label; 271 PetscSection coordSection,parentSection; 272 Vec coordinates; 273 IS globalNum = NULL; 274 const PetscScalar *array; 275 const PetscInt *gNum; 276 PetscInt bf,p,sdim,dim,depth,novl; 277 PetscInt cStart,cEnd,cEndInterior,fStart,fEnd,fEndInterior,vStart,vEnd; 278 PetscMPIInt commsize; 279 PetscBool ovl,isascii,enable_boundary,enable_ncmesh; 280 PetscBT pown; 281 PetscErrorCode ierr; 282 PetscInt *faces = NULL,fpc = 0,vpf = 0; 283 PetscInt fv1[] = {0,1}, 284 fv2tri[] = {0,1, 285 1,2, 286 2,0}, 287 fv2quad[] = {0,1, 288 1,2, 289 2,3, 290 3,0}, 291 fv3tet[] = {0,1,2, 292 0,3,1, 293 0,2,3, 294 2,1,3}, 295 fv3hex[] = {0,1,2,3, 296 4,5,6,7, 297 0,3,5,4, 298 2,1,7,6, 299 3,2,6,5, 300 0,4,7,1}; 301 PetscContainer glvis_container; 302 PetscBool enabled = PETSC_TRUE; 303 304 PetscFunctionBegin; 305 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 306 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 307 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 308 if (!isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERASCII"); 309 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&commsize);CHKERRQ(ierr); 310 if (commsize > 1) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Use single sequential viewers for parallel visualization"); 311 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 312 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 313 if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 314 315 /* get container: determines if a process visualizes is portion of the data or not */ 316 ierr = PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);CHKERRQ(ierr); 317 if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis container"); 318 { 319 PetscViewerGLVisInfo glvis_info; 320 ierr = PetscContainerGetPointer(glvis_container,(void**)&glvis_info);CHKERRQ(ierr); 321 enabled = glvis_info->enabled; 322 } 323 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, &fEndInterior, NULL, NULL);CHKERRQ(ierr); 324 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 325 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 326 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 327 328 /* 329 a couple of sections of the mesh specification are disabled 330 - boundary: unless we want to visualize boundary attributes, the boundary is not needed for proper mesh visualization 331 - 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 332 */ 333 enable_boundary = PETSC_FALSE; 334 enable_ncmesh = PETSC_FALSE; 335 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"GLVis PetscViewer DMPlex Options","PetscViewer");CHKERRQ(ierr); 336 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); 337 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); 338 ierr = PetscOptionsEnd();CHKERRQ(ierr); 339 340 /* Identify possible cells in the overlap */ 341 gNum = NULL; 342 novl = 0; 343 ovl = PETSC_FALSE; 344 pown = NULL; 345 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr); 346 if (commsize > 1) { 347 ierr = DMPlexGetCellNumbering(dm,&globalNum);CHKERRQ(ierr); 348 ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr); 349 for (p=cStart; p<cEnd; p++) { 350 if (gNum[p-cStart] < 0) { 351 ovl = PETSC_TRUE; 352 novl++; 353 } 354 } 355 if (ovl) { 356 /* it may happen that pown get not destroyed, if the user closes the window while this function is running. 357 TODO: garbage collector? attach pown to dm? */ 358 ierr = PetscBTCreate(PetscMax(cEnd-cStart,vEnd-vStart),&pown);CHKERRQ(ierr); 359 } 360 } 361 362 if (!enabled) { 363 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 364 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 365 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 366 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 367 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 368 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 369 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 370 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 371 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 372 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 373 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 374 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 375 if (globalNum) { 376 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 377 } 378 PetscFunctionReturn(0); 379 } 380 381 /* header */ 382 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 383 384 /* topological dimension */ 385 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 386 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 387 388 /* elements */ 389 /* TODO: is this the label we want to use for marking material IDs? 390 We should decide to have a single marker for these stuff 391 Proposal: DMSetMaterialLabel? 392 */ 393 ierr = DMGetLabel(dm,"Cell sets",&label);CHKERRQ(ierr); 394 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 395 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);CHKERRQ(ierr); 396 for (p=cStart;p<cEnd;p++) { 397 int vids[8]; 398 PetscInt i,nv = 0,cid = -1,mid = 1; 399 400 if (ovl) { 401 if (gNum[p-cStart] < 0) continue; 402 else { 403 ierr = PetscBTSet(pown,p-cStart);CHKERRQ(ierr); 404 } 405 } 406 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);CHKERRQ(ierr); 407 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,p,&nv,vids);CHKERRQ(ierr); 408 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 409 for (i=0;i<nv;i++) { 410 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[i]);CHKERRQ(ierr); 411 } 412 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 413 } 414 415 /* determine orientation of boundary mesh */ 416 if (cEnd-cStart) { 417 ierr = DMPlexGetConeSize(dm,cStart,&fpc);CHKERRQ(ierr); 418 switch(dim) { 419 case 1: 420 if (fpc != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case faces per cell %D",fpc); 421 faces = fv1; 422 vpf = 1; 423 break; 424 case 2: 425 switch (fpc) { 426 case 3: 427 faces = fv2tri; 428 vpf = 2; 429 break; 430 case 4: 431 faces = fv2quad; 432 vpf = 2; 433 break; 434 default: 435 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 436 break; 437 } 438 break; 439 case 3: 440 switch (fpc) { 441 case 4: 442 faces = fv3tet; 443 vpf = 3; 444 break; 445 case 6: 446 faces = fv3hex; 447 vpf = 4; 448 break; 449 default: 450 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 451 break; 452 } 453 break; 454 default: 455 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim"); 456 break; 457 } 458 } 459 460 /* boundary */ 461 ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);CHKERRQ(ierr); 462 fEnd = fEndInterior < 0 ? fEnd : fEndInterior; 463 if (!ovl) { 464 for (p=fStart,bf=0;p<fEnd;p++) { 465 PetscInt supportSize; 466 467 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 468 if (supportSize == 1) bf++; 469 } 470 } else { 471 for (p=fStart,bf=0;p<fEnd;p++) { 472 const PetscInt *support; 473 PetscInt i,supportSize; 474 PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE; 475 476 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 477 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 478 for (i=0;i<supportSize;i++) { 479 if (PetscBTLookup(pown,support[i]-cStart)) has_owned = PETSC_TRUE; 480 else has_ghost = PETSC_TRUE; 481 } 482 if ((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost)) bf++; 483 } 484 } 485 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 486 if (!enable_boundary) { 487 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 488 } else { 489 /* TODO: is this the label we want to use for marking boundary facets? 490 We should decide to have a single marker for these stuff 491 Proposal: DMSetBoundaryLabel? 492 */ 493 ierr = DMGetLabel(dm,"marker",&label);CHKERRQ(ierr); 494 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",bf);CHKERRQ(ierr); 495 if (!ovl) { 496 for (p=fStart;p<fEnd;p++) { 497 PetscInt supportSize; 498 499 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 500 if (supportSize == 1) { 501 const PetscInt *support, *cone; 502 PetscInt i,c,v,cid = -1,mid = -1; 503 PetscInt cellClosureSize,*cellClosure = NULL; 504 505 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 506 ierr = DMPlexGetCone(dm,support[0],&cone);CHKERRQ(ierr); 507 for (c=0;c<fpc;c++) 508 if (cone[c] == p) 509 break; 510 511 ierr = DMPlexGetTransitiveClosure(dm,support[0],PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 512 for (v=0;v<cellClosureSize;v++) 513 if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) 514 break; 515 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);CHKERRQ(ierr); 516 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 517 for (i=0;i<vpf;i++) { 518 ierr = PetscViewerASCIIPrintf(viewer," %D",cellClosure[2*(v+faces[c*vpf+i])]-vStart);CHKERRQ(ierr); 519 } 520 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 521 ierr = DMPlexRestoreTransitiveClosure(dm,support[0],PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 522 } 523 } 524 } else { 525 for (p=fStart;p<fEnd;p++) { 526 const PetscInt *support; 527 PetscInt i,supportSize,validcell = -1; 528 PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE; 529 530 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 531 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 532 for (i=0;i<supportSize;i++) { 533 if (PetscBTLookup(pown,support[i]-cStart)) { 534 has_owned = PETSC_TRUE; 535 validcell = support[i]; 536 } else has_ghost = PETSC_TRUE; 537 } 538 if ((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost)) { 539 const PetscInt *support, *cone; 540 PetscInt i,c,v,cid = -1,mid = -1; 541 PetscInt cellClosureSize,*cellClosure = NULL; 542 543 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 544 ierr = DMPlexGetCone(dm,validcell,&cone);CHKERRQ(ierr); 545 for (c=0;c<fpc;c++) 546 if (cone[c] == p) 547 break; 548 549 ierr = DMPlexGetTransitiveClosure(dm,validcell,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 550 for (v=0;v<cellClosureSize;v++) 551 if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) 552 break; 553 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);CHKERRQ(ierr); 554 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 555 for (i=0;i<vpf;i++) { 556 ierr = PetscViewerASCIIPrintf(viewer," %D",cellClosure[2*(v+faces[c*vpf+i])]-vStart);CHKERRQ(ierr); 557 } 558 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 559 ierr = DMPlexRestoreTransitiveClosure(dm,validcell,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 560 } 561 } 562 } 563 } 564 565 /* mark owned vertices */ 566 if (ovl) { 567 ierr = PetscBTMemzero(vEnd-vStart,pown);CHKERRQ(ierr); 568 for (p=cStart;p<cEnd;p++) { 569 PetscInt i,closureSize,*closure = NULL; 570 571 if (gNum[p-cStart] < 0) continue; 572 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 573 for (i=0;i<closureSize;i++) { 574 const PetscInt pp = closure[2*i]; 575 576 if (pp >= vStart && pp < vEnd) { 577 ierr = PetscBTSet(pown,pp-vStart);CHKERRQ(ierr); 578 } 579 } 580 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 581 } 582 } 583 if (globalNum) { 584 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 585 } 586 587 /* vertex_parents (Non-conforming meshes) */ 588 parentSection = NULL; 589 if (enable_ncmesh) { 590 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 591 } 592 if (parentSection) { 593 PetscInt vp,gvp; 594 595 for (vp=0,p=vStart;p<vEnd;p++) { 596 DMLabel dlabel; 597 PetscInt parent,depth; 598 599 if (PetscUnlikely(ovl && !PetscBTLookup(pown,p-vStart))) continue; 600 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 601 ierr = DMLabelGetValue(dlabel,p,&depth);CHKERRQ(ierr); 602 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 603 if (parent != p) vp++; 604 } 605 ierr = MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 606 if (gvp) { 607 PetscInt maxsupp; 608 PetscBool *skip = NULL; 609 610 ierr = PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");CHKERRQ(ierr); 611 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vp);CHKERRQ(ierr); 612 ierr = DMPlexGetMaxSizes(dm,NULL,&maxsupp);CHKERRQ(ierr); 613 ierr = PetscMalloc1(maxsupp,&skip);CHKERRQ(ierr); 614 for (p=vStart;p<vEnd;p++) { 615 DMLabel dlabel; 616 PetscInt parent; 617 618 if (PetscUnlikely(ovl && !PetscBTLookup(pown,p-vStart))) continue; 619 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 620 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 621 if (parent != p) { 622 int vids[8]; 623 PetscInt i,nv,size,n,numChildren,depth = -1; 624 const PetscInt *children; 625 ierr = DMPlexGetConeSize(dm,parent,&size);CHKERRQ(ierr); 626 switch (size) { 627 case 2: /* edge */ 628 nv = 0; 629 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,&nv,vids);CHKERRQ(ierr); 630 ierr = PetscViewerASCIIPrintf(viewer,"%D",p-vStart);CHKERRQ(ierr); 631 for (i=0;i<nv;i++) { 632 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[i]);CHKERRQ(ierr); 633 } 634 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 635 vp--; 636 break; 637 case 4: /* face */ 638 ierr = PetscMemzero(skip,maxsupp*sizeof(PetscBool));CHKERRQ(ierr); 639 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 640 for (n=0;n<numChildren;n++) { 641 ierr = DMLabelGetValue(dlabel,children[n],&depth);CHKERRQ(ierr); 642 if (!depth) { 643 const PetscInt *hvsupp,*hesupp,*cone; 644 PetscInt hvsuppSize,hesuppSize,coneSize; 645 PetscInt hv = children[n],he,f; 646 PetscBool found = PETSC_FALSE; 647 648 ierr = DMPlexGetSupportSize(dm,hv,&hvsuppSize);CHKERRQ(ierr); 649 ierr = DMPlexGetSupport(dm,hv,&hvsupp);CHKERRQ(ierr); 650 for (i=0;i<hvsuppSize;i++) { 651 PetscInt ep; 652 ierr = DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);CHKERRQ(ierr); 653 if (ep != hvsupp[i]) { 654 he = hvsupp[i]; 655 found = PETSC_TRUE; 656 } else { 657 skip[i] = PETSC_TRUE; 658 } 659 } 660 if (!found) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize); 661 he = hvsupp[i]; 662 skip[i] = PETSC_TRUE; 663 ierr = DMPlexGetCone(dm,he,&cone);CHKERRQ(ierr); 664 vids[0] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 665 ierr = DMPlexGetSupportSize(dm,he,&hesuppSize);CHKERRQ(ierr); 666 ierr = DMPlexGetSupport(dm,he,&hesupp);CHKERRQ(ierr); 667 for (f=0;f<hesuppSize;f++) { 668 PetscInt j; 669 670 ierr = DMPlexGetCone(dm,hesupp[f],&cone);CHKERRQ(ierr); 671 ierr = DMPlexGetConeSize(dm,hesupp[f],&coneSize);CHKERRQ(ierr); 672 for (j=0;j<coneSize;j++) { 673 PetscInt k; 674 for (k=0;k<hvsuppSize;k++) { 675 if (hvsupp[k] == cone[j]) { 676 skip[k] = PETSC_TRUE; 677 break; 678 } 679 } 680 } 681 } 682 for (i=0;i<hvsuppSize;i++) { 683 if (!skip[i]) { 684 ierr = DMPlexGetCone(dm,hvsupp[i],&cone);CHKERRQ(ierr); 685 vids[1] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 686 } 687 } 688 ierr = PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);CHKERRQ(ierr); 689 for (i=0;i<2;i++) { 690 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[i]-vStart);CHKERRQ(ierr); 691 } 692 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 693 vp--; 694 } 695 } 696 break; 697 default: 698 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %d",size); 699 } 700 } 701 } 702 ierr = PetscFree(skip);CHKERRQ(ierr); 703 } 704 if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp); 705 } 706 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 707 708 /* vertices */ 709 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 710 ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr); 711 ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 712 ierr = DMGetCoordinatesLocal(dm,&coordinates);CHKERRQ(ierr); 713 ierr = DMPlexGetDepthStratum(cdm,0,&vStart,&vEnd);CHKERRQ(ierr); 714 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 715 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);CHKERRQ(ierr); 716 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 717 ierr = VecGetArrayRead(coordinates,&array);CHKERRQ(ierr); 718 for (p=vStart;p<vEnd;p++) { 719 PetscInt i,st,dof; 720 721 ierr = PetscSectionGetDof(coordSection,p,&dof);CHKERRQ(ierr); 722 ierr = PetscSectionGetOffset(coordSection,p,&st);CHKERRQ(ierr); 723 for (i=st;i<st+dof-1;i++) { 724 ierr = PetscViewerASCIIPrintf(viewer,"%g ",PetscRealPart(array[i]));CHKERRQ(ierr); 725 } 726 ierr = PetscViewerASCIIPrintf(viewer,"%g\n",PetscRealPart(array[i]));CHKERRQ(ierr); 727 } 728 ierr = VecRestoreArrayRead(coordinates,&array);CHKERRQ(ierr); 729 PetscFunctionReturn(0); 730 } 731 732 /* dispatching, prints through the socket by prepending the mesh keyword to the usual ASCII dump */ 733 PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer) 734 { 735 PetscErrorCode ierr; 736 PetscBool isglvis,isascii; 737 738 PetscFunctionBegin; 739 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 740 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 741 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);CHKERRQ(ierr); 742 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 743 if (!isglvis && !isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERGLVIS or VIEWERASCII"); 744 if (isglvis) { 745 PetscViewer view; 746 PetscViewerGLVisType type; 747 748 ierr = PetscViewerGLVisGetType_Private(viewer,&type);CHKERRQ(ierr); 749 ierr = PetscViewerGLVisGetDMWindow_Private(viewer,&view);CHKERRQ(ierr); 750 if (view) { /* in the socket case, it may happen that the connection failed */ 751 if (type == PETSC_VIEWER_GLVIS_SOCKET) { 752 PetscMPIInt size,rank; 753 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr); 754 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 755 ierr = PetscViewerASCIIPrintf(view,"parallel %D %D\nmesh\n",size,rank);CHKERRQ(ierr); 756 } 757 ierr = DMPlexView_GLVis_ASCII(dm,view);CHKERRQ(ierr); 758 ierr = PetscViewerFlush(view);CHKERRQ(ierr); 759 if (type == PETSC_VIEWER_GLVIS_SOCKET) { 760 PetscInt dim; 761 const char* name; 762 763 ierr = PetscObjectGetName((PetscObject)dm,&name);CHKERRQ(ierr); 764 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 765 ierr = PetscViewerGLVisInitWindow_Private(view,PETSC_TRUE,dim,name);CHKERRQ(ierr); 766 ierr = PetscBarrier((PetscObject)dm);CHKERRQ(ierr); 767 } 768 } 769 } else { 770 ierr = DMPlexView_GLVis_ASCII(dm,viewer);CHKERRQ(ierr); 771 } 772 PetscFunctionReturn(0); 773 } 774