1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/ 2 #include <petsc/private/petscfeimpl.h> /*I "petscdmfield.h" I*/ 3 #include <petscdmplex.h> 4 5 const char *const DMFieldContinuities[] = { 6 "VERTEX", 7 "EDGE", 8 "FACET", 9 "CELL", 10 0 11 }; 12 13 PETSC_INTERN PetscErrorCode DMFieldCreate(DM dm,PetscInt numComponents,DMFieldContinuity continuity,DMField *field) 14 { 15 PetscErrorCode ierr; 16 DMField b; 17 18 PetscFunctionBegin; 19 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 20 PetscValidPointer(field,2); 21 ierr = DMFieldInitializePackage();CHKERRQ(ierr); 22 23 ierr = PetscHeaderCreate(b,DMFIELD_CLASSID,"DMField","Field over DM","DM",PetscObjectComm((PetscObject)dm),DMFieldDestroy,DMFieldView);CHKERRQ(ierr); 24 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 25 b->dm = dm; 26 b->continuity = continuity; 27 b->numComponents = numComponents; 28 *field = b; 29 PetscFunctionReturn(0); 30 } 31 32 /*@ 33 DMFieldDestroy - destroy a DMField 34 35 Collective 36 37 Input Arguments: 38 . field - address of DMField 39 40 Level: advanced 41 42 .seealso: DMFieldCreate() 43 @*/ 44 PetscErrorCode DMFieldDestroy(DMField *field) 45 { 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 if (!*field) PetscFunctionReturn(0); 50 PetscValidHeaderSpecific((*field),DMFIELD_CLASSID,1); 51 if (--((PetscObject)(*field))->refct > 0) {*field = 0; PetscFunctionReturn(0);} 52 if ((*field)->ops->destroy) {ierr = (*(*field)->ops->destroy)(*field);CHKERRQ(ierr);} 53 ierr = DMDestroy(&((*field)->dm));CHKERRQ(ierr); 54 ierr = PetscHeaderDestroy(field);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 /*@C 59 DMFieldView - view a DMField 60 61 Collective 62 63 Input Arguments: 64 + field - DMField 65 - viewer - viewer to display field, for example PETSC_VIEWER_STDOUT_WORLD 66 67 Level: advanced 68 69 .seealso: DMFieldCreate() 70 @*/ 71 PetscErrorCode DMFieldView(DMField field,PetscViewer viewer) 72 { 73 PetscErrorCode ierr; 74 PetscBool iascii; 75 76 PetscFunctionBegin; 77 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 78 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field),&viewer);CHKERRQ(ierr);} 79 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 80 PetscCheckSameComm(field,1,viewer,2); 81 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 82 if (iascii) { 83 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)field,viewer);CHKERRQ(ierr); 84 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 85 ierr = PetscViewerASCIIPrintf(viewer,"%D components\n",field->numComponents);CHKERRQ(ierr); 86 ierr = PetscViewerASCIIPrintf(viewer,"%s continuity\n",DMFieldContinuities[field->continuity]);CHKERRQ(ierr); 87 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 88 ierr = DMView(field->dm,viewer);CHKERRQ(ierr); 89 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 90 } 91 if (field->ops->view) {ierr = (*field->ops->view)(field,viewer);CHKERRQ(ierr);} 92 if (iascii) { 93 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 94 } 95 PetscFunctionReturn(0); 96 } 97 98 /*@C 99 DMFieldSetType - set the DMField implementation 100 101 Collective on DMField 102 103 Input Parameters: 104 + field - the DMField context 105 - type - a known method 106 107 Notes: 108 See "include/petscvec.h" for available methods (for instance) 109 + DMFIELDDA - a field defined only by its values at the corners of a DMDA 110 . DMFIELDDS - a field defined by a discretization over a mesh set with DMSetField() 111 - DMFIELDSHELL - a field defined by arbitrary callbacks 112 113 Level: advanced 114 115 .keywords: DMField, set, type 116 117 .seealso: DMFieldType, 118 @*/ 119 PetscErrorCode DMFieldSetType(DMField field,DMFieldType type) 120 { 121 PetscErrorCode ierr,(*r)(DMField); 122 PetscBool match; 123 124 PetscFunctionBegin; 125 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 126 PetscValidCharPointer(type,2); 127 128 ierr = PetscObjectTypeCompare((PetscObject)field,type,&match);CHKERRQ(ierr); 129 if (match) PetscFunctionReturn(0); 130 131 ierr = PetscFunctionListFind(DMFieldList,type,&r);CHKERRQ(ierr); 132 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DMField type %s",type); 133 /* Destroy the previous private DMField context */ 134 if (field->ops->destroy) { 135 ierr = (*(field)->ops->destroy)(field);CHKERRQ(ierr); 136 } 137 ierr = PetscMemzero(field->ops,sizeof(*field->ops));CHKERRQ(ierr); 138 ierr = PetscObjectChangeTypeName((PetscObject)field,type);CHKERRQ(ierr); 139 field->ops->create = r; 140 ierr = (*r)(field);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 /*@C 145 DMFieldGetType - Gets the DMField type name (as a string) from the DMField. 146 147 Not Collective 148 149 Input Parameter: 150 . field - The DMField context 151 152 Output Parameter: 153 . type - The DMField type name 154 155 Level: advanced 156 157 .keywords: DMField, get, type, name 158 .seealso: DMFieldSetType() 159 @*/ 160 PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type) 161 { 162 PetscErrorCode ierr; 163 164 PetscFunctionBegin; 165 PetscValidHeaderSpecific(field, VEC_TAGGER_CLASSID,1); 166 PetscValidPointer(type,2); 167 ierr = DMFieldRegisterAll();CHKERRQ(ierr); 168 *type = ((PetscObject)field)->type_name; 169 PetscFunctionReturn(0); 170 } 171 172 PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc) 173 { 174 PetscFunctionBegin; 175 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 176 PetscValidIntPointer(nc,2); 177 *nc = field->numComponents; 178 PetscFunctionReturn(0); 179 } 180 181 PetscErrorCode DMFieldGetDM(DMField field, DM *dm) 182 { 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 185 PetscValidPointer(dm,2); 186 *dm = field->dm; 187 PetscFunctionReturn(0); 188 } 189 190 PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H) 191 { 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 196 PetscValidHeaderSpecific(points,VEC_CLASSID,2); 197 if (B) PetscValidPointer(B,3); 198 if (D) PetscValidPointer(D,4); 199 if (H) PetscValidPointer(H,5); 200 if (field->ops->evaluate) { 201 ierr = (*field->ops->evaluate) (field, points, datatype, B, D, H);CHKERRQ(ierr); 202 } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type"); 203 PetscFunctionReturn(0); 204 } 205 206 PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H) 207 { 208 PetscErrorCode ierr; 209 210 PetscFunctionBegin; 211 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 212 PetscValidHeaderSpecific(cellIS,IS_CLASSID,2); 213 PetscValidHeader(points,3); 214 if (B) PetscValidPointer(B,4); 215 if (D) PetscValidPointer(D,5); 216 if (H) PetscValidPointer(H,6); 217 if (field->ops->evaluateFE) { 218 ierr = (*field->ops->evaluateFE) (field, cellIS, points, datatype, B, D, H);CHKERRQ(ierr); 219 } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type"); 220 PetscFunctionReturn(0); 221 } 222 223 PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H) 224 { 225 PetscErrorCode ierr; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 229 PetscValidHeaderSpecific(cellIS,IS_CLASSID,2); 230 if (B) PetscValidPointer(B,3); 231 if (D) PetscValidPointer(D,4); 232 if (H) PetscValidPointer(H,5); 233 if (field->ops->evaluateFV) { 234 ierr = (*field->ops->evaluateFV) (field, cellIS, datatype, B, D, H);CHKERRQ(ierr); 235 } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type"); 236 PetscFunctionReturn(0); 237 } 238 239 PetscErrorCode DMFieldGetFEInvariance(DMField field, IS cellIS, PetscBool *isConstant, PetscBool *isAffine, PetscBool *isQuadratic) 240 { 241 PetscErrorCode ierr; 242 243 PetscFunctionBegin; 244 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 245 PetscValidHeaderSpecific(cellIS,IS_CLASSID,2); 246 if (isConstant) PetscValidPointer(isConstant,3); 247 if (isAffine) PetscValidPointer(isAffine,4); 248 if (isQuadratic) PetscValidPointer(isQuadratic,5); 249 250 if (isConstant) *isConstant = PETSC_FALSE; 251 if (isAffine) *isAffine = PETSC_FALSE; 252 if (isQuadratic) *isQuadratic = PETSC_FALSE; 253 254 if (field->ops->getFEInvariance) { 255 ierr = (*field->ops->getFEInvariance) (field,cellIS,isConstant,isAffine,isQuadratic);CHKERRQ(ierr); 256 } 257 PetscFunctionReturn(0); 258 } 259 260 PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad) 261 { 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 266 PetscValidHeaderSpecific(pointIS,IS_CLASSID,2); 267 PetscValidPointer(quad,3); 268 269 *quad = NULL; 270 if (field->ops->createDefaultQuadrature) { 271 ierr = (*field->ops->createDefaultQuadrature)(field, pointIS, quad);CHKERRQ(ierr); 272 } 273 PetscFunctionReturn(0); 274 } 275 276 PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 277 { 278 PetscInt dim, dE; 279 PetscInt nPoints; 280 PetscFEGeom *g; 281 PetscErrorCode ierr; 282 283 PetscFunctionBegin; 284 PetscValidHeaderSpecific(field,DMFIELD_CLASSID,1); 285 PetscValidHeaderSpecific(pointIS,IS_CLASSID,2); 286 PetscValidHeader(quad,3); 287 ierr = ISGetLocalSize(pointIS,&nPoints);CHKERRQ(ierr); 288 dE = field->numComponents; 289 ierr = PetscFEGeomCreate(quad,nPoints,dE,faceData,&g);CHKERRQ(ierr); 290 ierr = DMFieldEvaluateFE(field,pointIS,quad,PETSC_REAL,g->v,g->J,NULL);CHKERRQ(ierr); 291 dim = g->dim; 292 if (dE > dim) { 293 /* space out J and make square Jacobians */ 294 PetscInt i, j, k, N = g->numPoints * g->numCells; 295 296 for (i = N-1; i >= 0; i--) { 297 PetscReal J[9]; 298 299 for (j = 0; j < dE; j++) { 300 for (k = 0; k < dim; k++) { 301 J[j*dE + k] = g->J[i*dE*dim + j*dim + k]; 302 } 303 } 304 switch (dim) { 305 case 0: 306 for (j = 0; j < dE; j++) { 307 for (k = 0; k < dE; k++) { 308 J[j * dE + k] = (j == k) ? 1. : 0.; 309 } 310 } 311 break; 312 case 1: 313 if (dE == 2) { 314 PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]); 315 316 J[1] = -J[2] / norm; 317 J[3] = J[0] / norm; 318 } else { 319 PetscReal inorm = 1./PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]); 320 PetscReal x = J[0] * inorm; 321 PetscReal y = J[3] * inorm; 322 PetscReal z = J[6] * inorm; 323 324 if (x > 0.) { 325 PetscReal inv1pX = 1./ (1. + x); 326 327 J[1] = -y; J[2] = -z; 328 J[4] = 1. - y*y*inv1pX; J[5] = -y*z*inv1pX; 329 J[7] = -y*z*inv1pX; J[8] = 1. - z*z*inv1pX; 330 } else { 331 PetscReal inv1mX = 1./ (1. - x); 332 333 J[1] = z; J[2] = y; 334 J[4] = -y*z*inv1mX; J[5] = 1. - y*y*inv1mX; 335 J[7] = 1. - z*z*inv1mX; J[8] = -y*z*inv1mX; 336 } 337 } 338 break; 339 case 2: 340 { 341 PetscReal inorm; 342 343 J[2] = J[3] * J[7] - J[6] * J[4]; 344 J[5] = J[6] * J[1] - J[0] * J[7]; 345 J[8] = J[0] * J[4] - J[3] * J[1]; 346 347 inorm = 1./ PetscSqrtReal(J[2]*J[2] + J[5]*J[5] + J[8]*J[8]); 348 349 J[2] *= inorm; 350 J[5] *= inorm; 351 J[8] *= inorm; 352 } 353 break; 354 } 355 for (j = 0; j < dE*dE; j++) { 356 g->J[i*dE*dE + j] = J[j]; 357 } 358 } 359 } 360 ierr = PetscFEGeomComplete(g);CHKERRQ(ierr); 361 ierr = DMFieldGetFEInvariance(field,pointIS,NULL,&g->isAffine,NULL);CHKERRQ(ierr); 362 if (faceData) { 363 if (!field->ops->computeFaceData) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMField implementation does not compute face data\n"); 364 ierr = (*field->ops->computeFaceData) (field, pointIS, quad, g);CHKERRQ(ierr); 365 } 366 *geom = g; 367 PetscFunctionReturn(0); 368 } 369