Actual source code: dtfe.c

petsc-master 2016-07-26
Report Typos and Errors
  1: /* Basis Jet Tabulation

  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.

  8:   \psi_i = \sum_k \alpha_{ki} \phi_k

 10: Our nodal basis is defined in terms of the dual basis $n_j$

 12:   n_j \cdot \psi_i = \delta_{ji}

 14: and we may act on the first equation to obtain

 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

 20: so the coefficients of the nodal basis in the prime basis are

 22:    \alpha = V^{-1}

 24: We will define the dual basis vectors $n_j$ using a quadrature rule.

 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$.

 30: I will have to count the degrees correctly for the Legendre product when we are on simplices.

 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 <petsc/private/dtimpl.h>
 39: #include <petsc/private/dmpleximpl.h> /* For CellRefiner */
 40: #include <petscdmshell.h>
 41: #include <petscdmplex.h>
 42: #include <petscblaslapack.h>

 44: PetscBool FEcite = PETSC_FALSE;
 45: const char FECitation[] = "@article{kirby2004,\n"
 46:                           "  title   = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n"
 47:                           "  journal = {ACM Transactions on Mathematical Software},\n"
 48:                           "  author  = {Robert C. Kirby},\n"
 49:                           "  volume  = {30},\n"
 50:                           "  number  = {4},\n"
 51:                           "  pages   = {502--516},\n"
 52:                           "  doi     = {10.1145/1039813.1039820},\n"
 53:                           "  year    = {2004}\n}\n";

 55: PetscClassId PETSCSPACE_CLASSID = 0;

 57: PetscFunctionList PetscSpaceList              = NULL;
 58: PetscBool         PetscSpaceRegisterAllCalled = PETSC_FALSE;

 62: /*@C
 63:   PetscSpaceRegister - Adds a new PetscSpace implementation

 65:   Not Collective

 67:   Input Parameters:
 68: + name        - The name of a new user-defined creation routine
 69: - create_func - The creation routine itself

 71:   Notes:
 72:   PetscSpaceRegister() may be called multiple times to add several user-defined PetscSpaces

 74:   Sample usage:
 75: .vb
 76:     PetscSpaceRegister("my_space", MyPetscSpaceCreate);
 77: .ve

 79:   Then, your PetscSpace type can be chosen with the procedural interface via
 80: .vb
 81:     PetscSpaceCreate(MPI_Comm, PetscSpace *);
 82:     PetscSpaceSetType(PetscSpace, "my_space");
 83: .ve
 84:    or at runtime via the option
 85: .vb
 86:     -petscspace_type my_space
 87: .ve

 89:   Level: advanced

 91: .keywords: PetscSpace, register
 92: .seealso: PetscSpaceRegisterAll(), PetscSpaceRegisterDestroy()

 94: @*/
 95: PetscErrorCode PetscSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscSpace))
 96: {

100:   PetscFunctionListAdd(&PetscSpaceList, sname, function);
101:   return(0);
102: }

106: /*@C
107:   PetscSpaceSetType - Builds a particular PetscSpace

109:   Collective on PetscSpace

111:   Input Parameters:
112: + sp   - The PetscSpace object
113: - name - The kind of space

115:   Options Database Key:
116: . -petscspace_type <type> - Sets the PetscSpace type; use -help for a list of available types

118:   Level: intermediate

120: .keywords: PetscSpace, set, type
121: .seealso: PetscSpaceGetType(), PetscSpaceCreate()
122: @*/
123: PetscErrorCode PetscSpaceSetType(PetscSpace sp, PetscSpaceType name)
124: {
125:   PetscErrorCode (*r)(PetscSpace);
126:   PetscBool      match;

131:   PetscObjectTypeCompare((PetscObject) sp, name, &match);
132:   if (match) return(0);

134:   PetscSpaceRegisterAll();
135:   PetscFunctionListFind(PetscSpaceList, name, &r);
136:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSpace type: %s", name);

138:   if (sp->ops->destroy) {
139:     (*sp->ops->destroy)(sp);
140:     sp->ops->destroy = NULL;
141:   }
142:   (*r)(sp);
143:   PetscObjectChangeTypeName((PetscObject) sp, name);
144:   return(0);
145: }

149: /*@C
150:   PetscSpaceGetType - Gets the PetscSpace type name (as a string) from the object.

152:   Not Collective

154:   Input Parameter:
155: . sp  - The PetscSpace

157:   Output Parameter:
158: . name - The PetscSpace type name

160:   Level: intermediate

162: .keywords: PetscSpace, get, type, name
163: .seealso: PetscSpaceSetType(), PetscSpaceCreate()
164: @*/
165: PetscErrorCode PetscSpaceGetType(PetscSpace sp, PetscSpaceType *name)
166: {

172:   if (!PetscSpaceRegisterAllCalled) {
173:     PetscSpaceRegisterAll();
174:   }
175:   *name = ((PetscObject) sp)->type_name;
176:   return(0);
177: }

181: /*@C
182:   PetscSpaceView - Views a PetscSpace

184:   Collective on PetscSpace

186:   Input Parameter:
187: + sp - the PetscSpace object to view
188: - v  - the viewer

190:   Level: developer

192: .seealso PetscSpaceDestroy()
193: @*/
194: PetscErrorCode PetscSpaceView(PetscSpace sp, PetscViewer v)
195: {

200:   if (!v) {
201:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
202:   }
203:   if (sp->ops->view) {
204:     (*sp->ops->view)(sp, v);
205:   }
206:   return(0);
207: }

211: /*@
212:   PetscSpaceSetFromOptions - sets parameters in a PetscSpace from the options database

214:   Collective on PetscSpace

216:   Input Parameter:
217: . sp - the PetscSpace object to set options for

219:   Options Database:
220: . -petscspace_order the approximation order of the space

222:   Level: developer

224: .seealso PetscSpaceView()
225: @*/
226: PetscErrorCode PetscSpaceSetFromOptions(PetscSpace sp)
227: {
228:   const char    *defaultType;
229:   char           name[256];
230:   PetscBool      flg;

235:   if (!((PetscObject) sp)->type_name) {
236:     defaultType = PETSCSPACEPOLYNOMIAL;
237:   } else {
238:     defaultType = ((PetscObject) sp)->type_name;
239:   }
240:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

242:   PetscObjectOptionsBegin((PetscObject) sp);
243:   PetscOptionsFList("-petscspace_type", "Linear space", "PetscSpaceSetType", PetscSpaceList, defaultType, name, 256, &flg);
244:   if (flg) {
245:     PetscSpaceSetType(sp, name);
246:   } else if (!((PetscObject) sp)->type_name) {
247:     PetscSpaceSetType(sp, defaultType);
248:   }
249:   PetscOptionsInt("-petscspace_order", "The approximation order", "PetscSpaceSetOrder", sp->order, &sp->order, NULL);
250:   if (sp->ops->setfromoptions) {
251:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
252:   }
253:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
254:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
255:   PetscOptionsEnd();
256:   PetscSpaceViewFromOptions(sp, NULL, "-petscspace_view");
257:   return(0);
258: }

262: /*@C
263:   PetscSpaceSetUp - Construct data structures for the PetscSpace

265:   Collective on PetscSpace

267:   Input Parameter:
268: . sp - the PetscSpace object to setup

270:   Level: developer

272: .seealso PetscSpaceView(), PetscSpaceDestroy()
273: @*/
274: PetscErrorCode PetscSpaceSetUp(PetscSpace sp)
275: {

280:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
281:   return(0);
282: }

286: /*@
287:   PetscSpaceDestroy - Destroys a PetscSpace object

289:   Collective on PetscSpace

291:   Input Parameter:
292: . sp - the PetscSpace object to destroy

294:   Level: developer

296: .seealso PetscSpaceView()
297: @*/
298: PetscErrorCode PetscSpaceDestroy(PetscSpace *sp)
299: {

303:   if (!*sp) return(0);

306:   if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
307:   ((PetscObject) (*sp))->refct = 0;
308:   DMDestroy(&(*sp)->dm);

310:   (*(*sp)->ops->destroy)(*sp);
311:   PetscHeaderDestroy(sp);
312:   return(0);
313: }

317: /*@
318:   PetscSpaceCreate - Creates an empty PetscSpace object. The type can then be set with PetscSpaceSetType().

320:   Collective on MPI_Comm

322:   Input Parameter:
323: . comm - The communicator for the PetscSpace object

325:   Output Parameter:
326: . sp - The PetscSpace object

328:   Level: beginner

330: .seealso: PetscSpaceSetType(), PETSCSPACEPOLYNOMIAL
331: @*/
332: PetscErrorCode PetscSpaceCreate(MPI_Comm comm, PetscSpace *sp)
333: {
334:   PetscSpace     s;

339:   PetscCitationsRegister(FECitation,&FEcite);
340:   *sp  = NULL;
341:   PetscFEInitializePackage();

343:   PetscHeaderCreate(s, PETSCSPACE_CLASSID, "PetscSpace", "Linear Space", "PetscSpace", comm, PetscSpaceDestroy, PetscSpaceView);

345:   s->order = 0;
346:   DMShellCreate(comm, &s->dm);

348:   *sp = s;
349:   return(0);
350: }

354: /* Dimension of the space, i.e. number of basis vectors */
355: PetscErrorCode PetscSpaceGetDimension(PetscSpace sp, PetscInt *dim)
356: {

362:   *dim = 0;
363:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
364:   return(0);
365: }

369: /*@
370:   PetscSpaceGetOrder - Return the order of approximation for this space

372:   Input Parameter:
373: . sp - The PetscSpace

375:   Output Parameter:
376: . order - The approximation order

378:   Level: intermediate

380: .seealso: PetscSpaceSetOrder(), PetscSpaceCreate(), PetscSpace
381: @*/
382: PetscErrorCode PetscSpaceGetOrder(PetscSpace sp, PetscInt *order)
383: {
387:   *order = sp->order;
388:   return(0);
389: }

393: /*@
394:   PetscSpaceSetOrder - Set the order of approximation for this space

396:   Input Parameters:
397: + sp - The PetscSpace
398: - order - The approximation order

400:   Level: intermediate

402: .seealso: PetscSpaceGetOrder(), PetscSpaceCreate(), PetscSpace
403: @*/
404: PetscErrorCode PetscSpaceSetOrder(PetscSpace sp, PetscInt order)
405: {
408:   sp->order = order;
409:   return(0);
410: }

