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