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