414: /*@C
415:   PetscSpaceEvaluate - Evaluate the basis functions and their derivatives (jet) at each point

417:   Input Parameters:
418: + sp      - The PetscSpace
419: . npoints - The number of evaluation points
420: - points  - The point coordinates

422:   Output Parameters:
423: + B - The function evaluations in a npoints x nfuncs array
424: . D - The derivative evaluations in a npoints x nfuncs x dim array
425: - H - The second derivative evaluations in a npoints x nfuncs x dim x dim array

427:   Level: advanced

429: .seealso: PetscFEGetTabulation(), PetscFEGetDefaultTabulation(), PetscSpaceCreate()
430: @*/
431: PetscErrorCode PetscSpaceEvaluate(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
432: {

441:   if (sp->ops->evaluate) {(*sp->ops->evaluate)(sp, npoints, points, B, D, H);}
442:   return(0);
443: }

447: PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
448: {
449:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
450:   PetscErrorCode   ierr;

453:   PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
454:   PetscOptionsInt("-petscspace_poly_num_variables", "The number of different variables, e.g. x and y", "PetscSpacePolynomialSetNumVariables", poly->numVariables, &poly->numVariables, NULL);
455:   PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
456:   PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
457:   PetscOptionsTail();
458:   return(0);
459: }

463: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer viewer)
464: {
465:   PetscSpace_Poly  *poly = (PetscSpace_Poly *) sp->data;
466:   PetscViewerFormat format;
467:   PetscErrorCode    ierr;

470:   PetscViewerGetFormat(viewer, &format);
471:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
472:     if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
473:     else              {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
474:   } else {
475:     if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
476:     else              {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
477:   }
478:   return(0);
479: }

483: PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
484: {
485:   PetscBool      iascii;

491:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
492:   if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
493:   return(0);
494: }

498: PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
499: {
500:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
501:   PetscInt         ndegree = sp->order+1;
502:   PetscInt         deg;
503:   PetscErrorCode   ierr;

506:   PetscMalloc1(ndegree, &poly->degrees);
507:   for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
508:   return(0);
509: }

513: PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
514: {
515:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
516:   PetscErrorCode   ierr;

519:   PetscFree(poly->degrees);
520:   PetscFree(poly);
521:   return(0);
522: }

526: PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
527: {
528:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
529:   PetscInt         deg  = sp->order;
530:   PetscInt         n    = poly->numVariables, i;
531:   PetscReal        D    = 1.0;

534:   if (poly->tensor) {
535:     *dim = 1;
536:     for (i = 0; i < n; ++i) *dim *= (deg+1);
537:   } else {
538:     for (i = 1; i <= n; ++i) {
539:       D *= ((PetscReal) (deg+i))/i;
540:     }
541:     *dim = (PetscInt) (D + 0.5);
542:   }
543:   return(0);
544: }

548: /*
549:   LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'.

551:   Input Parameters:
552: + len - The length of the tuple
553: . sum - The sum of all entries in the tuple
554: - ind - The current multi-index of the tuple, initialized to the 0 tuple

556:   Output Parameter:
557: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
558: . tup - A tuple of len integers addig to sum

560:   Level: developer

562: .seealso: 
563: */
564: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
565: {
566:   PetscInt       i;

570:   if (len == 1) {
571:     ind[0] = -1;
572:     tup[0] = sum;
573:   } else if (sum == 0) {
574:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
575:   } else {
576:     tup[0] = sum - ind[0];
577:     LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
578:     if (ind[1] < 0) {
579:       if (ind[0] == sum) {ind[0] = -1;}
580:       else               {ind[1] = 0; ++ind[0];}
581:     }
582:   }
583:   return(0);
584: }

588: /*
589:   TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'.

591:   Input Parameters:
592: + len - The length of the tuple
593: . max - The max for all entries in the tuple
594: - ind - The current multi-index of the tuple, initialized to the 0 tuple

596:   Output Parameter:
597: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
598: . tup - A tuple of len integers less than max

600:   Level: developer

602: .seealso: 
603: */
604: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
605: {
606:   PetscInt       i;

610:   if (len == 1) {
611:     tup[0] = ind[0]++;
612:     ind[0] = ind[0] >= max ? -1 : ind[0];
613:   } else if (max == 0) {
614:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
615:   } else {
616:     tup[0] = ind[0];
617:     TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
618:     if (ind[1] < 0) {
619:       ind[1] = 0;
620:       if (ind[0] == max-1) {ind[0] = -1;}
621:       else                 {++ind[0];}
622:     }
623:   }
624:   return(0);
625: }

629: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
630: {
631:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
632:   DM               dm      = sp->dm;
633:   PetscInt         ndegree = sp->order+1;
634:   PetscInt        *degrees = poly->degrees;
635:   PetscInt         dim     = poly->numVariables;
636:   PetscReal       *lpoints, *tmp, *LB, *LD, *LH;
637:   PetscInt        *ind, *tup;
638:   PetscInt         pdim, d, der, i, p, deg, o;
639:   PetscErrorCode   ierr;

642:   PetscSpaceGetDimension(sp, &pdim);
643:   DMGetWorkArray(dm, npoints, PETSC_REAL, &lpoints);
644:   DMGetWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
645:   if (B) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
646:   if (D) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
647:   if (H) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
648:   for (d = 0; d < dim; ++d) {
649:     for (p = 0; p < npoints; ++p) {
650:       lpoints[p] = points[p*dim+d];
651:     }
652:     PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
653:     /* LB, LD, LH (ndegree * dim x npoints) */
654:     for (deg = 0; deg < ndegree; ++deg) {
655:       for (p = 0; p < npoints; ++p) {
656:         if (B) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
657:         if (D) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
658:         if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
659:       }
660:     }
661:   }
662:   /* Multiply by A (pdim x ndegree * dim) */
663:   PetscMalloc2(dim,&ind,dim,&tup);
664:   if (B) {
665:     /* B (npoints x pdim) */
666:     if (poly->tensor) {
667:       i = 0;
668:       PetscMemzero(ind, dim * sizeof(PetscInt));
669:       while (ind[0] >= 0) {
670:         TensorPoint_Internal(dim, sp->order+1, ind, tup);
671:         for (p = 0; p < npoints; ++p) {
672:           B[p*pdim + i] = 1.0;
673:           for (d = 0; d < dim; ++d) {
674:             B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
675:           }
676:         }
677:         ++i;
678:       }
679:     } else {
680:       i = 0;
681:       for (o = 0; o <= sp->order; ++o) {
682:         PetscMemzero(ind, dim * sizeof(PetscInt));
683:         while (ind[0] >= 0) {
684:           LatticePoint_Internal(dim, o, ind, tup);
685:           for (p = 0; p < npoints; ++p) {
686:             B[p*pdim + i] = 1.0;
687:             for (d = 0; d < dim; ++d) {
688:               B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
689:             }
690:           }
691:           ++i;
692:         }
693:       }
694:     }
695:   }
696:   if (D) {
697:     /* D (npoints x pdim x dim) */
698:     if (poly->tensor) {
699:       i = 0;
700:       PetscMemzero(ind, dim * sizeof(PetscInt));
701:       while (ind[0] >= 0) {
702:         TensorPoint_Internal(dim, sp->order+1, ind, tup);
703:         for (p = 0; p < npoints; ++p) {
704:           for (der = 0; der < dim; ++der) {
705:             D[(p*pdim + i)*dim + der] = 1.0;
706:             for (d = 0; d < dim; ++d) {
707:               if (d == der) {
708:                 D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
709:               } else {
710:                 D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
711:               }
712:             }
713:           }
714:         }
715:         ++i;
716:       }
717:     } else {
718:       i = 0;
719:       for (o = 0; o <= sp->order; ++o) {
720:         PetscMemzero(ind, dim * sizeof(PetscInt));
721:         while (ind[0] >= 0) {
722:           LatticePoint_Internal(dim, o, ind, tup);
723:           for (p = 0; p < npoints; ++p) {
724:             for (der = 0; der < dim; ++der) {
725:               D[(p*pdim + i)*dim + der] = 1.0;
726:               for (d = 0; d < dim; ++d) {
727:                 if (d == der) {
728:                   D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
729:                 } else {
730:                   D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
731:                 }
732:               }
733:             }
734:           }
735:           ++i;
736:         }
737:       }
738:     }
739:   }
740:   if (H) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to code second derivatives");
741:   PetscFree2(ind,tup);
742:   if (B) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
743:   if (D) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
744:   if (H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
745:   DMRestoreWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
746:   DMRestoreWorkArray(dm, npoints, PETSC_REAL, &lpoints);
747:   return(0);
748: }

752: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
753: {
755:   sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
756:   sp->ops->setup          = PetscSpaceSetUp_Polynomial;
757:   sp->ops->view           = PetscSpaceView_Polynomial;
758:   sp->ops->destroy        = PetscSpaceDestroy_Polynomial;
759:   sp->ops->getdimension   = PetscSpaceGetDimension_Polynomial;
760:   sp->ops->evaluate       = PetscSpaceEvaluate_Polynomial;
761:   return(0);
762: }

764: /*MC
765:   PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of linear polynomials.

767:   Level: intermediate

769: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
770: M*/

774: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
775: {
776:   PetscSpace_Poly *poly;
777:   PetscErrorCode   ierr;

781:   PetscNewLog(sp,&poly);
782:   sp->data = poly;

784:   poly->numVariables = 0;
785:   poly->symmetric    = PETSC_FALSE;
786:   poly->tensor       = PETSC_FALSE;
787:   poly->degrees      = NULL;

789:   PetscSpaceInitialize_Polynomial(sp);
790:   return(0);
791: }

795: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
796: {
797:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

801:   poly->symmetric = sym;
802:   return(0);
803: }

807: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
808: {
809:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

814:   *sym = poly->symmetric;
815:   return(0);
816: }

820: /*@
821:   PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
822:   by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
823:   spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).

825:   Input Parameters:
826: + sp     - the function space object
827: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

829:   Level: beginner

831: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetOrder(), PetscSpacePolynomialSetNumVariables()
832: @*/
833: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
834: {
835:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

839:   poly->tensor = tensor;
840:   return(0);
841: }

845: /*@
846:   PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
847:   by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
848:   spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).

850:   Input Parameters:
851: . sp     - the function space object

853:   Output Parameters:
854: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

856:   Level: beginner

858: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetOrder(), PetscSpacePolynomialSetNumVariables()
859: @*/
860: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
861: {
862:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

867:   *tensor = poly->tensor;
868:   return(0);
869: }

873: PetscErrorCode PetscSpacePolynomialSetNumVariables(PetscSpace sp, PetscInt n)
874: {
875:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

879:   poly->numVariables = n;
880:   return(0);
881: }

885: PetscErrorCode PetscSpacePolynomialGetNumVariables(PetscSpace sp, PetscInt *n)
886: {
887:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

892:   *n = poly->numVariables;
893:   return(0);
894: }

898: PetscErrorCode PetscSpaceSetFromOptions_DG(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
899: {
900:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

904:   PetscOptionsHead(PetscOptionsObject,"PetscSpace DG options");
905:   PetscOptionsInt("-petscspace_dg_num_variables", "The number of different variables, e.g. x and y", "PetscSpaceDGSetNumVariables", dg->numVariables, &dg->numVariables, NULL);
906:   PetscOptionsTail();
907:   return(0);
908: }

912: PetscErrorCode PetscSpaceDGView_Ascii(PetscSpace sp, PetscViewer viewer)
913: {
914:   PetscSpace_DG    *dg = (PetscSpace_DG *) sp->data;
915:   PetscViewerFormat format;
916:   PetscErrorCode    ierr;

919:   PetscViewerGetFormat(viewer, &format);
920:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
921:     PetscViewerASCIIPrintf(viewer, "DG space in dimension %d:\n", dg->numVariables);
922:     PetscViewerASCIIPushTab(viewer);
923:     PetscQuadratureView(dg->quad, viewer);
924:     PetscViewerASCIIPopTab(viewer);
925:   } else {
926:     PetscViewerASCIIPrintf(viewer, "DG space in dimension %d on %d points\n", dg->numVariables, dg->quad->numPoints);
927:   }
928:   return(0);
929: }

933: PetscErrorCode PetscSpaceView_DG(PetscSpace sp, PetscViewer viewer)
934: {
935:   PetscBool      iascii;

941:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
942:   if (iascii) {PetscSpaceDGView_Ascii(sp, viewer);}
943:   return(0);
944: }

948: PetscErrorCode PetscSpaceSetUp_DG(PetscSpace sp)
949: {
950:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

954:   if (!dg->quad->points && sp->order) {
955:     PetscDTGaussJacobiQuadrature(dg->numVariables, sp->order, -1.0, 1.0, &dg->quad);
956:   }
957:   return(0);
958: }

962: PetscErrorCode PetscSpaceDestroy_DG(PetscSpace sp)
963: {
964:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

968:   PetscQuadratureDestroy(&dg->quad);
969:   return(0);
970: }

974: PetscErrorCode PetscSpaceGetDimension_DG(PetscSpace sp, PetscInt *dim)
975: {
976:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

979:   *dim = dg->quad->numPoints;
980:   return(0);
981: }

985: PetscErrorCode PetscSpaceEvaluate_DG(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
986: {
987:   PetscSpace_DG *dg  = (PetscSpace_DG *) sp->data;
988:   PetscInt       dim = dg->numVariables, d, p;

992:   if (D || H) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Cannot calculate derivatives for a DG space");
993:   if (npoints != dg->quad->numPoints) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot evaluate DG space on %d points != %d size", npoints, dg->quad->numPoints);
994:   PetscMemzero(B, npoints*npoints * sizeof(PetscReal));
995:   for (p = 0; p < npoints; ++p) {
996:     for (d = 0; d < dim; ++d) {
997:       if (PetscAbsReal(points[p*dim+d] - dg->quad->points[p*dim+d]) > 1.0e-10) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot evaluate DG point (%d, %d) %g != %g", p, d, points[p*dim+d], dg->quad->points[p*dim+d]);
998:     }
999:     B[p*npoints+p] = 1.0;
1000:   }
1001:   return(0);
1002: }

1006: PetscErrorCode PetscSpaceInitialize_DG(PetscSpace sp)
1007: {
1009:   sp->ops->setfromoptions = PetscSpaceSetFromOptions_DG;
1010:   sp->ops->setup          = PetscSpaceSetUp_DG;
1011:   sp->ops->view           = PetscSpaceView_DG;
1012:   sp->ops->destroy        = PetscSpaceDestroy_DG;
1013:   sp->ops->getdimension   = PetscSpaceGetDimension_DG;
1014:   sp->ops->evaluate       = PetscSpaceEvaluate_DG;
1015:   return(0);
1016: }

1018: /*MC
1019:   PETSCSPACEDG = "dg" - A PetscSpace object that encapsulates functions defined on a set of quadrature points.

1021:   Level: intermediate

1023: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
1024: M*/

1028: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_DG(PetscSpace sp)
1029: {
1030:   PetscSpace_DG *dg;

1035:   PetscNewLog(sp,&dg);
1036:   sp->data = dg;

1038:   dg->numVariables    = 0;
1039:   dg->quad->dim       = 0;
1040:   dg->quad->numPoints = 0;
1041:   dg->quad->points    = NULL;
1042:   dg->quad->weights   = NULL;

1044:   PetscSpaceInitialize_DG(sp);
1045:   return(0);
1046: }


1049: PetscClassId PETSCDUALSPACE_CLASSID = 0;

1051: PetscFunctionList PetscDualSpaceList              = NULL;
1052: PetscBool         PetscDualSpaceRegisterAllCalled = PETSC_FALSE;

1056: /*@C
1057:   PetscDualSpaceRegister - Adds a new PetscDualSpace implementation

1059:   Not Collective

1061:   Input Parameters:
1062: + name        - The name of a new user-defined creation routine
1063: - create_func - The creation routine itself

1065:   Notes:
1066:   PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces

1068:   Sample usage:
1069: .vb
1070:     PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
1071: .ve

1073:   Then, your PetscDualSpace type can be chosen with the procedural interface via
1074: .vb
1075:     PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
1076:     PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
1077: .ve
1078:    or at runtime via the option
1079: .vb
1080:     -petscdualspace_type my_dual_space
1081: .ve

1083:   Level: advanced

1085: .keywords: PetscDualSpace, register
1086: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()

1088: @*/
1089: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
1090: {

1094:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
1095:   return(0);
1096: }

1100: /*@C
1101:   PetscDualSpaceSetType - Builds a particular PetscDualSpace

1103:   Collective on PetscDualSpace

1105:   Input Parameters:
1106: + sp   - The PetscDualSpace object
1107: - name - The kind of space

1109:   Options Database Key:
1110: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types

1112:   Level: intermediate

1114: .keywords: PetscDualSpace, set, type
1115: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
1116: @*/
1117: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
1118: {
1119:   PetscErrorCode (*r)(PetscDualSpace);
1120:   PetscBool      match;

1125:   PetscObjectTypeCompare((PetscObject) sp, name, &match);
1126:   if (match) return(0);

1128:   if (!PetscDualSpaceRegisterAllCalled) {PetscDualSpaceRegisterAll();}
1129:   PetscFunctionListFind(PetscDualSpaceList, name, &r);
1130:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);

1132:   if (sp->ops->destroy) {
1133:     (*sp->ops->destroy)(sp);
1134:     sp->ops->destroy = NULL;
1135:   }
1136:   (*r)(sp);
1137:   PetscObjectChangeTypeName((PetscObject) sp, name);
1138:   return(0);
1139: }

1143: /*@C
1144:   PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.

1146:   Not Collective

1148:   Input Parameter:
1149: . sp  - The PetscDualSpace

1151:   Output Parameter:
1152: . name - The PetscDualSpace type name

1154:   Level: intermediate

1156: .keywords: PetscDualSpace, get, type, name
1157: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
1158: @*/
1159: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
1160: {

1166:   if (!PetscDualSpaceRegisterAllCalled) {
1167:     PetscDualSpaceRegisterAll();
1168:   }
1169:   *name = ((PetscObject) sp)->type_name;
1170:   return(0);
1171: }

1175: /*@
1176:   PetscDualSpaceView - Views a PetscDualSpace

1178:   Collective on PetscDualSpace

1180:   Input Parameter:
1181: + sp - the PetscDualSpace object to view
1182: - v  - the viewer

1184:   Level: developer

1186: .seealso PetscDualSpaceDestroy()
1187: @*/
1188: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
1189: {

1194:   if (!v) {
1195:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
1196:   }
1197:   if (sp->ops->view) {
1198:     (*sp->ops->view)(sp, v);
1199:   }
1200:   return(0);
1201: }

1205: /*@
1206:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

1208:   Collective on PetscDualSpace

1210:   Input Parameter:
1211: . sp - the PetscDualSpace object to set options for

1213:   Options Database:
1214: . -petscspace_order the approximation order of the space

1216:   Level: developer

1218: .seealso PetscDualSpaceView()
1219: @*/
1220: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
1221: {
1222:   const char    *defaultType;
1223:   char           name[256];
1224:   PetscBool      flg;

1229:   if (!((PetscObject) sp)->type_name) {
1230:     defaultType = PETSCDUALSPACELAGRANGE;
1231:   } else {
1232:     defaultType = ((PetscObject) sp)->type_name;
1233:   }
1234:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

1236:   PetscObjectOptionsBegin((PetscObject) sp);
1237:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
1238:   if (flg) {
1239:     PetscDualSpaceSetType(sp, name);
1240:   } else if (!((PetscObject) sp)->type_name) {
1241:     PetscDualSpaceSetType(sp, defaultType);
1242:   }
1243:   PetscOptionsInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
1244:   if (sp->ops->setfromoptions) {
1245:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
1246:   }
1247:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1248:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
1249:   PetscOptionsEnd();
1250:   PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
1251:   return(0);
1252: }

1256: /*@
1257:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

1259:   Collective on PetscDualSpace

1261:   Input Parameter:
1262: . sp - the PetscDualSpace object to setup

1264:   Level: developer

1266: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
1267: @*/
1268: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
1269: {

1274:   if (sp->setupcalled) return(0);
1275:   sp->setupcalled = PETSC_TRUE;
1276:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
1277:   return(0);
1278: }

1282: /*@
1283:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

1285:   Collective on PetscDualSpace

1287:   Input Parameter:
1288: . sp - the PetscDualSpace object to destroy

1290:   Level: developer

1292: .seealso PetscDualSpaceView()
1293: @*/
1294: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
1295: {
1296:   PetscInt       dim, f;

1300:   if (!*sp) return(0);

1303:   if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
1304:   ((PetscObject) (*sp))->refct = 0;

1306:   PetscDualSpaceGetDimension(*sp, &dim);
1307:   for (f = 0; f < dim; ++f) {
1308:     PetscQuadratureDestroy(&(*sp)->functional[f]);
1309:   }
1310:   PetscFree((*sp)->functional);
1311:   DMDestroy(&(*sp)->dm);

1313:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
1314:   PetscHeaderDestroy(sp);
1315:   return(0);
1316: }

1320: /*@
1321:   PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().

1323:   Collective on MPI_Comm

1325:   Input Parameter:
1326: . comm - The communicator for the PetscDualSpace object

1328:   Output Parameter:
1329: . sp - The PetscDualSpace object

1331:   Level: beginner

1333: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
1334: @*/
1335: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
1336: {
1337:   PetscDualSpace s;

1342:   PetscCitationsRegister(FECitation,&FEcite);
1343:   *sp  = NULL;
1344:   PetscFEInitializePackage();

1346:   PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);

1348:   s->order = 0;
1349:   s->setupcalled = PETSC_FALSE;

1351:   *sp = s;
1352:   return(0);
1353: }

1357: /*@
1358:   PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.

1360:   Collective on PetscDualSpace

1362:   Input Parameter:
1363: . sp - The original PetscDualSpace

1365:   Output Parameter:
1366: . spNew - The duplicate PetscDualSpace

1368:   Level: beginner

1370: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
1371: @*/
1372: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
1373: {

1379:   (*sp->ops->duplicate)(sp, spNew);
1380:   return(0);
1381: }

1385: /*@
1386:   PetscDualSpaceGetDM - Get the DM representing the reference cell

1388:   Not collective

1390:   Input Parameter:
1391: . sp - The PetscDualSpace

1393:   Output Parameter:
1394: . dm - The reference cell

1396:   Level: intermediate

1398: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
1399: @*/
1400: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
1401: {
1405:   *dm = sp->dm;
1406:   return(0);
1407: }

1411: /*@
1412:   PetscDualSpaceSetDM - Get the DM representing the reference cell

1414:   Not collective

1416:   Input Parameters:
1417: + sp - The PetscDualSpace
1418: - dm - The reference cell

1420:   Level: intermediate

1422: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
1423: @*/
1424: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
1425: {

1431:   DMDestroy(&sp->dm);
1432:   PetscObjectReference((PetscObject) dm);
1433:   sp->dm = dm;
1434:   return(0);
1435: }

1439: /*@
1440:   PetscDualSpaceGetOrder - Get the order of the dual space

1442:   Not collective

1444:   Input Parameter:
1445: . sp - The PetscDualSpace

1447:   Output Parameter:
1448: . order - The order

1450:   Level: intermediate

1452: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
1453: @*/
1454: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
1455: {
1459:   *order = sp->order;
1460:   return(0);
1461: }

1465: /*@
1466:   PetscDualSpaceSetOrder - Set the order of the dual space

1468:   Not collective

1470:   Input Parameters:
1471: + sp - The PetscDualSpace
1472: - order - The order

1474:   Level: intermediate

1476: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
1477: @*/
1478: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
1479: {
1482:   sp->order = order;
1483:   return(0);
1484: }

1488: /*@
1489:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

1491:   Not collective

1493:   Input Parameters:
1494: + sp - The PetscDualSpace
1495: - i  - The basis number

1497:   Output Parameter:
1498: . functional - The basis functional

1500:   Level: intermediate

1502: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
1503: @*/
1504: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
1505: {
1506:   PetscInt       dim;

1512:   PetscDualSpaceGetDimension(sp, &dim);
1513:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
1514:   *functional = sp->functional[i];
1515:   return(0);
1516: }

1520: /*@
1521:   PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals

1523:   Not collective

1525:   Input Parameter:
1526: . sp - The PetscDualSpace

1528:   Output Parameter:
1529: . dim - The dimension

1531:   Level: intermediate

1533: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1534: @*/
1535: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
1536: {

1542:   *dim = 0;
1543:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
1544:   return(0);
1545: }

1549: /*@C
1550:   PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension

1552:   Not collective

1554:   Input Parameter:
1555: . sp - The PetscDualSpace

1557:   Output Parameter:
1558: . numDof - An array of length dim+1 which holds the number of dofs for each dimension

1560:   Level: intermediate

1562: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1563: @*/
1564: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
1565: {

1571:   (*sp->ops->getnumdof)(sp, numDof);
1572:   if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
1573:   return(0);
1574: }

1578: /*@
1579:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

1581:   Collective on PetscDualSpace

1583:   Input Parameters:
1584: + sp      - The PetscDualSpace
1585: . dim     - The spatial dimension
1586: - simplex - Flag for simplex, otherwise use a tensor-product cell

1588:   Output Parameter:
1589: . refdm - The reference cell

1591:   Level: advanced

1593: .keywords: PetscDualSpace, reference cell
1594: .seealso: PetscDualSpaceCreate(), DMPLEX
1595: @*/
1596: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1597: {

1601:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
1602:   return(0);
1603: }

1607: /*@C
1608:   PetscDualSpaceApply - Apply a functional from the dual space basis to an input function

1610:   Input Parameters:
1611: + sp      - The PetscDualSpace object
1612: . f       - The basis functional index
1613: . time    - The time
1614: . geom    - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1615: . numComp - The number of components for the function
1616: . func    - The input function
1617: - ctx     - A context for the function

1619:   Output Parameter:
1620: . value   - numComp output values

1622:   Note: The calling sequence for the callback func is given by:

1624: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1625: $      PetscInt numComponents, PetscScalar values[], void *ctx)

1627:   Level: developer

1629: .seealso: PetscDualSpaceCreate()
1630: @*/
1631: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFECellGeom *geom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1632: {
1633:   DM               dm;
1634:   PetscQuadrature  quad;
1635:   PetscReal        x[3];
1636:   PetscScalar     *val;
1637:   PetscInt         dim, q, c;
1638:   PetscErrorCode   ierr;

1643:   dim  = geom->dim;
1644:   PetscDualSpaceGetDM(sp, &dm);
1645:   PetscDualSpaceGetFunctional(sp, f, &quad);
1646:   DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1647:   for (c = 0; c < numComp; ++c) value[c] = 0.0;
1648:   for (q = 0; q < quad->numPoints; ++q) {
1649:     CoordinatesRefToReal(geom->dimEmbed, dim, geom->v0, geom->J, &quad->points[q*dim], x);
1650:     (*func)(geom->dimEmbed, time, x, numComp, val, ctx);
1651:     for (c = 0; c < numComp; ++c) {
1652:       value[c] += val[c]*quad->weights[q];
1653:     }
1654:   }
1655:   DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1656:   return(0);
1657: }

1661: /*@C
1662:   PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function

1664:   Input Parameters:
1665: + sp      - The PetscDualSpace object
1666: . f       - The basis functional index
1667: . time    - The time
1668: . geom    - A context with geometric information for this cell, we currently just use the centroid
1669: . numComp - The number of components for the function
1670: . func    - The input function
1671: - ctx     - A context for the function

1673:   Output Parameter:
1674: . value   - numComp output values

1676:   Note: The calling sequence for the callback func is given by:

1678: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1679: $      PetscInt numComponents, PetscScalar values[], void *ctx)

1681:   Level: developer

1683: .seealso: PetscDualSpaceCreate()
1684: @*/
1685: PetscErrorCode PetscDualSpaceApplyFVM(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFVCellGeom *geom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1686: {
1687:   DM               dm;
1688:   PetscQuadrature  quad;
1689:   PetscScalar     *val;
1690:   PetscInt         dimEmbed, q, c;
1691:   PetscErrorCode   ierr;

1696:   PetscDualSpaceGetDM(sp, &dm);
1697:   DMGetCoordinateDim(dm, &dimEmbed);
1698:   PetscDualSpaceGetFunctional(sp, f, &quad);
1699:   DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1700:   for (c = 0; c < numComp; ++c) value[c] = 0.0;
1701:   for (q = 0; q < quad->numPoints; ++q) {
1702:     (*func)(dimEmbed, time, geom->centroid, numComp, val, ctx);
1703:     for (c = 0; c < numComp; ++c) {
1704:       value[c] += val[c]*quad->weights[q];
1705:     }
1706:   }
1707:   DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1708:   return(0);
1709: }

1713: /*@
1714:   PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a given height.

1716:   If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1717:   pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
1718:   support extracting subspaces, then NULL is returned.

1720:   This does not increment the reference count on the returned dual space, and the user should not destroy it.

1722:   Not collective

1724:   Input Parameters:
1725: + sp - the PetscDualSpace object
1726: - height - the height of the mesh point for which the subspace is desired

1728:   Output Parameters:
1729:   bdsp - the subspace: must be destroyed by the user

1731:   Level: advanced

1733: .seealso: PetscDualSpace
1734: @*/
1735: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
1736: {

1742:   *bdsp = NULL;
1743:   if (sp->ops->getheightsubspace) {
1744:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1745:   }
1746:   return(0);
1747: }

1751: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
1752: {
1753:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1754:   PetscReal           D   = 1.0;
1755:   PetscInt            n, i;
1756:   PetscErrorCode      ierr;

1759:   *dim = -1;                    /* Ensure that the compiler knows *dim is set. */
1760:   DMGetDimension(sp->dm, &n);
1761:   if (lag->simplex || !lag->continuous) {
1762:     for (i = 1; i <= n; ++i) {
1763:       D *= ((PetscReal) (order+i))/i;
1764:     }
1765:     *dim = (PetscInt) (D + 0.5);
1766:   } else {
1767:     *dim = 1;
1768:     for (i = 0; i < n; ++i) *dim *= (order+1);
1769:   }
1770:   return(0);
1771: }

1775: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
1776: {
1777:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1778:   PetscBool          continuous;
1779:   PetscInt           order;
1780:   PetscErrorCode     ierr;

1785:   PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
1786:   PetscDualSpaceGetOrder(sp,&order);
1787:   if (height == 0) {
1788:     PetscObjectReference((PetscObject)sp);
1789:     *bdsp = sp;
1790:   }
1791:   else if (continuous == PETSC_FALSE || !order) {
1792:     *bdsp = NULL;
1793:   }
1794:   else {
1795:     DM dm, K;
1796:     PetscInt dim;

1798:     PetscDualSpaceGetDM(sp,&dm);
1799:     DMGetDimension(dm,&dim);
1800:     if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim);}
1801:     PetscDualSpaceDuplicate(sp,bdsp);
1802:     PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplex, &K);
1803:     PetscDualSpaceSetDM(*bdsp, K);
1804:     DMDestroy(&K);
1805:     PetscDualSpaceSetUp(*bdsp);
1806:   }
1807:   return(0);
1808: }

