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