Actual source code: ex50.c
petsc-3.13.6 2020-09-29
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 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: 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];
72: PetscOptionsGetInt(NULL,NULL,"-order",&dof,NULL);
73: dof += 1;
75: PetscMalloc1(dof,&nodes);
76: PetscMalloc1(dof,&weights);
77: PetscDTGaussLobattoLegendreQuadrature(dof,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,nodes,weights);
78: DMGetCoordinatesLocal(dm,&cv);
79: DMGetCoordinateDM(dm,&cdm);
80: if (plex) {
81: PetscInt cEnd,cStart;
83: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
84: DMGetCoordinateSection(dm,&cSec);
85: nel = cEnd - cStart;
86: nex = nel;
87: ney = 1;
88: nez = 1;
89: } else {
90: DMDAVecGetArray(cdm,cv,&_coords);
91: DMDAGetElementsCorners(dm,&gx,&gy,&gz);
92: DMDAGetElementsSizes(dm,&nex,&ney,&nez);
93: nel = nex*ney*nez;
94: }
95: VecCreate(PETSC_COMM_WORLD,&v);
96: VecSetSizes(v,3*dof*dof*dof*nel,PETSC_DECIDE);
97: VecSetType(v,VECSTANDARD);
98: VecGetArray(v,&c);
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;
106: DMPlexVecGetClosure(dm,cSec,cv,i,NULL,&t);
107: shift[0] = t[0];
108: shift[1] = t[1];
109: shift[2] = t[2];
110: DMPlexVecRestoreClosure(dm,cSec,cv,i,NULL,&t);
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];
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: MapPoint(xyz,cptr);
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: DMDAVecRestoreArray(cdm,cv,&_coords);
137: }
138: VecRestoreArray(v,&c);
139: PetscSNPrintf(name,sizeof(name),"FiniteElementCollection: L2_T1_3D_P%D",dof-1);
140: PetscObjectSetName((PetscObject)v,name);
141: PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)v);
142: VecDestroy(&v);
143: PetscFree(nodes);
144: PetscFree(weights);
145: DMViewFromOptions(dm,NULL,"-view");
146: } else { /* map the whole domain to a sphere */
147: DMGetCoordinates(dm,&v);
148: VecGetLocalSize(v,&nl);
149: VecGetArray(v,&c);
150: for (i=0;i<nl/3;i++) {
151: MapPoint(c+3*i,c+3*i);
152: }
153: VecRestoreArray(v,&c);
154: DMSetCoordinates(dm,v);
155: DMViewFromOptions(dm,NULL,"-view");
156: }
157: DMDestroy(&dm);
158: return(0);
159: }
161: int main(int argc, char *argv[])
162: {
164: PetscBool ho = PETSC_TRUE;
165: PetscBool plex = PETSC_FALSE;
166: PetscInt cells[3] = {2,2,2};
168: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
169: PetscOptionsGetBool(NULL,NULL,"-ho",&ho,NULL);
170: PetscOptionsGetBool(NULL,NULL,"-plex",&plex,NULL);
171: PetscOptionsGetInt(NULL,NULL,"-nex",&cells[0],NULL);
172: PetscOptionsGetInt(NULL,NULL,"-ney",&cells[1],NULL);
173: PetscOptionsGetInt(NULL,NULL,"-nez",&cells[2],NULL);
174: test_3d(cells,plex,ho);
175: PetscFinalize();
176: return ierr;
177: }
179: /*TEST
181: testset:
182: nsize: 1
183: args: -view glvis:
185: test:
186: suffix: dmda_glvis_ho
187: args: -nex 1 -ney 1 -nez 1
189: test:
190: suffix: dmda_glvis_lo
191: args: -ho 0
193: test:
194: suffix: dmplex_glvis_ho
195: args: -nex 1 -ney 1 -nez 1
197: test:
198: suffix: dmplex_glvis_lo
199: args: -ho 0
201: TEST*/