1812: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1813: {
1814:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1815:   DM                  dm    = sp->dm;
1816:   PetscInt            order = sp->order;
1817:   PetscBool           disc  = lag->continuous ? PETSC_FALSE : PETSC_TRUE;
1818:   PetscSection        csection;
1819:   Vec                 coordinates;
1820:   PetscReal          *qpoints, *qweights;
1821:   PetscInt           *closure = NULL, closureSize, c;
1822:   PetscInt            depth, dim, pdimMax, pMax = 0, *pStart, *pEnd, cell, coneSize, d, n, f = 0;
1823:   PetscBool           simplex;
1824:   PetscErrorCode      ierr;

1827:   /* Classify element type */
1828:   DMGetDimension(dm, &dim);
1829:   DMPlexGetDepth(dm, &depth);
1830:   PetscCalloc1(dim+1, &lag->numDof);
1831:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
1832:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1833:   DMPlexGetConeSize(dm, pStart[depth], &coneSize);
1834:   DMGetCoordinateSection(dm, &csection);
1835:   DMGetCoordinatesLocal(dm, &coordinates);
1836:   if (depth == 1) {
1837:     if      (coneSize == dim+1)    simplex = PETSC_TRUE;
1838:     else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
1839:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1840:   }
1841:   else if (depth == dim) {
1842:     if      (coneSize == dim+1)   simplex = PETSC_TRUE;
1843:     else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
1844:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1845:   }
1846:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
1847:   lag->simplex = simplex;
1848:   PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
1849:   pdimMax *= (pEnd[dim] - pStart[dim]);
1850:   PetscMalloc1(pdimMax, &sp->functional);
1851:   for (d = 0; d <= depth; d++) {
1852:     pMax = PetscMax(pMax,pEnd[d]);
1853:   }
1854:   if (!dim) {
1855:     PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1856:     PetscMalloc1(1, &qpoints);
1857:     PetscMalloc1(1, &qweights);
1858:     PetscQuadratureSetOrder(sp->functional[f], 0);
1859:     PetscQuadratureSetData(sp->functional[f], PETSC_DETERMINE, 1, qpoints, qweights);
1860:     qpoints[0]  = 0.0;
1861:     qweights[0] = 1.0;
1862:     ++f;
1863:     lag->numDof[0] = 1;
1864:   } else {
1865:     PetscBT seen;

1867:     PetscBTCreate(pMax, &seen);
1868:     PetscBTMemzero(pMax, seen);
1869:     for (cell = pStart[depth]; cell < pEnd[depth]; ++cell) {
1870:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1871:       for (c = 0; c < closureSize*2; c += 2) {
1872:         const PetscInt p = closure[c];

1874:         if (PetscBTLookup(seen, p)) continue;
1875:         PetscBTSet(seen, p);
1876:         if ((p >= pStart[0]) && (p < pEnd[0])) {
1877:           /* Vertices */
1878:           const PetscScalar *coords;
1879:           PetscInt           dof, off, d;

1881:           if (order < 1 || disc) continue;
1882:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1883:           PetscMalloc1(dim, &qpoints);
1884:           PetscMalloc1(1, &qweights);
1885:           PetscQuadratureSetOrder(sp->functional[f], 0);
1886:           PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1887:           VecGetArrayRead(coordinates, &coords);
1888:           PetscSectionGetDof(csection, p, &dof);
1889:           PetscSectionGetOffset(csection, p, &off);
1890:           if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim);
1891:           for (d = 0; d < dof; ++d) {qpoints[d] = PetscRealPart(coords[off+d]);}
1892:           qweights[0] = 1.0;
1893:           ++f;
1894:           VecRestoreArrayRead(coordinates, &coords);
1895:           lag->numDof[0] = 1;
1896:         } else if ((p >= pStart[1]) && (p < pEnd[1])) {
1897:           /* Edges */
1898:           PetscScalar *coords;
1899:           PetscInt     num = ((dim == 1) && !order) ? 1 : order-1, k;

1901:           if (num < 1 || disc) continue;
1902:           coords = NULL;
1903:           DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);
1904:           if (n != dim*2) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d has %d coordinate values instead of %d", p, n, dim*2);
1905:           for (k = 1; k <= num; ++k) {
1906:             PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1907:             PetscMalloc1(dim, &qpoints);
1908:             PetscMalloc1(1, &qweights);
1909:             PetscQuadratureSetOrder(sp->functional[f], 0);
1910:             PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1911:             for (d = 0; d < dim; ++d) {qpoints[d] = k*PetscRealPart(coords[1*dim+d] - coords[0*dim+d])/order + PetscRealPart(coords[0*dim+d]);}
1912:             qweights[0] = 1.0;
1913:             ++f;
1914:           }
1915:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);
1916:           lag->numDof[1] = num;
1917:         } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) {
1918:           /* Faces */

1920:           if (disc) continue;
1921:           if ( simplex && (order < 3)) continue;
1922:           if (!simplex && (order < 2)) continue;
1923:           lag->numDof[depth-1] = 0;
1924:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces");
1925:         } else if ((p >= pStart[depth]) && (p < pEnd[depth])) {
1926:           /* Cells */
1927:           PetscInt     orderEff = lag->continuous && order ? (simplex ? order-3 : order-2) : order;
1928:           PetscReal    denom    = order ? (lag->continuous ? order : (simplex ? order+3 : order+2)) : (simplex ? 3 : 2);
1929:           PetscScalar *coords   = NULL;
1930:           PetscReal    dx = 2.0/denom, *v0, *J, *invJ, detJ;
1931:           PetscInt    *ind, *tup;
1932:           PetscInt     cdim, csize, d, d2, o;

1934:           lag->numDof[depth] = 0;
1935:           if (orderEff < 0) continue;
1936:           PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
1937:           DMPlexVecGetClosure(dm, csection, coordinates, p, &csize, &coords);
1938:           if (csize%dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate size %d is not divisible by spatial dimension %d", csize, dim);

1940:           PetscCalloc5(dim,&ind,dim,&tup,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1941:           DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
1942:           if (simplex || disc) {
1943:             for (o = 0; o <= orderEff; ++o) {
1944:               PetscMemzero(ind, dim*sizeof(PetscInt));
1945:               while (ind[0] >= 0) {
1946:                 PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1947:                 PetscMalloc1(dim, &qpoints);
1948:                 PetscMalloc1(1,   &qweights);
1949:                 PetscQuadratureSetOrder(sp->functional[f], 0);
1950:                 PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1951:                 LatticePoint_Internal(dim, o, ind, tup);
1952:                 for (d = 0; d < dim; ++d) {
1953:                   qpoints[d] = v0[d];
1954:                   for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
1955:                 }
1956:                 qweights[0] = 1.0;
1957:                 ++f;
1958:               }
1959:             }
1960:           } else {
1961:             while (ind[0] >= 0) {
1962:               PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1963:               PetscMalloc1(dim, &qpoints);
1964:               PetscMalloc1(1,   &qweights);
1965:               PetscQuadratureSetOrder(sp->functional[f], 0);
1966:               PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1967:               TensorPoint_Internal(dim, orderEff+1, ind, tup);
1968:               for (d = 0; d < dim; ++d) {
1969:                 qpoints[d] = v0[d];
1970:                 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d]+1)*dx);
1971:               }
1972:               qweights[0] = 1.0;
1973:               ++f;
1974:             }
1975:           }
1976:           PetscFree5(ind,tup,v0,J,invJ);
1977:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &csize, &coords);
1978:           lag->numDof[depth] = cdim;
1979:         }
1980:       }
1981:       DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);
1982:     }
1983:     PetscBTDestroy(&seen);
1984:   }
1985:   if (pEnd[dim] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d not equal to dimension %d", f, pdimMax);
1986:   PetscFree2(pStart,pEnd);
1987:   if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
1988:   lag->height    = 0;
1989:   lag->subspaces = NULL;
1990:   if (lag->continuous && sp->order > 0 && dim > 0) {
1991:     PetscInt i;

1993:     lag->height = dim;
1994:     PetscMalloc1(dim,&lag->subspaces);
1995:     PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);
1996:     PetscDualSpaceSetUp(lag->subspaces[0]);
1997:     for (i = 1; i < dim; i++) {
1998:       PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);
1999:       PetscObjectReference((PetscObject)(lag->subspaces[i]));
2000:     }
2001:   }
2002:   return(0);
2003: }

2007: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
2008: {
2009:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2010:   PetscInt            i;
2011:   PetscErrorCode      ierr;

2014:   for (i = 0; i < lag->height; i++) {
2015:     PetscDualSpaceDestroy(&lag->subspaces[i]);
2016:   }
2017:   PetscFree(lag->subspaces);
2018:   PetscFree(lag->numDof);
2019:   PetscFree(lag);
2020:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
2021:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
2022:   return(0);
2023: }

2027: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
2028: {
2029:   PetscInt       order;
2030:   PetscBool      cont;

2034:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2035:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
2036:   PetscDualSpaceGetOrder(sp, &order);
2037:   PetscDualSpaceSetOrder(*spNew, order);
2038:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
2039:   PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
2040:   return(0);
2041: }

2045: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
2046: {
2047:   PetscBool      continuous, flg;

2051:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
2052:   PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
2053:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
2054:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
2055:   PetscOptionsTail();
2056:   return(0);
2057: }

2061: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
2062: {
2063:   DM              K;
2064:   const PetscInt *numDof;
2065:   PetscInt        spatialDim, Nc, size = 0, d;
2066:   PetscErrorCode  ierr;

2069:   PetscDualSpaceGetDM(sp, &K);
2070:   PetscDualSpaceGetNumDof(sp, &numDof);
2071:   DMGetDimension(K, &spatialDim);
2072:   DMPlexGetHeightStratum(K, 0, NULL, &Nc);
2073:   if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
2074:   for (d = 0; d <= spatialDim; ++d) {
2075:     PetscInt pStart, pEnd;

2077:     DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
2078:     size += (pEnd-pStart)*numDof[d];
2079:   }
2080:   *dim = size;
2081:   return(0);
2082: }

2086: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
2087: {
2088:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2091:   *numDof = lag->numDof;
2092:   return(0);
2093: }

2097: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2098: {
2099:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2104:   *continuous = lag->continuous;
2105:   return(0);
2106: }

2110: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2111: {
2112:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2116:   lag->continuous = continuous;
2117:   return(0);
2118: }

2122: /*@
2123:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2125:   Not Collective

2127:   Input Parameter:
2128: . sp         - the PetscDualSpace

2130:   Output Parameter:
2131: . continuous - flag for element continuity

2133:   Level: intermediate

2135: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2136: .seealso: PetscDualSpaceLagrangeSetContinuity()
2137: @*/
2138: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2139: {

2145:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2146:   return(0);
2147: }

2151: /*@
2152:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2154:   Logically Collective on PetscDualSpace

2156:   Input Parameters:
2157: + sp         - the PetscDualSpace
2158: - continuous - flag for element continuity

2160:   Options Database:
2161: . -petscdualspace_lagrange_continuity <bool>

2163:   Level: intermediate

2165: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2166: .seealso: PetscDualSpaceLagrangeGetContinuity()
2167: @*/
2168: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2169: {

2175:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2176:   return(0);
2177: }

2181: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
2182: {
2183:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2184:   PetscErrorCode     ierr;

2189:   if (height == 0) {
2190:     *bdsp = sp;
2191:   }
2192:   else {
2193:     DM dm;
2194:     PetscInt dim;

2196:     PetscDualSpaceGetDM(sp,&dm);
2197:     DMGetDimension(dm,&dim);
2198:     if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim);}
2199:     if (height <= lag->height) {
2200:       *bdsp = lag->subspaces[height-1];
2201:     }
2202:     else {
2203:       *bdsp = NULL;
2204:     }
2205:   }
2206:   return(0);
2207: }

2211: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2212: {
2214:   sp->ops->setfromoptions    = PetscDualSpaceSetFromOptions_Lagrange;
2215:   sp->ops->setup             = PetscDualSpaceSetUp_Lagrange;
2216:   sp->ops->view              = NULL;
2217:   sp->ops->destroy           = PetscDualSpaceDestroy_Lagrange;
2218:   sp->ops->duplicate         = PetscDualSpaceDuplicate_Lagrange;
2219:   sp->ops->getdimension      = PetscDualSpaceGetDimension_Lagrange;
2220:   sp->ops->getnumdof         = PetscDualSpaceGetNumDof_Lagrange;
2221:   sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
2222:   return(0);
2223: }

2225: /*MC
2226:   PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals

2228:   Level: intermediate

2230: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2231: M*/

2235: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2236: {
2237:   PetscDualSpace_Lag *lag;
2238:   PetscErrorCode      ierr;

2242:   PetscNewLog(sp,&lag);
2243:   sp->data = lag;

2245:   lag->numDof     = NULL;
2246:   lag->simplex    = PETSC_TRUE;
2247:   lag->continuous = PETSC_TRUE;

2249:   PetscDualSpaceInitialize_Lagrange(sp);
2250:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2251:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2252:   return(0);
2253: }

2257: PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
2258: {
2259:   PetscDualSpace_Simple *s  = (PetscDualSpace_Simple *) sp->data;
2260:   DM                     dm = sp->dm;
2261:   PetscInt               dim;
2262:   PetscErrorCode         ierr;

2265:   DMGetDimension(dm, &dim);
2266:   PetscCalloc1(dim+1, &s->numDof);
2267:   return(0);
2268: }

2272: PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
2273: {
2274:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2275:   PetscErrorCode         ierr;

2278:   PetscFree(s->numDof);
2279:   PetscFree(s);
2280:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
2281:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
2282:   return(0);
2283: }

2287: PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace *spNew)
2288: {
2289:   PetscInt       dim, d;

2293:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2294:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACESIMPLE);
2295:   PetscDualSpaceGetDimension(sp, &dim);
2296:   PetscDualSpaceSimpleSetDimension(*spNew, dim);
2297:   for (d = 0; d < dim; ++d) {
2298:     PetscQuadrature q;

2300:     PetscDualSpaceGetFunctional(sp, d, &q);
2301:     PetscDualSpaceSimpleSetFunctional(*spNew, d, q);
2302:   }
2303:   return(0);
2304: }

2308: PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
2309: {
2311:   return(0);
2312: }

2316: PetscErrorCode PetscDualSpaceGetDimension_Simple(PetscDualSpace sp, PetscInt *dim)
2317: {
2318:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

2321:   *dim = s->dim;
2322:   return(0);
2323: }

2327: PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
2328: {
2329:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2330:   DM                     dm;
2331:   PetscInt               spatialDim, f;
2332:   PetscErrorCode         ierr;

2335:   for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
2336:   PetscFree(sp->functional);
2337:   s->dim = dim;
2338:   PetscCalloc1(s->dim, &sp->functional);
2339:   PetscFree(s->numDof);
2340:   PetscDualSpaceGetDM(sp, &dm);
2341:   DMGetCoordinateDim(dm, &spatialDim);
2342:   PetscCalloc1(spatialDim+1, &s->numDof);
2343:   s->numDof[spatialDim] = dim;
2344:   return(0);
2345: }

2349: PetscErrorCode PetscDualSpaceGetNumDof_Simple(PetscDualSpace sp, const PetscInt **numDof)
2350: {
2351:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

2354:   *numDof = s->numDof;
2355:   return(0);
2356: }

2360: PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
2361: {
2362:   PetscDualSpace_Simple *s   = (PetscDualSpace_Simple *) sp->data;
2363:   PetscReal              vol = 0.0;
2364:   PetscReal             *weights;
2365:   PetscInt               Nq, p;
2366:   PetscErrorCode         ierr;

2369:   if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim);
2370:   PetscQuadratureDuplicate(q, &sp->functional[f]);
2371:   /* Reweight so that it has unit volume */
2372:   PetscQuadratureGetData(sp->functional[f], NULL, &Nq, NULL, (const PetscReal **) &weights);
2373:   for (p = 0; p < Nq; ++p) vol += weights[p];
2374:   for (p = 0; p < Nq; ++p) weights[p] /= vol;
2375:   return(0);
2376: }

2380: /*@
2381:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

2383:   Logically Collective on PetscDualSpace

2385:   Input Parameters:
2386: + sp  - the PetscDualSpace
2387: - dim - the basis dimension

2389:   Level: intermediate

2391: .keywords: PetscDualSpace, dimension
2392: .seealso: PetscDualSpaceSimpleSetFunctional()
2393: @*/
2394: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
2395: {

2401:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
2402:   return(0);
2403: }

2407: /*@
2408:   PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space

2410:   Not Collective

2412:   Input Parameters:
2413: + sp  - the PetscDualSpace
2414: . f - the basis index
2415: - q - the basis functional

2417:   Level: intermediate

2419:   Note: The quadrature will be reweighted so that it has unit volume.

2421: .keywords: PetscDualSpace, functional
2422: .seealso: PetscDualSpaceSimpleSetDimension()
2423: @*/
2424: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
2425: {

2430:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
2431:   return(0);
2432: }

2436: PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
2437: {
2439:   sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple;
2440:   sp->ops->setup          = PetscDualSpaceSetUp_Simple;
2441:   sp->ops->view           = NULL;
2442:   sp->ops->destroy        = PetscDualSpaceDestroy_Simple;
2443:   sp->ops->duplicate      = PetscDualSpaceDuplicate_Simple;
2444:   sp->ops->getdimension   = PetscDualSpaceGetDimension_Simple;
2445:   sp->ops->getnumdof      = PetscDualSpaceGetNumDof_Simple;
2446:   return(0);
2447: }

2449: /*MC
2450:   PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals

2452:   Level: intermediate

2454: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2455: M*/

2459: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
2460: {
2461:   PetscDualSpace_Simple *s;
2462:   PetscErrorCode         ierr;

2466:   PetscNewLog(sp,&s);
2467:   sp->data = s;

2469:   s->dim    = 0;
2470:   s->numDof = NULL;

2472:   PetscDualSpaceInitialize_Simple(sp);
2473:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
2474:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
2475:   return(0);
2476: }


2479: PetscClassId PETSCFE_CLASSID = 0;

2481: PetscFunctionList PetscFEList              = NULL;
2482: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

2486: /*@C
2487:   PetscFERegister - Adds a new PetscFE implementation

2489:   Not Collective

2491:   Input Parameters:
2492: + name        - The name of a new user-defined creation routine
2493: - create_func - The creation routine itself

2495:   Notes:
2496:   PetscFERegister() may be called multiple times to add several user-defined PetscFEs

2498:   Sample usage:
2499: .vb
2500:     PetscFERegister("my_fe", MyPetscFECreate);
2501: .ve

2503:   Then, your PetscFE type can be chosen with the procedural interface via
2504: .vb
2505:     PetscFECreate(MPI_Comm, PetscFE *);
2506:     PetscFESetType(PetscFE, "my_fe");
2507: .ve
2508:    or at runtime via the option
2509: .vb
2510:     -petscfe_type my_fe
2511: .ve

2513:   Level: advanced

2515: .keywords: PetscFE, register
2516: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()

2518: @*/
2519: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2520: {

2524:   PetscFunctionListAdd(&PetscFEList, sname, function);
2525:   return(0);
2526: }

2530: /*@C
2531:   PetscFESetType - Builds a particular PetscFE

2533:   Collective on PetscFE

2535:   Input Parameters:
2536: + fem  - The PetscFE object
2537: - name - The kind of FEM space

2539:   Options Database Key:
2540: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types

2542:   Level: intermediate

2544: .keywords: PetscFE, set, type
2545: .seealso: PetscFEGetType(), PetscFECreate()
2546: @*/
2547: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2548: {
2549:   PetscErrorCode (*r)(PetscFE);
2550:   PetscBool      match;

2555:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
2556:   if (match) return(0);

2558:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
2559:   PetscFunctionListFind(PetscFEList, name, &r);
2560:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);

2562:   if (fem->ops->destroy) {
2563:     (*fem->ops->destroy)(fem);
2564:     fem->ops->destroy = NULL;
2565:   }
2566:   (*r)(fem);
2567:   PetscObjectChangeTypeName((PetscObject) fem, name);
2568:   return(0);
2569: }

2573: /*@C
2574:   PetscFEGetType - Gets the PetscFE type name (as a string) from the object.

2576:   Not Collective

2578:   Input Parameter:
2579: . fem  - The PetscFE

2581:   Output Parameter:
2582: . name - The PetscFE type name

2584:   Level: intermediate

2586: .keywords: PetscFE, get, type, name
2587: .seealso: PetscFESetType(), PetscFECreate()
2588: @*/
2589: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2590: {

2596:   if (!PetscFERegisterAllCalled) {
2597:     PetscFERegisterAll();
2598:   }
2599:   *name = ((PetscObject) fem)->type_name;
2600:   return(0);
2601: }

