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