1 2 #include <petsc.h> 3 #include <petscdm.h> 4 #include <petscdmplex.h> 5 #include <petscdmswarm.h> 6 7 8 PetscErrorCode subdivide_triangle(PetscReal v1[2],PetscReal v2[2],PetscReal v3[3],PetscInt depth,PetscInt max,PetscReal xi[],PetscInt *np) 9 { 10 PetscReal v12[2],v23[2],v31[2]; 11 PetscInt i; 12 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 if (depth == max) { 16 PetscReal cx[2]; 17 18 //printf("centroid \n"); 19 cx[0] = (v1[0] + v2[0] + v3[0])/3.0; 20 cx[1] = (v1[1] + v2[1] + v3[1])/3.0; 21 22 xi[2*(*np)+0] = cx[0]; 23 xi[2*(*np)+1] = cx[1]; 24 (*np)++; 25 PetscFunctionReturn(0); 26 } 27 28 /* calculate midpoints of each side */ 29 for (i = 0; i < 2; i++) { 30 v12[i] = (v1[i]+v2[i])/2.0; 31 v23[i] = (v2[i]+v3[i])/2.0; 32 v31[i] = (v3[i]+v1[i])/2.0; 33 } 34 35 /* recursively subdivide new triangles */ 36 ierr = subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);CHKERRQ(ierr); 37 ierr = subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);CHKERRQ(ierr); 38 ierr = subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);CHKERRQ(ierr); 39 ierr = subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 44 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm,DM dmc,PetscInt nsub) 45 { 46 PetscErrorCode ierr; 47 const PetscInt dim = 2; 48 PetscInt q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,depth; 49 PetscReal *xi; 50 PetscReal **basis; 51 Vec coorlocal; 52 PetscSection coordSection; 53 PetscScalar *elcoor = NULL; 54 PetscReal *swarm_coor; 55 PetscInt *swarm_cellid; 56 PetscReal v1[2],v2[2],v3[2]; 57 58 PetscFunctionBegin; 59 60 npoints_q = 1; 61 for (d=0; d<nsub; d++) { npoints_q *= 4; } 62 ierr = PetscMalloc1(dim*npoints_q,&xi);CHKERRQ(ierr); 63 64 v1[0] = 0.0; v1[1] = 0.0; 65 v2[0] = 1.0; v2[1] = 0.0; 66 v3[0] = 0.0; v3[1] = 1.0; 67 depth = 0; 68 pcnt = 0; 69 ierr = subdivide_triangle(v1,v2,v3,depth,nsub,xi,&pcnt);CHKERRQ(ierr); 70 71 npe = 3; /* nodes per element (triangle) */ 72 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 73 for (q=0; q<npoints_q; q++) { 74 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 75 76 basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 77 basis[q][1] = xi[dim*q+0]; 78 basis[q][2] = xi[dim*q+1]; 79 } 80 81 /* 0->cell, 1->edge, 2->vert */ 82 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 83 nel = pe - ps; 84 85 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 86 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 87 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 88 89 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 90 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 91 92 pcnt = 0; 93 for (e=0; e<nel; e++) { 94 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 95 96 for (q=0; q<npoints_q; q++) { 97 for (d=0; d<dim; d++) { 98 swarm_coor[dim*pcnt+d] = 0.0; 99 for (k=0; k<npe; k++) { 100 swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d]; 101 } 102 } 103 swarm_cellid[pcnt] = e; 104 pcnt++; 105 } 106 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 107 } 108 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 109 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 110 111 ierr = PetscFree(xi);CHKERRQ(ierr); 112 for (q=0; q<npoints_q; q++) { 113 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 114 } 115 ierr = PetscFree(basis);CHKERRQ(ierr); 116 117 PetscFunctionReturn(0); 118 } 119 120 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints) 121 { 122 PetscErrorCode ierr; 123 const PetscInt dim = 2; 124 PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k; 125 PetscReal *xi,ds,ds2; 126 PetscReal **basis; 127 Vec coorlocal; 128 PetscSection coordSection; 129 PetscScalar *elcoor = NULL; 130 PetscReal *swarm_coor; 131 PetscInt *swarm_cellid; 132 133 PetscFunctionBegin; 134 ierr = PetscMalloc1(dim*npoints*npoints,&xi);CHKERRQ(ierr); 135 pcnt = 0; 136 ds = 1.0/((PetscReal)(npoints-1)); 137 ds2 = 1.0/((PetscReal)(npoints)); 138 for (jj = 0; jj<npoints; jj++) { 139 for (ii=0; ii<npoints-jj; ii++) { 140 xi[dim*pcnt+0] = ii * ds; 141 xi[dim*pcnt+1] = jj * ds; 142 143 xi[dim*pcnt+0] *= (1.0 - 1.2*ds2); 144 xi[dim*pcnt+1] *= (1.0 - 1.2*ds2); 145 146 xi[dim*pcnt+0] += 0.35*ds2; 147 xi[dim*pcnt+1] += 0.35*ds2; 148 149 pcnt++; 150 } 151 } 152 npoints_q = pcnt; 153 154 npe = 3; /* nodes per element (triangle) */ 155 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 156 for (q=0; q<npoints_q; q++) { 157 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 158 159 basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 160 basis[q][1] = xi[dim*q+0]; 161 basis[q][2] = xi[dim*q+1]; 162 } 163 164 /* 0->cell, 1->edge, 2->vert */ 165 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 166 nel = pe - ps; 167 168 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 169 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 170 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 171 172 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 173 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 174 175 pcnt = 0; 176 for (e=0; e<nel; e++) { 177 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 178 179 for (q=0; q<npoints_q; q++) { 180 for (d=0; d<dim; d++) { 181 swarm_coor[dim*pcnt+d] = 0.0; 182 for (k=0; k<npe; k++) { 183 swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d]; 184 } 185 } 186 swarm_cellid[pcnt] = e; 187 pcnt++; 188 } 189 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 190 } 191 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 192 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 193 194 ierr = PetscFree(xi);CHKERRQ(ierr); 195 for (q=0; q<npoints_q; q++) { 196 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 197 } 198 ierr = PetscFree(basis);CHKERRQ(ierr); 199 200 PetscFunctionReturn(0); 201 } 202 203 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 204 { 205 PetscErrorCode ierr; 206 PetscInt dim; 207 208 PetscFunctionBegin; 209 ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr); 210 if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No 3D support for PLEX"); 211 switch (layout) { 212 case DMSWARMPIC_LAYOUT_REGULAR: 213 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);CHKERRQ(ierr); 214 break; 215 case DMSWARMPIC_LAYOUT_GAUSS: 216 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Gauss layout not supported for PLEX"); 217 break; 218 case DMSWARMPIC_LAYOUT_SUBDIVISION: 219 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(dm,celldm,layout_param);CHKERRQ(ierr); 220 break; 221 } 222 PetscFunctionReturn(0); 223 } 224