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