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