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