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