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 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 43 PetscCall(DMSetType(dm, DMPLEX)); 44 PetscCall(DMSetFromOptions(dm)); 45 } else { 46 PetscCall(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 PetscCall(DMSetUp(dm)); 49 if (!plex) { 50 PetscCall(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 PetscCall(PetscOptionsGetInt(NULL,NULL,"-order",&dof,NULL)); 64 dof += 1; 65 66 PetscCall(PetscMalloc1(dof,&nodes)); 67 PetscCall(PetscMalloc1(dof,&weights)); 68 PetscCall(PetscDTGaussLobattoLegendreQuadrature(dof,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,nodes,weights)); 69 PetscCall(DMGetCoordinatesLocal(dm,&cv)); 70 PetscCall(DMGetCoordinateDM(dm,&cdm)); 71 if (plex) { 72 PetscInt cEnd,cStart; 73 74 PetscCall(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 75 PetscCall(DMGetCoordinateSection(dm,&cSec)); 76 nel = cEnd - cStart; 77 nex = nel; 78 ney = 1; 79 nez = 1; 80 } else { 81 PetscCall(DMDAVecGetArray(cdm,cv,&_coords)); 82 PetscCall(DMDAGetElementsCorners(dm,&gx,&gy,&gz)); 83 PetscCall(DMDAGetElementsSizes(dm,&nex,&ney,&nez)); 84 nel = nex*ney*nez; 85 } 86 PetscCall(VecCreate(PETSC_COMM_WORLD,&v)); 87 PetscCall(VecSetSizes(v,3*dof*dof*dof*nel,PETSC_DECIDE)); 88 PetscCall(VecSetType(v,VECSTANDARD)); 89 PetscCall(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 PetscCall(DMPlexVecGetClosure(dm,cSec,cv,i,NULL,&t)); 98 shift[0] = t[0]; 99 shift[1] = t[1]; 100 shift[2] = t[2]; 101 PetscCall(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 PetscCall(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 PetscCall(DMDAVecRestoreArray(cdm,cv,&_coords)); 128 } 129 PetscCall(VecRestoreArray(v,&c)); 130 PetscCall(PetscSNPrintf(name,sizeof(name),"FiniteElementCollection: L2_T1_3D_P%D",dof-1)); 131 PetscCall(PetscObjectSetName((PetscObject)v,name)); 132 PetscCall(PetscObjectCompose((PetscObject)dm,"_glvis_mesh_coords",(PetscObject)v)); 133 PetscCall(VecDestroy(&v)); 134 PetscCall(PetscFree(nodes)); 135 PetscCall(PetscFree(weights)); 136 PetscCall(DMViewFromOptions(dm,NULL,"-view")); 137 } else { /* map the whole domain to a sphere */ 138 PetscCall(DMGetCoordinates(dm,&v)); 139 PetscCall(VecGetLocalSize(v,&nl)); 140 PetscCall(VecGetArray(v,&c)); 141 for (i=0;i<nl/3;i++) { 142 PetscCall(MapPoint(c+3*i,c+3*i)); 143 } 144 PetscCall(VecRestoreArray(v,&c)); 145 PetscCall(DMSetCoordinates(dm,v)); 146 PetscCall(DMViewFromOptions(dm,NULL,"-view")); 147 } 148 PetscCall(DMDestroy(&dm)); 149 PetscFunctionReturn(0); 150 } 151 152 int main(int argc, char *argv[]) 153 { 154 PetscBool ho = PETSC_TRUE; 155 PetscBool plex = PETSC_FALSE; 156 PetscInt cells[3] = {2,2,2}; 157 158 PetscCall(PetscInitialize(&argc,&argv,0,help)); 159 PetscCall(PetscOptionsGetBool(NULL,NULL,"-ho",&ho,NULL)); 160 PetscCall(PetscOptionsGetBool(NULL,NULL,"-plex",&plex,NULL)); 161 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nex",&cells[0],NULL)); 162 PetscCall(PetscOptionsGetInt(NULL,NULL,"-ney",&cells[1],NULL)); 163 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nez",&cells[2],NULL)); 164 PetscCall(test_3d(cells,plex,ho)); 165 PetscCall(PetscFinalize()); 166 return 0; 167 } 168 169 /*TEST 170 171 testset: 172 nsize: 1 173 args: -view glvis: 174 175 test: 176 suffix: dmda_glvis_ho 177 args: -nex 1 -ney 1 -nez 1 178 179 test: 180 suffix: dmda_glvis_lo 181 args: -ho 0 182 183 test: 184 suffix: dmplex_glvis_ho 185 args: -nex 1 -ney 1 -nez 1 186 187 test: 188 suffix: dmplex_glvis_lo 189 args: -ho 0 190 191 TEST*/ 192