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 PetscErrorCode ierr; 163 164 PetscFunctionBegin; 165 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 166 ierr = DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);CHKERRQ(ierr); 167 if (!quad) SETERRQ(PetscObjectComm((PetscObject) pointIS), PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation"); 168 ierr = DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);CHKERRQ(ierr); 169 ierr = DMGetCoordinateDim(dm, &dimC);CHKERRQ(ierr); 170 ierr = PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 171 if (Nq != 1) SETERRQ(PetscObjectComm((PetscObject) quad), PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point"); 172 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 173 ierr = PetscMalloc1(dimC * numPoints, &pfArray);CHKERRQ(ierr); 174 for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar) geom->v[p]; 175 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward);CHKERRQ(ierr); 176 ierr = DMFieldEvaluate(field, pushforward, type, B, D, H);CHKERRQ(ierr); 177 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 178 ierr = VecDestroy(&pushforward);CHKERRQ(ierr); 179 ierr = PetscFree(pfArray);CHKERRQ(ierr); 180 ierr = PetscFEGeomDestroy(&geom);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField)) 185 { 186 DMField_Shell *shell = (DMField_Shell *) field->data; 187 188 PetscFunctionBegin; 189 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 190 shell->destroy = destroy; 191 PetscFunctionReturn(0); 192 } 193 194 PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField,Vec,PetscDataType,void*,void*,void*)) 195 { 196 PetscFunctionBegin; 197 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 198 field->ops->evaluate = evaluate; 199 PetscFunctionReturn(0); 200 } 201 202 PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField,IS,PetscQuadrature,PetscDataType,void*,void*,void*)) 203 { 204 PetscFunctionBegin; 205 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 206 field->ops->evaluateFE = evaluateFE; 207 PetscFunctionReturn(0); 208 } 209 210 PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField,IS,PetscDataType,void*,void*,void*)) 211 { 212 PetscFunctionBegin; 213 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 214 field->ops->evaluateFV = evaluateFV; 215 PetscFunctionReturn(0); 216 } 217 218 PetscErrorCode DMFieldShellSetGetFEInvariance(DMField field, PetscErrorCode (*getFEInvariance)(DMField,IS,PetscBool*,PetscBool*,PetscBool*)) 219 { 220 PetscFunctionBegin; 221 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 222 field->ops->getFEInvariance = getFEInvariance; 223 PetscFunctionReturn(0); 224 } 225 226 PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField,IS,PetscQuadrature*)) 227 { 228 PetscFunctionBegin; 229 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 230 field->ops->createDefaultQuadrature = createDefaultQuadrature; 231 PetscFunctionReturn(0); 232 } 233 234 static PetscErrorCode DMFieldInitialize_Shell(DMField field) 235 { 236 PetscFunctionBegin; 237 field->ops->destroy = DMFieldDestroy_Shell; 238 field->ops->evaluate = NULL; 239 field->ops->evaluateFE = DMFieldShellEvaluateFEDefault; 240 field->ops->evaluateFV = DMFieldShellEvaluateFVDefault; 241 field->ops->getFEInvariance = NULL; 242 field->ops->createDefaultQuadrature = NULL; 243 field->ops->view = NULL; 244 PetscFunctionReturn(0); 245 } 246 247 PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field) 248 { 249 DMField_Shell *shell; 250 PetscErrorCode ierr; 251 252 PetscFunctionBegin; 253 ierr = PetscNewLog(field,&shell);CHKERRQ(ierr); 254 field->data = shell; 255 ierr = DMFieldInitialize_Shell(field);CHKERRQ(ierr); 256 PetscFunctionReturn(0); 257 } 258 259 PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field) 260 { 261 DMField b; 262 DMField_Shell *shell; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 267 if (ctx) PetscValidPointer(ctx, 4); 268 PetscValidPointer(field, 5); 269 ierr = DMFieldCreate(dm,numComponents,continuity,&b);CHKERRQ(ierr); 270 ierr = DMFieldSetType(b,DMFIELDSHELL);CHKERRQ(ierr); 271 shell = (DMField_Shell *) b->data; 272 shell->ctx = ctx; 273 *field = b; 274 PetscFunctionReturn(0); 275 } 276