Actual source code: ex50.c
petsc-3.11.1 2019-04-12
1: static char help[] = "Test GLVis high-order support with DMDAs\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
5: #include <petscdmplex.h>
6: #include <petscgll.h>
8: static PetscErrorCode MapPoint(PetscScalar xyz[],PetscScalar mxyz[])
9: {
10: PetscScalar x,y,z,x2,y2,z2;
12: x = xyz[0];
13: y = xyz[1];
14: z = xyz[2];
15: x2 = x*x;
16: y2 = y*y;
17: z2 = z*z;
18: mxyz[0] = x*PetscSqrtScalar(1.0-y2/2.0-z2/2.0 + y2*z2/3.0);
19: mxyz[1] = y*PetscSqrtScalar(1.0-z2/2.0-x2/2.0 + z2*x2/3.0);
20: mxyz[2] = z*PetscSqrtScalar(1.0-x2/2.0-y2/2.0 + x2*y2/3.0);
21: return 0;
22: }
24: static PetscErrorCode test_3d(PetscInt cells[], PetscBool plex, PetscBool ho)
25: {
26: DM dm;
27: Vec v;
28: PetscScalar *c;
29: PetscInt nl,i;
30: PetscReal u[3] = {1.0,1.0,1.0}, l[3] = {-1.0,-1.0,-1.0};
34: if (ho) {
35: u[0] = 2.*cells[0];
36: u[1] = 2.*cells[1];
37: u[2] = 2.*cells[2];
38: l[0] = 0.;
39: l[1] = 0.;
40: l[2] = 0.;
41: }
42: if (plex) {
43: DM dm2;
44: PetscPartitioner part;
46: DMPlexCreateBoxMesh(PETSC_COMM_WORLD,3,PETSC_FALSE,cells,l,u,NULL,PETSC_FALSE,&dm);
47: DMPlexGetPartitioner(dm,&part);
48: PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
49: DMPlexDistribute(dm,0,NULL,&dm2);
50: if (dm2) {
51: DMDestroy(&dm);
52: dm = dm2;
53: }
54: } else {
55: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,cells[0]+1,cells[1]+1,cells[2]+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&dm);
56: }
57: DMSetUp(dm);
58: if (!plex) {
59: DMDASetUniformCoordinates(dm,l[0],u[0],l[1],u[1],l[2],u[2]);
60: }
61: if (ho) { /* each element mapped to a sphere */
62: DM cdm;
63: Vec cv;
64: PetscGLL gll;
65: PetscSection cSec;
66: DMDACoor3d ***_coords;
67: PetscScalar shift[3],*cptr;
68: PetscInt nel,dof = 3,nex,ney,nez,gx = 0,gy = 0,gz = 0;
69: PetscInt i,j,k,pi,pj,pk;
70: char name[256];
72: PetscOptionsGetInt(NULL,NULL,"-order",&dof,NULL);
73: dof += 1;
74: PetscGLLCreate(dof,PETSCGLL_VIA_LINEARALGEBRA,&gll);
75: DMGetCoordinatesLocal(dm,&cv);
76: DMGetCoordinateDM(dm,&cdm);
77: if (plex) {
78: PetscInt cEnd,cStart;
80: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
81: DMGetCoordinateSection(dm,&cSec);
82: nel = cEnd - cStart;
83: nex = nel;
84: ney = 1;
85: nez = 1;
86: } else {
87: DMDAVecGetArray(cdm,cv,&_coords);
88: DMDAGetElementsCorners(dm,&gx,&gy,&gz);
89: DMDAGetElementsSizes(dm,&nex,&ney,&nez);
90: nel = nex*ney*nez;
91: }
92: VecCreate(PETSC_COMM_WORLD,&v);
93: VecSetSizes(v,3*dof*dof*dof*nel,PETSC_DECIDE);
94: VecSetType(v,VECSTANDARD);
95: VecGetArray(v,&c);
96: cptr = c;
97: for (k=gz;k<gz+nez;k++) {
98: for (j=gy;j<gy+ney;j++) {
99: for (i=gx;i<gx+nex;i++) {
100: if (plex) {
101: PetscScalar *t = NULL;
103: DMPlexVecGetClosure(dm,cSec,cv,i,NULL,&t);
104: shift[0] = t[0];
105: shift[1] = t[1];
106: shift[2] = t[2];
107: DMPlexVecRestoreClosure(dm,cSec,cv,i,NULL,&t);
108: } else {
109: shift[0] = _coords[k][j][i].x;
110: shift[1] = _coords[k][j][i].y;
111: shift[2] = _coords[k][j][i].z;
112: }
113: for (pk=0;pk<dof;pk++) {
114: PetscScalar xyz[3];
116: xyz[2] = gll.nodes[pk];
117: for (pj=0;pj<dof;pj++) {
118: xyz[1] = gll.nodes[pj];
119: for (pi=0;pi<dof;pi++) {
120: xyz[0] = gll.nodes[pi];
121: MapPoint(xyz,cptr);
122: cptr[0] += shift[0];
123: cptr[1] += shift[1];
124: cptr[2] += shift[2];
125: cptr += 3;
126: }
127: }
128: }
129: }
130: }
131: }
132: PetscGLLDestroy(&gll);
133: if (!plex) {
134: DMDAVecRestoreArray(cdm,cv,&_coords);
135: }
136: VecRestoreArray(v,&c);
137: PetscSNPrintf(name,sizeof(name),"FiniteElementCollection: L2_T1_3D_P%D",dof-1);
138: PetscObjectSetName((PetscObject)v,name);
139: PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)v);
140: VecDestroy(&v);
141: DMViewFromOptions(dm,NULL,"-view");
142: } else { /* map the whole domain to a sphere */
143: DMGetCoordinates(dm,&v);
144: VecGetLocalSize(v,&nl);
145: VecGetArray(v,&c);
146: for (i=0;i<nl/3;i++) {
147: MapPoint(c+3*i,c+3*i);
148: }
149: VecRestoreArray(v,&c);
150: DMSetCoordinates(dm,v);
151: DMViewFromOptions(dm,NULL,"-view");
152: }
153: DMDestroy(&dm);
154: return(0);
155: }
157: int main(int argc, char *argv[])
158: {
160: PetscBool ho = PETSC_TRUE;
161: PetscBool plex = PETSC_FALSE;
162: PetscInt cells[3] = {2,2,2};
164: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
165: PetscOptionsGetBool(NULL,NULL,"-ho",&ho,NULL);
166: PetscOptionsGetBool(NULL,NULL,"-plex",&plex,NULL);
167: PetscOptionsGetInt(NULL,NULL,"-nex",&cells[0],NULL);
168: PetscOptionsGetInt(NULL,NULL,"-ney",&cells[1],NULL);
169: PetscOptionsGetInt(NULL,NULL,"-nez",&cells[2],NULL);
170: test_3d(cells,plex,ho);
171: PetscFinalize();
172: return ierr;
173: }
175: /*TEST
177: testset:
178: nsize: 1
179: args: -view glvis:
181: test:
182: suffix: dmda_glvis_ho
183: args: -nex 1 -ney 1 -nez 1
185: test:
186: suffix: dmda_glvis_lo
187: args: -ho 0
189: test:
190: suffix: dmplex_glvis_ho
191: args: -nex 1 -ney 1 -nez 1
193: test:
194: suffix: dmplex_glvis_lo
195: args: -ho 0
197: TEST*/