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