Actual source code: grglvis.c
1: /* Routines to visualize DMDAs and fields through GLVis */
3: #include <petsc/private/dmdaimpl.h>
4: #include <petsc/private/glvisviewerimpl.h>
6: typedef struct {
7: PetscBool ll;
8: } DMDAGhostedGLVisViewerCtx;
10: static PetscErrorCode DMDAGhostedDestroyGLVisViewerCtx_Private(void **vctx)
11: {
12: PetscFree(*vctx);
13: return 0;
14: }
16: typedef struct {
17: Vec xlocal;
18: } DMDAFieldGLVisViewerCtx;
20: static PetscErrorCode DMDAFieldDestroyGLVisViewerCtx_Private(void *vctx)
21: {
22: DMDAFieldGLVisViewerCtx *ctx = (DMDAFieldGLVisViewerCtx*)vctx;
24: VecDestroy(&ctx->xlocal);
25: PetscFree(vctx);
26: return 0;
27: }
29: /*
30: dactx->ll is false -> all but the last proc per dimension claim the ghosted node on the right
31: dactx->ll is true -> all but the first proc per dimension claim the ghosted node on the left
32: */
33: static PetscErrorCode DMDAGetNumElementsGhosted(DM da, PetscInt *nex, PetscInt *ney, PetscInt *nez)
34: {
35: DMDAGhostedGLVisViewerCtx *dactx;
36: PetscInt sx,sy,sz,ien,jen,ken;
38: /* Appease -Wmaybe-uninitialized */
39: if (nex) *nex = -1;
40: if (ney) *ney = -1;
41: if (nez) *nez = -1;
42: DMDAGetCorners(da,&sx,&sy,&sz,&ien,&jen,&ken);
43: DMGetApplicationContext(da,&dactx);
44: if (dactx->ll) {
45: PetscInt dim;
47: DMGetDimension(da,&dim);
48: if (!sx) ien--;
49: if (!sy && dim > 1) jen--;
50: if (!sz && dim > 2) ken--;
51: } else {
52: PetscInt M,N,P;
54: DMDAGetInfo(da,NULL,&M,&N,&P,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
55: if (sx+ien == M) ien--;
56: if (sy+jen == N) jen--;
57: if (sz+ken == P) ken--;
58: }
59: if (nex) *nex = ien;
60: if (ney) *ney = jen;
61: if (nez) *nez = ken;
62: return 0;
63: }
65: /* inherits number of vertices from DMDAGetNumElementsGhosted */
66: static PetscErrorCode DMDAGetNumVerticesGhosted(DM da, PetscInt *nvx, PetscInt *nvy, PetscInt *nvz)
67: {
68: PetscInt ien = 0,jen = 0,ken = 0,dim;
69: PetscInt tote;
71: DMGetDimension(da,&dim);
72: DMDAGetNumElementsGhosted(da,&ien,&jen,&ken);
73: tote = ien * (dim > 1 ? jen : 1) * (dim > 2 ? ken : 1);
74: if (tote) {
75: ien = ien+1;
76: jen = dim > 1 ? jen+1 : 1;
77: ken = dim > 2 ? ken+1 : 1;
78: }
79: if (nvx) *nvx = ien;
80: if (nvy) *nvy = jen;
81: if (nvz) *nvz = ken;
82: return 0;
83: }
85: static PetscErrorCode DMDASampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXf[], void *vctx)
86: {
87: DM da;
88: DMDAFieldGLVisViewerCtx *ctx = (DMDAFieldGLVisViewerCtx*)vctx;
89: DMDAGhostedGLVisViewerCtx *dactx;
90: const PetscScalar *array;
91: PetscScalar **arrayf;
92: PetscInt i,f,ii,ien,jen,ken,ie,je,ke,bs,*bss;
93: PetscInt sx,sy,sz,gsx,gsy,gsz,ist,jst,kst,gm,gn,gp;
95: VecGetDM(ctx->xlocal,&da);
97: DMGetApplicationContext(da,&dactx);
98: VecGetBlockSize(ctx->xlocal,&bs);
99: DMGlobalToLocalBegin(da,(Vec)oX,INSERT_VALUES,ctx->xlocal);
100: DMGlobalToLocalEnd(da,(Vec)oX,INSERT_VALUES,ctx->xlocal);
101: DMDAGetNumVerticesGhosted(da,&ien,&jen,&ken);
102: DMDAGetGhostCorners(da,&gsx,&gsy,&gsz,&gm,&gn,&gp);
103: DMDAGetCorners(da,&sx,&sy,&sz,NULL,NULL,NULL);
104: if (dactx->ll) {
105: kst = jst = ist = 0;
106: } else {
107: kst = gsz != sz ? 1 : 0;
108: jst = gsy != sy ? 1 : 0;
109: ist = gsx != sx ? 1 : 0;
110: }
111: PetscMalloc2(nf,&arrayf,nf,&bss);
112: VecGetArrayRead(ctx->xlocal,&array);
113: for (f=0;f<nf;f++) {
114: VecGetBlockSize((Vec)oXf[f],&bss[f]);
115: VecGetArray((Vec)oXf[f],&arrayf[f]);
116: }
117: for (ke = kst, ii = 0; ke < kst + ken; ke++) {
118: for (je = jst; je < jst + jen; je++) {
119: for (ie = ist; ie < ist + ien; ie++) {
120: PetscInt cf,b;
121: i = ke * gm * gn + je * gm + ie;
122: for (f=0,cf=0;f<nf;f++)
123: for (b=0;b<bss[f];b++)
124: arrayf[f][bss[f]*ii+b] = array[i*bs+cf++];
125: ii++;
126: }
127: }
128: }
129: for (f=0;f<nf;f++) VecRestoreArray((Vec)oXf[f],&arrayf[f]);
130: VecRestoreArrayRead(ctx->xlocal,&array);
131: PetscFree2(arrayf,bss);
132: return 0;
133: }
135: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_DMDA(PetscObject oda, PetscViewer viewer)
136: {
137: DM da = (DM)oda,daview;
140: PetscObjectQuery(oda,"GLVisGraphicsDMDAGhosted",(PetscObject*)&daview);
141: if (!daview) {
142: DMDAGhostedGLVisViewerCtx *dactx;
143: DM dacoord = NULL;
144: Vec xcoor,xcoorl;
145: PetscBool hashocoord = PETSC_FALSE;
146: const PetscInt *lx,*ly,*lz;
147: PetscInt dim,M,N,P,m,n,p,dof,s,i;
149: PetscNew(&dactx);
150: dactx->ll = PETSC_TRUE; /* default to match elements layout obtained by DMDAGetElements */
151: PetscOptionsBegin(PetscObjectComm(oda),oda->prefix,"GLVis PetscViewer DMDA Options","PetscViewer");
152: PetscOptionsBool("-viewer_glvis_dm_da_ll","Left-looking subdomain view",NULL,dactx->ll,&dactx->ll,NULL);
153: PetscOptionsEnd();
154: /* Create a properly ghosted DMDA to visualize the mesh and the fields associated with */
155: DMDAGetInfo(da,&dim,&M,&N,&P,&m,&n,&p,&dof,&s,NULL,NULL,NULL,NULL);
156: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
157: PetscObjectQuery((PetscObject)da,"_glvis_mesh_coords",(PetscObject*)&xcoor);
158: if (!xcoor) {
159: DMGetCoordinates(da,&xcoor);
160: } else {
161: hashocoord = PETSC_TRUE;
162: }
163: PetscInfo(da,"Creating auxilary DMDA for managing GLVis graphics\n");
164: switch (dim) {
165: case 1:
166: DMDACreate1d(PetscObjectComm((PetscObject)da),DM_BOUNDARY_NONE,M,dof,1,lx,&daview);
167: if (!hashocoord) {
168: DMDACreate1d(PetscObjectComm((PetscObject)da),DM_BOUNDARY_NONE,M,1,1,lx,&dacoord);
169: }
170: break;
171: case 2:
172: DMDACreate2d(PetscObjectComm((PetscObject)da),DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,m,n,dof,1,lx,ly,&daview);
173: if (!hashocoord) {
174: DMDACreate2d(PetscObjectComm((PetscObject)da),DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,m,n,2,1,lx,ly,&dacoord);
175: }
176: break;
177: case 3:
178: DMDACreate3d(PetscObjectComm((PetscObject)da),DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,P,m,n,p,dof,1,lx,ly,lz,&daview);
179: if (!hashocoord) {
180: DMDACreate3d(PetscObjectComm((PetscObject)da),DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,P,m,n,p,3,1,lx,ly,lz,&dacoord);
181: }
182: break;
183: default:
184: SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Unsupported dimension %D",dim);
185: }
186: DMSetApplicationContext(daview,dactx);
187: DMSetApplicationContextDestroy(daview,DMDAGhostedDestroyGLVisViewerCtx_Private);
188: DMSetUp(daview);
189: if (!xcoor) {
190: DMDASetUniformCoordinates(daview,0.0,1.0,0.0,1.0,0.0,1.0);
191: DMGetCoordinates(daview,&xcoor);
192: }
193: if (dacoord) {
194: DMSetUp(dacoord);
195: DMCreateLocalVector(dacoord,&xcoorl);
196: DMGlobalToLocalBegin(dacoord,xcoor,INSERT_VALUES,xcoorl);
197: DMGlobalToLocalEnd(dacoord,xcoor,INSERT_VALUES,xcoorl);
198: DMDestroy(&dacoord);
199: } else {
200: PetscInt ien,jen,ken,nc,nl,cdof,deg;
201: char fecmesh[64];
202: const char *name;
203: PetscBool flg;
205: DMDAGetNumElementsGhosted(daview,&ien,&jen,&ken);
206: nc = ien*(jen>0 ? jen : 1)*(ken>0 ? ken : 1);
208: VecGetLocalSize(xcoor,&nl);
210: VecDuplicate(xcoor,&xcoorl);
211: VecCopy(xcoor,xcoorl);
212: VecSetDM(xcoorl,NULL);
213: PetscObjectGetName((PetscObject)xcoor,&name);
214: PetscStrbeginswith(name,"FiniteElementCollection:",&flg);
215: if (!flg) {
216: deg = 0;
217: if (nc && nl) {
218: cdof = nl/(nc*dim);
219: deg = 1;
220: while (1) {
221: PetscInt degd = 1;
222: for (i=0;i<dim;i++) degd *= (deg+1);
224: if (degd == cdof) break;
225: deg++;
226: }
227: }
228: PetscSNPrintf(fecmesh,sizeof(fecmesh),"FiniteElementCollection: L2_T1_%DD_P%D",dim,deg);
229: PetscObjectSetName((PetscObject)xcoorl,fecmesh);
230: } else {
231: PetscObjectSetName((PetscObject)xcoorl,name);
232: }
233: }
235: /* xcoorl is composed with the ghosted DMDA, the ghosted coordinate DMDA (if present) is only available through this vector */
236: PetscObjectCompose((PetscObject)daview,"GLVisGraphicsCoordsGhosted",(PetscObject)xcoorl);
237: PetscObjectDereference((PetscObject)xcoorl);
239: /* daview is composed with the original DMDA */
240: PetscObjectCompose(oda,"GLVisGraphicsDMDAGhosted",(PetscObject)daview);
241: PetscObjectDereference((PetscObject)daview);
242: }
244: /* customize the viewer if present */
245: if (viewer) {
246: DMDAFieldGLVisViewerCtx *ctx;
247: DMDAGhostedGLVisViewerCtx *dactx;
248: char fec[64];
249: Vec xlocal,*Ufield;
250: const char **dafieldname;
251: char **fec_type,**fieldname;
252: PetscInt *nlocal,*bss,*dims;
253: PetscInt dim,M,N,P,dof,s,i,nf;
254: PetscBool bsset;
256: DMDAGetInfo(daview,&dim,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
257: DMGetApplicationContext(daview,&dactx);
258: DMCreateLocalVector(daview,&xlocal);
259: DMDAGetFieldNames(da,(const char * const **)&dafieldname);
260: DMDAGetNumVerticesGhosted(daview,&M,&N,&P);
261: PetscSNPrintf(fec,sizeof(fec),"FiniteElementCollection: H1_%DD_P1",dim);
262: PetscMalloc6(dof,&fec_type,dof,&nlocal,dof,&bss,dof,&dims,dof,&fieldname,dof,&Ufield);
263: for (i=0;i<dof;i++) bss[i] = 1;
264: nf = dof;
266: PetscOptionsBegin(PetscObjectComm(oda),oda->prefix,"GLVis PetscViewer DMDA Field options","PetscViewer");
267: PetscOptionsIntArray("-viewer_glvis_dm_da_bs","Block sizes for subfields; enable vector representation",NULL,bss,&nf,&bsset);
268: PetscOptionsEnd();
269: if (bsset) {
270: PetscInt t;
271: for (i=0,t=0;i<nf;i++) t += bss[i];
273: } else nf = dof;
275: for (i=0,s=0;i<nf;i++) {
276: PetscStrallocpy(fec,&fec_type[i]);
277: if (bss[i] == 1) {
278: PetscStrallocpy(dafieldname[s],&fieldname[i]);
279: } else {
280: PetscInt b;
281: size_t tlen = 9; /* "Vector-" + end */
282: for (b=0;b<bss[i];b++) {
283: size_t len;
284: PetscStrlen(dafieldname[s+b],&len);
285: tlen += len + 1; /* field + "-" */
286: }
287: PetscMalloc1(tlen,&fieldname[i]);
288: PetscStrcpy(fieldname[i],"Vector-");
289: for (b=0;b<bss[i]-1;b++) {
290: PetscStrcat(fieldname[i],dafieldname[s+b]);
291: PetscStrcat(fieldname[i],"-");
292: }
293: PetscStrcat(fieldname[i],dafieldname[s+b]);
294: }
295: dims[i] = dim;
296: nlocal[i] = M*N*P*bss[i];
297: s += bss[i];
298: }
300: /* the viewer context takes ownership of xlocal and destroys it in DMDAFieldDestroyGLVisViewerCtx_Private */
301: PetscNew(&ctx);
302: ctx->xlocal = xlocal;
304: /* create work vectors */
305: for (i=0;i<nf;i++) {
306: VecCreateMPI(PetscObjectComm((PetscObject)da),nlocal[i],PETSC_DECIDE,&Ufield[i]);
307: PetscObjectSetName((PetscObject)Ufield[i],fieldname[i]);
308: VecSetBlockSize(Ufield[i],bss[i]);
309: VecSetDM(Ufield[i],da);
310: }
312: PetscViewerGLVisSetFields(viewer,nf,(const char**)fec_type,dims,DMDASampleGLVisFields_Private,(PetscObject*)Ufield,ctx,DMDAFieldDestroyGLVisViewerCtx_Private);
313: for (i=0;i<nf;i++) {
314: PetscFree(fec_type[i]);
315: PetscFree(fieldname[i]);
316: VecDestroy(&Ufield[i]);
317: }
318: PetscFree6(fec_type,nlocal,bss,dims,fieldname,Ufield);
319: }
320: return 0;
321: }
323: static PetscErrorCode DMDAView_GLVis_ASCII(DM dm, PetscViewer viewer)
324: {
325: DM da,cda;
326: Vec xcoorl;
327: PetscMPIInt size;
328: const PetscScalar *array;
329: PetscContainer glvis_container;
330: PetscInt dim,sdim,i,vid[8],mid,cid,cdof;
331: PetscInt sx,sy,sz,ie,je,ke,ien,jen,ken,nel;
332: PetscInt gsx,gsy,gsz,gm,gn,gp,kst,jst,ist;
333: PetscBool enabled = PETSC_TRUE, isascii;
334: const char *fmt;
338: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
340: MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&size);
342: DMGetDimension(dm,&dim);
344: /* get container: determines if a process visualizes is portion of the data or not */
345: PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);
347: {
348: PetscViewerGLVisInfo glvis_info;
349: PetscContainerGetPointer(glvis_container,(void**)&glvis_info);
350: enabled = glvis_info->enabled;
351: fmt = glvis_info->fmt;
352: }
353: /* this can happen if we are calling DMView outside of VecView_GLVis */
354: PetscObjectQuery((PetscObject)dm,"GLVisGraphicsDMDAGhosted",(PetscObject*)&da);
355: if (!da) DMSetUpGLVisViewer_DMDA((PetscObject)dm,NULL);
356: PetscObjectQuery((PetscObject)dm,"GLVisGraphicsDMDAGhosted",(PetscObject*)&da);
358: DMGetCoordinateDim(da,&sdim);
360: PetscViewerASCIIPrintf(viewer,"MFEM mesh v1.0\n");
361: PetscViewerASCIIPrintf(viewer,"\ndimension\n");
362: PetscViewerASCIIPrintf(viewer,"%D\n",dim);
364: if (!enabled) {
365: PetscViewerASCIIPrintf(viewer,"\nelements\n");
366: PetscViewerASCIIPrintf(viewer,"%D\n",0);
367: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
368: PetscViewerASCIIPrintf(viewer,"%D\n",0);
369: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
370: PetscViewerASCIIPrintf(viewer,"%D\n",0);
371: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
372: return 0;
373: }
375: DMDAGetNumElementsGhosted(da,&ien,&jen,&ken);
376: nel = ien;
377: if (dim > 1) nel *= jen;
378: if (dim > 2) nel *= ken;
379: PetscViewerASCIIPrintf(viewer,"\nelements\n");
380: PetscViewerASCIIPrintf(viewer,"%D\n",nel);
381: switch (dim) {
382: case 1:
383: for (ie = 0; ie < ien; ie++) {
384: vid[0] = ie;
385: vid[1] = ie+1;
386: mid = 1; /* material id */
387: cid = 1; /* segment */
388: PetscViewerASCIIPrintf(viewer,"%D %D %D %D\n",mid,cid,vid[0],vid[1]);
389: }
390: break;
391: case 2:
392: for (je = 0; je < jen; je++) {
393: for (ie = 0; ie < ien; ie++) {
394: vid[0] = je*(ien+1) + ie;
395: vid[1] = je*(ien+1) + ie+1;
396: vid[2] = (je+1)*(ien+1) + ie+1;
397: vid[3] = (je+1)*(ien+1) + ie;
398: mid = 1; /* material id */
399: cid = 3; /* quad */
400: PetscViewerASCIIPrintf(viewer,"%D %D %D %D %D %D\n",mid,cid,vid[0],vid[1],vid[2],vid[3]);
401: }
402: }
403: break;
404: case 3:
405: for (ke = 0; ke < ken; ke++) {
406: for (je = 0; je < jen; je++) {
407: for (ie = 0; ie < ien; ie++) {
408: vid[0] = ke*(jen+1)*(ien+1) + je*(ien+1) + ie;
409: vid[1] = ke*(jen+1)*(ien+1) + je*(ien+1) + ie+1;
410: vid[2] = ke*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie+1;
411: vid[3] = ke*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie;
412: vid[4] = (ke+1)*(jen+1)*(ien+1) + je*(ien+1) + ie;
413: vid[5] = (ke+1)*(jen+1)*(ien+1) + je*(ien+1) + ie+1;
414: vid[6] = (ke+1)*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie+1;
415: vid[7] = (ke+1)*(jen+1)*(ien+1) + (je+1)*(ien+1) + ie;
416: mid = 1; /* material id */
417: cid = 5; /* hex */
418: PetscViewerASCIIPrintf(viewer,"%D %D %D %D %D %D %D %D %D %D\n",mid,cid,vid[0],vid[1],vid[2],vid[3],vid[4],vid[5],vid[6],vid[7]);
419: }
420: }
421: }
422: break;
423: default:
424: SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Unsupported dimension %D",dim);
425: }
426: PetscViewerASCIIPrintf(viewer,"\nboundary\n");
427: PetscViewerASCIIPrintf(viewer,"%D\n",0);
429: /* vertex coordinates */
430: PetscObjectQuery((PetscObject)da,"GLVisGraphicsCoordsGhosted",(PetscObject*)&xcoorl);
432: DMDAGetNumVerticesGhosted(da,&ien,&jen,&ken);
433: PetscViewerASCIIPrintf(viewer,"\nvertices\n");
434: PetscViewerASCIIPrintf(viewer,"%D\n",ien*jen*ken);
435: if (nel) {
436: VecGetDM(xcoorl,&cda);
437: VecGetArrayRead(xcoorl,&array);
438: if (!cda) { /* HO viz */
439: const char *fecname;
440: PetscInt nc,nl;
442: PetscObjectGetName((PetscObject)xcoorl,&fecname);
443: PetscViewerASCIIPrintf(viewer,"nodes\n");
444: PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
445: PetscViewerASCIIPrintf(viewer,"%s\n",fecname);
446: PetscViewerASCIIPrintf(viewer,"VDim: %D\n",sdim);
447: PetscViewerASCIIPrintf(viewer,"Ordering: 1\n\n"); /*Ordering::byVDIM*/
448: /* L2 coordinates */
449: DMDAGetNumElementsGhosted(da,&ien,&jen,&ken);
450: VecGetLocalSize(xcoorl,&nl);
451: nc = ien*(jen>0 ? jen : 1)*(ken>0 ? ken : 1);
452: cdof = nc ? nl/nc : 0;
453: if (!ien) ien++;
454: if (!jen) jen++;
455: if (!ken) ken++;
456: ist = jst = kst = 0;
457: gm = ien;
458: gn = jen;
459: gp = ken;
460: } else {
461: DMDAGhostedGLVisViewerCtx *dactx;
463: DMGetApplicationContext(da,&dactx);
464: PetscViewerASCIIPrintf(viewer,"%D\n",sdim);
465: cdof = sdim;
466: DMDAGetCorners(da,&sx,&sy,&sz,NULL,NULL,NULL);
467: DMDAGetGhostCorners(da,&gsx,&gsy,&gsz,&gm,&gn,&gp);
468: if (dactx->ll) {
469: kst = jst = ist = 0;
470: } else {
471: kst = gsz != sz ? 1 : 0;
472: jst = gsy != sy ? 1 : 0;
473: ist = gsx != sx ? 1 : 0;
474: }
475: }
476: for (ke = kst; ke < kst + ken; ke++) {
477: for (je = jst; je < jst + jen; je++) {
478: for (ie = ist; ie < ist + ien; ie++) {
479: PetscInt c;
481: i = ke * gm * gn + je * gm + ie;
482: for (c=0;c<cdof/sdim;c++) {
483: PetscInt d;
484: for (d=0;d<sdim;d++) {
485: PetscViewerASCIIPrintf(viewer,fmt,PetscRealPart(array[cdof*i+c*sdim+d]));
486: }
487: PetscViewerASCIIPrintf(viewer,"\n");
488: }
489: }
490: }
491: }
492: VecRestoreArrayRead(xcoorl,&array);
493: }
494: return 0;
495: }
497: PetscErrorCode DMView_DA_GLVis(DM dm, PetscViewer viewer)
498: {
499: DMView_GLVis(dm,viewer,DMDAView_GLVis_ASCII);
500: return 0;
501: }