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