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,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 = PetscObjectQuery((PetscObject)dm,"_glvis_plex_gnum",(PetscObject*)&globalNum);CHKERRQ(ierr); 65 if (!globalNum) { 66 ierr = DMPlexCreateCellNumbering_Internal(dm,PETSC_TRUE,&globalNum);CHKERRQ(ierr); 67 ierr = PetscObjectCompose((PetscObject)dm,"_glvis_plex_gnum",(PetscObject)globalNum);CHKERRQ(ierr); 68 ierr = PetscObjectDereference((PetscObject)globalNum);CHKERRQ(ierr); 69 } 70 ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr); 71 ierr = PetscBTCreate(vEnd-vStart,&vown);CHKERRQ(ierr); 72 for (c = cStart, totc = 0; c < cEnd; c++) { 73 if (gNum[c-cStart] >= 0) { 74 PetscInt i,numPoints,*points = NULL; 75 76 totc++; 77 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 78 for (i=0;i<numPoints*2;i+= 2) { 79 if ((points[i] >= vStart) && (points[i] < vEnd)) { 80 ierr = PetscBTSet(vown,points[i]-vStart);CHKERRQ(ierr); 81 } 82 } 83 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 84 } 85 } 86 for (f=0,Nv=0;f<vEnd-vStart;f++) if (PetscLikely(PetscBTLookup(vown,f))) Nv++; 87 88 ierr = DMCreateLocalVector(dm,&xlocal);CHKERRQ(ierr); 89 ierr = VecGetLocalSize(xlocal,&totdofs);CHKERRQ(ierr); 90 ierr = DMGetLocalSection(dm,&s);CHKERRQ(ierr); 91 ierr = PetscSectionGetNumFields(s,&nfields);CHKERRQ(ierr); 92 for (f=0,maxfields=0;f<nfields;f++) { 93 PetscInt bs; 94 95 ierr = PetscSectionGetFieldComponents(s,f,&bs);CHKERRQ(ierr); 96 maxfields += bs; 97 } 98 ierr = PetscCalloc7(maxfields,&fieldname,maxfields,&nlocal,maxfields,&bs,maxfields,&dims,maxfields,&fec_type,totdofs,&idxs,maxfields,&Ufield);CHKERRQ(ierr); 99 ierr = PetscNew(&ctx);CHKERRQ(ierr); 100 ierr = PetscCalloc1(maxfields,&ctx->scctx);CHKERRQ(ierr); 101 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 102 if (ds) { 103 for (f=0;f<nfields;f++) { 104 const char* fname; 105 char name[256]; 106 PetscObject disc; 107 size_t len; 108 109 ierr = PetscSectionGetFieldName(s,f,&fname);CHKERRQ(ierr); 110 ierr = PetscStrlen(fname,&len);CHKERRQ(ierr); 111 if (len) { 112 ierr = PetscStrcpy(name,fname);CHKERRQ(ierr); 113 } else { 114 ierr = PetscSNPrintf(name,256,"Field%D",f);CHKERRQ(ierr); 115 } 116 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 117 if (disc) { 118 PetscClassId id; 119 PetscInt Nc; 120 char fec[64]; 121 122 ierr = PetscObjectGetClassId(disc, &id);CHKERRQ(ierr); 123 if (id == PETSCFE_CLASSID) { 124 PetscFE fem = (PetscFE)disc; 125 PetscDualSpace sp; 126 PetscDualSpaceType spname; 127 PetscInt order; 128 PetscBool islag,continuous,H1 = PETSC_TRUE; 129 130 ierr = PetscFEGetNumComponents(fem,&Nc);CHKERRQ(ierr); 131 ierr = PetscFEGetDualSpace(fem,&sp);CHKERRQ(ierr); 132 ierr = PetscDualSpaceGetType(sp,&spname);CHKERRQ(ierr); 133 ierr = PetscStrcmp(spname,PETSCDUALSPACELAGRANGE,&islag);CHKERRQ(ierr); 134 if (!islag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dual space"); 135 ierr = PetscDualSpaceLagrangeGetContinuity(sp,&continuous);CHKERRQ(ierr); 136 ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr); 137 if (continuous && order > 0) { /* no support for high-order viz, still have to figure out the numbering */ 138 ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: H1_%DD_P1",dim);CHKERRQ(ierr); 139 } else { 140 if (!continuous && order) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Discontinuous space visualization currently unsupported for order %D",order); 141 H1 = PETSC_FALSE; 142 ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P%D",dim,order);CHKERRQ(ierr); 143 } 144 ierr = PetscStrallocpy(name,&fieldname[ctx->nf]);CHKERRQ(ierr); 145 bs[ctx->nf] = Nc; 146 dims[ctx->nf] = dim; 147 if (H1) { 148 nlocal[ctx->nf] = Nc * Nv; 149 ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr); 150 ierr = VecCreateSeq(PETSC_COMM_SELF,Nv*Nc,&xfield);CHKERRQ(ierr); 151 for (i=0,cum=0;i<vEnd-vStart;i++) { 152 PetscInt j,off; 153 154 if (PetscUnlikely(!PetscBTLookup(vown,i))) continue; 155 ierr = PetscSectionGetFieldOffset(s,i+vStart,f,&off);CHKERRQ(ierr); 156 for (j=0;j<Nc;j++) idxs[cum++] = off + j; 157 } 158 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),Nv*Nc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr); 159 } else { 160 nlocal[ctx->nf] = Nc * totc; 161 ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr); 162 ierr = VecCreateSeq(PETSC_COMM_SELF,Nc*totc,&xfield);CHKERRQ(ierr); 163 for (i=0,cum=0;i<cEnd-cStart;i++) { 164 PetscInt j,off; 165 166 if (PetscUnlikely(gNum[i] < 0)) continue; 167 ierr = PetscSectionGetFieldOffset(s,i+cStart,f,&off);CHKERRQ(ierr); 168 for (j=0;j<Nc;j++) idxs[cum++] = off + j; 169 } 170 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc*Nc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr); 171 } 172 ierr = VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);CHKERRQ(ierr); 173 ierr = VecDestroy(&xfield);CHKERRQ(ierr); 174 ierr = ISDestroy(&isfield);CHKERRQ(ierr); 175 ctx->nf++; 176 } else if (id == PETSCFV_CLASSID) { 177 PetscInt c; 178 179 ierr = PetscFVGetNumComponents((PetscFV)disc,&Nc);CHKERRQ(ierr); 180 ierr = PetscSNPrintf(fec,64,"FiniteElementCollection: L2_%DD_P0",dim);CHKERRQ(ierr); 181 for (c = 0; c < Nc; c++) { 182 char comp[256]; 183 ierr = PetscSNPrintf(comp,256,"%s-Comp%D",name,c);CHKERRQ(ierr); 184 ierr = PetscStrallocpy(comp,&fieldname[ctx->nf]);CHKERRQ(ierr); 185 bs[ctx->nf] = 1; /* Does PetscFV support components with different block size? */ 186 nlocal[ctx->nf] = totc; 187 dims[ctx->nf] = dim; 188 ierr = PetscStrallocpy(fec,&fec_type[ctx->nf]);CHKERRQ(ierr); 189 ierr = VecCreateSeq(PETSC_COMM_SELF,totc,&xfield);CHKERRQ(ierr); 190 for (i=0,cum=0;i<cEnd-cStart;i++) { 191 PetscInt off; 192 193 if (PetscUnlikely(gNum[i])<0) continue; 194 ierr = PetscSectionGetFieldOffset(s,i+cStart,f,&off);CHKERRQ(ierr); 195 idxs[cum++] = off + c; 196 } 197 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)xlocal),totc,idxs,PETSC_USE_POINTER,&isfield);CHKERRQ(ierr); 198 ierr = VecScatterCreate(xlocal,isfield,xfield,NULL,&ctx->scctx[ctx->nf]);CHKERRQ(ierr); 199 ierr = VecDestroy(&xfield);CHKERRQ(ierr); 200 ierr = ISDestroy(&isfield);CHKERRQ(ierr); 201 ctx->nf++; 202 } 203 } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",f); 204 } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing discretization for field %D",f); 205 } 206 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Needs a DS attached to the DM"); 207 ierr = PetscBTDestroy(&vown);CHKERRQ(ierr); 208 ierr = VecDestroy(&xlocal);CHKERRQ(ierr); 209 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 210 211 /* create work vectors */ 212 for (f=0;f<ctx->nf;f++) { 213 ierr = VecCreateMPI(PetscObjectComm((PetscObject)dm),nlocal[f],PETSC_DECIDE,&Ufield[f]);CHKERRQ(ierr); 214 ierr = PetscObjectSetName((PetscObject)Ufield[f],fieldname[f]);CHKERRQ(ierr); 215 ierr = VecSetBlockSize(Ufield[f],bs[f]);CHKERRQ(ierr); 216 ierr = VecSetDM(Ufield[f],dm);CHKERRQ(ierr); 217 } 218 219 /* customize the viewer */ 220 ierr = PetscViewerGLVisSetFields(viewer,ctx->nf,(const char**)fec_type,dims,DMPlexSampleGLVisFields_Private,(PetscObject*)Ufield,ctx,DestroyGLVisViewerCtx_Private);CHKERRQ(ierr); 221 for (f=0;f<ctx->nf;f++) { 222 ierr = PetscFree(fieldname[f]);CHKERRQ(ierr); 223 ierr = PetscFree(fec_type[f]);CHKERRQ(ierr); 224 ierr = VecDestroy(&Ufield[f]);CHKERRQ(ierr); 225 } 226 ierr = PetscFree7(fieldname,nlocal,bs,dims,fec_type,idxs,Ufield);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 typedef enum {MFEM_POINT=0,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE,MFEM_TETRAHEDRON,MFEM_CUBE,MFEM_PRISM,MFEM_UNDEF} MFEM_cid; 231 232 MFEM_cid mfem_table_cid[4][7] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF}, 233 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_UNDEF}, 234 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE ,MFEM_UNDEF,MFEM_UNDEF}, 235 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_TETRAHEDRON,MFEM_PRISM,MFEM_CUBE } }; 236 237 MFEM_cid mfem_table_cid_unint[4][9] = { {MFEM_POINT,MFEM_UNDEF,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF}, 238 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_UNDEF ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF}, 239 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_TRIANGLE,MFEM_SQUARE ,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_UNDEF}, 240 {MFEM_POINT,MFEM_UNDEF,MFEM_SEGMENT,MFEM_UNDEF ,MFEM_TETRAHEDRON,MFEM_UNDEF,MFEM_PRISM,MFEM_UNDEF,MFEM_CUBE } }; 241 242 static PetscErrorCode DMPlexGetPointMFEMCellID_Internal(DM dm, DMLabel label, PetscInt minl, PetscInt p, PetscInt *mid, PetscInt *cid) 243 { 244 DMLabel dlabel; 245 PetscInt depth,csize,pdepth,dim; 246 PetscErrorCode ierr; 247 248 PetscFunctionBegin; 249 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 250 ierr = DMLabelGetValue(dlabel,p,&pdepth);CHKERRQ(ierr); 251 ierr = DMPlexGetConeSize(dm,p,&csize);CHKERRQ(ierr); 252 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 253 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 254 if (label) { 255 ierr = DMLabelGetValue(label,p,mid);CHKERRQ(ierr); 256 *mid = *mid - minl + 1; /* MFEM does not like negative markers */ 257 } else *mid = 1; 258 if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */ 259 #if defined PETSC_USE_DEBUG 260 if (dim < 0 || dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D",dim); 261 if (csize > 8) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found cone size %D for point %D",csize,p); 262 if (depth != 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Found depth %D for point %D. You should interpolate the mesh first",depth,p); 263 #endif 264 *cid = mfem_table_cid_unint[dim][csize]; 265 } else { 266 #if defined PETSC_USE_DEBUG 267 if (csize > 6) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cone size %D for point %D",csize,p); 268 if (pdepth < 0 || pdepth > 3) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Depth %D for point %D",csize,p); 269 #endif 270 *cid = mfem_table_cid[pdepth][csize]; 271 } 272 PetscFunctionReturn(0); 273 } 274 275 static PetscErrorCode DMPlexGetPointMFEMVertexIDs_Internal(DM dm, PetscInt p, PetscSection csec, PetscInt *nv, int vids[]) 276 { 277 PetscInt dim,sdim,dof = 0,off = 0,i,q,vStart,vEnd,numPoints,*points = NULL; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 282 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 283 sdim = dim; 284 if (csec) { 285 PetscInt sStart,sEnd; 286 287 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 288 ierr = PetscSectionGetChart(csec,&sStart,&sEnd);CHKERRQ(ierr); 289 ierr = PetscSectionGetOffset(csec,vStart,&off);CHKERRQ(ierr); 290 off = off/sdim; 291 if (p >= sStart && p < sEnd) { 292 ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr); 293 } 294 } 295 if (!dof) { 296 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 297 for (i=0,q=0;i<numPoints*2;i+= 2) 298 if ((points[i] >= vStart) && (points[i] < vEnd)) 299 vids[q++] = (int)(points[i]-vStart+off); 300 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&numPoints,&points);CHKERRQ(ierr); 301 } else { 302 ierr = PetscSectionGetOffset(csec,p,&off);CHKERRQ(ierr); 303 ierr = PetscSectionGetDof(csec,p,&dof);CHKERRQ(ierr); 304 for (q=0;q<dof/sdim;q++) vids[q] = (int)(off/sdim + q); 305 } 306 *nv = q; 307 PetscFunctionReturn(0); 308 } 309 310 static PetscErrorCode DMPlexGlvisInvertHybrid(PetscInt cid,int vids[]) 311 { 312 int tmp; 313 314 PetscFunctionBegin; 315 if (cid == MFEM_SQUARE) { /* PETSc stores hybrid quads not as counter-clockwise quad */ 316 tmp = vids[2]; 317 vids[2] = vids[3]; 318 vids[3] = tmp; 319 } else if (cid == MFEM_PRISM) { /* MFEM uses a different orientation for the base and top triangles of the wedge */ 320 tmp = vids[1]; 321 vids[1] = vids[2]; 322 vids[2] = tmp; 323 tmp = vids[4]; 324 vids[4] = vids[5]; 325 vids[5] = tmp; 326 } 327 PetscFunctionReturn(0); 328 } 329 330 /* 331 ASCII visualization/dump: full support for simplices and tensor product cells. It supports AMR 332 Higher order meshes are also supported 333 */ 334 static PetscErrorCode DMPlexView_GLVis_ASCII(DM dm, PetscViewer viewer) 335 { 336 DMLabel label; 337 PetscSection coordSection,parentSection; 338 Vec coordinates,hovec; 339 const PetscScalar *array; 340 PetscInt bf,p,sdim,dim,depth,novl,minl; 341 PetscInt cStart,cEnd,cEndInterior,vStart,vEnd,nvert; 342 PetscMPIInt size; 343 PetscBool localized,isascii; 344 PetscBool enable_mfem,enable_boundary,enable_ncmesh,view_ovl = PETSC_FALSE; 345 PetscBT pown,vown; 346 PetscErrorCode ierr; 347 PetscContainer glvis_container; 348 PetscBool cellvertex = PETSC_FALSE, periodic, enabled = PETSC_TRUE; 349 PetscBool enable_emark,enable_bmark; 350 const char *fmt; 351 char emark[64] = "",bmark[64] = ""; 352 353 PetscFunctionBegin; 354 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 355 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 356 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 357 if (!isascii) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer must be of type VIEWERASCII"); 358 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&size);CHKERRQ(ierr); 359 if (size > 1) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Use single sequential viewers for parallel visualization"); 360 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 361 362 /* get container: determines if a process visualizes is portion of the data or not */ 363 ierr = PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);CHKERRQ(ierr); 364 if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing GLVis container"); 365 { 366 PetscViewerGLVisInfo glvis_info; 367 ierr = PetscContainerGetPointer(glvis_container,(void**)&glvis_info);CHKERRQ(ierr); 368 enabled = glvis_info->enabled; 369 fmt = glvis_info->fmt; 370 } 371 372 /* Users can attach a coordinate vector to the DM in case they have a higher-order mesh 373 DMPlex does not currently support HO meshes, so there's no API for this */ 374 ierr = PetscObjectQuery((PetscObject)dm,"_glvis_mesh_coords",(PetscObject*)&hovec);CHKERRQ(ierr); 375 376 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 377 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 378 cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior; 379 ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr); 380 ierr = DMGetPeriodicity(dm,&periodic,NULL,NULL,NULL);CHKERRQ(ierr); 381 ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr); 382 if (periodic && !localized && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Coordinates need to be localized"); 383 ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 384 ierr = DMGetCoordinateDim(dm,&sdim);CHKERRQ(ierr); 385 ierr = DMGetCoordinatesLocal(dm,&coordinates);CHKERRQ(ierr); 386 if (!coordinates && !hovec) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector"); 387 388 /* 389 a couple of sections of the mesh specification are disabled 390 - 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) 391 - vertex_parents: used for non-conforming meshes only when we want to use MFEM as a discretization package 392 and be able to derefine the mesh (MFEM does not currently have to ability to read ncmeshes in parallel) 393 */ 394 enable_boundary = PETSC_FALSE; 395 enable_ncmesh = PETSC_FALSE; 396 enable_mfem = PETSC_FALSE; 397 enable_emark = PETSC_FALSE; 398 enable_bmark = PETSC_FALSE; 399 /* I'm tired of problems with negative values in the markers, disable them */ 400 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"GLVis PetscViewer DMPlex Options","PetscViewer");CHKERRQ(ierr); 401 ierr = PetscOptionsBool("-viewer_glvis_dm_plex_enable_boundary","Enable boundary section in mesh representation",NULL,enable_boundary,&enable_boundary,NULL);CHKERRQ(ierr); 402 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); 403 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); 404 ierr = PetscOptionsBool("-viewer_glvis_dm_plex_overlap","Include overlap region in local meshes",NULL,view_ovl,&view_ovl,NULL);CHKERRQ(ierr); 405 ierr = PetscOptionsString("-viewer_glvis_dm_plex_emarker","String for the material id label",NULL,emark,emark,sizeof(emark),&enable_emark);CHKERRQ(ierr); 406 ierr = PetscOptionsString("-viewer_glvis_dm_plex_bmarker","String for the boundary id label",NULL,bmark,bmark,sizeof(bmark),&enable_bmark);CHKERRQ(ierr); 407 ierr = PetscOptionsEnd();CHKERRQ(ierr); 408 if (enable_bmark) enable_boundary = PETSC_TRUE; 409 410 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr); 411 if (enable_ncmesh && size > 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not supported in parallel"); 412 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 413 if (enable_boundary && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. " 414 "Alternatively, run with -viewer_glvis_dm_plex_enable_boundary 0"); 415 if (enable_ncmesh && depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated. " 416 "Alternatively, run with -viewer_glvis_dm_plex_enable_ncmesh 0"); 417 if (depth >=0 && dim != depth) { /* not interpolated, it assumes cell-vertex mesh */ 418 if (depth != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported depth %D. You should interpolate the mesh first",depth); 419 cellvertex = PETSC_TRUE; 420 } 421 422 /* Identify possible cells in the overlap */ 423 novl = 0; 424 pown = NULL; 425 if (size > 1) { 426 IS globalNum = NULL; 427 const PetscInt *gNum; 428 PetscBool ovl = PETSC_FALSE; 429 430 ierr = PetscObjectQuery((PetscObject)dm,"_glvis_plex_gnum",(PetscObject*)&globalNum);CHKERRQ(ierr); 431 if (!globalNum) { 432 if (view_ovl) { 433 ierr = ISCreateStride(PetscObjectComm((PetscObject)dm),cEnd-cStart,0,1,&globalNum);CHKERRQ(ierr); 434 } else { 435 ierr = DMPlexCreateCellNumbering_Internal(dm,PETSC_TRUE,&globalNum);CHKERRQ(ierr); 436 } 437 ierr = PetscObjectCompose((PetscObject)dm,"_glvis_plex_gnum",(PetscObject)globalNum);CHKERRQ(ierr); 438 ierr = PetscObjectDereference((PetscObject)globalNum);CHKERRQ(ierr); 439 } 440 ierr = ISGetIndices(globalNum,&gNum);CHKERRQ(ierr); 441 for (p=cStart; p<cEnd; p++) { 442 if (gNum[p-cStart] < 0) { 443 ovl = PETSC_TRUE; 444 novl++; 445 } 446 } 447 if (ovl) { 448 /* it may happen that pown get not destroyed, if the user closes the window while this function is running. 449 TODO: garbage collector? attach pown to dm? */ 450 ierr = PetscBTCreate(cEnd-cStart,&pown);CHKERRQ(ierr); 451 for (p=cStart; p<cEnd; p++) { 452 if (gNum[p-cStart] < 0) continue; 453 else { 454 ierr = PetscBTSet(pown,p-cStart);CHKERRQ(ierr); 455 } 456 } 457 } 458 ierr = ISRestoreIndices(globalNum,&gNum);CHKERRQ(ierr); 459 } 460 461 /* return if this process is disabled */ 462 if (!enabled) { 463 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 464 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 465 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 466 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 467 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 468 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 469 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 470 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 471 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 472 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 473 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 474 PetscFunctionReturn(0); 475 } 476 477 if (enable_mfem) { 478 if (periodic && !hovec) { /* we need to generate a vector of L2 coordinates, as this is how MFEM handles periodic meshes */ 479 PetscInt vpc = 0; 480 char fec[64]; 481 int vids[8] = {0,1,2,3,4,5,6,7}; 482 int hexv[8] = {0,1,3,2,4,5,7,6}, tetv[4] = {0,1,2,3}; 483 int quadv[8] = {0,1,3,2}, triv[3] = {0,1,2}; 484 int *dof = NULL; 485 PetscScalar *array,*ptr; 486 487 ierr = PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: L2_T1_%DD_P1",dim);CHKERRQ(ierr); 488 if (cEndInterior < cEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Support for hybrid meshed not currently implemented"); 489 if (cEnd-cStart) { 490 PetscInt fpc; 491 492 ierr = DMPlexGetConeSize(dm,cStart,&fpc);CHKERRQ(ierr); 493 switch(dim) { 494 case 1: 495 vpc = 2; 496 dof = hexv; 497 break; 498 case 2: 499 switch (fpc) { 500 case 3: 501 vpc = 3; 502 dof = triv; 503 break; 504 case 4: 505 vpc = 4; 506 dof = quadv; 507 break; 508 default: 509 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 510 break; 511 } 512 break; 513 case 3: 514 switch (fpc) { 515 case 4: /* TODO: still need to understand L2 ordering for tets */ 516 vpc = 4; 517 dof = tetv; 518 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled tethraedral case"); 519 break; 520 case 6: 521 if (cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: vertices per cell %D",fpc); 522 vpc = 8; 523 dof = hexv; 524 break; 525 case 8: 526 if (!cellvertex) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 527 vpc = 8; 528 dof = hexv; 529 break; 530 default: 531 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 532 break; 533 } 534 break; 535 default: 536 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim"); 537 break; 538 } 539 ierr = DMPlexInvertCell(dim,vpc,vids);CHKERRQ(ierr); 540 } 541 if (!dof) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Missing dofs"); 542 ierr = VecCreateSeq(PETSC_COMM_SELF,(cEnd-cStart-novl)*vpc*sdim,&hovec);CHKERRQ(ierr); 543 ierr = PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)hovec);CHKERRQ(ierr); 544 ierr = PetscObjectDereference((PetscObject)hovec);CHKERRQ(ierr); 545 ierr = PetscObjectSetName((PetscObject)hovec,fec);CHKERRQ(ierr); 546 ierr = VecGetArray(hovec,&array);CHKERRQ(ierr); 547 ptr = array; 548 for (p=cStart;p<cEnd;p++) { 549 PetscInt csize,v,d; 550 PetscScalar *vals = NULL; 551 552 if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue; 553 ierr = DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 554 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); 555 for (v=0;v<vpc;v++) { 556 for (d=0;d<sdim;d++) { 557 ptr[sdim*dof[v]+d] = vals[sdim*vids[v]+d]; 558 } 559 } 560 ptr += vpc*sdim; 561 ierr = DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 562 } 563 ierr = VecRestoreArray(hovec,&array);CHKERRQ(ierr); 564 } 565 } 566 /* if we have high-order coordinates in 3D, we need to specify the boundary */ 567 if (hovec && dim == 3) enable_boundary = PETSC_TRUE; 568 569 /* header */ 570 ierr = PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.1\n");CHKERRQ(ierr); 571 572 /* topological dimension */ 573 ierr = PetscViewerASCIIPrintf(viewer,"\ndimension\n");CHKERRQ(ierr); 574 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",dim);CHKERRQ(ierr); 575 576 /* elements */ 577 minl = 1; 578 label = NULL; 579 if (enable_emark) { 580 PetscInt lminl = PETSC_MAX_INT; 581 582 ierr = DMGetLabel(dm,emark,&label);CHKERRQ(ierr); 583 if (label) { 584 IS vals; 585 PetscInt ldef; 586 587 ierr = DMLabelGetDefaultValue(label,&ldef);CHKERRQ(ierr); 588 ierr = DMLabelGetValueIS(label,&vals);CHKERRQ(ierr); 589 ierr = ISGetMinMax(vals,&lminl,NULL);CHKERRQ(ierr); 590 ierr = ISDestroy(&vals);CHKERRQ(ierr); 591 lminl = PetscMin(ldef,lminl); 592 } 593 ierr = MPIU_Allreduce(&lminl,&minl,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 594 if (minl == PETSC_MAX_INT) minl = 1; 595 } 596 ierr = PetscViewerASCIIPrintf(viewer,"\nelements\n");CHKERRQ(ierr); 597 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",cEnd-cStart-novl);CHKERRQ(ierr); 598 for (p=cStart;p<cEnd;p++) { 599 int vids[8]; 600 PetscInt i,nv = 0,cid = -1,mid = 1; 601 602 if (PetscUnlikely(pown && !PetscBTLookup(pown,p-cStart))) continue; 603 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,minl,p,&mid,&cid);CHKERRQ(ierr); 604 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,p,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 605 ierr = DMPlexInvertCell(dim,nv,vids);CHKERRQ(ierr); 606 if (p >= cEndInterior) { 607 ierr = DMPlexGlvisInvertHybrid(cid,vids);CHKERRQ(ierr); 608 } 609 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 610 for (i=0;i<nv;i++) { 611 ierr = PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);CHKERRQ(ierr); 612 } 613 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 614 } 615 616 /* boundary */ 617 ierr = PetscViewerASCIIPrintf(viewer,"\nboundary\n");CHKERRQ(ierr); 618 if (!enable_boundary) { 619 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",0);CHKERRQ(ierr); 620 } else { 621 DMLabel perLabel; 622 PetscBT bfaces; 623 PetscInt fStart,fEnd,*fcells; 624 PetscInt *faces = NULL,fpc = 0,vpf = 0, vpc = 0; 625 PetscInt *facesH = NULL,fpcH = 0,vpfH = 0, vpcH = 0; 626 PetscInt fv1[] = {0,1}, 627 fv2tri[] = {0,1, 628 1,2, 629 2,0}, 630 fv2quad[] = {0,1, 631 1,2, 632 2,3, 633 3,0}, 634 fv2quadH[] = {0,1, 635 2,3, 636 0,2, 637 1,3}, 638 fv3tet[] = {0,1,2, 639 0,3,1, 640 0,2,3, 641 2,1,3}, 642 fv3wedge[] = {0,1,2,-1, 643 3,4,5,-1, 644 0,1,3,4, 645 1,2,4,5, 646 2,0,5,3}, 647 fv3hex[] = {0,1,2,3, 648 4,5,6,7, 649 0,3,5,4, 650 2,1,7,6, 651 3,2,6,5, 652 0,4,7,1}; 653 654 /* determine orientation of boundary mesh */ 655 if (cEnd-cStart) { 656 if (cEndInterior < cEnd) { 657 ierr = DMPlexGetConeSize(dm,cEndInterior,&fpcH);CHKERRQ(ierr); 658 } 659 if (cEndInterior > cStart) { 660 ierr = DMPlexGetConeSize(dm,cStart,&fpc);CHKERRQ(ierr); 661 } 662 switch(dim) { 663 case 1: 664 if (fpc != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case faces per cell %D",fpc); 665 faces = fv1; 666 vpf = 1; 667 vpc = 2; 668 break; 669 case 2: 670 switch (fpc) { 671 case 0: 672 case 3: 673 faces = fv2tri; 674 vpf = 2; 675 vpc = 3; 676 if (fpcH == 4) { 677 facesH = fv2quadH; 678 vpfH = 2; 679 vpcH = 4; 680 } else if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH); 681 break; 682 case 4: 683 faces = fv2quad; 684 vpf = 2; 685 vpc = 4; 686 if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH); 687 break; 688 default: 689 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 690 break; 691 } 692 break; 693 case 3: 694 switch (fpc) { 695 case 0: 696 case 4: 697 faces = fv3tet; 698 vpf = 3; 699 vpc = 4; 700 if (fpcH == 5) { 701 facesH = fv3wedge; 702 vpfH = -4; 703 vpcH = 6; 704 } else if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH); 705 break; 706 case 6: 707 faces = fv3hex; 708 vpf = 4; 709 vpc = 8; 710 if (fpcH) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case hybrid faces per cell %D",fpcH); 711 break; 712 default: 713 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unhandled case: faces per cell %D",fpc); 714 break; 715 } 716 break; 717 default: 718 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dim"); 719 break; 720 } 721 } 722 ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);CHKERRQ(ierr); 723 ierr = PetscBTCreate(fEnd-fStart,&bfaces);CHKERRQ(ierr); 724 ierr = DMPlexGetMaxSizes(dm,NULL,&p);CHKERRQ(ierr); 725 ierr = PetscMalloc1(p,&fcells);CHKERRQ(ierr); 726 ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr); 727 if (!perLabel && localized) { /* this periodic cut can be moved up to DMPlex setup */ 728 ierr = DMCreateLabel(dm,"glvis_periodic_cut");CHKERRQ(ierr); 729 ierr = DMGetLabel(dm,"glvis_periodic_cut",&perLabel);CHKERRQ(ierr); 730 ierr = DMLabelSetDefaultValue(perLabel,1);CHKERRQ(ierr); 731 for (p=cStart;p<cEnd;p++) { 732 PetscInt dof, uvpc; 733 734 ierr = PetscSectionGetDof(coordSection,p,&dof);CHKERRQ(ierr); 735 if (dof) { 736 PetscInt v,csize,cellClosureSize,*cellClosure = NULL,*vidxs = NULL; 737 PetscScalar *vals = NULL; 738 if (p < cEndInterior) uvpc = vpc; 739 else uvpc = vpcH; 740 if (dof%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Incompatible number of cell dofs %D and space dimension %D",dof,sdim); 741 if (dof/sdim != uvpc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Incompatible number of cell dofs %D, vertices %D and space dim %D",dof/sdim,uvpc,sdim); 742 ierr = DMPlexVecGetClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 743 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 744 for (v=0;v<cellClosureSize;v++) 745 if (cellClosure[2*v] >= vStart && cellClosure[2*v] < vEnd) { 746 vidxs = cellClosure + 2*v; 747 break; 748 } 749 if (!vidxs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing vertices"); 750 for (v=0;v<uvpc;v++) { 751 PetscInt s; 752 753 for (s=0;s<sdim;s++) { 754 if (PetscAbsScalar(vals[v*sdim+s]-vals[v*sdim+s+uvpc*sdim])>PETSC_MACHINE_EPSILON) { 755 ierr = DMLabelSetValue(perLabel,vidxs[2*v],2);CHKERRQ(ierr); 756 } 757 } 758 } 759 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&cellClosureSize,&cellClosure);CHKERRQ(ierr); 760 ierr = DMPlexVecRestoreClosure(dm,coordSection,coordinates,p,&csize,&vals);CHKERRQ(ierr); 761 } 762 } 763 if (dim > 1) { 764 PetscInt eEnd,eStart; 765 766 ierr = DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);CHKERRQ(ierr); 767 for (p=eStart;p<eEnd;p++) { 768 const PetscInt *cone; 769 PetscInt coneSize,i; 770 PetscBool ispe = PETSC_TRUE; 771 772 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 773 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 774 for (i=0;i<coneSize;i++) { 775 PetscInt v; 776 777 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 778 ispe = (PetscBool)(ispe && (v==2)); 779 } 780 if (ispe && coneSize) { 781 PetscInt ch, numChildren; 782 const PetscInt *children; 783 784 ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr); 785 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 786 for (ch = 0; ch < numChildren; ch++) { 787 ierr = DMLabelSetValue(perLabel,children[ch],2);CHKERRQ(ierr); 788 } 789 } 790 } 791 if (dim > 2) { 792 for (p=fStart;p<fEnd;p++) { 793 const PetscInt *cone; 794 PetscInt coneSize,i; 795 PetscBool ispe = PETSC_TRUE; 796 797 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 798 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 799 for (i=0;i<coneSize;i++) { 800 PetscInt v; 801 802 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 803 ispe = (PetscBool)(ispe && (v==2)); 804 } 805 if (ispe && coneSize) { 806 PetscInt ch, numChildren; 807 const PetscInt *children; 808 809 ierr = DMLabelSetValue(perLabel,p,2);CHKERRQ(ierr); 810 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 811 for (ch = 0; ch < numChildren; ch++) { 812 ierr = DMLabelSetValue(perLabel,children[ch],2);CHKERRQ(ierr); 813 } 814 } 815 } 816 } 817 } 818 } 819 for (p=fStart;p<fEnd;p++) { 820 const PetscInt *support; 821 PetscInt supportSize; 822 PetscBool isbf = PETSC_FALSE; 823 824 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 825 if (pown) { 826 PetscBool has_owned = PETSC_FALSE, has_ghost = PETSC_FALSE; 827 PetscInt i; 828 829 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 830 for (i=0;i<supportSize;i++) { 831 if (PetscLikely(PetscBTLookup(pown,support[i]-cStart))) has_owned = PETSC_TRUE; 832 else has_ghost = PETSC_TRUE; 833 } 834 isbf = (PetscBool)((supportSize == 1 && has_owned) || (supportSize > 1 && has_owned && has_ghost)); 835 } else { 836 isbf = (PetscBool)(supportSize == 1); 837 } 838 if (!isbf && perLabel) { 839 const PetscInt *cone; 840 PetscInt coneSize,i; 841 842 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 843 ierr = DMPlexGetConeSize(dm,p,&coneSize);CHKERRQ(ierr); 844 isbf = PETSC_TRUE; 845 for (i=0;i<coneSize;i++) { 846 PetscInt v,d; 847 848 ierr = DMLabelGetValue(perLabel,cone[i],&v);CHKERRQ(ierr); 849 ierr = DMLabelGetDefaultValue(perLabel,&d);CHKERRQ(ierr); 850 isbf = (PetscBool)(isbf && v != d); 851 } 852 } 853 if (isbf) { 854 ierr = PetscBTSet(bfaces,p-fStart);CHKERRQ(ierr); 855 } 856 } 857 /* count boundary faces */ 858 for (p=fStart,bf=0;p<fEnd;p++) { 859 if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) { 860 const PetscInt *support; 861 PetscInt supportSize,c; 862 863 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 864 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 865 for (c=0;c<supportSize;c++) { 866 const PetscInt *cone; 867 PetscInt cell,cl,coneSize; 868 869 cell = support[c]; 870 if (pown && PetscUnlikely(!PetscBTLookup(pown,cell-cStart))) continue; 871 ierr = DMPlexGetCone(dm,cell,&cone);CHKERRQ(ierr); 872 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 873 for (cl=0;cl<coneSize;cl++) { 874 if (cone[cl] == p) { 875 bf += 1; 876 break; 877 } 878 } 879 } 880 } 881 } 882 minl = 1; 883 label = NULL; 884 if (enable_bmark) { 885 PetscInt lminl = PETSC_MAX_INT; 886 887 ierr = DMGetLabel(dm,bmark,&label);CHKERRQ(ierr); 888 if (label) { 889 IS vals; 890 PetscInt ldef; 891 892 ierr = DMLabelGetDefaultValue(label,&ldef);CHKERRQ(ierr); 893 ierr = DMLabelGetValueIS(label,&vals);CHKERRQ(ierr); 894 ierr = ISGetMinMax(vals,&lminl,NULL);CHKERRQ(ierr); 895 ierr = ISDestroy(&vals);CHKERRQ(ierr); 896 lminl = PetscMin(ldef,lminl); 897 } 898 ierr = MPIU_Allreduce(&lminl,&minl,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 899 if (minl == PETSC_MAX_INT) minl = 1; 900 } 901 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",bf);CHKERRQ(ierr); 902 for (p=fStart;p<fEnd;p++) { 903 if (PetscUnlikely(PetscBTLookup(bfaces,p-fStart))) { 904 const PetscInt *support; 905 PetscInt supportSize,c,nc = 0; 906 907 ierr = DMPlexGetSupportSize(dm,p,&supportSize);CHKERRQ(ierr); 908 ierr = DMPlexGetSupport(dm,p,&support);CHKERRQ(ierr); 909 if (pown) { 910 for (c=0;c<supportSize;c++) { 911 if (PetscLikely(PetscBTLookup(pown,support[c]-cStart))) { 912 fcells[nc++] = support[c]; 913 } 914 } 915 } else for (c=0;c<supportSize;c++) fcells[nc++] = support[c]; 916 for (c=0;c<nc;c++) { 917 const PetscInt *cone; 918 int vids[8]; 919 PetscInt i,coneSize,cell,cl,nv,cid = -1,mid = -1; 920 921 cell = fcells[c]; 922 ierr = DMPlexGetCone(dm,cell,&cone);CHKERRQ(ierr); 923 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 924 for (cl=0;cl<coneSize;cl++) 925 if (cone[cl] == p) 926 break; 927 if (cl == coneSize) continue; 928 929 /* face material id and type */ 930 ierr = DMPlexGetPointMFEMCellID_Internal(dm,label,minl,p,&mid,&cid);CHKERRQ(ierr); 931 ierr = PetscViewerASCIIPrintf(viewer,"%D %D",mid,cid);CHKERRQ(ierr); 932 /* vertex ids */ 933 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,cell,(localized && !hovec) ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 934 if (cell >= cEndInterior) { 935 PetscInt nv = vpfH, inc = vpfH; 936 if (vpfH < 0) { /* Wedge */ 937 if (cl == 0 || cl == 1) nv = 3; 938 else nv = 4; 939 inc = -vpfH; 940 } 941 for (i=0;i<nv;i++) { 942 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[facesH[cl*inc+i]]);CHKERRQ(ierr); 943 } 944 } else { 945 for (i=0;i<vpf;i++) { 946 ierr = PetscViewerASCIIPrintf(viewer," %d",vids[faces[cl*vpf+i]]);CHKERRQ(ierr); 947 } 948 } 949 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 950 bf -= 1; 951 } 952 } 953 } 954 if (bf) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Remaining boundary faces %D",bf); 955 ierr = PetscBTDestroy(&bfaces);CHKERRQ(ierr); 956 ierr = PetscFree(fcells);CHKERRQ(ierr); 957 } 958 959 /* mark owned vertices */ 960 vown = NULL; 961 if (pown) { 962 ierr = PetscBTCreate(vEnd-vStart,&vown);CHKERRQ(ierr); 963 for (p=cStart;p<cEnd;p++) { 964 PetscInt i,closureSize,*closure = NULL; 965 966 if (PetscUnlikely(!PetscBTLookup(pown,p-cStart))) continue; 967 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 968 for (i=0;i<closureSize;i++) { 969 const PetscInt pp = closure[2*i]; 970 971 if (pp >= vStart && pp < vEnd) { 972 ierr = PetscBTSet(vown,pp-vStart);CHKERRQ(ierr); 973 } 974 } 975 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 976 } 977 } 978 979 /* vertex_parents (Non-conforming meshes) */ 980 parentSection = NULL; 981 if (enable_ncmesh) { 982 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 983 } 984 if (parentSection) { 985 PetscInt vp,gvp; 986 987 for (vp=0,p=vStart;p<vEnd;p++) { 988 DMLabel dlabel; 989 PetscInt parent,depth; 990 991 if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue; 992 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 993 ierr = DMLabelGetValue(dlabel,p,&depth);CHKERRQ(ierr); 994 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 995 if (parent != p) vp++; 996 } 997 ierr = MPIU_Allreduce(&vp,&gvp,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 998 if (gvp) { 999 PetscInt maxsupp; 1000 PetscBool *skip = NULL; 1001 1002 ierr = PetscViewerASCIIPrintf(viewer,"\nvertex_parents\n");CHKERRQ(ierr); 1003 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vp);CHKERRQ(ierr); 1004 ierr = DMPlexGetMaxSizes(dm,NULL,&maxsupp);CHKERRQ(ierr); 1005 ierr = PetscMalloc1(maxsupp,&skip);CHKERRQ(ierr); 1006 for (p=vStart;p<vEnd;p++) { 1007 DMLabel dlabel; 1008 PetscInt parent; 1009 1010 if (PetscUnlikely(vown && !PetscBTLookup(vown,p-vStart))) continue; 1011 ierr = DMPlexGetDepthLabel(dm,&dlabel);CHKERRQ(ierr); 1012 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 1013 if (parent != p) { 1014 int vids[8] = { -1, -1, -1, -1, -1, -1, -1, -1 }; /* silent overzealous clang static analyzer */ 1015 PetscInt i,nv,ssize,n,numChildren,depth = -1; 1016 const PetscInt *children; 1017 1018 ierr = DMPlexGetConeSize(dm,parent,&ssize);CHKERRQ(ierr); 1019 switch (ssize) { 1020 case 2: /* edge */ 1021 nv = 0; 1022 ierr = DMPlexGetPointMFEMVertexIDs_Internal(dm,parent,localized ? coordSection : NULL,&nv,vids);CHKERRQ(ierr); 1023 ierr = DMPlexInvertCell(dim,nv,vids);CHKERRQ(ierr); 1024 ierr = PetscViewerASCIIPrintf(viewer,"%D",p-vStart);CHKERRQ(ierr); 1025 for (i=0;i<nv;i++) { 1026 ierr = PetscViewerASCIIPrintf(viewer," %D",(PetscInt)vids[i]);CHKERRQ(ierr); 1027 } 1028 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1029 vp--; 1030 break; 1031 case 4: /* face */ 1032 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 1033 for (n=0;n<numChildren;n++) { 1034 ierr = DMLabelGetValue(dlabel,children[n],&depth);CHKERRQ(ierr); 1035 if (!depth) { 1036 const PetscInt *hvsupp,*hesupp,*cone; 1037 PetscInt hvsuppSize,hesuppSize,coneSize; 1038 PetscInt hv = children[n],he = -1,f; 1039 1040 ierr = PetscArrayzero(skip,maxsupp);CHKERRQ(ierr); 1041 ierr = DMPlexGetSupportSize(dm,hv,&hvsuppSize);CHKERRQ(ierr); 1042 ierr = DMPlexGetSupport(dm,hv,&hvsupp);CHKERRQ(ierr); 1043 for (i=0;i<hvsuppSize;i++) { 1044 PetscInt ep; 1045 ierr = DMPlexGetTreeParent(dm,hvsupp[i],&ep,NULL);CHKERRQ(ierr); 1046 if (ep != hvsupp[i]) { 1047 he = hvsupp[i]; 1048 } else { 1049 skip[i] = PETSC_TRUE; 1050 } 1051 } 1052 if (he == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertex %D support size %D: hanging edge not found",hv,hvsuppSize); 1053 ierr = DMPlexGetCone(dm,he,&cone);CHKERRQ(ierr); 1054 vids[0] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 1055 ierr = DMPlexGetSupportSize(dm,he,&hesuppSize);CHKERRQ(ierr); 1056 ierr = DMPlexGetSupport(dm,he,&hesupp);CHKERRQ(ierr); 1057 for (f=0;f<hesuppSize;f++) { 1058 PetscInt j; 1059 1060 ierr = DMPlexGetCone(dm,hesupp[f],&cone);CHKERRQ(ierr); 1061 ierr = DMPlexGetConeSize(dm,hesupp[f],&coneSize);CHKERRQ(ierr); 1062 for (j=0;j<coneSize;j++) { 1063 PetscInt k; 1064 for (k=0;k<hvsuppSize;k++) { 1065 if (hvsupp[k] == cone[j]) { 1066 skip[k] = PETSC_TRUE; 1067 break; 1068 } 1069 } 1070 } 1071 } 1072 for (i=0;i<hvsuppSize;i++) { 1073 if (!skip[i]) { 1074 ierr = DMPlexGetCone(dm,hvsupp[i],&cone);CHKERRQ(ierr); 1075 vids[1] = (int)((cone[0] == hv) ? cone[1] : cone[0]); 1076 } 1077 } 1078 ierr = PetscViewerASCIIPrintf(viewer,"%D",hv-vStart);CHKERRQ(ierr); 1079 for (i=0;i<2;i++) { 1080 ierr = PetscViewerASCIIPrintf(viewer," %D",(PetscInt)(vids[i]-vStart));CHKERRQ(ierr); 1081 } 1082 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1083 vp--; 1084 } 1085 } 1086 break; 1087 default: 1088 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Don't know how to deal with support size %D",ssize); 1089 } 1090 } 1091 } 1092 ierr = PetscFree(skip);CHKERRQ(ierr); 1093 } 1094 if (vp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D hanging vertices",vp); 1095 } 1096 ierr = PetscBTDestroy(&pown);CHKERRQ(ierr); 1097 ierr = PetscBTDestroy(&vown);CHKERRQ(ierr); 1098 1099 /* vertices */ 1100 if (hovec) { /* higher-order meshes */ 1101 const char *fec; 1102 PetscInt i,n,s; 1103 1104 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 1105 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",vEnd-vStart);CHKERRQ(ierr); 1106 ierr = PetscViewerASCIIPrintf(viewer,"nodes\n");CHKERRQ(ierr); 1107 ierr = PetscObjectGetName((PetscObject)hovec,&fec);CHKERRQ(ierr); 1108 ierr = PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");CHKERRQ(ierr); 1109 ierr = PetscViewerASCIIPrintf(viewer,"%s\n",fec);CHKERRQ(ierr); 1110 ierr = PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);CHKERRQ(ierr); 1111 ierr = PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n");CHKERRQ(ierr); /*Ordering::byVDIM*/ 1112 ierr = VecGetArrayRead(hovec,&array);CHKERRQ(ierr); 1113 ierr = VecGetLocalSize(hovec,&n);CHKERRQ(ierr); 1114 if (n%sdim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Size of local coordinate vector %D incompatible with space dimension %D",n,sdim); 1115 for (i=0;i<n/sdim;i++) { 1116 for (s=0;s<sdim;s++) { 1117 ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[i*sdim+s]));CHKERRQ(ierr); 1118 } 1119 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1120 } 1121 ierr = VecRestoreArrayRead(hovec,&array);CHKERRQ(ierr); 1122 } else { 1123 ierr = VecGetLocalSize(coordinates,&nvert);CHKERRQ(ierr); 1124 ierr = PetscViewerASCIIPrintf(viewer,"\nvertices\n");CHKERRQ(ierr); 1125 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",nvert/sdim);CHKERRQ(ierr); 1126 ierr = PetscViewerASCIIPrintf(viewer,"%D\n",sdim);CHKERRQ(ierr); 1127 ierr = VecGetArrayRead(coordinates,&array);CHKERRQ(ierr); 1128 for (p=0;p<nvert/sdim;p++) { 1129 PetscInt s; 1130 for (s=0;s<sdim;s++) { 1131 PetscReal v = PetscRealPart(array[p*sdim+s]); 1132 1133 ierr = PetscViewerASCIIPrintf(viewer,fmt,PetscIsInfOrNanReal(v) ? 0.0 : v);CHKERRQ(ierr); 1134 } 1135 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1136 } 1137 ierr = VecRestoreArrayRead(coordinates,&array);CHKERRQ(ierr); 1138 } 1139 PetscFunctionReturn(0); 1140 } 1141 1142 PetscErrorCode DMPlexView_GLVis(DM dm, PetscViewer viewer) 1143 { 1144 PetscErrorCode ierr; 1145 PetscFunctionBegin; 1146 ierr = DMView_GLVis(dm,viewer,DMPlexView_GLVis_ASCII);CHKERRQ(ierr); 1147 PetscFunctionReturn(0); 1148 } 1149