2605: /*@C
2606:   PetscFEView - Views a PetscFE

2608:   Collective on PetscFE

2610:   Input Parameter:
2611: + fem - the PetscFE object to view
2612: - v   - the viewer

2614:   Level: developer

2616: .seealso PetscFEDestroy()
2617: @*/
2618: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
2619: {

2624:   if (!v) {
2625:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
2626:   }
2627:   if (fem->ops->view) {
2628:     (*fem->ops->view)(fem, v);
2629:   }
2630:   return(0);
2631: }

2635: /*@
2636:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

2638:   Collective on PetscFE

2640:   Input Parameter:
2641: . fem - the PetscFE object to set options for

2643:   Options Database:
2644: . -petscfe_num_blocks  the number of cell blocks to integrate concurrently
2645: . -petscfe_num_batches the number of cell batches to integrate serially

2647:   Level: developer

2649: .seealso PetscFEView()
2650: @*/
2651: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
2652: {
2653:   const char    *defaultType;
2654:   char           name[256];
2655:   PetscBool      flg;

2660:   if (!((PetscObject) fem)->type_name) {
2661:     defaultType = PETSCFEBASIC;
2662:   } else {
2663:     defaultType = ((PetscObject) fem)->type_name;
2664:   }
2665:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

2667:   PetscObjectOptionsBegin((PetscObject) fem);
2668:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
2669:   if (flg) {
2670:     PetscFESetType(fem, name);
2671:   } else if (!((PetscObject) fem)->type_name) {
2672:     PetscFESetType(fem, defaultType);
2673:   }
2674:   PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
2675:   PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
2676:   if (fem->ops->setfromoptions) {
2677:     (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
2678:   }
2679:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
2680:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
2681:   PetscOptionsEnd();
2682:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
2683:   return(0);
2684: }

2688: /*@C
2689:   PetscFESetUp - Construct data structures for the PetscFE

2691:   Collective on PetscFE

2693:   Input Parameter:
2694: . fem - the PetscFE object to setup

2696:   Level: developer

2698: .seealso PetscFEView(), PetscFEDestroy()
2699: @*/
2700: PetscErrorCode PetscFESetUp(PetscFE fem)
2701: {

2706:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
2707:   return(0);
2708: }

2712: /*@
2713:   PetscFEDestroy - Destroys a PetscFE object

2715:   Collective on PetscFE

2717:   Input Parameter:
2718: . fem - the PetscFE object to destroy

2720:   Level: developer

2722: .seealso PetscFEView()
2723: @*/
2724: PetscErrorCode PetscFEDestroy(PetscFE *fem)
2725: {

2729:   if (!*fem) return(0);

2732:   if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
2733:   ((PetscObject) (*fem))->refct = 0;

2735:   PetscFree((*fem)->numDof);
2736:   PetscFree((*fem)->invV);
2737:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
2738:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
2739:   PetscSpaceDestroy(&(*fem)->basisSpace);
2740:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
2741:   PetscQuadratureDestroy(&(*fem)->quadrature);

2743:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
2744:   PetscHeaderDestroy(fem);
2745:   return(0);
2746: }

2750: /*@
2751:   PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().

2753:   Collective on MPI_Comm

2755:   Input Parameter:
2756: . comm - The communicator for the PetscFE object

2758:   Output Parameter:
2759: . fem - The PetscFE object

2761:   Level: beginner

2763: .seealso: PetscFESetType(), PETSCFEGALERKIN
2764: @*/
2765: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
2766: {
2767:   PetscFE        f;

2772:   PetscCitationsRegister(FECitation,&FEcite);
2773:   *fem = NULL;
2774:   PetscFEInitializePackage();

2776:   PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);

2778:   f->basisSpace    = NULL;
2779:   f->dualSpace     = NULL;
2780:   f->numComponents = 1;
2781:   f->numDof        = NULL;
2782:   f->invV          = NULL;
2783:   f->B             = NULL;
2784:   f->D             = NULL;
2785:   f->H             = NULL;
2786:   PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
2787:   f->blockSize     = 0;
2788:   f->numBlocks     = 1;
2789:   f->batchSize     = 0;
2790:   f->numBatches    = 1;

2792:   *fem = f;
2793:   return(0);
2794: }

2798: /*@
2799:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

2801:   Not collective

2803:   Input Parameter:
2804: . fem - The PetscFE object

2806:   Output Parameter:
2807: . dim - The spatial dimension

2809:   Level: intermediate

2811: .seealso: PetscFECreate()
2812: @*/
2813: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
2814: {
2815:   DM             dm;

2821:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2822:   DMGetDimension(dm, dim);
2823:   return(0);
2824: }

2828: /*@
2829:   PetscFESetNumComponents - Sets the number of components in the element

2831:   Not collective

2833:   Input Parameters:
2834: + fem - The PetscFE object
2835: - comp - The number of field components

2837:   Level: intermediate

2839: .seealso: PetscFECreate()
2840: @*/
2841: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
2842: {
2845:   fem->numComponents = comp;
2846:   return(0);
2847: }

2851: /*@
2852:   PetscFEGetNumComponents - Returns the number of components in the element

2854:   Not collective

2856:   Input Parameter:
2857: . fem - The PetscFE object

2859:   Output Parameter:
2860: . comp - The number of field components

2862:   Level: intermediate

2864: .seealso: PetscFECreate()
2865: @*/
2866: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
2867: {
2871:   *comp = fem->numComponents;
2872:   return(0);
2873: }

2877: /*@
2878:   PetscFESetTileSizes - Sets the tile sizes for evaluation

2880:   Not collective

2882:   Input Parameters:
2883: + fem - The PetscFE object
2884: . blockSize - The number of elements in a block
2885: . numBlocks - The number of blocks in a batch
2886: . batchSize - The number of elements in a batch
2887: - numBatches - The number of batches in a chunk

2889:   Level: intermediate

2891: .seealso: PetscFECreate()
2892: @*/
2893: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
2894: {
2897:   fem->blockSize  = blockSize;
2898:   fem->numBlocks  = numBlocks;
2899:   fem->batchSize  = batchSize;
2900:   fem->numBatches = numBatches;
2901:   return(0);
2902: }

2906: /*@
2907:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

2909:   Not collective

2911:   Input Parameter:
2912: . fem - The PetscFE object

2914:   Output Parameters:
2915: + blockSize - The number of elements in a block
2916: . numBlocks - The number of blocks in a batch
2917: . batchSize - The number of elements in a batch
2918: - numBatches - The number of batches in a chunk

2920:   Level: intermediate

2922: .seealso: PetscFECreate()
2923: @*/
2924: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
2925: {
2932:   if (blockSize)  *blockSize  = fem->blockSize;
2933:   if (numBlocks)  *numBlocks  = fem->numBlocks;
2934:   if (batchSize)  *batchSize  = fem->batchSize;
2935:   if (numBatches) *numBatches = fem->numBatches;
2936:   return(0);
2937: }

2941: /*@
2942:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

2944:   Not collective

2946:   Input Parameter:
2947: . fem - The PetscFE object

2949:   Output Parameter:
2950: . sp - The PetscSpace object

2952:   Level: intermediate

2954: .seealso: PetscFECreate()
2955: @*/
2956: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
2957: {
2961:   *sp = fem->basisSpace;
2962:   return(0);
2963: }

2967: /*@
2968:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

2970:   Not collective

2972:   Input Parameters:
2973: + fem - The PetscFE object
2974: - sp - The PetscSpace object

2976:   Level: intermediate

2978: .seealso: PetscFECreate()
2979: @*/
2980: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
2981: {

2987:   PetscSpaceDestroy(&fem->basisSpace);
2988:   fem->basisSpace = sp;
2989:   PetscObjectReference((PetscObject) fem->basisSpace);
2990:   return(0);
2991: }

2995: /*@
2996:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

2998:   Not collective

3000:   Input Parameter:
3001: . fem - The PetscFE object

3003:   Output Parameter:
3004: . sp - The PetscDualSpace object

3006:   Level: intermediate

3008: .seealso: PetscFECreate()
3009: @*/
3010: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
3011: {
3015:   *sp = fem->dualSpace;
3016:   return(0);
3017: }

3021: /*@
3022:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

3024:   Not collective

3026:   Input Parameters:
3027: + fem - The PetscFE object
3028: - sp - The PetscDualSpace object

3030:   Level: intermediate

3032: .seealso: PetscFECreate()
3033: @*/
3034: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
3035: {

3041:   PetscDualSpaceDestroy(&fem->dualSpace);
3042:   fem->dualSpace = sp;
3043:   PetscObjectReference((PetscObject) fem->dualSpace);
3044:   return(0);
3045: }

3049: /*@
3050:   PetscFEGetQuadrature - Returns the PetscQuadreture used to calculate inner products

3052:   Not collective

3054:   Input Parameter:
3055: . fem - The PetscFE object

3057:   Output Parameter:
3058: . q - The PetscQuadrature object

3060:   Level: intermediate

3062: .seealso: PetscFECreate()
3063: @*/
3064: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
3065: {
3069:   *q = fem->quadrature;
3070:   return(0);
3071: }

3075: /*@
3076:   PetscFESetQuadrature - Sets the PetscQuadreture used to calculate inner products

3078:   Not collective

3080:   Input Parameters:
3081: + fem - The PetscFE object
3082: - q - The PetscQuadrature object

3084:   Level: intermediate

3086: .seealso: PetscFECreate()
3087: @*/
3088: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
3089: {

3094:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
3095:   PetscQuadratureDestroy(&fem->quadrature);
3096:   fem->quadrature = q;
3097:   PetscObjectReference((PetscObject) q);
3098:   return(0);
3099: }

3103: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
3104: {
3105:   const PetscInt *numDofDual;
3106:   PetscErrorCode  ierr;

3111:   PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
3112:   if (!fem->numDof) {
3113:     DM       dm;
3114:     PetscInt dim, d;

3116:     PetscDualSpaceGetDM(fem->dualSpace, &dm);
3117:     DMGetDimension(dm, &dim);
3118:     PetscMalloc1(dim+1, &fem->numDof);
3119:     for (d = 0; d <= dim; ++d) {
3120:       fem->numDof[d] = fem->numComponents*numDofDual[d];
3121:     }
3122:   }
3123:   *numDof = fem->numDof;
3124:   return(0);
3125: }

3129: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
3130: {
3131:   PetscInt         npoints;
3132:   const PetscReal *points;
3133:   PetscErrorCode   ierr;

3140:   PetscQuadratureGetData(fem->quadrature, NULL, &npoints, &points, NULL);
3141:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
3142:   if (B) *B = fem->B;
3143:   if (D) *D = fem->D;
3144:   if (H) *H = fem->H;
3145:   return(0);
3146: }

3150: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **F)
3151: {
3152:   PetscErrorCode   ierr;

3157:   if (!fem->F) {
3158:     PetscDualSpace  sp;
3159:     DM              dm;
3160:     const PetscInt *cone;
3161:     PetscReal      *centroids;
3162:     PetscInt        dim, numFaces, f;

3164:     PetscFEGetDualSpace(fem, &sp);
3165:     PetscDualSpaceGetDM(sp, &dm);
3166:     DMGetDimension(dm, &dim);
3167:     DMPlexGetConeSize(dm, 0, &numFaces);
3168:     DMPlexGetCone(dm, 0, &cone);
3169:     PetscMalloc1(numFaces*dim, &centroids);
3170:     for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f*dim], NULL);}
3171:     PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
3172:     PetscFree(centroids);
3173:   }
3174:   *F = fem->F;
3175:   return(0);
3176: }

3180: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3181: {
3182:   DM               dm;
3183:   PetscInt         pdim; /* Dimension of FE space P */
3184:   PetscInt         dim;  /* Spatial dimension */
3185:   PetscInt         comp; /* Field components */
3186:   PetscErrorCode   ierr;

3194:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3195:   DMGetDimension(dm, &dim);
3196:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3197:   PetscFEGetNumComponents(fem, &comp);
3198:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
3199:   if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
3200:   if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
3201:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
3202:   return(0);
3203: }

3207: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3208: {
3209:   DM             dm;

3214:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3215:   if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
3216:   if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
3217:   if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
3218:   return(0);
3219: }

3223: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
3224: {
3225:   PetscFE_Basic *b = (PetscFE_Basic *) fem->data;

3229:   PetscFree(b);
3230:   return(0);
3231: }

3235: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
3236: {
3237:   PetscSpace        basis;
3238:   PetscDualSpace    dual;
3239:   PetscQuadrature   q = NULL;
3240:   PetscInt          dim, Nq;
3241:   PetscViewerFormat format;
3242:   PetscErrorCode    ierr;

3245:   PetscFEGetBasisSpace(fe, &basis);
3246:   PetscFEGetDualSpace(fe, &dual);
3247:   PetscFEGetQuadrature(fe, &q);
3248:   PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
3249:   PetscViewerGetFormat(viewer, &format);
3250:   PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
3251:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
3252:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
3253:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
3254:     PetscViewerASCIIPushTab(viewer);
3255:     PetscQuadratureView(q, viewer);
3256:     PetscViewerASCIIPopTab(viewer);
3257:   } else {
3258:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
3259:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
3260:   }
3261:   PetscViewerASCIIPushTab(viewer);
3262:   PetscSpaceView(basis, viewer);
3263:   PetscDualSpaceView(dual, viewer);
3264:   PetscViewerASCIIPopTab(viewer);
3265:   return(0);
3266: }

3270: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
3271: {
3272:   PetscBool      iascii;

3278:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
3279:   if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
3280:   return(0);
3281: }

3285: /* Construct the change of basis from prime basis to nodal basis */
3286: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
3287: {
3288:   PetscScalar   *work, *invVscalar;
3289:   PetscBLASInt  *pivots;
3290:   PetscBLASInt   n, info;
3291:   PetscInt       pdim, j;

3295:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3296:   PetscMalloc1(pdim*pdim,&fem->invV);
3297: #if defined(PETSC_USE_COMPLEX)
3298:   PetscMalloc1(pdim*pdim,&invVscalar);
3299: #else
3300:   invVscalar = fem->invV;
3301: #endif
3302:   for (j = 0; j < pdim; ++j) {
3303:     PetscReal      *Bf;
3304:     PetscQuadrature f;
3305:     PetscInt        q, k;

3307:     PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
3308:     PetscMalloc1(f->numPoints*pdim,&Bf);
3309:     PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
3310:     for (k = 0; k < pdim; ++k) {
3311:       /* n_j \cdot \phi_k */
3312:       invVscalar[j*pdim+k] = 0.0;
3313:       for (q = 0; q < f->numPoints; ++q) {
3314:         invVscalar[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
3315:       }
3316:     }
3317:     PetscFree(Bf);
3318:   }
3319:   PetscMalloc2(pdim,&pivots,pdim,&work);
3320:   n = pdim;
3321:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info));
3322:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info));
3323: #if defined(PETSC_USE_COMPLEX)
3324:   for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]);
3325:   PetscFree(invVscalar);
3326: #endif
3327:   PetscFree2(pivots,work);
3328:   return(0);
3329: }

3333: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
3334: {

3338:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
3339:   return(0);
3340: }

3344: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
3345: {
3346:   DM               dm;
3347:   PetscInt         pdim; /* Dimension of FE space P */
3348:   PetscInt         dim;  /* Spatial dimension */
3349:   PetscInt         comp; /* Field components */
3350:   PetscReal       *tmpB, *tmpD, *tmpH;
3351:   PetscInt         p, d, j, k;
3352:   PetscErrorCode   ierr;

3355:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3356:   DMGetDimension(dm, &dim);
3357:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3358:   PetscFEGetNumComponents(fem, &comp);
3359:   /* Evaluate the prime basis functions at all points */
3360:   if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3361:   if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3362:   if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3363:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
3364:   /* Translate to the nodal basis */
3365:   for (p = 0; p < npoints; ++p) {
3366:     if (B) {
3367:       /* Multiply by V^{-1} (pdim x pdim) */
3368:       for (j = 0; j < pdim; ++j) {
3369:         const PetscInt i = (p*pdim + j)*comp;
3370:         PetscInt       c;

3372:         B[i] = 0.0;
3373:         for (k = 0; k < pdim; ++k) {
3374:           B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
3375:         }
3376:         for (c = 1; c < comp; ++c) {
3377:           B[i+c] = B[i];
3378:         }
3379:       }
3380:     }
3381:     if (D) {
3382:       /* Multiply by V^{-1} (pdim x pdim) */
3383:       for (j = 0; j < pdim; ++j) {
3384:         for (d = 0; d < dim; ++d) {
3385:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
3386:           PetscInt       c;

3388:           D[i] = 0.0;
3389:           for (k = 0; k < pdim; ++k) {
3390:             D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
3391:           }
3392:           for (c = 1; c < comp; ++c) {
3393:             D[((p*pdim + j)*comp + c)*dim + d] = D[i];
3394:           }
3395:         }
3396:       }
3397:     }
3398:     if (H) {
3399:       /* Multiply by V^{-1} (pdim x pdim) */
3400:       for (j = 0; j < pdim; ++j) {
3401:         for (d = 0; d < dim*dim; ++d) {
3402:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
3403:           PetscInt       c;

3405:           H[i] = 0.0;
3406:           for (k = 0; k < pdim; ++k) {
3407:             H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
3408:           }
3409:           for (c = 1; c < comp; ++c) {
3410:             H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
3411:           }
3412:         }
3413:       }
3414:     }
3415:   }
3416:   if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3417:   if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3418:   if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3419:   return(0);
3420: }

3424: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3425:                                       const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
3426: {
3427:   const PetscInt  debug = 0;
3428:   PetscPointFunc  obj_func;
3429:   PetscQuadrature quad;
3430:   PetscScalar    *u, *u_x, *a, *a_x;
3431:   PetscReal      *x;
3432:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3433:   PetscInt        dim, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
3434:   PetscErrorCode  ierr;

3437:   PetscDSGetObjective(prob, field, &obj_func);
3438:   if (!obj_func) return(0);
3439:   PetscFEGetSpatialDimension(fem, &dim);
3440:   PetscFEGetQuadrature(fem, &quad);
3441:   PetscDSGetNumFields(prob, &Nf);
3442:   PetscDSGetTotalDimension(prob, &totDim);
3443:   PetscDSGetComponentOffsets(prob, &uOff);
3444:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3445:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
3446:   PetscDSGetRefCoordArrays(prob, &x, NULL);
3447:   if (probAux) {
3448:     PetscDSGetNumFields(probAux, &NfAux);
3449:     PetscDSGetTotalDimension(probAux, &totDimAux);
3450:     PetscDSGetComponentOffsets(probAux, &aOff);
3451:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3452:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3453:   }
3454:   for (e = 0; e < Ne; ++e) {
3455:     const PetscReal *v0   = geom[e].v0;
3456:     const PetscReal *J    = geom[e].J;
3457:     const PetscReal *invJ = geom[e].invJ;
3458:     const PetscReal  detJ = geom[e].detJ;
3459:     const PetscReal *quadPoints, *quadWeights;
3460:     PetscInt         Nq, q;

3462:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3463:     if (debug > 1) {
3464:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3465: #ifndef PETSC_USE_COMPLEX
3466:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3467: #endif
3468:     }
3469:     for (q = 0; q < Nq; ++q) {
3470:       PetscScalar integrand;

3472:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3473:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3474:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       NULL, u, u_x, NULL);
3475:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3476:       obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &integrand);
3477:       integrand *= detJ*quadWeights[q];
3478:       integral[field] += PetscRealPart(integrand);
3479:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", PetscRealPart(integrand), integral[field]);}
3480:     }
3481:     cOffset    += totDim;
3482:     cOffsetAux += totDimAux;
3483:   }
3484:   return(0);
3485: }

