1 /* Basis Jet Tabulation
2
3 We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We
4 follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can
5 be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis
6 as a prime basis.
7
8 \psi_i = \sum_k \alpha_{ki} \phi_k
9
10 Our nodal basis is defined in terms of the dual basis $n_j$
11
12 n_j \cdot \psi_i = \delta_{ji}
13
14 and we may act on the first equation to obtain
15
16 n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k
17 \delta_{ji} = \sum_k \alpha_{ki} V_{jk}
18 I = V \alpha
19
20 so the coefficients of the nodal basis in the prime basis are
21
22 \alpha = V^{-1}
23
24 We will define the dual basis vectors $n_j$ using a quadrature rule.
25
26 Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials
27 (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can
28 be implemented exactly as in FIAT using functionals $L_j$.
29
30 I will have to count the degrees correctly for the Legendre product when we are on simplices.
31
32 We will have three objects:
33 - Space, P: this just need point evaluation I think
34 - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q
35 - FEM: This keeps {P, P', Q}
36 */
37 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
38 #include <petscdmplex.h>
39
40 PetscBool FEcite = PETSC_FALSE;
41 const char FECitation[] = "@article{kirby2004,\n"
42 " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n"
43 " journal = {ACM Transactions on Mathematical Software},\n"
44 " author = {Robert C. Kirby},\n"
45 " volume = {30},\n"
46 " number = {4},\n"
47 " pages = {502--516},\n"
48 " doi = {10.1145/1039813.1039820},\n"
49 " year = {2004}\n}\n";
50
51 PetscClassId PETSCFE_CLASSID = 0;
52
53 PetscLogEvent PETSCFE_SetUp;
54
55 PetscFunctionList PetscFEList = NULL;
56 PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
57
58 /*@C
59 PetscFERegister - Adds a new PetscFE implementation
60
61 Not Collective
62
63 Input Parameters:
64 + name - The name of a new user-defined creation routine
65 - create_func - The creation routine itself
66
67 Notes:
68 PetscFERegister() may be called multiple times to add several user-defined PetscFEs
69
70 Sample usage:
71 .vb
72 PetscFERegister("my_fe", MyPetscFECreate);
73 .ve
74
75 Then, your PetscFE type can be chosen with the procedural interface via
76 .vb
77 PetscFECreate(MPI_Comm, PetscFE *);
78 PetscFESetType(PetscFE, "my_fe");
79 .ve
80 or at runtime via the option
81 .vb
82 -petscfe_type my_fe
83 .ve
84
85 Level: advanced
86
87 .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
88
89 @*/
PetscFERegister(const char sname[],PetscErrorCode (* function)(PetscFE))90 PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
91 {
92 PetscErrorCode ierr;
93
94 PetscFunctionBegin;
95 ierr = PetscFunctionListAdd(&PetscFEList, sname, function);CHKERRQ(ierr);
96 PetscFunctionReturn(0);
97 }
98
99 /*@C
100 PetscFESetType - Builds a particular PetscFE
101
102 Collective on fem
103
104 Input Parameters:
105 + fem - The PetscFE object
106 - name - The kind of FEM space
107
108 Options Database Key:
109 . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
110
111 Level: intermediate
112
113 .seealso: PetscFEGetType(), PetscFECreate()
114 @*/
PetscFESetType(PetscFE fem,PetscFEType name)115 PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
116 {
117 PetscErrorCode (*r)(PetscFE);
118 PetscBool match;
119 PetscErrorCode ierr;
120
121 PetscFunctionBegin;
122 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
123 ierr = PetscObjectTypeCompare((PetscObject) fem, name, &match);CHKERRQ(ierr);
124 if (match) PetscFunctionReturn(0);
125
126 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);}
127 ierr = PetscFunctionListFind(PetscFEList, name, &r);CHKERRQ(ierr);
128 if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
129
130 if (fem->ops->destroy) {
131 ierr = (*fem->ops->destroy)(fem);CHKERRQ(ierr);
132 fem->ops->destroy = NULL;
133 }
134 ierr = (*r)(fem);CHKERRQ(ierr);
135 ierr = PetscObjectChangeTypeName((PetscObject) fem, name);CHKERRQ(ierr);
136 PetscFunctionReturn(0);
137 }
138
139 /*@C
140 PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
141
142 Not Collective
143
144 Input Parameter:
145 . fem - The PetscFE
146
147 Output Parameter:
148 . name - The PetscFE type name
149
150 Level: intermediate
151
152 .seealso: PetscFESetType(), PetscFECreate()
153 @*/
PetscFEGetType(PetscFE fem,PetscFEType * name)154 PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
155 {
156 PetscErrorCode ierr;
157
158 PetscFunctionBegin;
159 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
160 PetscValidPointer(name, 2);
161 if (!PetscFERegisterAllCalled) {
162 ierr = PetscFERegisterAll();CHKERRQ(ierr);
163 }
164 *name = ((PetscObject) fem)->type_name;
165 PetscFunctionReturn(0);
166 }
167
168 /*@C
169 PetscFEViewFromOptions - View from Options
170
171 Collective on PetscFE
172
173 Input Parameters:
174 + A - the PetscFE object
175 . obj - Optional object
176 - name - command line option
177
178 Level: intermediate
179 .seealso: PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate()
180 @*/
PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])181 PetscErrorCode PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[])
182 {
183 PetscErrorCode ierr;
184
185 PetscFunctionBegin;
186 PetscValidHeaderSpecific(A,PETSCFE_CLASSID,1);
187 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr);
188 PetscFunctionReturn(0);
189 }
190
191 /*@C
192 PetscFEView - Views a PetscFE
193
194 Collective on fem
195
196 Input Parameter:
197 + fem - the PetscFE object to view
198 - viewer - the viewer
199
200 Level: beginner
201
202 .seealso PetscFEDestroy()
203 @*/
PetscFEView(PetscFE fem,PetscViewer viewer)204 PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer)
205 {
206 PetscBool iascii;
207 PetscErrorCode ierr;
208
209 PetscFunctionBegin;
210 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
211 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
212 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);CHKERRQ(ierr);}
213 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);CHKERRQ(ierr);
214 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
215 if (fem->ops->view) {ierr = (*fem->ops->view)(fem, viewer);CHKERRQ(ierr);}
216 PetscFunctionReturn(0);
217 }
218
219 /*@
220 PetscFESetFromOptions - sets parameters in a PetscFE from the options database
221
222 Collective on fem
223
224 Input Parameter:
225 . fem - the PetscFE object to set options for
226
227 Options Database:
228 + -petscfe_num_blocks - the number of cell blocks to integrate concurrently
229 - -petscfe_num_batches - the number of cell batches to integrate serially
230
231 Level: intermediate
232
233 .seealso PetscFEView()
234 @*/
PetscFESetFromOptions(PetscFE fem)235 PetscErrorCode PetscFESetFromOptions(PetscFE fem)
236 {
237 const char *defaultType;
238 char name[256];
239 PetscBool flg;
240 PetscErrorCode ierr;
241
242 PetscFunctionBegin;
243 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
244 if (!((PetscObject) fem)->type_name) {
245 defaultType = PETSCFEBASIC;
246 } else {
247 defaultType = ((PetscObject) fem)->type_name;
248 }
249 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);}
250
251 ierr = PetscObjectOptionsBegin((PetscObject) fem);CHKERRQ(ierr);
252 ierr = PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);CHKERRQ(ierr);
253 if (flg) {
254 ierr = PetscFESetType(fem, name);CHKERRQ(ierr);
255 } else if (!((PetscObject) fem)->type_name) {
256 ierr = PetscFESetType(fem, defaultType);CHKERRQ(ierr);
257 }
258 ierr = PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);CHKERRQ(ierr);
259 ierr = PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);CHKERRQ(ierr);
260 if (fem->ops->setfromoptions) {
261 ierr = (*fem->ops->setfromoptions)(PetscOptionsObject,fem);CHKERRQ(ierr);
262 }
263 /* process any options handlers added with PetscObjectAddOptionsHandler() */
264 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);CHKERRQ(ierr);
265 ierr = PetscOptionsEnd();CHKERRQ(ierr);
266 ierr = PetscFEViewFromOptions(fem, NULL, "-petscfe_view");CHKERRQ(ierr);
267 PetscFunctionReturn(0);
268 }
269
270 /*@C
271 PetscFESetUp - Construct data structures for the PetscFE
272
273 Collective on fem
274
275 Input Parameter:
276 . fem - the PetscFE object to setup
277
278 Level: intermediate
279
280 .seealso PetscFEView(), PetscFEDestroy()
281 @*/
PetscFESetUp(PetscFE fem)282 PetscErrorCode PetscFESetUp(PetscFE fem)
283 {
284 PetscErrorCode ierr;
285
286 PetscFunctionBegin;
287 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
288 if (fem->setupcalled) PetscFunctionReturn(0);
289 ierr = PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);CHKERRQ(ierr);
290 fem->setupcalled = PETSC_TRUE;
291 if (fem->ops->setup) {ierr = (*fem->ops->setup)(fem);CHKERRQ(ierr);}
292 ierr = PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);CHKERRQ(ierr);
293 PetscFunctionReturn(0);
294 }
295
296 /*@
297 PetscFEDestroy - Destroys a PetscFE object
298
299 Collective on fem
300
301 Input Parameter:
302 . fem - the PetscFE object to destroy
303
304 Level: beginner
305
306 .seealso PetscFEView()
307 @*/
PetscFEDestroy(PetscFE * fem)308 PetscErrorCode PetscFEDestroy(PetscFE *fem)
309 {
310 PetscErrorCode ierr;
311
312 PetscFunctionBegin;
313 if (!*fem) PetscFunctionReturn(0);
314 PetscValidHeaderSpecific((*fem), PETSCFE_CLASSID, 1);
315
316 if (--((PetscObject)(*fem))->refct > 0) {*fem = NULL; PetscFunctionReturn(0);}
317 ((PetscObject) (*fem))->refct = 0;
318
319 if ((*fem)->subspaces) {
320 PetscInt dim, d;
321
322 ierr = PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);CHKERRQ(ierr);
323 for (d = 0; d < dim; ++d) {ierr = PetscFEDestroy(&(*fem)->subspaces[d]);CHKERRQ(ierr);}
324 }
325 ierr = PetscFree((*fem)->subspaces);CHKERRQ(ierr);
326 ierr = PetscFree((*fem)->invV);CHKERRQ(ierr);
327 ierr = PetscTabulationDestroy(&(*fem)->T);CHKERRQ(ierr);
328 ierr = PetscTabulationDestroy(&(*fem)->Tf);CHKERRQ(ierr);
329 ierr = PetscTabulationDestroy(&(*fem)->Tc);CHKERRQ(ierr);
330 ierr = PetscSpaceDestroy(&(*fem)->basisSpace);CHKERRQ(ierr);
331 ierr = PetscDualSpaceDestroy(&(*fem)->dualSpace);CHKERRQ(ierr);
332 ierr = PetscQuadratureDestroy(&(*fem)->quadrature);CHKERRQ(ierr);
333 ierr = PetscQuadratureDestroy(&(*fem)->faceQuadrature);CHKERRQ(ierr);
334
335 if ((*fem)->ops->destroy) {ierr = (*(*fem)->ops->destroy)(*fem);CHKERRQ(ierr);}
336 ierr = PetscHeaderDestroy(fem);CHKERRQ(ierr);
337 PetscFunctionReturn(0);
338 }
339
340 /*@
341 PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
342
343 Collective
344
345 Input Parameter:
346 . comm - The communicator for the PetscFE object
347
348 Output Parameter:
349 . fem - The PetscFE object
350
351 Level: beginner
352
353 .seealso: PetscFESetType(), PETSCFEGALERKIN
354 @*/
PetscFECreate(MPI_Comm comm,PetscFE * fem)355 PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
356 {
357 PetscFE f;
358 PetscErrorCode ierr;
359
360 PetscFunctionBegin;
361 PetscValidPointer(fem, 2);
362 ierr = PetscCitationsRegister(FECitation,&FEcite);CHKERRQ(ierr);
363 *fem = NULL;
364 ierr = PetscFEInitializePackage();CHKERRQ(ierr);
365
366 ierr = PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);CHKERRQ(ierr);
367
368 f->basisSpace = NULL;
369 f->dualSpace = NULL;
370 f->numComponents = 1;
371 f->subspaces = NULL;
372 f->invV = NULL;
373 f->T = NULL;
374 f->Tf = NULL;
375 f->Tc = NULL;
376 ierr = PetscArrayzero(&f->quadrature, 1);CHKERRQ(ierr);
377 ierr = PetscArrayzero(&f->faceQuadrature, 1);CHKERRQ(ierr);
378 f->blockSize = 0;
379 f->numBlocks = 1;
380 f->batchSize = 0;
381 f->numBatches = 1;
382
383 *fem = f;
384 PetscFunctionReturn(0);
385 }
386
387 /*@
388 PetscFEGetSpatialDimension - Returns the spatial dimension of the element
389
390 Not collective
391
392 Input Parameter:
393 . fem - The PetscFE object
394
395 Output Parameter:
396 . dim - The spatial dimension
397
398 Level: intermediate
399
400 .seealso: PetscFECreate()
401 @*/
PetscFEGetSpatialDimension(PetscFE fem,PetscInt * dim)402 PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
403 {
404 DM dm;
405 PetscErrorCode ierr;
406
407 PetscFunctionBegin;
408 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
409 PetscValidPointer(dim, 2);
410 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr);
411 ierr = DMGetDimension(dm, dim);CHKERRQ(ierr);
412 PetscFunctionReturn(0);
413 }
414
415 /*@
416 PetscFESetNumComponents - Sets the number of components in the element
417
418 Not collective
419
420 Input Parameters:
421 + fem - The PetscFE object
422 - comp - The number of field components
423
424 Level: intermediate
425
426 .seealso: PetscFECreate()
427 @*/
PetscFESetNumComponents(PetscFE fem,PetscInt comp)428 PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
429 {
430 PetscFunctionBegin;
431 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
432 fem->numComponents = comp;
433 PetscFunctionReturn(0);
434 }
435
436 /*@
437 PetscFEGetNumComponents - Returns the number of components in the element
438
439 Not collective
440
441 Input Parameter:
442 . fem - The PetscFE object
443
444 Output Parameter:
445 . comp - The number of field components
446
447 Level: intermediate
448
449 .seealso: PetscFECreate()
450 @*/
PetscFEGetNumComponents(PetscFE fem,PetscInt * comp)451 PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
452 {
453 PetscFunctionBegin;
454 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
455 PetscValidPointer(comp, 2);
456 *comp = fem->numComponents;
457 PetscFunctionReturn(0);
458 }
459
460 /*@
461 PetscFESetTileSizes - Sets the tile sizes for evaluation
462
463 Not collective
464
465 Input Parameters:
466 + fem - The PetscFE object
467 . blockSize - The number of elements in a block
468 . numBlocks - The number of blocks in a batch
469 . batchSize - The number of elements in a batch
470 - numBatches - The number of batches in a chunk
471
472 Level: intermediate
473
474 .seealso: PetscFECreate()
475 @*/
PetscFESetTileSizes(PetscFE fem,PetscInt blockSize,PetscInt numBlocks,PetscInt batchSize,PetscInt numBatches)476 PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
477 {
478 PetscFunctionBegin;
479 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
480 fem->blockSize = blockSize;
481 fem->numBlocks = numBlocks;
482 fem->batchSize = batchSize;
483 fem->numBatches = numBatches;
484 PetscFunctionReturn(0);
485 }
486
487 /*@
488 PetscFEGetTileSizes - Returns the tile sizes for evaluation
489
490 Not collective
491
492 Input Parameter:
493 . fem - The PetscFE object
494
495 Output Parameters:
496 + blockSize - The number of elements in a block
497 . numBlocks - The number of blocks in a batch
498 . batchSize - The number of elements in a batch
499 - numBatches - The number of batches in a chunk
500
501 Level: intermediate
502
503 .seealso: PetscFECreate()
504 @*/
PetscFEGetTileSizes(PetscFE fem,PetscInt * blockSize,PetscInt * numBlocks,PetscInt * batchSize,PetscInt * numBatches)505 PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
506 {
507 PetscFunctionBegin;
508 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
509 if (blockSize) PetscValidPointer(blockSize, 2);
510 if (numBlocks) PetscValidPointer(numBlocks, 3);
511 if (batchSize) PetscValidPointer(batchSize, 4);
512 if (numBatches) PetscValidPointer(numBatches, 5);
513 if (blockSize) *blockSize = fem->blockSize;
514 if (numBlocks) *numBlocks = fem->numBlocks;
515 if (batchSize) *batchSize = fem->batchSize;
516 if (numBatches) *numBatches = fem->numBatches;
517 PetscFunctionReturn(0);
518 }
519
520 /*@
521 PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
522
523 Not collective
524
525 Input Parameter:
526 . fem - The PetscFE object
527
528 Output Parameter:
529 . sp - The PetscSpace object
530
531 Level: intermediate
532
533 .seealso: PetscFECreate()
534 @*/
PetscFEGetBasisSpace(PetscFE fem,PetscSpace * sp)535 PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
536 {
537 PetscFunctionBegin;
538 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
539 PetscValidPointer(sp, 2);
540 *sp = fem->basisSpace;
541 PetscFunctionReturn(0);
542 }
543
544 /*@
545 PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
546
547 Not collective
548
549 Input Parameters:
550 + fem - The PetscFE object
551 - sp - The PetscSpace object
552
553 Level: intermediate
554
555 .seealso: PetscFECreate()
556 @*/
PetscFESetBasisSpace(PetscFE fem,PetscSpace sp)557 PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
558 {
559 PetscErrorCode ierr;
560
561 PetscFunctionBegin;
562 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
563 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 2);
564 ierr = PetscSpaceDestroy(&fem->basisSpace);CHKERRQ(ierr);
565 fem->basisSpace = sp;
566 ierr = PetscObjectReference((PetscObject) fem->basisSpace);CHKERRQ(ierr);
567 PetscFunctionReturn(0);
568 }
569
570 /*@
571 PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
572
573 Not collective
574
575 Input Parameter:
576 . fem - The PetscFE object
577
578 Output Parameter:
579 . sp - The PetscDualSpace object
580
581 Level: intermediate
582
583 .seealso: PetscFECreate()
584 @*/
PetscFEGetDualSpace(PetscFE fem,PetscDualSpace * sp)585 PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
586 {
587 PetscFunctionBegin;
588 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
589 PetscValidPointer(sp, 2);
590 *sp = fem->dualSpace;
591 PetscFunctionReturn(0);
592 }
593
594 /*@
595 PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
596
597 Not collective
598
599 Input Parameters:
600 + fem - The PetscFE object
601 - sp - The PetscDualSpace object
602
603 Level: intermediate
604
605 .seealso: PetscFECreate()
606 @*/
PetscFESetDualSpace(PetscFE fem,PetscDualSpace sp)607 PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
608 {
609 PetscErrorCode ierr;
610
611 PetscFunctionBegin;
612 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
613 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 2);
614 ierr = PetscDualSpaceDestroy(&fem->dualSpace);CHKERRQ(ierr);
615 fem->dualSpace = sp;
616 ierr = PetscObjectReference((PetscObject) fem->dualSpace);CHKERRQ(ierr);
617 PetscFunctionReturn(0);
618 }
619
620 /*@
621 PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products
622
623 Not collective
624
625 Input Parameter:
626 . fem - The PetscFE object
627
628 Output Parameter:
629 . q - The PetscQuadrature object
630
631 Level: intermediate
632
633 .seealso: PetscFECreate()
634 @*/
PetscFEGetQuadrature(PetscFE fem,PetscQuadrature * q)635 PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
636 {
637 PetscFunctionBegin;
638 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
639 PetscValidPointer(q, 2);
640 *q = fem->quadrature;
641 PetscFunctionReturn(0);
642 }
643
644 /*@
645 PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products
646
647 Not collective
648
649 Input Parameters:
650 + fem - The PetscFE object
651 - q - The PetscQuadrature object
652
653 Level: intermediate
654
655 .seealso: PetscFECreate()
656 @*/
PetscFESetQuadrature(PetscFE fem,PetscQuadrature q)657 PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
658 {
659 PetscInt Nc, qNc;
660 PetscErrorCode ierr;
661
662 PetscFunctionBegin;
663 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
664 if (q == fem->quadrature) PetscFunctionReturn(0);
665 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr);
666 ierr = PetscQuadratureGetNumComponents(q, &qNc);CHKERRQ(ierr);
667 if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
668 ierr = PetscTabulationDestroy(&fem->T);CHKERRQ(ierr);
669 ierr = PetscTabulationDestroy(&fem->Tc);CHKERRQ(ierr);
670 ierr = PetscObjectReference((PetscObject) q);CHKERRQ(ierr);
671 ierr = PetscQuadratureDestroy(&fem->quadrature);CHKERRQ(ierr);
672 fem->quadrature = q;
673 PetscFunctionReturn(0);
674 }
675
676 /*@
677 PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces
678
679 Not collective
680
681 Input Parameter:
682 . fem - The PetscFE object
683
684 Output Parameter:
685 . q - The PetscQuadrature object
686
687 Level: intermediate
688
689 .seealso: PetscFECreate()
690 @*/
PetscFEGetFaceQuadrature(PetscFE fem,PetscQuadrature * q)691 PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q)
692 {
693 PetscFunctionBegin;
694 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
695 PetscValidPointer(q, 2);
696 *q = fem->faceQuadrature;
697 PetscFunctionReturn(0);
698 }
699
700 /*@
701 PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces
702
703 Not collective
704
705 Input Parameters:
706 + fem - The PetscFE object
707 - q - The PetscQuadrature object
708
709 Level: intermediate
710
711 .seealso: PetscFECreate()
712 @*/
PetscFESetFaceQuadrature(PetscFE fem,PetscQuadrature q)713 PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q)
714 {
715 PetscInt Nc, qNc;
716 PetscErrorCode ierr;
717
718 PetscFunctionBegin;
719 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
720 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr);
721 ierr = PetscQuadratureGetNumComponents(q, &qNc);CHKERRQ(ierr);
722 if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc);
723 ierr = PetscTabulationDestroy(&fem->Tf);CHKERRQ(ierr);
724 ierr = PetscQuadratureDestroy(&fem->faceQuadrature);CHKERRQ(ierr);
725 fem->faceQuadrature = q;
726 ierr = PetscObjectReference((PetscObject) q);CHKERRQ(ierr);
727 PetscFunctionReturn(0);
728 }
729
730 /*@
731 PetscFECopyQuadrature - Copy both volumetric and surface quadrature
732
733 Not collective
734
735 Input Parameters:
736 + sfe - The PetscFE source for the quadratures
737 - tfe - The PetscFE target for the quadratures
738
739 Level: intermediate
740
741 .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature()
742 @*/
PetscFECopyQuadrature(PetscFE sfe,PetscFE tfe)743 PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe)
744 {
745 PetscQuadrature q;
746 PetscErrorCode ierr;
747
748 PetscFunctionBegin;
749 PetscValidHeaderSpecific(sfe, PETSCFE_CLASSID, 1);
750 PetscValidHeaderSpecific(tfe, PETSCFE_CLASSID, 2);
751 ierr = PetscFEGetQuadrature(sfe, &q);CHKERRQ(ierr);
752 ierr = PetscFESetQuadrature(tfe, q);CHKERRQ(ierr);
753 ierr = PetscFEGetFaceQuadrature(sfe, &q);CHKERRQ(ierr);
754 ierr = PetscFESetFaceQuadrature(tfe, q);CHKERRQ(ierr);
755 PetscFunctionReturn(0);
756 }
757
758 /*@C
759 PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension
760
761 Not collective
762
763 Input Parameter:
764 . fem - The PetscFE object
765
766 Output Parameter:
767 . numDof - Array with the number of dofs per dimension
768
769 Level: intermediate
770
771 .seealso: PetscFECreate()
772 @*/
PetscFEGetNumDof(PetscFE fem,const PetscInt ** numDof)773 PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
774 {
775 PetscErrorCode ierr;
776
777 PetscFunctionBegin;
778 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
779 PetscValidPointer(numDof, 2);
780 ierr = PetscDualSpaceGetNumDof(fem->dualSpace, numDof);CHKERRQ(ierr);
781 PetscFunctionReturn(0);
782 }
783
784 /*@C
785 PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell
786
787 Not collective
788
789 Input Parameter:
790 . fem - The PetscFE object
791
792 Output Parameter:
793 . T - The basis function values and derivatives at quadrature points
794
795 Note:
796 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
797 $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
798 $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
799
800 Level: intermediate
801
802 .seealso: PetscFECreateTabulation(), PetscTabulationDestroy()
803 @*/
PetscFEGetCellTabulation(PetscFE fem,PetscTabulation * T)804 PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscTabulation *T)
805 {
806 PetscInt npoints;
807 const PetscReal *points;
808 PetscErrorCode ierr;
809
810 PetscFunctionBegin;
811 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
812 PetscValidPointer(T, 2);
813 ierr = PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr);
814 if (!fem->T) {ierr = PetscFECreateTabulation(fem, 1, npoints, points, 1, &fem->T);CHKERRQ(ierr);}
815 *T = fem->T;
816 PetscFunctionReturn(0);
817 }
818
819 /*@C
820 PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell
821
822 Not collective
823
824 Input Parameter:
825 . fem - The PetscFE object
826
827 Output Parameters:
828 . Tf - The basis function values and derviatives at face quadrature points
829
830 Note:
831 $ T->T[0] = Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c
832 $ T->T[1] = Df[(((f*Nq + q)*pdim + i)*Nc + c)*dim + d] is the derivative value at point f,q for basis function i, component c, in direction d
833 $ T->T[2] = Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f,q for basis function i, component c, in directions d and e
834
835 Level: intermediate
836
837 .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
838 @*/
PetscFEGetFaceTabulation(PetscFE fem,PetscTabulation * Tf)839 PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscTabulation *Tf)
840 {
841 PetscErrorCode ierr;
842
843 PetscFunctionBegin;
844 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
845 PetscValidPointer(Tf, 2);
846 if (!fem->Tf) {
847 const PetscReal xi0[3] = {-1., -1., -1.};
848 PetscReal v0[3], J[9], detJ;
849 PetscQuadrature fq;
850 PetscDualSpace sp;
851 DM dm;
852 const PetscInt *faces;
853 PetscInt dim, numFaces, f, npoints, q;
854 const PetscReal *points;
855 PetscReal *facePoints;
856
857 ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr);
858 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr);
859 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
860 ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr);
861 ierr = DMPlexGetCone(dm, 0, &faces);CHKERRQ(ierr);
862 ierr = PetscFEGetFaceQuadrature(fem, &fq);CHKERRQ(ierr);
863 if (fq) {
864 ierr = PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr);
865 ierr = PetscMalloc1(numFaces*npoints*dim, &facePoints);CHKERRQ(ierr);
866 for (f = 0; f < numFaces; ++f) {
867 ierr = DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);CHKERRQ(ierr);
868 for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]);
869 }
870 ierr = PetscFECreateTabulation(fem, numFaces, npoints, facePoints, 1, &fem->Tf);CHKERRQ(ierr);
871 ierr = PetscFree(facePoints);CHKERRQ(ierr);
872 }
873 }
874 *Tf = fem->Tf;
875 PetscFunctionReturn(0);
876 }
877
878 /*@C
879 PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points
880
881 Not collective
882
883 Input Parameter:
884 . fem - The PetscFE object
885
886 Output Parameters:
887 . Tc - The basis function values at face centroid points
888
889 Note:
890 $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c
891
892 Level: intermediate
893
894 .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy()
895 @*/
PetscFEGetFaceCentroidTabulation(PetscFE fem,PetscTabulation * Tc)896 PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc)
897 {
898 PetscErrorCode ierr;
899
900 PetscFunctionBegin;
901 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
902 PetscValidPointer(Tc, 2);
903 if (!fem->Tc) {
904 PetscDualSpace sp;
905 DM dm;
906 const PetscInt *cone;
907 PetscReal *centroids;
908 PetscInt dim, numFaces, f;
909
910 ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr);
911 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr);
912 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
913 ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr);
914 ierr = DMPlexGetCone(dm, 0, &cone);CHKERRQ(ierr);
915 ierr = PetscMalloc1(numFaces*dim, ¢roids);CHKERRQ(ierr);
916 for (f = 0; f < numFaces; ++f) {ierr = DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);CHKERRQ(ierr);}
917 ierr = PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);CHKERRQ(ierr);
918 ierr = PetscFree(centroids);CHKERRQ(ierr);
919 }
920 *Tc = fem->Tc;
921 PetscFunctionReturn(0);
922 }
923
924 /*@C
925 PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
926
927 Not collective
928
929 Input Parameters:
930 + fem - The PetscFE object
931 . nrepl - The number of replicas
932 . npoints - The number of tabulation points in a replica
933 . points - The tabulation point coordinates
934 - K - The number of derivatives calculated
935
936 Output Parameter:
937 . T - The basis function values and derivatives at tabulation points
938
939 Note:
940 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
941 $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
942 $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
943
944 Level: intermediate
945
946 .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
947 @*/
PetscFECreateTabulation(PetscFE fem,PetscInt nrepl,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation * T)948 PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
949 {
950 DM dm;
951 PetscDualSpace Q;
952 PetscInt Nb; /* Dimension of FE space P */
953 PetscInt Nc; /* Field components */
954 PetscInt cdim; /* Reference coordinate dimension */
955 PetscInt k;
956 PetscErrorCode ierr;
957
958 PetscFunctionBegin;
959 if (!npoints || !fem->dualSpace || K < 0) {
960 *T = NULL;
961 PetscFunctionReturn(0);
962 }
963 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
964 PetscValidPointer(points, 4);
965 PetscValidPointer(T, 6);
966 ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr);
967 ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr);
968 ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr);
969 ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr);
970 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr);
971 ierr = PetscMalloc1(1, T);CHKERRQ(ierr);
972 (*T)->K = !cdim ? 0 : K;
973 (*T)->Nr = nrepl;
974 (*T)->Np = npoints;
975 (*T)->Nb = Nb;
976 (*T)->Nc = Nc;
977 (*T)->cdim = cdim;
978 ierr = PetscMalloc1((*T)->K+1, &(*T)->T);CHKERRQ(ierr);
979 for (k = 0; k <= (*T)->K; ++k) {
980 ierr = PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);CHKERRQ(ierr);
981 }
982 ierr = (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);CHKERRQ(ierr);
983 PetscFunctionReturn(0);
984 }
985
986 /*@C
987 PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
988
989 Not collective
990
991 Input Parameters:
992 + fem - The PetscFE object
993 . npoints - The number of tabulation points
994 . points - The tabulation point coordinates
995 . K - The number of derivatives calculated
996 - T - An existing tabulation object with enough allocated space
997
998 Output Parameter:
999 . T - The basis function values and derivatives at tabulation points
1000
1001 Note:
1002 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1003 $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1004 $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
1005
1006 Level: intermediate
1007
1008 .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy()
1009 @*/
PetscFEComputeTabulation(PetscFE fem,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation T)1010 PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
1011 {
1012 PetscErrorCode ierr;
1013
1014 PetscFunctionBeginHot;
1015 if (!npoints || !fem->dualSpace || K < 0) PetscFunctionReturn(0);
1016 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
1017 PetscValidPointer(points, 3);
1018 PetscValidPointer(T, 5);
1019 if (PetscDefined(USE_DEBUG)) {
1020 DM dm;
1021 PetscDualSpace Q;
1022 PetscInt Nb; /* Dimension of FE space P */
1023 PetscInt Nc; /* Field components */
1024 PetscInt cdim; /* Reference coordinate dimension */
1025
1026 ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr);
1027 ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr);
1028 ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr);
1029 ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr);
1030 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr);
1031 if (T->K != (!cdim ? 0 : K)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %D must match requested K %D", T->K, !cdim ? 0 : K);
1032 if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb);
1033 if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc);
1034 if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim);
1035 }
1036 T->Nr = 1;
1037 T->Np = npoints;
1038 ierr = (*fem->ops->createtabulation)(fem, npoints, points, K, T);CHKERRQ(ierr);
1039 PetscFunctionReturn(0);
1040 }
1041
1042 /*@C
1043 PetscTabulationDestroy - Frees memory from the associated tabulation.
1044
1045 Not collective
1046
1047 Input Parameter:
1048 . T - The tabulation
1049
1050 Level: intermediate
1051
1052 .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation()
1053 @*/
PetscTabulationDestroy(PetscTabulation * T)1054 PetscErrorCode PetscTabulationDestroy(PetscTabulation *T)
1055 {
1056 PetscInt k;
1057 PetscErrorCode ierr;
1058
1059 PetscFunctionBegin;
1060 PetscValidPointer(T, 1);
1061 if (!T || !(*T)) PetscFunctionReturn(0);
1062 for (k = 0; k <= (*T)->K; ++k) {ierr = PetscFree((*T)->T[k]);CHKERRQ(ierr);}
1063 ierr = PetscFree((*T)->T);CHKERRQ(ierr);
1064 ierr = PetscFree(*T);CHKERRQ(ierr);
1065 *T = NULL;
1066 PetscFunctionReturn(0);
1067 }
1068
PetscFECreatePointTrace(PetscFE fe,PetscInt refPoint,PetscFE * trFE)1069 PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
1070 {
1071 PetscSpace bsp, bsubsp;
1072 PetscDualSpace dsp, dsubsp;
1073 PetscInt dim, depth, numComp, i, j, coneSize, order;
1074 PetscFEType type;
1075 DM dm;
1076 DMLabel label;
1077 PetscReal *xi, *v, *J, detJ;
1078 const char *name;
1079 PetscQuadrature origin, fullQuad, subQuad;
1080 PetscErrorCode ierr;
1081
1082 PetscFunctionBegin;
1083 PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1);
1084 PetscValidPointer(trFE,3);
1085 ierr = PetscFEGetBasisSpace(fe,&bsp);CHKERRQ(ierr);
1086 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr);
1087 ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr);
1088 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1089 ierr = DMPlexGetDepthLabel(dm,&label);CHKERRQ(ierr);
1090 ierr = DMLabelGetValue(label,refPoint,&depth);CHKERRQ(ierr);
1091 ierr = PetscCalloc1(depth,&xi);CHKERRQ(ierr);
1092 ierr = PetscMalloc1(dim,&v);CHKERRQ(ierr);
1093 ierr = PetscMalloc1(dim*dim,&J);CHKERRQ(ierr);
1094 for (i = 0; i < depth; i++) xi[i] = 0.;
1095 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,&origin);CHKERRQ(ierr);
1096 ierr = PetscQuadratureSetData(origin,depth,0,1,xi,NULL);CHKERRQ(ierr);
1097 ierr = DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);CHKERRQ(ierr);
1098 /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */
1099 for (i = 1; i < dim; i++) {
1100 for (j = 0; j < depth; j++) {
1101 J[i * depth + j] = J[i * dim + j];
1102 }
1103 }
1104 ierr = PetscQuadratureDestroy(&origin);CHKERRQ(ierr);
1105 ierr = PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);CHKERRQ(ierr);
1106 ierr = PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);CHKERRQ(ierr);
1107 ierr = PetscSpaceSetUp(bsubsp);CHKERRQ(ierr);
1108 ierr = PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);CHKERRQ(ierr);
1109 ierr = PetscFEGetType(fe,&type);CHKERRQ(ierr);
1110 ierr = PetscFESetType(*trFE,type);CHKERRQ(ierr);
1111 ierr = PetscFEGetNumComponents(fe,&numComp);CHKERRQ(ierr);
1112 ierr = PetscFESetNumComponents(*trFE,numComp);CHKERRQ(ierr);
1113 ierr = PetscFESetBasisSpace(*trFE,bsubsp);CHKERRQ(ierr);
1114 ierr = PetscFESetDualSpace(*trFE,dsubsp);CHKERRQ(ierr);
1115 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr);
1116 if (name) {ierr = PetscFESetName(*trFE, name);CHKERRQ(ierr);}
1117 ierr = PetscFEGetQuadrature(fe,&fullQuad);CHKERRQ(ierr);
1118 ierr = PetscQuadratureGetOrder(fullQuad,&order);CHKERRQ(ierr);
1119 ierr = DMPlexGetConeSize(dm,refPoint,&coneSize);CHKERRQ(ierr);
1120 if (coneSize == 2 * depth) {
1121 ierr = PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr);
1122 } else {
1123 ierr = PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr);
1124 }
1125 ierr = PetscFESetQuadrature(*trFE,subQuad);CHKERRQ(ierr);
1126 ierr = PetscFESetUp(*trFE);CHKERRQ(ierr);
1127 ierr = PetscQuadratureDestroy(&subQuad);CHKERRQ(ierr);
1128 ierr = PetscSpaceDestroy(&bsubsp);CHKERRQ(ierr);
1129 PetscFunctionReturn(0);
1130 }
1131
PetscFECreateHeightTrace(PetscFE fe,PetscInt height,PetscFE * trFE)1132 PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE)
1133 {
1134 PetscInt hStart, hEnd;
1135 PetscDualSpace dsp;
1136 DM dm;
1137 PetscErrorCode ierr;
1138
1139 PetscFunctionBegin;
1140 PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1);
1141 PetscValidPointer(trFE,3);
1142 *trFE = NULL;
1143 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr);
1144 ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr);
1145 ierr = DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);CHKERRQ(ierr);
1146 if (hEnd <= hStart) PetscFunctionReturn(0);
1147 ierr = PetscFECreatePointTrace(fe,hStart,trFE);CHKERRQ(ierr);
1148 PetscFunctionReturn(0);
1149 }
1150
1151
1152 /*@
1153 PetscFEGetDimension - Get the dimension of the finite element space on a cell
1154
1155 Not collective
1156
1157 Input Parameter:
1158 . fe - The PetscFE
1159
1160 Output Parameter:
1161 . dim - The dimension
1162
1163 Level: intermediate
1164
1165 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
1166 @*/
PetscFEGetDimension(PetscFE fem,PetscInt * dim)1167 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
1168 {
1169 PetscErrorCode ierr;
1170
1171 PetscFunctionBegin;
1172 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
1173 PetscValidPointer(dim, 2);
1174 if (fem->ops->getdimension) {ierr = (*fem->ops->getdimension)(fem, dim);CHKERRQ(ierr);}
1175 PetscFunctionReturn(0);
1176 }
1177
1178 /*@C
1179 PetscFEPushforward - Map the reference element function to real space
1180
1181 Input Parameters:
1182 + fe - The PetscFE
1183 . fegeom - The cell geometry
1184 . Nv - The number of function values
1185 - vals - The function values
1186
1187 Output Parameter:
1188 . vals - The transformed function values
1189
1190 Level: advanced
1191
1192 Note: This just forwards the call onto PetscDualSpacePushforward().
1193
1194 Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1195
1196 .seealso: PetscDualSpacePushforward()
1197 @*/
PetscFEPushforward(PetscFE fe,PetscFEGeom * fegeom,PetscInt Nv,PetscScalar vals[])1198 PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1199 {
1200 PetscErrorCode ierr;
1201
1202 PetscFunctionBeginHot;
1203 ierr = PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr);
1204 PetscFunctionReturn(0);
1205 }
1206
1207 /*@C
1208 PetscFEPushforwardGradient - Map the reference element function gradient to real space
1209
1210 Input Parameters:
1211 + fe - The PetscFE
1212 . fegeom - The cell geometry
1213 . Nv - The number of function gradient values
1214 - vals - The function gradient values
1215
1216 Output Parameter:
1217 . vals - The transformed function gradient values
1218
1219 Level: advanced
1220
1221 Note: This just forwards the call onto PetscDualSpacePushforwardGradient().
1222
1223 Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.
1224
1225 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward()
1226 @*/
PetscFEPushforwardGradient(PetscFE fe,PetscFEGeom * fegeom,PetscInt Nv,PetscScalar vals[])1227 PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[])
1228 {
1229 PetscErrorCode ierr;
1230
1231 PetscFunctionBeginHot;
1232 ierr = PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr);
1233 PetscFunctionReturn(0);
1234 }
1235
1236 /*
1237 Purpose: Compute element vector for chunk of elements
1238
1239 Input:
1240 Sizes:
1241 Ne: number of elements
1242 Nf: number of fields
1243 PetscFE
1244 dim: spatial dimension
1245 Nb: number of basis functions
1246 Nc: number of field components
1247 PetscQuadrature
1248 Nq: number of quadrature points
1249
1250 Geometry:
1251 PetscFEGeom[Ne] possibly *Nq
1252 PetscReal v0s[dim]
1253 PetscReal n[dim]
1254 PetscReal jacobians[dim*dim]
1255 PetscReal jacobianInverses[dim*dim]
1256 PetscReal jacobianDeterminants
1257 FEM:
1258 PetscFE
1259 PetscQuadrature
1260 PetscReal quadPoints[Nq*dim]
1261 PetscReal quadWeights[Nq]
1262 PetscReal basis[Nq*Nb*Nc]
1263 PetscReal basisDer[Nq*Nb*Nc*dim]
1264 PetscScalar coefficients[Ne*Nb*Nc]
1265 PetscScalar elemVec[Ne*Nb*Nc]
1266
1267 Problem:
1268 PetscInt f: the active field
1269 f0, f1
1270
1271 Work Space:
1272 PetscFE
1273 PetscScalar f0[Nq*dim];
1274 PetscScalar f1[Nq*dim*dim];
1275 PetscScalar u[Nc];
1276 PetscScalar gradU[Nc*dim];
1277 PetscReal x[dim];
1278 PetscScalar realSpaceDer[dim];
1279
1280 Purpose: Compute element vector for N_cb batches of elements
1281
1282 Input:
1283 Sizes:
1284 N_cb: Number of serial cell batches
1285
1286 Geometry:
1287 PetscReal v0s[Ne*dim]
1288 PetscReal jacobians[Ne*dim*dim] possibly *Nq
1289 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
1290 PetscReal jacobianDeterminants[Ne] possibly *Nq
1291 FEM:
1292 static PetscReal quadPoints[Nq*dim]
1293 static PetscReal quadWeights[Nq]
1294 static PetscReal basis[Nq*Nb*Nc]
1295 static PetscReal basisDer[Nq*Nb*Nc*dim]
1296 PetscScalar coefficients[Ne*Nb*Nc]
1297 PetscScalar elemVec[Ne*Nb*Nc]
1298
1299 ex62.c:
1300 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
1301 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
1302 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
1303 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
1304
1305 ex52.c:
1306 PetscErrorCode IntegrateLaplacianBatchCPU(PetscInt Ne, PetscInt Nb, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user)
1307 PetscErrorCode IntegrateElasticityBatchCPU(PetscInt Ne, PetscInt Nb, PetscInt Ncomp, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user)
1308
1309 ex52_integrateElement.cu
1310 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
1311
1312 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
1313 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1314 PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1315
1316 ex52_integrateElementOpenCL.c:
1317 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
1318 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
1319 PetscLogEvent event, PetscInt debug, PetscInt pde_op)
1320
1321 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
1322 */
1323
1324 /*@C
1325 PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
1326
1327 Not collective
1328
1329 Input Parameters:
1330 + prob - The PetscDS specifying the discretizations and continuum functions
1331 . field - The field being integrated
1332 . Ne - The number of elements in the chunk
1333 . cgeom - The cell geometry for each cell in the chunk
1334 . coefficients - The array of FEM basis coefficients for the elements
1335 . probAux - The PetscDS specifying the auxiliary discretizations
1336 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1337
1338 Output Parameter:
1339 . integral - the integral for this field
1340
1341 Level: intermediate
1342
1343 .seealso: PetscFEIntegrateResidual()
1344 @*/
PetscFEIntegrate(PetscDS prob,PetscInt field,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscScalar integral[])1345 PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1346 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1347 {
1348 PetscFE fe;
1349 PetscErrorCode ierr;
1350
1351 PetscFunctionBegin;
1352 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1353 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr);
1354 if (fe->ops->integrate) {ierr = (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);}
1355 PetscFunctionReturn(0);
1356 }
1357
1358 /*@C
1359 PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration
1360
1361 Not collective
1362
1363 Input Parameters:
1364 + prob - The PetscDS specifying the discretizations and continuum functions
1365 . field - The field being integrated
1366 . obj_func - The function to be integrated
1367 . Ne - The number of elements in the chunk
1368 . fgeom - The face geometry for each face in the chunk
1369 . coefficients - The array of FEM basis coefficients for the elements
1370 . probAux - The PetscDS specifying the auxiliary discretizations
1371 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1372
1373 Output Parameter:
1374 . integral - the integral for this field
1375
1376 Level: intermediate
1377
1378 .seealso: PetscFEIntegrateResidual()
1379 @*/
PetscFEIntegrateBd(PetscDS prob,PetscInt field,void (* obj_func)(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[]),PetscInt Ne,PetscFEGeom * geom,const PetscScalar coefficients[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscScalar integral[])1380 PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field,
1381 void (*obj_func)(PetscInt, PetscInt, PetscInt,
1382 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1383 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1384 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1385 PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
1386 {
1387 PetscFE fe;
1388 PetscErrorCode ierr;
1389
1390 PetscFunctionBegin;
1391 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1392 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr);
1393 if (fe->ops->integratebd) {ierr = (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);}
1394 PetscFunctionReturn(0);
1395 }
1396
1397 /*@C
1398 PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
1399
1400 Not collective
1401
1402 Input Parameters:
1403 + prob - The PetscDS specifying the discretizations and continuum functions
1404 . field - The field being integrated
1405 . Ne - The number of elements in the chunk
1406 . cgeom - The cell geometry for each cell in the chunk
1407 . coefficients - The array of FEM basis coefficients for the elements
1408 . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1409 . probAux - The PetscDS specifying the auxiliary discretizations
1410 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1411 - t - The time
1412
1413 Output Parameter:
1414 . elemVec - the element residual vectors from each element
1415
1416 Note:
1417 $ Loop over batch of elements (e):
1418 $ Loop over quadrature points (q):
1419 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
1420 $ Call f_0 and f_1
1421 $ Loop over element vector entries (f,fc --> i):
1422 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
1423
1424 Level: intermediate
1425
1426 .seealso: PetscFEIntegrateResidual()
1427 @*/
PetscFEIntegrateResidual(PetscDS prob,PetscInt field,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])1428 PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
1429 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1430 {
1431 PetscFE fe;
1432 PetscErrorCode ierr;
1433
1434 PetscFunctionBegin;
1435 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1436 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr);
1437 if (fe->ops->integrateresidual) {ierr = (*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);}
1438 PetscFunctionReturn(0);
1439 }
1440
1441 /*@C
1442 PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
1443
1444 Not collective
1445
1446 Input Parameters:
1447 + prob - The PetscDS specifying the discretizations and continuum functions
1448 . field - The field being integrated
1449 . Ne - The number of elements in the chunk
1450 . fgeom - The face geometry for each cell in the chunk
1451 . coefficients - The array of FEM basis coefficients for the elements
1452 . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1453 . probAux - The PetscDS specifying the auxiliary discretizations
1454 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1455 - t - The time
1456
1457 Output Parameter:
1458 . elemVec - the element residual vectors from each element
1459
1460 Level: intermediate
1461
1462 .seealso: PetscFEIntegrateResidual()
1463 @*/
PetscFEIntegrateBdResidual(PetscDS prob,PetscInt field,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])1464 PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1465 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1466 {
1467 PetscFE fe;
1468 PetscErrorCode ierr;
1469
1470 PetscFunctionBegin;
1471 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1472 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr);
1473 if (fe->ops->integratebdresidual) {ierr = (*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);}
1474 PetscFunctionReturn(0);
1475 }
1476
1477 /*@C
1478 PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration
1479
1480 Not collective
1481
1482 Input Parameters:
1483 + prob - The PetscDS specifying the discretizations and continuum functions
1484 . field - The field being integrated
1485 . Ne - The number of elements in the chunk
1486 . fgeom - The face geometry for each cell in the chunk
1487 . coefficients - The array of FEM basis coefficients for the elements
1488 . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1489 . probAux - The PetscDS specifying the auxiliary discretizations
1490 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1491 - t - The time
1492
1493 Output Parameter
1494 . elemVec - the element residual vectors from each element
1495
1496 Level: developer
1497
1498 .seealso: PetscFEIntegrateResidual()
1499 @*/
PetscFEIntegrateHybridResidual(PetscDS prob,PetscInt field,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])1500 PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
1501 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
1502 {
1503 PetscFE fe;
1504 PetscErrorCode ierr;
1505
1506 PetscFunctionBegin;
1507 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1508 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr);
1509 if (fe->ops->integratehybridresidual) {ierr = (*fe->ops->integratehybridresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);}
1510 PetscFunctionReturn(0);
1511 }
1512
1513 /*@C
1514 PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
1515
1516 Not collective
1517
1518 Input Parameters:
1519 + prob - The PetscDS specifying the discretizations and continuum functions
1520 . jtype - The type of matrix pointwise functions that should be used
1521 . fieldI - The test field being integrated
1522 . fieldJ - The basis field being integrated
1523 . Ne - The number of elements in the chunk
1524 . cgeom - The cell geometry for each cell in the chunk
1525 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1526 . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1527 . probAux - The PetscDS specifying the auxiliary discretizations
1528 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1529 . t - The time
1530 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1531
1532 Output Parameter:
1533 . elemMat - the element matrices for the Jacobian from each element
1534
1535 Note:
1536 $ Loop over batch of elements (e):
1537 $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1538 $ Loop over quadrature points (q):
1539 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1540 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1541 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1542 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1543 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1544 Level: intermediate
1545
1546 .seealso: PetscFEIntegrateResidual()
1547 @*/
PetscFEIntegrateJacobian(PetscDS prob,PetscFEJacobianType jtype,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])1548 PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
1549 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1550 {
1551 PetscFE fe;
1552 PetscErrorCode ierr;
1553
1554 PetscFunctionBegin;
1555 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1556 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1557 if (fe->ops->integratejacobian) {ierr = (*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);}
1558 PetscFunctionReturn(0);
1559 }
1560
1561 /*@C
1562 PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
1563
1564 Not collective
1565
1566 Input Parameters:
1567 + prob - The PetscDS specifying the discretizations and continuum functions
1568 . fieldI - The test field being integrated
1569 . fieldJ - The basis field being integrated
1570 . Ne - The number of elements in the chunk
1571 . fgeom - The face geometry for each cell in the chunk
1572 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1573 . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1574 . probAux - The PetscDS specifying the auxiliary discretizations
1575 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1576 . t - The time
1577 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1578
1579 Output Parameter:
1580 . elemMat - the element matrices for the Jacobian from each element
1581
1582 Note:
1583 $ Loop over batch of elements (e):
1584 $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1585 $ Loop over quadrature points (q):
1586 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1587 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1588 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1589 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1590 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1591 Level: intermediate
1592
1593 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1594 @*/
PetscFEIntegrateBdJacobian(PetscDS prob,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])1595 PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1596 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1597 {
1598 PetscFE fe;
1599 PetscErrorCode ierr;
1600
1601 PetscFunctionBegin;
1602 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1603 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1604 if (fe->ops->integratebdjacobian) {ierr = (*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);}
1605 PetscFunctionReturn(0);
1606 }
1607
1608 /*@C
1609 PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration
1610
1611 Not collective
1612
1613 Input Parameters:
1614 + prob - The PetscDS specifying the discretizations and continuum functions
1615 . jtype - The type of matrix pointwise functions that should be used
1616 . fieldI - The test field being integrated
1617 . fieldJ - The basis field being integrated
1618 . Ne - The number of elements in the chunk
1619 . fgeom - The face geometry for each cell in the chunk
1620 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
1621 . coefficients_t - The array of FEM basis time derivative coefficients for the elements
1622 . probAux - The PetscDS specifying the auxiliary discretizations
1623 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
1624 . t - The time
1625 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term)
1626
1627 Output Parameter
1628 . elemMat - the element matrices for the Jacobian from each element
1629
1630 Note:
1631 $ Loop over batch of elements (e):
1632 $ Loop over element matrix entries (f,fc,g,gc --> i,j):
1633 $ Loop over quadrature points (q):
1634 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
1635 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
1636 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1637 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
1638 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
1639 Level: developer
1640
1641 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual()
1642 @*/
PetscFEIntegrateHybridJacobian(PetscDS prob,PetscFEJacobianType jtype,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS probAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])1643 PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1644 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1645 {
1646 PetscFE fe;
1647 PetscErrorCode ierr;
1648
1649 PetscFunctionBegin;
1650 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1651 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1652 if (fe->ops->integratehybridjacobian) {ierr = (*fe->ops->integratehybridjacobian)(prob, jtype, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);}
1653 PetscFunctionReturn(0);
1654 }
1655
1656 /*@
1657 PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height
1658
1659 Input Parameters:
1660 + fe - The finite element space
1661 - height - The height of the Plex point
1662
1663 Output Parameter:
1664 . subfe - The subspace of this FE space
1665
1666 Note: For example, if we want the subspace of this space for a face, we would choose height = 1.
1667
1668 Level: advanced
1669
1670 .seealso: PetscFECreateDefault()
1671 @*/
PetscFEGetHeightSubspace(PetscFE fe,PetscInt height,PetscFE * subfe)1672 PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe)
1673 {
1674 PetscSpace P, subP;
1675 PetscDualSpace Q, subQ;
1676 PetscQuadrature subq;
1677 PetscFEType fetype;
1678 PetscInt dim, Nc;
1679 PetscErrorCode ierr;
1680
1681 PetscFunctionBegin;
1682 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1);
1683 PetscValidPointer(subfe, 3);
1684 if (height == 0) {
1685 *subfe = fe;
1686 PetscFunctionReturn(0);
1687 }
1688 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr);
1689 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr);
1690 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1691 ierr = PetscFEGetFaceQuadrature(fe, &subq);CHKERRQ(ierr);
1692 ierr = PetscDualSpaceGetDimension(Q, &dim);CHKERRQ(ierr);
1693 if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);
1694 if (!fe->subspaces) {ierr = PetscCalloc1(dim, &fe->subspaces);CHKERRQ(ierr);}
1695 if (height <= dim) {
1696 if (!fe->subspaces[height-1]) {
1697 PetscFE sub = NULL;
1698 const char *name;
1699
1700 ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr);
1701 ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr);
1702 if (subQ) {
1703 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr);
1704 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr);
1705 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr);
1706 ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr);
1707 ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr);
1708 ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr);
1709 ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr);
1710 ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr);
1711 ierr = PetscFESetUp(sub);CHKERRQ(ierr);
1712 ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr);
1713 }
1714 fe->subspaces[height-1] = sub;
1715 }
1716 *subfe = fe->subspaces[height-1];
1717 } else {
1718 *subfe = NULL;
1719 }
1720 PetscFunctionReturn(0);
1721 }
1722
1723 /*@
1724 PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
1725 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1726 sparsity). It is also used to create an interpolation between regularly refined meshes.
1727
1728 Collective on fem
1729
1730 Input Parameter:
1731 . fe - The initial PetscFE
1732
1733 Output Parameter:
1734 . feRef - The refined PetscFE
1735
1736 Level: advanced
1737
1738 .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1739 @*/
PetscFERefine(PetscFE fe,PetscFE * feRef)1740 PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
1741 {
1742 PetscSpace P, Pref;
1743 PetscDualSpace Q, Qref;
1744 DM K, Kref;
1745 PetscQuadrature q, qref;
1746 const PetscReal *v0, *jac;
1747 PetscInt numComp, numSubelements;
1748 PetscInt cStart, cEnd, c;
1749 PetscDualSpace *cellSpaces;
1750 PetscErrorCode ierr;
1751
1752 PetscFunctionBegin;
1753 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr);
1754 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr);
1755 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1756 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr);
1757 /* Create space */
1758 ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr);
1759 Pref = P;
1760 /* Create dual space */
1761 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr);
1762 ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr);
1763 ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr);
1764 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr);
1765 ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr);
1766 ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr);
1767 /* TODO: fix for non-uniform refinement */
1768 for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q;
1769 ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr);
1770 ierr = PetscFree(cellSpaces);CHKERRQ(ierr);
1771 ierr = DMDestroy(&Kref);CHKERRQ(ierr);
1772 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr);
1773 /* Create element */
1774 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr);
1775 ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr);
1776 ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr);
1777 ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr);
1778 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr);
1779 ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr);
1780 ierr = PetscFESetUp(*feRef);CHKERRQ(ierr);
1781 ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr);
1782 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr);
1783 /* Create quadrature */
1784 ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr);
1785 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr);
1786 ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr);
1787 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr);
1788 PetscFunctionReturn(0);
1789 }
1790
1791 /*@C
1792 PetscFECreateDefault - Create a PetscFE for basic FEM computation
1793
1794 Collective
1795
1796 Input Parameters:
1797 + comm - The MPI comm
1798 . dim - The spatial dimension
1799 . Nc - The number of components
1800 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1801 . prefix - The options prefix, or NULL
1802 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1803
1804 Output Parameter:
1805 . fem - The PetscFE object
1806
1807 Note:
1808 Each object is SetFromOption() during creation, so that the object may be customized from the command line.
1809
1810 Level: beginner
1811
1812 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1813 @*/
PetscFECreateDefault(MPI_Comm comm,PetscInt dim,PetscInt Nc,PetscBool isSimplex,const char prefix[],PetscInt qorder,PetscFE * fem)1814 PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
1815 {
1816 PetscQuadrature q, fq;
1817 DM K;
1818 PetscSpace P;
1819 PetscDualSpace Q;
1820 PetscInt order, quadPointsPerEdge;
1821 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1822 PetscErrorCode ierr;
1823
1824 PetscFunctionBegin;
1825 /* Create space */
1826 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr);
1827 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr);
1828 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr);
1829 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr);
1830 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr);
1831 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr);
1832 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr);
1833 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr);
1834 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr);
1835 /* Create dual space */
1836 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr);
1837 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr);
1838 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr);
1839 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr);
1840 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr);
1841 ierr = DMDestroy(&K);CHKERRQ(ierr);
1842 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr);
1843 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr);
1844 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr);
1845 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr);
1846 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr);
1847 /* Create element */
1848 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr);
1849 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr);
1850 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr);
1851 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr);
1852 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr);
1853 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr);
1854 ierr = PetscFESetUp(*fem);CHKERRQ(ierr);
1855 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr);
1856 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr);
1857 /* Create quadrature (with specified order if given) */
1858 qorder = qorder >= 0 ? qorder : order;
1859 ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr);
1860 ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr);
1861 ierr = PetscOptionsEnd();CHKERRQ(ierr);
1862 quadPointsPerEdge = PetscMax(qorder + 1,1);
1863 if (isSimplex) {
1864 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr);
1865 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr);
1866 } else {
1867 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr);
1868 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr);
1869 }
1870 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr);
1871 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr);
1872 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr);
1873 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr);
1874 PetscFunctionReturn(0);
1875 }
1876
1877 /*@
1878 PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k
1879
1880 Collective
1881
1882 Input Parameters:
1883 + comm - The MPI comm
1884 . dim - The spatial dimension
1885 . Nc - The number of components
1886 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
1887 . k - The degree k of the space
1888 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree
1889
1890 Output Parameter:
1891 . fem - The PetscFE object
1892
1893 Level: beginner
1894
1895 Notes:
1896 For simplices, this element is the space of maximum polynomial degree k, otherwise it is a tensor product of 1D polynomials, each with maximal degree k.
1897
1898 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1899 @*/
PetscFECreateLagrange(MPI_Comm comm,PetscInt dim,PetscInt Nc,PetscBool isSimplex,PetscInt k,PetscInt qorder,PetscFE * fem)1900 PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem)
1901 {
1902 PetscQuadrature q, fq;
1903 DM K;
1904 PetscSpace P;
1905 PetscDualSpace Q;
1906 PetscInt quadPointsPerEdge;
1907 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1908 char name[64];
1909 PetscErrorCode ierr;
1910
1911 PetscFunctionBegin;
1912 /* Create space */
1913 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr);
1914 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr);
1915 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr);
1916 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr);
1917 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr);
1918 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr);
1919 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr);
1920 /* Create dual space */
1921 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr);
1922 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr);
1923 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr);
1924 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr);
1925 ierr = DMDestroy(&K);CHKERRQ(ierr);
1926 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr);
1927 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr);
1928 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr);
1929 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr);
1930 /* Create finite element */
1931 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr);
1932 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr);
1933 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr);
1934 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr);
1935 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr);
1936 ierr = PetscFESetUp(*fem);CHKERRQ(ierr);
1937 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr);
1938 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr);
1939 /* Create quadrature (with specified order if given) */
1940 qorder = qorder >= 0 ? qorder : k;
1941 quadPointsPerEdge = PetscMax(qorder + 1,1);
1942 if (isSimplex) {
1943 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr);
1944 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr);
1945 } else {
1946 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr);
1947 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr);
1948 }
1949 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr);
1950 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr);
1951 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr);
1952 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr);
1953 /* Set finite element name */
1954 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr);
1955 ierr = PetscFESetName(*fem, name);CHKERRQ(ierr);
1956 PetscFunctionReturn(0);
1957 }
1958
1959 /*@C
1960 PetscFESetName - Names the FE and its subobjects
1961
1962 Not collective
1963
1964 Input Parameters:
1965 + fe - The PetscFE
1966 - name - The name
1967
1968 Level: intermediate
1969
1970 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
1971 @*/
PetscFESetName(PetscFE fe,const char name[])1972 PetscErrorCode PetscFESetName(PetscFE fe, const char name[])
1973 {
1974 PetscSpace P;
1975 PetscDualSpace Q;
1976 PetscErrorCode ierr;
1977
1978 PetscFunctionBegin;
1979 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr);
1980 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr);
1981 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr);
1982 ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr);
1983 ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr);
1984 PetscFunctionReturn(0);
1985 }
1986
PetscFEEvaluateFieldJets_Internal(PetscDS ds,PetscInt Nf,PetscInt r,PetscInt q,PetscTabulation T[],PetscFEGeom * fegeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscScalar u[],PetscScalar u_x[],PetscScalar u_t[])1987 PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
1988 {
1989 PetscInt dOffset = 0, fOffset = 0, f;
1990 PetscErrorCode ierr;
1991
1992 for (f = 0; f < Nf; ++f) {
1993 PetscFE fe;
1994 const PetscInt cdim = T[f]->cdim;
1995 const PetscInt Nq = T[f]->Np;
1996 const PetscInt Nbf = T[f]->Nb;
1997 const PetscInt Ncf = T[f]->Nc;
1998 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
1999 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2000 PetscInt b, c, d;
2001
2002 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr);
2003 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2004 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2005 for (b = 0; b < Nbf; ++b) {
2006 for (c = 0; c < Ncf; ++c) {
2007 const PetscInt cidx = b*Ncf+c;
2008
2009 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2010 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2011 }
2012 }
2013 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr);
2014 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr);
2015 if (u_t) {
2016 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2017 for (b = 0; b < Nbf; ++b) {
2018 for (c = 0; c < Ncf; ++c) {
2019 const PetscInt cidx = b*Ncf+c;
2020
2021 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2022 }
2023 }
2024 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr);
2025 }
2026 fOffset += Ncf;
2027 dOffset += Nbf;
2028 }
2029 return 0;
2030 }
2031
PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds,PetscInt Nf,PetscInt r,PetscInt q,PetscTabulation T[],PetscFEGeom * fegeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscScalar u[],PetscScalar u_x[],PetscScalar u_t[])2032 PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
2033 {
2034 PetscInt dOffset = 0, fOffset = 0, g;
2035 PetscErrorCode ierr;
2036
2037 for (g = 0; g < 2*Nf-1; ++g) {
2038 if (!T[g/2]) continue;
2039 {
2040 PetscFE fe;
2041 const PetscInt f = g/2;
2042 const PetscInt cdim = T[f]->cdim;
2043 const PetscInt Nq = T[f]->Np;
2044 const PetscInt Nbf = T[f]->Nb;
2045 const PetscInt Ncf = T[f]->Nc;
2046 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf];
2047 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim];
2048 PetscInt b, c, d;
2049
2050 fe = (PetscFE) ds->disc[f];
2051 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
2052 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0;
2053 for (b = 0; b < Nbf; ++b) {
2054 for (c = 0; c < Ncf; ++c) {
2055 const PetscInt cidx = b*Ncf+c;
2056
2057 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
2058 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b];
2059 }
2060 }
2061 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr);
2062 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr);
2063 if (u_t) {
2064 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
2065 for (b = 0; b < Nbf; ++b) {
2066 for (c = 0; c < Ncf; ++c) {
2067 const PetscInt cidx = b*Ncf+c;
2068
2069 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
2070 }
2071 }
2072 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr);
2073 }
2074 fOffset += Ncf;
2075 dOffset += Nbf;
2076 }
2077 }
2078 return 0;
2079 }
2080
PetscFEEvaluateFaceFields_Internal(PetscDS prob,PetscInt field,PetscInt faceLoc,const PetscScalar coefficients[],PetscScalar u[])2081 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
2082 {
2083 PetscFE fe;
2084 PetscTabulation Tc;
2085 PetscInt b, c;
2086 PetscErrorCode ierr;
2087
2088 if (!prob) return 0;
2089 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr);
2090 ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr);
2091 {
2092 const PetscReal *faceBasis = Tc->T[0];
2093 const PetscInt Nb = Tc->Nb;
2094 const PetscInt Nc = Tc->Nc;
2095
2096 for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
2097 for (b = 0; b < Nb; ++b) {
2098 for (c = 0; c < Nc; ++c) {
2099 u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c];
2100 }
2101 }
2102 }
2103 return 0;
2104 }
2105
PetscFEUpdateElementVec_Internal(PetscFE fe,PetscTabulation T,PetscInt r,PetscScalar tmpBasis[],PetscScalar tmpBasisDer[],PetscFEGeom * fegeom,PetscScalar f0[],PetscScalar f1[],PetscScalar elemVec[])2106 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2107 {
2108 const PetscInt dE = T->cdim; /* fegeom->dimEmbed */
2109 const PetscInt Nq = T->Np;
2110 const PetscInt Nb = T->Nb;
2111 const PetscInt Nc = T->Nc;
2112 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc];
2113 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2114 PetscInt q, b, c, d;
2115 PetscErrorCode ierr;
2116
2117 for (b = 0; b < Nb; ++b) elemVec[b] = 0.0;
2118 for (q = 0; q < Nq; ++q) {
2119 for (b = 0; b < Nb; ++b) {
2120 for (c = 0; c < Nc; ++c) {
2121 const PetscInt bcidx = b*Nc+c;
2122
2123 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2124 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2125 }
2126 }
2127 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr);
2128 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr);
2129 for (b = 0; b < Nb; ++b) {
2130 for (c = 0; c < Nc; ++c) {
2131 const PetscInt bcidx = b*Nc+c;
2132 const PetscInt qcidx = q*Nc+c;
2133
2134 elemVec[b] += tmpBasis[bcidx]*f0[qcidx];
2135 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2136 }
2137 }
2138 }
2139 return(0);
2140 }
2141
PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe,PetscTabulation T,PetscInt r,PetscScalar tmpBasis[],PetscScalar tmpBasisDer[],PetscFEGeom * fegeom,PetscScalar f0[],PetscScalar f1[],PetscScalar elemVec[])2142 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
2143 {
2144 const PetscInt dE = T->cdim;
2145 const PetscInt Nq = T->Np;
2146 const PetscInt Nb = T->Nb;
2147 const PetscInt Nc = T->Nc;
2148 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc];
2149 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE];
2150 PetscInt q, b, c, d, s;
2151 PetscErrorCode ierr;
2152
2153 for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0;
2154 for (q = 0; q < Nq; ++q) {
2155 for (b = 0; b < Nb; ++b) {
2156 for (c = 0; c < Nc; ++c) {
2157 const PetscInt bcidx = b*Nc+c;
2158
2159 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx];
2160 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d];
2161 }
2162 }
2163 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr);
2164 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr);
2165 for (s = 0; s < 2; ++s) {
2166 for (b = 0; b < Nb; ++b) {
2167 for (c = 0; c < Nc; ++c) {
2168 const PetscInt bcidx = b*Nc+c;
2169 const PetscInt qcidx = (q*2+s)*Nc+c;
2170
2171 elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx];
2172 for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d];
2173 }
2174 }
2175 }
2176 }
2177 return(0);
2178 }
2179
PetscFEUpdateElementMat_Internal(PetscFE feI,PetscFE feJ,PetscInt r,PetscInt q,PetscTabulation TI,PetscScalar tmpBasisI[],PetscScalar tmpBasisDerI[],PetscTabulation TJ,PetscScalar tmpBasisJ[],PetscScalar tmpBasisDerJ[],PetscFEGeom * fegeom,const PetscScalar g0[],const PetscScalar g1[],const PetscScalar g2[],const PetscScalar g3[],PetscInt eOffset,PetscInt totDim,PetscInt offsetI,PetscInt offsetJ,PetscScalar elemMat[])2180 PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2181 {
2182 const PetscInt dE = TI->cdim;
2183 const PetscInt NqI = TI->Np;
2184 const PetscInt NbI = TI->Nb;
2185 const PetscInt NcI = TI->Nc;
2186 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI];
2187 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2188 const PetscInt NqJ = TJ->Np;
2189 const PetscInt NbJ = TJ->Nb;
2190 const PetscInt NcJ = TJ->Nc;
2191 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2192 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2193 PetscInt f, fc, g, gc, df, dg;
2194 PetscErrorCode ierr;
2195
2196 for (f = 0; f < NbI; ++f) {
2197 for (fc = 0; fc < NcI; ++fc) {
2198 const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2199
2200 tmpBasisI[fidx] = basisI[fidx];
2201 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2202 }
2203 }
2204 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr);
2205 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr);
2206 for (g = 0; g < NbJ; ++g) {
2207 for (gc = 0; gc < NcJ; ++gc) {
2208 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2209
2210 tmpBasisJ[gidx] = basisJ[gidx];
2211 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2212 }
2213 }
2214 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr);
2215 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr);
2216 for (f = 0; f < NbI; ++f) {
2217 for (fc = 0; fc < NcI; ++fc) {
2218 const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2219 const PetscInt i = offsetI+f; /* Element matrix row */
2220 for (g = 0; g < NbJ; ++g) {
2221 for (gc = 0; gc < NcJ; ++gc) {
2222 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2223 const PetscInt j = offsetJ+g; /* Element matrix column */
2224 const PetscInt fOff = eOffset+i*totDim+j;
2225
2226 elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx];
2227 for (df = 0; df < dE; ++df) {
2228 elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2229 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx];
2230 for (dg = 0; dg < dE; ++dg) {
2231 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2232 }
2233 }
2234 }
2235 }
2236 }
2237 }
2238 return(0);
2239 }
2240
PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI,PetscBool isHybridI,PetscFE feJ,PetscBool isHybridJ,PetscInt r,PetscInt q,PetscTabulation TI,PetscScalar tmpBasisI[],PetscScalar tmpBasisDerI[],PetscTabulation TJ,PetscScalar tmpBasisJ[],PetscScalar tmpBasisDerJ[],PetscFEGeom * fegeom,const PetscScalar g0[],const PetscScalar g1[],const PetscScalar g2[],const PetscScalar g3[],PetscInt eOffset,PetscInt totDim,PetscInt offsetI,PetscInt offsetJ,PetscScalar elemMat[])2241 PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[])
2242 {
2243 const PetscInt dE = TI->cdim;
2244 const PetscInt NqI = TI->Np;
2245 const PetscInt NbI = TI->Nb;
2246 const PetscInt NcI = TI->Nc;
2247 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI];
2248 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE];
2249 const PetscInt NqJ = TJ->Np;
2250 const PetscInt NbJ = TJ->Nb;
2251 const PetscInt NcJ = TJ->Nc;
2252 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ];
2253 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE];
2254 const PetscInt Ns = isHybridI ? 1 : 2;
2255 const PetscInt Nt = isHybridJ ? 1 : 2;
2256 PetscInt f, fc, g, gc, df, dg, s, t;
2257 PetscErrorCode ierr;
2258
2259 for (f = 0; f < NbI; ++f) {
2260 for (fc = 0; fc < NcI; ++fc) {
2261 const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2262
2263 tmpBasisI[fidx] = basisI[fidx];
2264 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df];
2265 }
2266 }
2267 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr);
2268 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr);
2269 for (g = 0; g < NbJ; ++g) {
2270 for (gc = 0; gc < NcJ; ++gc) {
2271 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2272
2273 tmpBasisJ[gidx] = basisJ[gidx];
2274 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg];
2275 }
2276 }
2277 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr);
2278 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr);
2279 for (s = 0; s < Ns; ++s) {
2280 for (f = 0; f < NbI; ++f) {
2281 for (fc = 0; fc < NcI; ++fc) {
2282 const PetscInt sc = NcI*s+fc; /* components from each side of the surface */
2283 const PetscInt fidx = f*NcI+fc; /* Test function basis index */
2284 const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */
2285 for (t = 0; t < Nt; ++t) {
2286 for (g = 0; g < NbJ; ++g) {
2287 for (gc = 0; gc < NcJ; ++gc) {
2288 const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */
2289 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
2290 const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */
2291 const PetscInt fOff = eOffset+i*totDim+j;
2292
2293 elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx];
2294 for (df = 0; df < dE; ++df) {
2295 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df];
2296 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx];
2297 for (dg = 0; dg < dE; ++dg) {
2298 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg];
2299 }
2300 }
2301 }
2302 }
2303 }
2304 }
2305 }
2306 }
2307 return(0);
2308 }
2309
PetscFECreateCellGeometry(PetscFE fe,PetscQuadrature quad,PetscFEGeom * cgeom)2310 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom)
2311 {
2312 PetscDualSpace dsp;
2313 DM dm;
2314 PetscQuadrature quadDef;
2315 PetscInt dim, cdim, Nq;
2316 PetscErrorCode ierr;
2317
2318 PetscFunctionBegin;
2319 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr);
2320 ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr);
2321 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
2322 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr);
2323 ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr);
2324 quad = quad ? quad : quadDef;
2325 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr);
2326 ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr);
2327 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr);
2328 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr);
2329 ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr);
2330 cgeom->dim = dim;
2331 cgeom->dimEmbed = cdim;
2332 cgeom->numCells = 1;
2333 cgeom->numPoints = Nq;
2334 ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr);
2335 PetscFunctionReturn(0);
2336 }
2337
PetscFEDestroyCellGeometry(PetscFE fe,PetscFEGeom * cgeom)2338 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom)
2339 {
2340 PetscErrorCode ierr;
2341
2342 PetscFunctionBegin;
2343 ierr = PetscFree(cgeom->v);CHKERRQ(ierr);
2344 ierr = PetscFree(cgeom->J);CHKERRQ(ierr);
2345 ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr);
2346 ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr);
2347 PetscFunctionReturn(0);
2348 }
2349