Actual source code: ex50.c

petsc-3.11.1 2019-04-12
Report Typos and Errors
  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 <petscgll.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:     PetscGLL     gll;
 65:     PetscSection cSec;
 66:     DMDACoor3d   ***_coords;
 67:     PetscScalar  shift[3],*cptr;
 68:     PetscInt     nel,dof = 3,nex,ney,nez,gx = 0,gy = 0,gz = 0;
 69:     PetscInt     i,j,k,pi,pj,pk;
 70:     char         name[256];

 72:     PetscOptionsGetInt(NULL,NULL,"-order",&dof,NULL);
 73:     dof += 1;
 74:     PetscGLLCreate(dof,PETSCGLL_VIA_LINEARALGEBRA,&gll);
 75:     DMGetCoordinatesLocal(dm,&cv);
 76:     DMGetCoordinateDM(dm,&cdm);
 77:     if (plex) {
 78:       PetscInt cEnd,cStart;

 80:       DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 81:       DMGetCoordinateSection(dm,&cSec);
 82:       nel  = cEnd - cStart;
 83:       nex  = nel;
 84:       ney  = 1;
 85:       nez  = 1;
 86:     } else {
 87:       DMDAVecGetArray(cdm,cv,&_coords);
 88:       DMDAGetElementsCorners(dm,&gx,&gy,&gz);
 89:       DMDAGetElementsSizes(dm,&nex,&ney,&nez);
 90:       nel  = nex*ney*nez;
 91:     }
 92:     VecCreate(PETSC_COMM_WORLD,&v);
 93:     VecSetSizes(v,3*dof*dof*dof*nel,PETSC_DECIDE);
 94:     VecSetType(v,VECSTANDARD);
 95:     VecGetArray(v,&c);
 96:     cptr = c;
 97:     for (k=gz;k<gz+nez;k++) {
 98:       for (j=gy;j<gy+ney;j++) {
 99:         for (i=gx;i<gx+nex;i++) {
100:           if (plex) {
101:             PetscScalar *t = NULL;

103:             DMPlexVecGetClosure(dm,cSec,cv,i,NULL,&t);
104:             shift[0] = t[0];
105:             shift[1] = t[1];
106:             shift[2] = t[2];
107:             DMPlexVecRestoreClosure(dm,cSec,cv,i,NULL,&t);
108:           } else {
109:             shift[0] = _coords[k][j][i].x;
110:             shift[1] = _coords[k][j][i].y;
111:             shift[2] = _coords[k][j][i].z;
112:           }
113:           for (pk=0;pk<dof;pk++) {
114:             PetscScalar xyz[3];

116:             xyz[2] = gll.nodes[pk];
117:             for (pj=0;pj<dof;pj++) {
118:               xyz[1] = gll.nodes[pj];
119:               for (pi=0;pi<dof;pi++) {
120:                 xyz[0] = gll.nodes[pi];
121:                 MapPoint(xyz,cptr);
122:                 cptr[0] += shift[0];
123:                 cptr[1] += shift[1];
124:                 cptr[2] += shift[2];
125:                 cptr += 3;
126:               }
127:             }
128:           }
129:         }
130:       }
131:     }
132:     PetscGLLDestroy(&gll);
133:     if (!plex) {
134:       DMDAVecRestoreArray(cdm,cv,&_coords);
135:     }
136:     VecRestoreArray(v,&c);
137:     PetscSNPrintf(name,sizeof(name),"FiniteElementCollection: L2_T1_3D_P%D",dof-1);
138:     PetscObjectSetName((PetscObject)v,name);
139:     PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)v);
140:     VecDestroy(&v);
141:     DMViewFromOptions(dm,NULL,"-view");
142:   } else { /* map the whole domain to a sphere */
143:     DMGetCoordinates(dm,&v);
144:     VecGetLocalSize(v,&nl);
145:     VecGetArray(v,&c);
146:     for (i=0;i<nl/3;i++) {
147:       MapPoint(c+3*i,c+3*i);
148:     }
149:     VecRestoreArray(v,&c);
150:     DMSetCoordinates(dm,v);
151:     DMViewFromOptions(dm,NULL,"-view");
152:   }
153:   DMDestroy(&dm);
154:   return(0);
155: }

157: int main(int argc, char *argv[])
158: {
160:   PetscBool      ho = PETSC_TRUE;
161:   PetscBool      plex = PETSC_FALSE;
162:   PetscInt       cells[3] = {2,2,2};

164:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
165:   PetscOptionsGetBool(NULL,NULL,"-ho",&ho,NULL);
166:   PetscOptionsGetBool(NULL,NULL,"-plex",&plex,NULL);
167:   PetscOptionsGetInt(NULL,NULL,"-nex",&cells[0],NULL);
168:   PetscOptionsGetInt(NULL,NULL,"-ney",&cells[1],NULL);
169:   PetscOptionsGetInt(NULL,NULL,"-nez",&cells[2],NULL);
170:   test_3d(cells,plex,ho);
171:   PetscFinalize();
172:   return ierr;
173: }

175: /*TEST

177:    testset:
178:      nsize: 1
179:      args: -view glvis:

181:      test:
182:         suffix: dmda_glvis_ho
183:         args: -nex 1 -ney 1 -nez 1

185:      test:
186:         suffix: dmda_glvis_lo
187:         args: -ho 0

189:      test:
190:         suffix: dmplex_glvis_ho
191:         args: -nex 1 -ney 1 -nez 1

193:      test:
194:         suffix: dmplex_glvis_lo
195:         args: -ho 0

197: TEST*/