xref: /petsc/src/dm/tutorials/ex20.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce) 
1 
2 static char help[] = "DMSwarm-PIC demonstator of inserting points into a cell DM \n\
3 Options: \n\
4 -mode {0,1} : 0 ==> DMDA, 1 ==> DMPLEX cell DM \n\
5 -dim {2,3}  : spatial dimension\n";
6 
7 #include <petsc.h>
8 #include <petscdm.h>
9 #include <petscdmda.h>
10 #include <petscdmplex.h>
11 #include <petscdmswarm.h>
12 
13 PetscErrorCode pic_insert_DMDA(PetscInt dim)
14 {
15   PetscErrorCode ierr;
16   DM             celldm = NULL,swarm;
17   PetscInt       dof,stencil_width;
18   PetscReal      min[3],max[3];
19   PetscInt       ndir[3];
20 
21   PetscFunctionBegin;
22   /* Create the background cell DM */
23   dof = 1;
24   stencil_width = 1;
25   if (dim == 2) {
26     ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,25,13,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&celldm);CHKERRQ(ierr);
27   }
28   if (dim == 3) {
29     ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,25,13,19,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&celldm);CHKERRQ(ierr);
30   }
31 
32   ierr = DMDASetElementType(celldm,DMDA_ELEMENT_Q1);CHKERRQ(ierr);
33   ierr = DMSetFromOptions(celldm);CHKERRQ(ierr);
34   ierr = DMSetUp(celldm);CHKERRQ(ierr);
35 
36   ierr = DMDASetUniformCoordinates(celldm,0.0,2.0,0.0,1.0,0.0,1.5);CHKERRQ(ierr);
37 
38   /* Create the DMSwarm */
39   ierr = DMCreate(PETSC_COMM_WORLD,&swarm);CHKERRQ(ierr);
40   ierr = DMSetType(swarm,DMSWARM);CHKERRQ(ierr);
41   ierr = DMSetDimension(swarm,dim);CHKERRQ(ierr);
42 
43   /* Configure swarm to be of type PIC */
44   ierr = DMSwarmSetType(swarm,DMSWARM_PIC);CHKERRQ(ierr);
45   ierr = DMSwarmSetCellDM(swarm,celldm);CHKERRQ(ierr);
46 
47   /* Register two scalar fields within the DMSwarm */
48   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);CHKERRQ(ierr);
49   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);CHKERRQ(ierr);
50   ierr = DMSwarmFinalizeFieldRegister(swarm);CHKERRQ(ierr);
51 
52   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
53   ierr = DMSwarmSetLocalSizes(swarm,4,0);CHKERRQ(ierr);
54 
55   /* Insert swarm coordinates cell-wise */
56   ierr = DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,3);CHKERRQ(ierr);
57   min[0] = 0.5; max[0] = 0.7;
58   min[1] = 0.5; max[1] = 0.8;
59   min[2] = 0.5; max[2] = 0.9;
60   ndir[0] = ndir[1] = ndir[2] = 30;
61   ierr = DMSwarmSetPointsUniformCoordinates(swarm,min,max,ndir,ADD_VALUES);CHKERRQ(ierr);
62 
63   /* This should be dispatched from a regular DMView() */
64   ierr = DMSwarmViewXDMF(swarm,"ex20.xmf");CHKERRQ(ierr);
65   ierr = DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
66   ierr = DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
67 
68   {
69     PetscInt    npoints,*list;
70     PetscMPIInt rank;
71 
72     ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
73     ierr = DMSwarmSortGetAccess(swarm);CHKERRQ(ierr);
74     ierr = DMSwarmSortGetNumberOfPointsPerCell(swarm,0,&npoints);CHKERRQ(ierr);
75     ierr = DMSwarmSortGetPointsPerCell(swarm,rank,&npoints,&list);CHKERRQ(ierr);
76     ierr = PetscFree(list);CHKERRQ(ierr);
77     ierr = DMSwarmSortRestoreAccess(swarm);CHKERRQ(ierr);
78   }
79   ierr = DMSwarmMigrate(swarm,PETSC_FALSE);CHKERRQ(ierr);
80   ierr = DMDestroy(&celldm);CHKERRQ(ierr);
81   ierr = DMDestroy(&swarm);CHKERRQ(ierr);
82   PetscFunctionReturn(0);
83 }
84 
85 PetscErrorCode pic_insert_DMPLEX_with_cell_list(PetscInt dim)
86 {
87   PetscErrorCode ierr;
88   DM             celldm = NULL,swarm,distributedMesh = NULL;
89   const  char    *fieldnames[] = {"viscosity"};
90 
91   PetscFunctionBegin;
92   /* Create the background cell DM */
93   if (dim == 2) {
94     PetscInt   cells_per_dim[2],nx[2];
95     PetscInt   n_tricells;
96     PetscInt   n_trivert;
97     PetscInt   *tricells;
98     PetscReal  *trivert,dx,dy;
99     PetscInt   ii,jj,cnt;
100 
101     cells_per_dim[0] = 4;
102     cells_per_dim[1] = 4;
103     n_tricells = cells_per_dim[0] * cells_per_dim[1] * 2;
104     nx[0] = cells_per_dim[0] + 1;
105     nx[1] = cells_per_dim[1] + 1;
106     n_trivert = nx[0] * nx[1];
107 
108     ierr = PetscMalloc1(n_tricells*3,&tricells);CHKERRQ(ierr);
109     ierr = PetscMalloc1(nx[0]*nx[1]*2,&trivert);CHKERRQ(ierr);
110 
111     /* verts */
112     cnt = 0;
113     dx = 2.0/((PetscReal)cells_per_dim[0]);
114     dy = 1.0/((PetscReal)cells_per_dim[1]);
115     for (jj=0; jj<nx[1]; jj++) {
116       for (ii=0; ii<nx[0]; ii++) {
117         trivert[2*cnt+0] = 0.0 + ii * dx;
118         trivert[2*cnt+1] = 0.0 + jj * dy;
119         cnt++;
120       }
121     }
122 
123     /* connectivity */
124     cnt = 0;
125     for (jj=0; jj<cells_per_dim[1]; jj++) {
126       for (ii=0; ii<cells_per_dim[0]; ii++) {
127         PetscInt idx,idx0,idx1,idx2,idx3;
128 
129         idx = (ii) + (jj) * nx[0];
130         idx0 = idx;
131         idx1 = idx0 + 1;
132         idx2 = idx1 + nx[0];
133         idx3 = idx0 + nx[0];
134 
135         tricells[3*cnt+0] = idx0;
136         tricells[3*cnt+1] = idx1;
137         tricells[3*cnt+2] = idx2;
138         cnt++;
139 
140         tricells[3*cnt+0] = idx0;
141         tricells[3*cnt+1] = idx2;
142         tricells[3*cnt+2] = idx3;
143         cnt++;
144       }
145     }
146     ierr = DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD,dim,n_tricells,n_trivert,3,PETSC_TRUE,tricells,dim,trivert,&celldm);CHKERRQ(ierr);
147     ierr = PetscFree(trivert);CHKERRQ(ierr);
148     ierr = PetscFree(tricells);CHKERRQ(ierr);
149   }
150   if (dim == 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only 2D PLEX example supported");
151 
152   /* Distribute mesh over processes */
153   ierr = DMPlexDistribute(celldm,0,NULL,&distributedMesh);CHKERRQ(ierr);
154   if (distributedMesh) {
155     ierr = DMDestroy(&celldm);CHKERRQ(ierr);
156     celldm = distributedMesh;
157   }
158   ierr = DMSetFromOptions(celldm);CHKERRQ(ierr);
159   {
160     PetscInt     numComp[] = {1};
161     PetscInt     numDof[] = {1,0,0}; /* vert, edge, cell */
162     PetscInt     numBC = 0;
163     PetscSection section;
164 
165     ierr = DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,&section);CHKERRQ(ierr);
166     ierr = DMSetLocalSection(celldm,section);CHKERRQ(ierr);
167     ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
168   }
169   ierr = DMSetUp(celldm);CHKERRQ(ierr);
170   {
171     PetscViewer viewer;
172 
173     ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr);
174     ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
175     ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
176     ierr = PetscViewerFileSetName(viewer,"ex20plex.vtk");CHKERRQ(ierr);
177     ierr = DMView(celldm,viewer);CHKERRQ(ierr);
178     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
179   }
180 
181   /* Create the DMSwarm */
182   ierr = DMCreate(PETSC_COMM_WORLD,&swarm);CHKERRQ(ierr);
183   ierr = DMSetType(swarm,DMSWARM);CHKERRQ(ierr);
184   ierr = DMSetDimension(swarm,dim);CHKERRQ(ierr);
185 
186   ierr = DMSwarmSetType(swarm,DMSWARM_PIC);CHKERRQ(ierr);
187   ierr = DMSwarmSetCellDM(swarm,celldm);CHKERRQ(ierr);
188 
189   /* Register two scalar fields within the DMSwarm */
190   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);CHKERRQ(ierr);
191   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);CHKERRQ(ierr);
192   ierr = DMSwarmFinalizeFieldRegister(swarm);CHKERRQ(ierr);
193 
194   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
195   ierr = DMSwarmSetLocalSizes(swarm,4,0);CHKERRQ(ierr);
196 
197   /* Insert swarm coordinates cell-wise */
198   ierr = DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_SUBDIVISION,2);CHKERRQ(ierr);
199   ierr = DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",1,fieldnames);CHKERRQ(ierr);
200   ierr = DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
201   ierr = DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
202   ierr = DMDestroy(&celldm);CHKERRQ(ierr);
203   ierr = DMDestroy(&swarm);CHKERRQ(ierr);
204   PetscFunctionReturn(0);
205 }
206 
207 PetscErrorCode pic_insert_DMPLEX(PetscBool is_simplex,PetscInt dim)
208 {
209   PetscErrorCode ierr;
210   DM             celldm,swarm,distributedMesh = NULL;
211   const char     *fieldnames[] = {"viscosity","DMSwarm_rank"};
212 
213   PetscFunctionBegin;
214 
215   /* Create the background cell DM */
216   {
217     PetscInt faces[3] = {4, 2, 4};
218     ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, &celldm);CHKERRQ(ierr);
219   }
220 
221   /* Distribute mesh over processes */
222   ierr = DMPlexDistribute(celldm,0,NULL,&distributedMesh);CHKERRQ(ierr);
223   if (distributedMesh) {
224     ierr = DMDestroy(&celldm);CHKERRQ(ierr);
225     celldm = distributedMesh;
226   }
227   ierr = DMSetFromOptions(celldm);CHKERRQ(ierr);
228   {
229     PetscInt     numComp[] = {1};
230     PetscInt     numDof[] = {1,0,0}; /* vert, edge, cell */
231     PetscInt     numBC = 0;
232     PetscSection section;
233 
234     ierr = DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,&section);CHKERRQ(ierr);
235     ierr = DMSetLocalSection(celldm,section);CHKERRQ(ierr);
236     ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
237   }
238   ierr = DMSetUp(celldm);CHKERRQ(ierr);
239   {
240     PetscViewer viewer;
241 
242     ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr);
243     ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
244     ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
245     ierr = PetscViewerFileSetName(viewer,"ex20plex.vtk");CHKERRQ(ierr);
246     ierr = DMView(celldm,viewer);CHKERRQ(ierr);
247     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
248   }
249 
250   ierr = DMCreate(PETSC_COMM_WORLD,&swarm);CHKERRQ(ierr);
251   ierr = DMSetType(swarm,DMSWARM);CHKERRQ(ierr);
252   ierr = DMSetDimension(swarm,dim);CHKERRQ(ierr);
253 
254   ierr = DMSwarmSetType(swarm,DMSWARM_PIC);CHKERRQ(ierr);
255   ierr = DMSwarmSetCellDM(swarm,celldm);CHKERRQ(ierr);
256 
257   /* Register two scalar fields within the DMSwarm */
258   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);CHKERRQ(ierr);
259   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);CHKERRQ(ierr);
260   ierr = DMSwarmFinalizeFieldRegister(swarm);CHKERRQ(ierr);
261 
262   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
263   ierr = DMSwarmSetLocalSizes(swarm,4,0);CHKERRQ(ierr);
264 
265   /* Insert swarm coordinates cell-wise */
266   ierr = DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_GAUSS,3);CHKERRQ(ierr);
267   ierr = DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",2,fieldnames);CHKERRQ(ierr);
268   ierr = DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
269   ierr = DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
270   ierr = DMDestroy(&celldm);CHKERRQ(ierr);
271   ierr = DMDestroy(&swarm);CHKERRQ(ierr);
272   PetscFunctionReturn(0);
273 }
274 
275 int main(int argc,char **args)
276 {
277   PetscErrorCode ierr;
278   PetscInt       mode = 0;
279   PetscInt       dim = 2;
280 
281   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
282   ierr = PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL);CHKERRQ(ierr);
283   ierr = PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);CHKERRQ(ierr);
284   switch (mode) {
285     case 0:
286       ierr = pic_insert_DMDA(dim);CHKERRQ(ierr);
287       break;
288     case 1:
289       /* tri / tet */
290       ierr = pic_insert_DMPLEX(PETSC_TRUE,dim);CHKERRQ(ierr);
291       break;
292     case 2:
293       /* quad / hex */
294       ierr = pic_insert_DMPLEX(PETSC_FALSE,dim);CHKERRQ(ierr);
295       break;
296     default:
297       ierr = pic_insert_DMDA(dim);CHKERRQ(ierr);
298       break;
299   }
300   ierr = PetscFinalize();
301   return ierr;
302 }
303 
304 /*TEST
305 
306    test:
307       args:
308       requires: !complex double
309       filter: grep -v DM_ | grep -v atomic
310       filter_output: grep -v atomic
311 
312    test:
313       suffix: 2
314       requires: triangle double !complex
315       args: -mode 1
316       filter: grep -v DM_ | grep -v atomic
317       filter_output: grep -v atomic
318 
319 TEST*/
320