1 #include <petscdm.h> 2 #include <petscdmda.h> /*I "petscdmda.h" I*/ 3 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 4 #include "../src/dm/impls/swarm/data_bucket.h" 5 6 static PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim, PetscInt np[], PetscInt *_npoints, PetscReal **_xi) 7 { 8 PetscReal *xi; 9 PetscInt d, npoints = 0, cnt; 10 PetscReal ds[] = {0.0, 0.0, 0.0}; 11 PetscInt ii, jj, kk; 12 13 PetscFunctionBegin; 14 switch (dim) { 15 case 1: 16 npoints = np[0]; 17 break; 18 case 2: 19 npoints = np[0] * np[1]; 20 break; 21 case 3: 22 npoints = np[0] * np[1] * np[2]; 23 break; 24 } 25 for (d = 0; d < dim; d++) ds[d] = 2.0 / ((PetscReal)np[d]); 26 27 PetscCall(PetscMalloc1(dim * npoints, &xi)); 28 switch (dim) { 29 case 1: 30 cnt = 0; 31 for (ii = 0; ii < np[0]; ii++) { 32 xi[dim * cnt + 0] = -1.0 + 0.5 * ds[d] + ii * ds[0]; 33 cnt++; 34 } 35 break; 36 37 case 2: 38 cnt = 0; 39 for (jj = 0; jj < np[1]; jj++) { 40 for (ii = 0; ii < np[0]; ii++) { 41 xi[dim * cnt + 0] = -1.0 + 0.5 * ds[0] + ii * ds[0]; 42 xi[dim * cnt + 1] = -1.0 + 0.5 * ds[1] + jj * ds[1]; 43 cnt++; 44 } 45 } 46 break; 47 48 case 3: 49 cnt = 0; 50 for (kk = 0; kk < np[2]; kk++) { 51 for (jj = 0; jj < np[1]; jj++) { 52 for (ii = 0; ii < np[0]; ii++) { 53 xi[dim * cnt + 0] = -1.0 + 0.5 * ds[0] + ii * ds[0]; 54 xi[dim * cnt + 1] = -1.0 + 0.5 * ds[1] + jj * ds[1]; 55 xi[dim * cnt + 2] = -1.0 + 0.5 * ds[2] + kk * ds[2]; 56 cnt++; 57 } 58 } 59 } 60 break; 61 } 62 *_npoints = npoints; 63 *_xi = xi; 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 static PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim, PetscInt np_1d, PetscInt *_npoints, PetscReal **_xi) 68 { 69 PetscQuadrature quadrature; 70 const PetscReal *quadrature_xi; 71 PetscReal *xi; 72 PetscInt d, q, npoints_q; 73 74 PetscFunctionBegin; 75 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, np_1d, -1.0, 1.0, &quadrature)); 76 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &quadrature_xi, NULL)); 77 PetscCall(PetscMalloc1(dim * npoints_q, &xi)); 78 for (q = 0; q < npoints_q; q++) { 79 for (d = 0; d < dim; d++) xi[dim * q + d] = quadrature_xi[dim * q + d]; 80 } 81 PetscCall(PetscQuadratureDestroy(&quadrature)); 82 *_npoints = npoints_q; 83 *_xi = xi; 84 PetscFunctionReturn(PETSC_SUCCESS); 85 } 86 87 static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm, DM dmc, PetscInt npoints, DMSwarmPICLayoutType layout) 88 { 89 DMSwarmCellDM celldm; 90 PetscInt dim, npoints_q; 91 PetscInt nel, npe, e, q, k, d; 92 const PetscInt *element_list; 93 PetscReal **basis; 94 PetscReal *xi; 95 Vec coor; 96 const PetscScalar *_coor; 97 PetscReal *elcoor; 98 PetscReal *swarm_coor; 99 PetscInt *swarm_cellid; 100 PetscInt pcnt, Nfc; 101 const char **coordFields, *cellid; 102 103 PetscFunctionBegin; 104 PetscCall(DMGetDimension(dm, &dim)); 105 switch (layout) { 106 case DMSWARMPIC_LAYOUT_REGULAR: { 107 PetscInt np_dir[3]; 108 np_dir[0] = np_dir[1] = np_dir[2] = npoints; 109 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi)); 110 } break; 111 case DMSWARMPIC_LAYOUT_GAUSS: 112 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim, npoints, &npoints_q, &xi)); 113 break; 114 115 case DMSWARMPIC_LAYOUT_SUBDIVISION: { 116 PetscInt s, nsub; 117 PetscInt np_dir[3]; 118 nsub = npoints; 119 np_dir[0] = 1; 120 for (s = 0; s < nsub; s++) np_dir[0] *= 2; 121 np_dir[1] = np_dir[0]; 122 np_dir[2] = np_dir[0]; 123 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi)); 124 } break; 125 default: 126 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "A valid DMSwarmPIC layout must be provided"); 127 } 128 129 PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list)); 130 PetscCall(PetscMalloc1(dim * npe, &elcoor)); 131 PetscCall(PetscMalloc1(npoints_q, &basis)); 132 for (q = 0; q < npoints_q; q++) { 133 PetscCall(PetscMalloc1(npe, &basis[q])); 134 135 switch (dim) { 136 case 1: 137 basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]); 138 basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]); 139 break; 140 case 2: 141 basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]); 142 basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]); 143 basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]); 144 basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]); 145 break; 146 147 case 3: 148 basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 149 basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 150 basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 151 basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 152 basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 153 basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 154 basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 155 basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 156 break; 157 } 158 } 159 160 PetscCall(DMSwarmGetCellDMActive(dm, &celldm)); 161 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields)); 162 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc); 163 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid)); 164 165 PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1)); 166 PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 167 PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 168 169 PetscCall(DMGetCoordinatesLocal(dmc, &coor)); 170 PetscCall(VecGetArrayRead(coor, &_coor)); 171 pcnt = 0; 172 for (e = 0; e < nel; e++) { 173 const PetscInt *element = &element_list[npe * e]; 174 175 for (k = 0; k < npe; k++) { 176 for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]); 177 } 178 179 for (q = 0; q < npoints_q; q++) { 180 for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] = 0.0; 181 for (k = 0; k < npe; k++) { 182 for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] += basis[q][k] * elcoor[dim * k + d]; 183 } 184 swarm_cellid[pcnt] = e; 185 pcnt++; 186 } 187 } 188 PetscCall(VecRestoreArrayRead(coor, &_coor)); 189 PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 190 PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 191 PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list)); 192 193 PetscCall(PetscFree(xi)); 194 PetscCall(PetscFree(elcoor)); 195 for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q])); 196 PetscCall(PetscFree(basis)); 197 PetscFunctionReturn(PETSC_SUCCESS); 198 } 199 200 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param) 201 { 202 DMDAElementType etype; 203 PetscInt dim; 204 205 PetscFunctionBegin; 206 PetscCall(DMDAGetElementType(celldm, &etype)); 207 PetscCall(DMGetDimension(celldm, &dim)); 208 switch (etype) { 209 case DMDA_ELEMENT_P1: 210 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DA support is not currently available for DMDA_ELEMENT_P1"); 211 case DMDA_ELEMENT_Q1: 212 PetscCheck(dim != 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Support only available for dim = 2, 3"); 213 PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm, celldm, layout_param, layout)); 214 break; 215 } 216 PetscFunctionReturn(PETSC_SUCCESS); 217 } 218