1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 4 PetscClassId PETSCDUALSPACE_CLASSID = 0; 5 6 PetscLogEvent PETSCDUALSPACE_SetUp; 7 8 PetscFunctionList PetscDualSpaceList = NULL; 9 PetscBool PetscDualSpaceRegisterAllCalled = PETSC_FALSE; 10 11 /* 12 PetscDualSpaceLatticePointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to at most 'max'. 13 Ordering is lexicographic with lowest index as least significant in ordering. 14 e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {0,2}. 15 16 Input Parameters: 17 + len - The length of the tuple 18 . max - The maximum sum 19 - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition 20 21 Output Parameter: 22 . tup - A tuple of len integers whos sum is at most 'max' 23 24 Level: developer 25 26 .seealso: PetscDualSpaceTensorPointLexicographic_Internal() 27 */ 28 PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[]) 29 { 30 PetscFunctionBegin; 31 while (len--) { 32 max -= tup[len]; 33 if (!max) { 34 tup[len] = 0; 35 break; 36 } 37 } 38 tup[++len]++; 39 PetscFunctionReturn(0); 40 } 41 42 /* 43 PetscDualSpaceTensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'. 44 Ordering is lexicographic with lowest index as least significant in ordering. 45 e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,1}, {0,2}, {1,2}, {2,2}. 46 47 Input Parameters: 48 + len - The length of the tuple 49 . max - The maximum value 50 - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition 51 52 Output Parameter: 53 . tup - A tuple of len integers whos sum is at most 'max' 54 55 Level: developer 56 57 .seealso: PetscDualSpaceLatticePointLexicographic_Internal() 58 */ 59 PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[]) 60 { 61 PetscInt i; 62 63 PetscFunctionBegin; 64 for (i = 0; i < len; i++) { 65 if (tup[i] < max) { 66 break; 67 } else { 68 tup[i] = 0; 69 } 70 } 71 tup[i]++; 72 PetscFunctionReturn(0); 73 } 74 75 /*@C 76 PetscDualSpaceRegister - Adds a new PetscDualSpace implementation 77 78 Not Collective 79 80 Input Parameters: 81 + name - The name of a new user-defined creation routine 82 - create_func - The creation routine itself 83 84 Notes: 85 PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces 86 87 Sample usage: 88 .vb 89 PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate); 90 .ve 91 92 Then, your PetscDualSpace type can be chosen with the procedural interface via 93 .vb 94 PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *); 95 PetscDualSpaceSetType(PetscDualSpace, "my_dual_space"); 96 .ve 97 or at runtime via the option 98 .vb 99 -petscdualspace_type my_dual_space 100 .ve 101 102 Level: advanced 103 104 .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy() 105 106 @*/ 107 PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace)) 108 { 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 ierr = PetscFunctionListAdd(&PetscDualSpaceList, sname, function);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 /*@C 117 PetscDualSpaceSetType - Builds a particular PetscDualSpace 118 119 Collective on sp 120 121 Input Parameters: 122 + sp - The PetscDualSpace object 123 - name - The kind of space 124 125 Options Database Key: 126 . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types 127 128 Level: intermediate 129 130 .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate() 131 @*/ 132 PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name) 133 { 134 PetscErrorCode (*r)(PetscDualSpace); 135 PetscBool match; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 140 ierr = PetscObjectTypeCompare((PetscObject) sp, name, &match);CHKERRQ(ierr); 141 if (match) PetscFunctionReturn(0); 142 143 if (!PetscDualSpaceRegisterAllCalled) {ierr = PetscDualSpaceRegisterAll();CHKERRQ(ierr);} 144 ierr = PetscFunctionListFind(PetscDualSpaceList, name, &r);CHKERRQ(ierr); 145 PetscCheckFalse(!r,PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name); 146 147 if (sp->ops->destroy) { 148 ierr = (*sp->ops->destroy)(sp);CHKERRQ(ierr); 149 sp->ops->destroy = NULL; 150 } 151 ierr = (*r)(sp);CHKERRQ(ierr); 152 ierr = PetscObjectChangeTypeName((PetscObject) sp, name);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 /*@C 157 PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object. 158 159 Not Collective 160 161 Input Parameter: 162 . sp - The PetscDualSpace 163 164 Output Parameter: 165 . name - The PetscDualSpace type name 166 167 Level: intermediate 168 169 .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate() 170 @*/ 171 PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name) 172 { 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 177 PetscValidPointer(name, 2); 178 if (!PetscDualSpaceRegisterAllCalled) { 179 ierr = PetscDualSpaceRegisterAll();CHKERRQ(ierr); 180 } 181 *name = ((PetscObject) sp)->type_name; 182 PetscFunctionReturn(0); 183 } 184 185 static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v) 186 { 187 PetscViewerFormat format; 188 PetscInt pdim, f; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 ierr = PetscDualSpaceGetDimension(sp, &pdim);CHKERRQ(ierr); 193 ierr = PetscObjectPrintClassNamePrefixType((PetscObject) sp, v);CHKERRQ(ierr); 194 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 195 if (sp->k) { 196 ierr = PetscViewerASCIIPrintf(v, "Dual space for %D-forms %swith %D components, size %D\n", PetscAbsInt(sp->k), sp->k < 0 ? "(stored in dual form) ": "", sp->Nc, pdim);CHKERRQ(ierr); 197 } else { 198 ierr = PetscViewerASCIIPrintf(v, "Dual space with %D components, size %D\n", sp->Nc, pdim);CHKERRQ(ierr); 199 } 200 if (sp->ops->view) {ierr = (*sp->ops->view)(sp, v);CHKERRQ(ierr);} 201 ierr = PetscViewerGetFormat(v, &format);CHKERRQ(ierr); 202 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 203 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 204 for (f = 0; f < pdim; ++f) { 205 ierr = PetscViewerASCIIPrintf(v, "Dual basis vector %D\n", f);CHKERRQ(ierr); 206 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 207 ierr = PetscQuadratureView(sp->functional[f], v);CHKERRQ(ierr); 208 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 209 } 210 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 211 } 212 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 /*@C 217 PetscDualSpaceViewFromOptions - View from Options 218 219 Collective on PetscDualSpace 220 221 Input Parameters: 222 + A - the PetscDualSpace object 223 . obj - Optional object, proivides prefix 224 - name - command line option 225 226 Level: intermediate 227 .seealso: PetscDualSpace, PetscDualSpaceView(), PetscObjectViewFromOptions(), PetscDualSpaceCreate() 228 @*/ 229 PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace A,PetscObject obj,const char name[]) 230 { 231 PetscErrorCode ierr; 232 233 PetscFunctionBegin; 234 PetscValidHeaderSpecific(A,PETSCDUALSPACE_CLASSID,1); 235 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 239 /*@ 240 PetscDualSpaceView - Views a PetscDualSpace 241 242 Collective on sp 243 244 Input Parameters: 245 + sp - the PetscDualSpace object to view 246 - v - the viewer 247 248 Level: beginner 249 250 .seealso PetscDualSpaceDestroy(), PetscDualSpace 251 @*/ 252 PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v) 253 { 254 PetscBool iascii; 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 259 if (v) PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2); 260 if (!v) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);CHKERRQ(ierr);} 261 ierr = PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 262 if (iascii) {ierr = PetscDualSpaceView_ASCII(sp, v);CHKERRQ(ierr);} 263 PetscFunctionReturn(0); 264 } 265 266 /*@ 267 PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database 268 269 Collective on sp 270 271 Input Parameter: 272 . sp - the PetscDualSpace object to set options for 273 274 Options Database: 275 + -petscdualspace_order <order> - the approximation order of the space 276 . -petscdualspace_form_degree <deg> - the form degree, say 0 for point evaluations, or 2 for area integrals 277 . -petscdualspace_components <c> - the number of components, say d for a vector field 278 - -petscdualspace_refcell <celltype> - Reference cell type name 279 280 Level: intermediate 281 282 .seealso PetscDualSpaceView(), PetscDualSpace, PetscObjectSetFromOptions() 283 @*/ 284 PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp) 285 { 286 DMPolytopeType refCell = DM_POLYTOPE_TRIANGLE; 287 const char *defaultType; 288 char name[256]; 289 PetscBool flg; 290 PetscErrorCode ierr; 291 292 PetscFunctionBegin; 293 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 294 if (!((PetscObject) sp)->type_name) { 295 defaultType = PETSCDUALSPACELAGRANGE; 296 } else { 297 defaultType = ((PetscObject) sp)->type_name; 298 } 299 if (!PetscSpaceRegisterAllCalled) {ierr = PetscSpaceRegisterAll();CHKERRQ(ierr);} 300 301 ierr = PetscObjectOptionsBegin((PetscObject) sp);CHKERRQ(ierr); 302 ierr = PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);CHKERRQ(ierr); 303 if (flg) { 304 ierr = PetscDualSpaceSetType(sp, name);CHKERRQ(ierr); 305 } else if (!((PetscObject) sp)->type_name) { 306 ierr = PetscDualSpaceSetType(sp, defaultType);CHKERRQ(ierr); 307 } 308 ierr = PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL,0);CHKERRQ(ierr); 309 ierr = PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL);CHKERRQ(ierr); 310 ierr = PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL,1);CHKERRQ(ierr); 311 if (sp->ops->setfromoptions) { 312 ierr = (*sp->ops->setfromoptions)(PetscOptionsObject,sp);CHKERRQ(ierr); 313 } 314 ierr = PetscOptionsEnum("-petscdualspace_refcell", "Reference cell shape", "PetscDualSpaceSetReferenceCell", DMPolytopeTypes, (PetscEnum) refCell, (PetscEnum *) &refCell, &flg);CHKERRQ(ierr); 315 if (flg) { 316 DM K; 317 318 ierr = DMPlexCreateReferenceCell(PETSC_COMM_SELF, refCell, &K);CHKERRQ(ierr); 319 ierr = PetscDualSpaceSetDM(sp, K);CHKERRQ(ierr); 320 ierr = DMDestroy(&K);CHKERRQ(ierr); 321 } 322 323 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 324 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);CHKERRQ(ierr); 325 ierr = PetscOptionsEnd();CHKERRQ(ierr); 326 sp->setfromoptionscalled = PETSC_TRUE; 327 PetscFunctionReturn(0); 328 } 329 330 /*@ 331 PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace 332 333 Collective on sp 334 335 Input Parameter: 336 . sp - the PetscDualSpace object to setup 337 338 Level: intermediate 339 340 .seealso PetscDualSpaceView(), PetscDualSpaceDestroy(), PetscDualSpace 341 @*/ 342 PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp) 343 { 344 PetscErrorCode ierr; 345 346 PetscFunctionBegin; 347 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 348 if (sp->setupcalled) PetscFunctionReturn(0); 349 ierr = PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);CHKERRQ(ierr); 350 sp->setupcalled = PETSC_TRUE; 351 if (sp->ops->setup) {ierr = (*sp->ops->setup)(sp);CHKERRQ(ierr);} 352 ierr = PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);CHKERRQ(ierr); 353 if (sp->setfromoptionscalled) {ierr = PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");CHKERRQ(ierr);} 354 PetscFunctionReturn(0); 355 } 356 357 static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm) 358 { 359 PetscInt pStart = -1, pEnd = -1, depth = -1; 360 PetscErrorCode ierr; 361 362 PetscFunctionBegin; 363 if (!dm) PetscFunctionReturn(0); 364 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 365 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 366 367 if (sp->pointSpaces) { 368 PetscInt i; 369 370 for (i = 0; i < pEnd - pStart; i++) { 371 ierr = PetscDualSpaceDestroy(&(sp->pointSpaces[i]));CHKERRQ(ierr); 372 } 373 } 374 ierr = PetscFree(sp->pointSpaces);CHKERRQ(ierr); 375 376 if (sp->heightSpaces) { 377 PetscInt i; 378 379 for (i = 0; i <= depth; i++) { 380 ierr = PetscDualSpaceDestroy(&(sp->heightSpaces[i]));CHKERRQ(ierr); 381 } 382 } 383 ierr = PetscFree(sp->heightSpaces);CHKERRQ(ierr); 384 385 ierr = PetscSectionDestroy(&(sp->pointSection));CHKERRQ(ierr); 386 ierr = PetscQuadratureDestroy(&(sp->intNodes));CHKERRQ(ierr); 387 ierr = VecDestroy(&(sp->intDofValues));CHKERRQ(ierr); 388 ierr = VecDestroy(&(sp->intNodeValues));CHKERRQ(ierr); 389 ierr = MatDestroy(&(sp->intMat));CHKERRQ(ierr); 390 ierr = PetscQuadratureDestroy(&(sp->allNodes));CHKERRQ(ierr); 391 ierr = VecDestroy(&(sp->allDofValues));CHKERRQ(ierr); 392 ierr = VecDestroy(&(sp->allNodeValues));CHKERRQ(ierr); 393 ierr = MatDestroy(&(sp->allMat));CHKERRQ(ierr); 394 ierr = PetscFree(sp->numDof);CHKERRQ(ierr); 395 PetscFunctionReturn(0); 396 } 397 398 /*@ 399 PetscDualSpaceDestroy - Destroys a PetscDualSpace object 400 401 Collective on sp 402 403 Input Parameter: 404 . sp - the PetscDualSpace object to destroy 405 406 Level: beginner 407 408 .seealso PetscDualSpaceView(), PetscDualSpace(), PetscDualSpaceCreate() 409 @*/ 410 PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp) 411 { 412 PetscInt dim, f; 413 DM dm; 414 PetscErrorCode ierr; 415 416 PetscFunctionBegin; 417 if (!*sp) PetscFunctionReturn(0); 418 PetscValidHeaderSpecific((*sp), PETSCDUALSPACE_CLASSID, 1); 419 420 if (--((PetscObject)(*sp))->refct > 0) {*sp = NULL; PetscFunctionReturn(0);} 421 ((PetscObject) (*sp))->refct = 0; 422 423 ierr = PetscDualSpaceGetDimension(*sp, &dim);CHKERRQ(ierr); 424 dm = (*sp)->dm; 425 426 if ((*sp)->ops->destroy) {ierr = (*(*sp)->ops->destroy)(*sp);CHKERRQ(ierr);} 427 ierr = PetscDualSpaceClearDMData_Internal(*sp, dm);CHKERRQ(ierr); 428 429 for (f = 0; f < dim; ++f) { 430 ierr = PetscQuadratureDestroy(&(*sp)->functional[f]);CHKERRQ(ierr); 431 } 432 ierr = PetscFree((*sp)->functional);CHKERRQ(ierr); 433 ierr = DMDestroy(&(*sp)->dm);CHKERRQ(ierr); 434 ierr = PetscHeaderDestroy(sp);CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437 438 /*@ 439 PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType(). 440 441 Collective 442 443 Input Parameter: 444 . comm - The communicator for the PetscDualSpace object 445 446 Output Parameter: 447 . sp - The PetscDualSpace object 448 449 Level: beginner 450 451 .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE 452 @*/ 453 PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp) 454 { 455 PetscDualSpace s; 456 PetscErrorCode ierr; 457 458 PetscFunctionBegin; 459 PetscValidPointer(sp, 2); 460 ierr = PetscCitationsRegister(FECitation,&FEcite);CHKERRQ(ierr); 461 *sp = NULL; 462 ierr = PetscFEInitializePackage();CHKERRQ(ierr); 463 464 ierr = PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);CHKERRQ(ierr); 465 466 s->order = 0; 467 s->Nc = 1; 468 s->k = 0; 469 s->spdim = -1; 470 s->spintdim = -1; 471 s->uniform = PETSC_TRUE; 472 s->setupcalled = PETSC_FALSE; 473 474 *sp = s; 475 PetscFunctionReturn(0); 476 } 477 478 /*@ 479 PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup. 480 481 Collective on sp 482 483 Input Parameter: 484 . sp - The original PetscDualSpace 485 486 Output Parameter: 487 . spNew - The duplicate PetscDualSpace 488 489 Level: beginner 490 491 .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType() 492 @*/ 493 PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew) 494 { 495 DM dm; 496 PetscDualSpaceType type; 497 const char *name; 498 PetscErrorCode ierr; 499 500 PetscFunctionBegin; 501 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 502 PetscValidPointer(spNew, 2); 503 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew);CHKERRQ(ierr); 504 ierr = PetscObjectGetName((PetscObject) sp, &name);CHKERRQ(ierr); 505 ierr = PetscObjectSetName((PetscObject) *spNew, name);CHKERRQ(ierr); 506 ierr = PetscDualSpaceGetType(sp, &type);CHKERRQ(ierr); 507 ierr = PetscDualSpaceSetType(*spNew, type);CHKERRQ(ierr); 508 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 509 ierr = PetscDualSpaceSetDM(*spNew, dm);CHKERRQ(ierr); 510 511 (*spNew)->order = sp->order; 512 (*spNew)->k = sp->k; 513 (*spNew)->Nc = sp->Nc; 514 (*spNew)->uniform = sp->uniform; 515 if (sp->ops->duplicate) {ierr = (*sp->ops->duplicate)(sp, *spNew);CHKERRQ(ierr);} 516 PetscFunctionReturn(0); 517 } 518 519 /*@ 520 PetscDualSpaceGetDM - Get the DM representing the reference cell 521 522 Not collective 523 524 Input Parameter: 525 . sp - The PetscDualSpace 526 527 Output Parameter: 528 . dm - The reference cell 529 530 Level: intermediate 531 532 .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate() 533 @*/ 534 PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm) 535 { 536 PetscFunctionBegin; 537 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 538 PetscValidPointer(dm, 2); 539 *dm = sp->dm; 540 PetscFunctionReturn(0); 541 } 542 543 /*@ 544 PetscDualSpaceSetDM - Get the DM representing the reference cell 545 546 Not collective 547 548 Input Parameters: 549 + sp - The PetscDualSpace 550 - dm - The reference cell 551 552 Level: intermediate 553 554 .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate() 555 @*/ 556 PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm) 557 { 558 PetscErrorCode ierr; 559 560 PetscFunctionBegin; 561 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 562 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 563 PetscCheckFalse(sp->setupcalled,PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change DM after dualspace is set up"); 564 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 565 if (sp->dm && sp->dm != dm) { 566 ierr = PetscDualSpaceClearDMData_Internal(sp, sp->dm);CHKERRQ(ierr); 567 } 568 ierr = DMDestroy(&sp->dm);CHKERRQ(ierr); 569 sp->dm = dm; 570 PetscFunctionReturn(0); 571 } 572 573 /*@ 574 PetscDualSpaceGetOrder - Get the order of the dual space 575 576 Not collective 577 578 Input Parameter: 579 . sp - The PetscDualSpace 580 581 Output Parameter: 582 . order - The order 583 584 Level: intermediate 585 586 .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate() 587 @*/ 588 PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order) 589 { 590 PetscFunctionBegin; 591 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 592 PetscValidPointer(order, 2); 593 *order = sp->order; 594 PetscFunctionReturn(0); 595 } 596 597 /*@ 598 PetscDualSpaceSetOrder - Set the order of the dual space 599 600 Not collective 601 602 Input Parameters: 603 + sp - The PetscDualSpace 604 - order - The order 605 606 Level: intermediate 607 608 .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate() 609 @*/ 610 PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order) 611 { 612 PetscFunctionBegin; 613 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 614 PetscCheckFalse(sp->setupcalled,PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change order after dualspace is set up"); 615 sp->order = order; 616 PetscFunctionReturn(0); 617 } 618 619 /*@ 620 PetscDualSpaceGetNumComponents - Return the number of components for this space 621 622 Input Parameter: 623 . sp - The PetscDualSpace 624 625 Output Parameter: 626 . Nc - The number of components 627 628 Note: A vector space, for example, will have d components, where d is the spatial dimension 629 630 Level: intermediate 631 632 .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace 633 @*/ 634 PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc) 635 { 636 PetscFunctionBegin; 637 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 638 PetscValidPointer(Nc, 2); 639 *Nc = sp->Nc; 640 PetscFunctionReturn(0); 641 } 642 643 /*@ 644 PetscDualSpaceSetNumComponents - Set the number of components for this space 645 646 Input Parameters: 647 + sp - The PetscDualSpace 648 - order - The number of components 649 650 Level: intermediate 651 652 .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace 653 @*/ 654 PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc) 655 { 656 PetscFunctionBegin; 657 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 658 PetscCheckFalse(sp->setupcalled,PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up"); 659 sp->Nc = Nc; 660 PetscFunctionReturn(0); 661 } 662 663 /*@ 664 PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space 665 666 Not collective 667 668 Input Parameters: 669 + sp - The PetscDualSpace 670 - i - The basis number 671 672 Output Parameter: 673 . functional - The basis functional 674 675 Level: intermediate 676 677 .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate() 678 @*/ 679 PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional) 680 { 681 PetscInt dim; 682 PetscErrorCode ierr; 683 684 PetscFunctionBegin; 685 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 686 PetscValidPointer(functional, 3); 687 ierr = PetscDualSpaceGetDimension(sp, &dim);CHKERRQ(ierr); 688 PetscCheckFalse((i < 0) || (i >= dim),PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim); 689 *functional = sp->functional[i]; 690 PetscFunctionReturn(0); 691 } 692 693 /*@ 694 PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals 695 696 Not collective 697 698 Input Parameter: 699 . sp - The PetscDualSpace 700 701 Output Parameter: 702 . dim - The dimension 703 704 Level: intermediate 705 706 .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate() 707 @*/ 708 PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim) 709 { 710 PetscErrorCode ierr; 711 712 PetscFunctionBegin; 713 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 714 PetscValidPointer(dim, 2); 715 if (sp->spdim < 0) { 716 PetscSection section; 717 718 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 719 if (section) { 720 ierr = PetscSectionGetStorageSize(section, &(sp->spdim));CHKERRQ(ierr); 721 } else sp->spdim = 0; 722 } 723 *dim = sp->spdim; 724 PetscFunctionReturn(0); 725 } 726 727 /*@ 728 PetscDualSpaceGetInteriorDimension - Get the interior dimension of the dual space, i.e. the number of basis functionals assigned to the interior of the reference domain 729 730 Not collective 731 732 Input Parameter: 733 . sp - The PetscDualSpace 734 735 Output Parameter: 736 . dim - The dimension 737 738 Level: intermediate 739 740 .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate() 741 @*/ 742 PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim) 743 { 744 PetscErrorCode ierr; 745 746 PetscFunctionBegin; 747 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 748 PetscValidPointer(intdim, 2); 749 if (sp->spintdim < 0) { 750 PetscSection section; 751 752 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 753 if (section) { 754 ierr = PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim));CHKERRQ(ierr); 755 } else sp->spintdim = 0; 756 } 757 *intdim = sp->spintdim; 758 PetscFunctionReturn(0); 759 } 760 761 /*@ 762 PetscDualSpaceGetUniform - Whether this dual space is uniform 763 764 Not collective 765 766 Input Parameters: 767 . sp - A dual space 768 769 Output Parameters: 770 . uniform - PETSC_TRUE if (a) the dual space is the same for each point in a stratum of the reference DMPlex, and 771 (b) every symmetry of each point in the reference DMPlex is also a symmetry of the point's dual space. 772 773 Level: advanced 774 775 Note: all of the usual spaces on simplex or tensor-product elements will be uniform, only reference cells 776 with non-uniform strata (like trianguar-prisms) or anisotropic hp dual spaces will not be uniform. 777 778 .seealso: PetscDualSpaceGetPointSubspace(), PetscDualSpaceGetSymmetries() 779 @*/ 780 PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform) 781 { 782 PetscFunctionBegin; 783 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 784 PetscValidPointer(uniform, 2); 785 *uniform = sp->uniform; 786 PetscFunctionReturn(0); 787 } 788 789 /*@C 790 PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension 791 792 Not collective 793 794 Input Parameter: 795 . sp - The PetscDualSpace 796 797 Output Parameter: 798 . numDof - An array of length dim+1 which holds the number of dofs for each dimension 799 800 Level: intermediate 801 802 .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate() 803 @*/ 804 PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof) 805 { 806 PetscErrorCode ierr; 807 808 PetscFunctionBegin; 809 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 810 PetscValidPointer(numDof, 2); 811 PetscCheckFalse(!sp->uniform,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform space does not have a fixed number of dofs for each height"); 812 if (!sp->numDof) { 813 DM dm; 814 PetscInt depth, d; 815 PetscSection section; 816 817 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 818 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 819 ierr = PetscCalloc1(depth+1,&(sp->numDof));CHKERRQ(ierr); 820 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 821 for (d = 0; d <= depth; d++) { 822 PetscInt dStart, dEnd; 823 824 ierr = DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);CHKERRQ(ierr); 825 if (dEnd <= dStart) continue; 826 ierr = PetscSectionGetDof(section, dStart, &(sp->numDof[d]));CHKERRQ(ierr); 827 828 } 829 } 830 *numDof = sp->numDof; 831 PetscCheckFalse(!*numDof,PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation"); 832 PetscFunctionReturn(0); 833 } 834 835 /* create the section of the right size and set a permutation for topological ordering */ 836 PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection) 837 { 838 DM dm; 839 PetscInt pStart, pEnd, cStart, cEnd, c, depth, count, i; 840 PetscInt *seen, *perm; 841 PetscSection section; 842 PetscErrorCode ierr; 843 844 PetscFunctionBegin; 845 dm = sp->dm; 846 ierr = PetscSectionCreate(PETSC_COMM_SELF, §ion);CHKERRQ(ierr); 847 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 848 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 849 ierr = PetscCalloc1(pEnd - pStart, &seen);CHKERRQ(ierr); 850 ierr = PetscMalloc1(pEnd - pStart, &perm);CHKERRQ(ierr); 851 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 852 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 853 for (c = cStart, count = 0; c < cEnd; c++) { 854 PetscInt closureSize = -1, e; 855 PetscInt *closure = NULL; 856 857 perm[count++] = c; 858 seen[c-pStart] = 1; 859 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 860 for (e = 0; e < closureSize; e++) { 861 PetscInt point = closure[2*e]; 862 863 if (seen[point-pStart]) continue; 864 perm[count++] = point; 865 seen[point-pStart] = 1; 866 } 867 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 868 } 869 PetscCheckFalse(count != pEnd - pStart,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad topological ordering"); 870 for (i = 0; i < pEnd - pStart; i++) if (perm[i] != i) break; 871 if (i < pEnd - pStart) { 872 IS permIS; 873 874 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS);CHKERRQ(ierr); 875 ierr = ISSetPermutation(permIS);CHKERRQ(ierr); 876 ierr = PetscSectionSetPermutation(section, permIS);CHKERRQ(ierr); 877 ierr = ISDestroy(&permIS);CHKERRQ(ierr); 878 } else { 879 ierr = PetscFree(perm);CHKERRQ(ierr); 880 } 881 ierr = PetscFree(seen);CHKERRQ(ierr); 882 *topSection = section; 883 PetscFunctionReturn(0); 884 } 885 886 /* mark boundary points and set up */ 887 PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section) 888 { 889 DM dm; 890 DMLabel boundary; 891 PetscInt pStart, pEnd, p; 892 PetscErrorCode ierr; 893 894 PetscFunctionBegin; 895 dm = sp->dm; 896 ierr = DMLabelCreate(PETSC_COMM_SELF,"boundary",&boundary);CHKERRQ(ierr); 897 ierr = PetscDualSpaceGetDM(sp,&dm);CHKERRQ(ierr); 898 ierr = DMPlexMarkBoundaryFaces(dm,1,boundary);CHKERRQ(ierr); 899 ierr = DMPlexLabelComplete(dm,boundary);CHKERRQ(ierr); 900 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 901 for (p = pStart; p < pEnd; p++) { 902 PetscInt bval; 903 904 ierr = DMLabelGetValue(boundary, p, &bval);CHKERRQ(ierr); 905 if (bval == 1) { 906 PetscInt dof; 907 908 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 909 ierr = PetscSectionSetConstraintDof(section, p, dof);CHKERRQ(ierr); 910 } 911 } 912 ierr = DMLabelDestroy(&boundary);CHKERRQ(ierr); 913 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 /*@ 918 PetscDualSpaceGetSection - Create a PetscSection over the reference cell with the layout from this space 919 920 Collective on sp 921 922 Input Parameters: 923 . sp - The PetscDualSpace 924 925 Output Parameter: 926 . section - The section 927 928 Level: advanced 929 930 .seealso: PetscDualSpaceCreate(), DMPLEX 931 @*/ 932 PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section) 933 { 934 PetscInt pStart, pEnd, p; 935 PetscErrorCode ierr; 936 937 PetscFunctionBegin; 938 if (!sp->pointSection) { 939 /* mark the boundary */ 940 ierr = PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection));CHKERRQ(ierr); 941 ierr = DMPlexGetChart(sp->dm,&pStart,&pEnd);CHKERRQ(ierr); 942 for (p = pStart; p < pEnd; p++) { 943 PetscDualSpace psp; 944 945 ierr = PetscDualSpaceGetPointSubspace(sp, p, &psp);CHKERRQ(ierr); 946 if (psp) { 947 PetscInt dof; 948 949 ierr = PetscDualSpaceGetInteriorDimension(psp, &dof);CHKERRQ(ierr); 950 ierr = PetscSectionSetDof(sp->pointSection,p,dof);CHKERRQ(ierr); 951 } 952 } 953 ierr = PetscDualSpaceSectionSetUp_Internal(sp,sp->pointSection);CHKERRQ(ierr); 954 } 955 *section = sp->pointSection; 956 PetscFunctionReturn(0); 957 } 958 959 /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs 960 * have one cell */ 961 PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd) 962 { 963 PetscReal *sv0, *v0, *J; 964 PetscSection section; 965 PetscInt dim, s, k; 966 DM dm; 967 PetscErrorCode ierr; 968 969 PetscFunctionBegin; 970 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 971 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 972 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 973 ierr = PetscMalloc3(dim, &v0, dim, &sv0, dim*dim, &J);CHKERRQ(ierr); 974 ierr = PetscDualSpaceGetFormDegree(sp, &k);CHKERRQ(ierr); 975 for (s = sStart; s < sEnd; s++) { 976 PetscReal detJ, hdetJ; 977 PetscDualSpace ssp; 978 PetscInt dof, off, f, sdim; 979 PetscInt i, j; 980 DM sdm; 981 982 ierr = PetscDualSpaceGetPointSubspace(sp, s, &ssp);CHKERRQ(ierr); 983 if (!ssp) continue; 984 ierr = PetscSectionGetDof(section, s, &dof);CHKERRQ(ierr); 985 ierr = PetscSectionGetOffset(section, s, &off);CHKERRQ(ierr); 986 /* get the first vertex of the reference cell */ 987 ierr = PetscDualSpaceGetDM(ssp, &sdm);CHKERRQ(ierr); 988 ierr = DMGetDimension(sdm, &sdim);CHKERRQ(ierr); 989 ierr = DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ);CHKERRQ(ierr); 990 ierr = DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ);CHKERRQ(ierr); 991 /* compactify Jacobian */ 992 for (i = 0; i < dim; i++) for (j = 0; j < sdim; j++) J[i* sdim + j] = J[i * dim + j]; 993 for (f = 0; f < dof; f++) { 994 PetscQuadrature fn; 995 996 ierr = PetscDualSpaceGetFunctional(ssp, f, &fn);CHKERRQ(ierr); 997 ierr = PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off+f]));CHKERRQ(ierr); 998 } 999 } 1000 ierr = PetscFree3(v0, sv0, J);CHKERRQ(ierr); 1001 PetscFunctionReturn(0); 1002 } 1003 1004 /*@C 1005 PetscDualSpaceApply - Apply a functional from the dual space basis to an input function 1006 1007 Input Parameters: 1008 + sp - The PetscDualSpace object 1009 . f - The basis functional index 1010 . time - The time 1011 . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian) (or evaluated at the coordinates of the functional) 1012 . numComp - The number of components for the function 1013 . func - The input function 1014 - ctx - A context for the function 1015 1016 Output Parameter: 1017 . value - numComp output values 1018 1019 Note: The calling sequence for the callback func is given by: 1020 1021 $ func(PetscInt dim, PetscReal time, const PetscReal x[], 1022 $ PetscInt numComponents, PetscScalar values[], void *ctx) 1023 1024 Level: beginner 1025 1026 .seealso: PetscDualSpaceCreate() 1027 @*/ 1028 PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value) 1029 { 1030 PetscErrorCode ierr; 1031 1032 PetscFunctionBegin; 1033 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1034 PetscValidPointer(cgeom, 4); 1035 PetscValidPointer(value, 8); 1036 ierr = (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);CHKERRQ(ierr); 1037 PetscFunctionReturn(0); 1038 } 1039 1040 /*@C 1041 PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData() 1042 1043 Input Parameters: 1044 + sp - The PetscDualSpace object 1045 - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData() 1046 1047 Output Parameter: 1048 . spValue - The values of all dual space functionals 1049 1050 Level: beginner 1051 1052 .seealso: PetscDualSpaceCreate() 1053 @*/ 1054 PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue) 1055 { 1056 PetscErrorCode ierr; 1057 1058 PetscFunctionBegin; 1059 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1060 ierr = (*sp->ops->applyall)(sp, pointEval, spValue);CHKERRQ(ierr); 1061 PetscFunctionReturn(0); 1062 } 1063 1064 /*@C 1065 PetscDualSpaceApplyInterior - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData() 1066 1067 Input Parameters: 1068 + sp - The PetscDualSpace object 1069 - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData() 1070 1071 Output Parameter: 1072 . spValue - The values of interior dual space functionals 1073 1074 Level: beginner 1075 1076 .seealso: PetscDualSpaceCreate() 1077 @*/ 1078 PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue) 1079 { 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1084 ierr = (*sp->ops->applyint)(sp, pointEval, spValue);CHKERRQ(ierr); 1085 PetscFunctionReturn(0); 1086 } 1087 1088 /*@C 1089 PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional. 1090 1091 Input Parameters: 1092 + sp - The PetscDualSpace object 1093 . f - The basis functional index 1094 . time - The time 1095 . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian) 1096 . Nc - The number of components for the function 1097 . func - The input function 1098 - ctx - A context for the function 1099 1100 Output Parameter: 1101 . value - The output value 1102 1103 Note: The calling sequence for the callback func is given by: 1104 1105 $ func(PetscInt dim, PetscReal time, const PetscReal x[], 1106 $ PetscInt numComponents, PetscScalar values[], void *ctx) 1107 1108 and the idea is to evaluate the functional as an integral 1109 1110 $ n(f) = int dx n(x) . f(x) 1111 1112 where both n and f have Nc components. 1113 1114 Level: beginner 1115 1116 .seealso: PetscDualSpaceCreate() 1117 @*/ 1118 PetscErrorCode PetscDualSpaceApplyDefault(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value) 1119 { 1120 DM dm; 1121 PetscQuadrature n; 1122 const PetscReal *points, *weights; 1123 PetscReal x[3]; 1124 PetscScalar *val; 1125 PetscInt dim, dE, qNc, c, Nq, q; 1126 PetscBool isAffine; 1127 PetscErrorCode ierr; 1128 1129 PetscFunctionBegin; 1130 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1131 PetscValidPointer(value, 8); 1132 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1133 ierr = PetscDualSpaceGetFunctional(sp, f, &n);CHKERRQ(ierr); 1134 ierr = PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);CHKERRQ(ierr); 1135 PetscCheckFalse(dim != cgeom->dim,PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim); 1136 PetscCheckFalse(qNc != Nc,PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc); 1137 ierr = DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);CHKERRQ(ierr); 1138 *value = 0.0; 1139 isAffine = cgeom->isAffine; 1140 dE = cgeom->dimEmbed; 1141 for (q = 0; q < Nq; ++q) { 1142 if (isAffine) { 1143 CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x); 1144 ierr = (*func)(dE, time, x, Nc, val, ctx);CHKERRQ(ierr); 1145 } else { 1146 ierr = (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);CHKERRQ(ierr); 1147 } 1148 for (c = 0; c < Nc; ++c) { 1149 *value += val[c]*weights[q*Nc+c]; 1150 } 1151 } 1152 ierr = DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);CHKERRQ(ierr); 1153 PetscFunctionReturn(0); 1154 } 1155 1156 /*@C 1157 PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllData() 1158 1159 Input Parameters: 1160 + sp - The PetscDualSpace object 1161 - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData() 1162 1163 Output Parameter: 1164 . spValue - The values of all dual space functionals 1165 1166 Level: beginner 1167 1168 .seealso: PetscDualSpaceCreate() 1169 @*/ 1170 PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue) 1171 { 1172 Vec pointValues, dofValues; 1173 Mat allMat; 1174 PetscErrorCode ierr; 1175 1176 PetscFunctionBegin; 1177 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1178 PetscValidScalarPointer(pointEval, 2); 1179 PetscValidScalarPointer(spValue, 3); 1180 ierr = PetscDualSpaceGetAllData(sp, NULL, &allMat);CHKERRQ(ierr); 1181 if (!(sp->allNodeValues)) { 1182 ierr = MatCreateVecs(allMat, &(sp->allNodeValues), NULL);CHKERRQ(ierr); 1183 } 1184 pointValues = sp->allNodeValues; 1185 if (!(sp->allDofValues)) { 1186 ierr = MatCreateVecs(allMat, NULL, &(sp->allDofValues));CHKERRQ(ierr); 1187 } 1188 dofValues = sp->allDofValues; 1189 ierr = VecPlaceArray(pointValues, pointEval);CHKERRQ(ierr); 1190 ierr = VecPlaceArray(dofValues, spValue);CHKERRQ(ierr); 1191 ierr = MatMult(allMat, pointValues, dofValues);CHKERRQ(ierr); 1192 ierr = VecResetArray(dofValues);CHKERRQ(ierr); 1193 ierr = VecResetArray(pointValues);CHKERRQ(ierr); 1194 PetscFunctionReturn(0); 1195 } 1196 1197 /*@C 1198 PetscDualSpaceApplyInteriorDefault - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData() 1199 1200 Input Parameters: 1201 + sp - The PetscDualSpace object 1202 - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData() 1203 1204 Output Parameter: 1205 . spValue - The values of interior dual space functionals 1206 1207 Level: beginner 1208 1209 .seealso: PetscDualSpaceCreate() 1210 @*/ 1211 PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue) 1212 { 1213 Vec pointValues, dofValues; 1214 Mat intMat; 1215 PetscErrorCode ierr; 1216 1217 PetscFunctionBegin; 1218 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1219 PetscValidScalarPointer(pointEval, 2); 1220 PetscValidScalarPointer(spValue, 3); 1221 ierr = PetscDualSpaceGetInteriorData(sp, NULL, &intMat);CHKERRQ(ierr); 1222 if (!(sp->intNodeValues)) { 1223 ierr = MatCreateVecs(intMat, &(sp->intNodeValues), NULL);CHKERRQ(ierr); 1224 } 1225 pointValues = sp->intNodeValues; 1226 if (!(sp->intDofValues)) { 1227 ierr = MatCreateVecs(intMat, NULL, &(sp->intDofValues));CHKERRQ(ierr); 1228 } 1229 dofValues = sp->intDofValues; 1230 ierr = VecPlaceArray(pointValues, pointEval);CHKERRQ(ierr); 1231 ierr = VecPlaceArray(dofValues, spValue);CHKERRQ(ierr); 1232 ierr = MatMult(intMat, pointValues, dofValues);CHKERRQ(ierr); 1233 ierr = VecResetArray(dofValues);CHKERRQ(ierr); 1234 ierr = VecResetArray(pointValues);CHKERRQ(ierr); 1235 PetscFunctionReturn(0); 1236 } 1237 1238 /*@ 1239 PetscDualSpaceGetAllData - Get all quadrature nodes from this space, and the matrix that sends quadrature node values to degree-of-freedom values 1240 1241 Input Parameter: 1242 . sp - The dualspace 1243 1244 Output Parameters: 1245 + allNodes - A PetscQuadrature object containing all evaluation nodes 1246 - allMat - A Mat for the node-to-dof transformation 1247 1248 Level: advanced 1249 1250 .seealso: PetscDualSpaceCreate() 1251 @*/ 1252 PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat) 1253 { 1254 PetscErrorCode ierr; 1255 1256 PetscFunctionBegin; 1257 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1258 if (allNodes) PetscValidPointer(allNodes,2); 1259 if (allMat) PetscValidPointer(allMat,3); 1260 if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) { 1261 PetscQuadrature qpoints; 1262 Mat amat; 1263 1264 ierr = (*sp->ops->createalldata)(sp,&qpoints,&amat);CHKERRQ(ierr); 1265 ierr = PetscQuadratureDestroy(&(sp->allNodes));CHKERRQ(ierr); 1266 ierr = MatDestroy(&(sp->allMat));CHKERRQ(ierr); 1267 sp->allNodes = qpoints; 1268 sp->allMat = amat; 1269 } 1270 if (allNodes) *allNodes = sp->allNodes; 1271 if (allMat) *allMat = sp->allMat; 1272 PetscFunctionReturn(0); 1273 } 1274 1275 /*@ 1276 PetscDualSpaceCreateAllDataDefault - Create all evaluation nodes and the node-to-dof matrix by examining functionals 1277 1278 Input Parameter: 1279 . sp - The dualspace 1280 1281 Output Parameters: 1282 + allNodes - A PetscQuadrature object containing all evaluation nodes 1283 - allMat - A Mat for the node-to-dof transformation 1284 1285 Level: advanced 1286 1287 .seealso: PetscDualSpaceCreate() 1288 @*/ 1289 PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat) 1290 { 1291 PetscInt spdim; 1292 PetscInt numPoints, offset; 1293 PetscReal *points; 1294 PetscInt f, dim; 1295 PetscInt Nc, nrows, ncols; 1296 PetscInt maxNumPoints; 1297 PetscQuadrature q; 1298 Mat A; 1299 PetscErrorCode ierr; 1300 1301 PetscFunctionBegin; 1302 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 1303 ierr = PetscDualSpaceGetDimension(sp,&spdim);CHKERRQ(ierr); 1304 if (!spdim) { 1305 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);CHKERRQ(ierr); 1306 ierr = PetscQuadratureSetData(*allNodes,0,0,0,NULL,NULL);CHKERRQ(ierr); 1307 } 1308 nrows = spdim; 1309 ierr = PetscDualSpaceGetFunctional(sp,0,&q);CHKERRQ(ierr); 1310 ierr = PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);CHKERRQ(ierr); 1311 maxNumPoints = numPoints; 1312 for (f = 1; f < spdim; f++) { 1313 PetscInt Np; 1314 1315 ierr = PetscDualSpaceGetFunctional(sp,f,&q);CHKERRQ(ierr); 1316 ierr = PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);CHKERRQ(ierr); 1317 numPoints += Np; 1318 maxNumPoints = PetscMax(maxNumPoints,Np); 1319 } 1320 ncols = numPoints * Nc; 1321 ierr = PetscMalloc1(dim*numPoints,&points);CHKERRQ(ierr); 1322 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A);CHKERRQ(ierr); 1323 for (f = 0, offset = 0; f < spdim; f++) { 1324 const PetscReal *p, *w; 1325 PetscInt Np, i; 1326 PetscInt fnc; 1327 1328 ierr = PetscDualSpaceGetFunctional(sp,f,&q);CHKERRQ(ierr); 1329 ierr = PetscQuadratureGetData(q,NULL,&fnc,&Np,&p,&w);CHKERRQ(ierr); 1330 PetscCheckFalse(fnc != Nc,PETSC_COMM_SELF, PETSC_ERR_PLIB, "functional component mismatch"); 1331 for (i = 0; i < Np * dim; i++) { 1332 points[offset* dim + i] = p[i]; 1333 } 1334 for (i = 0; i < Np * Nc; i++) { 1335 ierr = MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES);CHKERRQ(ierr); 1336 } 1337 offset += Np; 1338 } 1339 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1340 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1341 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);CHKERRQ(ierr); 1342 ierr = PetscQuadratureSetData(*allNodes,dim,0,numPoints,points,NULL);CHKERRQ(ierr); 1343 *allMat = A; 1344 PetscFunctionReturn(0); 1345 } 1346 1347 /*@ 1348 PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from 1349 this space, as well as the matrix that computes the degrees of freedom from the quadrature values. Degrees of 1350 freedom are interior degrees of freedom if they belong (by PetscDualSpaceGetSection()) to interior points in the 1351 reference DMPlex: complementary boundary degrees of freedom are marked as constrained in the section returned by 1352 PetscDualSpaceGetSection()). 1353 1354 Input Parameter: 1355 . sp - The dualspace 1356 1357 Output Parameters: 1358 + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom 1359 - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is 1360 the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section, 1361 npoints is the number of points in intNodes and nc is PetscDualSpaceGetNumComponents(). 1362 1363 Level: advanced 1364 1365 .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetDimension(), PetscDualSpaceGetNumComponents(), PetscQuadratureGetData() 1366 @*/ 1367 PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat) 1368 { 1369 PetscErrorCode ierr; 1370 1371 PetscFunctionBegin; 1372 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1373 if (intNodes) PetscValidPointer(intNodes,2); 1374 if (intMat) PetscValidPointer(intMat,3); 1375 if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) { 1376 PetscQuadrature qpoints; 1377 Mat imat; 1378 1379 ierr = (*sp->ops->createintdata)(sp,&qpoints,&imat);CHKERRQ(ierr); 1380 ierr = PetscQuadratureDestroy(&(sp->intNodes));CHKERRQ(ierr); 1381 ierr = MatDestroy(&(sp->intMat));CHKERRQ(ierr); 1382 sp->intNodes = qpoints; 1383 sp->intMat = imat; 1384 } 1385 if (intNodes) *intNodes = sp->intNodes; 1386 if (intMat) *intMat = sp->intMat; 1387 PetscFunctionReturn(0); 1388 } 1389 1390 /*@ 1391 PetscDualSpaceCreateInteriorDataDefault - Create quadrature points by examining interior functionals and create the matrix mapping quadrature point values to interior dual space values 1392 1393 Input Parameter: 1394 . sp - The dualspace 1395 1396 Output Parameters: 1397 + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom 1398 - intMat - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is 1399 the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section, 1400 npoints is the number of points in allNodes and nc is PetscDualSpaceGetNumComponents(). 1401 1402 Level: advanced 1403 1404 .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetInteriorData() 1405 @*/ 1406 PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat) 1407 { 1408 DM dm; 1409 PetscInt spdim0; 1410 PetscInt Nc; 1411 PetscInt pStart, pEnd, p, f; 1412 PetscSection section; 1413 PetscInt numPoints, offset, matoffset; 1414 PetscReal *points; 1415 PetscInt dim; 1416 PetscInt *nnz; 1417 PetscQuadrature q; 1418 Mat imat; 1419 PetscErrorCode ierr; 1420 1421 PetscFunctionBegin; 1422 PetscValidHeaderSpecific(sp,PETSCDUALSPACE_CLASSID,1); 1423 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 1424 ierr = PetscSectionGetConstrainedStorageSize(section, &spdim0);CHKERRQ(ierr); 1425 if (!spdim0) { 1426 *intNodes = NULL; 1427 *intMat = NULL; 1428 PetscFunctionReturn(0); 1429 } 1430 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 1431 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 1432 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1433 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1434 ierr = PetscMalloc1(spdim0, &nnz);CHKERRQ(ierr); 1435 for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) { 1436 PetscInt dof, cdof, off, d; 1437 1438 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 1439 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr); 1440 if (!(dof - cdof)) continue; 1441 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 1442 for (d = 0; d < dof; d++, off++, f++) { 1443 PetscInt Np; 1444 1445 ierr = PetscDualSpaceGetFunctional(sp,off,&q);CHKERRQ(ierr); 1446 ierr = PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);CHKERRQ(ierr); 1447 nnz[f] = Np * Nc; 1448 numPoints += Np; 1449 } 1450 } 1451 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat);CHKERRQ(ierr); 1452 ierr = PetscFree(nnz);CHKERRQ(ierr); 1453 ierr = PetscMalloc1(dim*numPoints,&points);CHKERRQ(ierr); 1454 for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) { 1455 PetscInt dof, cdof, off, d; 1456 1457 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 1458 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr); 1459 if (!(dof - cdof)) continue; 1460 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 1461 for (d = 0; d < dof; d++, off++, f++) { 1462 const PetscReal *p; 1463 const PetscReal *w; 1464 PetscInt Np, i; 1465 1466 ierr = PetscDualSpaceGetFunctional(sp,off,&q);CHKERRQ(ierr); 1467 ierr = PetscQuadratureGetData(q,NULL,NULL,&Np,&p,&w);CHKERRQ(ierr); 1468 for (i = 0; i < Np * dim; i++) { 1469 points[offset + i] = p[i]; 1470 } 1471 for (i = 0; i < Np * Nc; i++) { 1472 ierr = MatSetValue(imat, f, matoffset + i, w[i],INSERT_VALUES);CHKERRQ(ierr); 1473 } 1474 offset += Np * dim; 1475 matoffset += Np * Nc; 1476 } 1477 } 1478 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,intNodes);CHKERRQ(ierr); 1479 ierr = PetscQuadratureSetData(*intNodes,dim,0,numPoints,points,NULL);CHKERRQ(ierr); 1480 ierr = MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1481 ierr = MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1482 *intMat = imat; 1483 PetscFunctionReturn(0); 1484 } 1485 1486 /*@C 1487 PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid. 1488 1489 Input Parameters: 1490 + sp - The PetscDualSpace object 1491 . f - The basis functional index 1492 . time - The time 1493 . cgeom - A context with geometric information for this cell, we currently just use the centroid 1494 . Nc - The number of components for the function 1495 . func - The input function 1496 - ctx - A context for the function 1497 1498 Output Parameter: 1499 . value - The output value (scalar) 1500 1501 Note: The calling sequence for the callback func is given by: 1502 1503 $ func(PetscInt dim, PetscReal time, const PetscReal x[], 1504 $ PetscInt numComponents, PetscScalar values[], void *ctx) 1505 1506 and the idea is to evaluate the functional as an integral 1507 1508 $ n(f) = int dx n(x) . f(x) 1509 1510 where both n and f have Nc components. 1511 1512 Level: beginner 1513 1514 .seealso: PetscDualSpaceCreate() 1515 @*/ 1516 PetscErrorCode PetscDualSpaceApplyFVM(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFVCellGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value) 1517 { 1518 DM dm; 1519 PetscQuadrature n; 1520 const PetscReal *points, *weights; 1521 PetscScalar *val; 1522 PetscInt dimEmbed, qNc, c, Nq, q; 1523 PetscErrorCode ierr; 1524 1525 PetscFunctionBegin; 1526 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1527 PetscValidPointer(value, 8); 1528 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1529 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1530 ierr = PetscDualSpaceGetFunctional(sp, f, &n);CHKERRQ(ierr); 1531 ierr = PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);CHKERRQ(ierr); 1532 PetscCheckFalse(qNc != Nc,PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc); 1533 ierr = DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);CHKERRQ(ierr); 1534 *value = 0.; 1535 for (q = 0; q < Nq; ++q) { 1536 ierr = (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);CHKERRQ(ierr); 1537 for (c = 0; c < Nc; ++c) { 1538 *value += val[c]*weights[q*Nc+c]; 1539 } 1540 } 1541 ierr = DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);CHKERRQ(ierr); 1542 PetscFunctionReturn(0); 1543 } 1544 1545 /*@ 1546 PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a 1547 given height. This assumes that the reference cell is symmetric over points of this height. 1548 1549 If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and 1550 pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not 1551 support extracting subspaces, then NULL is returned. 1552 1553 This does not increment the reference count on the returned dual space, and the user should not destroy it. 1554 1555 Not collective 1556 1557 Input Parameters: 1558 + sp - the PetscDualSpace object 1559 - height - the height of the mesh point for which the subspace is desired 1560 1561 Output Parameter: 1562 . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the 1563 point, which will be of lesser dimension if height > 0. 1564 1565 Level: advanced 1566 1567 .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace 1568 @*/ 1569 PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp) 1570 { 1571 PetscInt depth = -1, cStart, cEnd; 1572 DM dm; 1573 PetscErrorCode ierr; 1574 1575 PetscFunctionBegin; 1576 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1577 PetscValidPointer(subsp,3); 1578 PetscCheckFalse(!(sp->uniform),PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform dual space does not have a single dual space at each height"); 1579 *subsp = NULL; 1580 dm = sp->dm; 1581 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1582 PetscCheckFalse(height < 0 || height > depth,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid height"); 1583 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1584 if (height == 0 && cEnd == cStart + 1) { 1585 *subsp = sp; 1586 PetscFunctionReturn(0); 1587 } 1588 if (!sp->heightSpaces) { 1589 PetscInt h; 1590 ierr = PetscCalloc1(depth+1, &(sp->heightSpaces));CHKERRQ(ierr); 1591 1592 for (h = 0; h <= depth; h++) { 1593 if (h == 0 && cEnd == cStart + 1) continue; 1594 if (sp->ops->createheightsubspace) {ierr = (*sp->ops->createheightsubspace)(sp,height,&(sp->heightSpaces[h]));CHKERRQ(ierr);} 1595 else if (sp->pointSpaces) { 1596 PetscInt hStart, hEnd; 1597 1598 ierr = DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 1599 if (hEnd > hStart) { 1600 const char *name; 1601 1602 ierr = PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));CHKERRQ(ierr); 1603 if (sp->pointSpaces[hStart]) { 1604 ierr = PetscObjectGetName((PetscObject) sp, &name);CHKERRQ(ierr); 1605 ierr = PetscObjectSetName((PetscObject) sp->pointSpaces[hStart], name);CHKERRQ(ierr); 1606 } 1607 sp->heightSpaces[h] = sp->pointSpaces[hStart]; 1608 } 1609 } 1610 } 1611 } 1612 *subsp = sp->heightSpaces[height]; 1613 PetscFunctionReturn(0); 1614 } 1615 1616 /*@ 1617 PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point. 1618 1619 If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not 1620 defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting 1621 subspaces, then NULL is returned. 1622 1623 This does not increment the reference count on the returned dual space, and the user should not destroy it. 1624 1625 Not collective 1626 1627 Input Parameters: 1628 + sp - the PetscDualSpace object 1629 - point - the point (in the dual space's DM) for which the subspace is desired 1630 1631 Output Parameters: 1632 bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the 1633 point, which will be of lesser dimension if height > 0. 1634 1635 Level: advanced 1636 1637 .seealso: PetscDualSpace 1638 @*/ 1639 PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp) 1640 { 1641 PetscInt pStart = 0, pEnd = 0, cStart, cEnd; 1642 DM dm; 1643 PetscErrorCode ierr; 1644 1645 PetscFunctionBegin; 1646 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1647 PetscValidPointer(bdsp,3); 1648 *bdsp = NULL; 1649 dm = sp->dm; 1650 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1651 PetscCheckFalse(point < pStart || point > pEnd,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point"); 1652 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1653 if (point == cStart && cEnd == cStart + 1) { /* the dual space is only equivalent to the dual space on a cell if the reference mesh has just one cell */ 1654 *bdsp = sp; 1655 PetscFunctionReturn(0); 1656 } 1657 if (!sp->pointSpaces) { 1658 PetscInt p; 1659 ierr = PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));CHKERRQ(ierr); 1660 1661 for (p = 0; p < pEnd - pStart; p++) { 1662 if (p + pStart == cStart && cEnd == cStart + 1) continue; 1663 if (sp->ops->createpointsubspace) {ierr = (*sp->ops->createpointsubspace)(sp,p+pStart,&(sp->pointSpaces[p]));CHKERRQ(ierr);} 1664 else if (sp->heightSpaces || sp->ops->createheightsubspace) { 1665 PetscInt dim, depth, height; 1666 DMLabel label; 1667 1668 ierr = DMPlexGetDepth(dm,&dim);CHKERRQ(ierr); 1669 ierr = DMPlexGetDepthLabel(dm,&label);CHKERRQ(ierr); 1670 ierr = DMLabelGetValue(label,p+pStart,&depth);CHKERRQ(ierr); 1671 height = dim - depth; 1672 ierr = PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));CHKERRQ(ierr); 1673 ierr = PetscObjectReference((PetscObject)sp->pointSpaces[p]);CHKERRQ(ierr); 1674 } 1675 } 1676 } 1677 *bdsp = sp->pointSpaces[point - pStart]; 1678 PetscFunctionReturn(0); 1679 } 1680 1681 /*@C 1682 PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis 1683 1684 Not collective 1685 1686 Input Parameter: 1687 . sp - the PetscDualSpace object 1688 1689 Output Parameters: 1690 + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation 1691 - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation 1692 1693 Note: The permutation and flip arrays are organized in the following way 1694 $ perms[p][ornt][dof # on point] = new local dof # 1695 $ flips[p][ornt][dof # on point] = reversal or not 1696 1697 Level: developer 1698 1699 @*/ 1700 PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips) 1701 { 1702 PetscErrorCode ierr; 1703 1704 PetscFunctionBegin; 1705 PetscValidHeaderSpecific(sp,PETSCDUALSPACE_CLASSID,1); 1706 if (perms) {PetscValidPointer(perms,2); *perms = NULL;} 1707 if (flips) {PetscValidPointer(flips,3); *flips = NULL;} 1708 if (sp->ops->getsymmetries) {ierr = (sp->ops->getsymmetries)(sp,perms,flips);CHKERRQ(ierr);} 1709 PetscFunctionReturn(0); 1710 } 1711 1712 /*@ 1713 PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this 1714 dual space's functionals. 1715 1716 Input Parameter: 1717 . dsp - The PetscDualSpace 1718 1719 Output Parameter: 1720 . k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored 1721 in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example, 1722 the 1-form basis in 3-D is (dx, dy, dz), and the 2-form basis in 3-D is (dx wedge dy, dx wedge dz, dy wedge dz). 1723 If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the 1724 Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms 1725 but are stored as 1-forms. 1726 1727 Level: developer 1728 1729 .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType 1730 @*/ 1731 PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k) 1732 { 1733 PetscFunctionBeginHot; 1734 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1735 PetscValidPointer(k, 2); 1736 *k = dsp->k; 1737 PetscFunctionReturn(0); 1738 } 1739 1740 /*@ 1741 PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this 1742 dual space's functionals. 1743 1744 Input Parameters: 1745 + dsp - The PetscDualSpace 1746 - k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored 1747 in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example, 1748 the 1-form basis in 3-D is (dx, dy, dz), and the 2-form basis in 3-D is (dx wedge dy, dx wedge dz, dy wedge dz). 1749 If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the 1750 Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms 1751 but are stored as 1-forms. 1752 1753 Level: developer 1754 1755 .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType 1756 @*/ 1757 PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k) 1758 { 1759 PetscInt dim; 1760 1761 PetscFunctionBeginHot; 1762 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1763 PetscCheckFalse(dsp->setupcalled,PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up"); 1764 dim = dsp->dm->dim; 1765 PetscCheckFalse(k < -dim || k > dim,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported %D-form on %D-dimensional reference cell", PetscAbsInt(k), dim); 1766 dsp->k = k; 1767 PetscFunctionReturn(0); 1768 } 1769 1770 /*@ 1771 PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space 1772 1773 Input Parameter: 1774 . dsp - The PetscDualSpace 1775 1776 Output Parameter: 1777 . k - The simplex dimension 1778 1779 Level: developer 1780 1781 Note: Currently supported values are 1782 $ 0: These are H_1 methods that only transform coordinates 1783 $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM) 1784 $ 2: These are the same as 1 1785 $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM) 1786 1787 .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType 1788 @*/ 1789 PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k) 1790 { 1791 PetscInt dim; 1792 1793 PetscFunctionBeginHot; 1794 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1795 PetscValidPointer(k, 2); 1796 dim = dsp->dm->dim; 1797 if (!dsp->k) *k = IDENTITY_TRANSFORM; 1798 else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM; 1799 else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM; 1800 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation"); 1801 PetscFunctionReturn(0); 1802 } 1803 1804 /*@C 1805 PetscDualSpaceTransform - Transform the function values 1806 1807 Input Parameters: 1808 + dsp - The PetscDualSpace 1809 . trans - The type of transform 1810 . isInverse - Flag to invert the transform 1811 . fegeom - The cell geometry 1812 . Nv - The number of function samples 1813 . Nc - The number of function components 1814 - vals - The function values 1815 1816 Output Parameter: 1817 . vals - The transformed function values 1818 1819 Level: intermediate 1820 1821 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1822 1823 .seealso: PetscDualSpaceTransformGradient(), PetscDualSpaceTransformHessian(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType 1824 @*/ 1825 PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[]) 1826 { 1827 PetscReal Jstar[9] = {0}; 1828 PetscInt dim, v, c, Nk; 1829 PetscErrorCode ierr; 1830 1831 PetscFunctionBeginHot; 1832 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1833 PetscValidPointer(fegeom, 4); 1834 PetscValidPointer(vals, 7); 1835 /* TODO: not handling dimEmbed != dim right now */ 1836 dim = dsp->dm->dim; 1837 /* No change needed for 0-forms */ 1838 if (!dsp->k) PetscFunctionReturn(0); 1839 ierr = PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);CHKERRQ(ierr); 1840 /* TODO: use fegeom->isAffine */ 1841 ierr = PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);CHKERRQ(ierr); 1842 for (v = 0; v < Nv; ++v) { 1843 switch (Nk) { 1844 case 1: 1845 for (c = 0; c < Nc; c++) vals[v*Nc + c] *= Jstar[0]; 1846 break; 1847 case 2: 1848 for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]); 1849 break; 1850 case 3: 1851 for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]); 1852 break; 1853 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %D for transformation", Nk); 1854 } 1855 } 1856 PetscFunctionReturn(0); 1857 } 1858 1859 /*@C 1860 PetscDualSpaceTransformGradient - Transform the function gradient values 1861 1862 Input Parameters: 1863 + dsp - The PetscDualSpace 1864 . trans - The type of transform 1865 . isInverse - Flag to invert the transform 1866 . fegeom - The cell geometry 1867 . Nv - The number of function gradient samples 1868 . Nc - The number of function components 1869 - vals - The function gradient values 1870 1871 Output Parameter: 1872 . vals - The transformed function gradient values 1873 1874 Level: intermediate 1875 1876 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1877 1878 .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType 1879 @*/ 1880 PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[]) 1881 { 1882 const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed; 1883 PetscInt v, c, d; 1884 1885 PetscFunctionBeginHot; 1886 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1887 PetscValidPointer(fegeom, 4); 1888 PetscValidPointer(vals, 7); 1889 #ifdef PETSC_USE_DEBUG 1890 PetscCheckFalse(dE <= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE); 1891 #endif 1892 /* Transform gradient */ 1893 if (dim == dE) { 1894 for (v = 0; v < Nv; ++v) { 1895 for (c = 0; c < Nc; ++c) { 1896 switch (dim) 1897 { 1898 case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];break; 1899 case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break; 1900 case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break; 1901 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1902 } 1903 } 1904 } 1905 } else { 1906 for (v = 0; v < Nv; ++v) { 1907 for (c = 0; c < Nc; ++c) { 1908 DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v*Nc+c)*dE], &vals[(v*Nc+c)*dE]); 1909 } 1910 } 1911 } 1912 /* Assume its a vector, otherwise assume its a bunch of scalars */ 1913 if (Nc == 1 || Nc != dim) PetscFunctionReturn(0); 1914 switch (trans) { 1915 case IDENTITY_TRANSFORM: break; 1916 case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */ 1917 if (isInverse) { 1918 for (v = 0; v < Nv; ++v) { 1919 for (d = 0; d < dim; ++d) { 1920 switch (dim) 1921 { 1922 case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1923 case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1924 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1925 } 1926 } 1927 } 1928 } else { 1929 for (v = 0; v < Nv; ++v) { 1930 for (d = 0; d < dim; ++d) { 1931 switch (dim) 1932 { 1933 case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1934 case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1935 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1936 } 1937 } 1938 } 1939 } 1940 break; 1941 case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */ 1942 if (isInverse) { 1943 for (v = 0; v < Nv; ++v) { 1944 for (d = 0; d < dim; ++d) { 1945 switch (dim) 1946 { 1947 case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1948 case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1949 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1950 } 1951 for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0]; 1952 } 1953 } 1954 } else { 1955 for (v = 0; v < Nv; ++v) { 1956 for (d = 0; d < dim; ++d) { 1957 switch (dim) 1958 { 1959 case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1960 case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1961 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1962 } 1963 for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0]; 1964 } 1965 } 1966 } 1967 break; 1968 } 1969 PetscFunctionReturn(0); 1970 } 1971 1972 /*@C 1973 PetscDualSpaceTransformHessian - Transform the function Hessian values 1974 1975 Input Parameters: 1976 + dsp - The PetscDualSpace 1977 . trans - The type of transform 1978 . isInverse - Flag to invert the transform 1979 . fegeom - The cell geometry 1980 . Nv - The number of function Hessian samples 1981 . Nc - The number of function components 1982 - vals - The function gradient values 1983 1984 Output Parameter: 1985 . vals - The transformed function Hessian values 1986 1987 Level: intermediate 1988 1989 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1990 1991 .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType 1992 @*/ 1993 PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[]) 1994 { 1995 const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed; 1996 PetscInt v, c; 1997 1998 PetscFunctionBeginHot; 1999 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2000 PetscValidPointer(fegeom, 4); 2001 PetscValidPointer(vals, 7); 2002 #ifdef PETSC_USE_DEBUG 2003 PetscCheckFalse(dE <= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE); 2004 #endif 2005 /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */ 2006 if (dim == dE) { 2007 for (v = 0; v < Nv; ++v) { 2008 for (c = 0; c < Nc; ++c) { 2009 switch (dim) 2010 { 2011 case 1: vals[(v*Nc+c)*dim*dim] *= PetscSqr(fegeom->invJ[0]);break; 2012 case 2: DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break; 2013 case 3: DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break; 2014 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 2015 } 2016 } 2017 } 2018 } else { 2019 for (v = 0; v < Nv; ++v) { 2020 for (c = 0; c < Nc; ++c) { 2021 DMPlex_PTAPReal_Internal(fegeom->invJ, dim, dE, &vals[(v*Nc+c)*dE*dE], &vals[(v*Nc+c)*dE*dE]); 2022 } 2023 } 2024 } 2025 /* Assume its a vector, otherwise assume its a bunch of scalars */ 2026 if (Nc == 1 || Nc != dim) PetscFunctionReturn(0); 2027 switch (trans) { 2028 case IDENTITY_TRANSFORM: break; 2029 case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */ 2030 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported"); 2031 case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */ 2032 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported"); 2033 } 2034 PetscFunctionReturn(0); 2035 } 2036 2037 /*@C 2038 PetscDualSpacePullback - Transform the given functional so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method. 2039 2040 Input Parameters: 2041 + dsp - The PetscDualSpace 2042 . fegeom - The geometry for this cell 2043 . Nq - The number of function samples 2044 . Nc - The number of function components 2045 - pointEval - The function values 2046 2047 Output Parameter: 2048 . pointEval - The transformed function values 2049 2050 Level: advanced 2051 2052 Note: Functions transform in a complementary way (pushforward) to functionals, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex. 2053 2054 Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2055 2056 .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2057 @*/ 2058 PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2059 { 2060 PetscDualSpaceTransformType trans; 2061 PetscInt k; 2062 PetscErrorCode ierr; 2063 2064 PetscFunctionBeginHot; 2065 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2066 PetscValidPointer(fegeom, 2); 2067 PetscValidPointer(pointEval, 5); 2068 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2069 This determines their transformation properties. */ 2070 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2071 switch (k) 2072 { 2073 case 0: /* H^1 point evaluations */ 2074 trans = IDENTITY_TRANSFORM;break; 2075 case 1: /* Hcurl preserves tangential edge traces */ 2076 trans = COVARIANT_PIOLA_TRANSFORM;break; 2077 case 2: 2078 case 3: /* Hdiv preserve normal traces */ 2079 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2080 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2081 } 2082 ierr = PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2083 PetscFunctionReturn(0); 2084 } 2085 2086 /*@C 2087 PetscDualSpacePushforward - Transform the given function so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method. 2088 2089 Input Parameters: 2090 + dsp - The PetscDualSpace 2091 . fegeom - The geometry for this cell 2092 . Nq - The number of function samples 2093 . Nc - The number of function components 2094 - pointEval - The function values 2095 2096 Output Parameter: 2097 . pointEval - The transformed function values 2098 2099 Level: advanced 2100 2101 Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex. 2102 2103 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2104 2105 .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2106 @*/ 2107 PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2108 { 2109 PetscDualSpaceTransformType trans; 2110 PetscInt k; 2111 PetscErrorCode ierr; 2112 2113 PetscFunctionBeginHot; 2114 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2115 PetscValidPointer(fegeom, 2); 2116 PetscValidPointer(pointEval, 5); 2117 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2118 This determines their transformation properties. */ 2119 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2120 switch (k) 2121 { 2122 case 0: /* H^1 point evaluations */ 2123 trans = IDENTITY_TRANSFORM;break; 2124 case 1: /* Hcurl preserves tangential edge traces */ 2125 trans = COVARIANT_PIOLA_TRANSFORM;break; 2126 case 2: 2127 case 3: /* Hdiv preserve normal traces */ 2128 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2129 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2130 } 2131 ierr = PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 /*@C 2136 PetscDualSpacePushforwardGradient - Transform the given function gradient so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method. 2137 2138 Input Parameters: 2139 + dsp - The PetscDualSpace 2140 . fegeom - The geometry for this cell 2141 . Nq - The number of function gradient samples 2142 . Nc - The number of function components 2143 - pointEval - The function gradient values 2144 2145 Output Parameter: 2146 . pointEval - The transformed function gradient values 2147 2148 Level: advanced 2149 2150 Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex. 2151 2152 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2153 2154 .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2155 @*/ 2156 PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2157 { 2158 PetscDualSpaceTransformType trans; 2159 PetscInt k; 2160 PetscErrorCode ierr; 2161 2162 PetscFunctionBeginHot; 2163 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2164 PetscValidPointer(fegeom, 2); 2165 PetscValidPointer(pointEval, 5); 2166 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2167 This determines their transformation properties. */ 2168 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2169 switch (k) 2170 { 2171 case 0: /* H^1 point evaluations */ 2172 trans = IDENTITY_TRANSFORM;break; 2173 case 1: /* Hcurl preserves tangential edge traces */ 2174 trans = COVARIANT_PIOLA_TRANSFORM;break; 2175 case 2: 2176 case 3: /* Hdiv preserve normal traces */ 2177 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2178 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2179 } 2180 ierr = PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2181 PetscFunctionReturn(0); 2182 } 2183 2184 /*@C 2185 PetscDualSpacePushforwardHessian - Transform the given function Hessian so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method. 2186 2187 Input Parameters: 2188 + dsp - The PetscDualSpace 2189 . fegeom - The geometry for this cell 2190 . Nq - The number of function Hessian samples 2191 . Nc - The number of function components 2192 - pointEval - The function gradient values 2193 2194 Output Parameter: 2195 . pointEval - The transformed function Hessian values 2196 2197 Level: advanced 2198 2199 Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex. 2200 2201 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2202 2203 .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2204 @*/ 2205 PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2206 { 2207 PetscDualSpaceTransformType trans; 2208 PetscInt k; 2209 PetscErrorCode ierr; 2210 2211 PetscFunctionBeginHot; 2212 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2213 PetscValidPointer(fegeom, 2); 2214 PetscValidPointer(pointEval, 5); 2215 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2216 This determines their transformation properties. */ 2217 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2218 switch (k) 2219 { 2220 case 0: /* H^1 point evaluations */ 2221 trans = IDENTITY_TRANSFORM;break; 2222 case 1: /* Hcurl preserves tangential edge traces */ 2223 trans = COVARIANT_PIOLA_TRANSFORM;break; 2224 case 2: 2225 case 3: /* Hdiv preserve normal traces */ 2226 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2227 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2228 } 2229 ierr = PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2230 PetscFunctionReturn(0); 2231 } 2232