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 PetscErrorCode ierr; 32 33 PetscFunctionBeginUser; 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; 45 46 ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,3,PETSC_FALSE,cells,l,u,NULL,PETSC_FALSE,&dm);CHKERRQ(ierr); 47 ierr = DMPlexGetPartitioner(dm,&part);CHKERRQ(ierr); 48 ierr = PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);CHKERRQ(ierr); 49 ierr = DMPlexDistribute(dm,0,NULL,&dm2);CHKERRQ(ierr); 50 if (dm2) { 51 ierr = DMDestroy(&dm);CHKERRQ(ierr); 52 dm = dm2; 53 } 54 } else { 55 ierr = 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);CHKERRQ(ierr); 56 } 57 ierr = DMSetUp(dm);CHKERRQ(ierr); 58 if (!plex) { 59 ierr = DMDASetUniformCoordinates(dm,l[0],u[0],l[1],u[1],l[2],u[2]);CHKERRQ(ierr); 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]; 71 72 ierr = PetscOptionsGetInt(NULL,NULL,"-order",&dof,NULL);CHKERRQ(ierr); 73 dof += 1; 74 75 ierr = PetscMalloc1(dof,&nodes);CHKERRQ(ierr); 76 ierr = PetscMalloc1(dof,&weights);CHKERRQ(ierr); 77 ierr = PetscDTGaussLobattoLegendreQuadrature(dof,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,nodes,weights);CHKERRQ(ierr); 78 ierr = DMGetCoordinatesLocal(dm,&cv);CHKERRQ(ierr); 79 ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr); 80 if (plex) { 81 PetscInt cEnd,cStart; 82 83 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 84 ierr = DMGetCoordinateSection(dm,&cSec);CHKERRQ(ierr); 85 nel = cEnd - cStart; 86 nex = nel; 87 ney = 1; 88 nez = 1; 89 } else { 90 ierr = DMDAVecGetArray(cdm,cv,&_coords);CHKERRQ(ierr); 91 ierr = DMDAGetElementsCorners(dm,&gx,&gy,&gz);CHKERRQ(ierr); 92 ierr = DMDAGetElementsSizes(dm,&nex,&ney,&nez);CHKERRQ(ierr); 93 nel = nex*ney*nez; 94 } 95 ierr = VecCreate(PETSC_COMM_WORLD,&v);CHKERRQ(ierr); 96 ierr = VecSetSizes(v,3*dof*dof*dof*nel,PETSC_DECIDE);CHKERRQ(ierr); 97 ierr = VecSetType(v,VECSTANDARD);CHKERRQ(ierr); 98 ierr = VecGetArray(v,&c);CHKERRQ(ierr); 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; 105 106 ierr = DMPlexVecGetClosure(dm,cSec,cv,i,NULL,&t);CHKERRQ(ierr); 107 shift[0] = t[0]; 108 shift[1] = t[1]; 109 shift[2] = t[2]; 110 ierr = DMPlexVecRestoreClosure(dm,cSec,cv,i,NULL,&t);CHKERRQ(ierr); 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]; 118 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 ierr = MapPoint(xyz,cptr);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(cdm,cv,&_coords);CHKERRQ(ierr); 137 } 138 ierr = VecRestoreArray(v,&c);CHKERRQ(ierr); 139 ierr = PetscSNPrintf(name,sizeof(name),"FiniteElementCollection: L2_T1_3D_P%D",dof-1);CHKERRQ(ierr); 140 ierr = PetscObjectSetName((PetscObject)v,name);CHKERRQ(ierr); 141 ierr = PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)v);CHKERRQ(ierr); 142 ierr = VecDestroy(&v);CHKERRQ(ierr); 143 ierr = PetscFree(nodes);CHKERRQ(ierr); 144 ierr = PetscFree(weights);CHKERRQ(ierr); 145 ierr = DMViewFromOptions(dm,NULL,"-view");CHKERRQ(ierr); 146 } else { /* map the whole domain to a sphere */ 147 ierr = DMGetCoordinates(dm,&v);CHKERRQ(ierr); 148 ierr = VecGetLocalSize(v,&nl);CHKERRQ(ierr); 149 ierr = VecGetArray(v,&c);CHKERRQ(ierr); 150 for (i=0;i<nl/3;i++) { 151 ierr = MapPoint(c+3*i,c+3*i);CHKERRQ(ierr); 152 } 153 ierr = VecRestoreArray(v,&c);CHKERRQ(ierr); 154 ierr = DMSetCoordinates(dm,v);CHKERRQ(ierr); 155 ierr = DMViewFromOptions(dm,NULL,"-view");CHKERRQ(ierr); 156 } 157 ierr = DMDestroy(&dm);CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 int main(int argc, char *argv[]) 162 { 163 PetscErrorCode ierr; 164 PetscBool ho = PETSC_TRUE; 165 PetscBool plex = PETSC_FALSE; 166 PetscInt cells[3] = {2,2,2}; 167 168 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 169 ierr = PetscOptionsGetBool(NULL,NULL,"-ho",&ho,NULL);CHKERRQ(ierr); 170 ierr = PetscOptionsGetBool(NULL,NULL,"-plex",&plex,NULL);CHKERRQ(ierr); 171 ierr = PetscOptionsGetInt(NULL,NULL,"-nex",&cells[0],NULL);CHKERRQ(ierr); 172 ierr = PetscOptionsGetInt(NULL,NULL,"-ney",&cells[1],NULL);CHKERRQ(ierr); 173 ierr = PetscOptionsGetInt(NULL,NULL,"-nez",&cells[2],NULL);CHKERRQ(ierr); 174 ierr = test_3d(cells,plex,ho);CHKERRQ(ierr); 175 ierr = PetscFinalize(); 176 return ierr; 177 } 178 179 /*TEST 180 181 testset: 182 nsize: 1 183 args: -view glvis: 184 185 test: 186 suffix: dmda_glvis_ho 187 args: -nex 1 -ney 1 -nez 1 188 189 test: 190 suffix: dmda_glvis_lo 191 args: -ho 0 192 193 test: 194 suffix: dmplex_glvis_ho 195 args: -nex 1 -ney 1 -nez 1 196 197 test: 198 suffix: dmplex_glvis_lo 199 args: -ho 0 200 201 TEST*/ 202