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 /*@ 1487 PetscDualSpaceEqual - Determine if a dual space is equivalent 1488 1489 Input Parameters: 1490 + A - A PetscDualSpace object 1491 - B - Another PetscDualSpace object 1492 1493 Output Parameter: 1494 . equal - PETSC_TRUE if the dual spaces are equivalent 1495 1496 Level: advanced 1497 1498 .seealso: PetscDualSpaceCreate() 1499 @*/ 1500 PetscErrorCode PetscDualSpaceEqual(PetscDualSpace A, PetscDualSpace B, PetscBool *equal) 1501 { 1502 PetscErrorCode ierr; 1503 PetscInt sizeA, sizeB, dimA, dimB; 1504 const PetscInt *dofA, *dofB; 1505 PetscQuadrature quadA, quadB; 1506 Mat matA, matB; 1507 1508 PetscFunctionBegin; 1509 PetscValidHeaderSpecific(A,PETSCDUALSPACE_CLASSID,1); 1510 PetscValidHeaderSpecific(B,PETSCDUALSPACE_CLASSID,2); 1511 PetscValidBoolPointer(equal,3); 1512 *equal = PETSC_FALSE; 1513 ierr = PetscDualSpaceGetDimension(A, &sizeA);CHKERRQ(ierr); 1514 ierr = PetscDualSpaceGetDimension(B, &sizeB);CHKERRQ(ierr); 1515 if (sizeB != sizeA) { 1516 PetscFunctionReturn(0); 1517 } 1518 ierr = DMGetDimension(A->dm, &dimA);CHKERRQ(ierr); 1519 ierr = DMGetDimension(B->dm, &dimB);CHKERRQ(ierr); 1520 if (dimA != dimB) { 1521 PetscFunctionReturn(0); 1522 } 1523 1524 ierr = PetscDualSpaceGetNumDof(A, &dofA);CHKERRQ(ierr); 1525 ierr = PetscDualSpaceGetNumDof(B, &dofB);CHKERRQ(ierr); 1526 for (PetscInt d=0; d<dimA; d++) { 1527 if (dofA[d] != dofB[d]) { 1528 PetscFunctionReturn(0); 1529 } 1530 } 1531 1532 ierr = PetscDualSpaceGetInteriorData(A, &quadA, &matA);CHKERRQ(ierr); 1533 ierr = PetscDualSpaceGetInteriorData(B, &quadB, &matB);CHKERRQ(ierr); 1534 if (!quadA && !quadB) { 1535 *equal = PETSC_TRUE; 1536 } else if (quadA && quadB) { 1537 ierr = PetscQuadratureEqual(quadA, quadB, equal);CHKERRQ(ierr); 1538 if (*equal == PETSC_FALSE) PetscFunctionReturn(0); 1539 if (!matA && !matB) PetscFunctionReturn(0); 1540 if (matA && matB) {ierr = MatEqual(matA, matB, equal);CHKERRQ(ierr);} 1541 else *equal = PETSC_FALSE; 1542 } 1543 PetscFunctionReturn(0); 1544 } 1545 1546 /*@C 1547 PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid. 1548 1549 Input Parameters: 1550 + sp - The PetscDualSpace object 1551 . f - The basis functional index 1552 . time - The time 1553 . cgeom - A context with geometric information for this cell, we currently just use the centroid 1554 . Nc - The number of components for the function 1555 . func - The input function 1556 - ctx - A context for the function 1557 1558 Output Parameter: 1559 . value - The output value (scalar) 1560 1561 Note: The calling sequence for the callback func is given by: 1562 1563 $ func(PetscInt dim, PetscReal time, const PetscReal x[], 1564 $ PetscInt numComponents, PetscScalar values[], void *ctx) 1565 1566 and the idea is to evaluate the functional as an integral 1567 1568 $ n(f) = int dx n(x) . f(x) 1569 1570 where both n and f have Nc components. 1571 1572 Level: beginner 1573 1574 .seealso: PetscDualSpaceCreate() 1575 @*/ 1576 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) 1577 { 1578 DM dm; 1579 PetscQuadrature n; 1580 const PetscReal *points, *weights; 1581 PetscScalar *val; 1582 PetscInt dimEmbed, qNc, c, Nq, q; 1583 PetscErrorCode ierr; 1584 1585 PetscFunctionBegin; 1586 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1587 PetscValidPointer(value, 8); 1588 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1589 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1590 ierr = PetscDualSpaceGetFunctional(sp, f, &n);CHKERRQ(ierr); 1591 ierr = PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);CHKERRQ(ierr); 1592 PetscCheckFalse(qNc != Nc,PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc); 1593 ierr = DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);CHKERRQ(ierr); 1594 *value = 0.; 1595 for (q = 0; q < Nq; ++q) { 1596 ierr = (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);CHKERRQ(ierr); 1597 for (c = 0; c < Nc; ++c) { 1598 *value += val[c]*weights[q*Nc+c]; 1599 } 1600 } 1601 ierr = DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);CHKERRQ(ierr); 1602 PetscFunctionReturn(0); 1603 } 1604 1605 /*@ 1606 PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a 1607 given height. This assumes that the reference cell is symmetric over points of this height. 1608 1609 If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and 1610 pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not 1611 support extracting subspaces, then NULL is returned. 1612 1613 This does not increment the reference count on the returned dual space, and the user should not destroy it. 1614 1615 Not collective 1616 1617 Input Parameters: 1618 + sp - the PetscDualSpace object 1619 - height - the height of the mesh point for which the subspace is desired 1620 1621 Output Parameter: 1622 . subsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the 1623 point, which will be of lesser dimension if height > 0. 1624 1625 Level: advanced 1626 1627 .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace 1628 @*/ 1629 PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp) 1630 { 1631 PetscInt depth = -1, cStart, cEnd; 1632 DM dm; 1633 PetscErrorCode ierr; 1634 1635 PetscFunctionBegin; 1636 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1637 PetscValidPointer(subsp,3); 1638 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"); 1639 *subsp = NULL; 1640 dm = sp->dm; 1641 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1642 PetscCheckFalse(height < 0 || height > depth,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid height"); 1643 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1644 if (height == 0 && cEnd == cStart + 1) { 1645 *subsp = sp; 1646 PetscFunctionReturn(0); 1647 } 1648 if (!sp->heightSpaces) { 1649 PetscInt h; 1650 ierr = PetscCalloc1(depth+1, &(sp->heightSpaces));CHKERRQ(ierr); 1651 1652 for (h = 0; h <= depth; h++) { 1653 if (h == 0 && cEnd == cStart + 1) continue; 1654 if (sp->ops->createheightsubspace) {ierr = (*sp->ops->createheightsubspace)(sp,height,&(sp->heightSpaces[h]));CHKERRQ(ierr);} 1655 else if (sp->pointSpaces) { 1656 PetscInt hStart, hEnd; 1657 1658 ierr = DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 1659 if (hEnd > hStart) { 1660 const char *name; 1661 1662 ierr = PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));CHKERRQ(ierr); 1663 if (sp->pointSpaces[hStart]) { 1664 ierr = PetscObjectGetName((PetscObject) sp, &name);CHKERRQ(ierr); 1665 ierr = PetscObjectSetName((PetscObject) sp->pointSpaces[hStart], name);CHKERRQ(ierr); 1666 } 1667 sp->heightSpaces[h] = sp->pointSpaces[hStart]; 1668 } 1669 } 1670 } 1671 } 1672 *subsp = sp->heightSpaces[height]; 1673 PetscFunctionReturn(0); 1674 } 1675 1676 /*@ 1677 PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point. 1678 1679 If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not 1680 defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting 1681 subspaces, then NULL is returned. 1682 1683 This does not increment the reference count on the returned dual space, and the user should not destroy it. 1684 1685 Not collective 1686 1687 Input Parameters: 1688 + sp - the PetscDualSpace object 1689 - point - the point (in the dual space's DM) for which the subspace is desired 1690 1691 Output Parameters: 1692 bdsp - the subspace. Note that the functionals in the subspace are with respect to the intrinsic geometry of the 1693 point, which will be of lesser dimension if height > 0. 1694 1695 Level: advanced 1696 1697 .seealso: PetscDualSpace 1698 @*/ 1699 PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp) 1700 { 1701 PetscInt pStart = 0, pEnd = 0, cStart, cEnd; 1702 DM dm; 1703 PetscErrorCode ierr; 1704 1705 PetscFunctionBegin; 1706 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1707 PetscValidPointer(bdsp,3); 1708 *bdsp = NULL; 1709 dm = sp->dm; 1710 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1711 PetscCheckFalse(point < pStart || point > pEnd,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point"); 1712 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1713 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 */ 1714 *bdsp = sp; 1715 PetscFunctionReturn(0); 1716 } 1717 if (!sp->pointSpaces) { 1718 PetscInt p; 1719 ierr = PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));CHKERRQ(ierr); 1720 1721 for (p = 0; p < pEnd - pStart; p++) { 1722 if (p + pStart == cStart && cEnd == cStart + 1) continue; 1723 if (sp->ops->createpointsubspace) {ierr = (*sp->ops->createpointsubspace)(sp,p+pStart,&(sp->pointSpaces[p]));CHKERRQ(ierr);} 1724 else if (sp->heightSpaces || sp->ops->createheightsubspace) { 1725 PetscInt dim, depth, height; 1726 DMLabel label; 1727 1728 ierr = DMPlexGetDepth(dm,&dim);CHKERRQ(ierr); 1729 ierr = DMPlexGetDepthLabel(dm,&label);CHKERRQ(ierr); 1730 ierr = DMLabelGetValue(label,p+pStart,&depth);CHKERRQ(ierr); 1731 height = dim - depth; 1732 ierr = PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));CHKERRQ(ierr); 1733 ierr = PetscObjectReference((PetscObject)sp->pointSpaces[p]);CHKERRQ(ierr); 1734 } 1735 } 1736 } 1737 *bdsp = sp->pointSpaces[point - pStart]; 1738 PetscFunctionReturn(0); 1739 } 1740 1741 /*@C 1742 PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis 1743 1744 Not collective 1745 1746 Input Parameter: 1747 . sp - the PetscDualSpace object 1748 1749 Output Parameters: 1750 + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation 1751 - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation 1752 1753 Note: The permutation and flip arrays are organized in the following way 1754 $ perms[p][ornt][dof # on point] = new local dof # 1755 $ flips[p][ornt][dof # on point] = reversal or not 1756 1757 Level: developer 1758 1759 @*/ 1760 PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips) 1761 { 1762 PetscErrorCode ierr; 1763 1764 PetscFunctionBegin; 1765 PetscValidHeaderSpecific(sp,PETSCDUALSPACE_CLASSID,1); 1766 if (perms) {PetscValidPointer(perms,2); *perms = NULL;} 1767 if (flips) {PetscValidPointer(flips,3); *flips = NULL;} 1768 if (sp->ops->getsymmetries) {ierr = (sp->ops->getsymmetries)(sp,perms,flips);CHKERRQ(ierr);} 1769 PetscFunctionReturn(0); 1770 } 1771 1772 /*@ 1773 PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this 1774 dual space's functionals. 1775 1776 Input Parameter: 1777 . dsp - The PetscDualSpace 1778 1779 Output Parameter: 1780 . k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored 1781 in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example, 1782 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). 1783 If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the 1784 Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms 1785 but are stored as 1-forms. 1786 1787 Level: developer 1788 1789 .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType 1790 @*/ 1791 PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k) 1792 { 1793 PetscFunctionBeginHot; 1794 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1795 PetscValidPointer(k, 2); 1796 *k = dsp->k; 1797 PetscFunctionReturn(0); 1798 } 1799 1800 /*@ 1801 PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this 1802 dual space's functionals. 1803 1804 Input Parameters: 1805 + dsp - The PetscDualSpace 1806 - k - The *signed* degree k of the k. If k >= 0, this means that the degrees of freedom are k-forms, and are stored 1807 in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms. So for example, 1808 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). 1809 If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the 1810 Hodge star map. So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms 1811 but are stored as 1-forms. 1812 1813 Level: developer 1814 1815 .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType 1816 @*/ 1817 PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k) 1818 { 1819 PetscInt dim; 1820 1821 PetscFunctionBeginHot; 1822 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1823 PetscCheckFalse(dsp->setupcalled,PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up"); 1824 dim = dsp->dm->dim; 1825 PetscCheckFalse(k < -dim || k > dim,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported %D-form on %D-dimensional reference cell", PetscAbsInt(k), dim); 1826 dsp->k = k; 1827 PetscFunctionReturn(0); 1828 } 1829 1830 /*@ 1831 PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space 1832 1833 Input Parameter: 1834 . dsp - The PetscDualSpace 1835 1836 Output Parameter: 1837 . k - The simplex dimension 1838 1839 Level: developer 1840 1841 Note: Currently supported values are 1842 $ 0: These are H_1 methods that only transform coordinates 1843 $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM) 1844 $ 2: These are the same as 1 1845 $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM) 1846 1847 .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType 1848 @*/ 1849 PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k) 1850 { 1851 PetscInt dim; 1852 1853 PetscFunctionBeginHot; 1854 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1855 PetscValidPointer(k, 2); 1856 dim = dsp->dm->dim; 1857 if (!dsp->k) *k = IDENTITY_TRANSFORM; 1858 else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM; 1859 else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM; 1860 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation"); 1861 PetscFunctionReturn(0); 1862 } 1863 1864 /*@C 1865 PetscDualSpaceTransform - Transform the function values 1866 1867 Input Parameters: 1868 + dsp - The PetscDualSpace 1869 . trans - The type of transform 1870 . isInverse - Flag to invert the transform 1871 . fegeom - The cell geometry 1872 . Nv - The number of function samples 1873 . Nc - The number of function components 1874 - vals - The function values 1875 1876 Output Parameter: 1877 . vals - The transformed function values 1878 1879 Level: intermediate 1880 1881 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1882 1883 .seealso: PetscDualSpaceTransformGradient(), PetscDualSpaceTransformHessian(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType 1884 @*/ 1885 PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[]) 1886 { 1887 PetscReal Jstar[9] = {0}; 1888 PetscInt dim, v, c, Nk; 1889 PetscErrorCode ierr; 1890 1891 PetscFunctionBeginHot; 1892 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1893 PetscValidPointer(fegeom, 4); 1894 PetscValidPointer(vals, 7); 1895 /* TODO: not handling dimEmbed != dim right now */ 1896 dim = dsp->dm->dim; 1897 /* No change needed for 0-forms */ 1898 if (!dsp->k) PetscFunctionReturn(0); 1899 ierr = PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);CHKERRQ(ierr); 1900 /* TODO: use fegeom->isAffine */ 1901 ierr = PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);CHKERRQ(ierr); 1902 for (v = 0; v < Nv; ++v) { 1903 switch (Nk) { 1904 case 1: 1905 for (c = 0; c < Nc; c++) vals[v*Nc + c] *= Jstar[0]; 1906 break; 1907 case 2: 1908 for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]); 1909 break; 1910 case 3: 1911 for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]); 1912 break; 1913 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %D for transformation", Nk); 1914 } 1915 } 1916 PetscFunctionReturn(0); 1917 } 1918 1919 /*@C 1920 PetscDualSpaceTransformGradient - Transform the function gradient values 1921 1922 Input Parameters: 1923 + dsp - The PetscDualSpace 1924 . trans - The type of transform 1925 . isInverse - Flag to invert the transform 1926 . fegeom - The cell geometry 1927 . Nv - The number of function gradient samples 1928 . Nc - The number of function components 1929 - vals - The function gradient values 1930 1931 Output Parameter: 1932 . vals - The transformed function gradient values 1933 1934 Level: intermediate 1935 1936 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1937 1938 .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType 1939 @*/ 1940 PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[]) 1941 { 1942 const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed; 1943 PetscInt v, c, d; 1944 1945 PetscFunctionBeginHot; 1946 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 1947 PetscValidPointer(fegeom, 4); 1948 PetscValidPointer(vals, 7); 1949 #ifdef PETSC_USE_DEBUG 1950 PetscCheckFalse(dE <= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE); 1951 #endif 1952 /* Transform gradient */ 1953 if (dim == dE) { 1954 for (v = 0; v < Nv; ++v) { 1955 for (c = 0; c < Nc; ++c) { 1956 switch (dim) 1957 { 1958 case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];break; 1959 case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break; 1960 case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break; 1961 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1962 } 1963 } 1964 } 1965 } else { 1966 for (v = 0; v < Nv; ++v) { 1967 for (c = 0; c < Nc; ++c) { 1968 DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v*Nc+c)*dE], &vals[(v*Nc+c)*dE]); 1969 } 1970 } 1971 } 1972 /* Assume its a vector, otherwise assume its a bunch of scalars */ 1973 if (Nc == 1 || Nc != dim) PetscFunctionReturn(0); 1974 switch (trans) { 1975 case IDENTITY_TRANSFORM: break; 1976 case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */ 1977 if (isInverse) { 1978 for (v = 0; v < Nv; ++v) { 1979 for (d = 0; d < dim; ++d) { 1980 switch (dim) 1981 { 1982 case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1983 case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1984 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1985 } 1986 } 1987 } 1988 } else { 1989 for (v = 0; v < Nv; ++v) { 1990 for (d = 0; d < dim; ++d) { 1991 switch (dim) 1992 { 1993 case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1994 case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 1995 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 1996 } 1997 } 1998 } 1999 } 2000 break; 2001 case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */ 2002 if (isInverse) { 2003 for (v = 0; v < Nv; ++v) { 2004 for (d = 0; d < dim; ++d) { 2005 switch (dim) 2006 { 2007 case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 2008 case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 2009 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 2010 } 2011 for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0]; 2012 } 2013 } 2014 } else { 2015 for (v = 0; v < Nv; ++v) { 2016 for (d = 0; d < dim; ++d) { 2017 switch (dim) 2018 { 2019 case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 2020 case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break; 2021 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 2022 } 2023 for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0]; 2024 } 2025 } 2026 } 2027 break; 2028 } 2029 PetscFunctionReturn(0); 2030 } 2031 2032 /*@C 2033 PetscDualSpaceTransformHessian - Transform the function Hessian values 2034 2035 Input Parameters: 2036 + dsp - The PetscDualSpace 2037 . trans - The type of transform 2038 . isInverse - Flag to invert the transform 2039 . fegeom - The cell geometry 2040 . Nv - The number of function Hessian samples 2041 . Nc - The number of function components 2042 - vals - The function gradient values 2043 2044 Output Parameter: 2045 . vals - The transformed function Hessian values 2046 2047 Level: intermediate 2048 2049 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2050 2051 .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType 2052 @*/ 2053 PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[]) 2054 { 2055 const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed; 2056 PetscInt v, c; 2057 2058 PetscFunctionBeginHot; 2059 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2060 PetscValidPointer(fegeom, 4); 2061 PetscValidPointer(vals, 7); 2062 #ifdef PETSC_USE_DEBUG 2063 PetscCheckFalse(dE <= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE); 2064 #endif 2065 /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */ 2066 if (dim == dE) { 2067 for (v = 0; v < Nv; ++v) { 2068 for (c = 0; c < Nc; ++c) { 2069 switch (dim) 2070 { 2071 case 1: vals[(v*Nc+c)*dim*dim] *= PetscSqr(fegeom->invJ[0]);break; 2072 case 2: DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break; 2073 case 3: DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break; 2074 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim); 2075 } 2076 } 2077 } 2078 } else { 2079 for (v = 0; v < Nv; ++v) { 2080 for (c = 0; c < Nc; ++c) { 2081 DMPlex_PTAPReal_Internal(fegeom->invJ, dim, dE, &vals[(v*Nc+c)*dE*dE], &vals[(v*Nc+c)*dE*dE]); 2082 } 2083 } 2084 } 2085 /* Assume its a vector, otherwise assume its a bunch of scalars */ 2086 if (Nc == 1 || Nc != dim) PetscFunctionReturn(0); 2087 switch (trans) { 2088 case IDENTITY_TRANSFORM: break; 2089 case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */ 2090 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported"); 2091 case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */ 2092 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported"); 2093 } 2094 PetscFunctionReturn(0); 2095 } 2096 2097 /*@C 2098 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. 2099 2100 Input Parameters: 2101 + dsp - The PetscDualSpace 2102 . fegeom - The geometry for this cell 2103 . Nq - The number of function samples 2104 . Nc - The number of function components 2105 - pointEval - The function values 2106 2107 Output Parameter: 2108 . pointEval - The transformed function values 2109 2110 Level: advanced 2111 2112 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. 2113 2114 Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2115 2116 .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2117 @*/ 2118 PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2119 { 2120 PetscDualSpaceTransformType trans; 2121 PetscInt k; 2122 PetscErrorCode ierr; 2123 2124 PetscFunctionBeginHot; 2125 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2126 PetscValidPointer(fegeom, 2); 2127 PetscValidPointer(pointEval, 5); 2128 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2129 This determines their transformation properties. */ 2130 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2131 switch (k) 2132 { 2133 case 0: /* H^1 point evaluations */ 2134 trans = IDENTITY_TRANSFORM;break; 2135 case 1: /* Hcurl preserves tangential edge traces */ 2136 trans = COVARIANT_PIOLA_TRANSFORM;break; 2137 case 2: 2138 case 3: /* Hdiv preserve normal traces */ 2139 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2140 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2141 } 2142 ierr = PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2143 PetscFunctionReturn(0); 2144 } 2145 2146 /*@C 2147 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. 2148 2149 Input Parameters: 2150 + dsp - The PetscDualSpace 2151 . fegeom - The geometry for this cell 2152 . Nq - The number of function samples 2153 . Nc - The number of function components 2154 - pointEval - The function values 2155 2156 Output Parameter: 2157 . pointEval - The transformed function values 2158 2159 Level: advanced 2160 2161 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. 2162 2163 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2164 2165 .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2166 @*/ 2167 PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2168 { 2169 PetscDualSpaceTransformType trans; 2170 PetscInt k; 2171 PetscErrorCode ierr; 2172 2173 PetscFunctionBeginHot; 2174 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2175 PetscValidPointer(fegeom, 2); 2176 PetscValidPointer(pointEval, 5); 2177 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2178 This determines their transformation properties. */ 2179 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2180 switch (k) 2181 { 2182 case 0: /* H^1 point evaluations */ 2183 trans = IDENTITY_TRANSFORM;break; 2184 case 1: /* Hcurl preserves tangential edge traces */ 2185 trans = COVARIANT_PIOLA_TRANSFORM;break; 2186 case 2: 2187 case 3: /* Hdiv preserve normal traces */ 2188 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2189 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2190 } 2191 ierr = PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2192 PetscFunctionReturn(0); 2193 } 2194 2195 /*@C 2196 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. 2197 2198 Input Parameters: 2199 + dsp - The PetscDualSpace 2200 . fegeom - The geometry for this cell 2201 . Nq - The number of function gradient samples 2202 . Nc - The number of function components 2203 - pointEval - The function gradient values 2204 2205 Output Parameter: 2206 . pointEval - The transformed function gradient values 2207 2208 Level: advanced 2209 2210 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. 2211 2212 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2213 2214 .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2215 @*/ 2216 PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2217 { 2218 PetscDualSpaceTransformType trans; 2219 PetscInt k; 2220 PetscErrorCode ierr; 2221 2222 PetscFunctionBeginHot; 2223 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2224 PetscValidPointer(fegeom, 2); 2225 PetscValidPointer(pointEval, 5); 2226 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2227 This determines their transformation properties. */ 2228 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2229 switch (k) 2230 { 2231 case 0: /* H^1 point evaluations */ 2232 trans = IDENTITY_TRANSFORM;break; 2233 case 1: /* Hcurl preserves tangential edge traces */ 2234 trans = COVARIANT_PIOLA_TRANSFORM;break; 2235 case 2: 2236 case 3: /* Hdiv preserve normal traces */ 2237 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2238 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2239 } 2240 ierr = PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2241 PetscFunctionReturn(0); 2242 } 2243 2244 /*@C 2245 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. 2246 2247 Input Parameters: 2248 + dsp - The PetscDualSpace 2249 . fegeom - The geometry for this cell 2250 . Nq - The number of function Hessian samples 2251 . Nc - The number of function components 2252 - pointEval - The function gradient values 2253 2254 Output Parameter: 2255 . pointEval - The transformed function Hessian values 2256 2257 Level: advanced 2258 2259 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. 2260 2261 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 2262 2263 .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm() 2264 @*/ 2265 PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[]) 2266 { 2267 PetscDualSpaceTransformType trans; 2268 PetscInt k; 2269 PetscErrorCode ierr; 2270 2271 PetscFunctionBeginHot; 2272 PetscValidHeaderSpecific(dsp, PETSCDUALSPACE_CLASSID, 1); 2273 PetscValidPointer(fegeom, 2); 2274 PetscValidPointer(pointEval, 5); 2275 /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k. 2276 This determines their transformation properties. */ 2277 ierr = PetscDualSpaceGetDeRahm(dsp, &k);CHKERRQ(ierr); 2278 switch (k) 2279 { 2280 case 0: /* H^1 point evaluations */ 2281 trans = IDENTITY_TRANSFORM;break; 2282 case 1: /* Hcurl preserves tangential edge traces */ 2283 trans = COVARIANT_PIOLA_TRANSFORM;break; 2284 case 2: 2285 case 3: /* Hdiv preserve normal traces */ 2286 trans = CONTRAVARIANT_PIOLA_TRANSFORM;break; 2287 default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k); 2288 } 2289 ierr = PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);CHKERRQ(ierr); 2290 PetscFunctionReturn(0); 2291 } 2292