3489: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3490:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
3491: {
3492:   const PetscInt  debug = 0;
3493:   PetscPointFunc  f0_func;
3494:   PetscPointFunc  f1_func;
3495:   PetscQuadrature quad;
3496:   PetscReal     **basisField, **basisFieldDer;
3497:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3498:   PetscReal      *x;
3499:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3500:   PetscInt        dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3501:   PetscErrorCode  ierr;

3504:   PetscFEGetSpatialDimension(fem, &dim);
3505:   PetscFEGetQuadrature(fem, &quad);
3506:   PetscFEGetDimension(fem, &Nb);
3507:   PetscFEGetNumComponents(fem, &Nc);
3508:   PetscDSGetNumFields(prob, &Nf);
3509:   PetscDSGetTotalDimension(prob, &totDim);
3510:   PetscDSGetComponentOffsets(prob, &uOff);
3511:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3512:   PetscDSGetFieldOffset(prob, field, &fOffset);
3513:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3514:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3515:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3516:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3517:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3518:   if (probAux) {
3519:     PetscDSGetNumFields(probAux, &NfAux);
3520:     PetscDSGetTotalDimension(probAux, &totDimAux);
3521:     PetscDSGetComponentOffsets(probAux, &aOff);
3522:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3523:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3524:   }
3525:   for (e = 0; e < Ne; ++e) {
3526:     const PetscReal *v0   = geom[e].v0;
3527:     const PetscReal *J    = geom[e].J;
3528:     const PetscReal *invJ = geom[e].invJ;
3529:     const PetscReal  detJ = geom[e].detJ;
3530:     const PetscReal *quadPoints, *quadWeights;
3531:     PetscInt         Nq, q;

3533:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3534:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3535:     PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
3536:     if (debug > 1) {
3537:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3538: #ifndef PETSC_USE_COMPLEX
3539:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3540: #endif
3541:     }
3542:     for (q = 0; q < Nq; ++q) {
3543:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3544:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3545:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3546:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3547:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, &f0[q*Nc]);
3548:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, refSpaceDer);
3549:       TransformF(dim, dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3550:     }
3551:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3552:     cOffset    += totDim;
3553:     cOffsetAux += totDimAux;
3554:   }
3555:   return(0);
3556: }

3560: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3561:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
3562: {
3563:   const PetscInt  debug = 0;
3564:   PetscBdPointFunc f0_func;
3565:   PetscBdPointFunc f1_func;
3566:   PetscQuadrature quad;
3567:   PetscReal     **basisField, **basisFieldDer;
3568:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3569:   PetscReal      *x;
3570:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3571:   PetscInt        dim, Nf, NfAux = 0, bdim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3572:   PetscErrorCode  ierr;

3575:   PetscFEGetSpatialDimension(fem, &dim);
3576:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3577:   bdim = dim-1;
3578:   PetscFEGetQuadrature(fem, &quad);
3579:   PetscFEGetDimension(fem, &Nb);
3580:   PetscFEGetNumComponents(fem, &Nc);
3581:   PetscDSGetNumFields(prob, &Nf);
3582:   PetscDSGetTotalBdDimension(prob, &totDim);
3583:   PetscDSGetComponentBdOffsets(prob, &uOff);
3584:   PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
3585:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3586:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3587:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3588:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3589:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3590:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3591:   if (probAux) {
3592:     PetscDSGetNumFields(probAux, &NfAux);
3593:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3594:     PetscDSGetComponentBdOffsets(probAux, &aOff);
3595:     PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
3596:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3597:   }
3598:   for (e = 0; e < Ne; ++e) {
3599:     const PetscReal *v0          = geom[e].v0;
3600:     const PetscReal *n           = geom[e].n;
3601:     const PetscReal *J           = geom[e].J;
3602:     const PetscReal *invJ        = geom[e].invJ;
3603:     const PetscReal  detJ        = geom[e].detJ;
3604:     const PetscReal *quadPoints, *quadWeights;
3605:     PetscInt         Nq, q;

3607:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3608:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3609:     PetscMemzero(f1, Nq*Nc*bdim * sizeof(PetscScalar));
3610:     if (debug > 1) {
3611:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3612: #ifndef PETSC_USE_COMPLEX
3613:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3614: #endif
3615:     }
3616:     for (q = 0; q < Nq; ++q) {
3617:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3618:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3619:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3620:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3621:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, &f0[q*Nc]);
3622:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, refSpaceDer);
3623:       TransformF(dim, bdim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3624:     }
3625:     UpdateElementVec(bdim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3626:     cOffset    += totDim;
3627:     cOffsetAux += totDimAux;
3628:   }
3629:   return(0);
3630: }

3634: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3635:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
3636: {
3637:   const PetscInt  debug      = 0;
3638:   PetscPointJac   g0_func;
3639:   PetscPointJac   g1_func;
3640:   PetscPointJac   g2_func;
3641:   PetscPointJac   g3_func;
3642:   PetscFE         fe;
3643:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3644:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3645:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3646:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3647:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3648:   PetscQuadrature quad;
3649:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3650:   PetscReal      *x;
3651:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3652:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3653:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3654:   PetscInt        dim, Nf, NfAux = 0, totDim, totDimAux = 0, e;
3655:   PetscErrorCode  ierr;

3658:   PetscFEGetSpatialDimension(fem, &dim);
3659:   PetscFEGetQuadrature(fem, &quad);
3660:   PetscDSGetNumFields(prob, &Nf);
3661:   PetscDSGetTotalDimension(prob, &totDim);
3662:   PetscDSGetComponentOffsets(prob, &uOff);
3663:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3664:   switch(jtype) {
3665:   case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
3666:   case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
3667:   case PETSCFE_JACOBIAN:     PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
3668:   }
3669:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
3670:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3671:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3672:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3673:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3674:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3675:   PetscFEGetDimension(fe, &NbI);
3676:   PetscFEGetNumComponents(fe, &NcI);
3677:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3678:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3679:   PetscFEGetDimension(fe, &NbJ);
3680:   PetscFEGetNumComponents(fe, &NcJ);
3681:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3682:   if (probAux) {
3683:     PetscDSGetNumFields(probAux, &NfAux);
3684:     PetscDSGetTotalDimension(probAux, &totDimAux);
3685:     PetscDSGetComponentOffsets(probAux, &aOff);
3686:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3687:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3688:   }
3689:   basisI    = basisField[fieldI];
3690:   basisJ    = basisField[fieldJ];
3691:   basisDerI = basisFieldDer[fieldI];
3692:   basisDerJ = basisFieldDer[fieldJ];
3693:   /* Initialize here in case the function is not defined */
3694:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3695:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3696:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3697:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3698:   for (e = 0; e < Ne; ++e) {
3699:     const PetscReal *v0   = geom[e].v0;
3700:     const PetscReal *J    = geom[e].J;
3701:     const PetscReal *invJ = geom[e].invJ;
3702:     const PetscReal  detJ = geom[e].detJ;
3703:     const PetscReal *quadPoints, *quadWeights;
3704:     PetscInt         Nq, q;

3706:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3707:     for (q = 0; q < Nq; ++q) {
3708:       PetscInt f, g, fc, gc, c;

3710:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3711:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3712:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3713:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3714:       if (g0_func) {
3715:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3716:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, g0);
3717:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3718:       }
3719:       if (g1_func) {
3720:         PetscInt d, d2;
3721:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3722:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
3723:         for (fc = 0; fc < NcI; ++fc) {
3724:           for (gc = 0; gc < NcJ; ++gc) {
3725:             for (d = 0; d < dim; ++d) {
3726:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3727:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3728:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3729:             }
3730:           }
3731:         }
3732:       }
3733:       if (g2_func) {
3734:         PetscInt d, d2;
3735:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3736:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
3737:         for (fc = 0; fc < NcI; ++fc) {
3738:           for (gc = 0; gc < NcJ; ++gc) {
3739:             for (d = 0; d < dim; ++d) {
3740:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3741:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3742:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3743:             }
3744:           }
3745:         }
3746:       }
3747:       if (g3_func) {
3748:         PetscInt d, d2, dp, d3;
3749:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3750:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
3751:         for (fc = 0; fc < NcI; ++fc) {
3752:           for (gc = 0; gc < NcJ; ++gc) {
3753:             for (d = 0; d < dim; ++d) {
3754:               for (dp = 0; dp < dim; ++dp) {
3755:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3756:                 for (d2 = 0; d2 < dim; ++d2) {
3757:                   for (d3 = 0; d3 < dim; ++d3) {
3758:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3759:                   }
3760:                 }
3761:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3762:               }
3763:             }
3764:           }
3765:         }
3766:       }

3768:       for (f = 0; f < NbI; ++f) {
3769:         for (fc = 0; fc < NcI; ++fc) {
3770:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3771:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3772:           for (g = 0; g < NbJ; ++g) {
3773:             for (gc = 0; gc < NcJ; ++gc) {
3774:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3775:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3776:               const PetscInt fOff = eOffset+i*totDim+j;
3777:               PetscInt       d, d2;

3779:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3780:               for (d = 0; d < dim; ++d) {
3781:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3782:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3783:                 for (d2 = 0; d2 < dim; ++d2) {
3784:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3785:                 }
3786:               }
3787:             }
3788:           }
3789:         }
3790:       }
3791:     }
3792:     if (debug > 1) {
3793:       PetscInt fc, f, gc, g;

3795:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3796:       for (fc = 0; fc < NcI; ++fc) {
3797:         for (f = 0; f < NbI; ++f) {
3798:           const PetscInt i = offsetI + f*NcI+fc;
3799:           for (gc = 0; gc < NcJ; ++gc) {
3800:             for (g = 0; g < NbJ; ++g) {
3801:               const PetscInt j = offsetJ + g*NcJ+gc;
3802:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3803:             }
3804:           }
3805:           PetscPrintf(PETSC_COMM_SELF, "\n");
3806:         }
3807:       }
3808:     }
3809:     cOffset    += totDim;
3810:     cOffsetAux += totDimAux;
3811:     eOffset    += PetscSqr(totDim);
3812:   }
3813:   return(0);
3814: }

3818: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3819:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
3820: {
3821:   const PetscInt  debug      = 0;
3822:   PetscBdPointJac g0_func;
3823:   PetscBdPointJac g1_func;
3824:   PetscBdPointJac g2_func;
3825:   PetscBdPointJac g3_func;
3826:   PetscFE         fe;
3827:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3828:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3829:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3830:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3831:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3832:   PetscQuadrature quad;
3833:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3834:   PetscReal      *x;
3835:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3836:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3837:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3838:   PetscInt        dim, Nf, NfAux = 0, bdim, totDim, totDimAux = 0, e;
3839:   PetscErrorCode  ierr;

3842:   PetscFEGetQuadrature(fem, &quad);
3843:   PetscFEGetSpatialDimension(fem, &dim);
3844:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3845:   bdim = dim-1;
3846:   PetscDSGetNumFields(prob, &Nf);
3847:   PetscDSGetTotalBdDimension(prob, &totDim);
3848:   PetscDSGetComponentBdOffsets(prob, &uOff);
3849:   PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
3850:   PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3851:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3852:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3853:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3854:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3855:   PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
3856:   PetscFEGetDimension(fe, &NbI);
3857:   PetscFEGetNumComponents(fe, &NcI);
3858:   PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
3859:   PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
3860:   PetscFEGetDimension(fe, &NbJ);
3861:   PetscFEGetNumComponents(fe, &NcJ);
3862:   PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
3863:   if (probAux) {
3864:     PetscDSGetNumFields(probAux, &NfAux);
3865:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3866:     PetscDSGetComponentBdOffsets(probAux, &aOff);
3867:     PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
3868:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3869:   }
3870:   basisI    = basisField[fieldI];
3871:   basisJ    = basisField[fieldJ];
3872:   basisDerI = basisFieldDer[fieldI];
3873:   basisDerJ = basisFieldDer[fieldJ];
3874:   /* Initialize here in case the function is not defined */
3875:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3876:   PetscMemzero(g1, NcI*NcJ*bdim * sizeof(PetscScalar));
3877:   PetscMemzero(g2, NcI*NcJ*bdim * sizeof(PetscScalar));
3878:   PetscMemzero(g3, NcI*NcJ*bdim*bdim * sizeof(PetscScalar));
3879:   for (e = 0; e < Ne; ++e) {
3880:     const PetscReal *v0   = geom[e].v0;
3881:     const PetscReal *n    = geom[e].n;
3882:     const PetscReal *J    = geom[e].J;
3883:     const PetscReal *invJ = geom[e].invJ;
3884:     const PetscReal  detJ = geom[e].detJ;
3885:     const PetscReal *quadPoints, *quadWeights;
3886:     PetscInt         Nq, q;

3888:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3889:     for (q = 0; q < Nq; ++q) {
3890:       PetscInt f, g, fc, gc, c;

3892:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3893:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3894:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3895:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3896:       if (g0_func) {
3897:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3898:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, n, g0);
3899:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3900:       }
3901:       if (g1_func) {
3902:         PetscInt d, d2;
3903:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3904:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, n, refSpaceDer);
3905:         for (fc = 0; fc < NcI; ++fc) {
3906:           for (gc = 0; gc < NcJ; ++gc) {
3907:             for (d = 0; d < bdim; ++d) {
3908:               g1[(fc*NcJ+gc)*bdim+d] = 0.0;
3909:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*bdim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3910:               g1[(fc*NcJ+gc)*bdim+d] *= detJ*quadWeights[q];
3911:             }
3912:           }
3913:         }
3914:       }
3915:       if (g2_func) {
3916:         PetscInt d, d2;
3917:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3918:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, n, refSpaceDer);
3919:         for (fc = 0; fc < NcI; ++fc) {
3920:           for (gc = 0; gc < NcJ; ++gc) {
3921:             for (d = 0; d < bdim; ++d) {
3922:               g2[(fc*NcJ+gc)*bdim+d] = 0.0;
3923:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*bdim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3924:               g2[(fc*NcJ+gc)*bdim+d] *= detJ*quadWeights[q];
3925:             }
3926:           }
3927:         }
3928:       }
3929:       if (g3_func) {
3930:         PetscInt d, d2, dp, d3;
3931:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3932:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, n, g3);
3933:         for (fc = 0; fc < NcI; ++fc) {
3934:           for (gc = 0; gc < NcJ; ++gc) {
3935:             for (d = 0; d < bdim; ++d) {
3936:               for (dp = 0; dp < bdim; ++dp) {
3937:                 g3[((fc*NcJ+gc)*bdim+d)*bdim+dp] = 0.0;
3938:                 for (d2 = 0; d2 < dim; ++d2) {
3939:                   for (d3 = 0; d3 < dim; ++d3) {
3940:                     g3[((fc*NcJ+gc)*bdim+d)*bdim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3941:                   }
3942:                 }
3943:                 g3[((fc*NcJ+gc)*bdim+d)*bdim+dp] *= detJ*quadWeights[q];
3944:               }
3945:             }
3946:           }
3947:         }
3948:       }

3950:       for (f = 0; f < NbI; ++f) {
3951:         for (fc = 0; fc < NcI; ++fc) {
3952:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3953:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3954:           for (g = 0; g < NbJ; ++g) {
3955:             for (gc = 0; gc < NcJ; ++gc) {
3956:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3957:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3958:               const PetscInt fOff = eOffset+i*totDim+j;
3959:               PetscInt       d, d2;

3961:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3962:               for (d = 0; d < bdim; ++d) {
3963:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
3964:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
3965:                 for (d2 = 0; d2 < bdim; ++d2) {
3966:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
3967:                 }
3968:               }
3969:             }
3970:           }
3971:         }
3972:       }
3973:     }
3974:     if (debug > 1) {
3975:       PetscInt fc, f, gc, g;

3977:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3978:       for (fc = 0; fc < NcI; ++fc) {
3979:         for (f = 0; f < NbI; ++f) {
3980:           const PetscInt i = offsetI + f*NcI+fc;
3981:           for (gc = 0; gc < NcJ; ++gc) {
3982:             for (g = 0; g < NbJ; ++g) {
3983:               const PetscInt j = offsetJ + g*NcJ+gc;
3984:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3985:             }
3986:           }
3987:           PetscPrintf(PETSC_COMM_SELF, "\n");
3988:         }
3989:       }
3990:     }
3991:     cOffset    += totDim;
3992:     cOffsetAux += totDimAux;
3993:     eOffset    += PetscSqr(totDim);
3994:   }
3995:   return(0);
3996: }

4000: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
4001: {
4003:   fem->ops->setfromoptions          = NULL;
4004:   fem->ops->setup                   = PetscFESetUp_Basic;
4005:   fem->ops->view                    = PetscFEView_Basic;
4006:   fem->ops->destroy                 = PetscFEDestroy_Basic;
4007:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4008:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
4009:   fem->ops->integrate               = PetscFEIntegrate_Basic;
4010:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
4011:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
4012:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
4013:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
4014:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
4015:   return(0);
4016: }

4018: /*MC
4019:   PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization

4021:   Level: intermediate

4023: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4024: M*/

4028: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
4029: {
4030:   PetscFE_Basic *b;

4035:   PetscNewLog(fem,&b);
4036:   fem->data = b;

4038:   PetscFEInitialize_Basic(fem);
4039:   return(0);
4040: }

4044: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
4045: {
4046:   PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;

4050:   PetscFree(na);
4051:   return(0);
4052: }

4056: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4057:                                                   const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
4058: {
4059:   const PetscInt  debug = 0;
4060:   PetscPointFunc  f0_func;
4061:   PetscPointFunc  f1_func;
4062:   PetscQuadrature quad;
4063:   PetscReal     **basisField, **basisFieldDer;
4064:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
4065:   PetscReal      *x;
4066:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4067:   PetscInt        dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
4068:   PetscErrorCode  ierr;

4071:   PetscFEGetSpatialDimension(fem, &dim);
4072:   PetscFEGetQuadrature(fem, &quad);
4073:   PetscFEGetDimension(fem, &Nb);
4074:   PetscFEGetNumComponents(fem, &Nc);
4075:   PetscDSGetNumFields(prob, &Nf);
4076:   PetscDSGetTotalDimension(prob, &totDim);
4077:   PetscDSGetComponentOffsets(prob, &uOff);
4078:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4079:   PetscDSGetFieldOffset(prob, field, &fOffset);
4080:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4081:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
4082:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4083:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
4084:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4085:   if (probAux) {
4086:     PetscDSGetNumFields(probAux, &NfAux);
4087:     PetscDSGetTotalDimension(probAux, &totDimAux);
4088:     PetscDSGetComponentOffsets(probAux, &aOff);
4089:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4090:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4091:   }
4092:   for (e = 0; e < Ne; ++e) {
4093:     const PetscReal *quadPoints, *quadWeights;
4094:     PetscInt         Nq, q;

4096:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4097:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
4098:     PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
4099:     for (q = 0; q < Nq; ++q) {
4100:       const PetscReal *v0   = geom[e*Nq+q].v0;
4101:       const PetscReal *J    = geom[e*Nq+q].J;
4102:       const PetscReal *invJ = geom[e*Nq+q].invJ;
4103:       const PetscReal  detJ = geom[e*Nq+q].detJ;

4105:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4106:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4107:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4108:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4109:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, &f0[q*Nc]);
4110:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, refSpaceDer);
4111:       TransformF(dim, dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4112:     }
4113:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4114:     cOffset    += totDim;
4115:     cOffsetAux += totDimAux;
4116:   }
4117:   return(0);
4118: }

4122: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4123:                                                     const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
4124: {
4125:   const PetscInt   debug = 0;
4126:   PetscBdPointFunc f0_func;
4127:   PetscBdPointFunc f1_func;
4128:   PetscQuadrature  quad;
4129:   PetscReal      **basisField, **basisFieldDer;
4130:   PetscScalar     *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4131:   PetscReal       *x;
4132:   PetscInt        *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4133:   PetscInt         dim, Nf, NfAux = 0, bdim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
4134:   PetscErrorCode   ierr;

4137:   PetscFEGetSpatialDimension(fem, &dim);
4138:   dim += 1; /* Spatial dimension is one higher than topological dimension */
4139:   bdim = dim-1;
4140:   PetscFEGetQuadrature(fem, &quad);
4141:   PetscFEGetDimension(fem, &Nb);
4142:   PetscFEGetNumComponents(fem, &Nc);
4143:   PetscDSGetNumFields(prob, &Nf);
4144:   PetscDSGetTotalBdDimension(prob, &totDim);
4145:   PetscDSGetComponentBdOffsets(prob, &uOff);
4146:   PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
4147:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
4148:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
4149:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4150:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4151:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
4152:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
4153:   if (probAux) {
4154:     PetscDSGetNumFields(probAux, &NfAux);
4155:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
4156:     PetscDSGetComponentBdOffsets(probAux, &aOff);
4157:     PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
4158:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4159:   }
4160:   for (e = 0; e < Ne; ++e) {
4161:     const PetscReal *quadPoints, *quadWeights;
4162:     PetscInt         Nq, q;

4164:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4165:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
4166:     PetscMemzero(f1, Nq*Nc*bdim * sizeof(PetscScalar));
4167:     for (q = 0; q < Nq; ++q) {
4168:       const PetscReal *v0   = geom[e*Nq+q].v0;
4169:       const PetscReal *n    = geom[e*Nq+q].n;
4170:       const PetscReal *J    = geom[e*Nq+q].J;
4171:       const PetscReal *invJ = geom[e*Nq+q].invJ;
4172:       const PetscReal  detJ = geom[e*Nq+q].detJ;

4174:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4175:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
4176:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4177:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4178:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, &f0[q*Nc]);
4179:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, refSpaceDer);
4180:       TransformF(dim, bdim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4181:     }
4182:     UpdateElementVec(bdim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4183:     cOffset    += totDim;
4184:     cOffsetAux += totDimAux;
4185:   }
4186:   return(0);
4187: }

