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