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