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,PetscInt 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 default: 128 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"A valid DMSwarmPIC layout must be provided"); 129 break; 130 } 131 132 ierr = DMDAGetElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr); 133 134 ierr = PetscMalloc1(dim*npe,&elcoor);CHKERRQ(ierr); 135 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 136 for (q=0; q<npoints_q; q++) { 137 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 138 139 switch (dim) { 140 case 1: 141 basis[q][0] = 0.5*(1.0 - xi[dim*q+0]); 142 basis[q][1] = 0.5*(1.0 + xi[dim*q+0]); 143 break; 144 case 2: 145 basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]); 146 basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]); 147 basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]); 148 basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]); 149 break; 150 151 case 3: 152 basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]); 153 basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]); 154 basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]); 155 basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]); 156 basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]); 157 basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]); 158 basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]); 159 basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]); 160 break; 161 } 162 } 163 164 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 165 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 166 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 167 168 ierr = DMGetCoordinatesLocal(dmc,&coor);CHKERRQ(ierr); 169 ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr); 170 pcnt = 0; 171 for (e=0; e<nel; e++) { 172 const PetscInt *element = &element_list[npe*e]; 173 174 for (k=0; k<npe; k++) { 175 for (d=0; d<dim; d++) { 176 elcoor[dim*k+d] = _coor[ dim*element[k] + d ]; 177 } 178 } 179 180 for (q=0; q<npoints_q; q++) { 181 for (d=0; d<dim; d++) { 182 swarm_coor[dim*pcnt+d] = 0.0; 183 } 184 for (k=0; k<npe; k++) { 185 for (d=0; d<dim; d++) { 186 swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d]; 187 } 188 } 189 swarm_cellid[pcnt] = e; 190 pcnt++; 191 } 192 } 193 ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr); 194 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 195 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 196 ierr = DMDARestoreElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr); 197 198 ierr = PetscFree(xi);CHKERRQ(ierr); 199 ierr = PetscFree(elcoor);CHKERRQ(ierr); 200 for (q=0; q<npoints_q; q++) { 201 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 202 } 203 ierr = PetscFree(basis);CHKERRQ(ierr); 204 205 PetscFunctionReturn(0); 206 } 207 208 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 209 { 210 PetscErrorCode ierr; 211 DMDAElementType etype; 212 PetscInt dim; 213 214 PetscFunctionBegin; 215 ierr = DMDAGetElementType(celldm,&etype);CHKERRQ(ierr); 216 ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr); 217 switch (etype) { 218 case DMDA_ELEMENT_P1: 219 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DA support is not currently available for DMDA_ELEMENT_P1"); 220 break; 221 case DMDA_ELEMENT_Q1: 222 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only available for dim = 2, 3"); 223 ierr = private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);CHKERRQ(ierr); 224 break; 225 } 226 PetscFunctionReturn(0); 227 } 228