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