4191: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
4192:                                                   const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
4193: {
4194:   const PetscInt  debug      = 0;
4195:   PetscPointJac   g0_func;
4196:   PetscPointJac   g1_func;
4197:   PetscPointJac   g2_func;
4198:   PetscPointJac   g3_func;
4199:   PetscFE         fe;
4200:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
4201:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
4202:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
4203:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
4204:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
4205:   PetscQuadrature quad;
4206:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4207:   PetscReal      *x;
4208:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
4209:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
4210:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4211:   PetscInt        dim, Nf, NfAux = 0, totDim, totDimAux = 0, e;
4212:   PetscErrorCode  ierr;

4215:   PetscFEGetSpatialDimension(fem, &dim);
4216:   PetscFEGetQuadrature(fem, &quad);
4217:   PetscDSGetNumFields(prob, &Nf);
4218:   PetscDSGetTotalDimension(prob, &totDim);
4219:   PetscDSGetComponentOffsets(prob, &uOff);
4220:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4221:   switch(jtype) {
4222:   case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4223:   case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4224:   case PETSCFE_JACOBIAN:     PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4225:   }
4226:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4227:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4228:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
4229:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4230:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4231:   PetscFEGetDimension(fe, &NbI);
4232:   PetscFEGetNumComponents(fe, &NcI);
4233:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
4234:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
4235:   PetscFEGetDimension(fe, &NbJ);
4236:   PetscFEGetNumComponents(fe, &NcJ);
4237:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
4238:   if (probAux) {
4239:     PetscDSGetNumFields(probAux, &NfAux);
4240:     PetscDSGetTotalDimension(probAux, &totDimAux);
4241:     PetscDSGetComponentOffsets(probAux, &aOff);
4242:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4243:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4244:   }
4245:   basisI    = basisField[fieldI];
4246:   basisJ    = basisField[fieldJ];
4247:   basisDerI = basisFieldDer[fieldI];
4248:   basisDerJ = basisFieldDer[fieldJ];
4249:   /* Initialize here in case the function is not defined */
4250:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4251:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
4252:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
4253:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4254:   for (e = 0; e < Ne; ++e) {
4255:     const PetscReal *quadPoints, *quadWeights;
4256:     PetscInt         Nq, q;

4258:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4259:     for (q = 0; q < Nq; ++q) {
4260:       const PetscReal *v0   = geom[e*Nq+q].v0;
4261:       const PetscReal *J    = geom[e*Nq+q].J;
4262:       const PetscReal *invJ = geom[e*Nq+q].invJ;
4263:       const PetscReal  detJ = geom[e*Nq+q].detJ;
4264:       PetscInt         f, g, fc, gc, c;

4266:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4267:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4268:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4269:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4270:       if (g0_func) {
4271:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4272:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, g0);
4273:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
4274:       }
4275:       if (g1_func) {
4276:         PetscInt d, d2;
4277:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4278:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4279:         for (fc = 0; fc < NcI; ++fc) {
4280:           for (gc = 0; gc < NcJ; ++gc) {
4281:             for (d = 0; d < dim; ++d) {
4282:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
4283:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4284:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4285:             }
4286:           }
4287:         }
4288:       }
4289:       if (g2_func) {
4290:         PetscInt d, d2;
4291:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4292:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4293:         for (fc = 0; fc < NcI; ++fc) {
4294:           for (gc = 0; gc < NcJ; ++gc) {
4295:             for (d = 0; d < dim; ++d) {
4296:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
4297:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4298:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4299:             }
4300:           }
4301:         }
4302:       }
4303:       if (g3_func) {
4304:         PetscInt d, d2, dp, d3;
4305:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4306:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4307:         for (fc = 0; fc < NcI; ++fc) {
4308:           for (gc = 0; gc < NcJ; ++gc) {
4309:             for (d = 0; d < dim; ++d) {
4310:               for (dp = 0; dp < dim; ++dp) {
4311:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
4312:                 for (d2 = 0; d2 < dim; ++d2) {
4313:                   for (d3 = 0; d3 < dim; ++d3) {
4314:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
4315:                   }
4316:                 }
4317:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
4318:               }
4319:             }
4320:           }
4321:         }
4322:       }

4324:       for (f = 0; f < NbI; ++f) {
4325:         for (fc = 0; fc < NcI; ++fc) {
4326:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
4327:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
4328:           for (g = 0; g < NbJ; ++g) {
4329:             for (gc = 0; gc < NcJ; ++gc) {
4330:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
4331:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
4332:               const PetscInt fOff = eOffset+i*totDim+j;
4333:               PetscInt       d, d2;

4335:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
4336:               for (d = 0; d < dim; ++d) {
4337:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
4338:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
4339:                 for (d2 = 0; d2 < dim; ++d2) {
4340:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
4341:                 }
4342:               }
4343:             }
4344:           }
4345:         }
4346:       }
4347:     }
4348:     if (debug > 1) {
4349:       PetscInt fc, f, gc, g;

4351:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
4352:       for (fc = 0; fc < NcI; ++fc) {
4353:         for (f = 0; f < NbI; ++f) {
4354:           const PetscInt i = offsetI + f*NcI+fc;
4355:           for (gc = 0; gc < NcJ; ++gc) {
4356:             for (g = 0; g < NbJ; ++g) {
4357:               const PetscInt j = offsetJ + g*NcJ+gc;
4358:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
4359:             }
4360:           }
4361:           PetscPrintf(PETSC_COMM_SELF, "\n");
4362:         }
4363:       }
4364:     }
4365:     cOffset    += totDim;
4366:     cOffsetAux += totDimAux;
4367:     eOffset    += PetscSqr(totDim);
4368:   }
4369:   return(0);
4370: }

4374: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
4375: {
4377:   fem->ops->setfromoptions          = NULL;
4378:   fem->ops->setup                   = PetscFESetUp_Basic;
4379:   fem->ops->view                    = NULL;
4380:   fem->ops->destroy                 = PetscFEDestroy_Nonaffine;
4381:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4382:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
4383:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Nonaffine;
4384:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Nonaffine;
4385:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
4386:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Nonaffine;
4387:   return(0);
4388: }

4390: /*MC
4391:   PETSCFENONAFFINE = "nonaffine" - A PetscFE object that integrates with basic tiling and no vectorization for non-affine mappings

4393:   Level: intermediate

4395: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4396: M*/

4400: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
4401: {
4402:   PetscFE_Nonaffine *na;
4403:   PetscErrorCode     ierr;

4407:   PetscNewLog(fem, &na);
4408:   fem->data = na;

4410:   PetscFEInitialize_Nonaffine(fem);
4411:   return(0);
4412: }

4414: #ifdef PETSC_HAVE_OPENCL

4418: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
4419: {
4420:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4421:   PetscErrorCode  ierr;

4424:   clReleaseCommandQueue(ocl->queue_id);
4425:   ocl->queue_id = 0;
4426:   clReleaseContext(ocl->ctx_id);
4427:   ocl->ctx_id = 0;
4428:   PetscFree(ocl);
4429:   return(0);
4430: }

4432: #define STRING_ERROR_CHECK(MSG) do { string_tail += count; if (string_tail == end_of_buffer) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, MSG);} while(0)
4433: enum {LAPLACIAN = 0, ELASTICITY = 1};

4437: /* dim     Number of spatial dimensions:          2                   */
4438: /* N_b     Number of basis functions:             generated           */
4439: /* N_{bt}  Number of total basis functions:       N_b * N_{comp}      */
4440: /* N_q     Number of quadrature points:           generated           */
4441: /* N_{bs}  Number of block cells                  LCM(N_b, N_q)       */
4442: /* N_{bst} Number of block cell components        LCM(N_{bt}, N_q)    */
4443: /* N_{bl}  Number of concurrent blocks            generated           */
4444: /* N_t     Number of threads:                     N_{bl} * N_{bs}     */
4445: /* N_{cbc} Number of concurrent basis      cells: N_{bl} * N_q        */
4446: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b        */
4447: /* N_{sbc} Number of serial     basis      cells: N_{bs} / N_q        */
4448: /* N_{sqc} Number of serial     quadrature cells: N_{bs} / N_b        */
4449: /* N_{cb}  Number of serial cell batches:         input               */
4450: /* N_c     Number of total cells:                 N_{cb}*N_{t}/N_{comp} */
4451: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
4452: {
4453:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4454:   PetscQuadrature q;
4455:   char           *string_tail   = *string_buffer;
4456:   char           *end_of_buffer = *string_buffer + buffer_length;
4457:   char            float_str[]   = "float", double_str[]  = "double";
4458:   char           *numeric_str   = &(float_str[0]);
4459:   PetscInt        op            = ocl->op;
4460:   PetscBool       useField      = PETSC_FALSE;
4461:   PetscBool       useFieldDer   = PETSC_TRUE;
4462:   PetscBool       useFieldAux   = useAux;
4463:   PetscBool       useFieldDerAux= PETSC_FALSE;
4464:   PetscBool       useF0         = PETSC_TRUE;
4465:   PetscBool       useF1         = PETSC_TRUE;
4466:   PetscReal      *basis, *basisDer;
4467:   PetscInt        dim, N_b, N_c, N_q, N_t, p, d, b, c;
4468:   size_t          count;
4469:   PetscErrorCode  ierr;

4472:   PetscFEGetSpatialDimension(fem, &dim);
4473:   PetscFEGetDimension(fem, &N_b);
4474:   PetscFEGetNumComponents(fem, &N_c);
4475:   PetscFEGetQuadrature(fem, &q);
4476:   N_q  = q->numPoints;
4477:   N_t  = N_b * N_c * N_q * N_bl;
4478:   /* Enable device extension for double precision */
4479:   if (ocl->realType == PETSC_DOUBLE) {
4480:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4481: "#if defined(cl_khr_fp64)\n"
4482: "#  pragma OPENCL EXTENSION cl_khr_fp64: enable\n"
4483: "#elif defined(cl_amd_fp64)\n"
4484: "#  pragma OPENCL EXTENSION cl_amd_fp64: enable\n"
4485: "#endif\n",
4486:                               &count);STRING_ERROR_CHECK("Message to short");
4487:     numeric_str  = &(double_str[0]);
4488:   }
4489:   /* Kernel API */
4490:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4491: "\n"
4492: "__kernel void integrateElementQuadrature(int N_cb, __global %s *coefficients, __global %s *coefficientsAux, __global %s *jacobianInverses, __global %s *jacobianDeterminants, __global %s *elemVec)\n"
4493: "{\n",
4494:                        &count, numeric_str, numeric_str, numeric_str, numeric_str, numeric_str);STRING_ERROR_CHECK("Message to short");
4495:   /* Quadrature */
4496:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4497: "  /* Quadrature points\n"
4498: "   - (x1,y1,x2,y2,...) */\n"
4499: "  const %s points[%d] = {\n",
4500:                        &count, numeric_str, N_q*dim);STRING_ERROR_CHECK("Message to short");
4501:   for (p = 0; p < N_q; ++p) {
4502:     for (d = 0; d < dim; ++d) {
4503:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->points[p*dim+d]);STRING_ERROR_CHECK("Message to short");
4504:     }
4505:   }
4506:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4507:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4508: "  /* Quadrature weights\n"
4509: "   - (v1,v2,...) */\n"
4510: "  const %s weights[%d] = {\n",
4511:                        &count, numeric_str, N_q);STRING_ERROR_CHECK("Message to short");
4512:   for (p = 0; p < N_q; ++p) {
4513:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->weights[p]);STRING_ERROR_CHECK("Message to short");
4514:   }
4515:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4516:   /* Basis Functions */
4517:   PetscFEGetDefaultTabulation(fem, &basis, &basisDer, NULL);
4518:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4519: "  /* Nodal basis function evaluations\n"
4520: "    - basis component is fastest varying, the basis function, then point */\n"
4521: "  const %s Basis[%d] = {\n",
4522:                        &count, numeric_str, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4523:   for (p = 0; p < N_q; ++p) {
4524:     for (b = 0; b < N_b; ++b) {
4525:       for (c = 0; c < N_c; ++c) {
4526:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, basis[(p*N_b + b)*N_c + c]);STRING_ERROR_CHECK("Message to short");
4527:       }
4528:     }
4529:   }
4530:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4531:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4532: "\n"
4533: "  /* Nodal basis function derivative evaluations,\n"
4534: "      - derivative direction is fastest varying, then basis component, then basis function, then point */\n"
4535: "  const %s%d BasisDerivatives[%d] = {\n",
4536:                        &count, numeric_str, dim, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4537:   for (p = 0; p < N_q; ++p) {
4538:     for (b = 0; b < N_b; ++b) {
4539:       for (c = 0; c < N_c; ++c) {
4540:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "(%s%d)(", &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4541:         for (d = 0; d < dim; ++d) {
4542:           if (d > 0) {
4543:             PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, ", %g", &count, basisDer[((p*N_b + b)*dim + d)*N_c + c]);STRING_ERROR_CHECK("Message to short");
4544:           } else {
4545:             PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g", &count, basisDer[((p*N_b + b)*dim + d)*N_c + c]);STRING_ERROR_CHECK("Message to short");
4546:           }
4547:         }
4548:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "),\n", &count);STRING_ERROR_CHECK("Message to short");
4549:       }
4550:     }
4551:   }
4552:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4553:   /* Sizes */
4554:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4555: "  const int dim    = %d;                           // The spatial dimension\n"
4556: "  const int N_bl   = %d;                           // The number of concurrent blocks\n"
4557: "  const int N_b    = %d;                           // The number of basis functions\n"
4558: "  const int N_comp = %d;                           // The number of basis function components\n"
4559: "  const int N_bt   = N_b*N_comp;                    // The total number of scalar basis functions\n"
4560: "  const int N_q    = %d;                           // The number of quadrature points\n"
4561: "  const int N_bst  = N_bt*N_q;                      // The block size, LCM(N_b*N_comp, N_q), Notice that a block is not processed simultaneously\n"
4562: "  const int N_t    = N_bst*N_bl;                    // The number of threads, N_bst * N_bl\n"
4563: "  const int N_bc   = N_t/N_comp;                    // The number of cells per batch (N_b*N_q*N_bl)\n"
4564: "  const int N_sbc  = N_bst / (N_q * N_comp);\n"
4565: "  const int N_sqc  = N_bst / N_bt;\n"
4566: "  /*const int N_c    = N_cb * N_bc;*/\n"
4567: "\n"
4568: "  /* Calculated indices */\n"
4569: "  /*const int tidx    = get_local_id(0) + get_local_size(0)*get_local_id(1);*/\n"
4570: "  const int tidx    = get_local_id(0);\n"
4571: "  const int blidx   = tidx / N_bst;                  // Block number for this thread\n"
4572: "  const int bidx    = tidx %% N_bt;                   // Basis function mapped to this thread\n"
4573: "  const int cidx    = tidx %% N_comp;                 // Basis component mapped to this thread\n"
4574: "  const int qidx    = tidx %% N_q;                    // Quadrature point mapped to this thread\n"
4575: "  const int blbidx  = tidx %% N_q + blidx*N_q;        // Cell mapped to this thread in the basis phase\n"
4576: "  const int blqidx  = tidx %% N_b + blidx*N_b;        // Cell mapped to this thread in the quadrature phase\n"
4577: "  const int gidx    = get_group_id(1)*get_num_groups(0) + get_group_id(0);\n"
4578: "  const int Goffset = gidx*N_cb*N_bc;\n",
4579:                             &count, dim, N_bl, N_b, N_c, N_q);STRING_ERROR_CHECK("Message to short");
4580:   /* Local memory */
4581:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4582: "\n"
4583: "  /* Quadrature data */\n"
4584: "  %s                w;                   // $w_q$, Quadrature weight at $x_q$\n"
4585: "  __local %s         phi_i[%d];    //[N_bt*N_q];  // $\\phi_i(x_q)$, Value of the basis function $i$ at $x_q$\n"
4586: "  __local %s%d       phiDer_i[%d]; //[N_bt*N_q];  // $\\frac{\\partial\\phi_i(x_q)}{\\partial x_d}$, Value of the derivative of basis function $i$ in direction $x_d$ at $x_q$\n"
4587: "  /* Geometric data */\n"
4588: "  __local %s        detJ[%d]; //[N_t];           // $|J(x_q)|$, Jacobian determinant at $x_q$\n"
4589: "  __local %s        invJ[%d];//[N_t*dim*dim];   // $J^{-1}(x_q)$, Jacobian inverse at $x_q$\n",
4590:                             &count, numeric_str, numeric_str, N_b*N_c*N_q, numeric_str, dim, N_b*N_c*N_q, numeric_str, N_t,
4591:                             numeric_str, N_t*dim*dim, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4592:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4593: "  /* FEM data */\n"
4594: "  __local %s        u_i[%d]; //[N_t*N_bt];       // Coefficients $u_i$ of the field $u|_{\\mathcal{T}} = \\sum_i u_i \\phi_i$\n",
4595:                             &count, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4596:   if (useAux) {
4597:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4598: "  __local %s        a_i[%d]; //[N_t];            // Coefficients $a_i$ of the auxiliary field $a|_{\\mathcal{T}} = \\sum_i a_i \\phi^R_i$\n",
4599:                             &count, numeric_str, N_t);STRING_ERROR_CHECK("Message to short");
4600:   }
4601:   if (useF0) {
4602:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4603: "  /* Intermediate calculations */\n"
4604: "  __local %s         f_0[%d]; //[N_t*N_sqc];      // $f_0(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4605:                               &count, numeric_str, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4606:   }
4607:   if (useF1) {
4608:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4609: "  __local %s%d       f_1[%d]; //[N_t*N_sqc];      // $f_1(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4610:                               &count, numeric_str, dim, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4611:   }
4612:   /* TODO: If using elasticity, put in mu/lambda coefficients */
4613:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4614: "  /* Output data */\n"
4615: "  %s                e_i;                 // Coefficient $e_i$ of the residual\n\n",
4616:                             &count, numeric_str);STRING_ERROR_CHECK("Message to short");
4617:   /* One-time loads */
4618:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4619: "  /* These should be generated inline */\n"
4620: "  /* Load quadrature weights */\n"
4621: "  w = weights[qidx];\n"
4622: "  /* Load basis tabulation \\phi_i for this cell */\n"
4623: "  if (tidx < N_bt*N_q) {\n"
4624: "    phi_i[tidx]    = Basis[tidx];\n"
4625: "    phiDer_i[tidx] = BasisDerivatives[tidx];\n"
4626: "  }\n\n",
4627:                        &count);STRING_ERROR_CHECK("Message to short");
4628:   /* Batch loads */
4629:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4630: "  for (int batch = 0; batch < N_cb; ++batch) {\n"
4631: "    /* Load geometry */\n"
4632: "    detJ[tidx] = jacobianDeterminants[Goffset+batch*N_bc+tidx];\n"
4633: "    for (int n = 0; n < dim*dim; ++n) {\n"
4634: "      const int offset = n*N_t;\n"
4635: "      invJ[offset+tidx] = jacobianInverses[(Goffset+batch*N_bc)*dim*dim+offset+tidx];\n"
4636: "    }\n"
4637: "    /* Load coefficients u_i for this cell */\n"
4638: "    for (int n = 0; n < N_bt; ++n) {\n"
4639: "      const int offset = n*N_t;\n"
4640: "      u_i[offset+tidx] = coefficients[(Goffset*N_bt)+batch*N_t*N_b+offset+tidx];\n"
4641: "    }\n",
4642:                        &count);STRING_ERROR_CHECK("Message to short");
4643:   if (useAux) {
4644:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4645: "    /* Load coefficients a_i for this cell */\n"
4646: "    /* TODO: This should not be N_t here, it should be N_bc*N_comp_aux */\n"
4647: "    a_i[tidx] = coefficientsAux[Goffset+batch*N_t+tidx];\n",
4648:                             &count);STRING_ERROR_CHECK("Message to short");
4649:   }
4650:   /* Quadrature phase */
4651:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4652: "    barrier(CLK_LOCAL_MEM_FENCE);\n"
4653: "\n"
4654: "    /* Map coefficients to values at quadrature points */\n"
4655: "    for (int c = 0; c < N_sqc; ++c) {\n"
4656: "      const int cell          = c*N_bl*N_b + blqidx;\n"
4657: "      const int fidx          = (cell*N_q + qidx)*N_comp + cidx;\n",
4658:                        &count);STRING_ERROR_CHECK("Message to short");
4659:   if (useField) {
4660:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4661: "      %s  u[%d]; //[N_comp];     // $u(x_q)$, Value of the field at $x_q$\n",
4662:                               &count, numeric_str, N_c);STRING_ERROR_CHECK("Message to short");
4663:   }
4664:   if (useFieldDer) {
4665:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4666: "      %s%d   gradU[%d]; //[N_comp]; // $\\nabla u(x_q)$, Value of the field gradient at $x_q$\n",
4667:                               &count, numeric_str, dim, N_c);STRING_ERROR_CHECK("Message to short");
4668:   }
4669:   if (useFieldAux) {
4670:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4671: "      %s  a[%d]; //[1];     // $a(x_q)$, Value of the auxiliary fields at $x_q$\n",
4672:                               &count, numeric_str, 1);STRING_ERROR_CHECK("Message to short");
4673:   }
4674:   if (useFieldDerAux) {
4675:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4676: "      %s%d   gradA[%d]; //[1]; // $\\nabla a(x_q)$, Value of the auxiliary field gradient at $x_q$\n",
4677:                               &count, numeric_str, dim, 1);STRING_ERROR_CHECK("Message to short");
4678:   }
4679:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4680: "\n"
4681: "      for (int comp = 0; comp < N_comp; ++comp) {\n",
4682:                             &count);STRING_ERROR_CHECK("Message to short");
4683:   if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        u[comp] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4684:   if (useFieldDer) {
4685:     switch (dim) {
4686:     case 1:
4687:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        gradU[comp].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4688:     case 2:
4689:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        gradU[comp].x = 0.0; gradU[comp].y = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4690:     case 3:
4691:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        gradU[comp].x = 0.0; gradU[comp].y = 0.0; gradU[comp].z = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4692:     }
4693:   }
4694:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4695: "      }\n",
4696:                             &count);STRING_ERROR_CHECK("Message to short");
4697:   if (useFieldAux) {
4698:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      a[0] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");
4699:   }
4700:   if (useFieldDerAux) {
4701:     switch (dim) {
4702:     case 1:
4703:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      gradA[0].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4704:     case 2:
4705:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      gradA[0].x = 0.0; gradA[0].y = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4706:     case 3:
4707:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      gradA[0].x = 0.0; gradA[0].y = 0.0; gradA[0].z = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4708:     }
4709:   }
4710:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4711: "      /* Get field and derivatives at this quadrature point */\n"
4712: "      for (int i = 0; i < N_b; ++i) {\n"
4713: "        for (int comp = 0; comp < N_comp; ++comp) {\n"
4714: "          const int b    = i*N_comp+comp;\n"
4715: "          const int pidx = qidx*N_bt + b;\n"
4716: "          const int uidx = cell*N_bt + b;\n"
4717: "          %s%d   realSpaceDer;\n\n",
4718:                             &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4719:   if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"          u[comp] += u_i[uidx]*phi_i[pidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4720:   if (useFieldDer) {
4721:     switch (dim) {
4722:     case 2:
4723:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4724: "          realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y;\n"
4725: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4726: "          realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y;\n"
4727: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n",
4728:                            &count);STRING_ERROR_CHECK("Message to short");break;
4729:     case 3:
4730:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4731: "          realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+0]*phiDer_i[pidx].z;\n"
4732: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4733: "          realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+1]*phiDer_i[pidx].z;\n"
4734: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n"
4735: "          realSpaceDer.z = invJ[cell*dim*dim+0*dim+2]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+2]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+2]*phiDer_i[pidx].z;\n"
4736: "          gradU[comp].z += u_i[uidx]*realSpaceDer.z;\n",
4737:                            &count);STRING_ERROR_CHECK("Message to short");break;
4738:     }
4739:   }
4740:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4741: "        }\n"
4742: "      }\n",
4743:                             &count);STRING_ERROR_CHECK("Message to short");
4744:   if (useFieldAux) {
4745:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"          a[0] += a_i[cell];\n", &count);STRING_ERROR_CHECK("Message to short");
4746:   }
4747:   /* Calculate residual at quadrature points: Should be generated by an weak form egine */
4748:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4749: "      /* Process values at quadrature points */\n",
4750:                             &count);STRING_ERROR_CHECK("Message to short");
4751:   switch (op) {
4752:   case LAPLACIAN:
4753:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4754:     if (useF1) {
4755:       if (useAux) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_1[fidx] = a[0]*gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4756:       else        {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_1[fidx] = gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4757:     }
4758:     break;
4759:   case ELASTICITY:
4760:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4761:     if (useF1) {
4762:     switch (dim) {
4763:     case 2:
4764:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4765: "      switch (cidx) {\n"
4766: "      case 0:\n"
4767: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].x + gradU[0].x);\n"
4768: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].y + gradU[1].x);\n"
4769: "        break;\n"
4770: "      case 1:\n"
4771: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].x + gradU[0].y);\n"
4772: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].y + gradU[1].y);\n"
4773: "      }\n",
4774:                            &count);STRING_ERROR_CHECK("Message to short");break;
4775:     case 3:
4776:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4777: "      switch (cidx) {\n"
4778: "      case 0:\n"
4779: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].x + gradU[0].x);\n"
4780: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].y + gradU[1].x);\n"
4781: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].z + gradU[2].x);\n"
4782: "        break;\n"
4783: "      case 1:\n"
4784: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].x + gradU[0].y);\n"
4785: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[1].y);\n"
4786: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[2].y);\n"
4787: "        break;\n"
4788: "      case 2:\n"
4789: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].x + gradU[0].z);\n"
4790: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[1].z);\n"
4791: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[2].z);\n"
4792: "      }\n",
4793:                            &count);STRING_ERROR_CHECK("Message to short");break;
4794:     }}
4795:     break;
4796:   default:
4797:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PDE operator %d is not supported", op);
4798:   }
4799:   if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_0[fidx] *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");}
4800:   if (useF1) {
4801:     switch (dim) {
4802:     case 1:
4803:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_1[fidx].x *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4804:     case 2:
4805:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_1[fidx].x *= detJ[cell]*w; f_1[fidx].y *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4806:     case 3:
4807:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_1[fidx].x *= detJ[cell]*w; f_1[fidx].y *= detJ[cell]*w; f_1[fidx].z *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4808:     }
4809:   }
4810:   /* Thread transpose */
4811:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4812: "    }\n\n"
4813: "    /* ==== TRANSPOSE THREADS ==== */\n"
4814: "    barrier(CLK_LOCAL_MEM_FENCE);\n\n",
4815:                        &count);STRING_ERROR_CHECK("Message to short");
4816:   /* Basis phase */
4817:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4818: "    /* Map values at quadrature points to coefficients */\n"
4819: "    for (int c = 0; c < N_sbc; ++c) {\n"
4820: "      const int cell = c*N_bl*N_q + blbidx; /* Cell number in batch */\n"
4821: "\n"
4822: "      e_i = 0.0;\n"
4823: "      for (int q = 0; q < N_q; ++q) {\n"
4824: "        const int pidx = q*N_bt + bidx;\n"
4825: "        const int fidx = (cell*N_q + q)*N_comp + cidx;\n"
4826: "        %s%d   realSpaceDer;\n\n",
4827:                        &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");

