1 #include <petscdm.h> 2 #include <petscdmda.h> 3 #include <petscdmswarm.h> 4 #include <petsc/private/dmswarmimpl.h> 5 #include "../src/dm/impls/swarm/data_bucket.h" 6 7 PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim,PetscInt np[],PetscInt *_npoints,PetscReal **_xi) 8 { 9 PetscReal *xi; 10 PetscInt d,npoints=0,cnt; 11 PetscReal ds[] = {0.0,0.0,0.0}; 12 PetscInt ii,jj,kk; 13 14 PetscFunctionBegin; 15 switch (dim) { 16 case 1: 17 npoints = np[0]; 18 break; 19 case 2: 20 npoints = np[0]*np[1]; 21 break; 22 case 3: 23 npoints = np[0]*np[1]*np[2]; 24 break; 25 } 26 for (d=0; d<dim; d++) { 27 ds[d] = 2.0 / ((PetscReal)np[d]); 28 } 29 30 PetscCall(PetscMalloc1(dim*npoints,&xi)); 31 switch (dim) { 32 case 1: 33 cnt = 0; 34 for (ii=0; ii<np[0]; ii++) { 35 xi[dim*cnt+0] = -1.0 + 0.5*ds[d] + ii*ds[0]; 36 cnt++; 37 } 38 break; 39 40 case 2: 41 cnt = 0; 42 for (jj=0; jj<np[1]; jj++) { 43 for (ii=0; ii<np[0]; ii++) { 44 xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0]; 45 xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1]; 46 cnt++; 47 } 48 } 49 break; 50 51 case 3: 52 cnt = 0; 53 for (kk=0; kk<np[2]; kk++) { 54 for (jj=0; jj<np[1]; jj++) { 55 for (ii=0; ii<np[0]; ii++) { 56 xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0]; 57 xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1]; 58 xi[dim*cnt+2] = -1.0 + 0.5*ds[2] + kk*ds[2]; 59 cnt++; 60 } 61 } 62 } 63 break; 64 } 65 *_npoints = npoints; 66 *_xi = xi; 67 PetscFunctionReturn(0); 68 } 69 70 PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt *_npoints,PetscReal **_xi) 71 { 72 PetscQuadrature quadrature; 73 const PetscReal *quadrature_xi; 74 PetscReal *xi; 75 PetscInt d,q,npoints_q; 76 77 PetscFunctionBegin; 78 PetscCall(PetscDTGaussTensorQuadrature(dim,1,np_1d,-1.0,1.0,&quadrature)); 79 PetscCall(PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&quadrature_xi,NULL)); 80 PetscCall(PetscMalloc1(dim*npoints_q,&xi)); 81 for (q=0; q<npoints_q; q++) { 82 for (d=0; d<dim; d++) { 83 xi[dim*q+d] = quadrature_xi[dim*q+d]; 84 } 85 } 86 PetscCall(PetscQuadratureDestroy(&quadrature)); 87 *_npoints = npoints_q; 88 *_xi = xi; 89 PetscFunctionReturn(0); 90 } 91 92 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout) 93 { 94 PetscInt dim,npoints_q; 95 PetscInt nel,npe,e,q,k,d; 96 const PetscInt *element_list; 97 PetscReal **basis; 98 PetscReal *xi; 99 Vec coor; 100 const PetscScalar *_coor; 101 PetscReal *elcoor; 102 PetscReal *swarm_coor; 103 PetscInt *swarm_cellid; 104 PetscInt pcnt; 105 106 PetscFunctionBegin; 107 PetscCall(DMGetDimension(dm,&dim)); 108 switch (layout) { 109 case DMSWARMPIC_LAYOUT_REGULAR: 110 { 111 PetscInt np_dir[3]; 112 np_dir[0] = np_dir[1] = np_dir[2] = npoints; 113 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi)); 114 } 115 break; 116 case DMSWARMPIC_LAYOUT_GAUSS: 117 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim,npoints,&npoints_q,&xi)); 118 break; 119 120 case DMSWARMPIC_LAYOUT_SUBDIVISION: 121 { 122 PetscInt s,nsub; 123 PetscInt np_dir[3]; 124 nsub = npoints; 125 np_dir[0] = 1; 126 for (s=0; s<nsub; s++) { 127 np_dir[0] *= 2; 128 } 129 np_dir[1] = np_dir[0]; 130 np_dir[2] = np_dir[0]; 131 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi)); 132 } 133 break; 134 default: 135 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"A valid DMSwarmPIC layout must be provided"); 136 } 137 138 PetscCall(DMDAGetElements(dmc,&nel,&npe,&element_list)); 139 PetscCall(PetscMalloc1(dim*npe,&elcoor)); 140 PetscCall(PetscMalloc1(npoints_q,&basis)); 141 for (q=0; q<npoints_q; q++) { 142 PetscCall(PetscMalloc1(npe,&basis[q])); 143 144 switch (dim) { 145 case 1: 146 basis[q][0] = 0.5*(1.0 - xi[dim*q+0]); 147 basis[q][1] = 0.5*(1.0 + xi[dim*q+0]); 148 break; 149 case 2: 150 basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]); 151 basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]); 152 basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]); 153 basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]); 154 break; 155 156 case 3: 157 basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]); 158 basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]); 159 basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]); 160 basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]); 161 basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]); 162 basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]); 163 basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]); 164 basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]); 165 break; 166 } 167 } 168 169 PetscCall(DMSwarmSetLocalSizes(dm,npoints_q*nel,-1)); 170 PetscCall(DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor)); 171 PetscCall(DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid)); 172 173 PetscCall(DMGetCoordinatesLocal(dmc,&coor)); 174 PetscCall(VecGetArrayRead(coor,&_coor)); 175 pcnt = 0; 176 for (e=0; e<nel; e++) { 177 const PetscInt *element = &element_list[npe*e]; 178 179 for (k=0; k<npe; k++) { 180 for (d=0; d<dim; d++) { 181 elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]); 182 } 183 } 184 185 for (q=0; q<npoints_q; q++) { 186 for (d=0; d<dim; d++) { 187 swarm_coor[dim*pcnt+d] = 0.0; 188 } 189 for (k=0; k<npe; k++) { 190 for (d=0; d<dim; d++) { 191 swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d]; 192 } 193 } 194 swarm_cellid[pcnt] = e; 195 pcnt++; 196 } 197 } 198 PetscCall(VecRestoreArrayRead(coor,&_coor)); 199 PetscCall(DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid)); 200 PetscCall(DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor)); 201 PetscCall(DMDARestoreElements(dmc,&nel,&npe,&element_list)); 202 203 PetscCall(PetscFree(xi)); 204 PetscCall(PetscFree(elcoor)); 205 for (q=0; q<npoints_q; q++) { 206 PetscCall(PetscFree(basis[q])); 207 } 208 PetscCall(PetscFree(basis)); 209 PetscFunctionReturn(0); 210 } 211 212 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 213 { 214 DMDAElementType etype; 215 PetscInt dim; 216 217 PetscFunctionBegin; 218 PetscCall(DMDAGetElementType(celldm,&etype)); 219 PetscCall(DMGetDimension(celldm,&dim)); 220 switch (etype) { 221 case DMDA_ELEMENT_P1: 222 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DA support is not currently available for DMDA_ELEMENT_P1"); 223 case DMDA_ELEMENT_Q1: 224 PetscCheckFalse(dim == 1,PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only available for dim = 2, 3"); 225 PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout)); 226 break; 227 } 228 PetscFunctionReturn(0); 229 } 230 231 PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field) 232 { 233 Vec v_field_l,denom_l,coor_l,denom; 234 PetscScalar *_field_l,*_denom_l; 235 PetscInt k,p,e,npoints,nel,npe; 236 PetscInt *mpfield_cell; 237 PetscReal *mpfield_coor; 238 const PetscInt *element_list; 239 const PetscInt *element; 240 PetscScalar xi_p[2],Ni[4]; 241 const PetscScalar *_coor; 242 243 PetscFunctionBegin; 244 PetscCall(VecZeroEntries(v_field)); 245 246 PetscCall(DMGetLocalVector(dm,&v_field_l)); 247 PetscCall(DMGetGlobalVector(dm,&denom)); 248 PetscCall(DMGetLocalVector(dm,&denom_l)); 249 PetscCall(VecZeroEntries(v_field_l)); 250 PetscCall(VecZeroEntries(denom)); 251 PetscCall(VecZeroEntries(denom_l)); 252 253 PetscCall(VecGetArray(v_field_l,&_field_l)); 254 PetscCall(VecGetArray(denom_l,&_denom_l)); 255 256 PetscCall(DMGetCoordinatesLocal(dm,&coor_l)); 257 PetscCall(VecGetArrayRead(coor_l,&_coor)); 258 259 PetscCall(DMDAGetElements(dm,&nel,&npe,&element_list)); 260 PetscCall(DMSwarmGetLocalSize(swarm,&npoints)); 261 PetscCall(DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor)); 262 PetscCall(DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell)); 263 264 for (p=0; p<npoints; p++) { 265 PetscReal *coor_p; 266 const PetscScalar *x0; 267 const PetscScalar *x2; 268 PetscScalar dx[2]; 269 270 e = mpfield_cell[p]; 271 coor_p = &mpfield_coor[2*p]; 272 element = &element_list[npe*e]; 273 274 /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */ 275 x0 = &_coor[2*element[0]]; 276 x2 = &_coor[2*element[2]]; 277 278 dx[0] = x2[0] - x0[0]; 279 dx[1] = x2[1] - x0[1]; 280 281 xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0; 282 xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0; 283 284 /* evaluate basis functions */ 285 Ni[0] = 0.25*(1.0 - xi_p[0])*(1.0 - xi_p[1]); 286 Ni[1] = 0.25*(1.0 + xi_p[0])*(1.0 - xi_p[1]); 287 Ni[2] = 0.25*(1.0 + xi_p[0])*(1.0 + xi_p[1]); 288 Ni[3] = 0.25*(1.0 - xi_p[0])*(1.0 + xi_p[1]); 289 290 for (k=0; k<npe; k++) { 291 _field_l[ element[k] ] += Ni[k] * swarm_field[p]; 292 _denom_l[ element[k] ] += Ni[k]; 293 } 294 } 295 296 PetscCall(DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell)); 297 PetscCall(DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor)); 298 PetscCall(DMDARestoreElements(dm,&nel,&npe,&element_list)); 299 PetscCall(VecRestoreArrayRead(coor_l,&_coor)); 300 PetscCall(VecRestoreArray(v_field_l,&_field_l)); 301 PetscCall(VecRestoreArray(denom_l,&_denom_l)); 302 303 PetscCall(DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field)); 304 PetscCall(DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field)); 305 PetscCall(DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom)); 306 PetscCall(DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom)); 307 308 PetscCall(VecPointwiseDivide(v_field,v_field,denom)); 309 310 PetscCall(DMRestoreLocalVector(dm,&v_field_l)); 311 PetscCall(DMRestoreLocalVector(dm,&denom_l)); 312 PetscCall(DMRestoreGlobalVector(dm,&denom)); 313 PetscFunctionReturn(0); 314 } 315 316 PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]) 317 { 318 PetscInt f,dim; 319 DMDAElementType etype; 320 321 PetscFunctionBegin; 322 PetscCall(DMDAGetElementType(celldm,&etype)); 323 PetscCheckFalse(etype == DMDA_ELEMENT_P1,PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"Only Q1 DMDA supported"); 324 325 PetscCall(DMGetDimension(swarm,&dim)); 326 switch (dim) { 327 case 2: 328 for (f=0; f<nfields; f++) { 329 PetscReal *swarm_field; 330 331 PetscCall(DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field)); 332 PetscCall(DMSwarmProjectField_ApproxQ1_DA_2D(swarm,swarm_field,celldm,vecs[f])); 333 } 334 break; 335 case 3: 336 SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D"); 337 default: 338 break; 339 } 340 PetscFunctionReturn(0); 341 } 342