1 static char help[] = "Test GLVis high-order support with DMDAs\n\n";
2 
3 #include <petscdm.h>
4 #include <petscdmda.h>
5 #include <petscdmplex.h>
6 #include <petscdt.h>
7 
MapPoint(PetscScalar xyz[],PetscScalar mxyz[])8 static PetscErrorCode MapPoint(PetscScalar xyz[],PetscScalar mxyz[])
9 {
10   PetscScalar x,y,z,x2,y2,z2;
11 
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 }
23 
test_3d(PetscInt cells[],PetscBool plex,PetscBool ho)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};
31   PetscErrorCode ierr;
32 
33   PetscFunctionBeginUser;
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;
45 
46     ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,3,PETSC_FALSE,cells,l,u,NULL,PETSC_FALSE,&dm);CHKERRQ(ierr);
47     ierr = DMPlexGetPartitioner(dm,&part);CHKERRQ(ierr);
48     ierr = PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);CHKERRQ(ierr);
49     ierr = DMPlexDistribute(dm,0,NULL,&dm2);CHKERRQ(ierr);
50     if (dm2) {
51       ierr = DMDestroy(&dm);CHKERRQ(ierr);
52       dm   = dm2;
53     }
54   } else {
55     ierr = 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);CHKERRQ(ierr);
56   }
57   ierr = DMSetUp(dm);CHKERRQ(ierr);
58   if (!plex) {
59     ierr = DMDASetUniformCoordinates(dm,l[0],u[0],l[1],u[1],l[2],u[2]);CHKERRQ(ierr);
60   }
61   if (ho) { /* each element mapped to a sphere */
62     DM           cdm;
63     Vec          cv;
64     PetscSection cSec;
65     DMDACoor3d   ***_coords;
66     PetscScalar  shift[3],*cptr;
67     PetscInt     nel,dof = 3,nex,ney,nez,gx = 0,gy = 0,gz = 0;
68     PetscInt     i,j,k,pi,pj,pk;
69     PetscReal    *nodes,*weights;
70     char         name[256];
71 
72     ierr = PetscOptionsGetInt(NULL,NULL,"-order",&dof,NULL);CHKERRQ(ierr);
73     dof += 1;
74 
75     ierr = PetscMalloc1(dof,&nodes);CHKERRQ(ierr);
76     ierr = PetscMalloc1(dof,&weights);CHKERRQ(ierr);
77     ierr = PetscDTGaussLobattoLegendreQuadrature(dof,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,nodes,weights);CHKERRQ(ierr);
78     ierr = DMGetCoordinatesLocal(dm,&cv);CHKERRQ(ierr);
79     ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
80     if (plex) {
81       PetscInt cEnd,cStart;
82 
83       ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
84       ierr = DMGetCoordinateSection(dm,&cSec);CHKERRQ(ierr);
85       nel  = cEnd - cStart;
86       nex  = nel;
87       ney  = 1;
88       nez  = 1;
89     } else {
90       ierr = DMDAVecGetArray(cdm,cv,&_coords);CHKERRQ(ierr);
91       ierr = DMDAGetElementsCorners(dm,&gx,&gy,&gz);CHKERRQ(ierr);
92       ierr = DMDAGetElementsSizes(dm,&nex,&ney,&nez);CHKERRQ(ierr);
93       nel  = nex*ney*nez;
94     }
95     ierr = VecCreate(PETSC_COMM_WORLD,&v);CHKERRQ(ierr);
96     ierr = VecSetSizes(v,3*dof*dof*dof*nel,PETSC_DECIDE);CHKERRQ(ierr);
97     ierr = VecSetType(v,VECSTANDARD);CHKERRQ(ierr);
98     ierr = VecGetArray(v,&c);CHKERRQ(ierr);
99     cptr = c;
100     for (k=gz;k<gz+nez;k++) {
101       for (j=gy;j<gy+ney;j++) {
102         for (i=gx;i<gx+nex;i++) {
103           if (plex) {
104             PetscScalar *t = NULL;
105 
106             ierr = DMPlexVecGetClosure(dm,cSec,cv,i,NULL,&t);CHKERRQ(ierr);
107             shift[0] = t[0];
108             shift[1] = t[1];
109             shift[2] = t[2];
110             ierr = DMPlexVecRestoreClosure(dm,cSec,cv,i,NULL,&t);CHKERRQ(ierr);
111           } else {
112             shift[0] = _coords[k][j][i].x;
113             shift[1] = _coords[k][j][i].y;
114             shift[2] = _coords[k][j][i].z;
115           }
116           for (pk=0;pk<dof;pk++) {
117             PetscScalar xyz[3];
118 
119             xyz[2] = nodes[pk];
120             for (pj=0;pj<dof;pj++) {
121               xyz[1] = nodes[pj];
122               for (pi=0;pi<dof;pi++) {
123                 xyz[0] = nodes[pi];
124                 ierr = MapPoint(xyz,cptr);CHKERRQ(ierr);
125                 cptr[0] += shift[0];
126                 cptr[1] += shift[1];
127                 cptr[2] += shift[2];
128                 cptr += 3;
129               }
130             }
131           }
132         }
133       }
134     }
135     if (!plex) {
136       ierr = DMDAVecRestoreArray(cdm,cv,&_coords);CHKERRQ(ierr);
137     }
138     ierr = VecRestoreArray(v,&c);CHKERRQ(ierr);
139     ierr = PetscSNPrintf(name,sizeof(name),"FiniteElementCollection: L2_T1_3D_P%D",dof-1);CHKERRQ(ierr);
140     ierr = PetscObjectSetName((PetscObject)v,name);CHKERRQ(ierr);
141     ierr = PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)v);CHKERRQ(ierr);
142     ierr = VecDestroy(&v);CHKERRQ(ierr);
143     ierr = PetscFree(nodes);CHKERRQ(ierr);
144     ierr = PetscFree(weights);CHKERRQ(ierr);
145     ierr = DMViewFromOptions(dm,NULL,"-view");CHKERRQ(ierr);
146   } else { /* map the whole domain to a sphere */
147     ierr = DMGetCoordinates(dm,&v);CHKERRQ(ierr);
148     ierr = VecGetLocalSize(v,&nl);CHKERRQ(ierr);
149     ierr = VecGetArray(v,&c);CHKERRQ(ierr);
150     for (i=0;i<nl/3;i++) {
151       ierr = MapPoint(c+3*i,c+3*i);CHKERRQ(ierr);
152     }
153     ierr = VecRestoreArray(v,&c);CHKERRQ(ierr);
154     ierr = DMSetCoordinates(dm,v);CHKERRQ(ierr);
155     ierr = DMViewFromOptions(dm,NULL,"-view");CHKERRQ(ierr);
156   }
157   ierr = DMDestroy(&dm);CHKERRQ(ierr);
158   PetscFunctionReturn(0);
159 }
160 
main(int argc,char * argv[])161 int main(int argc, char *argv[])
162 {
163   PetscErrorCode ierr;
164   PetscBool      ho = PETSC_TRUE;
165   PetscBool      plex = PETSC_FALSE;
166   PetscInt       cells[3] = {2,2,2};
167 
168   ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
169   ierr = PetscOptionsGetBool(NULL,NULL,"-ho",&ho,NULL);CHKERRQ(ierr);
170   ierr = PetscOptionsGetBool(NULL,NULL,"-plex",&plex,NULL);CHKERRQ(ierr);
171   ierr = PetscOptionsGetInt(NULL,NULL,"-nex",&cells[0],NULL);CHKERRQ(ierr);
172   ierr = PetscOptionsGetInt(NULL,NULL,"-ney",&cells[1],NULL);CHKERRQ(ierr);
173   ierr = PetscOptionsGetInt(NULL,NULL,"-nez",&cells[2],NULL);CHKERRQ(ierr);
174   ierr = test_3d(cells,plex,ho);CHKERRQ(ierr);
175   ierr = PetscFinalize();
176   return ierr;
177 }
178 
179 /*TEST
180 
181    testset:
182      nsize: 1
183      args: -view glvis:
184 
185      test:
186         suffix: dmda_glvis_ho
187         args: -nex 1 -ney 1 -nez 1
188 
189      test:
190         suffix: dmda_glvis_lo
191         args: -ho 0
192 
193      test:
194         suffix: dmplex_glvis_ho
195         args: -nex 1 -ney 1 -nez 1
196 
197      test:
198         suffix: dmplex_glvis_lo
199         args: -ho 0
200 
201 TEST*/
202