1 #include <petsc/private/petscdsimpl.h> /*I "petscds.h" I*/ 2 3 PetscClassId PETSCDS_CLASSID = 0; 4 5 PetscFunctionList PetscDSList = NULL; 6 PetscBool PetscDSRegisterAllCalled = PETSC_FALSE; 7 8 /* A PetscDS (Discrete System) encodes a set of equations posed in a discrete space, which represents a set of 9 nonlinear continuum equations. The equations can have multiple fields, each field having a different 10 discretization. In addition, different pieces of the domain can have different field combinations and equations. 11 12 The DS provides the user a description of the approximation space on any given cell. It also gives pointwise 13 functions representing the equations. 14 15 Each field is associated with a label, marking the cells on which it is supported. Note that a field can be 16 supported on the closure of a cell not in the label due to overlap of the boundary of neighboring cells. The DM 17 then creates a DS for each set of cells with identical approximation spaces. When assembling, the user asks for 18 the space associated with a given cell. DMPlex uses the labels associated with each DS in the default integration loop. 19 */ 20 21 /*@C 22 PetscDSRegister - Adds a new `PetscDS` implementation 23 24 Not Collective; No Fortran Support 25 26 Input Parameters: 27 + sname - The name of a new user-defined creation routine 28 - function - The creation routine itself 29 30 Example Usage: 31 .vb 32 PetscDSRegister("my_ds", MyPetscDSCreate); 33 .ve 34 35 Then, your PetscDS type can be chosen with the procedural interface via 36 .vb 37 PetscDSCreate(MPI_Comm, PetscDS *); 38 PetscDSSetType(PetscDS, "my_ds"); 39 .ve 40 or at runtime via the option 41 .vb 42 -petscds_type my_ds 43 .ve 44 45 Level: advanced 46 47 Note: 48 `PetscDSRegister()` may be called multiple times to add several user-defined `PetscDSs` 49 50 .seealso: `PetscDSType`, `PetscDS`, `PetscDSRegisterAll()`, `PetscDSRegisterDestroy()` 51 @*/ 52 PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS)) 53 { 54 PetscFunctionBegin; 55 PetscCall(PetscFunctionListAdd(&PetscDSList, sname, function)); 56 PetscFunctionReturn(PETSC_SUCCESS); 57 } 58 59 /*@ 60 PetscDSSetType - Builds a particular `PetscDS` 61 62 Collective; No Fortran Support 63 64 Input Parameters: 65 + prob - The `PetscDS` object 66 - name - The `PetscDSType` 67 68 Options Database Key: 69 . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types 70 71 Level: intermediate 72 73 .seealso: `PetscDSType`, `PetscDS`, `PetscDSGetType()`, `PetscDSCreate()` 74 @*/ 75 PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name) 76 { 77 PetscErrorCode (*r)(PetscDS); 78 PetscBool match; 79 80 PetscFunctionBegin; 81 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 82 PetscCall(PetscObjectTypeCompare((PetscObject)prob, name, &match)); 83 if (match) PetscFunctionReturn(PETSC_SUCCESS); 84 85 PetscCall(PetscDSRegisterAll()); 86 PetscCall(PetscFunctionListFind(PetscDSList, name, &r)); 87 PetscCheck(r, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name); 88 89 PetscTryTypeMethod(prob, destroy); 90 prob->ops->destroy = NULL; 91 92 PetscCall((*r)(prob)); 93 PetscCall(PetscObjectChangeTypeName((PetscObject)prob, name)); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 /*@ 98 PetscDSGetType - Gets the `PetscDSType` name (as a string) from the `PetscDS` 99 100 Not Collective; No Fortran Support 101 102 Input Parameter: 103 . prob - The `PetscDS` 104 105 Output Parameter: 106 . name - The `PetscDSType` name 107 108 Level: intermediate 109 110 .seealso: `PetscDSType`, `PetscDS`, `PetscDSSetType()`, `PetscDSCreate()` 111 @*/ 112 PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name) 113 { 114 PetscFunctionBegin; 115 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 116 PetscAssertPointer(name, 2); 117 PetscCall(PetscDSRegisterAll()); 118 *name = ((PetscObject)prob)->type_name; 119 PetscFunctionReturn(PETSC_SUCCESS); 120 } 121 122 static PetscErrorCode PetscDSView_Ascii(PetscDS ds, PetscViewer viewer) 123 { 124 PetscViewerFormat format; 125 const PetscScalar *constants; 126 PetscInt Nf, numConstants, f; 127 128 PetscFunctionBegin; 129 PetscCall(PetscDSGetNumFields(ds, &Nf)); 130 PetscCall(PetscViewerGetFormat(viewer, &format)); 131 PetscCall(PetscViewerASCIIPrintf(viewer, "Discrete System with %" PetscInt_FMT " fields\n", Nf)); 132 PetscCall(PetscViewerASCIIPushTab(viewer)); 133 PetscCall(PetscViewerASCIIPrintf(viewer, " cell total dim %" PetscInt_FMT " total comp %" PetscInt_FMT "\n", ds->totDim, ds->totComp)); 134 if (ds->isCohesive) PetscCall(PetscViewerASCIIPrintf(viewer, " cohesive cell\n")); 135 for (f = 0; f < Nf; ++f) { 136 DSBoundary b; 137 PetscObject obj; 138 PetscClassId id; 139 PetscQuadrature q; 140 const char *name; 141 PetscInt Nc, Nq, Nqc; 142 143 PetscCall(PetscDSGetDiscretization(ds, f, &obj)); 144 PetscCall(PetscObjectGetClassId(obj, &id)); 145 PetscCall(PetscObjectGetName(obj, &name)); 146 PetscCall(PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>")); 147 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 148 if (id == PETSCFE_CLASSID) { 149 PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc)); 150 PetscCall(PetscFEGetQuadrature((PetscFE)obj, &q)); 151 PetscCall(PetscViewerASCIIPrintf(viewer, " FEM")); 152 } else if (id == PETSCFV_CLASSID) { 153 PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc)); 154 PetscCall(PetscFVGetQuadrature((PetscFV)obj, &q)); 155 PetscCall(PetscViewerASCIIPrintf(viewer, " FVM")); 156 } else SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f); 157 if (Nc > 1) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " components", Nc)); 158 else PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " component ", Nc)); 159 if (ds->implicit[f]) PetscCall(PetscViewerASCIIPrintf(viewer, " (implicit)")); 160 else PetscCall(PetscViewerASCIIPrintf(viewer, " (explicit)")); 161 if (q) { 162 PetscCall(PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL)); 163 PetscCall(PetscViewerASCIIPrintf(viewer, " (Nq %" PetscInt_FMT " Nqc %" PetscInt_FMT ")", Nq, Nqc)); 164 } 165 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "-jet", ds->jetDegree[f])); 166 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 167 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 168 PetscCall(PetscViewerASCIIPushTab(viewer)); 169 if (id == PETSCFE_CLASSID) PetscCall(PetscFEView((PetscFE)obj, viewer)); 170 else if (id == PETSCFV_CLASSID) PetscCall(PetscFVView((PetscFV)obj, viewer)); 171 PetscCall(PetscViewerASCIIPopTab(viewer)); 172 173 for (b = ds->boundary; b; b = b->next) { 174 char *name; 175 PetscInt c, i; 176 177 if (b->field != f) continue; 178 PetscCall(PetscViewerASCIIPushTab(viewer)); 179 PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->lname, DMBoundaryConditionTypes[b->type])); 180 if (!b->Nc) { 181 PetscCall(PetscViewerASCIIPrintf(viewer, " all components\n")); 182 } else { 183 PetscCall(PetscViewerASCIIPrintf(viewer, " components: ")); 184 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 185 for (c = 0; c < b->Nc; ++c) { 186 if (c > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", ")); 187 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->comps[c])); 188 } 189 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 190 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 191 } 192 PetscCall(PetscViewerASCIIPrintf(viewer, " values: ")); 193 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 194 for (i = 0; i < b->Nv; ++i) { 195 if (i > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", ")); 196 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->values[i])); 197 } 198 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 199 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 200 #if defined(__clang__) 201 PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic") 202 #elif defined(__GNUC__) || defined(__GNUG__) 203 PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat") 204 #endif 205 if (b->func) { 206 PetscCall(PetscDLAddr(b->func, &name)); 207 if (name) PetscCall(PetscViewerASCIIPrintf(viewer, " func: %s\n", name)); 208 else PetscCall(PetscViewerASCIIPrintf(viewer, " func: %p\n", b->func)); 209 PetscCall(PetscFree(name)); 210 } 211 if (b->func_t) { 212 PetscCall(PetscDLAddr(b->func_t, &name)); 213 if (name) PetscCall(PetscViewerASCIIPrintf(viewer, " func_t: %s\n", name)); 214 else PetscCall(PetscViewerASCIIPrintf(viewer, " func_t: %p\n", b->func_t)); 215 PetscCall(PetscFree(name)); 216 } 217 PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END() 218 PetscCall(PetscWeakFormView(b->wf, viewer)); 219 PetscCall(PetscViewerASCIIPopTab(viewer)); 220 } 221 } 222 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 223 if (numConstants) { 224 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " constants\n", numConstants)); 225 PetscCall(PetscViewerASCIIPushTab(viewer)); 226 for (f = 0; f < numConstants; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(constants[f]))); 227 PetscCall(PetscViewerASCIIPopTab(viewer)); 228 } 229 PetscCall(PetscWeakFormView(ds->wf, viewer)); 230 PetscCall(PetscViewerASCIIPopTab(viewer)); 231 PetscFunctionReturn(PETSC_SUCCESS); 232 } 233 234 /*@ 235 PetscDSViewFromOptions - View a `PetscDS` based on values in the options database 236 237 Collective 238 239 Input Parameters: 240 + A - the `PetscDS` object 241 . obj - Optional object that provides the options prefix used in the search 242 - name - command line option 243 244 Level: intermediate 245 246 .seealso: `PetscDSType`, `PetscDS`, `PetscDSView()`, `PetscObjectViewFromOptions()`, `PetscDSCreate()` 247 @*/ 248 PetscErrorCode PetscDSViewFromOptions(PetscDS A, PetscObject obj, const char name[]) 249 { 250 PetscFunctionBegin; 251 PetscValidHeaderSpecific(A, PETSCDS_CLASSID, 1); 252 PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name)); 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255 256 /*@ 257 PetscDSView - Views a `PetscDS` 258 259 Collective 260 261 Input Parameters: 262 + prob - the `PetscDS` object to view 263 - v - the viewer 264 265 Level: developer 266 267 .seealso: `PetscDSType`, `PetscDS`, `PetscViewer`, `PetscDSDestroy()`, `PetscDSViewFromOptions()` 268 @*/ 269 PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v) 270 { 271 PetscBool iascii; 272 273 PetscFunctionBegin; 274 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 275 if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)prob), &v)); 276 else PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2); 277 PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii)); 278 if (iascii) PetscCall(PetscDSView_Ascii(prob, v)); 279 PetscTryTypeMethod(prob, view, v); 280 PetscFunctionReturn(PETSC_SUCCESS); 281 } 282 283 /*@ 284 PetscDSSetFromOptions - sets parameters in a `PetscDS` from the options database 285 286 Collective 287 288 Input Parameter: 289 . prob - the `PetscDS` object to set options for 290 291 Options Database Keys: 292 + -petscds_type <type> - Set the `PetscDS` type 293 . -petscds_view <view opt> - View the `PetscDS` 294 . -petscds_jac_pre - Turn formation of a separate Jacobian preconditioner on or off 295 . -bc_<name> <ids> - Specify a list of label ids for a boundary condition 296 - -bc_<name>_comp <comps> - Specify a list of field components to constrain for a boundary condition 297 298 Level: intermediate 299 300 .seealso: `PetscDS`, `PetscDSView()` 301 @*/ 302 PetscErrorCode PetscDSSetFromOptions(PetscDS prob) 303 { 304 DSBoundary b; 305 const char *defaultType; 306 char name[256]; 307 PetscBool flg; 308 309 PetscFunctionBegin; 310 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 311 if (!((PetscObject)prob)->type_name) { 312 defaultType = PETSCDSBASIC; 313 } else { 314 defaultType = ((PetscObject)prob)->type_name; 315 } 316 PetscCall(PetscDSRegisterAll()); 317 318 PetscObjectOptionsBegin((PetscObject)prob); 319 for (b = prob->boundary; b; b = b->next) { 320 char optname[1024]; 321 PetscInt ids[1024], len = 1024; 322 PetscBool flg; 323 324 PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name)); 325 PetscCall(PetscMemzero(ids, sizeof(ids))); 326 PetscCall(PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg)); 327 if (flg) { 328 b->Nv = len; 329 PetscCall(PetscFree(b->values)); 330 PetscCall(PetscMalloc1(len, &b->values)); 331 PetscCall(PetscArraycpy(b->values, ids, len)); 332 PetscCall(PetscWeakFormRewriteKeys(b->wf, b->label, len, b->values)); 333 } 334 len = 1024; 335 PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name)); 336 PetscCall(PetscMemzero(ids, sizeof(ids))); 337 PetscCall(PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg)); 338 if (flg) { 339 b->Nc = len; 340 PetscCall(PetscFree(b->comps)); 341 PetscCall(PetscMalloc1(len, &b->comps)); 342 PetscCall(PetscArraycpy(b->comps, ids, len)); 343 } 344 } 345 PetscCall(PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg)); 346 if (flg) { 347 PetscCall(PetscDSSetType(prob, name)); 348 } else if (!((PetscObject)prob)->type_name) { 349 PetscCall(PetscDSSetType(prob, defaultType)); 350 } 351 PetscCall(PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg)); 352 PetscCall(PetscOptionsBool("-petscds_force_quad", "Discrete System", "PetscDSSetForceQuad", prob->forceQuad, &prob->forceQuad, &flg)); 353 PetscCall(PetscOptionsInt("-petscds_print_integrate", "Discrete System", "", prob->printIntegrate, &prob->printIntegrate, NULL)); 354 PetscTryTypeMethod(prob, setfromoptions); 355 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 356 PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)prob, PetscOptionsObject)); 357 PetscOptionsEnd(); 358 if (prob->Nf) PetscCall(PetscDSViewFromOptions(prob, NULL, "-petscds_view")); 359 PetscFunctionReturn(PETSC_SUCCESS); 360 } 361 362 /*@ 363 PetscDSSetUp - Construct data structures for the `PetscDS` 364 365 Collective 366 367 Input Parameter: 368 . prob - the `PetscDS` object to setup 369 370 Level: developer 371 372 .seealso: `PetscDS`, `PetscDSView()`, `PetscDSDestroy()` 373 @*/ 374 PetscErrorCode PetscDSSetUp(PetscDS prob) 375 { 376 const PetscInt Nf = prob->Nf; 377 PetscBool hasH = PETSC_FALSE; 378 PetscInt maxOrder[4] = {-2, -2, -2, -2}; 379 PetscInt dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f; 380 381 PetscFunctionBegin; 382 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 383 if (prob->setup) PetscFunctionReturn(PETSC_SUCCESS); 384 /* Calculate sizes */ 385 PetscCall(PetscDSGetSpatialDimension(prob, &dim)); 386 PetscCall(PetscDSGetCoordinateDimension(prob, &dimEmbed)); 387 prob->totDim = prob->totComp = 0; 388 PetscCall(PetscMalloc2(Nf, &prob->Nc, Nf, &prob->Nb)); 389 PetscCall(PetscCalloc2(Nf + 1, &prob->off, Nf + 1, &prob->offDer)); 390 PetscCall(PetscCalloc6(Nf + 1, &prob->offCohesive[0], Nf + 1, &prob->offCohesive[1], Nf + 1, &prob->offCohesive[2], Nf + 1, &prob->offDerCohesive[0], Nf + 1, &prob->offDerCohesive[1], Nf + 1, &prob->offDerCohesive[2])); 391 PetscCall(PetscMalloc2(Nf, &prob->T, Nf, &prob->Tf)); 392 if (prob->forceQuad) { 393 // Note: This assumes we have one kind of cell at each dimension. 394 // We can fix this by having quadrature hold the celltype 395 PetscQuadrature maxQuad[4] = {NULL, NULL, NULL, NULL}; 396 397 for (f = 0; f < Nf; ++f) { 398 PetscObject obj; 399 PetscClassId id; 400 PetscQuadrature q = NULL, fq = NULL; 401 PetscInt dim = -1, order = -1, forder = -1; 402 403 PetscCall(PetscDSGetDiscretization(prob, f, &obj)); 404 if (!obj) continue; 405 PetscCall(PetscObjectGetClassId(obj, &id)); 406 if (id == PETSCFE_CLASSID) { 407 PetscFE fe = (PetscFE)obj; 408 409 PetscCall(PetscFEGetQuadrature(fe, &q)); 410 PetscCall(PetscFEGetFaceQuadrature(fe, &fq)); 411 } else if (id == PETSCFV_CLASSID) { 412 PetscFV fv = (PetscFV)obj; 413 414 PetscCall(PetscFVGetQuadrature(fv, &q)); 415 } 416 if (q) { 417 PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL)); 418 PetscCall(PetscQuadratureGetOrder(q, &order)); 419 if (order > maxOrder[dim]) { 420 maxOrder[dim] = order; 421 maxQuad[dim] = q; 422 } 423 } 424 if (fq) { 425 PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL)); 426 PetscCall(PetscQuadratureGetOrder(fq, &forder)); 427 if (forder > maxOrder[dim]) { 428 maxOrder[dim] = forder; 429 maxQuad[dim] = fq; 430 } 431 } 432 } 433 for (f = 0; f < Nf; ++f) { 434 PetscObject obj; 435 PetscClassId id; 436 PetscQuadrature q; 437 PetscInt dim; 438 439 PetscCall(PetscDSGetDiscretization(prob, f, &obj)); 440 if (!obj) continue; 441 PetscCall(PetscObjectGetClassId(obj, &id)); 442 if (id == PETSCFE_CLASSID) { 443 PetscFE fe = (PetscFE)obj; 444 445 PetscCall(PetscFEGetQuadrature(fe, &q)); 446 PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL)); 447 PetscCall(PetscFESetQuadrature(fe, maxQuad[dim])); 448 PetscCall(PetscFESetFaceQuadrature(fe, dim ? maxQuad[dim - 1] : NULL)); 449 } else if (id == PETSCFV_CLASSID) { 450 PetscFV fv = (PetscFV)obj; 451 452 PetscCall(PetscFVGetQuadrature(fv, &q)); 453 PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL)); 454 PetscCall(PetscFVSetQuadrature(fv, maxQuad[dim])); 455 } 456 } 457 } 458 for (f = 0; f < Nf; ++f) { 459 PetscObject obj; 460 PetscClassId id; 461 PetscQuadrature q = NULL; 462 PetscInt Nq = 0, Nb, Nc; 463 464 PetscCall(PetscDSGetDiscretization(prob, f, &obj)); 465 if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE; 466 if (!obj) { 467 /* Empty mesh */ 468 Nb = Nc = 0; 469 prob->T[f] = prob->Tf[f] = NULL; 470 } else { 471 PetscCall(PetscObjectGetClassId(obj, &id)); 472 if (id == PETSCFE_CLASSID) { 473 PetscFE fe = (PetscFE)obj; 474 475 PetscCall(PetscFEGetQuadrature(fe, &q)); 476 { 477 PetscQuadrature fq; 478 PetscInt dim, order; 479 480 PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL)); 481 PetscCall(PetscQuadratureGetOrder(q, &order)); 482 if (maxOrder[dim] < 0) maxOrder[dim] = order; 483 PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " cell quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS cell quadrature order", f, order, maxOrder[dim]); 484 PetscCall(PetscFEGetFaceQuadrature(fe, &fq)); 485 if (fq) { 486 PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL)); 487 PetscCall(PetscQuadratureGetOrder(fq, &order)); 488 if (maxOrder[dim] < 0) maxOrder[dim] = order; 489 PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " face quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS face quadrature order", f, order, maxOrder[dim]); 490 } 491 } 492 PetscCall(PetscFEGetDimension(fe, &Nb)); 493 PetscCall(PetscFEGetNumComponents(fe, &Nc)); 494 PetscCall(PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f])); 495 PetscCall(PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f])); 496 } else if (id == PETSCFV_CLASSID) { 497 PetscFV fv = (PetscFV)obj; 498 499 PetscCall(PetscFVGetQuadrature(fv, &q)); 500 PetscCall(PetscFVGetNumComponents(fv, &Nc)); 501 Nb = Nc; 502 PetscCall(PetscFVGetCellTabulation(fv, &prob->T[f])); 503 /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */ 504 } else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f); 505 } 506 prob->Nc[f] = Nc; 507 prob->Nb[f] = Nb; 508 prob->off[f + 1] = Nc + prob->off[f]; 509 prob->offDer[f + 1] = Nc * dim + prob->offDer[f]; 510 prob->offCohesive[0][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[0][f]; 511 prob->offDerCohesive[0][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[0][f]; 512 prob->offCohesive[1][f] = (prob->cohesive[f] ? 0 : Nc) + prob->offCohesive[0][f]; 513 prob->offDerCohesive[1][f] = (prob->cohesive[f] ? 0 : Nc) * dimEmbed + prob->offDerCohesive[0][f]; 514 prob->offCohesive[2][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[2][f]; 515 prob->offDerCohesive[2][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[2][f]; 516 if (q) PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL)); 517 NqMax = PetscMax(NqMax, Nq); 518 NbMax = PetscMax(NbMax, Nb); 519 NcMax = PetscMax(NcMax, Nc); 520 prob->totDim += Nb; 521 prob->totComp += Nc; 522 /* There are two faces for all fields on a cohesive cell, except for cohesive fields */ 523 if (prob->isCohesive && !prob->cohesive[f]) prob->totDim += Nb; 524 } 525 prob->offCohesive[1][Nf] = prob->offCohesive[0][Nf]; 526 prob->offDerCohesive[1][Nf] = prob->offDerCohesive[0][Nf]; 527 /* Allocate works space */ 528 NsMax = 2; /* A non-cohesive discretizations can be used on a cohesive cell, so we need this extra workspace for all DS */ 529 PetscCall(PetscMalloc3(NsMax * prob->totComp, &prob->u, NsMax * prob->totComp, &prob->u_t, NsMax * prob->totComp * dimEmbed + (hasH ? NsMax * prob->totComp * dimEmbed * dimEmbed : 0), &prob->u_x)); 530 PetscCall(PetscMalloc5(dimEmbed, &prob->x, NbMax * NcMax, &prob->basisReal, NbMax * NcMax * dimEmbed, &prob->basisDerReal, NbMax * NcMax, &prob->testReal, NbMax * NcMax * dimEmbed, &prob->testDerReal)); 531 PetscCall(PetscMalloc6(NsMax * NqMax * NcMax, &prob->f0, NsMax * NqMax * NcMax * dimEmbed, &prob->f1, NsMax * NsMax * NqMax * NcMax * NcMax, &prob->g0, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed, &prob->g1, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed, 532 &prob->g2, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed * dimEmbed, &prob->g3)); 533 PetscTryTypeMethod(prob, setup); 534 prob->setup = PETSC_TRUE; 535 PetscFunctionReturn(PETSC_SUCCESS); 536 } 537 538 static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob) 539 { 540 PetscFunctionBegin; 541 PetscCall(PetscFree2(prob->Nc, prob->Nb)); 542 PetscCall(PetscFree2(prob->off, prob->offDer)); 543 PetscCall(PetscFree6(prob->offCohesive[0], prob->offCohesive[1], prob->offCohesive[2], prob->offDerCohesive[0], prob->offDerCohesive[1], prob->offDerCohesive[2])); 544 PetscCall(PetscFree2(prob->T, prob->Tf)); 545 PetscCall(PetscFree3(prob->u, prob->u_t, prob->u_x)); 546 PetscCall(PetscFree5(prob->x, prob->basisReal, prob->basisDerReal, prob->testReal, prob->testDerReal)); 547 PetscCall(PetscFree6(prob->f0, prob->f1, prob->g0, prob->g1, prob->g2, prob->g3)); 548 PetscFunctionReturn(PETSC_SUCCESS); 549 } 550 551 static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew) 552 { 553 PetscObject *tmpd; 554 PetscBool *tmpi; 555 PetscInt *tmpk; 556 PetscBool *tmpc; 557 PetscPointFunc *tmpup; 558 PetscSimplePointFn **tmpexactSol, **tmpexactSol_t; 559 void **tmpexactCtx, **tmpexactCtx_t; 560 void **tmpctx; 561 PetscInt Nf = prob->Nf, f; 562 563 PetscFunctionBegin; 564 if (Nf >= NfNew) PetscFunctionReturn(PETSC_SUCCESS); 565 prob->setup = PETSC_FALSE; 566 PetscCall(PetscDSDestroyStructs_Static(prob)); 567 PetscCall(PetscMalloc4(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpc, NfNew, &tmpk)); 568 for (f = 0; f < Nf; ++f) { 569 tmpd[f] = prob->disc[f]; 570 tmpi[f] = prob->implicit[f]; 571 tmpc[f] = prob->cohesive[f]; 572 tmpk[f] = prob->jetDegree[f]; 573 } 574 for (f = Nf; f < NfNew; ++f) { 575 tmpd[f] = NULL; 576 tmpi[f] = PETSC_TRUE, tmpc[f] = PETSC_FALSE; 577 tmpk[f] = 1; 578 } 579 PetscCall(PetscFree4(prob->disc, prob->implicit, prob->cohesive, prob->jetDegree)); 580 PetscCall(PetscWeakFormSetNumFields(prob->wf, NfNew)); 581 prob->Nf = NfNew; 582 prob->disc = tmpd; 583 prob->implicit = tmpi; 584 prob->cohesive = tmpc; 585 prob->jetDegree = tmpk; 586 PetscCall(PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx)); 587 for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f]; 588 for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f]; 589 for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL; 590 for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL; 591 PetscCall(PetscFree2(prob->update, prob->ctx)); 592 prob->update = tmpup; 593 prob->ctx = tmpctx; 594 PetscCall(PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t)); 595 for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f]; 596 for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f]; 597 for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f]; 598 for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f]; 599 for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL; 600 for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL; 601 for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL; 602 for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL; 603 PetscCall(PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t)); 604 prob->exactSol = tmpexactSol; 605 prob->exactCtx = tmpexactCtx; 606 prob->exactSol_t = tmpexactSol_t; 607 prob->exactCtx_t = tmpexactCtx_t; 608 PetscFunctionReturn(PETSC_SUCCESS); 609 } 610 611 /*@ 612 PetscDSDestroy - Destroys a `PetscDS` object 613 614 Collective 615 616 Input Parameter: 617 . ds - the `PetscDS` object to destroy 618 619 Level: developer 620 621 .seealso: `PetscDSView()` 622 @*/ 623 PetscErrorCode PetscDSDestroy(PetscDS *ds) 624 { 625 PetscInt f; 626 627 PetscFunctionBegin; 628 if (!*ds) PetscFunctionReturn(PETSC_SUCCESS); 629 PetscValidHeaderSpecific(*ds, PETSCDS_CLASSID, 1); 630 631 if (--((PetscObject)*ds)->refct > 0) { 632 *ds = NULL; 633 PetscFunctionReturn(PETSC_SUCCESS); 634 } 635 ((PetscObject)*ds)->refct = 0; 636 if ((*ds)->subprobs) { 637 PetscInt dim, d; 638 639 PetscCall(PetscDSGetSpatialDimension(*ds, &dim)); 640 for (d = 0; d < dim; ++d) PetscCall(PetscDSDestroy(&(*ds)->subprobs[d])); 641 } 642 PetscCall(PetscFree((*ds)->subprobs)); 643 PetscCall(PetscDSDestroyStructs_Static(*ds)); 644 for (f = 0; f < (*ds)->Nf; ++f) PetscCall(PetscObjectDereference((*ds)->disc[f])); 645 PetscCall(PetscFree4((*ds)->disc, (*ds)->implicit, (*ds)->cohesive, (*ds)->jetDegree)); 646 PetscCall(PetscWeakFormDestroy(&(*ds)->wf)); 647 PetscCall(PetscFree2((*ds)->update, (*ds)->ctx)); 648 PetscCall(PetscFree4((*ds)->exactSol, (*ds)->exactCtx, (*ds)->exactSol_t, (*ds)->exactCtx_t)); 649 PetscTryTypeMethod(*ds, destroy); 650 PetscCall(PetscDSDestroyBoundary(*ds)); 651 PetscCall(PetscFree((*ds)->constants)); 652 for (PetscInt c = 0; c < DM_NUM_POLYTOPES; ++c) { 653 const PetscInt Na = DMPolytopeTypeGetNumArrangements((DMPolytopeType)c); 654 if ((*ds)->quadPerm[c]) 655 for (PetscInt o = 0; o < Na; ++o) PetscCall(ISDestroy(&(*ds)->quadPerm[c][o])); 656 PetscCall(PetscFree((*ds)->quadPerm[c])); 657 (*ds)->quadPerm[c] = NULL; 658 } 659 PetscCall(PetscHeaderDestroy(ds)); 660 PetscFunctionReturn(PETSC_SUCCESS); 661 } 662 663 /*@ 664 PetscDSCreate - Creates an empty `PetscDS` object. The type can then be set with `PetscDSSetType()`. 665 666 Collective 667 668 Input Parameter: 669 . comm - The communicator for the `PetscDS` object 670 671 Output Parameter: 672 . ds - The `PetscDS` object 673 674 Level: beginner 675 676 .seealso: `PetscDS`, `PetscDSSetType()`, `PETSCDSBASIC`, `PetscDSType` 677 @*/ 678 PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds) 679 { 680 PetscDS p; 681 682 PetscFunctionBegin; 683 PetscAssertPointer(ds, 2); 684 PetscCall(PetscDSInitializePackage()); 685 686 PetscCall(PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView)); 687 p->Nf = 0; 688 p->setup = PETSC_FALSE; 689 p->numConstants = 0; 690 p->numFuncConstants = 3; // Row and col fields, cell size 691 p->dimEmbed = -1; 692 p->useJacPre = PETSC_TRUE; 693 p->forceQuad = PETSC_TRUE; 694 PetscCall(PetscMalloc1(p->numConstants + p->numFuncConstants, &p->constants)); 695 PetscCall(PetscWeakFormCreate(comm, &p->wf)); 696 PetscCall(PetscArrayzero(p->quadPerm, DM_NUM_POLYTOPES)); 697 *ds = p; 698 PetscFunctionReturn(PETSC_SUCCESS); 699 } 700 701 /*@ 702 PetscDSGetNumFields - Returns the number of fields in the `PetscDS` 703 704 Not Collective 705 706 Input Parameter: 707 . prob - The `PetscDS` object 708 709 Output Parameter: 710 . Nf - The number of fields 711 712 Level: beginner 713 714 .seealso: `PetscDS`, `PetscDSGetSpatialDimension()`, `PetscDSCreate()` 715 @*/ 716 PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf) 717 { 718 PetscFunctionBegin; 719 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 720 PetscAssertPointer(Nf, 2); 721 *Nf = prob->Nf; 722 PetscFunctionReturn(PETSC_SUCCESS); 723 } 724 725 /*@ 726 PetscDSGetSpatialDimension - Returns the spatial dimension of the `PetscDS`, meaning the topological dimension of the discretizations 727 728 Not Collective 729 730 Input Parameter: 731 . prob - The `PetscDS` object 732 733 Output Parameter: 734 . dim - The spatial dimension 735 736 Level: beginner 737 738 .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 739 @*/ 740 PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim) 741 { 742 PetscFunctionBegin; 743 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 744 PetscAssertPointer(dim, 2); 745 *dim = 0; 746 if (prob->Nf) { 747 PetscObject obj; 748 PetscClassId id; 749 750 PetscCall(PetscDSGetDiscretization(prob, 0, &obj)); 751 if (obj) { 752 PetscCall(PetscObjectGetClassId(obj, &id)); 753 if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetSpatialDimension((PetscFE)obj, dim)); 754 else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetSpatialDimension((PetscFV)obj, dim)); 755 else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0); 756 } 757 } 758 PetscFunctionReturn(PETSC_SUCCESS); 759 } 760 761 /*@ 762 PetscDSGetCoordinateDimension - Returns the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded 763 764 Not Collective 765 766 Input Parameter: 767 . prob - The `PetscDS` object 768 769 Output Parameter: 770 . dimEmbed - The coordinate dimension 771 772 Level: beginner 773 774 .seealso: `PetscDS`, `PetscDSSetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 775 @*/ 776 PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed) 777 { 778 PetscFunctionBegin; 779 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 780 PetscAssertPointer(dimEmbed, 2); 781 PetscCheck(prob->dimEmbed >= 0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS"); 782 *dimEmbed = prob->dimEmbed; 783 PetscFunctionReturn(PETSC_SUCCESS); 784 } 785 786 /*@ 787 PetscDSSetCoordinateDimension - Set the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded 788 789 Logically Collective 790 791 Input Parameters: 792 + prob - The `PetscDS` object 793 - dimEmbed - The coordinate dimension 794 795 Level: beginner 796 797 .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 798 @*/ 799 PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed) 800 { 801 PetscFunctionBegin; 802 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 803 PetscCheck(dimEmbed >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %" PetscInt_FMT, dimEmbed); 804 prob->dimEmbed = dimEmbed; 805 PetscFunctionReturn(PETSC_SUCCESS); 806 } 807 808 /*@ 809 PetscDSGetForceQuad - Returns the flag to force matching quadratures among the field discretizations 810 811 Not collective 812 813 Input Parameter: 814 . ds - The `PetscDS` object 815 816 Output Parameter: 817 . forceQuad - The flag 818 819 Level: intermediate 820 821 .seealso: `PetscDS`, `PetscDSSetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 822 @*/ 823 PetscErrorCode PetscDSGetForceQuad(PetscDS ds, PetscBool *forceQuad) 824 { 825 PetscFunctionBegin; 826 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 827 PetscAssertPointer(forceQuad, 2); 828 *forceQuad = ds->forceQuad; 829 PetscFunctionReturn(PETSC_SUCCESS); 830 } 831 832 /*@ 833 PetscDSSetForceQuad - Set the flag to force matching quadratures among the field discretizations 834 835 Logically collective on ds 836 837 Input Parameters: 838 + ds - The `PetscDS` object 839 - forceQuad - The flag 840 841 Level: intermediate 842 843 .seealso: `PetscDS`, `PetscDSGetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 844 @*/ 845 PetscErrorCode PetscDSSetForceQuad(PetscDS ds, PetscBool forceQuad) 846 { 847 PetscFunctionBegin; 848 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 849 ds->forceQuad = forceQuad; 850 PetscFunctionReturn(PETSC_SUCCESS); 851 } 852 853 /*@ 854 PetscDSIsCohesive - Returns the flag indicating that this `PetscDS` is for a cohesive cell 855 856 Not Collective 857 858 Input Parameter: 859 . ds - The `PetscDS` object 860 861 Output Parameter: 862 . isCohesive - The flag 863 864 Level: developer 865 866 .seealso: `PetscDS`, `PetscDSGetNumCohesive()`, `PetscDSGetCohesive()`, `PetscDSSetCohesive()`, `PetscDSCreate()` 867 @*/ 868 PetscErrorCode PetscDSIsCohesive(PetscDS ds, PetscBool *isCohesive) 869 { 870 PetscFunctionBegin; 871 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 872 PetscAssertPointer(isCohesive, 2); 873 *isCohesive = ds->isCohesive; 874 PetscFunctionReturn(PETSC_SUCCESS); 875 } 876 877 /*@ 878 PetscDSGetNumCohesive - Returns the number of cohesive fields, meaning those defined on the interior of a cohesive cell 879 880 Not Collective 881 882 Input Parameter: 883 . ds - The `PetscDS` object 884 885 Output Parameter: 886 . numCohesive - The number of cohesive fields 887 888 Level: developer 889 890 .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSCreate()` 891 @*/ 892 PetscErrorCode PetscDSGetNumCohesive(PetscDS ds, PetscInt *numCohesive) 893 { 894 PetscInt f; 895 896 PetscFunctionBegin; 897 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 898 PetscAssertPointer(numCohesive, 2); 899 *numCohesive = 0; 900 for (f = 0; f < ds->Nf; ++f) *numCohesive += ds->cohesive[f] ? 1 : 0; 901 PetscFunctionReturn(PETSC_SUCCESS); 902 } 903 904 /*@ 905 PetscDSGetCohesive - Returns the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell 906 907 Not Collective 908 909 Input Parameters: 910 + ds - The `PetscDS` object 911 - f - The field index 912 913 Output Parameter: 914 . isCohesive - The flag 915 916 Level: developer 917 918 .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()` 919 @*/ 920 PetscErrorCode PetscDSGetCohesive(PetscDS ds, PetscInt f, PetscBool *isCohesive) 921 { 922 PetscFunctionBegin; 923 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 924 PetscAssertPointer(isCohesive, 3); 925 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 926 *isCohesive = ds->cohesive[f]; 927 PetscFunctionReturn(PETSC_SUCCESS); 928 } 929 930 /*@ 931 PetscDSSetCohesive - Set the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell 932 933 Not Collective 934 935 Input Parameters: 936 + ds - The `PetscDS` object 937 . f - The field index 938 - isCohesive - The flag for a cohesive field 939 940 Level: developer 941 942 .seealso: `PetscDS`, `PetscDSGetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()` 943 @*/ 944 PetscErrorCode PetscDSSetCohesive(PetscDS ds, PetscInt f, PetscBool isCohesive) 945 { 946 PetscInt i; 947 948 PetscFunctionBegin; 949 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 950 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 951 ds->cohesive[f] = isCohesive; 952 ds->isCohesive = PETSC_FALSE; 953 for (i = 0; i < ds->Nf; ++i) ds->isCohesive = ds->isCohesive || ds->cohesive[f] ? PETSC_TRUE : PETSC_FALSE; 954 PetscFunctionReturn(PETSC_SUCCESS); 955 } 956 957 /*@ 958 PetscDSGetTotalDimension - Returns the total size of the approximation space for this system 959 960 Not Collective 961 962 Input Parameter: 963 . prob - The `PetscDS` object 964 965 Output Parameter: 966 . dim - The total problem dimension 967 968 Level: beginner 969 970 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()` 971 @*/ 972 PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim) 973 { 974 PetscFunctionBegin; 975 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 976 PetscCall(PetscDSSetUp(prob)); 977 PetscAssertPointer(dim, 2); 978 *dim = prob->totDim; 979 PetscFunctionReturn(PETSC_SUCCESS); 980 } 981 982 /*@ 983 PetscDSGetTotalComponents - Returns the total number of components in this system 984 985 Not Collective 986 987 Input Parameter: 988 . prob - The `PetscDS` object 989 990 Output Parameter: 991 . Nc - The total number of components 992 993 Level: beginner 994 995 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()` 996 @*/ 997 PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc) 998 { 999 PetscFunctionBegin; 1000 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1001 PetscCall(PetscDSSetUp(prob)); 1002 PetscAssertPointer(Nc, 2); 1003 *Nc = prob->totComp; 1004 PetscFunctionReturn(PETSC_SUCCESS); 1005 } 1006 1007 /*@ 1008 PetscDSGetDiscretization - Returns the discretization object for the given field 1009 1010 Not Collective 1011 1012 Input Parameters: 1013 + prob - The `PetscDS` object 1014 - f - The field number 1015 1016 Output Parameter: 1017 . disc - The discretization object 1018 1019 Level: beginner 1020 1021 .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1022 @*/ 1023 PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc) 1024 { 1025 PetscFunctionBeginHot; 1026 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1027 PetscAssertPointer(disc, 3); 1028 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 1029 *disc = prob->disc[f]; 1030 PetscFunctionReturn(PETSC_SUCCESS); 1031 } 1032 1033 /*@ 1034 PetscDSSetDiscretization - Sets the discretization object for the given field 1035 1036 Not Collective 1037 1038 Input Parameters: 1039 + prob - The `PetscDS` object 1040 . f - The field number 1041 - disc - The discretization object 1042 1043 Level: beginner 1044 1045 .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSGetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1046 @*/ 1047 PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc) 1048 { 1049 PetscFunctionBegin; 1050 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1051 if (disc) PetscAssertPointer(disc, 3); 1052 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 1053 PetscCall(PetscDSEnlarge_Static(prob, f + 1)); 1054 PetscCall(PetscObjectDereference(prob->disc[f])); 1055 prob->disc[f] = disc; 1056 PetscCall(PetscObjectReference(disc)); 1057 if (disc) { 1058 PetscClassId id; 1059 1060 PetscCall(PetscObjectGetClassId(disc, &id)); 1061 if (id == PETSCFE_CLASSID) { 1062 PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE)); 1063 } else if (id == PETSCFV_CLASSID) { 1064 PetscCall(PetscDSSetImplicit(prob, f, PETSC_FALSE)); 1065 } 1066 PetscCall(PetscDSSetJetDegree(prob, f, 1)); 1067 } 1068 PetscFunctionReturn(PETSC_SUCCESS); 1069 } 1070 1071 /*@ 1072 PetscDSGetWeakForm - Returns the weak form object 1073 1074 Not Collective 1075 1076 Input Parameter: 1077 . ds - The `PetscDS` object 1078 1079 Output Parameter: 1080 . wf - The weak form object 1081 1082 Level: beginner 1083 1084 .seealso: `PetscWeakForm`, `PetscDSSetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1085 @*/ 1086 PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf) 1087 { 1088 PetscFunctionBegin; 1089 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1090 PetscAssertPointer(wf, 2); 1091 *wf = ds->wf; 1092 PetscFunctionReturn(PETSC_SUCCESS); 1093 } 1094 1095 /*@ 1096 PetscDSSetWeakForm - Sets the weak form object 1097 1098 Not Collective 1099 1100 Input Parameters: 1101 + ds - The `PetscDS` object 1102 - wf - The weak form object 1103 1104 Level: beginner 1105 1106 .seealso: `PetscWeakForm`, `PetscDSGetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1107 @*/ 1108 PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf) 1109 { 1110 PetscFunctionBegin; 1111 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1112 PetscValidHeaderSpecific(wf, PETSCWEAKFORM_CLASSID, 2); 1113 PetscCall(PetscObjectDereference((PetscObject)ds->wf)); 1114 ds->wf = wf; 1115 PetscCall(PetscObjectReference((PetscObject)wf)); 1116 PetscCall(PetscWeakFormSetNumFields(wf, ds->Nf)); 1117 PetscFunctionReturn(PETSC_SUCCESS); 1118 } 1119 1120 /*@ 1121 PetscDSAddDiscretization - Adds a discretization object 1122 1123 Not Collective 1124 1125 Input Parameters: 1126 + prob - The `PetscDS` object 1127 - disc - The boundary discretization object 1128 1129 Level: beginner 1130 1131 .seealso: `PetscWeakForm`, `PetscDSGetDiscretization()`, `PetscDSSetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1132 @*/ 1133 PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc) 1134 { 1135 PetscFunctionBegin; 1136 PetscCall(PetscDSSetDiscretization(prob, prob->Nf, disc)); 1137 PetscFunctionReturn(PETSC_SUCCESS); 1138 } 1139 1140 /*@ 1141 PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the `PetscDS` 1142 1143 Not Collective 1144 1145 Input Parameter: 1146 . prob - The `PetscDS` object 1147 1148 Output Parameter: 1149 . q - The quadrature object 1150 1151 Level: intermediate 1152 1153 .seealso: `PetscDS`, `PetscQuadrature`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1154 @*/ 1155 PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q) 1156 { 1157 PetscObject obj; 1158 PetscClassId id; 1159 1160 PetscFunctionBegin; 1161 *q = NULL; 1162 if (!prob->Nf) PetscFunctionReturn(PETSC_SUCCESS); 1163 PetscCall(PetscDSGetDiscretization(prob, 0, &obj)); 1164 PetscCall(PetscObjectGetClassId(obj, &id)); 1165 if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetQuadrature((PetscFE)obj, q)); 1166 else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetQuadrature((PetscFV)obj, q)); 1167 else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0); 1168 PetscFunctionReturn(PETSC_SUCCESS); 1169 } 1170 1171 /*@ 1172 PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for `TSARKIMEX` 1173 1174 Not Collective 1175 1176 Input Parameters: 1177 + prob - The `PetscDS` object 1178 - f - The field number 1179 1180 Output Parameter: 1181 . implicit - The flag indicating what kind of solve to use for this field 1182 1183 Level: developer 1184 1185 .seealso: `TSARKIMEX`, `PetscDS`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1186 @*/ 1187 PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit) 1188 { 1189 PetscFunctionBegin; 1190 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1191 PetscAssertPointer(implicit, 3); 1192 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 1193 *implicit = prob->implicit[f]; 1194 PetscFunctionReturn(PETSC_SUCCESS); 1195 } 1196 1197 /*@ 1198 PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for `TSARKIMEX` 1199 1200 Not Collective 1201 1202 Input Parameters: 1203 + prob - The `PetscDS` object 1204 . f - The field number 1205 - implicit - The flag indicating what kind of solve to use for this field 1206 1207 Level: developer 1208 1209 .seealso: `TSARKIMEX`, `PetscDSGetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1210 @*/ 1211 PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit) 1212 { 1213 PetscFunctionBegin; 1214 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1215 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 1216 prob->implicit[f] = implicit; 1217 PetscFunctionReturn(PETSC_SUCCESS); 1218 } 1219 1220 /*@ 1221 PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate. 1222 1223 Not Collective 1224 1225 Input Parameters: 1226 + ds - The `PetscDS` object 1227 - f - The field number 1228 1229 Output Parameter: 1230 . k - The highest derivative we need to tabulate 1231 1232 Level: developer 1233 1234 .seealso: `PetscDS`, `PetscDSSetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1235 @*/ 1236 PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k) 1237 { 1238 PetscFunctionBegin; 1239 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1240 PetscAssertPointer(k, 3); 1241 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 1242 *k = ds->jetDegree[f]; 1243 PetscFunctionReturn(PETSC_SUCCESS); 1244 } 1245 1246 /*@ 1247 PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate. 1248 1249 Not Collective 1250 1251 Input Parameters: 1252 + ds - The `PetscDS` object 1253 . f - The field number 1254 - k - The highest derivative we need to tabulate 1255 1256 Level: developer 1257 1258 .seealso: `PetscDS`, `PetscDSGetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 1259 @*/ 1260 PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k) 1261 { 1262 PetscFunctionBegin; 1263 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1264 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 1265 ds->jetDegree[f] = k; 1266 PetscFunctionReturn(PETSC_SUCCESS); 1267 } 1268 1269 /*@C 1270 PetscDSGetObjective - Get the pointwise objective function for a given test field 1271 1272 Not Collective 1273 1274 Input Parameters: 1275 + ds - The `PetscDS` 1276 - f - The test field number 1277 1278 Output Parameter: 1279 . obj - integrand for the test function term 1280 1281 Calling sequence of `obj`: 1282 + dim - the spatial dimension 1283 . Nf - the number of fields 1284 . NfAux - the number of auxiliary fields 1285 . uOff - the offset into u[] and u_t[] for each field 1286 . uOff_x - the offset into u_x[] for each field 1287 . u - each field evaluated at the current point 1288 . u_t - the time derivative of each field evaluated at the current point 1289 . u_x - the gradient of each field evaluated at the current point 1290 . aOff - the offset into a[] and a_t[] for each auxiliary field 1291 . aOff_x - the offset into a_x[] for each auxiliary field 1292 . a - each auxiliary field evaluated at the current point 1293 . a_t - the time derivative of each auxiliary field evaluated at the current point 1294 . a_x - the gradient of auxiliary each field evaluated at the current point 1295 . t - current time 1296 . x - coordinates of the current point 1297 . numConstants - number of constant parameters 1298 . constants - constant parameters 1299 - obj - output values at the current point 1300 1301 Level: intermediate 1302 1303 Note: 1304 We are using a first order FEM model for the weak form\: $ \int_\Omega \phi\,\mathrm{obj}(u, u_t, \nabla u, x, t)$ 1305 1306 .seealso: `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()` 1307 @*/ 1308 PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f, void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])) 1309 { 1310 PetscPointFunc *tmp; 1311 PetscInt n; 1312 1313 PetscFunctionBegin; 1314 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1315 PetscAssertPointer(obj, 3); 1316 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 1317 PetscCall(PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp)); 1318 *obj = tmp ? tmp[0] : NULL; 1319 PetscFunctionReturn(PETSC_SUCCESS); 1320 } 1321 1322 /*@C 1323 PetscDSSetObjective - Set the pointwise objective function for a given test field 1324 1325 Not Collective 1326 1327 Input Parameters: 1328 + ds - The `PetscDS` 1329 . f - The test field number 1330 - obj - integrand for the test function term 1331 1332 Calling sequence of `obj`: 1333 + dim - the spatial dimension 1334 . Nf - the number of fields 1335 . NfAux - the number of auxiliary fields 1336 . uOff - the offset into u[] and u_t[] for each field 1337 . uOff_x - the offset into u_x[] for each field 1338 . u - each field evaluated at the current point 1339 . u_t - the time derivative of each field evaluated at the current point 1340 . u_x - the gradient of each field evaluated at the current point 1341 . aOff - the offset into a[] and a_t[] for each auxiliary field 1342 . aOff_x - the offset into a_x[] for each auxiliary field 1343 . a - each auxiliary field evaluated at the current point 1344 . a_t - the time derivative of each auxiliary field evaluated at the current point 1345 . a_x - the gradient of auxiliary each field evaluated at the current point 1346 . t - current time 1347 . x - coordinates of the current point 1348 . numConstants - number of constant parameters 1349 . constants - constant parameters 1350 - obj - output values at the current point 1351 1352 Level: intermediate 1353 1354 Note: 1355 We are using a first order FEM model for the weak form\: $ \int_\Omega \phi\,\mathrm{obj}(u, u_t, \nabla u, x, t)$ 1356 1357 .seealso: `PetscDS`, `PetscDSGetObjective()`, `PetscDSSetResidual()` 1358 @*/ 1359 PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f, void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])) 1360 { 1361 PetscFunctionBegin; 1362 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1363 if (obj) PetscValidFunction(obj, 3); 1364 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 1365 PetscCall(PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj)); 1366 PetscFunctionReturn(PETSC_SUCCESS); 1367 } 1368 1369 /*@C 1370 PetscDSGetResidual - Get the pointwise residual function for a given test field 1371 1372 Not Collective 1373 1374 Input Parameters: 1375 + ds - The `PetscDS` 1376 - f - The test field number 1377 1378 Output Parameters: 1379 + f0 - integrand for the test function term 1380 - f1 - integrand for the test function gradient term 1381 1382 Calling sequence of `f0`: 1383 + dim - the spatial dimension 1384 . Nf - the number of fields 1385 . NfAux - the number of auxiliary fields 1386 . uOff - the offset into u[] and u_t[] for each field 1387 . uOff_x - the offset into u_x[] for each field 1388 . u - each field evaluated at the current point 1389 . u_t - the time derivative of each field evaluated at the current point 1390 . u_x - the gradient of each field evaluated at the current point 1391 . aOff - the offset into a[] and a_t[] for each auxiliary field 1392 . aOff_x - the offset into a_x[] for each auxiliary field 1393 . a - each auxiliary field evaluated at the current point 1394 . a_t - the time derivative of each auxiliary field evaluated at the current point 1395 . a_x - the gradient of auxiliary each field evaluated at the current point 1396 . t - current time 1397 . x - coordinates of the current point 1398 . numConstants - number of constant parameters 1399 . constants - constant parameters 1400 - f0 - output values at the current point 1401 1402 Level: intermediate 1403 1404 Note: 1405 `f1` has an identical form and is omitted for brevity. 1406 1407 We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$ 1408 1409 .seealso: `PetscDS`, `PetscDSSetResidual()` 1410 @*/ 1411 PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1412 { 1413 PetscPointFunc *tmp0, *tmp1; 1414 PetscInt n0, n1; 1415 1416 PetscFunctionBegin; 1417 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1418 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 1419 PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1)); 1420 *f0 = tmp0 ? tmp0[0] : NULL; 1421 *f1 = tmp1 ? tmp1[0] : NULL; 1422 PetscFunctionReturn(PETSC_SUCCESS); 1423 } 1424 1425 /*@C 1426 PetscDSSetResidual - Set the pointwise residual function for a given test field 1427 1428 Not Collective 1429 1430 Input Parameters: 1431 + ds - The `PetscDS` 1432 . f - The test field number 1433 . f0 - integrand for the test function term 1434 - f1 - integrand for the test function gradient term 1435 1436 Calling sequence of `f0`: 1437 + dim - the spatial dimension 1438 . Nf - the number of fields 1439 . NfAux - the number of auxiliary fields 1440 . uOff - the offset into u[] and u_t[] for each field 1441 . uOff_x - the offset into u_x[] for each field 1442 . u - each field evaluated at the current point 1443 . u_t - the time derivative of each field evaluated at the current point 1444 . u_x - the gradient of each field evaluated at the current point 1445 . aOff - the offset into a[] and a_t[] for each auxiliary field 1446 . aOff_x - the offset into a_x[] for each auxiliary field 1447 . a - each auxiliary field evaluated at the current point 1448 . a_t - the time derivative of each auxiliary field evaluated at the current point 1449 . a_x - the gradient of auxiliary each field evaluated at the current point 1450 . t - current time 1451 . x - coordinates of the current point 1452 . numConstants - number of constant parameters 1453 . constants - constant parameters 1454 - f0 - output values at the current point 1455 1456 Level: intermediate 1457 1458 Note: 1459 `f1` has an identical form and is omitted for brevity. 1460 1461 We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$ 1462 1463 .seealso: `PetscDS`, `PetscDSGetResidual()` 1464 @*/ 1465 PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1466 { 1467 PetscFunctionBegin; 1468 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1469 if (f0) PetscValidFunction(f0, 3); 1470 if (f1) PetscValidFunction(f1, 4); 1471 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 1472 PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1)); 1473 PetscFunctionReturn(PETSC_SUCCESS); 1474 } 1475 1476 /*@C 1477 PetscDSGetRHSResidual - Get the pointwise RHS residual function for explicit timestepping for a given test field 1478 1479 Not Collective 1480 1481 Input Parameters: 1482 + ds - The `PetscDS` 1483 - f - The test field number 1484 1485 Output Parameters: 1486 + f0 - integrand for the test function term 1487 - f1 - integrand for the test function gradient term 1488 1489 Calling sequence of `f0`: 1490 + dim - the spatial dimension 1491 . Nf - the number of fields 1492 . NfAux - the number of auxiliary fields 1493 . uOff - the offset into u[] and u_t[] for each field 1494 . uOff_x - the offset into u_x[] for each field 1495 . u - each field evaluated at the current point 1496 . u_t - the time derivative of each field evaluated at the current point 1497 . u_x - the gradient of each field evaluated at the current point 1498 . aOff - the offset into a[] and a_t[] for each auxiliary field 1499 . aOff_x - the offset into a_x[] for each auxiliary field 1500 . a - each auxiliary field evaluated at the current point 1501 . a_t - the time derivative of each auxiliary field evaluated at the current point 1502 . a_x - the gradient of auxiliary each field evaluated at the current point 1503 . t - current time 1504 . x - coordinates of the current point 1505 . numConstants - number of constant parameters 1506 . constants - constant parameters 1507 - f0 - output values at the current point 1508 1509 Level: intermediate 1510 1511 Note: 1512 `f1` has an identical form and is omitted for brevity. 1513 1514 We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$ 1515 1516 .seealso: `PetscDS`, `PetscDSSetRHSResidual()` 1517 @*/ 1518 PetscErrorCode PetscDSGetRHSResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1519 { 1520 PetscPointFunc *tmp0, *tmp1; 1521 PetscInt n0, n1; 1522 1523 PetscFunctionBegin; 1524 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1525 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 1526 PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 100, &n0, &tmp0, &n1, &tmp1)); 1527 *f0 = tmp0 ? tmp0[0] : NULL; 1528 *f1 = tmp1 ? tmp1[0] : NULL; 1529 PetscFunctionReturn(PETSC_SUCCESS); 1530 } 1531 1532 /*@C 1533 PetscDSSetRHSResidual - Set the pointwise residual function for explicit timestepping for a given test field 1534 1535 Not Collective 1536 1537 Input Parameters: 1538 + ds - The `PetscDS` 1539 . f - The test field number 1540 . f0 - integrand for the test function term 1541 - f1 - integrand for the test function gradient term 1542 1543 Calling sequence for the callbacks `f0`: 1544 + dim - the spatial dimension 1545 . Nf - the number of fields 1546 . NfAux - the number of auxiliary fields 1547 . uOff - the offset into u[] and u_t[] for each field 1548 . uOff_x - the offset into u_x[] for each field 1549 . u - each field evaluated at the current point 1550 . u_t - the time derivative of each field evaluated at the current point 1551 . u_x - the gradient of each field evaluated at the current point 1552 . aOff - the offset into a[] and a_t[] for each auxiliary field 1553 . aOff_x - the offset into a_x[] for each auxiliary field 1554 . a - each auxiliary field evaluated at the current point 1555 . a_t - the time derivative of each auxiliary field evaluated at the current point 1556 . a_x - the gradient of auxiliary each field evaluated at the current point 1557 . t - current time 1558 . x - coordinates of the current point 1559 . numConstants - number of constant parameters 1560 . constants - constant parameters 1561 - f0 - output values at the current point 1562 1563 Level: intermediate 1564 1565 Note: 1566 `f1` has an identical form and is omitted for brevity. 1567 1568 We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$ 1569 1570 .seealso: `PetscDS`, `PetscDSGetResidual()` 1571 @*/ 1572 PetscErrorCode PetscDSSetRHSResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1573 { 1574 PetscFunctionBegin; 1575 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1576 if (f0) PetscValidFunction(f0, 3); 1577 if (f1) PetscValidFunction(f1, 4); 1578 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 1579 PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1)); 1580 PetscFunctionReturn(PETSC_SUCCESS); 1581 } 1582 1583 /*@ 1584 PetscDSHasJacobian - Checks that the Jacobian functions have been set 1585 1586 Not Collective 1587 1588 Input Parameter: 1589 . ds - The `PetscDS` 1590 1591 Output Parameter: 1592 . hasJac - flag that pointwise function for the Jacobian has been set 1593 1594 Level: intermediate 1595 1596 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()` 1597 @*/ 1598 PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac) 1599 { 1600 PetscFunctionBegin; 1601 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1602 PetscCall(PetscWeakFormHasJacobian(ds->wf, hasJac)); 1603 PetscFunctionReturn(PETSC_SUCCESS); 1604 } 1605 1606 /*@C 1607 PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field 1608 1609 Not Collective 1610 1611 Input Parameters: 1612 + ds - The `PetscDS` 1613 . f - The test field number 1614 - g - The field number 1615 1616 Output Parameters: 1617 + g0 - integrand for the test and basis function term 1618 . g1 - integrand for the test function and basis function gradient term 1619 . g2 - integrand for the test function gradient and basis function term 1620 - g3 - integrand for the test function gradient and basis function gradient term 1621 1622 Calling sequence of `g0`: 1623 + dim - the spatial dimension 1624 . Nf - the number of fields 1625 . NfAux - the number of auxiliary fields 1626 . uOff - the offset into u[] and u_t[] for each field 1627 . uOff_x - the offset into u_x[] for each field 1628 . u - each field evaluated at the current point 1629 . u_t - the time derivative of each field evaluated at the current point 1630 . u_x - the gradient of each field evaluated at the current point 1631 . aOff - the offset into a[] and a_t[] for each auxiliary field 1632 . aOff_x - the offset into a_x[] for each auxiliary field 1633 . a - each auxiliary field evaluated at the current point 1634 . a_t - the time derivative of each auxiliary field evaluated at the current point 1635 . a_x - the gradient of auxiliary each field evaluated at the current point 1636 . t - current time 1637 . u_tShift - the multiplier a for dF/dU_t 1638 . x - coordinates of the current point 1639 . numConstants - number of constant parameters 1640 . constants - constant parameters 1641 - g0 - output values at the current point 1642 1643 Level: intermediate 1644 1645 Note: 1646 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 1647 1648 We are using a first order FEM model for the weak form\: 1649 1650 $$ 1651 \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi 1652 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi 1653 $$ 1654 1655 .seealso: `PetscDS`, `PetscDSSetJacobian()` 1656 @*/ 1657 PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1658 { 1659 PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3; 1660 PetscInt n0, n1, n2, n3; 1661 1662 PetscFunctionBegin; 1663 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1664 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 1665 PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf); 1666 PetscCall(PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3)); 1667 *g0 = tmp0 ? tmp0[0] : NULL; 1668 *g1 = tmp1 ? tmp1[0] : NULL; 1669 *g2 = tmp2 ? tmp2[0] : NULL; 1670 *g3 = tmp3 ? tmp3[0] : NULL; 1671 PetscFunctionReturn(PETSC_SUCCESS); 1672 } 1673 1674 /*@C 1675 PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields 1676 1677 Not Collective 1678 1679 Input Parameters: 1680 + ds - The `PetscDS` 1681 . f - The test field number 1682 . g - The field number 1683 . g0 - integrand for the test and basis function term 1684 . g1 - integrand for the test function and basis function gradient term 1685 . g2 - integrand for the test function gradient and basis function term 1686 - g3 - integrand for the test function gradient and basis function gradient term 1687 1688 Calling sequence of `g0`: 1689 + dim - the spatial dimension 1690 . Nf - the number of fields 1691 . NfAux - the number of auxiliary fields 1692 . uOff - the offset into u[] and u_t[] for each field 1693 . uOff_x - the offset into u_x[] for each field 1694 . u - each field evaluated at the current point 1695 . u_t - the time derivative of each field evaluated at the current point 1696 . u_x - the gradient of each field evaluated at the current point 1697 . aOff - the offset into a[] and a_t[] for each auxiliary field 1698 . aOff_x - the offset into a_x[] for each auxiliary field 1699 . a - each auxiliary field evaluated at the current point 1700 . a_t - the time derivative of each auxiliary field evaluated at the current point 1701 . a_x - the gradient of auxiliary each field evaluated at the current point 1702 . t - current time 1703 . u_tShift - the multiplier a for dF/dU_t 1704 . x - coordinates of the current point 1705 . numConstants - number of constant parameters 1706 . constants - constant parameters 1707 - g0 - output values at the current point 1708 1709 Level: intermediate 1710 1711 Note: 1712 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 1713 1714 We are using a first order FEM model for the weak form\: 1715 1716 $$ 1717 \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi 1718 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi 1719 $$ 1720 1721 .seealso: `PetscDS`, `PetscDSGetJacobian()` 1722 @*/ 1723 PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1724 { 1725 PetscFunctionBegin; 1726 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1727 if (g0) PetscValidFunction(g0, 4); 1728 if (g1) PetscValidFunction(g1, 5); 1729 if (g2) PetscValidFunction(g2, 6); 1730 if (g3) PetscValidFunction(g3, 7); 1731 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 1732 PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g); 1733 PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3)); 1734 PetscFunctionReturn(PETSC_SUCCESS); 1735 } 1736 1737 /*@ 1738 PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner 1739 1740 Not Collective 1741 1742 Input Parameters: 1743 + prob - The `PetscDS` 1744 - useJacPre - flag that enables construction of a Jacobian preconditioner 1745 1746 Level: intermediate 1747 1748 Developer Note: 1749 Should be called `PetscDSSetUseJacobianPreconditioner()` 1750 1751 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()` 1752 @*/ 1753 PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre) 1754 { 1755 PetscFunctionBegin; 1756 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1757 prob->useJacPre = useJacPre; 1758 PetscFunctionReturn(PETSC_SUCCESS); 1759 } 1760 1761 /*@ 1762 PetscDSHasJacobianPreconditioner - Checks if a Jacobian preconditioner matrix has been set 1763 1764 Not Collective 1765 1766 Input Parameter: 1767 . ds - The `PetscDS` 1768 1769 Output Parameter: 1770 . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set 1771 1772 Level: intermediate 1773 1774 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()` 1775 @*/ 1776 PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre) 1777 { 1778 PetscFunctionBegin; 1779 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1780 *hasJacPre = PETSC_FALSE; 1781 if (!ds->useJacPre) PetscFunctionReturn(PETSC_SUCCESS); 1782 PetscCall(PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre)); 1783 PetscFunctionReturn(PETSC_SUCCESS); 1784 } 1785 1786 /*@C 1787 PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian preconditioner function for given test and basis field. If this is missing, 1788 the system matrix is used to build the preconditioner. 1789 1790 Not Collective 1791 1792 Input Parameters: 1793 + ds - The `PetscDS` 1794 . f - The test field number 1795 - g - The field number 1796 1797 Output Parameters: 1798 + g0 - integrand for the test and basis function term 1799 . g1 - integrand for the test function and basis function gradient term 1800 . g2 - integrand for the test function gradient and basis function term 1801 - g3 - integrand for the test function gradient and basis function gradient term 1802 1803 Calling sequence of `g0`: 1804 + dim - the spatial dimension 1805 . Nf - the number of fields 1806 . NfAux - the number of auxiliary fields 1807 . uOff - the offset into u[] and u_t[] for each field 1808 . uOff_x - the offset into u_x[] for each field 1809 . u - each field evaluated at the current point 1810 . u_t - the time derivative of each field evaluated at the current point 1811 . u_x - the gradient of each field evaluated at the current point 1812 . aOff - the offset into a[] and a_t[] for each auxiliary field 1813 . aOff_x - the offset into a_x[] for each auxiliary field 1814 . a - each auxiliary field evaluated at the current point 1815 . a_t - the time derivative of each auxiliary field evaluated at the current point 1816 . a_x - the gradient of auxiliary each field evaluated at the current point 1817 . t - current time 1818 . u_tShift - the multiplier a for dF/dU_t 1819 . x - coordinates of the current point 1820 . numConstants - number of constant parameters 1821 . constants - constant parameters 1822 - g0 - output values at the current point 1823 1824 Level: intermediate 1825 1826 Note: 1827 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 1828 We are using a first order FEM model for the weak form\: 1829 1830 $$ 1831 \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi 1832 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi 1833 $$ 1834 1835 .seealso: `PetscDS`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()` 1836 @*/ 1837 PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1838 { 1839 PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3; 1840 PetscInt n0, n1, n2, n3; 1841 1842 PetscFunctionBegin; 1843 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1844 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 1845 PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf); 1846 PetscCall(PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3)); 1847 *g0 = tmp0 ? tmp0[0] : NULL; 1848 *g1 = tmp1 ? tmp1[0] : NULL; 1849 *g2 = tmp2 ? tmp2[0] : NULL; 1850 *g3 = tmp3 ? tmp3[0] : NULL; 1851 PetscFunctionReturn(PETSC_SUCCESS); 1852 } 1853 1854 /*@C 1855 PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian preconditioner function for given test and basis fields. 1856 If this is missing, the system matrix is used to build the preconditioner. 1857 1858 Not Collective 1859 1860 Input Parameters: 1861 + ds - The `PetscDS` 1862 . f - The test field number 1863 . g - The field number 1864 . g0 - integrand for the test and basis function term 1865 . g1 - integrand for the test function and basis function gradient term 1866 . g2 - integrand for the test function gradient and basis function term 1867 - g3 - integrand for the test function gradient and basis function gradient term 1868 1869 Calling sequence of `g0`: 1870 + dim - the spatial dimension 1871 . Nf - the number of fields 1872 . NfAux - the number of auxiliary fields 1873 . uOff - the offset into u[] and u_t[] for each field 1874 . uOff_x - the offset into u_x[] for each field 1875 . u - each field evaluated at the current point 1876 . u_t - the time derivative of each field evaluated at the current point 1877 . u_x - the gradient of each field evaluated at the current point 1878 . aOff - the offset into a[] and a_t[] for each auxiliary field 1879 . aOff_x - the offset into a_x[] for each auxiliary field 1880 . a - each auxiliary field evaluated at the current point 1881 . a_t - the time derivative of each auxiliary field evaluated at the current point 1882 . a_x - the gradient of auxiliary each field evaluated at the current point 1883 . t - current time 1884 . u_tShift - the multiplier a for dF/dU_t 1885 . x - coordinates of the current point 1886 . numConstants - number of constant parameters 1887 . constants - constant parameters 1888 - g0 - output values at the current point 1889 1890 Level: intermediate 1891 1892 Note: 1893 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 1894 1895 We are using a first order FEM model for the weak form\: 1896 1897 $$ 1898 \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi 1899 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi 1900 $$ 1901 1902 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()` 1903 @*/ 1904 PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1905 { 1906 PetscFunctionBegin; 1907 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1908 if (g0) PetscValidFunction(g0, 4); 1909 if (g1) PetscValidFunction(g1, 5); 1910 if (g2) PetscValidFunction(g2, 6); 1911 if (g3) PetscValidFunction(g3, 7); 1912 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 1913 PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g); 1914 PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3)); 1915 PetscFunctionReturn(PETSC_SUCCESS); 1916 } 1917 1918 /*@ 1919 PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set 1920 1921 Not Collective 1922 1923 Input Parameter: 1924 . ds - The `PetscDS` 1925 1926 Output Parameter: 1927 . hasDynJac - flag that pointwise function for dynamic Jacobian has been set 1928 1929 Level: intermediate 1930 1931 .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscDSGetJacobian()` 1932 @*/ 1933 PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac) 1934 { 1935 PetscFunctionBegin; 1936 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1937 PetscCall(PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac)); 1938 PetscFunctionReturn(PETSC_SUCCESS); 1939 } 1940 1941 /*@C 1942 PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field 1943 1944 Not Collective 1945 1946 Input Parameters: 1947 + ds - The `PetscDS` 1948 . f - The test field number 1949 - g - The field number 1950 1951 Output Parameters: 1952 + g0 - integrand for the test and basis function term 1953 . g1 - integrand for the test function and basis function gradient term 1954 . g2 - integrand for the test function gradient and basis function term 1955 - g3 - integrand for the test function gradient and basis function gradient term 1956 1957 Calling sequence of `g0`: 1958 + dim - the spatial dimension 1959 . Nf - the number of fields 1960 . NfAux - the number of auxiliary fields 1961 . uOff - the offset into u[] and u_t[] for each field 1962 . uOff_x - the offset into u_x[] for each field 1963 . u - each field evaluated at the current point 1964 . u_t - the time derivative of each field evaluated at the current point 1965 . u_x - the gradient of each field evaluated at the current point 1966 . aOff - the offset into a[] and a_t[] for each auxiliary field 1967 . aOff_x - the offset into a_x[] for each auxiliary field 1968 . a - each auxiliary field evaluated at the current point 1969 . a_t - the time derivative of each auxiliary field evaluated at the current point 1970 . a_x - the gradient of auxiliary each field evaluated at the current point 1971 . t - current time 1972 . u_tShift - the multiplier a for dF/dU_t 1973 . x - coordinates of the current point 1974 . numConstants - number of constant parameters 1975 . constants - constant parameters 1976 - g0 - output values at the current point 1977 1978 Level: intermediate 1979 1980 Note: 1981 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 1982 1983 We are using a first order FEM model for the weak form\: 1984 1985 $$ 1986 \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi 1987 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi 1988 $$ 1989 1990 .seealso: `PetscDS`, `PetscDSSetJacobian()` 1991 @*/ 1992 PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 1993 { 1994 PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3; 1995 PetscInt n0, n1, n2, n3; 1996 1997 PetscFunctionBegin; 1998 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1999 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 2000 PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf); 2001 PetscCall(PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3)); 2002 *g0 = tmp0 ? tmp0[0] : NULL; 2003 *g1 = tmp1 ? tmp1[0] : NULL; 2004 *g2 = tmp2 ? tmp2[0] : NULL; 2005 *g3 = tmp3 ? tmp3[0] : NULL; 2006 PetscFunctionReturn(PETSC_SUCCESS); 2007 } 2008 2009 /*@C 2010 PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields 2011 2012 Not Collective 2013 2014 Input Parameters: 2015 + ds - The `PetscDS` 2016 . f - The test field number 2017 . g - The field number 2018 . g0 - integrand for the test and basis function term 2019 . g1 - integrand for the test function and basis function gradient term 2020 . g2 - integrand for the test function gradient and basis function term 2021 - g3 - integrand for the test function gradient and basis function gradient term 2022 2023 Calling sequence of `g0`: 2024 + dim - the spatial dimension 2025 . Nf - the number of fields 2026 . NfAux - the number of auxiliary fields 2027 . uOff - the offset into u[] and u_t[] for each field 2028 . uOff_x - the offset into u_x[] for each field 2029 . u - each field evaluated at the current point 2030 . u_t - the time derivative of each field evaluated at the current point 2031 . u_x - the gradient of each field evaluated at the current point 2032 . aOff - the offset into a[] and a_t[] for each auxiliary field 2033 . aOff_x - the offset into a_x[] for each auxiliary field 2034 . a - each auxiliary field evaluated at the current point 2035 . a_t - the time derivative of each auxiliary field evaluated at the current point 2036 . a_x - the gradient of auxiliary each field evaluated at the current point 2037 . t - current time 2038 . u_tShift - the multiplier a for dF/dU_t 2039 . x - coordinates of the current point 2040 . numConstants - number of constant parameters 2041 . constants - constant parameters 2042 - g0 - output values at the current point 2043 2044 Level: intermediate 2045 2046 Note: 2047 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 2048 2049 We are using a first order FEM model for the weak form\: 2050 2051 $$ 2052 \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi 2053 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi 2054 $$ 2055 2056 .seealso: `PetscDS`, `PetscDSGetJacobian()` 2057 @*/ 2058 PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 2059 { 2060 PetscFunctionBegin; 2061 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2062 if (g0) PetscValidFunction(g0, 4); 2063 if (g1) PetscValidFunction(g1, 5); 2064 if (g2) PetscValidFunction(g2, 6); 2065 if (g3) PetscValidFunction(g3, 7); 2066 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2067 PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g); 2068 PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3)); 2069 PetscFunctionReturn(PETSC_SUCCESS); 2070 } 2071 2072 /*@C 2073 PetscDSGetRiemannSolver - Returns the Riemann solver for the given field 2074 2075 Not Collective 2076 2077 Input Parameters: 2078 + ds - The `PetscDS` object 2079 - f - The field number 2080 2081 Output Parameter: 2082 . r - Riemann solver 2083 2084 Calling sequence of `r`: 2085 + dim - The spatial dimension 2086 . Nf - The number of fields 2087 . x - The coordinates at a point on the interface 2088 . n - The normal vector to the interface 2089 . uL - The state vector to the left of the interface 2090 . uR - The state vector to the right of the interface 2091 . flux - output array of flux through the interface 2092 . numConstants - number of constant parameters 2093 . constants - constant parameters 2094 - ctx - optional user context 2095 2096 Level: intermediate 2097 2098 .seealso: `PetscDS`, `PetscDSSetRiemannSolver()` 2099 @*/ 2100 PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f, void (**r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx)) 2101 { 2102 PetscRiemannFunc *tmp; 2103 PetscInt n; 2104 2105 PetscFunctionBegin; 2106 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2107 PetscAssertPointer(r, 3); 2108 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 2109 PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp)); 2110 *r = tmp ? tmp[0] : NULL; 2111 PetscFunctionReturn(PETSC_SUCCESS); 2112 } 2113 2114 /*@C 2115 PetscDSSetRiemannSolver - Sets the Riemann solver for the given field 2116 2117 Not Collective 2118 2119 Input Parameters: 2120 + ds - The `PetscDS` object 2121 . f - The field number 2122 - r - Riemann solver 2123 2124 Calling sequence of `r`: 2125 + dim - The spatial dimension 2126 . Nf - The number of fields 2127 . x - The coordinates at a point on the interface 2128 . n - The normal vector to the interface 2129 . uL - The state vector to the left of the interface 2130 . uR - The state vector to the right of the interface 2131 . flux - output array of flux through the interface 2132 . numConstants - number of constant parameters 2133 . constants - constant parameters 2134 - ctx - optional user context 2135 2136 Level: intermediate 2137 2138 .seealso: `PetscDS`, `PetscDSGetRiemannSolver()` 2139 @*/ 2140 PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f, void (*r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx)) 2141 { 2142 PetscFunctionBegin; 2143 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2144 if (r) PetscValidFunction(r, 3); 2145 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2146 PetscCall(PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r)); 2147 PetscFunctionReturn(PETSC_SUCCESS); 2148 } 2149 2150 /*@C 2151 PetscDSGetUpdate - Get the pointwise update function for a given field 2152 2153 Not Collective 2154 2155 Input Parameters: 2156 + ds - The `PetscDS` 2157 - f - The field number 2158 2159 Output Parameter: 2160 . update - update function 2161 2162 Calling sequence of `update`: 2163 + dim - the spatial dimension 2164 . Nf - the number of fields 2165 . NfAux - the number of auxiliary fields 2166 . uOff - the offset into u[] and u_t[] for each field 2167 . uOff_x - the offset into u_x[] for each field 2168 . u - each field evaluated at the current point 2169 . u_t - the time derivative of each field evaluated at the current point 2170 . u_x - the gradient of each field evaluated at the current point 2171 . aOff - the offset into a[] and a_t[] for each auxiliary field 2172 . aOff_x - the offset into a_x[] for each auxiliary field 2173 . a - each auxiliary field evaluated at the current point 2174 . a_t - the time derivative of each auxiliary field evaluated at the current point 2175 . a_x - the gradient of auxiliary each field evaluated at the current point 2176 . t - current time 2177 . x - coordinates of the current point 2178 . numConstants - number of constant parameters 2179 . constants - constant parameters 2180 - uNew - new value for field at the current point 2181 2182 Level: intermediate 2183 2184 .seealso: `PetscDS`, `PetscDSSetUpdate()`, `PetscDSSetResidual()` 2185 @*/ 2186 PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f, void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[])) 2187 { 2188 PetscFunctionBegin; 2189 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2190 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 2191 if (update) { 2192 PetscAssertPointer(update, 3); 2193 *update = ds->update[f]; 2194 } 2195 PetscFunctionReturn(PETSC_SUCCESS); 2196 } 2197 2198 /*@C 2199 PetscDSSetUpdate - Set the pointwise update function for a given field 2200 2201 Not Collective 2202 2203 Input Parameters: 2204 + ds - The `PetscDS` 2205 . f - The field number 2206 - update - update function 2207 2208 Calling sequence of `update`: 2209 + dim - the spatial dimension 2210 . Nf - the number of fields 2211 . NfAux - the number of auxiliary fields 2212 . uOff - the offset into u[] and u_t[] for each field 2213 . uOff_x - the offset into u_x[] for each field 2214 . u - each field evaluated at the current point 2215 . u_t - the time derivative of each field evaluated at the current point 2216 . u_x - the gradient of each field evaluated at the current point 2217 . aOff - the offset into a[] and a_t[] for each auxiliary field 2218 . aOff_x - the offset into a_x[] for each auxiliary field 2219 . a - each auxiliary field evaluated at the current point 2220 . a_t - the time derivative of each auxiliary field evaluated at the current point 2221 . a_x - the gradient of auxiliary each field evaluated at the current point 2222 . t - current time 2223 . x - coordinates of the current point 2224 . numConstants - number of constant parameters 2225 . constants - constant parameters 2226 - uNew - new field values at the current point 2227 2228 Level: intermediate 2229 2230 .seealso: `PetscDS`, `PetscDSGetResidual()` 2231 @*/ 2232 PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f, void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[])) 2233 { 2234 PetscFunctionBegin; 2235 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2236 if (update) PetscValidFunction(update, 3); 2237 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2238 PetscCall(PetscDSEnlarge_Static(ds, f + 1)); 2239 ds->update[f] = update; 2240 PetscFunctionReturn(PETSC_SUCCESS); 2241 } 2242 2243 PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void *ctx) 2244 { 2245 PetscFunctionBegin; 2246 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2247 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 2248 PetscAssertPointer(ctx, 3); 2249 *(void **)ctx = ds->ctx[f]; 2250 PetscFunctionReturn(PETSC_SUCCESS); 2251 } 2252 2253 PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx) 2254 { 2255 PetscFunctionBegin; 2256 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2257 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2258 PetscCall(PetscDSEnlarge_Static(ds, f + 1)); 2259 ds->ctx[f] = ctx; 2260 PetscFunctionReturn(PETSC_SUCCESS); 2261 } 2262 2263 /*@C 2264 PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field 2265 2266 Not Collective 2267 2268 Input Parameters: 2269 + ds - The PetscDS 2270 - f - The test field number 2271 2272 Output Parameters: 2273 + f0 - boundary integrand for the test function term 2274 - f1 - boundary integrand for the test function gradient term 2275 2276 Calling sequence of `f0`: 2277 + dim - the spatial dimension 2278 . Nf - the number of fields 2279 . NfAux - the number of auxiliary fields 2280 . uOff - the offset into u[] and u_t[] for each field 2281 . uOff_x - the offset into u_x[] for each field 2282 . u - each field evaluated at the current point 2283 . u_t - the time derivative of each field evaluated at the current point 2284 . u_x - the gradient of each field evaluated at the current point 2285 . aOff - the offset into a[] and a_t[] for each auxiliary field 2286 . aOff_x - the offset into a_x[] for each auxiliary field 2287 . a - each auxiliary field evaluated at the current point 2288 . a_t - the time derivative of each auxiliary field evaluated at the current point 2289 . a_x - the gradient of auxiliary each field evaluated at the current point 2290 . t - current time 2291 . x - coordinates of the current point 2292 . n - unit normal at the current point 2293 . numConstants - number of constant parameters 2294 . constants - constant parameters 2295 - f0 - output values at the current point 2296 2297 Level: intermediate 2298 2299 Note: 2300 The calling sequence of `f1` is identical, and therefore omitted for brevity. 2301 2302 We are using a first order FEM model for the weak form\: 2303 2304 $$ 2305 \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n 2306 $$ 2307 2308 .seealso: `PetscDS`, `PetscDSSetBdResidual()` 2309 @*/ 2310 PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 2311 { 2312 PetscBdPointFunc *tmp0, *tmp1; 2313 PetscInt n0, n1; 2314 2315 PetscFunctionBegin; 2316 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2317 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 2318 PetscCall(PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1)); 2319 *f0 = tmp0 ? tmp0[0] : NULL; 2320 *f1 = tmp1 ? tmp1[0] : NULL; 2321 PetscFunctionReturn(PETSC_SUCCESS); 2322 } 2323 2324 /*@C 2325 PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field 2326 2327 Not Collective 2328 2329 Input Parameters: 2330 + ds - The `PetscDS` 2331 . f - The test field number 2332 . f0 - boundary integrand for the test function term 2333 - f1 - boundary integrand for the test function gradient term 2334 2335 Calling sequence of `f0`: 2336 + dim - the spatial dimension 2337 . Nf - the number of fields 2338 . NfAux - the number of auxiliary fields 2339 . uOff - the offset into u[] and u_t[] for each field 2340 . uOff_x - the offset into u_x[] for each field 2341 . u - each field evaluated at the current point 2342 . u_t - the time derivative of each field evaluated at the current point 2343 . u_x - the gradient of each field evaluated at the current point 2344 . aOff - the offset into a[] and a_t[] for each auxiliary field 2345 . aOff_x - the offset into a_x[] for each auxiliary field 2346 . a - each auxiliary field evaluated at the current point 2347 . a_t - the time derivative of each auxiliary field evaluated at the current point 2348 . a_x - the gradient of auxiliary each field evaluated at the current point 2349 . t - current time 2350 . x - coordinates of the current point 2351 . n - unit normal at the current point 2352 . numConstants - number of constant parameters 2353 . constants - constant parameters 2354 - f0 - output values at the current point 2355 2356 Level: intermediate 2357 2358 Note: 2359 The calling sequence of `f1` is identical, and therefore omitted for brevity. 2360 2361 We are using a first order FEM model for the weak form\: 2362 2363 $$ 2364 \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n 2365 $$ 2366 2367 .seealso: `PetscDS`, `PetscDSGetBdResidual()` 2368 @*/ 2369 PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 2370 { 2371 PetscFunctionBegin; 2372 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2373 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2374 PetscCall(PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1)); 2375 PetscFunctionReturn(PETSC_SUCCESS); 2376 } 2377 2378 /*@ 2379 PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set 2380 2381 Not Collective 2382 2383 Input Parameter: 2384 . ds - The `PetscDS` 2385 2386 Output Parameter: 2387 . hasBdJac - flag that pointwise function for the boundary Jacobian has been set 2388 2389 Level: intermediate 2390 2391 .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()` 2392 @*/ 2393 PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac) 2394 { 2395 PetscFunctionBegin; 2396 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2397 PetscAssertPointer(hasBdJac, 2); 2398 PetscCall(PetscWeakFormHasBdJacobian(ds->wf, hasBdJac)); 2399 PetscFunctionReturn(PETSC_SUCCESS); 2400 } 2401 2402 /*@C 2403 PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field 2404 2405 Not Collective 2406 2407 Input Parameters: 2408 + ds - The `PetscDS` 2409 . f - The test field number 2410 - g - The field number 2411 2412 Output Parameters: 2413 + g0 - integrand for the test and basis function term 2414 . g1 - integrand for the test function and basis function gradient term 2415 . g2 - integrand for the test function gradient and basis function term 2416 - g3 - integrand for the test function gradient and basis function gradient term 2417 2418 Calling sequence of `g0`: 2419 + dim - the spatial dimension 2420 . Nf - the number of fields 2421 . NfAux - the number of auxiliary fields 2422 . uOff - the offset into u[] and u_t[] for each field 2423 . uOff_x - the offset into u_x[] for each field 2424 . u - each field evaluated at the current point 2425 . u_t - the time derivative of each field evaluated at the current point 2426 . u_x - the gradient of each field evaluated at the current point 2427 . aOff - the offset into a[] and a_t[] for each auxiliary field 2428 . aOff_x - the offset into a_x[] for each auxiliary field 2429 . a - each auxiliary field evaluated at the current point 2430 . a_t - the time derivative of each auxiliary field evaluated at the current point 2431 . a_x - the gradient of auxiliary each field evaluated at the current point 2432 . t - current time 2433 . u_tShift - the multiplier a for dF/dU_t 2434 . x - coordinates of the current point 2435 . n - normal at the current point 2436 . numConstants - number of constant parameters 2437 . constants - constant parameters 2438 - g0 - output values at the current point 2439 2440 Level: intermediate 2441 2442 Note: 2443 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 2444 2445 We are using a first order FEM model for the weak form\: 2446 2447 $$ 2448 \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi 2449 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi 2450 $$ 2451 2452 .seealso: `PetscDS`, `PetscDSSetBdJacobian()` 2453 @*/ 2454 PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 2455 { 2456 PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3; 2457 PetscInt n0, n1, n2, n3; 2458 2459 PetscFunctionBegin; 2460 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2461 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 2462 PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf); 2463 PetscCall(PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3)); 2464 *g0 = tmp0 ? tmp0[0] : NULL; 2465 *g1 = tmp1 ? tmp1[0] : NULL; 2466 *g2 = tmp2 ? tmp2[0] : NULL; 2467 *g3 = tmp3 ? tmp3[0] : NULL; 2468 PetscFunctionReturn(PETSC_SUCCESS); 2469 } 2470 2471 /*@C 2472 PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field 2473 2474 Not Collective 2475 2476 Input Parameters: 2477 + ds - The PetscDS 2478 . f - The test field number 2479 . g - The field number 2480 . g0 - integrand for the test and basis function term 2481 . g1 - integrand for the test function and basis function gradient term 2482 . g2 - integrand for the test function gradient and basis function term 2483 - g3 - integrand for the test function gradient and basis function gradient term 2484 2485 Calling sequence of `g0`: 2486 + dim - the spatial dimension 2487 . Nf - the number of fields 2488 . NfAux - the number of auxiliary fields 2489 . uOff - the offset into u[] and u_t[] for each field 2490 . uOff_x - the offset into u_x[] for each field 2491 . u - each field evaluated at the current point 2492 . u_t - the time derivative of each field evaluated at the current point 2493 . u_x - the gradient of each field evaluated at the current point 2494 . aOff - the offset into a[] and a_t[] for each auxiliary field 2495 . aOff_x - the offset into a_x[] for each auxiliary field 2496 . a - each auxiliary field evaluated at the current point 2497 . a_t - the time derivative of each auxiliary field evaluated at the current point 2498 . a_x - the gradient of auxiliary each field evaluated at the current point 2499 . t - current time 2500 . u_tShift - the multiplier a for dF/dU_t 2501 . x - coordinates of the current point 2502 . n - normal at the current point 2503 . numConstants - number of constant parameters 2504 . constants - constant parameters 2505 - g0 - output values at the current point 2506 2507 Level: intermediate 2508 2509 Note: 2510 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 2511 2512 We are using a first order FEM model for the weak form\: 2513 2514 $$ 2515 \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi 2516 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi 2517 $$ 2518 2519 .seealso: `PetscDS`, `PetscDSGetBdJacobian()` 2520 @*/ 2521 PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 2522 { 2523 PetscFunctionBegin; 2524 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2525 if (g0) PetscValidFunction(g0, 4); 2526 if (g1) PetscValidFunction(g1, 5); 2527 if (g2) PetscValidFunction(g2, 6); 2528 if (g3) PetscValidFunction(g3, 7); 2529 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2530 PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g); 2531 PetscCall(PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3)); 2532 PetscFunctionReturn(PETSC_SUCCESS); 2533 } 2534 2535 /*@ 2536 PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set 2537 2538 Not Collective 2539 2540 Input Parameter: 2541 . ds - The `PetscDS` 2542 2543 Output Parameter: 2544 . hasBdJacPre - flag that pointwise function for the boundary Jacobian preconditioner has been set 2545 2546 Level: intermediate 2547 2548 .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()` 2549 @*/ 2550 PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre) 2551 { 2552 PetscFunctionBegin; 2553 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2554 PetscAssertPointer(hasBdJacPre, 2); 2555 PetscCall(PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre)); 2556 PetscFunctionReturn(PETSC_SUCCESS); 2557 } 2558 2559 /*@C 2560 PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian preconditioner function for given test and basis field 2561 2562 Not Collective; No Fortran Support 2563 2564 Input Parameters: 2565 + ds - The `PetscDS` 2566 . f - The test field number 2567 - g - The field number 2568 2569 Output Parameters: 2570 + g0 - integrand for the test and basis function term 2571 . g1 - integrand for the test function and basis function gradient term 2572 . g2 - integrand for the test function gradient and basis function term 2573 - g3 - integrand for the test function gradient and basis function gradient term 2574 2575 Calling sequence of `g0`: 2576 + dim - the spatial dimension 2577 . Nf - the number of fields 2578 . NfAux - the number of auxiliary fields 2579 . uOff - the offset into u[] and u_t[] for each field 2580 . uOff_x - the offset into u_x[] for each field 2581 . u - each field evaluated at the current point 2582 . u_t - the time derivative of each field evaluated at the current point 2583 . u_x - the gradient of each field evaluated at the current point 2584 . aOff - the offset into a[] and a_t[] for each auxiliary field 2585 . aOff_x - the offset into a_x[] for each auxiliary field 2586 . a - each auxiliary field evaluated at the current point 2587 . a_t - the time derivative of each auxiliary field evaluated at the current point 2588 . a_x - the gradient of auxiliary each field evaluated at the current point 2589 . t - current time 2590 . u_tShift - the multiplier a for dF/dU_t 2591 . x - coordinates of the current point 2592 . n - normal at the current point 2593 . numConstants - number of constant parameters 2594 . constants - constant parameters 2595 - g0 - output values at the current point 2596 2597 Level: intermediate 2598 2599 Note: 2600 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 2601 2602 We are using a first order FEM model for the weak form\: 2603 2604 $$ 2605 \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi 2606 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi 2607 $$ 2608 2609 .seealso: `PetscDS`, `PetscDSSetBdJacobianPreconditioner()` 2610 @*/ 2611 PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 2612 { 2613 PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3; 2614 PetscInt n0, n1, n2, n3; 2615 2616 PetscFunctionBegin; 2617 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2618 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 2619 PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf); 2620 PetscCall(PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3)); 2621 *g0 = tmp0 ? tmp0[0] : NULL; 2622 *g1 = tmp1 ? tmp1[0] : NULL; 2623 *g2 = tmp2 ? tmp2[0] : NULL; 2624 *g3 = tmp3 ? tmp3[0] : NULL; 2625 PetscFunctionReturn(PETSC_SUCCESS); 2626 } 2627 2628 /*@C 2629 PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field 2630 2631 Not Collective; No Fortran Support 2632 2633 Input Parameters: 2634 + ds - The `PetscDS` 2635 . f - The test field number 2636 . g - The field number 2637 . g0 - integrand for the test and basis function term 2638 . g1 - integrand for the test function and basis function gradient term 2639 . g2 - integrand for the test function gradient and basis function term 2640 - g3 - integrand for the test function gradient and basis function gradient term 2641 2642 Calling sequence of `g0': 2643 + dim - the spatial dimension 2644 . Nf - the number of fields 2645 . NfAux - the number of auxiliary fields 2646 . uOff - the offset into u[] and u_t[] for each field 2647 . uOff_x - the offset into u_x[] for each field 2648 . u - each field evaluated at the current point 2649 . u_t - the time derivative of each field evaluated at the current point 2650 . u_x - the gradient of each field evaluated at the current point 2651 . aOff - the offset into a[] and a_t[] for each auxiliary field 2652 . aOff_x - the offset into a_x[] for each auxiliary field 2653 . a - each auxiliary field evaluated at the current point 2654 . a_t - the time derivative of each auxiliary field evaluated at the current point 2655 . a_x - the gradient of auxiliary each field evaluated at the current point 2656 . t - current time 2657 . u_tShift - the multiplier a for dF/dU_t 2658 . x - coordinates of the current point 2659 . n - normal at the current point 2660 . numConstants - number of constant parameters 2661 . constants - constant parameters 2662 - g0 - output values at the current point 2663 2664 Level: intermediate 2665 2666 Note: 2667 `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity. 2668 2669 We are using a first order FEM model for the weak form\: 2670 2671 $$ 2672 \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi 2673 + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi 2674 $$ 2675 2676 .seealso: `PetscDS`, `PetscDSGetBdJacobianPreconditioner()` 2677 @*/ 2678 PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) 2679 { 2680 PetscFunctionBegin; 2681 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2682 if (g0) PetscValidFunction(g0, 4); 2683 if (g1) PetscValidFunction(g1, 5); 2684 if (g2) PetscValidFunction(g2, 6); 2685 if (g3) PetscValidFunction(g3, 7); 2686 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2687 PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g); 2688 PetscCall(PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3)); 2689 PetscFunctionReturn(PETSC_SUCCESS); 2690 } 2691 2692 /*@C 2693 PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field 2694 2695 Not Collective 2696 2697 Input Parameters: 2698 + prob - The PetscDS 2699 - f - The test field number 2700 2701 Output Parameters: 2702 + sol - exact solution for the test field 2703 - ctx - exact solution context 2704 2705 Calling sequence of `exactSol`: 2706 + dim - the spatial dimension 2707 . t - current time 2708 . x - coordinates of the current point 2709 . Nc - the number of field components 2710 . u - the solution field evaluated at the current point 2711 - ctx - a user context 2712 2713 Level: intermediate 2714 2715 .seealso: `PetscDS`, `PetscDSSetExactSolution()`, `PetscDSGetExactSolutionTimeDerivative()` 2716 @*/ 2717 PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx) 2718 { 2719 PetscFunctionBegin; 2720 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 2721 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 2722 if (sol) { 2723 PetscAssertPointer(sol, 3); 2724 *sol = prob->exactSol[f]; 2725 } 2726 if (ctx) { 2727 PetscAssertPointer(ctx, 4); 2728 *ctx = prob->exactCtx[f]; 2729 } 2730 PetscFunctionReturn(PETSC_SUCCESS); 2731 } 2732 2733 /*@C 2734 PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field 2735 2736 Not Collective 2737 2738 Input Parameters: 2739 + prob - The `PetscDS` 2740 . f - The test field number 2741 . sol - solution function for the test fields 2742 - ctx - solution context or `NULL` 2743 2744 Calling sequence of `sol`: 2745 + dim - the spatial dimension 2746 . t - current time 2747 . x - coordinates of the current point 2748 . Nc - the number of field components 2749 . u - the solution field evaluated at the current point 2750 - ctx - a user context 2751 2752 Level: intermediate 2753 2754 .seealso: `PetscDS`, `PetscDSGetExactSolution()` 2755 @*/ 2756 PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx) 2757 { 2758 PetscFunctionBegin; 2759 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 2760 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2761 PetscCall(PetscDSEnlarge_Static(prob, f + 1)); 2762 if (sol) { 2763 PetscValidFunction(sol, 3); 2764 prob->exactSol[f] = sol; 2765 } 2766 if (ctx) { 2767 PetscValidFunction(ctx, 4); 2768 prob->exactCtx[f] = ctx; 2769 } 2770 PetscFunctionReturn(PETSC_SUCCESS); 2771 } 2772 2773 /*@C 2774 PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field 2775 2776 Not Collective 2777 2778 Input Parameters: 2779 + prob - The `PetscDS` 2780 - f - The test field number 2781 2782 Output Parameters: 2783 + sol - time derivative of the exact solution for the test field 2784 - ctx - time derivative of the exact solution context 2785 2786 Calling sequence of `exactSol`: 2787 + dim - the spatial dimension 2788 . t - current time 2789 . x - coordinates of the current point 2790 . Nc - the number of field components 2791 . u - the solution field evaluated at the current point 2792 - ctx - a user context 2793 2794 Level: intermediate 2795 2796 .seealso: `PetscDS`, `PetscDSSetExactSolutionTimeDerivative()`, `PetscDSGetExactSolution()` 2797 @*/ 2798 PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx) 2799 { 2800 PetscFunctionBegin; 2801 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 2802 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 2803 if (sol) { 2804 PetscAssertPointer(sol, 3); 2805 *sol = prob->exactSol_t[f]; 2806 } 2807 if (ctx) { 2808 PetscAssertPointer(ctx, 4); 2809 *ctx = prob->exactCtx_t[f]; 2810 } 2811 PetscFunctionReturn(PETSC_SUCCESS); 2812 } 2813 2814 /*@C 2815 PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field 2816 2817 Not Collective 2818 2819 Input Parameters: 2820 + prob - The `PetscDS` 2821 . f - The test field number 2822 . sol - time derivative of the solution function for the test fields 2823 - ctx - time derivative of the solution context or `NULL` 2824 2825 Calling sequence of `sol`: 2826 + dim - the spatial dimension 2827 . t - current time 2828 . x - coordinates of the current point 2829 . Nc - the number of field components 2830 . u - the solution field evaluated at the current point 2831 - ctx - a user context 2832 2833 Level: intermediate 2834 2835 .seealso: `PetscDS`, `PetscDSGetExactSolutionTimeDerivative()`, `PetscDSSetExactSolution()` 2836 @*/ 2837 PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx) 2838 { 2839 PetscFunctionBegin; 2840 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 2841 PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f); 2842 PetscCall(PetscDSEnlarge_Static(prob, f + 1)); 2843 if (sol) { 2844 PetscValidFunction(sol, 3); 2845 prob->exactSol_t[f] = sol; 2846 } 2847 if (ctx) { 2848 PetscValidFunction(ctx, 4); 2849 prob->exactCtx_t[f] = ctx; 2850 } 2851 PetscFunctionReturn(PETSC_SUCCESS); 2852 } 2853 2854 /*@C 2855 PetscDSGetConstants - Returns the array of constants passed to point functions 2856 2857 Not Collective 2858 2859 Input Parameter: 2860 . ds - The `PetscDS` object 2861 2862 Output Parameters: 2863 + numConstants - The number of constants 2864 - constants - The array of constants, NULL if there are none 2865 2866 Level: intermediate 2867 2868 .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSCreate()` 2869 @*/ 2870 PetscErrorCode PetscDSGetConstants(PetscDS ds, PetscInt *numConstants, const PetscScalar *constants[]) 2871 { 2872 PetscFunctionBegin; 2873 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2874 if (numConstants) { 2875 PetscAssertPointer(numConstants, 2); 2876 *numConstants = ds->numConstants; 2877 } 2878 if (constants) { 2879 PetscAssertPointer(constants, 3); 2880 *constants = ds->constants; 2881 } 2882 PetscFunctionReturn(PETSC_SUCCESS); 2883 } 2884 2885 /*@C 2886 PetscDSSetConstants - Set the array of constants passed to point functions 2887 2888 Not Collective 2889 2890 Input Parameters: 2891 + ds - The `PetscDS` object 2892 . numConstants - The number of constants 2893 - constants - The array of constants, `NULL` if there are none 2894 2895 Level: intermediate 2896 2897 .seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()` 2898 @*/ 2899 PetscErrorCode PetscDSSetConstants(PetscDS ds, PetscInt numConstants, PetscScalar constants[]) 2900 { 2901 PetscFunctionBegin; 2902 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2903 if (numConstants != ds->numConstants) { 2904 PetscCall(PetscFree(ds->constants)); 2905 ds->numConstants = numConstants; 2906 PetscCall(PetscMalloc1(ds->numConstants + ds->numFuncConstants, &ds->constants)); 2907 } 2908 if (ds->numConstants) { 2909 PetscAssertPointer(constants, 3); 2910 PetscCall(PetscArraycpy(ds->constants, constants, ds->numConstants)); 2911 } 2912 PetscFunctionReturn(PETSC_SUCCESS); 2913 } 2914 2915 /*@C 2916 PetscDSSetIntegrationParameters - Set the parameters for a particular integration 2917 2918 Not Collective 2919 2920 Input Parameters: 2921 + ds - The `PetscDS` object 2922 . fieldI - The test field for a given point function, or PETSC_DETERMINE 2923 - fieldJ - The basis field for a given point function, or PETSC_DETERMINE 2924 2925 Level: intermediate 2926 2927 .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()` 2928 @*/ 2929 PetscErrorCode PetscDSSetIntegrationParameters(PetscDS ds, PetscInt fieldI, PetscInt fieldJ) 2930 { 2931 PetscFunctionBegin; 2932 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2933 ds->constants[ds->numConstants] = fieldI; 2934 ds->constants[ds->numConstants + 1] = fieldJ; 2935 PetscFunctionReturn(PETSC_SUCCESS); 2936 } 2937 2938 /*@C 2939 PetscDSSetCellParameters - Set the parameters for a particular cell 2940 2941 Not Collective 2942 2943 Input Parameters: 2944 + ds - The `PetscDS` object 2945 - volume - The cell volume 2946 2947 Level: intermediate 2948 2949 .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()` 2950 @*/ 2951 PetscErrorCode PetscDSSetCellParameters(PetscDS ds, PetscReal volume) 2952 { 2953 PetscFunctionBegin; 2954 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 2955 ds->constants[ds->numConstants + 2] = volume; 2956 PetscFunctionReturn(PETSC_SUCCESS); 2957 } 2958 2959 /*@ 2960 PetscDSGetFieldIndex - Returns the index of the given field 2961 2962 Not Collective 2963 2964 Input Parameters: 2965 + prob - The `PetscDS` object 2966 - disc - The discretization object 2967 2968 Output Parameter: 2969 . f - The field number 2970 2971 Level: beginner 2972 2973 .seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 2974 @*/ 2975 PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f) 2976 { 2977 PetscInt g; 2978 2979 PetscFunctionBegin; 2980 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 2981 PetscAssertPointer(f, 3); 2982 *f = -1; 2983 for (g = 0; g < prob->Nf; ++g) { 2984 if (disc == prob->disc[g]) break; 2985 } 2986 PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS."); 2987 *f = g; 2988 PetscFunctionReturn(PETSC_SUCCESS); 2989 } 2990 2991 /*@ 2992 PetscDSGetFieldSize - Returns the size of the given field in the full space basis 2993 2994 Not Collective 2995 2996 Input Parameters: 2997 + prob - The `PetscDS` object 2998 - f - The field number 2999 3000 Output Parameter: 3001 . size - The size 3002 3003 Level: beginner 3004 3005 .seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3006 @*/ 3007 PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size) 3008 { 3009 PetscFunctionBegin; 3010 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3011 PetscAssertPointer(size, 3); 3012 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 3013 PetscCall(PetscDSSetUp(prob)); 3014 *size = prob->Nb[f]; 3015 PetscFunctionReturn(PETSC_SUCCESS); 3016 } 3017 3018 /*@ 3019 PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis 3020 3021 Not Collective 3022 3023 Input Parameters: 3024 + prob - The `PetscDS` object 3025 - f - The field number 3026 3027 Output Parameter: 3028 . off - The offset 3029 3030 Level: beginner 3031 3032 .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3033 @*/ 3034 PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off) 3035 { 3036 PetscInt size, g; 3037 3038 PetscFunctionBegin; 3039 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3040 PetscAssertPointer(off, 3); 3041 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 3042 *off = 0; 3043 for (g = 0; g < f; ++g) { 3044 PetscCall(PetscDSGetFieldSize(prob, g, &size)); 3045 *off += size; 3046 } 3047 PetscFunctionReturn(PETSC_SUCCESS); 3048 } 3049 3050 /*@ 3051 PetscDSGetFieldOffsetCohesive - Returns the offset of the given field in the full space basis on a cohesive cell 3052 3053 Not Collective 3054 3055 Input Parameters: 3056 + ds - The `PetscDS` object 3057 - f - The field number 3058 3059 Output Parameter: 3060 . off - The offset 3061 3062 Level: beginner 3063 3064 .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3065 @*/ 3066 PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off) 3067 { 3068 PetscInt size, g; 3069 3070 PetscFunctionBegin; 3071 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3072 PetscAssertPointer(off, 3); 3073 PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf); 3074 *off = 0; 3075 for (g = 0; g < f; ++g) { 3076 PetscBool cohesive; 3077 3078 PetscCall(PetscDSGetCohesive(ds, g, &cohesive)); 3079 PetscCall(PetscDSGetFieldSize(ds, g, &size)); 3080 *off += cohesive ? size : size * 2; 3081 } 3082 PetscFunctionReturn(PETSC_SUCCESS); 3083 } 3084 3085 /*@ 3086 PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point 3087 3088 Not Collective 3089 3090 Input Parameter: 3091 . prob - The `PetscDS` object 3092 3093 Output Parameter: 3094 . dimensions - The number of dimensions 3095 3096 Level: beginner 3097 3098 .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3099 @*/ 3100 PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[]) 3101 { 3102 PetscFunctionBegin; 3103 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3104 PetscCall(PetscDSSetUp(prob)); 3105 PetscAssertPointer(dimensions, 2); 3106 *dimensions = prob->Nb; 3107 PetscFunctionReturn(PETSC_SUCCESS); 3108 } 3109 3110 /*@ 3111 PetscDSGetComponents - Returns the number of components for each field on an evaluation point 3112 3113 Not Collective 3114 3115 Input Parameter: 3116 . prob - The `PetscDS` object 3117 3118 Output Parameter: 3119 . components - The number of components 3120 3121 Level: beginner 3122 3123 .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3124 @*/ 3125 PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[]) 3126 { 3127 PetscFunctionBegin; 3128 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3129 PetscCall(PetscDSSetUp(prob)); 3130 PetscAssertPointer(components, 2); 3131 *components = prob->Nc; 3132 PetscFunctionReturn(PETSC_SUCCESS); 3133 } 3134 3135 /*@ 3136 PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point 3137 3138 Not Collective 3139 3140 Input Parameters: 3141 + prob - The `PetscDS` object 3142 - f - The field number 3143 3144 Output Parameter: 3145 . off - The offset 3146 3147 Level: beginner 3148 3149 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3150 @*/ 3151 PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off) 3152 { 3153 PetscFunctionBegin; 3154 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3155 PetscAssertPointer(off, 3); 3156 PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf); 3157 PetscCall(PetscDSSetUp(prob)); 3158 *off = prob->off[f]; 3159 PetscFunctionReturn(PETSC_SUCCESS); 3160 } 3161 3162 /*@ 3163 PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point 3164 3165 Not Collective 3166 3167 Input Parameter: 3168 . prob - The `PetscDS` object 3169 3170 Output Parameter: 3171 . offsets - The offsets 3172 3173 Level: beginner 3174 3175 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3176 @*/ 3177 PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[]) 3178 { 3179 PetscFunctionBegin; 3180 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3181 PetscAssertPointer(offsets, 2); 3182 PetscCall(PetscDSSetUp(prob)); 3183 *offsets = prob->off; 3184 PetscFunctionReturn(PETSC_SUCCESS); 3185 } 3186 3187 /*@ 3188 PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point 3189 3190 Not Collective 3191 3192 Input Parameter: 3193 . prob - The `PetscDS` object 3194 3195 Output Parameter: 3196 . offsets - The offsets 3197 3198 Level: beginner 3199 3200 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3201 @*/ 3202 PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[]) 3203 { 3204 PetscFunctionBegin; 3205 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3206 PetscAssertPointer(offsets, 2); 3207 PetscCall(PetscDSSetUp(prob)); 3208 *offsets = prob->offDer; 3209 PetscFunctionReturn(PETSC_SUCCESS); 3210 } 3211 3212 /*@ 3213 PetscDSGetComponentOffsetsCohesive - Returns the offset of each field on an evaluation point 3214 3215 Not Collective 3216 3217 Input Parameters: 3218 + ds - The `PetscDS` object 3219 - s - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive 3220 3221 Output Parameter: 3222 . offsets - The offsets 3223 3224 Level: beginner 3225 3226 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3227 @*/ 3228 PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[]) 3229 { 3230 PetscFunctionBegin; 3231 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3232 PetscAssertPointer(offsets, 3); 3233 PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS"); 3234 PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s); 3235 PetscCall(PetscDSSetUp(ds)); 3236 *offsets = ds->offCohesive[s]; 3237 PetscFunctionReturn(PETSC_SUCCESS); 3238 } 3239 3240 /*@ 3241 PetscDSGetComponentDerivativeOffsetsCohesive - Returns the offset of each field derivative on an evaluation point 3242 3243 Not Collective 3244 3245 Input Parameters: 3246 + ds - The `PetscDS` object 3247 - s - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive 3248 3249 Output Parameter: 3250 . offsets - The offsets 3251 3252 Level: beginner 3253 3254 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()` 3255 @*/ 3256 PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[]) 3257 { 3258 PetscFunctionBegin; 3259 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3260 PetscAssertPointer(offsets, 3); 3261 PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS"); 3262 PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s); 3263 PetscCall(PetscDSSetUp(ds)); 3264 *offsets = ds->offDerCohesive[s]; 3265 PetscFunctionReturn(PETSC_SUCCESS); 3266 } 3267 3268 /*@C 3269 PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization 3270 3271 Not Collective 3272 3273 Input Parameter: 3274 . prob - The `PetscDS` object 3275 3276 Output Parameter: 3277 . T - The basis function and derivatives tabulation at quadrature points for each field 3278 3279 Level: intermediate 3280 3281 .seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()` 3282 @*/ 3283 PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[]) 3284 { 3285 PetscFunctionBegin; 3286 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3287 PetscAssertPointer(T, 2); 3288 PetscCall(PetscDSSetUp(prob)); 3289 *T = prob->T; 3290 PetscFunctionReturn(PETSC_SUCCESS); 3291 } 3292 3293 /*@C 3294 PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces 3295 3296 Not Collective 3297 3298 Input Parameter: 3299 . prob - The `PetscDS` object 3300 3301 Output Parameter: 3302 . Tf - The basis function and derivative tabulation on each local face at quadrature points for each and field 3303 3304 Level: intermediate 3305 3306 .seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()` 3307 @*/ 3308 PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[]) 3309 { 3310 PetscFunctionBegin; 3311 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3312 PetscAssertPointer(Tf, 2); 3313 PetscCall(PetscDSSetUp(prob)); 3314 *Tf = prob->Tf; 3315 PetscFunctionReturn(PETSC_SUCCESS); 3316 } 3317 3318 PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x) 3319 { 3320 PetscFunctionBegin; 3321 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3322 PetscCall(PetscDSSetUp(prob)); 3323 if (u) { 3324 PetscAssertPointer(u, 2); 3325 *u = prob->u; 3326 } 3327 if (u_t) { 3328 PetscAssertPointer(u_t, 3); 3329 *u_t = prob->u_t; 3330 } 3331 if (u_x) { 3332 PetscAssertPointer(u_x, 4); 3333 *u_x = prob->u_x; 3334 } 3335 PetscFunctionReturn(PETSC_SUCCESS); 3336 } 3337 3338 PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3) 3339 { 3340 PetscFunctionBegin; 3341 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3342 PetscCall(PetscDSSetUp(prob)); 3343 if (f0) { 3344 PetscAssertPointer(f0, 2); 3345 *f0 = prob->f0; 3346 } 3347 if (f1) { 3348 PetscAssertPointer(f1, 3); 3349 *f1 = prob->f1; 3350 } 3351 if (g0) { 3352 PetscAssertPointer(g0, 4); 3353 *g0 = prob->g0; 3354 } 3355 if (g1) { 3356 PetscAssertPointer(g1, 5); 3357 *g1 = prob->g1; 3358 } 3359 if (g2) { 3360 PetscAssertPointer(g2, 6); 3361 *g2 = prob->g2; 3362 } 3363 if (g3) { 3364 PetscAssertPointer(g3, 7); 3365 *g3 = prob->g3; 3366 } 3367 PetscFunctionReturn(PETSC_SUCCESS); 3368 } 3369 3370 PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal) 3371 { 3372 PetscFunctionBegin; 3373 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3374 PetscCall(PetscDSSetUp(prob)); 3375 if (x) { 3376 PetscAssertPointer(x, 2); 3377 *x = prob->x; 3378 } 3379 if (basisReal) { 3380 PetscAssertPointer(basisReal, 3); 3381 *basisReal = prob->basisReal; 3382 } 3383 if (basisDerReal) { 3384 PetscAssertPointer(basisDerReal, 4); 3385 *basisDerReal = prob->basisDerReal; 3386 } 3387 if (testReal) { 3388 PetscAssertPointer(testReal, 5); 3389 *testReal = prob->testReal; 3390 } 3391 if (testDerReal) { 3392 PetscAssertPointer(testDerReal, 6); 3393 *testDerReal = prob->testDerReal; 3394 } 3395 PetscFunctionReturn(PETSC_SUCCESS); 3396 } 3397 3398 /*@C 3399 PetscDSAddBoundary - Add a boundary condition to the model. 3400 3401 Collective 3402 3403 Input Parameters: 3404 + ds - The PetscDS object 3405 . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann) 3406 . name - The BC name 3407 . label - The label defining constrained points 3408 . Nv - The number of `DMLabel` values for constrained points 3409 . values - An array of label values for constrained points 3410 . field - The field to constrain 3411 . Nc - The number of constrained field components (0 will constrain all fields) 3412 . comps - An array of constrained component numbers 3413 . bcFunc - A pointwise function giving boundary values 3414 . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL 3415 - ctx - An optional user context for bcFunc 3416 3417 Output Parameter: 3418 . bd - The boundary number 3419 3420 Options Database Keys: 3421 + -bc_<boundary name> <num> - Overrides the boundary ids 3422 - -bc_<boundary name>_comp <num> - Overrides the boundary components 3423 3424 Level: developer 3425 3426 Note: 3427 Both `bcFunc` and `bcFunc_t` will depend on the boundary condition type. If the type if `DM_BC_ESSENTIAL`, then the calling sequence is\: 3428 3429 $ void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[]) 3430 3431 If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is\: 3432 .vb 3433 void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux, 3434 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 3435 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 3436 PetscReal time, const PetscReal x[], PetscScalar bcval[]) 3437 .ve 3438 + dim - the spatial dimension 3439 . Nf - the number of fields 3440 . uOff - the offset into u[] and u_t[] for each field 3441 . uOff_x - the offset into u_x[] for each field 3442 . u - each field evaluated at the current point 3443 . u_t - the time derivative of each field evaluated at the current point 3444 . u_x - the gradient of each field evaluated at the current point 3445 . aOff - the offset into a[] and a_t[] for each auxiliary field 3446 . aOff_x - the offset into a_x[] for each auxiliary field 3447 . a - each auxiliary field evaluated at the current point 3448 . a_t - the time derivative of each auxiliary field evaluated at the current point 3449 . a_x - the gradient of auxiliary each field evaluated at the current point 3450 . t - current time 3451 . x - coordinates of the current point 3452 . numConstants - number of constant parameters 3453 . constants - constant parameters 3454 - bcval - output values at the current point 3455 3456 Notes: 3457 The pointwise functions are used to provide boundary values for essential boundary 3458 conditions. In FEM, they are acting upon by dual basis functionals to generate FEM 3459 coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary 3460 integrals should be performed, using the kernels from `PetscDSSetBdResidual()`. 3461 3462 .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()` 3463 @*/ 3464 PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd) 3465 { 3466 DSBoundary head = ds->boundary, b; 3467 PetscInt n = 0; 3468 const char *lname; 3469 3470 PetscFunctionBegin; 3471 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3472 PetscValidLogicalCollectiveEnum(ds, type, 2); 3473 PetscAssertPointer(name, 3); 3474 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4); 3475 PetscValidLogicalCollectiveInt(ds, Nv, 5); 3476 PetscValidLogicalCollectiveInt(ds, field, 7); 3477 PetscValidLogicalCollectiveInt(ds, Nc, 8); 3478 PetscCheck(field >= 0 && field < ds->Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", field, ds->Nf); 3479 if (Nc > 0) { 3480 PetscInt *fcomps; 3481 PetscInt c; 3482 3483 PetscCall(PetscDSGetComponents(ds, &fcomps)); 3484 PetscCheck(Nc <= fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Number of constrained components %" PetscInt_FMT " > %" PetscInt_FMT " components for field %" PetscInt_FMT, Nc, fcomps[field], field); 3485 for (c = 0; c < Nc; ++c) { 3486 PetscCheck(comps[c] >= 0 && comps[c] < fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Constrained component[%" PetscInt_FMT "] %" PetscInt_FMT " not in [0, %" PetscInt_FMT ") components for field %" PetscInt_FMT, c, comps[c], fcomps[field], field); 3487 } 3488 } 3489 PetscCall(PetscNew(&b)); 3490 PetscCall(PetscStrallocpy(name, (char **)&b->name)); 3491 PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf)); 3492 PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf)); 3493 PetscCall(PetscMalloc1(Nv, &b->values)); 3494 if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv)); 3495 PetscCall(PetscMalloc1(Nc, &b->comps)); 3496 if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc)); 3497 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 3498 PetscCall(PetscStrallocpy(lname, (char **)&b->lname)); 3499 b->type = type; 3500 b->label = label; 3501 b->Nv = Nv; 3502 b->field = field; 3503 b->Nc = Nc; 3504 b->func = bcFunc; 3505 b->func_t = bcFunc_t; 3506 b->ctx = ctx; 3507 b->next = NULL; 3508 /* Append to linked list so that we can preserve the order */ 3509 if (!head) ds->boundary = b; 3510 while (head) { 3511 if (!head->next) { 3512 head->next = b; 3513 head = b; 3514 } 3515 head = head->next; 3516 ++n; 3517 } 3518 if (bd) { 3519 PetscAssertPointer(bd, 13); 3520 *bd = n; 3521 } 3522 PetscFunctionReturn(PETSC_SUCCESS); 3523 } 3524 3525 // PetscClangLinter pragma ignore: -fdoc-section-header-unknown 3526 /*@C 3527 PetscDSAddBoundaryByName - Add a boundary condition to the model. 3528 3529 Collective 3530 3531 Input Parameters: 3532 + ds - The `PetscDS` object 3533 . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann) 3534 . name - The BC name 3535 . lname - The naem of the label defining constrained points 3536 . Nv - The number of `DMLabel` values for constrained points 3537 . values - An array of label values for constrained points 3538 . field - The field to constrain 3539 . Nc - The number of constrained field components (0 will constrain all fields) 3540 . comps - An array of constrained component numbers 3541 . bcFunc - A pointwise function giving boundary values 3542 . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL 3543 - ctx - An optional user context for bcFunc 3544 3545 Output Parameter: 3546 . bd - The boundary number 3547 3548 Options Database Keys: 3549 + -bc_<boundary name> <num> - Overrides the boundary ids 3550 - -bc_<boundary name>_comp <num> - Overrides the boundary components 3551 3552 Calling Sequence of `bcFunc` and `bcFunc_t`: 3553 If the type is `DM_BC_ESSENTIAL` 3554 .vb 3555 void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[]) 3556 .ve 3557 If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, 3558 .vb 3559 void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux, 3560 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 3561 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 3562 PetscReal time, const PetscReal x[], PetscScalar bcval[]) 3563 .ve 3564 + dim - the spatial dimension 3565 . Nf - the number of fields 3566 . uOff - the offset into u[] and u_t[] for each field 3567 . uOff_x - the offset into u_x[] for each field 3568 . u - each field evaluated at the current point 3569 . u_t - the time derivative of each field evaluated at the current point 3570 . u_x - the gradient of each field evaluated at the current point 3571 . aOff - the offset into a[] and a_t[] for each auxiliary field 3572 . aOff_x - the offset into a_x[] for each auxiliary field 3573 . a - each auxiliary field evaluated at the current point 3574 . a_t - the time derivative of each auxiliary field evaluated at the current point 3575 . a_x - the gradient of auxiliary each field evaluated at the current point 3576 . t - current time 3577 . x - coordinates of the current point 3578 . numConstants - number of constant parameters 3579 . constants - constant parameters 3580 - bcval - output values at the current point 3581 3582 Level: developer 3583 3584 Notes: 3585 The pointwise functions are used to provide boundary values for essential boundary 3586 conditions. In FEM, they are acting upon by dual basis functionals to generate FEM 3587 coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary 3588 integrals should be performed, using the kernels from `PetscDSSetBdResidual()`. 3589 3590 This function should only be used with `DMFOREST` currently, since labels cannot be defined before the underlying `DMPLEX` is built. 3591 3592 .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()` 3593 @*/ 3594 PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd) 3595 { 3596 DSBoundary head = ds->boundary, b; 3597 PetscInt n = 0; 3598 3599 PetscFunctionBegin; 3600 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3601 PetscValidLogicalCollectiveEnum(ds, type, 2); 3602 PetscAssertPointer(name, 3); 3603 PetscAssertPointer(lname, 4); 3604 PetscValidLogicalCollectiveInt(ds, Nv, 5); 3605 PetscValidLogicalCollectiveInt(ds, field, 7); 3606 PetscValidLogicalCollectiveInt(ds, Nc, 8); 3607 PetscCall(PetscNew(&b)); 3608 PetscCall(PetscStrallocpy(name, (char **)&b->name)); 3609 PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf)); 3610 PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf)); 3611 PetscCall(PetscMalloc1(Nv, &b->values)); 3612 if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv)); 3613 PetscCall(PetscMalloc1(Nc, &b->comps)); 3614 if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc)); 3615 PetscCall(PetscStrallocpy(lname, (char **)&b->lname)); 3616 b->type = type; 3617 b->label = NULL; 3618 b->Nv = Nv; 3619 b->field = field; 3620 b->Nc = Nc; 3621 b->func = bcFunc; 3622 b->func_t = bcFunc_t; 3623 b->ctx = ctx; 3624 b->next = NULL; 3625 /* Append to linked list so that we can preserve the order */ 3626 if (!head) ds->boundary = b; 3627 while (head) { 3628 if (!head->next) { 3629 head->next = b; 3630 head = b; 3631 } 3632 head = head->next; 3633 ++n; 3634 } 3635 if (bd) { 3636 PetscAssertPointer(bd, 13); 3637 *bd = n; 3638 } 3639 PetscFunctionReturn(PETSC_SUCCESS); 3640 } 3641 3642 /*@C 3643 PetscDSUpdateBoundary - Change a boundary condition for the model. 3644 3645 Input Parameters: 3646 + ds - The `PetscDS` object 3647 . bd - The boundary condition number 3648 . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann) 3649 . name - The BC name 3650 . label - The label defining constrained points 3651 . Nv - The number of `DMLabel` ids for constrained points 3652 . values - An array of ids for constrained points 3653 . field - The field to constrain 3654 . Nc - The number of constrained field components 3655 . comps - An array of constrained component numbers 3656 . bcFunc - A pointwise function giving boundary values 3657 . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL 3658 - ctx - An optional user context for bcFunc 3659 3660 Level: developer 3661 3662 Notes: 3663 The pointwise functions are used to provide boundary values for essential boundary 3664 conditions. In FEM, they are acting upon by dual basis functionals to generate FEM 3665 coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary 3666 integrals should be performed, using the kernels from `PetscDSSetBdResidual()`. 3667 3668 The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from `PetscDSGetNumBoundary()`. 3669 See `PetscDSAddBoundary()` for a description of the calling sequences for the callbacks. 3670 3671 .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel` 3672 @*/ 3673 PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx) 3674 { 3675 DSBoundary b = ds->boundary; 3676 PetscInt n = 0; 3677 3678 PetscFunctionBegin; 3679 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3680 while (b) { 3681 if (n == bd) break; 3682 b = b->next; 3683 ++n; 3684 } 3685 PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n); 3686 if (name) { 3687 PetscCall(PetscFree(b->name)); 3688 PetscCall(PetscStrallocpy(name, (char **)&b->name)); 3689 } 3690 b->type = type; 3691 if (label) { 3692 const char *name; 3693 3694 b->label = label; 3695 PetscCall(PetscFree(b->lname)); 3696 PetscCall(PetscObjectGetName((PetscObject)label, &name)); 3697 PetscCall(PetscStrallocpy(name, (char **)&b->lname)); 3698 } 3699 if (Nv >= 0) { 3700 b->Nv = Nv; 3701 PetscCall(PetscFree(b->values)); 3702 PetscCall(PetscMalloc1(Nv, &b->values)); 3703 if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv)); 3704 } 3705 if (field >= 0) b->field = field; 3706 if (Nc >= 0) { 3707 b->Nc = Nc; 3708 PetscCall(PetscFree(b->comps)); 3709 PetscCall(PetscMalloc1(Nc, &b->comps)); 3710 if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc)); 3711 } 3712 if (bcFunc) b->func = bcFunc; 3713 if (bcFunc_t) b->func_t = bcFunc_t; 3714 if (ctx) b->ctx = ctx; 3715 PetscFunctionReturn(PETSC_SUCCESS); 3716 } 3717 3718 /*@ 3719 PetscDSGetNumBoundary - Get the number of registered BC 3720 3721 Input Parameter: 3722 . ds - The `PetscDS` object 3723 3724 Output Parameter: 3725 . numBd - The number of BC 3726 3727 Level: intermediate 3728 3729 .seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()` 3730 @*/ 3731 PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd) 3732 { 3733 DSBoundary b = ds->boundary; 3734 3735 PetscFunctionBegin; 3736 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3737 PetscAssertPointer(numBd, 2); 3738 *numBd = 0; 3739 while (b) { 3740 ++(*numBd); 3741 b = b->next; 3742 } 3743 PetscFunctionReturn(PETSC_SUCCESS); 3744 } 3745 3746 /*@C 3747 PetscDSGetBoundary - Gets a boundary condition to the model 3748 3749 Input Parameters: 3750 + ds - The `PetscDS` object 3751 - bd - The BC number 3752 3753 Output Parameters: 3754 + wf - The `PetscWeakForm` holding the pointwise functions 3755 . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann) 3756 . name - The BC name 3757 . label - The label defining constrained points 3758 . Nv - The number of `DMLabel` ids for constrained points 3759 . values - An array of ids for constrained points 3760 . field - The field to constrain 3761 . Nc - The number of constrained field components 3762 . comps - An array of constrained component numbers 3763 . func - A pointwise function giving boundary values 3764 . func_t - A pointwise function giving the time derivative of the boundary values 3765 - ctx - An optional user context for bcFunc 3766 3767 Options Database Keys: 3768 + -bc_<boundary name> <num> - Overrides the boundary ids 3769 - -bc_<boundary name>_comp <num> - Overrides the boundary components 3770 3771 Level: developer 3772 3773 .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel` 3774 @*/ 3775 PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, PetscWeakForm *wf, DMBoundaryConditionType *type, const char *name[], DMLabel *label, PetscInt *Nv, const PetscInt *values[], PetscInt *field, PetscInt *Nc, const PetscInt *comps[], void (**func)(void), void (**func_t)(void), void **ctx) 3776 { 3777 DSBoundary b = ds->boundary; 3778 PetscInt n = 0; 3779 3780 PetscFunctionBegin; 3781 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3782 while (b) { 3783 if (n == bd) break; 3784 b = b->next; 3785 ++n; 3786 } 3787 PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n); 3788 if (wf) { 3789 PetscAssertPointer(wf, 3); 3790 *wf = b->wf; 3791 } 3792 if (type) { 3793 PetscAssertPointer(type, 4); 3794 *type = b->type; 3795 } 3796 if (name) { 3797 PetscAssertPointer(name, 5); 3798 *name = b->name; 3799 } 3800 if (label) { 3801 PetscAssertPointer(label, 6); 3802 *label = b->label; 3803 } 3804 if (Nv) { 3805 PetscAssertPointer(Nv, 7); 3806 *Nv = b->Nv; 3807 } 3808 if (values) { 3809 PetscAssertPointer(values, 8); 3810 *values = b->values; 3811 } 3812 if (field) { 3813 PetscAssertPointer(field, 9); 3814 *field = b->field; 3815 } 3816 if (Nc) { 3817 PetscAssertPointer(Nc, 10); 3818 *Nc = b->Nc; 3819 } 3820 if (comps) { 3821 PetscAssertPointer(comps, 11); 3822 *comps = b->comps; 3823 } 3824 if (func) { 3825 PetscAssertPointer(func, 12); 3826 *func = b->func; 3827 } 3828 if (func_t) { 3829 PetscAssertPointer(func_t, 13); 3830 *func_t = b->func_t; 3831 } 3832 if (ctx) { 3833 PetscAssertPointer(ctx, 14); 3834 *ctx = b->ctx; 3835 } 3836 PetscFunctionReturn(PETSC_SUCCESS); 3837 } 3838 3839 /*@ 3840 PetscDSUpdateBoundaryLabels - Update `DMLabel` in each boundary condition using the label name and the input `DM` 3841 3842 Not Collective 3843 3844 Input Parameters: 3845 + ds - The source `PetscDS` object 3846 - dm - The `DM` holding labels 3847 3848 Level: intermediate 3849 3850 .seealso: `PetscDS`, `DMBoundary`, `DM`, `PetscDSCopyBoundary()`, `PetscDSCreate()`, `DMGetLabel()` 3851 @*/ 3852 PetscErrorCode PetscDSUpdateBoundaryLabels(PetscDS ds, DM dm) 3853 { 3854 DSBoundary b; 3855 3856 PetscFunctionBegin; 3857 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3858 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 3859 for (b = ds->boundary; b; b = b->next) { 3860 if (b->lname) PetscCall(DMGetLabel(dm, b->lname, &b->label)); 3861 } 3862 PetscFunctionReturn(PETSC_SUCCESS); 3863 } 3864 3865 static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew) 3866 { 3867 PetscFunctionBegin; 3868 PetscCall(PetscNew(bNew)); 3869 PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf)); 3870 PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf)); 3871 PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name))); 3872 PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname))); 3873 (*bNew)->type = b->type; 3874 (*bNew)->label = b->label; 3875 (*bNew)->Nv = b->Nv; 3876 PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values)); 3877 PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv)); 3878 (*bNew)->field = b->field; 3879 (*bNew)->Nc = b->Nc; 3880 PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps)); 3881 PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc)); 3882 (*bNew)->func = b->func; 3883 (*bNew)->func_t = b->func_t; 3884 (*bNew)->ctx = b->ctx; 3885 PetscFunctionReturn(PETSC_SUCCESS); 3886 } 3887 3888 /*@ 3889 PetscDSCopyBoundary - Copy all boundary condition objects to the new problem 3890 3891 Not Collective 3892 3893 Input Parameters: 3894 + ds - The source `PetscDS` object 3895 . numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields 3896 - fields - The selected fields, or NULL for all fields 3897 3898 Output Parameter: 3899 . newds - The target `PetscDS`, now with a copy of the boundary conditions 3900 3901 Level: intermediate 3902 3903 .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()` 3904 @*/ 3905 PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds) 3906 { 3907 DSBoundary b, *lastnext; 3908 3909 PetscFunctionBegin; 3910 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 3911 PetscValidHeaderSpecific(newds, PETSCDS_CLASSID, 4); 3912 if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS); 3913 PetscCall(PetscDSDestroyBoundary(newds)); 3914 lastnext = &newds->boundary; 3915 for (b = ds->boundary; b; b = b->next) { 3916 DSBoundary bNew; 3917 PetscInt fieldNew = -1; 3918 3919 if (numFields > 0 && fields) { 3920 PetscInt f; 3921 3922 for (f = 0; f < numFields; ++f) 3923 if (b->field == fields[f]) break; 3924 if (f == numFields) continue; 3925 fieldNew = f; 3926 } 3927 PetscCall(DSBoundaryDuplicate_Internal(b, &bNew)); 3928 bNew->field = fieldNew < 0 ? b->field : fieldNew; 3929 *lastnext = bNew; 3930 lastnext = &bNew->next; 3931 } 3932 PetscFunctionReturn(PETSC_SUCCESS); 3933 } 3934 3935 /*@ 3936 PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS` 3937 3938 Not Collective 3939 3940 Input Parameter: 3941 . ds - The `PetscDS` object 3942 3943 Level: intermediate 3944 3945 .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()` 3946 @*/ 3947 PetscErrorCode PetscDSDestroyBoundary(PetscDS ds) 3948 { 3949 DSBoundary next = ds->boundary; 3950 3951 PetscFunctionBegin; 3952 while (next) { 3953 DSBoundary b = next; 3954 3955 next = b->next; 3956 PetscCall(PetscWeakFormDestroy(&b->wf)); 3957 PetscCall(PetscFree(b->name)); 3958 PetscCall(PetscFree(b->lname)); 3959 PetscCall(PetscFree(b->values)); 3960 PetscCall(PetscFree(b->comps)); 3961 PetscCall(PetscFree(b)); 3962 } 3963 PetscFunctionReturn(PETSC_SUCCESS); 3964 } 3965 3966 /*@ 3967 PetscDSSelectDiscretizations - Copy discretizations to the new problem with different field layout 3968 3969 Not Collective 3970 3971 Input Parameters: 3972 + prob - The `PetscDS` object 3973 . numFields - Number of new fields 3974 . fields - Old field number for each new field 3975 . minDegree - Minimum degree for a discretization, or `PETSC_DETERMINE` for no limit 3976 - maxDegree - Maximum degree for a discretization, or `PETSC_DETERMINE` for no limit 3977 3978 Output Parameter: 3979 . newprob - The `PetscDS` copy 3980 3981 Level: intermediate 3982 3983 .seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()` 3984 @*/ 3985 PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscInt minDegree, PetscInt maxDegree, PetscDS newprob) 3986 { 3987 PetscInt Nf, Nfn, fn; 3988 3989 PetscFunctionBegin; 3990 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 3991 if (fields) PetscAssertPointer(fields, 3); 3992 PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 6); 3993 PetscCall(PetscDSGetNumFields(prob, &Nf)); 3994 PetscCall(PetscDSGetNumFields(newprob, &Nfn)); 3995 numFields = numFields < 0 ? Nf : numFields; 3996 for (fn = 0; fn < numFields; ++fn) { 3997 const PetscInt f = fields ? fields[fn] : fn; 3998 PetscObject disc; 3999 PetscClassId id; 4000 4001 if (f >= Nf) continue; 4002 PetscCall(PetscDSGetDiscretization(prob, f, &disc)); 4003 PetscCallContinue(PetscObjectGetClassId(disc, &id)); 4004 if (id == PETSCFE_CLASSID) { 4005 PetscFE fe; 4006 4007 PetscCall(PetscFELimitDegree((PetscFE)disc, minDegree, maxDegree, &fe)); 4008 PetscCall(PetscDSSetDiscretization(newprob, fn, (PetscObject)fe)); 4009 PetscCall(PetscFEDestroy(&fe)); 4010 } else { 4011 PetscCall(PetscDSSetDiscretization(newprob, fn, disc)); 4012 } 4013 } 4014 PetscFunctionReturn(PETSC_SUCCESS); 4015 } 4016 4017 /*@ 4018 PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout 4019 4020 Not Collective 4021 4022 Input Parameters: 4023 + prob - The `PetscDS` object 4024 . numFields - Number of new fields 4025 - fields - Old field number for each new field 4026 4027 Output Parameter: 4028 . newprob - The `PetscDS` copy 4029 4030 Level: intermediate 4031 4032 .seealso: `PetscDS`, `PetscDSSelectDiscretizations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()` 4033 @*/ 4034 PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob) 4035 { 4036 PetscInt Nf, Nfn, fn, gn; 4037 4038 PetscFunctionBegin; 4039 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 4040 if (fields) PetscAssertPointer(fields, 3); 4041 PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 4); 4042 PetscCall(PetscDSGetNumFields(prob, &Nf)); 4043 PetscCall(PetscDSGetNumFields(newprob, &Nfn)); 4044 PetscCheck(numFields <= Nfn, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields %" PetscInt_FMT " to transfer must not be greater then the total number of fields %" PetscInt_FMT, numFields, Nfn); 4045 for (fn = 0; fn < numFields; ++fn) { 4046 const PetscInt f = fields ? fields[fn] : fn; 4047 PetscPointFunc obj; 4048 PetscPointFunc f0, f1; 4049 PetscBdPointFunc f0Bd, f1Bd; 4050 PetscRiemannFunc r; 4051 4052 if (f >= Nf) continue; 4053 PetscCall(PetscDSGetObjective(prob, f, &obj)); 4054 PetscCall(PetscDSGetResidual(prob, f, &f0, &f1)); 4055 PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd)); 4056 PetscCall(PetscDSGetRiemannSolver(prob, f, &r)); 4057 PetscCall(PetscDSSetObjective(newprob, fn, obj)); 4058 PetscCall(PetscDSSetResidual(newprob, fn, f0, f1)); 4059 PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd)); 4060 PetscCall(PetscDSSetRiemannSolver(newprob, fn, r)); 4061 for (gn = 0; gn < numFields; ++gn) { 4062 const PetscInt g = fields ? fields[gn] : gn; 4063 PetscPointJac g0, g1, g2, g3; 4064 PetscPointJac g0p, g1p, g2p, g3p; 4065 PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd; 4066 4067 if (g >= Nf) continue; 4068 PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3)); 4069 PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p)); 4070 PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd)); 4071 PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3)); 4072 PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p)); 4073 PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd)); 4074 } 4075 } 4076 PetscFunctionReturn(PETSC_SUCCESS); 4077 } 4078 4079 /*@ 4080 PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS` 4081 4082 Not Collective 4083 4084 Input Parameter: 4085 . prob - The `PetscDS` object 4086 4087 Output Parameter: 4088 . newprob - The `PetscDS` copy 4089 4090 Level: intermediate 4091 4092 .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()` 4093 @*/ 4094 PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob) 4095 { 4096 PetscWeakForm wf, newwf; 4097 PetscInt Nf, Ng; 4098 4099 PetscFunctionBegin; 4100 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 4101 PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 2); 4102 PetscCall(PetscDSGetNumFields(prob, &Nf)); 4103 PetscCall(PetscDSGetNumFields(newprob, &Ng)); 4104 PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng); 4105 PetscCall(PetscDSGetWeakForm(prob, &wf)); 4106 PetscCall(PetscDSGetWeakForm(newprob, &newwf)); 4107 PetscCall(PetscWeakFormCopy(wf, newwf)); 4108 PetscFunctionReturn(PETSC_SUCCESS); 4109 } 4110 4111 /*@ 4112 PetscDSCopyConstants - Copy all constants to another `PetscDS` 4113 4114 Not Collective 4115 4116 Input Parameter: 4117 . prob - The `PetscDS` object 4118 4119 Output Parameter: 4120 . newprob - The `PetscDS` copy 4121 4122 Level: intermediate 4123 4124 .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()` 4125 @*/ 4126 PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob) 4127 { 4128 PetscInt Nc; 4129 const PetscScalar *constants; 4130 4131 PetscFunctionBegin; 4132 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 4133 PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 2); 4134 PetscCall(PetscDSGetConstants(prob, &Nc, &constants)); 4135 PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants)); 4136 PetscFunctionReturn(PETSC_SUCCESS); 4137 } 4138 4139 /*@ 4140 PetscDSCopyExactSolutions - Copy all exact solutions to another `PetscDS` 4141 4142 Not Collective 4143 4144 Input Parameter: 4145 . ds - The `PetscDS` object 4146 4147 Output Parameter: 4148 . newds - The `PetscDS` copy 4149 4150 Level: intermediate 4151 4152 .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()` 4153 @*/ 4154 PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds) 4155 { 4156 PetscSimplePointFn *sol; 4157 void *ctx; 4158 PetscInt Nf, f; 4159 4160 PetscFunctionBegin; 4161 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 4162 PetscValidHeaderSpecific(newds, PETSCDS_CLASSID, 2); 4163 PetscCall(PetscDSGetNumFields(ds, &Nf)); 4164 for (f = 0; f < Nf; ++f) { 4165 PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx)); 4166 PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx)); 4167 PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx)); 4168 PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx)); 4169 } 4170 PetscFunctionReturn(PETSC_SUCCESS); 4171 } 4172 4173 PetscErrorCode PetscDSCopy(PetscDS ds, PetscInt minDegree, PetscInt maxDegree, DM dmNew, PetscDS dsNew) 4174 { 4175 DSBoundary b; 4176 PetscInt cdim, Nf, f, d; 4177 PetscBool isCohesive; 4178 void *ctx; 4179 4180 PetscFunctionBegin; 4181 PetscCall(PetscDSCopyConstants(ds, dsNew)); 4182 PetscCall(PetscDSCopyExactSolutions(ds, dsNew)); 4183 PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, minDegree, maxDegree, dsNew)); 4184 PetscCall(PetscDSCopyEquations(ds, dsNew)); 4185 PetscCall(PetscDSGetNumFields(ds, &Nf)); 4186 for (f = 0; f < Nf; ++f) { 4187 PetscCall(PetscDSGetContext(ds, f, &ctx)); 4188 PetscCall(PetscDSSetContext(dsNew, f, ctx)); 4189 PetscCall(PetscDSGetCohesive(ds, f, &isCohesive)); 4190 PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive)); 4191 PetscCall(PetscDSGetJetDegree(ds, f, &d)); 4192 PetscCall(PetscDSSetJetDegree(dsNew, f, d)); 4193 } 4194 if (Nf) { 4195 PetscCall(PetscDSGetCoordinateDimension(ds, &cdim)); 4196 PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim)); 4197 } 4198 PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew)); 4199 for (b = dsNew->boundary; b; b = b->next) { 4200 PetscCall(DMGetLabel(dmNew, b->lname, &b->label)); 4201 /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */ 4202 //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name); 4203 } 4204 PetscFunctionReturn(PETSC_SUCCESS); 4205 } 4206 4207 PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob) 4208 { 4209 PetscInt dim, Nf, f; 4210 4211 PetscFunctionBegin; 4212 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 4213 PetscAssertPointer(subprob, 3); 4214 if (height == 0) { 4215 *subprob = prob; 4216 PetscFunctionReturn(PETSC_SUCCESS); 4217 } 4218 PetscCall(PetscDSGetNumFields(prob, &Nf)); 4219 PetscCall(PetscDSGetSpatialDimension(prob, &dim)); 4220 PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height); 4221 if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs)); 4222 if (!prob->subprobs[height - 1]) { 4223 PetscInt cdim; 4224 4225 PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1])); 4226 PetscCall(PetscDSGetCoordinateDimension(prob, &cdim)); 4227 PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim)); 4228 for (f = 0; f < Nf; ++f) { 4229 PetscFE subfe; 4230 PetscObject obj; 4231 PetscClassId id; 4232 4233 PetscCall(PetscDSGetDiscretization(prob, f, &obj)); 4234 PetscCall(PetscObjectGetClassId(obj, &id)); 4235 if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe)); 4236 else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f); 4237 PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe)); 4238 } 4239 } 4240 *subprob = prob->subprobs[height - 1]; 4241 PetscFunctionReturn(PETSC_SUCCESS); 4242 } 4243 4244 PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm) 4245 { 4246 IS permIS; 4247 PetscQuadrature quad; 4248 DMPolytopeType ct; 4249 const PetscInt *perm; 4250 PetscInt Na, Nq; 4251 4252 PetscFunctionBeginHot; 4253 PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad)); 4254 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL)); 4255 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 4256 PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq); 4257 Na = DMPolytopeTypeGetNumArrangements(ct) / 2; 4258 PetscCheck(ornt >= -Na && ornt < Na, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Orientation %" PetscInt_FMT " of %s is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", ornt, DMPolytopeTypes[ct], -Na, Na); 4259 if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct])); 4260 permIS = ds->quadPerm[(PetscInt)ct][ornt + Na]; 4261 PetscCall(ISGetIndices(permIS, &perm)); 4262 *qperm = perm[q]; 4263 PetscCall(ISRestoreIndices(permIS, &perm)); 4264 PetscFunctionReturn(PETSC_SUCCESS); 4265 } 4266 4267 PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype) 4268 { 4269 PetscObject obj; 4270 PetscClassId id; 4271 PetscInt Nf; 4272 4273 PetscFunctionBegin; 4274 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 4275 PetscAssertPointer(disctype, 3); 4276 *disctype = PETSC_DISC_NONE; 4277 PetscCall(PetscDSGetNumFields(ds, &Nf)); 4278 PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf); 4279 PetscCall(PetscDSGetDiscretization(ds, f, &obj)); 4280 if (obj) { 4281 PetscCall(PetscObjectGetClassId(obj, &id)); 4282 if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE; 4283 else *disctype = PETSC_DISC_FV; 4284 } 4285 PetscFunctionReturn(PETSC_SUCCESS); 4286 } 4287 4288 static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds) 4289 { 4290 PetscFunctionBegin; 4291 PetscCall(PetscFree(ds->data)); 4292 PetscFunctionReturn(PETSC_SUCCESS); 4293 } 4294 4295 static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds) 4296 { 4297 PetscFunctionBegin; 4298 ds->ops->setfromoptions = NULL; 4299 ds->ops->setup = NULL; 4300 ds->ops->view = NULL; 4301 ds->ops->destroy = PetscDSDestroy_Basic; 4302 PetscFunctionReturn(PETSC_SUCCESS); 4303 } 4304 4305 /*MC 4306 PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions 4307 4308 Level: intermediate 4309 4310 .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()` 4311 M*/ 4312 4313 PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds) 4314 { 4315 PetscDS_Basic *b; 4316 4317 PetscFunctionBegin; 4318 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 4319 PetscCall(PetscNew(&b)); 4320 ds->data = b; 4321 4322 PetscCall(PetscDSInitialize_Basic(ds)); 4323 PetscFunctionReturn(PETSC_SUCCESS); 4324 } 4325