1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/ 2 3 typedef struct _n_DMField_Shell { 4 void *ctx; 5 PetscErrorCode (*destroy)(DMField); 6 } DMField_Shell; 7 8 PetscErrorCode DMFieldShellGetContext(DMField field, void *ctx) { 9 PetscBool flg; 10 11 PetscFunctionBegin; 12 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1); 13 PetscValidPointer(ctx, 2); 14 PetscCall(PetscObjectTypeCompare((PetscObject)field, DMFIELDSHELL, &flg)); 15 if (flg) *(void **)ctx = ((DMField_Shell *)(field->data))->ctx; 16 else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Cannot get context from non-shell shield"); 17 PetscFunctionReturn(0); 18 } 19 20 static PetscErrorCode DMFieldDestroy_Shell(DMField field) { 21 DMField_Shell *shell = (DMField_Shell *)field->data; 22 23 PetscFunctionBegin; 24 if (shell->destroy) PetscCall((*(shell->destroy))(field)); 25 PetscCall(PetscFree(field->data)); 26 PetscFunctionReturn(0); 27 } 28 29 PetscErrorCode DMFieldShellEvaluateFEDefault(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H) { 30 DM dm = field->dm; 31 DMField coordField; 32 PetscFEGeom *geom; 33 Vec pushforward; 34 PetscInt dimC, dim, numPoints, Nq, p, Nc; 35 PetscScalar *pfArray; 36 37 PetscFunctionBegin; 38 Nc = field->numComponents; 39 PetscCall(DMGetCoordinateField(dm, &coordField)); 40 PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom)); 41 PetscCall(DMGetCoordinateDim(dm, &dimC)); 42 PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL)); 43 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 44 PetscCall(PetscMalloc1(dimC * Nq * numPoints, &pfArray)); 45 for (p = 0; p < numPoints * dimC * Nq; p++) pfArray[p] = (PetscScalar)geom->v[p]; 46 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * Nq * numPoints, PETSC_DETERMINE, pfArray, &pushforward)); 47 PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H)); 48 /* TODO: handle covariant/contravariant pullbacks */ 49 if (D) { 50 if (type == PETSC_SCALAR) { 51 PetscScalar *sD = (PetscScalar *)D; 52 PetscInt q; 53 54 for (p = 0; p < numPoints * Nq; p++) { 55 for (q = 0; q < Nc; q++) { 56 PetscScalar d[3]; 57 58 PetscInt i, j; 59 60 for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i]; 61 for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.; 62 for (i = 0; i < dimC; i++) { 63 for (j = 0; j < dimC; j++) { sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j]; } 64 } 65 } 66 } 67 } else { 68 PetscReal *rD = (PetscReal *)D; 69 PetscInt q; 70 71 for (p = 0; p < numPoints * Nq; p++) { 72 for (q = 0; q < Nc; q++) { 73 PetscReal d[3]; 74 75 PetscInt i, j; 76 77 for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i]; 78 for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.; 79 for (i = 0; i < dimC; i++) { 80 for (j = 0; j < dimC; j++) { rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j]; } 81 } 82 } 83 } 84 } 85 } 86 if (H) { 87 if (type == PETSC_SCALAR) { 88 PetscScalar *sH = (PetscScalar *)H; 89 PetscInt q; 90 91 for (p = 0; p < numPoints * Nq; p++) { 92 for (q = 0; q < Nc; q++) { 93 PetscScalar d[3][3]; 94 95 PetscInt i, j, k, l; 96 97 for (i = 0; i < dimC; i++) 98 for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j]; 99 for (i = 0; i < dimC; i++) 100 for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.; 101 for (i = 0; i < dimC; i++) { 102 for (j = 0; j < dimC; j++) { 103 for (k = 0; k < dimC; k++) { 104 for (l = 0; l < dimC; l++) { sH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l]; } 105 } 106 } 107 } 108 } 109 } 110 } else { 111 PetscReal *rH = (PetscReal *)H; 112 PetscInt q; 113 114 for (p = 0; p < numPoints * Nq; p++) { 115 for (q = 0; q < Nc; q++) { 116 PetscReal d[3][3]; 117 118 PetscInt i, j, k, l; 119 120 for (i = 0; i < dimC; i++) 121 for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j]; 122 for (i = 0; i < dimC; i++) 123 for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.; 124 for (i = 0; i < dimC; i++) { 125 for (j = 0; j < dimC; j++) { 126 for (k = 0; k < dimC; k++) { 127 for (l = 0; l < dimC; l++) { rH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l]; } 128 } 129 } 130 } 131 } 132 } 133 } 134 } 135 PetscCall(VecDestroy(&pushforward)); 136 PetscCall(PetscFree(pfArray)); 137 PetscCall(PetscFEGeomDestroy(&geom)); 138 PetscFunctionReturn(0); 139 } 140 141 PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H) { 142 DM dm = field->dm; 143 DMField coordField; 144 PetscFEGeom *geom; 145 Vec pushforward; 146 PetscInt dimC, dim, numPoints, Nq, p; 147 PetscScalar *pfArray; 148 PetscQuadrature quad; 149 MPI_Comm comm; 150 151 PetscFunctionBegin; 152 PetscCall(PetscObjectGetComm((PetscObject)field, &comm)); 153 PetscCall(DMGetDimension(dm, &dim)); 154 PetscCall(DMGetCoordinateDim(dm, &dimC)); 155 PetscCall(DMGetCoordinateField(dm, &coordField)); 156 PetscCall(DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad)); 157 PetscCheck(quad, comm, PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation"); 158 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL)); 159 PetscCheck(Nq == 1, comm, PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point"); 160 PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom)); 161 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 162 PetscCall(PetscMalloc1(dimC * numPoints, &pfArray)); 163 for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar)geom->v[p]; 164 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward)); 165 PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H)); 166 PetscCall(PetscQuadratureDestroy(&quad)); 167 PetscCall(VecDestroy(&pushforward)); 168 PetscCall(PetscFree(pfArray)); 169 PetscCall(PetscFEGeomDestroy(&geom)); 170 PetscFunctionReturn(0); 171 } 172 173 PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField)) { 174 DMField_Shell *shell = (DMField_Shell *)field->data; 175 176 PetscFunctionBegin; 177 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1); 178 shell->destroy = destroy; 179 PetscFunctionReturn(0); 180 } 181 182 PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField, Vec, PetscDataType, void *, void *, void *)) { 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1); 185 field->ops->evaluate = evaluate; 186 PetscFunctionReturn(0); 187 } 188 189 PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField, IS, PetscQuadrature, PetscDataType, void *, void *, void *)) { 190 PetscFunctionBegin; 191 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1); 192 field->ops->evaluateFE = evaluateFE; 193 PetscFunctionReturn(0); 194 } 195 196 PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField, IS, PetscDataType, void *, void *, void *)) { 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1); 199 field->ops->evaluateFV = evaluateFV; 200 PetscFunctionReturn(0); 201 } 202 203 PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField, IS, PetscInt *, PetscInt *)) { 204 PetscFunctionBegin; 205 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1); 206 field->ops->getDegree = getDegree; 207 PetscFunctionReturn(0); 208 } 209 210 PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField, IS, PetscQuadrature *)) { 211 PetscFunctionBegin; 212 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1); 213 field->ops->createDefaultQuadrature = createDefaultQuadrature; 214 PetscFunctionReturn(0); 215 } 216 217 static PetscErrorCode DMFieldInitialize_Shell(DMField field) { 218 PetscFunctionBegin; 219 field->ops->destroy = DMFieldDestroy_Shell; 220 field->ops->evaluate = NULL; 221 field->ops->evaluateFE = DMFieldShellEvaluateFEDefault; 222 field->ops->evaluateFV = DMFieldShellEvaluateFVDefault; 223 field->ops->getDegree = NULL; 224 field->ops->createDefaultQuadrature = NULL; 225 field->ops->view = NULL; 226 PetscFunctionReturn(0); 227 } 228 229 PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field) { 230 DMField_Shell *shell; 231 232 PetscFunctionBegin; 233 PetscCall(PetscNewLog(field, &shell)); 234 field->data = shell; 235 PetscCall(DMFieldInitialize_Shell(field)); 236 PetscFunctionReturn(0); 237 } 238 239 PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field) { 240 DMField b; 241 DMField_Shell *shell; 242 243 PetscFunctionBegin; 244 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 245 if (ctx) PetscValidPointer(ctx, 4); 246 PetscValidPointer(field, 5); 247 PetscCall(DMFieldCreate(dm, numComponents, continuity, &b)); 248 PetscCall(DMFieldSetType(b, DMFIELDSHELL)); 249 shell = (DMField_Shell *)b->data; 250 shell->ctx = ctx; 251 *field = b; 252 PetscFunctionReturn(0); 253 } 254