4829:   if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"        e_i += phi_i[pidx]*f_0[fidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4830:   if (useF1) {
4831:     switch (dim) {
4832:     case 2:
4833:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4834: "        realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y;\n"
4835: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4836: "        realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y;\n"
4837: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n",
4838:                            &count);STRING_ERROR_CHECK("Message to short");break;
4839:     case 3:
4840:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4841: "        realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+0]*phiDer_i[pidx].z;\n"
4842: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4843: "        realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+1]*phiDer_i[pidx].z;\n"
4844: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n"
4845: "        realSpaceDer.z = invJ[cell*dim*dim+0*dim+2]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+2]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+2]*phiDer_i[pidx].z;\n"
4846: "        e_i           += realSpaceDer.z*f_1[fidx].z;\n",
4847:                            &count);STRING_ERROR_CHECK("Message to short");break;
4848:     }
4849:   }
4850:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4851: "      }\n"
4852: "      /* Write element vector for N_{cbc} cells at a time */\n"
4853: "      elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
4854: "    }\n"
4855: "    /* ==== Could do one write per batch ==== */\n"
4856: "  }\n"
4857: "  return;\n"
4858: "}\n",
4859:                        &count);STRING_ERROR_CHECK("Message to short");
4860:   return(0);
4861: }

4865: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
4866: {
4867:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4868:   PetscInt        dim, N_bl;
4869:   PetscBool       flg;
4870:   char           *buffer;
4871:   size_t          len;
4872:   char            errMsg[8192];
4873:   cl_int          ierr2;
4874:   PetscErrorCode  ierr;

4877:   PetscFEGetSpatialDimension(fem, &dim);
4878:   PetscMalloc1(8192, &buffer);
4879:   PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
4880:   PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
4881:   PetscOptionsHasName(((PetscObject)fem)->options,((PetscObject)fem)->prefix, "-petscfe_opencl_kernel_print", &flg);
4882:   if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
4883:   len  = strlen(buffer);
4884:   *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
4885:   clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
4886:   if (ierr != CL_SUCCESS) {
4887:     clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
4888:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
4889:   }
4890:   PetscFree(buffer);
4891:   *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
4892:   return(0);
4893: }

4897: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
4898: {
4899:   const PetscInt Nblocks = N/blockSize;

4902:   if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
4903:   *z = 1;
4904:   for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
4905:     *y = Nblocks / *x;
4906:     if (*x * *y == Nblocks) break;
4907:   }
4908:   if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
4909:   return(0);
4910: }

4914: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
4915: {
4916:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4917:   PetscStageLog     stageLog;
4918:   PetscEventPerfLog eventLog = NULL;
4919:   PetscInt          stage;
4920:   PetscErrorCode    ierr;

4923:   PetscLogGetStageLog(&stageLog);
4924:   PetscStageLogGetCurrent(stageLog, &stage);
4925:   PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
4926:     /* Log performance info */
4927:   eventLog->eventInfo[ocl->residualEvent].count++;
4928:   eventLog->eventInfo[ocl->residualEvent].time  += time;
4929:   eventLog->eventInfo[ocl->residualEvent].flops += flops;
4930:   return(0);
4931: }

4935: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4936:                                                const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
4937: {
4938:   /* Nbc = batchSize */
4939:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4940:   PetscPointFunc    f0_func;
4941:   PetscPointFunc    f1_func;
4942:   PetscQuadrature   q;
4943:   PetscInt          dim;
4944:   PetscInt          N_b;    /* The number of basis functions */
4945:   PetscInt          N_comp; /* The number of basis function components */
4946:   PetscInt          N_bt;   /* The total number of scalar basis functions */
4947:   PetscInt          N_q;    /* The number of quadrature points */
4948:   PetscInt          N_bst;  /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
4949:   PetscInt          N_t;    /* The number of threads, N_bst * N_bl */
4950:   PetscInt          N_bl;   /* The number of blocks */
4951:   PetscInt          N_bc;   /* The batch size, N_bl*N_q*N_b */
4952:   PetscInt          N_cb;   /* The number of batches */
4953:   PetscInt          numFlops, f0Flops = 0, f1Flops = 0;
4954:   PetscBool         useAux      = probAux ? PETSC_TRUE : PETSC_FALSE;
4955:   PetscBool         useField    = PETSC_FALSE;
4956:   PetscBool         useFieldDer = PETSC_TRUE;
4957:   PetscBool         useF0       = PETSC_TRUE;
4958:   PetscBool         useF1       = PETSC_TRUE;
4959:   /* OpenCL variables */
4960:   cl_program        ocl_prog;
4961:   cl_kernel         ocl_kernel;
4962:   cl_event          ocl_ev;         /* The event for tracking kernel execution */
4963:   cl_ulong          ns_start;       /* Nanoseconds counter on GPU at kernel start */
4964:   cl_ulong          ns_end;         /* Nanoseconds counter on GPU at kernel stop */
4965:   cl_mem            o_jacobianInverses, o_jacobianDeterminants;
4966:   cl_mem            o_coefficients, o_coefficientsAux, o_elemVec;
4967:   float            *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
4968:   double           *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
4969:   PetscReal        *r_invJ = NULL, *r_detJ = NULL;
4970:   void             *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
4971:   size_t            local_work_size[3], global_work_size[3];
4972:   size_t            realSize, x, y, z;
4973:   PetscErrorCode    ierr;

4976:   if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
4977:   PetscFEGetSpatialDimension(fem, &dim);
4978:   PetscFEGetQuadrature(fem, &q);
4979:   PetscFEGetDimension(fem, &N_b);
4980:   PetscFEGetNumComponents(fem, &N_comp);
4981:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4982:   PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
4983:   N_bt  = N_b*N_comp;
4984:   N_q   = q->numPoints;
4985:   N_bst = N_bt*N_q;
4986:   N_t   = N_bst*N_bl;
4987:   if (N_bc*N_comp != N_t) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of threads %d should be %d * %d", N_t, N_bc, N_comp);
4988:   /* Calculate layout */
4989:   if (Ne % (N_cb*N_bc)) { /* Remainder cells */
4990:     PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
4991:     return(0);
4992:   }
4993:   PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
4994:   local_work_size[0]  = N_bc*N_comp;
4995:   local_work_size[1]  = 1;
4996:   local_work_size[2]  = 1;
4997:   global_work_size[0] = x * local_work_size[0];
4998:   global_work_size[1] = y * local_work_size[1];
4999:   global_work_size[2] = z * local_work_size[2];
5000:   PetscInfo7(fem, "GPU layout grid(%d,%d,%d) block(%d,%d,%d) with %d batches\n", x, y, z, local_work_size[0], local_work_size[1], local_work_size[2], N_cb);
5001:   PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
5002:   /* Generate code */
5003:   if (probAux) {
5004:     PetscSpace P;
5005:     PetscInt   NfAux, order, f;

5007:     PetscDSGetNumFields(probAux, &NfAux);
5008:     for (f = 0; f < NfAux; ++f) {
5009:       PetscFE feAux;

5011:       PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
5012:       PetscFEGetBasisSpace(feAux, &P);
5013:       PetscSpaceGetOrder(P, &order);
5014:       if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
5015:     }
5016:   }
5017:   PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
5018:   /* Create buffers on the device and send data over */
5019:   PetscDataTypeGetSize(ocl->realType, &realSize);
5020:   if (sizeof(PetscReal) != realSize) {
5021:     switch (ocl->realType) {
5022:     case PETSC_FLOAT:
5023:     {
5024:       PetscInt c, b, d;

5026:       PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
5027:       for (c = 0; c < Ne; ++c) {
5028:         f_detJ[c] = (float) geom[c].detJ;
5029:         for (d = 0; d < dim*dim; ++d) {
5030:           f_invJ[c*dim*dim+d] = (float) geom[c].invJ[d];
5031:         }
5032:         for (b = 0; b < N_bt; ++b) {
5033:           f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
5034:         }
5035:       }
5036:       if (coefficientsAux) { /* Assume P0 */
5037:         for (c = 0; c < Ne; ++c) {
5038:           f_coeffAux[c] = (float) coefficientsAux[c];
5039:         }
5040:       }
5041:       oclCoeff      = (void *) f_coeff;
5042:       if (coefficientsAux) {
5043:         oclCoeffAux = (void *) f_coeffAux;
5044:       } else {
5045:         oclCoeffAux = NULL;
5046:       }
5047:       oclInvJ       = (void *) f_invJ;
5048:       oclDetJ       = (void *) f_detJ;
5049:     }
5050:     break;
5051:     case PETSC_DOUBLE:
5052:     {
5053:       PetscInt c, b, d;

5055:       PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
5056:       for (c = 0; c < Ne; ++c) {
5057:         d_detJ[c] = (double) geom[c].detJ;
5058:         for (d = 0; d < dim*dim; ++d) {
5059:           d_invJ[c*dim*dim+d] = (double) geom[c].invJ[d];
5060:         }
5061:         for (b = 0; b < N_bt; ++b) {
5062:           d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
5063:         }
5064:       }
5065:       if (coefficientsAux) { /* Assume P0 */
5066:         for (c = 0; c < Ne; ++c) {
5067:           d_coeffAux[c] = (double) coefficientsAux[c];
5068:         }
5069:       }
5070:       oclCoeff      = (void *) d_coeff;
5071:       if (coefficientsAux) {
5072:         oclCoeffAux = (void *) d_coeffAux;
5073:       } else {
5074:         oclCoeffAux = NULL;
5075:       }
5076:       oclInvJ       = (void *) d_invJ;
5077:       oclDetJ       = (void *) d_detJ;
5078:     }
5079:     break;
5080:     default:
5081:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
5082:     }
5083:   } else {
5084:     PetscInt c, d;

5086:     PetscMalloc2(Ne*dim*dim,&r_invJ,Ne,&r_detJ);
5087:     for (c = 0; c < Ne; ++c) {
5088:       r_detJ[c] = geom[c].detJ;
5089:       for (d = 0; d < dim*dim; ++d) {
5090:         r_invJ[c*dim*dim+d] = geom[c].invJ[d];
5091:       }
5092:     }
5093:     oclCoeff    = (void *) coefficients;
5094:     oclCoeffAux = (void *) coefficientsAux;
5095:     oclInvJ     = (void *) r_invJ;
5096:     oclDetJ     = (void *) r_detJ;
5097:   }
5098:   o_coefficients         = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt    * realSize, oclCoeff,    &ierr);
5099:   if (coefficientsAux) {
5100:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclCoeffAux, &ierr);
5101:   } else {
5102:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY,                        Ne         * realSize, oclCoeffAux, &ierr);
5103:   }
5104:   o_jacobianInverses     = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ,     &ierr);
5105:   o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclDetJ,     &ierr);
5106:   o_elemVec              = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY,                       Ne*N_bt    * realSize, NULL,        &ierr);
5107:   /* Kernel launch */
5108:   clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
5109:   clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
5110:   clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
5111:   clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
5112:   clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
5113:   clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
5114:   clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
5115:   /* Read data back from device */
5116:   if (sizeof(PetscReal) != realSize) {
5117:     switch (ocl->realType) {
5118:     case PETSC_FLOAT:
5119:     {
5120:       float   *elem;
5121:       PetscInt c, b;

5123:       PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
5124:       PetscMalloc1(Ne*N_bt, &elem);
5125:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
5126:       for (c = 0; c < Ne; ++c) {
5127:         for (b = 0; b < N_bt; ++b) {
5128:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
5129:         }
5130:       }
5131:       PetscFree(elem);
5132:     }
5133:     break;
5134:     case PETSC_DOUBLE:
5135:     {
5136:       double  *elem;
5137:       PetscInt c, b;

5139:       PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
5140:       PetscMalloc1(Ne*N_bt, &elem);
5141:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
5142:       for (c = 0; c < Ne; ++c) {
5143:         for (b = 0; b < N_bt; ++b) {
5144:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
5145:         }
5146:       }
5147:       PetscFree(elem);
5148:     }
5149:     break;
5150:     default:
5151:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
5152:     }
5153:   } else {
5154:     PetscFree2(r_invJ,r_detJ);
5155:     clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
5156:   }
5157:   /* Log performance */
5158:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
5159:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END,   sizeof(cl_ulong), &ns_end,   NULL);
5160:   f0Flops = 0;
5161:   switch (ocl->op) {
5162:   case LAPLACIAN:
5163:     f1Flops = useAux ? dim : 0;break;
5164:   case ELASTICITY:
5165:     f1Flops = 2*dim*dim;break;
5166:   }
5167:   numFlops = Ne*(
5168:     N_q*(
5169:       N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
5170:       /*+
5171:        N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
5172:       +
5173:       N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
5174:     +
5175:     N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
5176:   PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
5177:   /* Cleanup */
5178:   clReleaseMemObject(o_coefficients);
5179:   clReleaseMemObject(o_coefficientsAux);
5180:   clReleaseMemObject(o_jacobianInverses);
5181:   clReleaseMemObject(o_jacobianDeterminants);
5182:   clReleaseMemObject(o_elemVec);
5183:   clReleaseKernel(ocl_kernel);
5184:   clReleaseProgram(ocl_prog);
5185:   return(0);
5186: }

5190: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
5191: {
5193:   fem->ops->setfromoptions          = NULL;
5194:   fem->ops->setup                   = PetscFESetUp_Basic;
5195:   fem->ops->view                    = NULL;
5196:   fem->ops->destroy                 = PetscFEDestroy_OpenCL;
5197:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
5198:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
5199:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_OpenCL;
5200:   fem->ops->integratebdresidual     = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
5201:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
5202:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
5203:   return(0);
5204: }

5206: /*MC
5207:   PETSCFEOPENCL = "opencl" - A PetscFE object that integrates using a vectorized OpenCL implementation

5209:   Level: intermediate

5211: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5212: M*/

5216: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
5217: {
5218:   PetscFE_OpenCL *ocl;
5219:   cl_uint         num_platforms;
5220:   cl_platform_id  platform_ids[42];
5221:   cl_uint         num_devices;
5222:   cl_device_id    device_ids[42];
5223:   cl_int          ierr2;
5224:   PetscErrorCode  ierr;

5228:   PetscNewLog(fem,&ocl);
5229:   fem->data = ocl;

5231:   /* Init Platform */
5232:   clGetPlatformIDs(42, platform_ids, &num_platforms);
5233:   if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
5234:   ocl->pf_id = platform_ids[0];
5235:   /* Init Device */
5236:   clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
5237:   if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
5238:   ocl->dev_id = device_ids[0];
5239:   /* Create context with one command queue */
5240:   ocl->ctx_id   = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
5241:   ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
5242:   /* Types */
5243:   ocl->realType = PETSC_FLOAT;
5244:   /* Register events */
5245:   PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
5246:   /* Equation handling */
5247:   ocl->op = LAPLACIAN;

5249:   PetscFEInitialize_OpenCL(fem);
5250:   return(0);
5251: }

5255: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
5256: {
5257:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

5261:   ocl->realType = realType;
5262:   return(0);
5263: }

5267: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
5268: {
5269:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

5274:   *realType = ocl->realType;
5275:   return(0);
5276: }

5278: #endif /* PETSC_HAVE_OPENCL */

5282: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
5283: {
5284:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5285:   PetscErrorCode     ierr;

5288:   CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5289:   PetscFree(cmp->embedding);
5290:   PetscFree(cmp);
5291:   return(0);
5292: }

5296: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
5297: {
5298:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5299:   DM                 K;
5300:   PetscReal         *subpoint;
5301:   PetscBLASInt      *pivots;
5302:   PetscBLASInt       n, info;
5303:   PetscScalar       *work, *invVscalar;
5304:   PetscInt           dim, pdim, spdim, j, s;
5305:   PetscErrorCode     ierr;

5308:   /* Get affine mapping from reference cell to each subcell */
5309:   PetscDualSpaceGetDM(fem->dualSpace, &K);
5310:   DMGetDimension(K, &dim);
5311:   DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
5312:   CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5313:   /* Determine dof embedding into subelements */
5314:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5315:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
5316:   PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
5317:   DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
5318:   for (s = 0; s < cmp->numSubelements; ++s) {
5319:     PetscInt sd = 0;

5321:     for (j = 0; j < pdim; ++j) {
5322:       PetscBool       inside;
5323:       PetscQuadrature f;
5324:       PetscInt        d, e;

5326:       PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
5327:       /* Apply transform to first point, and check that point is inside subcell */
5328:       for (d = 0; d < dim; ++d) {
5329:         subpoint[d] = -1.0;
5330:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
5331:       }
5332:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5333:       if (inside) {cmp->embedding[s*spdim+sd++] = j;}
5334:     }
5335:     if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
5336:   }
5337:   DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
5338:   /* Construct the change of basis from prime basis to nodal basis for each subelement */
5339:   PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
5340:   PetscMalloc2(spdim,&pivots,spdim,&work);
5341: #if defined(PETSC_USE_COMPLEX)
5342:   PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);
5343: #else
5344:   invVscalar = fem->invV;
5345: #endif
5346:   for (s = 0; s < cmp->numSubelements; ++s) {
5347:     for (j = 0; j < spdim; ++j) {
5348:       PetscReal      *Bf;
5349:       PetscQuadrature f;
5350:       PetscInt        q, k;

5352:       PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
5353:       PetscMalloc1(f->numPoints*spdim,&Bf);
5354:       PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
5355:       for (k = 0; k < spdim; ++k) {
5356:         /* n_j \cdot \phi_k */
5357:         invVscalar[(s*spdim + j)*spdim+k] = 0.0;
5358:         for (q = 0; q < f->numPoints; ++q) {
5359:           invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
5360:         }
5361:       }
5362:       PetscFree(Bf);
5363:     }
5364:     n = spdim;
5365:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info));
5366:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info));
5367:   }
5368: #if defined(PETSC_USE_COMPLEX)
5369:   for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
5370:   PetscFree(invVscalar);
5371: #endif
5372:   PetscFree2(pivots,work);
5373:   return(0);
5374: }

