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