1 2 #include <petsc.h> 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 #include <petscdmswarm.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,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 33 switch (dim) { 34 case 1: 35 cnt = 0; 36 for (ii=0; ii<np[0]; ii++) { 37 xi[dim*cnt+0] = -1.0 + 0.5*ds[d] + ii*ds[0]; 38 cnt++; 39 } 40 break; 41 42 case 2: 43 cnt = 0; 44 for (jj=0; jj<np[1]; jj++) { 45 for (ii=0; ii<np[0]; ii++) { 46 xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0]; 47 xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1]; 48 cnt++; 49 } 50 } 51 break; 52 53 case 3: 54 cnt = 0; 55 for (kk=0; kk<np[2]; kk++) { 56 for (jj=0; jj<np[1]; jj++) { 57 for (ii=0; ii<np[0]; ii++) { 58 xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0]; 59 xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1]; 60 xi[dim*cnt+2] = -1.0 + 0.5*ds[2] + kk*ds[2]; 61 cnt++; 62 } 63 } 64 } 65 break; 66 } 67 68 *_npoints = npoints; 69 *_xi = xi; 70 71 PetscFunctionReturn(0); 72 } 73 74 PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt *_npoints,PetscReal **_xi) 75 { 76 PetscErrorCode ierr; 77 PetscQuadrature quadrature; 78 const PetscReal *quadrature_xi; 79 PetscReal *xi; 80 PetscInt d,q,npoints_q; 81 82 PetscFunctionBegin; 83 ierr = PetscDTGaussTensorQuadrature(dim,np_1d,-1.0,1.0,&quadrature);CHKERRQ(ierr); 84 ierr = PetscQuadratureGetData(quadrature,NULL,&npoints_q,&quadrature_xi,NULL);CHKERRQ(ierr); 85 ierr = PetscMalloc1(dim*npoints_q,&xi);CHKERRQ(ierr); 86 for (q=0; q<npoints_q; q++) { 87 for (d=0; d<dim; d++) { 88 xi[dim*q+d] = quadrature_xi[dim*q+d]; 89 } 90 } 91 ierr = PetscQuadratureDestroy(&quadrature);CHKERRQ(ierr); 92 93 *_npoints = npoints_q; 94 *_xi = xi; 95 96 PetscFunctionReturn(0); 97 } 98 99 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout) 100 { 101 PetscErrorCode ierr; 102 PetscInt dim,npoints_q; 103 PetscInt nel,npe,e,q,k,d; 104 const PetscInt *element_list; 105 PetscReal **basis; 106 PetscReal *xi; 107 Vec coor; 108 const PetscScalar *_coor; 109 PetscReal *elcoor; 110 PetscReal *swarm_coor; 111 PetscInt *swarm_cellid; 112 PetscInt pcnt; 113 114 PetscFunctionBegin; 115 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 116 switch (layout) { 117 case DMSWARMPIC_LAYOUT_REGULAR: 118 { 119 PetscInt np_dir[3]; 120 np_dir[0] = np_dir[1] = np_dir[2] = npoints; 121 ierr = private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);CHKERRQ(ierr); 122 } 123 break; 124 case DMSWARMPIC_LAYOUT_GAUSS: 125 ierr = private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim,npoints,&npoints_q,&xi);CHKERRQ(ierr); 126 break; 127 128 case DMSWARMPIC_LAYOUT_SUBDIVISION: 129 { 130 PetscInt s,nsub; 131 PetscInt np_dir[3]; 132 nsub = npoints; 133 np_dir[0] = 1; 134 for (s=0; s<nsub; s++) { 135 np_dir[0] *= 2; 136 } 137 np_dir[1] = np_dir[0]; 138 np_dir[2] = np_dir[0]; 139 ierr = private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);CHKERRQ(ierr); 140 } 141 break; 142 default: 143 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"A valid DMSwarmPIC layout must be provided"); 144 break; 145 } 146 147 ierr = DMDAGetElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr); 148 149 ierr = PetscMalloc1(dim*npe,&elcoor);CHKERRQ(ierr); 150 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 151 for (q=0; q<npoints_q; q++) { 152 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 153 154 switch (dim) { 155 case 1: 156 basis[q][0] = 0.5*(1.0 - xi[dim*q+0]); 157 basis[q][1] = 0.5*(1.0 + xi[dim*q+0]); 158 break; 159 case 2: 160 basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]); 161 basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]); 162 basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]); 163 basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]); 164 break; 165 166 case 3: 167 basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]); 168 basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]); 169 basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]); 170 basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]); 171 basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]); 172 basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]); 173 basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]); 174 basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]); 175 break; 176 } 177 } 178 179 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 180 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 181 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 182 183 ierr = DMGetCoordinatesLocal(dmc,&coor);CHKERRQ(ierr); 184 ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr); 185 pcnt = 0; 186 for (e=0; e<nel; e++) { 187 const PetscInt *element = &element_list[npe*e]; 188 189 for (k=0; k<npe; k++) { 190 for (d=0; d<dim; d++) { 191 elcoor[dim*k+d] = _coor[ dim*element[k] + d ]; 192 } 193 } 194 195 for (q=0; q<npoints_q; q++) { 196 for (d=0; d<dim; d++) { 197 swarm_coor[dim*pcnt+d] = 0.0; 198 } 199 for (k=0; k<npe; k++) { 200 for (d=0; d<dim; d++) { 201 swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d]; 202 } 203 } 204 swarm_cellid[pcnt] = e; 205 pcnt++; 206 } 207 } 208 ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr); 209 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 210 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 211 ierr = DMDARestoreElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr); 212 213 ierr = PetscFree(xi);CHKERRQ(ierr); 214 ierr = PetscFree(elcoor);CHKERRQ(ierr); 215 for (q=0; q<npoints_q; q++) { 216 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 217 } 218 ierr = PetscFree(basis);CHKERRQ(ierr); 219 220 PetscFunctionReturn(0); 221 } 222 223 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 224 { 225 PetscErrorCode ierr; 226 DMDAElementType etype; 227 PetscInt dim; 228 229 PetscFunctionBegin; 230 ierr = DMDAGetElementType(celldm,&etype);CHKERRQ(ierr); 231 ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr); 232 switch (etype) { 233 case DMDA_ELEMENT_P1: 234 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DA support is not currently available for DMDA_ELEMENT_P1"); 235 break; 236 case DMDA_ELEMENT_Q1: 237 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only available for dim = 2, 3"); 238 ierr = private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);CHKERRQ(ierr); 239 break; 240 } 241 PetscFunctionReturn(0); 242 } 243