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,*Ufield; 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 = DMGetSection(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 = PetscCalloc7(maxfields,&fieldname,maxfields,&nlocal,maxfields,&bs,maxfields,&dims,maxfields,&fec_type,totdofs,&idxs,maxfields,&Ufield);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 = VecScatterCreateWithData(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 = VecScatterCreateWithData(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 /* create work vectors */ 209 for (f=0;f<ctx->nf;f++) { 210 ierr = VecCreateMPI(PetscObjectComm((PetscObject)dm),nlocal[f],PETSC_DECIDE,&Ufield[f]);CHKERRQ(ierr); 211 ierr = PetscObjectSetName((PetscObject)Ufield[f],fieldname[f]);CHKERRQ(ierr); 212 ierr = VecSetBlockSize(Ufield[f],bs[f]);CHKERRQ(ierr); 213 ierr = VecSetDM(Ufield[f],dm);CHKERRQ(ierr); 214 } 215 216 /* customize the viewer */ 217 ierr = PetscViewerGLVisSetFields(viewer,ctx->nf,(const char**)fec_type,dims,DMPlexSampleGLVisFields_Private,(PetscObject*)Ufield,ctx,DestroyGLVisViewerCtx_Private);CHKERRQ(ierr); 218 for (f=0;f<ctx->nf;f++) { 219 ierr = PetscFree(fieldname[f]);CHKERRQ(ierr); 220 ierr = PetscFree(fec_type[f]);CHKERRQ(ierr); 221 ierr = VecDestroy(&Ufield[f]);CHKERRQ(ierr); 222 } 223 ierr = PetscFree7(fieldname,nlocal,bs,dims,fec_type,idxs,Ufield);CHKERRQ(ierr); 224 PetscFunctionReturn(0); 225 } 226 227 typedef enum {MFEM_POINT=0,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE,MFEM_TETRAHEDRON,MFEM_CUBE,MFEM_UNDEF} MFEM_cid; 228 229 MFEM_cid mfem_table_cid[4][7] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF}, 230 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF}, 231 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE ,MFEM_UNDEF,MFEM_UNDEF}, 232 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_TETRAHEDRON,MFEM_UNDEF,MFEM_CUBE } }; 233 234 MFEM_cid mfem_table_cid_unint[4][9] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF}, 235 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF}, 236 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE ,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF}, 237 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_TETRAHEDRON,MFEM_UNDEF,MFEM_UNDEF,MFEM_UNDEF,MFEM_CUBE } }; 238 239 static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt p, PetscInt *mid, PetscInt *cid) 240 { 241 DMLabel dlabel; 242 PetscInt depth,csize,pdepth,dim; 243 PetscErrorCode ierr; 244 245 PetscFunctionBegin; 246 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 247 ierr = DMLabelGetValue(dlabel,p,&pdepth);CHKERRQ(ierr); 248 ierr = DMPlexGetConeSize(dm,p,&csize);CHKERRQ(ierr); 249 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 250 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 251 if (label) { 252 ierr = DMLabelGetValue(label,p,mid);CHKERRQ(ierr); 253 } else *mid = 1; 254 if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */ 255 #if defined PETSC_USE_DEBUG 256 if (dim < 0 || dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D",dim); 257 if (csize > 8) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found cone size %D for point %D",csize,p); 258 if (depth != 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found depth %D for point %D. You should interpolate the mesh first",depth,p); 259 #endif 260 *cid = mfem_table_cid_unint[dim][csize]; 261 } else { 262 #if defined PETSC_USE_DEBUG 263 if (csize > 6) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cone size %D for point %D",csize,p); 264 if (pdepth < 0 || pdepth > 3) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Depth %D for point %D",csize,p); 265 #endif 266 *cid = mfem_table_cid[pdepth][csize]; 267 } 268 PetscFunctionReturn(0); 269 } 270 271 static PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscSection csec, PetscInt *nv, int vids[]) 272 { 273 PetscInt dim,sdim,dof = 0,off = 0,i,q,vStart,vEnd,numPoints,*points = NULL; 274 PetscErrorCode ierr; 275 276 PetscFunctionBegin; 277 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 278 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 279 sdim = dim; 280 if (csec) { 281 PetscInt sStart,sEnd; 282 283 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 284 ierr = PetscSectionGetChart(csec,&sStart,&sEnd);CHKERRQ(ierr); 285 ierr = PetscSectionGetOffset(csec,vStart,&off);CHKERRQ(ierr); 286 off = off/sdim; 287 if (p >= sStart && p < sEnd) { 288 ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr); 289 } 290 } 291 if (!dof) { 292 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 293 for (i=0,q=0;i<numPoints*2;i+= 2) 294 if ((points[i] >= vStart) && (points[i] < vEnd)) 295 vids[q++] = (int)(points[i]-vStart+off); 296 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 297 } else { 298 ierr = PetscSectionGetOffset(csec,p,&off);CHKERRQ(ierr); 299 ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr); 300 for (q=0;q<dof/sdim;q++) vids[q] = (int)(off/sdim + q); 301 } 302 *nv = q; 303 PetscFunctionReturn(0); 304 } 305 306 /* 307 ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR 308 Higher order meshes are also supported 309 */ 310 static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer) 311 { 312 DMLabel label; 313 PetscSection coordSection,parentSection; 314 Vec coordinates,hovec; 315 const PetscScalar *array; 316 PetscInt bf,p,sdim,dim,depth,novl; 317 PetscInt cStart,cEnd,cEndInterior,vStart,vEnd,nvert; 318 PetscMPIInt size; 319 PetscBool localized,isascii; 320 PetscBool enable_mfem,enable_boundary,enable_ncmesh; 321 PetscBT pown,vown; 322 PetscErrorCode ierr; 323 PetscContainer glvis_container; 324 PetscBool cellvertex = PETSC_FALSE, periodic, enabled = PETSC_TRUE; 325 const char *fmt; 326 char emark[64] = "",bmark[64] = ""; 327 328 PetscFunctionBegin; 329 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 330 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 331 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 332 if (!isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERASCII"); 333 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&size);CHKERRQ(ierr); 334 if (size > 1) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Use single sequential viewers for parallel visualization"); 335 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 336 337 /* get container: determines if a process visualizes is portion of the data or not */ 338 ierr = PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);CHKERRQ(ierr); 339 if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis container"); 340 { 341 PetscViewerGLVisInfo glvis_info; 342 ierr = PetscContainerGetPointer(glvis_container,(void**)&glvis_info);CHKERRQ(ierr); 343 enabled = glvis_info->enabled; 344 fmt = glvis_info->fmt; 345 } 346 347 /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh 348 DMPlex does not currently support HO meshes, so there's no API for this */ 349 ierr = PetscObjectQuery((PetscObject)dm,"_glvis_mesh_coords",(PetscObject*)&hovec);CHKERRQ(ierr); 350 351 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 352 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 353 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 354 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 355 ierr = DMGetPeriodicity(dm,&periodic,NULL,NULL,NULL);CHKERRQ(ierr); 356 ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr); 357 if (periodic && !localized && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Coordinates need to be localized"); 358 ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 359 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 360 ierr = DMGetCoordinatesLocal(dm,&coordinates);CHKERRQ(ierr); 361 if (!coordinates && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector"); 362 363 /* 364 a couple of sections of the mesh specification are disabled 365 - boundary: the boundary is not needed for proper mesh visualization unless we want to visualize boundary attributes or we have high-order coordinates in 3D (topologically) 366 - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package 367 and be able to derefine the mesh (MFEM does not currently have to ability to read ncmeshes in parallel) 368 */ 369 enable_boundary = PETSC_FALSE; 370 enable_ncmesh = PETSC_FALSE; 371 enable_mfem = PETSC_FALSE; 372 /* I'm tired of problems with negative values in the markers, disable them */ 373 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"GLVis PetscViewer DMPlex Options","PetscViewer");CHKERRQ(ierr); 374 ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary","Enable boundary section in mesh representation",NULL,enable_boundary,&enable_boundary,NULL);CHKERRQ(ierr); 375 ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_ncmesh","Enable vertex_parents section in mesh representation (allows derefinement)",NULL,enable_ncmesh,&enable_ncmesh,NULL);CHKERRQ(ierr); 376 ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_mfem","Dump a mesh that can be used with MFEM's FiniteElementSpaces",NULL,enable_mfem,&enable_mfem,NULL);CHKERRQ(ierr); 377 ierr = PetscOptionsString("-viewer_glvis_dm_plex_emarker","String for the material id label",NULL,emark,emark,sizeof(emark),NULL);CHKERRQ(ierr); 378 ierr = PetscOptionsString("-viewer_glvis_dm_plex_bmarker","String for the boundary id label",NULL,bmark,bmark,sizeof(bmark),NULL);CHKERRQ(ierr); 379 ierr = PetscOptionsEnd();CHKERRQ(ierr); 380 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr); 381 if (enable_ncmesh && size > 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not supported in parallel"); 382 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 383 if (enable_boundary && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. " 384 "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0"); 385 if (enable_ncmesh && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. " 386 "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0"); 387 if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */ 388 if (depth != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported depth %D. You should interpolate the mesh first",depth); 389 cellvertex = PETSC_TRUE; 390 } 391 392 /* Identify possible cells in the overlap */ 393 novl = 0; 394 pown = NULL; 395 if (size > 1) { 396 IS globalNum = NULL; 397 const PetscInt *gNum; 398 PetscBool ovl = PETSC_FALSE; 399 400 ierr = DMPlexGetCellNumbering(dm,&globalNum);CHKERRQ(ierr); 401 ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr); 402 for (p=cStart; p<cEnd; p++) { 403 if (gNum[p-cStart] < 0) { 404 ovl = PETSC_TRUE; 405 novl++; 406 } 407 } 408 if (ovl) { 409 /* it may happen that pown get not destroyed, if the user closes the window while this function is running. 410 TODO: garbage collector? attach pown to dm? */ 411 ierr = PetscBTCreate(cEnd-cStart,&pown);CHKERRQ(ierr); 412 for (p=cStart; p<cEnd; p++) { 413 if (gNum[p-cStart] < 0) continue; 414 else { 415 ierr = PetscBTSet(pown,p-cStart);CHKERRQ(ierr); 416 } 417 } 418 } 419 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 420 } 421 422 /* return if this process is disabled */ 423 if (!enabled) { 424 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 425 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 426 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 427 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 428 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 429 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 430 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 431 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 432 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 433 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 434 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437 438 if (enable_mfem) { 439 if (periodic && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */ 440 PetscInt vpc = 0; 441 char fec[64]; 442 int vids[8] = {0,1,2,3,4,5,6,7}; 443 int hexv[8] = {0,1,3,2,4,5,7,6}, tetv[4] = {0,1,2,3}; 444 int quadv[8] = {0,1,3,2}, triv[3] = {0,1,2}; 445 int *dof = NULL; 446 PetscScalar *array,*ptr; 447 448 ierr = PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: L2_T1_%DD_P1",dim);CHKERRQ(ierr); 449 if (cEnd-cStart) { 450 PetscInt fpc; 451 452 ierr = DMPlexGetConeSize(dm,cStart,&fpc);CHKERRQ(ierr); 453 switch(dim) { 454 case 1: 455 vpc = 2; 456 dof = hexv; 457 break; 458 case 2: 459 switch (fpc) { 460 case 3: 461 vpc = 3; 462 dof = triv; 463 break; 464 case 4: 465 vpc = 4; 466 dof = quadv; 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: /* TODO: still need to understand L2 ordering for tets */ 476 vpc = 4; 477 dof = tetv; 478 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled tethraedral case"); 479 break; 480 case 6: 481 if (cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: vertices per cell %D",fpc); 482 vpc = 8; 483 dof = hexv; 484 break; 485 case 8: 486 if (!cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 487 vpc = 8; 488 dof = hexv; 489 break; 490 default: 491 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 492 break; 493 } 494 break; 495 default: 496 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim"); 497 break; 498 } 499 ierr = DMPlexInvertCell(dim,vpc,vids);CHKERRQ(ierr); 500 } 501 if (!dof) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing dofs"); 502 ierr = VecCreateSeq(PETSC_COMM_SELF,(cEnd-cStart-novl)*vpc*sdim,&hovec);CHKERRQ(ierr); 503 ierr = PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)hovec);CHKERRQ(ierr); 504 ierr = PetscObjectDereference((PetscObject)hovec);CHKERRQ(ierr); 505 ierr = PetscObjectSetName((PetscObject)hovec,fec);CHKERRQ(ierr); 506 ierr = VecGetArray(hovec,&array);CHKERRQ(ierr); 507 ptr = array; 508 for (p=cStart;p<cEnd;p++) { 509 PetscInt csize,v,d; 510 PetscScalar *vals = NULL; 511 512 if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue; 513 ierr = DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 514 if (csize != vpc*sdim && csize != vpc*sdim*2) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported closure size %D (vpc %D, sdim %D)",csize,vpc,sdim); 515 for (v=0;v<vpc;v++) { 516 for (d=0;d<sdim;d++) { 517 ptr[sdim*dof[v]+d] = vals[sdim*vids[v]+d]; 518 } 519 } 520 ptr += vpc*sdim; 521 ierr = DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 522 } 523 ierr = VecRestoreArray(hovec,&array);CHKERRQ(ierr); 524 } 525 } 526 /* if we have high-order coordinates in 3D, we need to specify the boundary */ 527 if (hovec && dim == 3) enable_boundary = PETSC_TRUE; 528 529 /* header */ 530 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 531 532 /* topological dimension */ 533 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 534 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 535 536 /* elements */ 537 ierr = DMGetLabel(dm,emark,&label);CHKERRQ(ierr); 538 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 539 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);CHKERRQ(ierr); 540 for (p=cStart;p<cEnd;p++) { 541 int vids[8]; 542 PetscInt i,nv = 0,cid = -1,mid = 1; 543 544 if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue; 545 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);CHKERRQ(ierr); 546 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,p,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 547 ierr = DMPlexInvertCell(dim,nv,vids);CHKERRQ(ierr); 548 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 549 for (i=0;i<nv;i++) { 550 ierr = PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);CHKERRQ(ierr); 551 } 552 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 553 } 554 555 /* boundary */ 556 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 557 if (!enable_boundary) { 558 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 559 } else { 560 DMLabel perLabel; 561 PetscBT bfaces; 562 PetscInt fStart,fEnd,fEndInterior,*fcells; 563 PetscInt *faces = NULL,fpc = 0,vpf = 0, vpc = 0; 564 PetscInt fv1[] = {0,1}, 565 fv2tri[] = {0,1, 566 1,2, 567 2,0}, 568 fv2quad[] = {0,1, 569 1,2, 570 2,3, 571 3,0}, 572 fv3tet[] = {0,1,2, 573 0,3,1, 574 0,2,3, 575 2,1,3}, 576 fv3hex[] = {0,1,2,3, 577 4,5,6,7, 578 0,3,5,4, 579 2,1,7,6, 580 3,2,6,5, 581 0,4,7,1}; 582 583 /* determine orientation of boundary mesh */ 584 if (cEnd-cStart) { 585 ierr = DMPlexGetConeSize(dm,cStart,&fpc);CHKERRQ(ierr); 586 switch(dim) { 587 case 1: 588 if (fpc != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case faces per cell %D",fpc); 589 faces = fv1; 590 vpf = 1; 591 vpc = 2; 592 break; 593 case 2: 594 switch (fpc) { 595 case 3: 596 faces = fv2tri; 597 vpf = 2; 598 vpc = 3; 599 break; 600 case 4: 601 faces = fv2quad; 602 vpf = 2; 603 vpc = 4; 604 break; 605 default: 606 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 607 break; 608 } 609 break; 610 case 3: 611 switch (fpc) { 612 case 4: 613 faces = fv3tet; 614 vpf = 3; 615 vpc = 4; 616 break; 617 case 6: 618 faces = fv3hex; 619 vpf = 4; 620 vpc = 8; 621 break; 622 default: 623 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 624 break; 625 } 626 break; 627 default: 628 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim"); 629 break; 630 } 631 } 632 ierr = DMPlexGetHybridBounds(dm,NULL,&fEndInterior,NULL,NULL);CHKERRQ(ierr); 633 ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);CHKERRQ(ierr); 634 fEnd = fEndInterior < 0 ? fEnd : fEndInterior; 635 ierr = PetscBTCreate(fEnd-fStart,&bfaces);CHKERRQ(ierr); 636 ierr = DMPlexGetMaxSizes(dm,NULL,&p);CHKERRQ(ierr); 637 ierr = PetscMalloc1(p,&fcells);CHKERRQ(ierr); 638 ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr); 639 if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */ 640 ierr = DMCreateLabel(dm,"glvis_periodic_cut");CHKERRQ(ierr); 641 ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr); 642 ierr = DMLabelSetDefaultValue(perLabel,1);CHKERRQ(ierr); 643 for (p=cStart;p<cEnd;p++) { 644 PetscInt dof; 645 ierr = PetscSectionGetDof(coordSection,p,&dof);CHKERRQ(ierr); 646 if (dof) { 647 PetscInt v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL; 648 PetscScalar *vals = NULL; 649 650 if (dof%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Incompatible number of cell dofs %D and space dimension %D",dof,sdim); 651 if (dof/sdim != vpc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Incompatible number of cell dofs %D, vertices %D and space dim %D",dof/sdim,vpc,sdim); 652 ierr = DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 653 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 654 for (v=0;v<cellClosureSize;v++) 655 if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) { 656 vidxs = cellClosure + 2*v; 657 break; 658 } 659 if (!vidxs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing vertices"); 660 for (v=0;v<vpc;v++) { 661 PetscInt s; 662 663 for (s=0;s<sdim;s++) { 664 if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+vpc*sdim])>PETSC_MACHINE_EPSILON) { 665 ierr = DMLabelSetValue(perLabel,vidxs[2*v],2);CHKERRQ(ierr); 666 } 667 } 668 } 669 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 670 ierr = DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 671 } 672 } 673 if (dim > 1) { 674 PetscInt eEnd,eStart,eEndInterior; 675 676 ierr = DMPlexGetHybridBounds(dm,NULL,NULL,&eEndInterior,NULL);CHKERRQ(ierr); 677 ierr = DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);CHKERRQ(ierr); 678 eEnd = eEndInterior < 0 ? eEnd : eEndInterior; 679 for (p=eStart;p<eEnd;p++) { 680 const PetscInt *cone; 681 PetscInt coneSize,i; 682 PetscBool ispe = PETSC_TRUE; 683 684 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 685 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 686 for (i=0;i<coneSize;i++) { 687 PetscInt v; 688 689 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 690 ispe = (PetscBool)(ispe && (v==2)); 691 } 692 if (ispe && coneSize) { 693 PetscInt ch, numChildren; 694 const PetscInt *children; 695 696 ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr); 697 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 698 for (ch = 0; ch < numChildren; ch++) { 699 ierr = DMLabelSetValue(perLabel,children[ch],2);CHKERRQ(ierr); 700 } 701 } 702 } 703 if (dim > 2) { 704 for (p=fStart;p<fEnd;p++) { 705 const PetscInt *cone; 706 PetscInt coneSize,i; 707 PetscBool ispe = PETSC_TRUE; 708 709 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 710 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 711 for (i=0;i<coneSize;i++) { 712 PetscInt v; 713 714 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 715 ispe = (PetscBool)(ispe && (v==2)); 716 } 717 if (ispe && coneSize) { 718 PetscInt ch, numChildren; 719 const PetscInt *children; 720 721 ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr); 722 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 723 for (ch = 0; ch < numChildren; ch++) { 724 ierr = DMLabelSetValue(perLabel,children[ch],2);CHKERRQ(ierr); 725 } 726 } 727 } 728 } 729 } 730 } 731 for (p=fStart;p<fEnd;p++) { 732 const PetscInt *support; 733 PetscInt supportSize; 734 PetscBool isbf = PETSC_FALSE; 735 736 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 737 if (pown) { 738 PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE; 739 PetscInt i; 740 741 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 742 for (i=0;i<supportSize;i++) { 743 if (PetscLikely(PetscBTLookup(pown,support[i]-cStart))) has_owned = PETSC_TRUE; 744 else has_ghost = PETSC_TRUE; 745 } 746 isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost)); 747 } else { 748 isbf = (PetscBool)(supportSize == 1); 749 } 750 if (!isbf && perLabel) { 751 const PetscInt *cone; 752 PetscInt coneSize,i; 753 754 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 755 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 756 isbf = PETSC_TRUE; 757 for (i=0;i<coneSize;i++) { 758 PetscInt v,d; 759 760 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 761 ierr = DMLabelGetDefaultValue(perLabel,&d);CHKERRQ(ierr); 762 isbf = (PetscBool)(isbf && v != d); 763 } 764 } 765 if (isbf) { 766 ierr = PetscBTSet(bfaces,p-fStart);CHKERRQ(ierr); 767 } 768 } 769 /* count boundary faces */ 770 for (p=fStart,bf=0;p<fEnd;p++) { 771 if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) { 772 const PetscInt *support; 773 PetscInt supportSize,c; 774 775 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 776 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 777 for (c=0;c<supportSize;c++) { 778 const PetscInt *cone; 779 PetscInt cell,cl; 780 781 cell = support[c]; 782 if (pown && PetscUnlikely(!PetscBTLookup(pown,cell-cStart))) continue; 783 ierr = DMPlexGetCone(dm,cell,&cone);CHKERRQ(ierr); 784 for (cl=0;cl<fpc;cl++) { 785 if (cone[cl] == p) { 786 bf += 1; 787 break; 788 } 789 } 790 } 791 } 792 } 793 ierr = DMGetLabel(dm,bmark,&label);CHKERRQ(ierr); 794 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",bf);CHKERRQ(ierr); 795 for (p=fStart;p<fEnd;p++) { 796 if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) { 797 const PetscInt *support; 798 PetscInt supportSize,c,nc = 0; 799 800 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 801 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 802 if (pown) { 803 for (c=0;c<supportSize;c++) { 804 if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) { 805 fcells[nc++] = support[c]; 806 } 807 } 808 } else for (c=0;c<supportSize;c++) fcells[nc++] = support[c]; 809 for (c=0;c<nc;c++) { 810 const PetscInt *cone; 811 int vids[8]; 812 PetscInt i,cell,cl,nv,cid = -1,mid = -1; 813 814 cell = fcells[c]; 815 ierr = DMPlexGetCone(dm,cell,&cone);CHKERRQ(ierr); 816 for (cl=0;cl<fpc;cl++) 817 if (cone[cl] == p) 818 break; 819 if (cl == fpc) continue; 820 821 /* face material id and type */ 822 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,p,&mid,&cid);CHKERRQ(ierr); 823 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 824 /* vertex ids */ 825 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,cell,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 826 for (i=0;i<vpf;i++) { 827 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[faces[cl*vpf+i]]);CHKERRQ(ierr); 828 } 829 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 830 bf -= 1; 831 } 832 } 833 } 834 if (bf) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Remaining boundary faces %D",bf); 835 ierr = PetscBTDestroy(&bfaces);CHKERRQ(ierr); 836 ierr = PetscFree(fcells);CHKERRQ(ierr); 837 } 838 839 /* mark owned vertices */ 840 vown = NULL; 841 if (pown) { 842 ierr = PetscBTCreate(vEnd-vStart,&vown);CHKERRQ(ierr); 843 for (p=cStart;p<cEnd;p++) { 844 PetscInt i,closureSize,*closure = NULL; 845 846 if (PetscUnlikely(!PetscBTLookup(pown,p-cStart))) continue; 847 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 848 for (i=0;i<closureSize;i++) { 849 const PetscInt pp = closure[2*i]; 850 851 if (pp >= vStart && pp < vEnd) { 852 ierr = PetscBTSet(vown,pp-vStart);CHKERRQ(ierr); 853 } 854 } 855 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 856 } 857 } 858 859 /* vertex_parents (Non-conforming meshes) */ 860 parentSection = NULL; 861 if (enable_ncmesh) { 862 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 863 } 864 if (parentSection) { 865 PetscInt vp,gvp; 866 867 for (vp=0,p=vStart;p<vEnd;p++) { 868 DMLabel dlabel; 869 PetscInt parent,depth; 870 871 if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue; 872 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 873 ierr = DMLabelGetValue(dlabel,p,&depth);CHKERRQ(ierr); 874 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 875 if (parent != p) vp++; 876 } 877 ierr = MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 878 if (gvp) { 879 PetscInt maxsupp; 880 PetscBool *skip = NULL; 881 882 ierr = PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");CHKERRQ(ierr); 883 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vp);CHKERRQ(ierr); 884 ierr = DMPlexGetMaxSizes(dm,NULL,&maxsupp);CHKERRQ(ierr); 885 ierr = PetscMalloc1(maxsupp,&skip);CHKERRQ(ierr); 886 for (p=vStart;p<vEnd;p++) { 887 DMLabel dlabel; 888 PetscInt parent; 889 890 if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue; 891 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 892 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 893 if (parent != p) { 894 int vids[8] = { -1, -1, -1, -1, -1, -1, -1, -1 }; /* silent overzealous clang static analyzer */ 895 PetscInt i,nv,ssize,n,numChildren,depth = -1; 896 const PetscInt *children; 897 898 ierr = DMPlexGetConeSize(dm,parent,&ssize);CHKERRQ(ierr); 899 switch (ssize) { 900 case 2: /* edge */ 901 nv = 0; 902 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 903 ierr = DMPlexInvertCell(dim,nv,vids);CHKERRQ(ierr); 904 ierr = PetscViewerASCIIPrintf(viewer,"%D",p-vStart);CHKERRQ(ierr); 905 for (i=0;i<nv;i++) { 906 ierr = PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);CHKERRQ(ierr); 907 } 908 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 909 vp--; 910 break; 911 case 4: /* face */ 912 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 913 for (n=0;n<numChildren;n++) { 914 ierr = DMLabelGetValue(dlabel,children[n],&depth);CHKERRQ(ierr); 915 if (!depth) { 916 const PetscInt *hvsupp,*hesupp,*cone; 917 PetscInt hvsuppSize,hesuppSize,coneSize; 918 PetscInt hv = children[n],he = -1,f; 919 920 ierr = PetscMemzero(skip,maxsupp*sizeof(PetscBool));CHKERRQ(ierr); 921 ierr = DMPlexGetSupportSize(dm,hv,&hvsuppSize);CHKERRQ(ierr); 922 ierr = DMPlexGetSupport(dm,hv,&hvsupp);CHKERRQ(ierr); 923 for (i=0;i<hvsuppSize;i++) { 924 PetscInt ep; 925 ierr = DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);CHKERRQ(ierr); 926 if (ep != hvsupp[i]) { 927 he = hvsupp[i]; 928 } else { 929 skip[i] = PETSC_TRUE; 930 } 931 } 932 if (he == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize); 933 ierr = DMPlexGetCone(dm,he,&cone);CHKERRQ(ierr); 934 vids[0] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 935 ierr = DMPlexGetSupportSize(dm,he,&hesuppSize);CHKERRQ(ierr); 936 ierr = DMPlexGetSupport(dm,he,&hesupp);CHKERRQ(ierr); 937 for (f=0;f<hesuppSize;f++) { 938 PetscInt j; 939 940 ierr = DMPlexGetCone(dm,hesupp[f],&cone);CHKERRQ(ierr); 941 ierr = DMPlexGetConeSize(dm,hesupp[f],&coneSize);CHKERRQ(ierr); 942 for (j=0;j<coneSize;j++) { 943 PetscInt k; 944 for (k=0;k<hvsuppSize;k++) { 945 if (hvsupp[k] == cone[j]) { 946 skip[k] = PETSC_TRUE; 947 break; 948 } 949 } 950 } 951 } 952 for (i=0;i<hvsuppSize;i++) { 953 if (!skip[i]) { 954 ierr = DMPlexGetCone(dm,hvsupp[i],&cone);CHKERRQ(ierr); 955 vids[1] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 956 } 957 } 958 ierr = PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);CHKERRQ(ierr); 959 for (i=0;i<2;i++) { 960 ierr = PetscViewerASCIIPrintf(viewer," %D",(PetscInt)(vids[i]-vStart));CHKERRQ(ierr); 961 } 962 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 963 vp--; 964 } 965 } 966 break; 967 default: 968 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %D",ssize); 969 } 970 } 971 } 972 ierr = PetscFree(skip);CHKERRQ(ierr); 973 } 974 if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp); 975 } 976 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 977 ierr = PetscBTDestroy(&vown);CHKERRQ(ierr); 978 979 /* vertices */ 980 if (hovec) { /* higher-order meshes */ 981 const char *fec; 982 PetscInt i,n,s; 983 984 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 985 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);CHKERRQ(ierr); 986 ierr = PetscViewerASCIIPrintf(viewer,"nodes\n");CHKERRQ(ierr); 987 ierr = PetscObjectGetName((PetscObject)hovec,&fec);CHKERRQ(ierr); 988 ierr = PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");CHKERRQ(ierr); 989 ierr = PetscViewerASCIIPrintf(viewer,"%s\n",fec);CHKERRQ(ierr); 990 ierr = PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);CHKERRQ(ierr); 991 ierr = PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n");CHKERRQ(ierr); /*Ordering::byVDIM*/ 992 ierr = VecGetArrayRead(hovec,&array);CHKERRQ(ierr); 993 ierr = VecGetLocalSize(hovec,&n);CHKERRQ(ierr); 994 if (n%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Size of local coordinate vector %D incompatible with space dimension %D",n,sdim); 995 for (i=0;i<n/sdim;i++) { 996 for (s=0;s<sdim;s++) { 997 ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[i*sdim+s]));CHKERRQ(ierr); 998 } 999 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1000 } 1001 ierr = VecRestoreArrayRead(hovec,&array);CHKERRQ(ierr); 1002 } else { 1003 ierr = VecGetLocalSize(coordinates,&nvert);CHKERRQ(ierr); 1004 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 1005 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);CHKERRQ(ierr); 1006 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 1007 ierr = VecGetArrayRead(coordinates,&array);CHKERRQ(ierr); 1008 for (p=0;p<nvert/sdim;p++) { 1009 PetscInt s; 1010 for (s=0;s<sdim;s++) { 1011 PetscReal v = PetscRealPart(array[p*sdim+s]); 1012 1013 ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscIsInfOrNanReal(v) ? 0.0 : v);CHKERRQ(ierr); 1014 } 1015 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1016 } 1017 ierr = VecRestoreArrayRead(coordinates,&array);CHKERRQ(ierr); 1018 } 1019 PetscFunctionReturn(0); 1020 } 1021 1022 PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer) 1023 { 1024 PetscErrorCode ierr; 1025 PetscFunctionBegin; 1026 ierr = DMView_GLVis(dm,viewer,DMPlexView_GLVis_ASCII);CHKERRQ(ierr); 1027 PetscFunctionReturn(0); 1028 } 1029