5378: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
5379: {
5380:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5381:   DM                 dm;
5382:   PetscInt           pdim;  /* Dimension of FE space P */
5383:   PetscInt           spdim; /* Dimension of subelement FE space P */
5384:   PetscInt           dim;   /* Spatial dimension */
5385:   PetscInt           comp;  /* Field components */
5386:   PetscInt          *subpoints;
5387:   PetscReal         *tmpB, *tmpD, *tmpH, *subpoint;
5388:   PetscInt           p, s, d, e, j, k;
5389:   PetscErrorCode     ierr;

5392:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
5393:   DMGetDimension(dm, &dim);
5394:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
5395:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5396:   PetscFEGetNumComponents(fem, &comp);
5397:   /* Divide points into subelements */
5398:   DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
5399:   DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
5400:   for (p = 0; p < npoints; ++p) {
5401:     for (s = 0; s < cmp->numSubelements; ++s) {
5402:       PetscBool inside;

5404:       /* Apply transform, and check that point is inside cell */
5405:       for (d = 0; d < dim; ++d) {
5406:         subpoint[d] = -1.0;
5407:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
5408:       }
5409:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5410:       if (inside) {subpoints[p] = s; break;}
5411:     }
5412:     if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
5413:   }
5414:   DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
5415:   /* Evaluate the prime basis functions at all points */
5416:   if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5417:   if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5418:   if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5419:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
5420:   /* Translate to the nodal basis */
5421:   if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
5422:   if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
5423:   if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
5424:   for (p = 0; p < npoints; ++p) {
5425:     const PetscInt s = subpoints[p];

5427:     if (B) {
5428:       /* Multiply by V^{-1} (spdim x spdim) */
5429:       for (j = 0; j < spdim; ++j) {
5430:         const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
5431:         PetscInt       c;

5433:         B[i] = 0.0;
5434:         for (k = 0; k < spdim; ++k) {
5435:           B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
5436:         }
5437:         for (c = 1; c < comp; ++c) {
5438:           B[i+c] = B[i];
5439:         }
5440:       }
5441:     }
5442:     if (D) {
5443:       /* Multiply by V^{-1} (spdim x spdim) */
5444:       for (j = 0; j < spdim; ++j) {
5445:         for (d = 0; d < dim; ++d) {
5446:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
5447:           PetscInt       c;

5449:           D[i] = 0.0;
5450:           for (k = 0; k < spdim; ++k) {
5451:             D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
5452:           }
5453:           for (c = 1; c < comp; ++c) {
5454:             D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
5455:           }
5456:         }
5457:       }
5458:     }
5459:     if (H) {
5460:       /* Multiply by V^{-1} (pdim x pdim) */
5461:       for (j = 0; j < spdim; ++j) {
5462:         for (d = 0; d < dim*dim; ++d) {
5463:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
5464:           PetscInt       c;

5466:           H[i] = 0.0;
5467:           for (k = 0; k < spdim; ++k) {
5468:             H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
5469:           }
5470:           for (c = 1; c < comp; ++c) {
5471:             H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
5472:           }
5473:         }
5474:       }
5475:     }
5476:   }
5477:   DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
5478:   if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5479:   if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5480:   if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5481:   return(0);
5482: }

5486: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
5487: {
5489:   fem->ops->setfromoptions          = NULL;
5490:   fem->ops->setup                   = PetscFESetUp_Composite;
5491:   fem->ops->view                    = NULL;
5492:   fem->ops->destroy                 = PetscFEDestroy_Composite;
5493:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
5494:   fem->ops->gettabulation           = PetscFEGetTabulation_Composite;
5495:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
5496:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
5497:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
5498:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
5499:   return(0);
5500: }

5502: /*MC
5503:   PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element

5505:   Level: intermediate

5507: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5508: M*/

5512: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
5513: {
5514:   PetscFE_Composite *cmp;
5515:   PetscErrorCode     ierr;

5519:   PetscNewLog(fem, &cmp);
5520:   fem->data = cmp;

5522:   cmp->cellRefiner    = REFINER_NOOP;
5523:   cmp->numSubelements = -1;
5524:   cmp->v0             = NULL;
5525:   cmp->jac            = NULL;

5527:   PetscFEInitialize_Composite(fem);
5528:   return(0);
5529: }

5533: /*@C
5534:   PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement

5536:   Not collective

5538:   Input Parameter:
5539: . fem - The PetscFE object

5541:   Output Parameters:
5542: + blockSize - The number of elements in a block
5543: . numBlocks - The number of blocks in a batch
5544: . batchSize - The number of elements in a batch
5545: - numBatches - The number of batches in a chunk

5547:   Level: intermediate

5549: .seealso: PetscFECreate()
5550: @*/
5551: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
5552: {
5553:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;

5561:   return(0);
5562: }

5566: /*@
5567:   PetscFEGetDimension - Get the dimension of the finite element space on a cell

5569:   Not collective

5571:   Input Parameter:
5572: . fe - The PetscFE

5574:   Output Parameter:
5575: . dim - The dimension

5577:   Level: intermediate

5579: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
5580: @*/
5581: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
5582: {

5588:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
5589:   return(0);
5590: }

5592: /*
5593: Purpose: Compute element vector for chunk of elements

5595: Input:
5596:   Sizes:
5597:      Ne:  number of elements
5598:      Nf:  number of fields
5599:      PetscFE
5600:        dim: spatial dimension
5601:        Nb:  number of basis functions
5602:        Nc:  number of field components
5603:        PetscQuadrature
5604:          Nq:  number of quadrature points

5606:   Geometry:
5607:      PetscFECellGeom[Ne] possibly *Nq
5608:        PetscReal v0s[dim]
5609:        PetscReal n[dim]
5610:        PetscReal jacobians[dim*dim]
5611:        PetscReal jacobianInverses[dim*dim]
5612:        PetscReal jacobianDeterminants
5613:   FEM:
5614:      PetscFE
5615:        PetscQuadrature
5616:          PetscReal   quadPoints[Nq*dim]
5617:          PetscReal   quadWeights[Nq]
5618:        PetscReal   basis[Nq*Nb*Nc]
5619:        PetscReal   basisDer[Nq*Nb*Nc*dim]
5620:      PetscScalar coefficients[Ne*Nb*Nc]
5621:      PetscScalar elemVec[Ne*Nb*Nc]

5623:   Problem:
5624:      PetscInt f: the active field
5625:      f0, f1

5627:   Work Space:
5628:      PetscFE
5629:        PetscScalar f0[Nq*dim];
5630:        PetscScalar f1[Nq*dim*dim];
5631:        PetscScalar u[Nc];
5632:        PetscScalar gradU[Nc*dim];
5633:        PetscReal   x[dim];
5634:        PetscScalar realSpaceDer[dim];

5636: Purpose: Compute element vector for N_cb batches of elements

5638: Input:
5639:   Sizes:
5640:      N_cb: Number of serial cell batches

5642:   Geometry:
5643:      PetscReal v0s[Ne*dim]
5644:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
5645:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5646:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
5647:   FEM:
5648:      static PetscReal   quadPoints[Nq*dim]
5649:      static PetscReal   quadWeights[Nq]
5650:      static PetscReal   basis[Nq*Nb*Nc]
5651:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
5652:      PetscScalar coefficients[Ne*Nb*Nc]
5653:      PetscScalar elemVec[Ne*Nb*Nc]

5655: ex62.c:
5656:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
5657:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
5658:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
5659:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

5661: ex52.c:
5662:   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)
5663:   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)

5665: ex52_integrateElement.cu
5666: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)

5668: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
5669:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5670:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

5672: ex52_integrateElementOpenCL.c:
5673: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
5674:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5675:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

5677: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
5678: */

5682: /*@C
5683:   PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration

5685:   Not collective

5687:   Input Parameters:
5688: + fem          - The PetscFE object for the field being integrated
5689: . prob         - The PetscDS specifying the discretizations and continuum functions
5690: . field        - The field being integrated
5691: . Ne           - The number of elements in the chunk
5692: . geom         - The cell geometry for each cell in the chunk
5693: . coefficients - The array of FEM basis coefficients for the elements
5694: . probAux      - The PetscDS specifying the auxiliary discretizations
5695: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5697:   Output Parameter
5698: . integral     - the integral for this field

5700:   Level: developer

5702: .seealso: PetscFEIntegrateResidual()
5703: @*/
5704: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5705:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
5706: {

5712:   if (fem->ops->integrate) {(*fem->ops->integrate)(fem, prob, field, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
5713:   return(0);
5714: }

5718: /*@C
5719:   PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration

5721:   Not collective

5723:   Input Parameters:
5724: + fem          - The PetscFE object for the field being integrated
5725: . prob         - The PetscDS specifying the discretizations and continuum functions
5726: . field        - The field being integrated
5727: . Ne           - The number of elements in the chunk
5728: . geom         - The cell geometry for each cell in the chunk
5729: . coefficients - The array of FEM basis coefficients for the elements
5730: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5731: . probAux      - The PetscDS specifying the auxiliary discretizations
5732: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5733: - t            - The time

5735:   Output Parameter
5736: . elemVec      - the element residual vectors from each element

5738:   Note:
5739: $ Loop over batch of elements (e):
5740: $   Loop over quadrature points (q):
5741: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
5742: $     Call f_0 and f_1
5743: $   Loop over element vector entries (f,fc --> i):
5744: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

5746:   Level: developer

5748: .seealso: PetscFEIntegrateResidual()
5749: @*/
5750: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5751:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
5752: {

5758:   if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
5759:   return(0);
5760: }

5764: /*@C
5765:   PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary

5767:   Not collective

5769:   Input Parameters:
5770: + fem          - The PetscFE object for the field being integrated
5771: . prob         - The PetscDS specifying the discretizations and continuum functions
5772: . field        - The field being integrated
5773: . Ne           - The number of elements in the chunk
5774: . geom         - The cell geometry for each cell in the chunk
5775: . coefficients - The array of FEM basis coefficients for the elements
5776: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5777: . probAux      - The PetscDS specifying the auxiliary discretizations
5778: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5779: - t            - The time

5781:   Output Parameter
5782: . elemVec      - the element residual vectors from each element

5784:   Level: developer

5786: .seealso: PetscFEIntegrateResidual()
5787: @*/
5788: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5789:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
5790: {

5795:   if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);}
5796:   return(0);
5797: }

5801: /*@C
5802:   PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration

5804:   Not collective

5806:   Input Parameters:
5807: + fem          - The PetscFE object for the field being integrated
5808: . prob         - The PetscDS specifying the discretizations and continuum functions
5809: . jtype        - The type of matrix pointwise functions that should be used
5810: . fieldI       - The test field being integrated
5811: . fieldJ       - The basis field being integrated
5812: . Ne           - The number of elements in the chunk
5813: . geom         - The cell geometry for each cell in the chunk
5814: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5815: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5816: . probAux      - The PetscDS specifying the auxiliary discretizations
5817: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5818: . t            - The time
5819: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

5821:   Output Parameter
5822: . elemMat      - the element matrices for the Jacobian from each element

5824:   Note:
5825: $ Loop over batch of elements (e):
5826: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5827: $     Loop over quadrature points (q):
5828: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5829: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5830: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5831: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5832: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5833: */
5834: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5835:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
5836: {

5841:   if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, jtype, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
5842:   return(0);
5843: }

5847: /*C
5848:   PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration

5850:   Not collective

5852:   Input Parameters:
5853: + fem          = The PetscFE object for the field being integrated
5854: . prob         - The PetscDS specifying the discretizations and continuum functions
5855: . fieldI       - The test field being integrated
5856: . fieldJ       - The basis field being integrated
5857: . Ne           - The number of elements in the chunk
5858: . geom         - The cell geometry for each cell in the chunk
5859: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5860: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5861: . probAux      - The PetscDS specifying the auxiliary discretizations
5862: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5863: . t            - The time
5864: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

5866:   Output Parameter
5867: . elemMat              - the element matrices for the Jacobian from each element

5869:   Note:
5870: $ Loop over batch of elements (e):
5871: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5872: $     Loop over quadrature points (q):
5873: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5874: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5875: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5876: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5877: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5878: */
5879: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5880:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
5881: {

5886:   if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);}
5887:   return(0);
5888: }

5892: /*@
5893:   PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
5894:   to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
5895:   sparsity). It is also used to create an interpolation between regularly refined meshes.

5897:   Input Parameter:
5898: . fe - The initial PetscFE

5900:   Output Parameter:
5901: . feRef - The refined PetscFE

5903:   Level: developer

5905: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5906: @*/
5907: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
5908: {
5909:   PetscSpace       P, Pref;
5910:   PetscDualSpace   Q, Qref;
5911:   DM               K, Kref;
5912:   PetscQuadrature  q, qref;
5913:   const PetscReal *v0, *jac;
5914:   PetscInt         numComp, numSubelements;
5915:   PetscErrorCode   ierr;

5918:   PetscFEGetBasisSpace(fe, &P);
5919:   PetscFEGetDualSpace(fe, &Q);
5920:   PetscFEGetQuadrature(fe, &q);
5921:   PetscDualSpaceGetDM(Q, &K);
5922:   /* Create space */
5923:   PetscObjectReference((PetscObject) P);
5924:   Pref = P;
5925:   /* Create dual space */
5926:   PetscDualSpaceDuplicate(Q, &Qref);
5927:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
5928:   PetscDualSpaceSetDM(Qref, Kref);
5929:   DMDestroy(&Kref);
5930:   PetscDualSpaceSetUp(Qref);
5931:   /* Create element */
5932:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
5933:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
5934:   PetscFESetBasisSpace(*feRef, Pref);
5935:   PetscFESetDualSpace(*feRef, Qref);
5936:   PetscFEGetNumComponents(fe,    &numComp);
5937:   PetscFESetNumComponents(*feRef, numComp);
5938:   PetscFESetUp(*feRef);
5939:   PetscSpaceDestroy(&Pref);
5940:   PetscDualSpaceDestroy(&Qref);
5941:   /* Create quadrature */
5942:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
5943:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
5944:   PetscFESetQuadrature(*feRef, qref);
5945:   PetscQuadratureDestroy(&qref);
5946:   return(0);
5947: }

5951: /*@
5952:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

5954:   Collective on DM

5956:   Input Parameters:
5957: + dm         - The underlying DM for the domain
5958: . dim        - The spatial dimension
5959: . numComp    - The number of components
5960: . isSimplex  - Flag for simplex reference cell, otherwise its a tensor product
5961: . prefix     - The options prefix, or NULL
5962: - qorder     - The quadrature order

5964:   Output Parameter:
5965: . fem - The PetscFE object

5967:   Level: beginner

5969: .keywords: PetscFE, finite element
5970: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
5971: @*/
5972: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
5973: {
5974:   PetscQuadrature q;
5975:   DM              K;
5976:   PetscSpace      P;
5977:   PetscDualSpace  Q;
5978:   PetscInt        order;
5979:   PetscErrorCode  ierr;

5982:   /* Create space */
5983:   PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
5984:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
5985:   PetscSpaceSetFromOptions(P);
5986:   PetscSpacePolynomialSetNumVariables(P, dim);
5987:   PetscSpaceSetUp(P);
5988:   PetscSpaceGetOrder(P, &order);
5989:   /* Create dual space */
5990:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
5991:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
5992:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
5993:   PetscDualSpaceSetDM(Q, K);
5994:   DMDestroy(&K);
5995:   PetscDualSpaceSetOrder(Q, order);
5996:   PetscDualSpaceSetFromOptions(Q);
5997:   PetscDualSpaceSetUp(Q);
5998:   /* Create element */
5999:   PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
6000:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
6001:   PetscFESetFromOptions(*fem);
6002:   PetscFESetBasisSpace(*fem, P);
6003:   PetscFESetDualSpace(*fem, Q);
6004:   PetscFESetNumComponents(*fem, numComp);
6005:   PetscFESetUp(*fem);
6006:   PetscSpaceDestroy(&P);
6007:   PetscDualSpaceDestroy(&Q);
6008:   /* Create quadrature (with specified order if given) */
6009:   if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
6010:   else           {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
6011:   PetscFESetQuadrature(*fem, q);
6012:   PetscQuadratureDestroy(&q);
6013:   return(0);
6014: }