Actual source code: ex50.c
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 <petscdt.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 PETSC_SUCCESS;
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};
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];
61: PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &dof, NULL));
62: dof += 1;
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;
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;
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];
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: }
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};
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: }
164: /*TEST
166: testset:
167: nsize: 1
168: args: -view glvis:
170: test:
171: suffix: dmda_glvis_ho
172: args: -nex 1 -ney 1 -nez 1
174: test:
175: suffix: dmda_glvis_lo
176: args: -ho 0
178: test:
179: suffix: dmplex_glvis_ho
180: args: -nex 1 -ney 1 -nez 1
182: test:
183: suffix: dmplex_glvis_lo
184: args: -ho 0
186: TEST*/