Actual source code: ex50.c
petsc-3.12.4 2020-02-04
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 <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*/