Actual source code: dtfe.c

petsc-dev 2014-08-28
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>
 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:   if (!PetscSpaceRegisterAllCalled) {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:   PetscSpaceViewFromOptions - Processes command line options to determine if/how a PetscSpace is to be viewed.

214:   Collective on PetscSpace

216:   Input Parameters:
217: + sp   - the PetscSpace
218: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
219: - optionname - option to activate viewing

221:   Level: intermediate

223: .keywords: PetscSpace, view, options, database
224: .seealso: VecViewFromOptions(), MatViewFromOptions()
225: */
226: PetscErrorCode PetscSpaceViewFromOptions(PetscSpace sp, const char prefix[], const char optionname[])
227: {
228:   PetscViewer       viewer;
229:   PetscViewerFormat format;
230:   PetscBool         flg;
231:   PetscErrorCode    ierr;

234:   if (prefix) {
235:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), prefix, optionname, &viewer, &format, &flg);
236:   } else {
237:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), ((PetscObject) sp)->prefix, optionname, &viewer, &format, &flg);
238:   }
239:   if (flg) {
240:     PetscViewerPushFormat(viewer, format);
241:     PetscSpaceView(sp, viewer);
242:     PetscViewerPopFormat(viewer);
243:     PetscViewerDestroy(&viewer);
244:   }
245:   return(0);
246: }

250: /*@
251:   PetscSpaceSetFromOptions - sets parameters in a PetscSpace from the options database

253:   Collective on PetscSpace

255:   Input Parameter:
256: . sp - the PetscSpace object to set options for

258:   Options Database:
259: . -petscspace_order the approximation order of the space

261:   Level: developer

263: .seealso PetscSpaceView()
264: @*/
265: PetscErrorCode PetscSpaceSetFromOptions(PetscSpace sp)
266: {
267:   const char    *defaultType;
268:   char           name[256];
269:   PetscBool      flg;

274:   if (!((PetscObject) sp)->type_name) {
275:     defaultType = PETSCSPACEPOLYNOMIAL;
276:   } else {
277:     defaultType = ((PetscObject) sp)->type_name;
278:   }
279:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

281:   PetscObjectOptionsBegin((PetscObject) sp);
282:   PetscOptionsFList("-petscspace_type", "Linear space", "PetscSpaceSetType", PetscSpaceList, defaultType, name, 256, &flg);
283:   if (flg) {
284:     PetscSpaceSetType(sp, name);
285:   } else if (!((PetscObject) sp)->type_name) {
286:     PetscSpaceSetType(sp, defaultType);
287:   }
288:   PetscOptionsInt("-petscspace_order", "The approximation order", "PetscSpaceSetOrder", sp->order, &sp->order, NULL);
289:   if (sp->ops->setfromoptions) {
290:     (*sp->ops->setfromoptions)(sp);
291:   }
292:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
293:   PetscObjectProcessOptionsHandlers((PetscObject) sp);
294:   PetscOptionsEnd();
295:   PetscSpaceViewFromOptions(sp, NULL, "-petscspace_view");
296:   return(0);
297: }

301: /*@C
302:   PetscSpaceSetUp - Construct data structures for the PetscSpace

304:   Collective on PetscSpace

306:   Input Parameter:
307: . sp - the PetscSpace object to setup

309:   Level: developer

311: .seealso PetscSpaceView(), PetscSpaceDestroy()
312: @*/
313: PetscErrorCode PetscSpaceSetUp(PetscSpace sp)
314: {

319:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
320:   return(0);
321: }

325: /*@
326:   PetscSpaceDestroy - Destroys a PetscSpace object

328:   Collective on PetscSpace

330:   Input Parameter:
331: . sp - the PetscSpace object to destroy

333:   Level: developer

335: .seealso PetscSpaceView()
336: @*/
337: PetscErrorCode PetscSpaceDestroy(PetscSpace *sp)
338: {

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

345:   if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
346:   ((PetscObject) (*sp))->refct = 0;
347:   DMDestroy(&(*sp)->dm);

349:   (*(*sp)->ops->destroy)(*sp);
350:   PetscHeaderDestroy(sp);
351:   return(0);
352: }

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

359:   Collective on MPI_Comm

361:   Input Parameter:
362: . comm - The communicator for the PetscSpace object

364:   Output Parameter:
365: . sp - The PetscSpace object

367:   Level: beginner

369: .seealso: PetscSpaceSetType(), PETSCSPACEPOLYNOMIAL
370: @*/
371: PetscErrorCode PetscSpaceCreate(MPI_Comm comm, PetscSpace *sp)
372: {
373:   PetscSpace     s;

378:   PetscCitationsRegister(FECitation,&FEcite);
379:   *sp  = NULL;
380:   PetscFEInitializePackage();

382:   PetscHeaderCreate(s, _p_PetscSpace, struct _PetscSpaceOps, PETSCSPACE_CLASSID, "PetscSpace", "Linear Space", "PetscSpace", comm, PetscSpaceDestroy, PetscSpaceView);
383:   PetscMemzero(s->ops, sizeof(struct _PetscSpaceOps));

385:   s->order = 0;
386:   DMShellCreate(comm, &s->dm);

388:   *sp = s;
389:   return(0);
390: }

394: /* Dimension of the space, i.e. number of basis vectors */
395: PetscErrorCode PetscSpaceGetDimension(PetscSpace sp, PetscInt *dim)
396: {

402:   *dim = 0;
403:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
404:   return(0);
405: }

409: PetscErrorCode PetscSpaceGetOrder(PetscSpace sp, PetscInt *order)
410: {
414:   *order = sp->order;
415:   return(0);
416: }

420: PetscErrorCode PetscSpaceSetOrder(PetscSpace sp, PetscInt order)
421: {
424:   sp->order = order;
425:   return(0);
426: }

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

433:   Input Parameters:
434: + sp      - The PetscSpace
435: . npoints - The number of evaluation points
436: - points  - The point coordinates

438:   Output Parameters:
439: + B - The function evaluations in a npoints x nfuncs array
440: . D - The derivative evaluations in a npoints x nfuncs x dim array
441: - H - The second derivative evaluations in a npoints x nfuncs x dim x dim array

443:   Level: advanced

445: .seealso: PetscFEGetTabulation(), PetscFEGetDefaultTabulation(), PetscSpaceCreate()
446: @*/
447: PetscErrorCode PetscSpaceEvaluate(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
448: {

457:   if (sp->ops->evaluate) {(*sp->ops->evaluate)(sp, npoints, points, B, D, H);}
458:   return(0);
459: }

463: PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscSpace sp)
464: {
465:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
466:   PetscErrorCode   ierr;

469:   PetscOptionsHead("PetscSpace Polynomial Options");
470:   PetscOptionsInt("-petscspace_poly_num_variables", "The number of different variables, e.g. x and y", "PetscSpacePolynomialSetNumVariables", poly->numVariables, &poly->numVariables, NULL);
471:   PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
472:   PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
473:   PetscOptionsTail();
474:   return(0);
475: }

479: PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer viewer)
480: {
481:   PetscSpace_Poly  *poly = (PetscSpace_Poly *) sp->data;
482:   PetscViewerFormat format;
483:   PetscErrorCode    ierr;

486:   PetscViewerGetFormat(viewer, &format);
487:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
488:     if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
489:     else              {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
490:   } else {
491:     if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
492:     else              {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
493:   }
494:   return(0);
495: }

499: PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
500: {
501:   PetscBool      iascii;

507:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
508:   if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
509:   return(0);
510: }

514: PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
515: {
516:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
517:   PetscInt         ndegree = sp->order+1;
518:   PetscInt         deg;
519:   PetscErrorCode   ierr;

522:   PetscMalloc1(ndegree, &poly->degrees);
523:   for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
524:   return(0);
525: }

529: PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
530: {
531:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
532:   PetscErrorCode   ierr;

535:   PetscFree(poly->degrees);
536:   PetscFree(poly);
537:   return(0);
538: }

542: PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
543: {
544:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
545:   PetscInt         deg  = sp->order;
546:   PetscInt         n    = poly->numVariables, i;
547:   PetscReal        D    = 1.0;

550:   if (poly->tensor) {
551:     *dim = 1;
552:     for (i = 0; i < n; ++i) *dim *= (deg+1);
553:   } else {
554:     for (i = 1; i <= n; ++i) {
555:       D *= ((PetscReal) (deg+i))/i;
556:     }
557:     *dim = (PetscInt) (D + 0.5);
558:   }
559:   return(0);
560: }

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

567:   Input Parameters:
568: + len - The length of the tuple
569: . sum - The sum of all entries in the tuple
570: - ind - The current multi-index of the tuple, initialized to the 0 tuple

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

576:   Level: developer

578: .seealso: 
579: */
580: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
581: {
582:   PetscInt       i;

586:   if (len == 1) {
587:     ind[0] = -1;
588:     tup[0] = sum;
589:   } else if (sum == 0) {
590:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
591:   } else {
592:     tup[0] = sum - ind[0];
593:     LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
594:     if (ind[1] < 0) {
595:       if (ind[0] == sum) {ind[0] = -1;}
596:       else               {ind[1] = 0; ++ind[0];}
597:     }
598:   }
599:   return(0);
600: }

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

607:   Input Parameters:
608: + len - The length of the tuple
609: . max - The max for all entries in the tuple
610: - ind - The current multi-index of the tuple, initialized to the 0 tuple

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

616:   Level: developer

618: .seealso: 
619: */
620: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
621: {
622:   PetscInt       i;

626:   if (len == 1) {
627:     tup[0] = ind[0]++;
628:     ind[0] = ind[0] >= max ? -1 : ind[0];
629:   } else if (max == 0) {
630:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
631:   } else {
632:     tup[0] = ind[0];
633:     TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
634:     if (ind[1] < 0) {
635:       ind[1] = 0;
636:       if (ind[0] == max-1) {ind[0] = -1;}
637:       else                 {++ind[0];}
638:     }
639:   }
640:   return(0);
641: }

645: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
646: {
647:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
648:   DM               dm      = sp->dm;
649:   PetscInt         ndegree = sp->order+1;
650:   PetscInt        *degrees = poly->degrees;
651:   PetscInt         dim     = poly->numVariables;
652:   PetscReal       *lpoints, *tmp, *LB, *LD, *LH;
653:   PetscInt        *ind, *tup;
654:   PetscInt         pdim, d, der, i, p, deg, o;
655:   PetscErrorCode   ierr;

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

768: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
769: {
771:   sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
772:   sp->ops->setup          = PetscSpaceSetUp_Polynomial;
773:   sp->ops->view           = PetscSpaceView_Polynomial;
774:   sp->ops->destroy        = PetscSpaceDestroy_Polynomial;
775:   sp->ops->getdimension   = PetscSpaceGetDimension_Polynomial;
776:   sp->ops->evaluate       = PetscSpaceEvaluate_Polynomial;
777:   return(0);
778: }

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

783:   Level: intermediate

785: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
786: M*/

790: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
791: {
792:   PetscSpace_Poly *poly;
793:   PetscErrorCode   ierr;

797:   PetscNewLog(sp,&poly);
798:   sp->data = poly;

800:   poly->numVariables = 0;
801:   poly->symmetric    = PETSC_FALSE;
802:   poly->tensor       = PETSC_FALSE;
803:   poly->degrees      = NULL;

805:   PetscSpaceInitialize_Polynomial(sp);
806:   return(0);
807: }

811: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
812: {
813:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

817:   poly->symmetric = sym;
818:   return(0);
819: }

823: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
824: {
825:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

830:   *sym = poly->symmetric;
831:   return(0);
832: }

836: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
837: {
838:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

842:   poly->tensor = tensor;
843:   return(0);
844: }

848: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
849: {
850:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

855:   *tensor = poly->tensor;
856:   return(0);
857: }

861: PetscErrorCode PetscSpacePolynomialSetNumVariables(PetscSpace sp, PetscInt n)
862: {
863:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

867:   poly->numVariables = n;
868:   return(0);
869: }

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

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

886: PetscErrorCode PetscSpaceSetFromOptions_DG(PetscSpace sp)
887: {
888:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

892:   PetscOptionsHead("PetscSpace DG Options");
893:   PetscOptionsInt("-petscspace_dg_num_variables", "The number of different variables, e.g. x and y", "PetscSpaceDGSetNumVariables", dg->numVariables, &dg->numVariables, NULL);
894:   PetscOptionsTail();
895:   return(0);
896: }

900: PetscErrorCode PetscSpaceDGView_Ascii(PetscSpace sp, PetscViewer viewer)
901: {
902:   PetscSpace_DG    *dg = (PetscSpace_DG *) sp->data;
903:   PetscViewerFormat format;
904:   PetscErrorCode    ierr;

907:   PetscViewerGetFormat(viewer, &format);
908:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
909:     PetscViewerASCIIPrintf(viewer, "DG space in dimension %d:\n", dg->numVariables);
910:     PetscViewerASCIIPushTab(viewer);
911:     PetscQuadratureView(dg->quad, viewer);
912:     PetscViewerASCIIPopTab(viewer);
913:   } else {
914:     PetscViewerASCIIPrintf(viewer, "DG space in dimension %d on %d points\n", dg->numVariables, dg->quad->numPoints);
915:   }
916:   return(0);
917: }

921: PetscErrorCode PetscSpaceView_DG(PetscSpace sp, PetscViewer viewer)
922: {
923:   PetscBool      iascii;

929:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
930:   if (iascii) {PetscSpaceDGView_Ascii(sp, viewer);}
931:   return(0);
932: }

936: PetscErrorCode PetscSpaceSetUp_DG(PetscSpace sp)
937: {
938:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

942:   if (!dg->quad->points && sp->order) {
943:     PetscDTGaussJacobiQuadrature(dg->numVariables, sp->order, -1.0, 1.0, &dg->quad);
944:   }
945:   return(0);
946: }

950: PetscErrorCode PetscSpaceDestroy_DG(PetscSpace sp)
951: {
952:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

956:   PetscQuadratureDestroy(&dg->quad);
957:   return(0);
958: }

962: PetscErrorCode PetscSpaceGetDimension_DG(PetscSpace sp, PetscInt *dim)
963: {
964:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

967:   *dim = dg->quad->numPoints;
968:   return(0);
969: }

973: PetscErrorCode PetscSpaceEvaluate_DG(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
974: {
975:   PetscSpace_DG *dg  = (PetscSpace_DG *) sp->data;
976:   PetscInt       dim = dg->numVariables, d, p;

980:   if (D || H) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Cannot calculate derivatives for a DG space");
981:   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);
982:   PetscMemzero(B, npoints*npoints * sizeof(PetscReal));
983:   for (p = 0; p < npoints; ++p) {
984:     for (d = 0; d < dim; ++d) {
985:       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]);
986:     }
987:     B[p*npoints+p] = 1.0;
988:   }
989:   return(0);
990: }

994: PetscErrorCode PetscSpaceInitialize_DG(PetscSpace sp)
995: {
997:   sp->ops->setfromoptions = PetscSpaceSetFromOptions_DG;
998:   sp->ops->setup          = PetscSpaceSetUp_DG;
999:   sp->ops->view           = PetscSpaceView_DG;
1000:   sp->ops->destroy        = PetscSpaceDestroy_DG;
1001:   sp->ops->getdimension   = PetscSpaceGetDimension_DG;
1002:   sp->ops->evaluate       = PetscSpaceEvaluate_DG;
1003:   return(0);
1004: }

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

1009:   Level: intermediate

1011: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
1012: M*/

1016: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_DG(PetscSpace sp)
1017: {
1018:   PetscSpace_DG *dg;

1023:   PetscNewLog(sp,&dg);
1024:   sp->data = dg;

1026:   dg->numVariables    = 0;
1027:   dg->quad->dim       = 0;
1028:   dg->quad->numPoints = 0;
1029:   dg->quad->points    = NULL;
1030:   dg->quad->weights   = NULL;

1032:   PetscSpaceInitialize_DG(sp);
1033:   return(0);
1034: }


1037: PetscClassId PETSCDUALSPACE_CLASSID = 0;

1039: PetscFunctionList PetscDualSpaceList              = NULL;
1040: PetscBool         PetscDualSpaceRegisterAllCalled = PETSC_FALSE;

1044: /*@C
1045:   PetscDualSpaceRegister - Adds a new PetscDualSpace implementation

1047:   Not Collective

1049:   Input Parameters:
1050: + name        - The name of a new user-defined creation routine
1051: - create_func - The creation routine itself

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

1056:   Sample usage:
1057: .vb
1058:     PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
1059: .ve

1061:   Then, your PetscDualSpace type can be chosen with the procedural interface via
1062: .vb
1063:     PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
1064:     PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
1065: .ve
1066:    or at runtime via the option
1067: .vb
1068:     -petscdualspace_type my_dual_space
1069: .ve

1071:   Level: advanced

1073: .keywords: PetscDualSpace, register
1074: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()

1076: @*/
1077: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
1078: {

1082:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
1083:   return(0);
1084: }

1088: /*@C
1089:   PetscDualSpaceSetType - Builds a particular PetscDualSpace

1091:   Collective on PetscDualSpace

1093:   Input Parameters:
1094: + sp   - The PetscDualSpace object
1095: - name - The kind of space

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

1100:   Level: intermediate

1102: .keywords: PetscDualSpace, set, type
1103: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
1104: @*/
1105: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
1106: {
1107:   PetscErrorCode (*r)(PetscDualSpace);
1108:   PetscBool      match;

1113:   PetscObjectTypeCompare((PetscObject) sp, name, &match);
1114:   if (match) return(0);

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

1120:   if (sp->ops->destroy) {
1121:     (*sp->ops->destroy)(sp);
1122:     sp->ops->destroy = NULL;
1123:   }
1124:   (*r)(sp);
1125:   PetscObjectChangeTypeName((PetscObject) sp, name);
1126:   return(0);
1127: }

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

1134:   Not Collective

1136:   Input Parameter:
1137: . sp  - The PetscDualSpace

1139:   Output Parameter:
1140: . name - The PetscDualSpace type name

1142:   Level: intermediate

1144: .keywords: PetscDualSpace, get, type, name
1145: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
1146: @*/
1147: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
1148: {

1154:   if (!PetscDualSpaceRegisterAllCalled) {
1155:     PetscDualSpaceRegisterAll();
1156:   }
1157:   *name = ((PetscObject) sp)->type_name;
1158:   return(0);
1159: }

1163: /*@
1164:   PetscDualSpaceView - Views a PetscDualSpace

1166:   Collective on PetscDualSpace

1168:   Input Parameter:
1169: + sp - the PetscDualSpace object to view
1170: - v  - the viewer

1172:   Level: developer

1174: .seealso PetscDualSpaceDestroy()
1175: @*/
1176: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
1177: {

1182:   if (!v) {
1183:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
1184:   }
1185:   if (sp->ops->view) {
1186:     (*sp->ops->view)(sp, v);
1187:   }
1188:   return(0);
1189: }

1193: /*@C
1194:   PetscDualSpaceViewFromOptions - Processes command line options to determine if/how a PetscDualSpace is to be viewed.

1196:   Collective on PetscDualSpace

1198:   Input Parameters:
1199: + sp   - the PetscDualSpace
1200: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
1201: - optionname - option to activate viewing

1203:   Level: intermediate

1205: .keywords: PetscDualSpace, view, options, database
1206: .seealso: VecViewFromOptions(), MatViewFromOptions()
1207: @*/
1208: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace sp, const char prefix[], const char optionname[])
1209: {
1210:   PetscViewer       viewer;
1211:   PetscViewerFormat format;
1212:   PetscBool         flg;
1213:   PetscErrorCode    ierr;

1216:   if (prefix) {
1217:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), prefix, optionname, &viewer, &format, &flg);
1218:   } else {
1219:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), ((PetscObject) sp)->prefix, optionname, &viewer, &format, &flg);
1220:   }
1221:   if (flg) {
1222:     PetscViewerPushFormat(viewer, format);
1223:     PetscDualSpaceView(sp, viewer);
1224:     PetscViewerPopFormat(viewer);
1225:     PetscViewerDestroy(&viewer);
1226:   }
1227:   return(0);
1228: }

1232: /*@
1233:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

1235:   Collective on PetscDualSpace

1237:   Input Parameter:
1238: . sp - the PetscDualSpace object to set options for

1240:   Options Database:
1241: . -petscspace_order the approximation order of the space

1243:   Level: developer

1245: .seealso PetscDualSpaceView()
1246: @*/
1247: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
1248: {
1249:   const char    *defaultType;
1250:   char           name[256];
1251:   PetscBool      flg;

1256:   if (!((PetscObject) sp)->type_name) {
1257:     defaultType = PETSCDUALSPACELAGRANGE;
1258:   } else {
1259:     defaultType = ((PetscObject) sp)->type_name;
1260:   }
1261:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

1263:   PetscObjectOptionsBegin((PetscObject) sp);
1264:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
1265:   if (flg) {
1266:     PetscDualSpaceSetType(sp, name);
1267:   } else if (!((PetscObject) sp)->type_name) {
1268:     PetscDualSpaceSetType(sp, defaultType);
1269:   }
1270:   PetscOptionsInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
1271:   if (sp->ops->setfromoptions) {
1272:     (*sp->ops->setfromoptions)(sp);
1273:   }
1274:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1275:   PetscObjectProcessOptionsHandlers((PetscObject) sp);
1276:   PetscOptionsEnd();
1277:   PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
1278:   return(0);
1279: }

1283: /*@
1284:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

1286:   Collective on PetscDualSpace

1288:   Input Parameter:
1289: . sp - the PetscDualSpace object to setup

1291:   Level: developer

1293: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
1294: @*/
1295: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
1296: {

1301:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
1302:   return(0);
1303: }

1307: /*@
1308:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

1310:   Collective on PetscDualSpace

1312:   Input Parameter:
1313: . sp - the PetscDualSpace object to destroy

1315:   Level: developer

1317: .seealso PetscDualSpaceView()
1318: @*/
1319: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
1320: {
1321:   PetscInt       dim, f;

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

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

1331:   PetscDualSpaceGetDimension(*sp, &dim);
1332:   for (f = 0; f < dim; ++f) {
1333:     PetscQuadratureDestroy(&(*sp)->functional[f]);
1334:   }
1335:   PetscFree((*sp)->functional);
1336:   DMDestroy(&(*sp)->dm);

1338:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
1339:   PetscHeaderDestroy(sp);
1340:   return(0);
1341: }

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

1348:   Collective on MPI_Comm

1350:   Input Parameter:
1351: . comm - The communicator for the PetscDualSpace object

1353:   Output Parameter:
1354: . sp - The PetscDualSpace object

1356:   Level: beginner

1358: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
1359: @*/
1360: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
1361: {
1362:   PetscDualSpace s;

1367:   PetscCitationsRegister(FECitation,&FEcite);
1368:   *sp  = NULL;
1369:   PetscFEInitializePackage();

1371:   PetscHeaderCreate(s, _p_PetscDualSpace, struct _PetscDualSpaceOps, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
1372:   PetscMemzero(s->ops, sizeof(struct _PetscDualSpaceOps));

1374:   s->order = 0;

1376:   *sp = s;
1377:   return(0);
1378: }

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

1385:   Collective on PetscDualSpace

1387:   Input Parameter:
1388: . sp - The original PetscDualSpace

1390:   Output Parameter:
1391: . spNew - The duplicate PetscDualSpace

1393:   Level: beginner

1395: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
1396: @*/
1397: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
1398: {

1404:   (*sp->ops->duplicate)(sp, spNew);
1405:   return(0);
1406: }

1410: /*@
1411:   PetscDualSpaceGetDM - Get the DM representing the reference cell

1413:   Not collective

1415:   Input Parameter:
1416: . sp - The PetscDualSpace

1418:   Output Parameter:
1419: . dm - The reference cell

1421:   Level: intermediate

1423: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
1424: @*/
1425: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
1426: {
1430:   *dm = sp->dm;
1431:   return(0);
1432: }

1436: /*@
1437:   PetscDualSpaceSetDM - Get the DM representing the reference cell

1439:   Not collective

1441:   Input Parameters:
1442: + sp - The PetscDualSpace
1443: - dm - The reference cell

1445:   Level: intermediate

1447: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
1448: @*/
1449: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
1450: {

1456:   DMDestroy(&sp->dm);
1457:   PetscObjectReference((PetscObject) dm);
1458:   sp->dm = dm;
1459:   return(0);
1460: }

1464: /*@
1465:   PetscDualSpaceGetOrder - Get the order of the dual space

1467:   Not collective

1469:   Input Parameter:
1470: . sp - The PetscDualSpace

1472:   Output Parameter:
1473: . order - The order

1475:   Level: intermediate

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

1490: /*@
1491:   PetscDualSpaceSetOrder - Set the order of the dual space

1493:   Not collective

1495:   Input Parameters:
1496: + sp - The PetscDualSpace
1497: - order - The order

1499:   Level: intermediate

1501: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
1502: @*/
1503: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
1504: {
1507:   sp->order = order;
1508:   return(0);
1509: }

1513: /*@
1514:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

1516:   Not collective

1518:   Input Parameters:
1519: + sp - The PetscDualSpace
1520: - i  - The basis number

1522:   Output Parameter:
1523: . functional - The basis functional

1525:   Level: intermediate

1527: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
1528: @*/
1529: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
1530: {
1531:   PetscInt       dim;

1537:   PetscDualSpaceGetDimension(sp, &dim);
1538:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
1539:   *functional = sp->functional[i];
1540:   return(0);
1541: }

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

1548:   Not collective

1550:   Input Parameter:
1551: . sp - The PetscDualSpace

1553:   Output Parameter:
1554: . dim - The dimension

1556:   Level: intermediate

1558: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1559: @*/
1560: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
1561: {

1567:   *dim = 0;
1568:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
1569:   return(0);
1570: }

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

1577:   Not collective

1579:   Input Parameter:
1580: . sp - The PetscDualSpace

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

1585:   Level: intermediate

1587: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1588: @*/
1589: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
1590: {

1596:   *numDof = NULL;
1597:   if (sp->ops->getnumdof) {(*sp->ops->getnumdof)(sp, numDof);}
1598:   return(0);
1599: }

1603: /*@
1604:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

1606:   Collective on PetscDualSpace

1608:   Input Parameters:
1609: + sp      - The PetscDualSpace
1610: . dim     - The spatial dimension
1611: - simplex - Flag for simplex, otherwise use a tensor-product cell

1613:   Output Parameter:
1614: . refdm - The reference cell

1616:   Level: advanced

1618: .keywords: PetscDualSpace, reference cell
1619: .seealso: PetscDualSpaceCreate(), DMPLEX
1620: @*/
1621: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1622: {

1626:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
1627:   return(0);
1628: }

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

1635:   Input Parameters:
1636: + sp      - The PetscDualSpace object
1637: . f       - The basis functional index
1638: . geom    - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1639: . numComp - The number of components for the function
1640: . func    - The input function
1641: - ctx     - A context for the function

1643:   Output Parameter:
1644: . value   - numComp output values

1646:   Level: developer

1648: .seealso: PetscDualSpaceCreate()
1649: @*/
1650: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscCellGeometry geom, PetscInt numComp, void (*func)(const PetscReal [], PetscScalar *, void *), void *ctx, PetscScalar *value)
1651: {
1652:   DM               dm;
1653:   PetscQuadrature  quad;
1654:   const PetscReal *v0 = geom.v0;
1655:   const PetscReal *J  = geom.J;
1656:   PetscReal        x[3];
1657:   PetscScalar     *val;
1658:   PetscInt         dim, q, c;
1659:   PetscErrorCode   ierr;

1664:   PetscDualSpaceGetDM(sp, &dm);
1665:   DMGetDimension(dm, &dim);
1666:   PetscDualSpaceGetFunctional(sp, f, &quad);
1667:   DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1668:   for (c = 0; c < numComp; ++c) value[c] = 0.0;
1669:   for (q = 0; q < quad->numPoints; ++q) {
1670:     CoordinatesRefToReal(dim, dim, v0, J, &quad->points[q*dim], x);
1671:     (*func)(x, val, ctx);
1672:     for (c = 0; c < numComp; ++c) {
1673:       value[c] += val[c]*quad->weights[q];
1674:     }
1675:   }
1676:   DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1677:   return(0);
1678: }

1682: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
1683: {
1684:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1685:   PetscReal           D   = 1.0;
1686:   PetscInt            n, i;
1687:   PetscErrorCode      ierr;

1690:   DMGetDimension(sp->dm, &n);
1691:   if (lag->simplex || !lag->continuous) {
1692:     for (i = 1; i <= n; ++i) {
1693:       D *= ((PetscReal) (order+i))/i;
1694:     }
1695:     *dim = (PetscInt) (D + 0.5);
1696:   } else {
1697:     *dim = 1;
1698:     for (i = 0; i < n; ++i) *dim *= (order+1);
1699:   }
1700:   return(0);
1701: }

1705: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1706: {
1707:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1708:   DM                  dm    = sp->dm;
1709:   PetscInt            order = sp->order;
1710:   PetscBool           disc  = lag->continuous ? PETSC_FALSE : PETSC_TRUE;
1711:   PetscSection        csection;
1712:   Vec                 coordinates;
1713:   PetscReal          *qpoints, *qweights;
1714:   PetscInt           *closure = NULL, closureSize, c;
1715:   PetscInt            depth, dim, pdimMax, pMax = 0, *pStart, *pEnd, cell, coneSize, d, n, f = 0;
1716:   PetscBool           simplex;
1717:   PetscErrorCode      ierr;

1720:   /* Classify element type */
1721:   DMGetDimension(dm, &dim);
1722:   DMPlexGetDepth(dm, &depth);
1723:   PetscCalloc1(dim+1, &lag->numDof);
1724:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
1725:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1726:   DMPlexGetConeSize(dm, pStart[depth], &coneSize);
1727:   DMGetCoordinateSection(dm, &csection);
1728:   DMGetCoordinatesLocal(dm, &coordinates);
1729:   if (depth == 1) {
1730:     if      (coneSize == dim+1)    simplex = PETSC_TRUE;
1731:     else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
1732:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1733:   }
1734:   else if (depth == dim) {
1735:     if      (coneSize == dim+1)   simplex = PETSC_TRUE;
1736:     else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
1737:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1738:   }
1739:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
1740:   lag->simplex = simplex;
1741:   PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
1742:   pdimMax *= (pEnd[dim] - pStart[dim]);
1743:   PetscMalloc1(pdimMax, &sp->functional);
1744:   for (d = 0; d <= depth; d++) {
1745:     pMax = PetscMax(pMax,pEnd[d]);
1746:   }
1747:   if (!dim) {
1748:     PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1749:     PetscMalloc1(1, &qpoints);
1750:     PetscMalloc1(1, &qweights);
1751:     PetscQuadratureSetOrder(sp->functional[f], 0);
1752:     PetscQuadratureSetData(sp->functional[f], PETSC_DETERMINE, 1, qpoints, qweights);
1753:     qpoints[0]  = 0.0;
1754:     qweights[0] = 1.0;
1755:     ++f;
1756:     lag->numDof[0] = 1;
1757:   } else {
1758:     PetscBT seen;

1760:     PetscBTCreate(pMax, &seen);
1761:     PetscBTMemzero(pMax, seen);
1762:     for (cell = pStart[depth]; cell < pEnd[depth]; ++cell) {
1763:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1764:       for (c = 0; c < closureSize*2; c += 2) {
1765:         const PetscInt p = closure[c];

1767:         if (PetscBTLookup(seen, p)) continue;
1768:         PetscBTSet(seen, p);
1769:         if ((p >= pStart[0]) && (p < pEnd[0])) {
1770:           /* Vertices */
1771:           const PetscScalar *coords;
1772:           PetscInt           dof, off, d;

1774:           if (order < 1 || disc) continue;
1775:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1776:           PetscMalloc1(dim, &qpoints);
1777:           PetscMalloc1(1, &qweights);
1778:           PetscQuadratureSetOrder(sp->functional[f], 0);
1779:           PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1780:           VecGetArrayRead(coordinates, &coords);
1781:           PetscSectionGetDof(csection, p, &dof);
1782:           PetscSectionGetOffset(csection, p, &off);
1783:           if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim);
1784:           for (d = 0; d < dof; ++d) {qpoints[d] = PetscRealPart(coords[off+d]);}
1785:           qweights[0] = 1.0;
1786:           ++f;
1787:           VecRestoreArrayRead(coordinates, &coords);
1788:           lag->numDof[0] = 1;
1789:         } else if ((p >= pStart[1]) && (p < pEnd[1])) {
1790:           /* Edges */
1791:           PetscScalar *coords;
1792:           PetscInt     num = order-1, k;

1794:           if (order < 2 || disc) continue;
1795:           coords = NULL;
1796:           DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);
1797:           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);
1798:           for (k = 1; k <= num; ++k) {
1799:             PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1800:             PetscMalloc1(dim, &qpoints);
1801:             PetscMalloc1(1, &qweights);
1802:             PetscQuadratureSetOrder(sp->functional[f], 0);
1803:             PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1804:             for (d = 0; d < dim; ++d) {qpoints[d] = k*PetscRealPart(coords[1*dim+d] - coords[0*dim+d])/order + PetscRealPart(coords[0*dim+d]);}
1805:             qweights[0] = 1.0;
1806:             ++f;
1807:           }
1808:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);
1809:           lag->numDof[1] = num;
1810:         } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) {
1811:           /* Faces */

1813:           if (disc) continue;
1814:           if ( simplex && (order < 3)) continue;
1815:           if (!simplex && (order < 2)) continue;
1816:           lag->numDof[depth-1] = 0;
1817:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces");
1818:         } else if ((p >= pStart[depth]) && (p < pEnd[depth])) {
1819:           /* Cells */
1820:           PetscInt     orderEff = lag->continuous && order ? (simplex ? order-3 : order-2) : order;
1821:           PetscReal    denom    = order ? (lag->continuous ? order : (simplex ? order+3 : order+2)) : (simplex ? 3 : 2);
1822:           PetscScalar *coords   = NULL;
1823:           PetscReal    dx = 2.0/denom, *v0, *J, *invJ, detJ;
1824:           PetscInt    *ind, *tup;
1825:           PetscInt     cdim, csize, d, d2, o;

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

1833:           PetscCalloc5(dim,&ind,dim,&tup,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1834:           DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
1835:           if (simplex || !lag->continuous) {
1836:             for (o = 0; o <= orderEff; ++o) {
1837:               PetscMemzero(ind, dim*sizeof(PetscInt));
1838:               while (ind[0] >= 0) {
1839:                 PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1840:                 PetscMalloc1(dim, &qpoints);
1841:                 PetscMalloc1(1,   &qweights);
1842:                 PetscQuadratureSetOrder(sp->functional[f], 0);
1843:                 PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1844:                 LatticePoint_Internal(dim, o, ind, tup);
1845:                 for (d = 0; d < dim; ++d) {
1846:                   qpoints[d] = v0[d];
1847:                   for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
1848:                 }
1849:                 qweights[0] = 1.0;
1850:                 ++f;
1851:               }
1852:             }
1853:           } else {
1854:             while (ind[0] >= 0) {
1855:               PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1856:               PetscMalloc1(dim, &qpoints);
1857:               PetscMalloc1(1,   &qweights);
1858:               PetscQuadratureSetOrder(sp->functional[f], 0);
1859:               PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1860:               TensorPoint_Internal(dim, orderEff+1, ind, tup);
1861:               for (d = 0; d < dim; ++d) {
1862:                 qpoints[d] = v0[d];
1863:                 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d]+1)*dx);
1864:               }
1865:               qweights[0] = 1.0;
1866:               ++f;
1867:             }
1868:           }
1869:           PetscFree5(ind,tup,v0,J,invJ);
1870:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &csize, &coords);
1871:           lag->numDof[depth] = cdim;
1872:         }
1873:       }
1874:       DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);
1875:     }
1876:     PetscBTDestroy(&seen);
1877:   }
1878:   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);
1879:   PetscFree2(pStart,pEnd);
1880:   if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
1881:   return(0);
1882: }

1886: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
1887: {
1888:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1889:   PetscErrorCode      ierr;

1892:   PetscFree(lag->numDof);
1893:   PetscFree(lag);
1894:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
1895:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
1896:   return(0);
1897: }

1901: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
1902: {
1903:   PetscInt       order;
1904:   PetscBool      cont;

1908:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
1909:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
1910:   PetscDualSpaceGetOrder(sp, &order);
1911:   PetscDualSpaceSetOrder(*spNew, order);
1912:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
1913:   PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
1914:   return(0);
1915: }

1919: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscDualSpace sp)
1920: {
1921:   PetscBool      continuous, flg;

1925:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
1926:   PetscOptionsHead("PetscDualSpace Lagrange Options");
1927:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
1928:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
1929:   PetscOptionsTail();
1930:   return(0);
1931: }

1935: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
1936: {
1937:   DM              K;
1938:   const PetscInt *numDof;
1939:   PetscInt        spatialDim, Nc, size = 0, d;
1940:   PetscErrorCode  ierr;

1943:   PetscDualSpaceGetDM(sp, &K);
1944:   PetscDualSpaceGetNumDof(sp, &numDof);
1945:   DMGetDimension(K, &spatialDim);
1946:   DMPlexGetHeightStratum(K, 0, NULL, &Nc);
1947:   if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
1948:   for (d = 0; d <= spatialDim; ++d) {
1949:     PetscInt pStart, pEnd;

1951:     DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
1952:     size += (pEnd-pStart)*numDof[d];
1953:   }
1954:   *dim = size;
1955:   return(0);
1956: }

1960: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
1961: {
1962:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

1965:   *numDof = lag->numDof;
1966:   return(0);
1967: }

1971: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
1972: {
1973:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

1978:   *continuous = lag->continuous;
1979:   return(0);
1980: }

1984: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
1985: {
1986:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

1990:   lag->continuous = continuous;
1991:   return(0);
1992: }

1996: /*@
1997:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

1999:   Not Collective

2001:   Input Parameter:
2002: . sp         - the PetscDualSpace

2004:   Output Parameter:
2005: . continuous - flag for element continuity

2007:   Level: intermediate

2009: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2010: .seealso: PetscDualSpaceLagrangeSetContinuity()
2011: @*/
2012: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2013: {

2019:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2020:   return(0);
2021: }

2025: /*@
2026:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2028:   Logically Collective on PetscDualSpace

2030:   Input Parameters:
2031: + sp         - the PetscDualSpace
2032: - continuous - flag for element continuity

2034:   Options Database:
2035: . -petscdualspace_lagrange_continuity <bool>

2037:   Level: intermediate

2039: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2040: .seealso: PetscDualSpaceLagrangeGetContinuity()
2041: @*/
2042: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2043: {

2049:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2050:   return(0);
2051: }

2055: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2056: {
2058:   sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
2059:   sp->ops->setup          = PetscDualSpaceSetUp_Lagrange;
2060:   sp->ops->view           = NULL;
2061:   sp->ops->destroy        = PetscDualSpaceDestroy_Lagrange;
2062:   sp->ops->duplicate      = PetscDualSpaceDuplicate_Lagrange;
2063:   sp->ops->getdimension   = PetscDualSpaceGetDimension_Lagrange;
2064:   sp->ops->getnumdof      = PetscDualSpaceGetNumDof_Lagrange;
2065:   return(0);
2066: }

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

2071:   Level: intermediate

2073: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2074: M*/

2078: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2079: {
2080:   PetscDualSpace_Lag *lag;
2081:   PetscErrorCode      ierr;

2085:   PetscNewLog(sp,&lag);
2086:   sp->data = lag;

2088:   lag->numDof     = NULL;
2089:   lag->simplex    = PETSC_TRUE;
2090:   lag->continuous = PETSC_TRUE;

2092:   PetscDualSpaceInitialize_Lagrange(sp);
2093:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2094:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2095:   return(0);
2096: }


2099: PetscClassId PETSCFE_CLASSID = 0;

2101: PetscFunctionList PetscFEList              = NULL;
2102: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

2106: /*@C
2107:   PetscFERegister - Adds a new PetscFE implementation

2109:   Not Collective

2111:   Input Parameters:
2112: + name        - The name of a new user-defined creation routine
2113: - create_func - The creation routine itself

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

2118:   Sample usage:
2119: .vb
2120:     PetscFERegister("my_fe", MyPetscFECreate);
2121: .ve

2123:   Then, your PetscFE type can be chosen with the procedural interface via
2124: .vb
2125:     PetscFECreate(MPI_Comm, PetscFE *);
2126:     PetscFESetType(PetscFE, "my_fe");
2127: .ve
2128:    or at runtime via the option
2129: .vb
2130:     -petscfe_type my_fe
2131: .ve

2133:   Level: advanced

2135: .keywords: PetscFE, register
2136: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()

2138: @*/
2139: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2140: {

2144:   PetscFunctionListAdd(&PetscFEList, sname, function);
2145:   return(0);
2146: }

2150: /*@C
2151:   PetscFESetType - Builds a particular PetscFE

2153:   Collective on PetscFE

2155:   Input Parameters:
2156: + fem  - The PetscFE object
2157: - name - The kind of FEM space

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

2162:   Level: intermediate

2164: .keywords: PetscFE, set, type
2165: .seealso: PetscFEGetType(), PetscFECreate()
2166: @*/
2167: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2168: {
2169:   PetscErrorCode (*r)(PetscFE);
2170:   PetscBool      match;

2175:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
2176:   if (match) return(0);

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

2182:   if (fem->ops->destroy) {
2183:     (*fem->ops->destroy)(fem);
2184:     fem->ops->destroy = NULL;
2185:   }
2186:   (*r)(fem);
2187:   PetscObjectChangeTypeName((PetscObject) fem, name);
2188:   return(0);
2189: }

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

2196:   Not Collective

2198:   Input Parameter:
2199: . fem  - The PetscFE

2201:   Output Parameter:
2202: . name - The PetscFE type name

2204:   Level: intermediate

2206: .keywords: PetscFE, get, type, name
2207: .seealso: PetscFESetType(), PetscFECreate()
2208: @*/
2209: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2210: {

2216:   if (!PetscFERegisterAllCalled) {
2217:     PetscFERegisterAll();
2218:   }
2219:   *name = ((PetscObject) fem)->type_name;
2220:   return(0);
2221: }

2225: /*@C
2226:   PetscFEView - Views a PetscFE

2228:   Collective on PetscFE

2230:   Input Parameter:
2231: + fem - the PetscFE object to view
2232: - v   - the viewer

2234:   Level: developer

2236: .seealso PetscFEDestroy()
2237: @*/
2238: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
2239: {

2244:   if (!v) {
2245:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
2246:   }
2247:   if (fem->ops->view) {
2248:     (*fem->ops->view)(fem, v);
2249:   }
2250:   return(0);
2251: }

2255: /*
2256:   PetscFEViewFromOptions - Processes command line options to determine if/how a PetscFE is to be viewed.

2258:   Collective on PetscFE

2260:   Input Parameters:
2261: + fem    - the PetscFE
2262: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
2263: - optionname - option to activate viewing

2265:   Level: intermediate

2267: .keywords: PetscFE, view, options, database
2268: .seealso: VecViewFromOptions(), MatViewFromOptions()
2269: */
2270: PetscErrorCode PetscFEViewFromOptions(PetscFE fem, const char prefix[], const char optionname[])
2271: {
2272:   PetscViewer       viewer;
2273:   PetscViewerFormat format;
2274:   PetscBool         flg;
2275:   PetscErrorCode    ierr;

2278:   if (prefix) {
2279:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), prefix, optionname, &viewer, &format, &flg);
2280:   } else {
2281:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), ((PetscObject) fem)->prefix, optionname, &viewer, &format, &flg);
2282:   }
2283:   if (flg) {
2284:     PetscViewerPushFormat(viewer, format);
2285:     PetscFEView(fem, viewer);
2286:     PetscViewerPopFormat(viewer);
2287:     PetscViewerDestroy(&viewer);
2288:   }
2289:   return(0);
2290: }

2294: /*@
2295:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

2297:   Collective on PetscFE

2299:   Input Parameter:
2300: . fem - the PetscFE object to set options for

2302:   Options Database:
2303: . -petscfe_num_blocks  the number of cell blocks to integrate concurrently
2304: . -petscfe_num_batches the number of cell batches to integrate serially

2306:   Level: developer

2308: .seealso PetscFEView()
2309: @*/
2310: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
2311: {
2312:   const char    *defaultType;
2313:   char           name[256];
2314:   PetscBool      flg;

2319:   if (!((PetscObject) fem)->type_name) {
2320:     defaultType = PETSCFEBASIC;
2321:   } else {
2322:     defaultType = ((PetscObject) fem)->type_name;
2323:   }
2324:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

2326:   PetscObjectOptionsBegin((PetscObject) fem);
2327:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
2328:   if (flg) {
2329:     PetscFESetType(fem, name);
2330:   } else if (!((PetscObject) fem)->type_name) {
2331:     PetscFESetType(fem, defaultType);
2332:   }
2333:   PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
2334:   PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
2335:   if (fem->ops->setfromoptions) {
2336:     (*fem->ops->setfromoptions)(fem);
2337:   }
2338:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
2339:   PetscObjectProcessOptionsHandlers((PetscObject) fem);
2340:   PetscOptionsEnd();
2341:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
2342:   return(0);
2343: }

2347: /*@C
2348:   PetscFESetUp - Construct data structures for the PetscFE

2350:   Collective on PetscFE

2352:   Input Parameter:
2353: . fem - the PetscFE object to setup

2355:   Level: developer

2357: .seealso PetscFEView(), PetscFEDestroy()
2358: @*/
2359: PetscErrorCode PetscFESetUp(PetscFE fem)
2360: {

2365:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
2366:   return(0);
2367: }

2371: /*@
2372:   PetscFEDestroy - Destroys a PetscFE object

2374:   Collective on PetscFE

2376:   Input Parameter:
2377: . fem - the PetscFE object to destroy

2379:   Level: developer

2381: .seealso PetscFEView()
2382: @*/
2383: PetscErrorCode PetscFEDestroy(PetscFE *fem)
2384: {

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

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

2394:   PetscFree((*fem)->numDof);
2395:   PetscFree((*fem)->invV);
2396:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
2397:   PetscSpaceDestroy(&(*fem)->basisSpace);
2398:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
2399:   PetscQuadratureDestroy(&(*fem)->quadrature);

2401:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
2402:   PetscHeaderDestroy(fem);
2403:   return(0);
2404: }

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

2411:   Collective on MPI_Comm

2413:   Input Parameter:
2414: . comm - The communicator for the PetscFE object

2416:   Output Parameter:
2417: . fem - The PetscFE object

2419:   Level: beginner

2421: .seealso: PetscFESetType(), PETSCFEGALERKIN
2422: @*/
2423: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
2424: {
2425:   PetscFE        f;

2430:   PetscCitationsRegister(FECitation,&FEcite);
2431:   *fem = NULL;
2432:   PetscFEInitializePackage();

2434:   PetscHeaderCreate(f, _p_PetscFE, struct _PetscFEOps, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
2435:   PetscMemzero(f->ops, sizeof(struct _PetscFEOps));

2437:   f->basisSpace    = NULL;
2438:   f->dualSpace     = NULL;
2439:   f->numComponents = 1;
2440:   f->numDof        = NULL;
2441:   f->invV          = NULL;
2442:   f->B             = NULL;
2443:   f->D             = NULL;
2444:   f->H             = NULL;
2445:   PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
2446:   f->blockSize     = 0;
2447:   f->numBlocks     = 1;
2448:   f->batchSize     = 0;
2449:   f->numBatches    = 1;

2451:   *fem = f;
2452:   return(0);
2453: }

2457: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
2458: {
2459:   DM             dm;

2465:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2466:   DMGetDimension(dm, dim);
2467:   return(0);
2468: }

2472: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
2473: {
2476:   fem->numComponents = comp;
2477:   return(0);
2478: }

2482: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
2483: {
2487:   *comp = fem->numComponents;
2488:   return(0);
2489: }

2493: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
2494: {
2497:   fem->blockSize  = blockSize;
2498:   fem->numBlocks  = numBlocks;
2499:   fem->batchSize  = batchSize;
2500:   fem->numBatches = numBatches;
2501:   return(0);
2502: }

2506: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
2507: {
2514:   if (blockSize)  *blockSize  = fem->blockSize;
2515:   if (numBlocks)  *numBlocks  = fem->numBlocks;
2516:   if (batchSize)  *batchSize  = fem->batchSize;
2517:   if (numBatches) *numBatches = fem->numBatches;
2518:   return(0);
2519: }

2523: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
2524: {
2528:   *sp = fem->basisSpace;
2529:   return(0);
2530: }

2534: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
2535: {

2541:   PetscSpaceDestroy(&fem->basisSpace);
2542:   fem->basisSpace = sp;
2543:   PetscObjectReference((PetscObject) fem->basisSpace);
2544:   return(0);
2545: }

2549: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
2550: {
2554:   *sp = fem->dualSpace;
2555:   return(0);
2556: }

2560: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
2561: {

2567:   PetscDualSpaceDestroy(&fem->dualSpace);
2568:   fem->dualSpace = sp;
2569:   PetscObjectReference((PetscObject) fem->dualSpace);
2570:   return(0);
2571: }

2575: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
2576: {
2580:   *q = fem->quadrature;
2581:   return(0);
2582: }

2586: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
2587: {

2592:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
2593:   PetscQuadratureDestroy(&fem->quadrature);
2594:   fem->quadrature = q;
2595:   PetscObjectReference((PetscObject) q);
2596:   return(0);
2597: }

2601: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
2602: {
2603:   const PetscInt *numDofDual;
2604:   PetscErrorCode  ierr;

2609:   PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
2610:   if (!fem->numDof) {
2611:     DM       dm;
2612:     PetscInt dim, d;

2614:     PetscDualSpaceGetDM(fem->dualSpace, &dm);
2615:     DMGetDimension(dm, &dim);
2616:     PetscMalloc1((dim+1), &fem->numDof);
2617:     for (d = 0; d <= dim; ++d) {
2618:       fem->numDof[d] = fem->numComponents*numDofDual[d];
2619:     }
2620:   }
2621:   *numDof = fem->numDof;
2622:   return(0);
2623: }

2627: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
2628: {
2629:   PetscInt         npoints = fem->quadrature->numPoints;
2630:   const PetscReal *points  = fem->quadrature->points;
2631:   PetscErrorCode   ierr;

2638:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
2639:   if (B) *B = fem->B;
2640:   if (D) *D = fem->D;
2641:   if (H) *H = fem->H;
2642:   return(0);
2643: }

2647: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
2648: {
2649:   DM               dm;
2650:   PetscInt         pdim; /* Dimension of FE space P */
2651:   PetscInt         dim;  /* Spatial dimension */
2652:   PetscInt         comp; /* Field components */
2653:   PetscErrorCode   ierr;

2661:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2662:   DMGetDimension(dm, &dim);
2663:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2664:   PetscFEGetNumComponents(fem, &comp);
2665:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
2666:   if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
2667:   if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
2668:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
2669:   return(0);
2670: }

2674: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
2675: {
2676:   DM             dm;

2681:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2682:   if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
2683:   if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
2684:   if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
2685:   return(0);
2686: }

2690: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
2691: {
2692:   PetscFE_Basic *b = (PetscFE_Basic *) fem->data;

2696:   PetscFree(b);
2697:   return(0);
2698: }

2702: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
2703: {
2704:   PetscSpace        basis;
2705:   PetscDualSpace    dual;
2706:   PetscQuadrature   q;
2707:   PetscInt          dim, Nq;
2708:   PetscViewerFormat format;
2709:   PetscErrorCode    ierr;

2712:   PetscFEGetBasisSpace(fe, &basis);
2713:   PetscFEGetDualSpace(fe, &dual);
2714:   PetscFEGetQuadrature(fe, &q);
2715:   PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
2716:   PetscViewerGetFormat(viewer, &format);
2717:   PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
2718:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2719:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
2720:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
2721:     PetscViewerASCIIPushTab(viewer);
2722:     PetscQuadratureView(q, viewer);
2723:     PetscViewerASCIIPopTab(viewer);
2724:   } else {
2725:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
2726:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
2727:   }
2728:   PetscViewerASCIIPushTab(viewer);
2729:   PetscSpaceView(basis, viewer);
2730:   PetscDualSpaceView(dual, viewer);
2731:   PetscViewerASCIIPopTab(viewer);
2732:   return(0);
2733: }

2737: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
2738: {
2739:   PetscBool      iascii;

2745:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2746:   if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
2747:   return(0);
2748: }

2752: /* Construct the change of basis from prime basis to nodal basis */
2753: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
2754: {
2755:   PetscReal     *work;
2756:   PetscBLASInt  *pivots;
2757: #ifndef PETSC_USE_COMPLEX
2758:   PetscBLASInt   n, info;
2759: #endif
2760:   PetscInt       pdim, j;

2764:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2765:   PetscMalloc1(pdim*pdim,&fem->invV);
2766:   for (j = 0; j < pdim; ++j) {
2767:     PetscReal      *Bf;
2768:     PetscQuadrature f;
2769:     PetscInt        q, k;

2771:     PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
2772:     PetscMalloc1(f->numPoints*pdim,&Bf);
2773:     PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
2774:     for (k = 0; k < pdim; ++k) {
2775:       /* n_j \cdot \phi_k */
2776:       fem->invV[j*pdim+k] = 0.0;
2777:       for (q = 0; q < f->numPoints; ++q) {
2778:         fem->invV[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
2779:       }
2780:     }
2781:     PetscFree(Bf);
2782:   }
2783:   PetscMalloc2(pdim,&pivots,pdim,&work);
2784: #ifndef PETSC_USE_COMPLEX
2785:   n = pdim;
2786:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, fem->invV, &n, pivots, &info));
2787:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
2788: #endif
2789:   PetscFree2(pivots,work);
2790:   return(0);
2791: }

2795: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
2796: {

2800:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
2801:   return(0);
2802: }

2806: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
2807: {
2808:   DM               dm;
2809:   PetscInt         pdim; /* Dimension of FE space P */
2810:   PetscInt         dim;  /* Spatial dimension */
2811:   PetscInt         comp; /* Field components */
2812:   PetscReal       *tmpB, *tmpD, *tmpH;
2813:   PetscInt         p, d, j, k;
2814:   PetscErrorCode   ierr;

2817:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2818:   DMGetDimension(dm, &dim);
2819:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2820:   PetscFEGetNumComponents(fem, &comp);
2821:   /* Evaluate the prime basis functions at all points */
2822:   if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
2823:   if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
2824:   if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
2825:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
2826:   /* Translate to the nodal basis */
2827:   for (p = 0; p < npoints; ++p) {
2828:     if (B) {
2829:       /* Multiply by V^{-1} (pdim x pdim) */
2830:       for (j = 0; j < pdim; ++j) {
2831:         const PetscInt i = (p*pdim + j)*comp;
2832:         PetscInt       c;

2834:         B[i] = 0.0;
2835:         for (k = 0; k < pdim; ++k) {
2836:           B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
2837:         }
2838:         for (c = 1; c < comp; ++c) {
2839:           B[i+c] = B[i];
2840:         }
2841:       }
2842:     }
2843:     if (D) {
2844:       /* Multiply by V^{-1} (pdim x pdim) */
2845:       for (j = 0; j < pdim; ++j) {
2846:         for (d = 0; d < dim; ++d) {
2847:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
2848:           PetscInt       c;

2850:           D[i] = 0.0;
2851:           for (k = 0; k < pdim; ++k) {
2852:             D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
2853:           }
2854:           for (c = 1; c < comp; ++c) {
2855:             D[((p*pdim + j)*comp + c)*dim + d] = D[i];
2856:           }
2857:         }
2858:       }
2859:     }
2860:     if (H) {
2861:       /* Multiply by V^{-1} (pdim x pdim) */
2862:       for (j = 0; j < pdim; ++j) {
2863:         for (d = 0; d < dim*dim; ++d) {
2864:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
2865:           PetscInt       c;

2867:           H[i] = 0.0;
2868:           for (k = 0; k < pdim; ++k) {
2869:             H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
2870:           }
2871:           for (c = 1; c < comp; ++c) {
2872:             H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
2873:           }
2874:         }
2875:       }
2876:     }
2877:   }
2878:   if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
2879:   if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
2880:   if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
2881:   return(0);
2882: }

2886: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
2887:                                       const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
2888: {
2889:   const PetscInt  debug = 0;
2890:   void          (*obj_func)(const PetscScalar u[], const PetscScalar u_x[], const PetscScalar u_t[], const PetscScalar a[], const PetscScalar a_x[], const PetscScalar a_t[], const PetscReal x[], PetscScalar obj[]);
2891:   PetscQuadrature quad;
2892:   PetscScalar    *u, *u_x, *a, *a_x;
2893:   PetscReal      *x;
2894:   PetscInt        dim, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, e;
2895:   PetscErrorCode  ierr;

2898:   PetscDSGetObjective(prob, field, &obj_func);
2899:   if (!obj_func) return(0);
2900:   PetscFEGetSpatialDimension(fem, &dim);
2901:   PetscFEGetQuadrature(fem, &quad);
2902:   PetscDSGetTotalDimension(prob, &totDim);
2903:   PetscDSGetTotalDimension(probAux, &totDimAux);
2904:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
2905:   PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
2906:   PetscDSGetRefCoordArrays(prob, &x, NULL);
2907:   for (e = 0; e < Ne; ++e) {
2908:     const PetscReal *v0   = &geom.v0[e*dim];
2909:     const PetscReal *J    = &geom.J[e*dim*dim];
2910:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
2911:     const PetscReal  detJ = geom.detJ[e];
2912:     const PetscReal *quadPoints, *quadWeights;
2913:     PetscInt         Nq, q;

2915:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
2916:     if (debug > 1) {
2917:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
2918: #ifndef PETSC_USE_COMPLEX
2919:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
2920: #endif
2921:     }
2922:     for (q = 0; q < Nq; ++q) {
2923:       PetscScalar integrand;

2925:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
2926:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
2927:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       NULL, u, u_x, NULL);
2928:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
2929:       obj_func(u, NULL, u_x, a, NULL, a_x, x, &integrand);
2930:       integrand *= detJ*quadWeights[q];
2931:       integral[field] += PetscRealPart(integrand);
2932:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", PetscRealPart(integrand), integral[field]);}
2933:     }
2934:     cOffset    += totDim;
2935:     cOffsetAux += totDimAux;
2936:   }
2937:   return(0);
2938: }

2942: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
2943:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
2944: {
2945:   const PetscInt  debug = 0;
2946:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
2947:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
2948:   PetscQuadrature quad;
2949:   PetscReal     **basisField, **basisFieldDer;
2950:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
2951:   PetscReal      *x;
2952:   PetscInt        dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
2953:   PetscErrorCode  ierr;

2956:   PetscFEGetSpatialDimension(fem, &dim);
2957:   PetscFEGetQuadrature(fem, &quad);
2958:   PetscFEGetDimension(fem, &Nb);
2959:   PetscFEGetNumComponents(fem, &Nc);
2960:   PetscDSGetTotalDimension(prob, &totDim);
2961:   PetscDSGetFieldOffset(prob, field, &fOffset);
2962:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
2963:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
2964:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
2965:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
2966:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
2967:   if (probAux) {
2968:     PetscDSGetTotalDimension(probAux, &totDimAux);
2969:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
2970:   }
2971:   for (e = 0; e < Ne; ++e) {
2972:     const PetscReal *v0   = &geom.v0[e*dim];
2973:     const PetscReal *J    = &geom.J[e*dim*dim];
2974:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
2975:     const PetscReal  detJ = geom.detJ[e];
2976:     const PetscReal *quadPoints, *quadWeights;
2977:     PetscInt         Nq, q;

2979:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
2980:     if (debug > 1) {
2981:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
2982: #ifndef PETSC_USE_COMPLEX
2983:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
2984: #endif
2985:     }
2986:     for (q = 0; q < Nq; ++q) {
2987:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
2988:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
2989:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
2990:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
2991:       f0_func(u, u_t, u_x, a, NULL, a_x, x, &f0[q*Nc]);
2992:       f1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
2993:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
2994:     }
2995:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
2996:     cOffset    += totDim;
2997:     cOffsetAux += totDimAux;
2998:   }
2999:   return(0);
3000: }

3004: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3005:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3006: {
3007:   const PetscInt  debug = 0;
3008:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]);
3009:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]);
3010:   PetscQuadrature quad;
3011:   PetscReal     **basisField, **basisFieldDer;
3012:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3013:   PetscReal      *x;
3014:   PetscInt        dim, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3015:   PetscErrorCode  ierr;

3018:   PetscFEGetSpatialDimension(fem, &dim);
3019:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3020:   PetscFEGetQuadrature(fem, &quad);
3021:   PetscFEGetDimension(fem, &Nb);
3022:   PetscFEGetNumComponents(fem, &Nc);
3023:   PetscDSGetTotalBdDimension(prob, &totDim);
3024:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3025:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3026:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3027:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3028:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3029:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3030:   if (probAux) {
3031:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3032:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3033:   }
3034:   for (e = 0; e < Ne; ++e) {
3035:     const PetscReal *v0          = &geom.v0[e*dim];
3036:     const PetscReal *n           = &geom.n[e*dim];
3037:     const PetscReal *J           = &geom.J[e*dim*dim];
3038:     const PetscReal *invJ        = &geom.invJ[e*dim*dim];
3039:     const PetscReal  detJ        = geom.detJ[e];
3040:     const PetscReal *quadPoints, *quadWeights;
3041:     PetscInt         Nq, q;

3043:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3044:     if (debug > 1) {
3045:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3046: #ifndef PETSC_USE_COMPLEX
3047:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3048: #endif
3049:     }
3050:     for (q = 0; q < Nq; ++q) {
3051:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3052:       CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3053:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3054:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3055:       f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3056:       f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3057:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3058:     }
3059:     UpdateElementVec(dim-1, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3060:     cOffset    += totDim;
3061:     cOffsetAux += totDimAux;
3062:   }
3063:   return(0);
3064: }

3068: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3069:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3070: {
3071:   const PetscInt  debug      = 0;
3072:   void          (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g0[]);
3073:   void          (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[]);
3074:   void          (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[]);
3075:   void          (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[]);
3076:   PetscFE         fe;
3077:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3078:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3079:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3080:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3081:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3082:   PetscQuadrature quad;
3083:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3084:   PetscReal      *x;
3085:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3086:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3087:   PetscInt        dim, totDim, totDimAux, e;
3088:   PetscErrorCode  ierr;

3091:   PetscFEGetSpatialDimension(fem, &dim);
3092:   PetscFEGetQuadrature(fem, &quad);
3093:   PetscDSGetTotalDimension(prob, &totDim);
3094:   PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3095:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3096:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3097:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3098:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3099:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3100:   PetscFEGetDimension(fe, &NbI);
3101:   PetscFEGetNumComponents(fe, &NcI);
3102:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3103:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3104:   PetscFEGetDimension(fe, &NbJ);
3105:   PetscFEGetNumComponents(fe, &NcJ);
3106:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3107:   if (probAux) {
3108:     PetscDSGetTotalDimension(probAux, &totDimAux);
3109:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3110:   }
3111:   basisI    = basisField[fieldI];
3112:   basisJ    = basisField[fieldJ];
3113:   basisDerI = basisFieldDer[fieldI];
3114:   basisDerJ = basisFieldDer[fieldJ];
3115:   /* Initialize here in case the function is not defined */
3116:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3117:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3118:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3119:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3120:   for (e = 0; e < Ne; ++e) {
3121:     const PetscReal *v0   = &geom.v0[e*dim];
3122:     const PetscReal *J    = &geom.J[e*dim*dim];
3123:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
3124:     const PetscReal  detJ = geom.detJ[e];
3125:     const PetscReal *quadPoints, *quadWeights;
3126:     PetscInt         Nq, q;

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

3132:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3133:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3134:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3135:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3136:       if (g0_func) {
3137:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3138:         g0_func(u, u_t, u_x, a, NULL, a_x, x, g0);
3139:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3140:       }
3141:       if (g1_func) {
3142:         PetscInt d, d2;
3143:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3144:         g1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3145:         for (fc = 0; fc < NcI; ++fc) {
3146:           for (gc = 0; gc < NcJ; ++gc) {
3147:             for (d = 0; d < dim; ++d) {
3148:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3149:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3150:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3151:             }
3152:           }
3153:         }
3154:       }
3155:       if (g2_func) {
3156:         PetscInt d, d2;
3157:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3158:         g2_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3159:         for (fc = 0; fc < NcI; ++fc) {
3160:           for (gc = 0; gc < NcJ; ++gc) {
3161:             for (d = 0; d < dim; ++d) {
3162:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3163:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3164:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3165:             }
3166:           }
3167:         }
3168:       }
3169:       if (g3_func) {
3170:         PetscInt d, d2, dp, d3;
3171:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3172:         g3_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3173:         for (fc = 0; fc < NcI; ++fc) {
3174:           for (gc = 0; gc < NcJ; ++gc) {
3175:             for (d = 0; d < dim; ++d) {
3176:               for (dp = 0; dp < dim; ++dp) {
3177:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3178:                 for (d2 = 0; d2 < dim; ++d2) {
3179:                   for (d3 = 0; d3 < dim; ++d3) {
3180:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3181:                   }
3182:                 }
3183:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3184:               }
3185:             }
3186:           }
3187:         }
3188:       }

3190:       for (f = 0; f < NbI; ++f) {
3191:         for (fc = 0; fc < NcI; ++fc) {
3192:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3193:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3194:           for (g = 0; g < NbJ; ++g) {
3195:             for (gc = 0; gc < NcJ; ++gc) {
3196:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3197:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3198:               const PetscInt fOff = eOffset+i*totDim+j;
3199:               PetscInt       d, d2;

3201:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3202:               for (d = 0; d < dim; ++d) {
3203:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3204:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3205:                 for (d2 = 0; d2 < dim; ++d2) {
3206:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3207:                 }
3208:               }
3209:             }
3210:           }
3211:         }
3212:       }
3213:     }
3214:     if (debug > 1) {
3215:       PetscInt fc, f, gc, g;

3217:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3218:       for (fc = 0; fc < NcI; ++fc) {
3219:         for (f = 0; f < NbI; ++f) {
3220:           const PetscInt i = offsetI + f*NcI+fc;
3221:           for (gc = 0; gc < NcJ; ++gc) {
3222:             for (g = 0; g < NbJ; ++g) {
3223:               const PetscInt j = offsetJ + g*NcJ+gc;
3224:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3225:             }
3226:           }
3227:           PetscPrintf(PETSC_COMM_SELF, "\n");
3228:         }
3229:       }
3230:     }
3231:     cOffset    += totDim;
3232:     cOffsetAux += totDimAux;
3233:     eOffset    += PetscSqr(totDim);
3234:   }
3235:   return(0);
3236: }

3240: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3241:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3242: {
3243:   const PetscInt  debug      = 0;
3244:   void          (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g0[]);
3245:   void          (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g1[]);
3246:   void          (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g2[]);
3247:   void          (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g3[]);
3248:   PetscFE         fe;
3249:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3250:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3251:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3252:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3253:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3254:   PetscQuadrature quad;
3255:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3256:   PetscReal      *x;
3257:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3258:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3259:   PetscInt        dim, bdim, totDim, totDimAux, e;
3260:   PetscErrorCode  ierr;

3263:   PetscFEGetQuadrature(fem, &quad);
3264:   PetscFEGetSpatialDimension(fem, &dim);
3265:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3266:   bdim = dim-1;
3267:   PetscDSGetTotalBdDimension(prob, &totDim);
3268:   PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3269:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3270:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3271:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3272:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3273:   PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
3274:   PetscFEGetDimension(fe, &NbI);
3275:   PetscFEGetNumComponents(fe, &NcI);
3276:   PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
3277:   PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
3278:   PetscFEGetDimension(fe, &NbJ);
3279:   PetscFEGetNumComponents(fe, &NcJ);
3280:   PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
3281:   if (probAux) {
3282:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3283:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3284:   }
3285:   basisI    = basisField[fieldI];
3286:   basisJ    = basisField[fieldJ];
3287:   basisDerI = basisFieldDer[fieldI];
3288:   basisDerJ = basisFieldDer[fieldJ];
3289:   /* Initialize here in case the function is not defined */
3290:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3291:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3292:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3293:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3294:   for (e = 0; e < Ne; ++e) {
3295:     const PetscReal *v0   = &geom.v0[e*dim];
3296:     const PetscReal *n    = &geom.n[e*dim];
3297:     const PetscReal *J    = &geom.J[e*dim*dim];
3298:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
3299:     const PetscReal  detJ = geom.detJ[e];
3300:     const PetscReal *quadPoints, *quadWeights;
3301:     PetscInt         Nq, q;

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

3307:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3308:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3309:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3310:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3311:       /* TODO: I think I have a mistake in the dim loops here */
3312:       if (g0_func) {
3313:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3314:         g0_func(u, u_t, u_x, a, NULL, a_x, x, n, g0);
3315:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3316:       }
3317:       if (g1_func) {
3318:         PetscInt d, d2;
3319:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3320:         g1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3321:         for (fc = 0; fc < NcI; ++fc) {
3322:           for (gc = 0; gc < NcJ; ++gc) {
3323:             for (d = 0; d < dim; ++d) {
3324:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3325:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3326:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3327:             }
3328:           }
3329:         }
3330:       }
3331:       if (g2_func) {
3332:         PetscInt d, d2;
3333:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3334:         g2_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3335:         for (fc = 0; fc < NcI; ++fc) {
3336:           for (gc = 0; gc < NcJ; ++gc) {
3337:             for (d = 0; d < dim; ++d) {
3338:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3339:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3340:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3341:             }
3342:           }
3343:         }
3344:       }
3345:       if (g3_func) {
3346:         PetscInt d, d2, dp, d3;
3347:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3348:         g3_func(u, u_t, u_x, a, NULL, a_x, x, n, g3);
3349:         for (fc = 0; fc < NcI; ++fc) {
3350:           for (gc = 0; gc < NcJ; ++gc) {
3351:             for (d = 0; d < dim; ++d) {
3352:               for (dp = 0; dp < dim; ++dp) {
3353:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3354:                 for (d2 = 0; d2 < dim; ++d2) {
3355:                   for (d3 = 0; d3 < dim; ++d3) {
3356:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3357:                   }
3358:                 }
3359:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3360:               }
3361:             }
3362:           }
3363:         }
3364:       }

3366:       for (f = 0; f < NbI; ++f) {
3367:         for (fc = 0; fc < NcI; ++fc) {
3368:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3369:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3370:           for (g = 0; g < NbJ; ++g) {
3371:             for (gc = 0; gc < NcJ; ++gc) {
3372:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3373:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3374:               const PetscInt fOff = eOffset+i*totDim+j;
3375:               PetscInt       d, d2;

3377:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3378:               for (d = 0; d < bdim; ++d) {
3379:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
3380:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
3381:                 for (d2 = 0; d2 < bdim; ++d2) {
3382:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
3383:                 }
3384:               }
3385:             }
3386:           }
3387:         }
3388:       }
3389:     }
3390:     if (debug > 1) {
3391:       PetscInt fc, f, gc, g;

3393:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3394:       for (fc = 0; fc < NcI; ++fc) {
3395:         for (f = 0; f < NbI; ++f) {
3396:           const PetscInt i = offsetI + f*NcI+fc;
3397:           for (gc = 0; gc < NcJ; ++gc) {
3398:             for (g = 0; g < NbJ; ++g) {
3399:               const PetscInt j = offsetJ + g*NcJ+gc;
3400:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3401:             }
3402:           }
3403:           PetscPrintf(PETSC_COMM_SELF, "\n");
3404:         }
3405:       }
3406:     }
3407:     cOffset    += totDim;
3408:     cOffsetAux += totDimAux;
3409:     eOffset    += PetscSqr(totDim);
3410:   }
3411:   return(0);
3412: }

3416: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
3417: {
3419:   fem->ops->setfromoptions          = NULL;
3420:   fem->ops->setup                   = PetscFESetUp_Basic;
3421:   fem->ops->view                    = PetscFEView_Basic;
3422:   fem->ops->destroy                 = PetscFEDestroy_Basic;
3423:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
3424:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
3425:   fem->ops->integrate               = PetscFEIntegrate_Basic;
3426:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
3427:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
3428:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
3429:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
3430:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
3431:   return(0);
3432: }

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

3437:   Level: intermediate

3439: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3440: M*/

3444: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
3445: {
3446:   PetscFE_Basic *b;

3451:   PetscNewLog(fem,&b);
3452:   fem->data = b;

3454:   PetscFEInitialize_Basic(fem);
3455:   return(0);
3456: }

3460: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
3461: {
3462:   PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;

3466:   PetscFree(na);
3467:   return(0);
3468: }

3472: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3473:                                                   const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3474: {
3475:   const PetscInt  debug = 0;
3476:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
3477:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
3478:   PetscQuadrature quad;
3479:   PetscReal     **basisField, **basisFieldDer;
3480:   PetscScalar    *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3481:   PetscReal      *x;
3482:   PetscInt        dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3483:   PetscErrorCode  ierr;

3486:   PetscFEGetSpatialDimension(fem, &dim);
3487:   PetscFEGetQuadrature(fem, &quad);
3488:   PetscFEGetDimension(fem, &Nb);
3489:   PetscFEGetNumComponents(fem, &Nc);
3490:   PetscDSGetTotalDimension(prob, &totDim);
3491:   PetscDSGetFieldOffset(prob, field, &fOffset);
3492:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3493:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3494:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3495:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3496:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3497:   if (probAux) {
3498:     PetscDSGetTotalDimension(probAux, &totDimAux);
3499:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3500:   }
3501:   for (e = 0; e < Ne; ++e) {
3502:     const PetscReal *quadPoints, *quadWeights;
3503:     PetscInt         Nq, q;

3505:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3506:     for (q = 0; q < Nq; ++q) {
3507:       const PetscReal *v0   = &geom.v0[(e*Nq+q)*dim];
3508:       const PetscReal *J    = &geom.J[(e*Nq+q)*dim*dim];
3509:       const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3510:       const PetscReal  detJ = geom.detJ[e*Nq+q];

3512:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3513:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3514:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3515:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3516:       f0_func(u, u_t, u_x, a, NULL, a_x, x, &f0[q*Nc]);
3517:       f1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3518:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3519:     }
3520:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3521:     cOffset    += totDim;
3522:     cOffsetAux += totDimAux;
3523:   }
3524:   return(0);
3525: }

3529: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3530:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3531: {
3532:   const PetscInt  debug = 0;
3533:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]);
3534:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]);
3535:   PetscQuadrature quad;
3536:   PetscReal    **basisField, **basisFieldDer;
3537:   PetscScalar    *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3538:   PetscReal      *x;
3539:   PetscInt        dim, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3540:   PetscErrorCode  ierr;

3543:   PetscFEGetSpatialDimension(fem, &dim);
3544:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3545:   PetscFEGetQuadrature(fem, &quad);
3546:   PetscFEGetDimension(fem, &Nb);
3547:   PetscFEGetNumComponents(fem, &Nc);
3548:   PetscDSGetTotalBdDimension(prob, &totDim);
3549:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3550:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3551:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3552:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3553:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3554:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3555:   if (probAux) {
3556:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3557:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3558:   }
3559:   for (e = 0; e < Ne; ++e) {
3560:     const PetscReal *quadPoints, *quadWeights;
3561:     PetscInt         Nq, q;

3563:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3564:     for (q = 0; q < Nq; ++q) {
3565:       const PetscReal *v0   = &geom.v0[(e*Nq+q)*dim];
3566:       const PetscReal *n    = &geom.n[(e*Nq+q)*dim];
3567:       const PetscReal *J    = &geom.J[(e*Nq+q)*dim*dim];
3568:       const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3569:       const PetscReal  detJ = geom.detJ[e*Nq+q];

3571:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3572:       CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3573:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3574:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3575:       f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3576:       f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3577:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3578:     }
3579:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3580:     cOffset    += totDim;
3581:     cOffsetAux += totDimAux;
3582:   }
3583:   return(0);
3584: }

3588: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3589:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3590: {
3591:   const PetscInt  debug      = 0;
3592:   void          (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g0[]);
3593:   void          (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[]);
3594:   void          (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[]);
3595:   void          (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[]);
3596:   PetscFE         fe;
3597:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3598:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3599:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3600:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3601:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3602:   PetscQuadrature quad;
3603:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3604:   PetscReal      *x;
3605:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3606:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3607:   PetscInt        dim, totDim, totDimAux, e;
3608:   PetscErrorCode  ierr;

3611:   PetscFEGetSpatialDimension(fem, &dim);
3612:   PetscFEGetQuadrature(fem, &quad);
3613:   PetscDSGetTotalDimension(prob, &totDim);
3614:   PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3615:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3616:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3617:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3618:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3619:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3620:   PetscFEGetDimension(fe, &NbI);
3621:   PetscFEGetNumComponents(fe, &NcI);
3622:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3623:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3624:   PetscFEGetDimension(fe, &NbJ);
3625:   PetscFEGetNumComponents(fe, &NcJ);
3626:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3627:   if (probAux) {
3628:     PetscDSGetTotalDimension(probAux, &totDimAux);
3629:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3630:   }
3631:   basisI    = basisField[fieldI];
3632:   basisJ    = basisField[fieldJ];
3633:   basisDerI = basisFieldDer[fieldI];
3634:   basisDerJ = basisFieldDer[fieldJ];
3635:   /* Initialize here in case the function is not defined */
3636:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3637:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3638:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3639:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3640:   for (e = 0; e < Ne; ++e) {
3641:     const PetscReal *quadPoints, *quadWeights;
3642:     PetscInt         Nq, q;

3644:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3645:     for (q = 0; q < Nq; ++q) {
3646:       const PetscReal *v0   = &geom.v0[(e*Nq+q)*dim];
3647:       const PetscReal *J    = &geom.J[(e*Nq+q)*dim*dim];
3648:       const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3649:       const PetscReal  detJ = geom.detJ[e*Nq+q];
3650:       PetscInt         f, g, fc, gc, c;

3652:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3653:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3654:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3655:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3656:       if (g0_func) {
3657:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3658:         g0_func(u, u_t, u_x, a, NULL, a_x, x, g0);
3659:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3660:       }
3661:       if (g1_func) {
3662:         PetscInt d, d2;
3663:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3664:         g1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3665:         for (fc = 0; fc < NcI; ++fc) {
3666:           for (gc = 0; gc < NcJ; ++gc) {
3667:             for (d = 0; d < dim; ++d) {
3668:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3669:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3670:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3671:             }
3672:           }
3673:         }
3674:       }
3675:       if (g2_func) {
3676:         PetscInt d, d2;
3677:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3678:         g2_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3679:         for (fc = 0; fc < NcI; ++fc) {
3680:           for (gc = 0; gc < NcJ; ++gc) {
3681:             for (d = 0; d < dim; ++d) {
3682:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3683:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3684:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3685:             }
3686:           }
3687:         }
3688:       }
3689:       if (g3_func) {
3690:         PetscInt d, d2, dp, d3;
3691:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3692:         g3_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3693:         for (fc = 0; fc < NcI; ++fc) {
3694:           for (gc = 0; gc < NcJ; ++gc) {
3695:             for (d = 0; d < dim; ++d) {
3696:               for (dp = 0; dp < dim; ++dp) {
3697:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3698:                 for (d2 = 0; d2 < dim; ++d2) {
3699:                   for (d3 = 0; d3 < dim; ++d3) {
3700:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3701:                   }
3702:                 }
3703:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3704:               }
3705:             }
3706:           }
3707:         }
3708:       }

3710:       for (f = 0; f < NbI; ++f) {
3711:         for (fc = 0; fc < NcI; ++fc) {
3712:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3713:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3714:           for (g = 0; g < NbJ; ++g) {
3715:             for (gc = 0; gc < NcJ; ++gc) {
3716:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3717:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3718:               const PetscInt fOff = eOffset+i*totDim+j;
3719:               PetscInt       d, d2;

3721:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3722:               for (d = 0; d < dim; ++d) {
3723:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3724:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3725:                 for (d2 = 0; d2 < dim; ++d2) {
3726:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3727:                 }
3728:               }
3729:             }
3730:           }
3731:         }
3732:       }
3733:     }
3734:     if (debug > 1) {
3735:       PetscInt fc, f, gc, g;

3737:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3738:       for (fc = 0; fc < NcI; ++fc) {
3739:         for (f = 0; f < NbI; ++f) {
3740:           const PetscInt i = offsetI + f*NcI+fc;
3741:           for (gc = 0; gc < NcJ; ++gc) {
3742:             for (g = 0; g < NbJ; ++g) {
3743:               const PetscInt j = offsetJ + g*NcJ+gc;
3744:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3745:             }
3746:           }
3747:           PetscPrintf(PETSC_COMM_SELF, "\n");
3748:         }
3749:       }
3750:     }
3751:     cOffset    += totDim;
3752:     cOffsetAux += totDimAux;
3753:     eOffset    += PetscSqr(totDim);
3754:   }
3755:   return(0);
3756: }

3760: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
3761: {
3763:   fem->ops->setfromoptions          = NULL;
3764:   fem->ops->setup                   = PetscFESetUp_Basic;
3765:   fem->ops->view                    = NULL;
3766:   fem->ops->destroy                 = PetscFEDestroy_Nonaffine;
3767:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
3768:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
3769:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Nonaffine;
3770:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Nonaffine;
3771:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
3772:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Nonaffine;
3773:   return(0);
3774: }

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

3779:   Level: intermediate

3781: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3782: M*/

3786: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
3787: {
3788:   PetscFE_Nonaffine *na;
3789:   PetscErrorCode     ierr;

3793:   PetscNewLog(fem, &na);
3794:   fem->data = na;

3796:   PetscFEInitialize_Nonaffine(fem);
3797:   return(0);
3798: }

3800: #ifdef PETSC_HAVE_OPENCL

3804: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
3805: {
3806:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
3807:   PetscErrorCode  ierr;

3810:   clReleaseCommandQueue(ocl->queue_id);
3811:   ocl->queue_id = 0;
3812:   clReleaseContext(ocl->ctx_id);
3813:   ocl->ctx_id = 0;
3814:   PetscFree(ocl);
3815:   return(0);
3816: }

3818: #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)
3819: enum {LAPLACIAN = 0, ELASTICITY = 1};

3823: /* dim     Number of spatial dimensions:          2                   */
3824: /* N_b     Number of basis functions:             generated           */
3825: /* N_{bt}  Number of total basis functions:       N_b * N_{comp}      */
3826: /* N_q     Number of quadrature points:           generated           */
3827: /* N_{bs}  Number of block cells                  LCM(N_b, N_q)       */
3828: /* N_{bst} Number of block cell components        LCM(N_{bt}, N_q)    */
3829: /* N_{bl}  Number of concurrent blocks            generated           */
3830: /* N_t     Number of threads:                     N_{bl} * N_{bs}     */
3831: /* N_{cbc} Number of concurrent basis      cells: N_{bl} * N_q        */
3832: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b        */
3833: /* N_{sbc} Number of serial     basis      cells: N_{bs} / N_q        */
3834: /* N_{sqc} Number of serial     quadrature cells: N_{bs} / N_b        */
3835: /* N_{cb}  Number of serial cell batches:         input               */
3836: /* N_c     Number of total cells:                 N_{cb}*N_{t}/N_{comp} */
3837: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
3838: {
3839:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
3840:   PetscQuadrature q;
3841:   char           *string_tail   = *string_buffer;
3842:   char           *end_of_buffer = *string_buffer + buffer_length;
3843:   char            float_str[]   = "float", double_str[]  = "double";
3844:   char           *numeric_str   = &(float_str[0]);
3845:   PetscInt        op            = ocl->op;
3846:   PetscBool       useField      = PETSC_FALSE;
3847:   PetscBool       useFieldDer   = PETSC_TRUE;
3848:   PetscBool       useFieldAux   = useAux;
3849:   PetscBool       useFieldDerAux= PETSC_FALSE;
3850:   PetscBool       useF0         = PETSC_TRUE;
3851:   PetscBool       useF1         = PETSC_TRUE;
3852:   PetscReal      *basis, *basisDer;
3853:   PetscInt        dim, N_b, N_c, N_q, N_t, p, d, b, c;
3854:   size_t          count;
3855:   PetscErrorCode  ierr;

3858:   PetscFEGetSpatialDimension(fem, &dim);
3859:   PetscFEGetDimension(fem, &N_b);
3860:   PetscFEGetNumComponents(fem, &N_c);
3861:   PetscFEGetQuadrature(fem, &q);
3862:   N_q  = q->numPoints;
3863:   N_t  = N_b * N_c * N_q * N_bl;
3864:   /* Enable device extension for double precision */
3865:   if (ocl->realType == PETSC_DOUBLE) {
3866:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3867: "#if defined(cl_khr_fp64)\n"
3868: "#  pragma OPENCL EXTENSION cl_khr_fp64: enable\n"
3869: "#elif defined(cl_amd_fp64)\n"
3870: "#  pragma OPENCL EXTENSION cl_amd_fp64: enable\n"
3871: "#endif\n",
3872:                               &count);STRING_ERROR_CHECK("Message to short");
3873:     numeric_str  = &(double_str[0]);
3874:   }
3875:   /* Kernel API */
3876:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3877: "\n"
3878: "__kernel void integrateElementQuadrature(int N_cb, __global %s *coefficients, __global %s *coefficientsAux, __global %s *jacobianInverses, __global %s *jacobianDeterminants, __global %s *elemVec)\n"
3879: "{\n",
3880:                        &count, numeric_str, numeric_str, numeric_str, numeric_str, numeric_str);STRING_ERROR_CHECK("Message to short");
3881:   /* Quadrature */
3882:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3883: "  /* Quadrature points\n"
3884: "   - (x1,y1,x2,y2,...) */\n"
3885: "  const %s points[%d] = {\n",
3886:                        &count, numeric_str, N_q*dim);STRING_ERROR_CHECK("Message to short");
3887:   for (p = 0; p < N_q; ++p) {
3888:     for (d = 0; d < dim; ++d) {
3889:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->points[p*dim+d]);STRING_ERROR_CHECK("Message to short");
3890:     }
3891:   }
3892:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3893:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3894: "  /* Quadrature weights\n"
3895: "   - (v1,v2,...) */\n"
3896: "  const %s weights[%d] = {\n",
3897:                        &count, numeric_str, N_q);STRING_ERROR_CHECK("Message to short");
3898:   for (p = 0; p < N_q; ++p) {
3899:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->weights[p]);STRING_ERROR_CHECK("Message to short");
3900:   }
3901:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3902:   /* Basis Functions */
3903:   PetscFEGetDefaultTabulation(fem, &basis, &basisDer, NULL);
3904:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3905: "  /* Nodal basis function evaluations\n"
3906: "    - basis component is fastest varying, the basis function, then point */\n"
3907: "  const %s Basis[%d] = {\n",
3908:                        &count, numeric_str, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
3909:   for (p = 0; p < N_q; ++p) {
3910:     for (b = 0; b < N_b; ++b) {
3911:       for (c = 0; c < N_c; ++c) {
3912:         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");
3913:       }
3914:     }
3915:   }
3916:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3917:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3918: "\n"
3919: "  /* Nodal basis function derivative evaluations,\n"
3920: "      - derivative direction is fastest varying, then basis component, then basis function, then point */\n"
3921: "  const %s%d BasisDerivatives[%d] = {\n",
3922:                        &count, numeric_str, dim, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
3923:   for (p = 0; p < N_q; ++p) {
3924:     for (b = 0; b < N_b; ++b) {
3925:       for (c = 0; c < N_c; ++c) {
3926:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "(%s%d)(", &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
3927:         for (d = 0; d < dim; ++d) {
3928:           if (d > 0) {
3929:             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");
3930:           } else {
3931:             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");
3932:           }
3933:         }
3934:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "),\n", &count);STRING_ERROR_CHECK("Message to short");
3935:       }
3936:     }
3937:   }
3938:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3939:   /* Sizes */
3940:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3941: "  const int dim    = %d;                           // The spatial dimension\n"
3942: "  const int N_bl   = %d;                           // The number of concurrent blocks\n"
3943: "  const int N_b    = %d;                           // The number of basis functions\n"
3944: "  const int N_comp = %d;                           // The number of basis function components\n"
3945: "  const int N_bt   = N_b*N_comp;                    // The total number of scalar basis functions\n"
3946: "  const int N_q    = %d;                           // The number of quadrature points\n"
3947: "  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"
3948: "  const int N_t    = N_bst*N_bl;                    // The number of threads, N_bst * N_bl\n"
3949: "  const int N_bc   = N_t/N_comp;                    // The number of cells per batch (N_b*N_q*N_bl)\n"
3950: "  const int N_sbc  = N_bst / (N_q * N_comp);\n"
3951: "  const int N_sqc  = N_bst / N_bt;\n"
3952: "  /*const int N_c    = N_cb * N_bc;*/\n"
3953: "\n"
3954: "  /* Calculated indices */\n"
3955: "  /*const int tidx    = get_local_id(0) + get_local_size(0)*get_local_id(1);*/\n"
3956: "  const int tidx    = get_local_id(0);\n"
3957: "  const int blidx   = tidx / N_bst;                  // Block number for this thread\n"
3958: "  const int bidx    = tidx %% N_bt;                   // Basis function mapped to this thread\n"
3959: "  const int cidx    = tidx %% N_comp;                 // Basis component mapped to this thread\n"
3960: "  const int qidx    = tidx %% N_q;                    // Quadrature point mapped to this thread\n"
3961: "  const int blbidx  = tidx %% N_q + blidx*N_q;        // Cell mapped to this thread in the basis phase\n"
3962: "  const int blqidx  = tidx %% N_b + blidx*N_b;        // Cell mapped to this thread in the quadrature phase\n"
3963: "  const int gidx    = get_group_id(1)*get_num_groups(0) + get_group_id(0);\n"
3964: "  const int Goffset = gidx*N_cb*N_bc;\n",
3965:                             &count, dim, N_bl, N_b, N_c, N_q);STRING_ERROR_CHECK("Message to short");
3966:   /* Local memory */
3967:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3968: "\n"
3969: "  /* Quadrature data */\n"
3970: "  %s                w;                   // $w_q$, Quadrature weight at $x_q$\n"
3971: "  __local %s         phi_i[%d];    //[N_bt*N_q];  // $\\phi_i(x_q)$, Value of the basis function $i$ at $x_q$\n"
3972: "  __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"
3973: "  /* Geometric data */\n"
3974: "  __local %s        detJ[%d]; //[N_t];           // $|J(x_q)|$, Jacobian determinant at $x_q$\n"
3975: "  __local %s        invJ[%d];//[N_t*dim*dim];   // $J^{-1}(x_q)$, Jacobian inverse at $x_q$\n",
3976:                             &count, numeric_str, numeric_str, N_b*N_c*N_q, numeric_str, dim, N_b*N_c*N_q, numeric_str, N_t,
3977:                             numeric_str, N_t*dim*dim, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
3978:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3979: "  /* FEM data */\n"
3980: "  __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",
3981:                             &count, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
3982:   if (useAux) {
3983:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3984: "  __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",
3985:                             &count, numeric_str, N_t);STRING_ERROR_CHECK("Message to short");
3986:   }
3987:   if (useF0) {
3988:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3989: "  /* Intermediate calculations */\n"
3990: "  __local %s         f_0[%d]; //[N_t*N_sqc];      // $f_0(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
3991:                               &count, numeric_str, N_t*N_q);STRING_ERROR_CHECK("Message to short");
3992:   }
3993:   if (useF1) {
3994:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3995: "  __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",
3996:                               &count, numeric_str, dim, N_t*N_q);STRING_ERROR_CHECK("Message to short");
3997:   }
3998:   /* TODO: If using elasticity, put in mu/lambda coefficients */
3999:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4000: "  /* Output data */\n"
4001: "  %s                e_i;                 // Coefficient $e_i$ of the residual\n\n",
4002:                             &count, numeric_str);STRING_ERROR_CHECK("Message to short");
4003:   /* One-time loads */
4004:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4005: "  /* These should be generated inline */\n"
4006: "  /* Load quadrature weights */\n"
4007: "  w = weights[qidx];\n"
4008: "  /* Load basis tabulation \\phi_i for this cell */\n"
4009: "  if (tidx < N_bt*N_q) {\n"
4010: "    phi_i[tidx]    = Basis[tidx];\n"
4011: "    phiDer_i[tidx] = BasisDerivatives[tidx];\n"
4012: "  }\n\n",
4013:                        &count);STRING_ERROR_CHECK("Message to short");
4014:   /* Batch loads */
4015:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4016: "  for (int batch = 0; batch < N_cb; ++batch) {\n"
4017: "    /* Load geometry */\n"
4018: "    detJ[tidx] = jacobianDeterminants[Goffset+batch*N_bc+tidx];\n"
4019: "    for (int n = 0; n < dim*dim; ++n) {\n"
4020: "      const int offset = n*N_t;\n"
4021: "      invJ[offset+tidx] = jacobianInverses[(Goffset+batch*N_bc)*dim*dim+offset+tidx];\n"
4022: "    }\n"
4023: "    /* Load coefficients u_i for this cell */\n"
4024: "    for (int n = 0; n < N_bt; ++n) {\n"
4025: "      const int offset = n*N_t;\n"
4026: "      u_i[offset+tidx] = coefficients[(Goffset*N_bt)+batch*N_t*N_b+offset+tidx];\n"
4027: "    }\n",
4028:                        &count);STRING_ERROR_CHECK("Message to short");
4029:   if (useAux) {
4030:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4031: "    /* Load coefficients a_i for this cell */\n"
4032: "    /* TODO: This should not be N_t here, it should be N_bc*N_comp_aux */\n"
4033: "    a_i[tidx] = coefficientsAux[Goffset+batch*N_t+tidx];\n",
4034:                             &count);STRING_ERROR_CHECK("Message to short");
4035:   }
4036:   /* Quadrature phase */
4037:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4038: "    barrier(CLK_LOCAL_MEM_FENCE);\n"
4039: "\n"
4040: "    /* Map coefficients to values at quadrature points */\n"
4041: "    for (int c = 0; c < N_sqc; ++c) {\n"
4042: "      const int cell          = c*N_bl*N_b + blqidx;\n"
4043: "      const int fidx          = (cell*N_q + qidx)*N_comp + cidx;\n",
4044:                        &count);STRING_ERROR_CHECK("Message to short");
4045:   if (useField) {
4046:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4047: "      %s  u[%d]; //[N_comp];     // $u(x_q)$, Value of the field at $x_q$\n",
4048:                               &count, numeric_str, N_c);STRING_ERROR_CHECK("Message to short");
4049:   }
4050:   if (useFieldDer) {
4051:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4052: "      %s%d   gradU[%d]; //[N_comp]; // $\\nabla u(x_q)$, Value of the field gradient at $x_q$\n",
4053:                               &count, numeric_str, dim, N_c);STRING_ERROR_CHECK("Message to short");
4054:   }
4055:   if (useFieldAux) {
4056:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4057: "      %s  a[%d]; //[1];     // $a(x_q)$, Value of the auxiliary fields at $x_q$\n",
4058:                               &count, numeric_str, 1);STRING_ERROR_CHECK("Message to short");
4059:   }
4060:   if (useFieldDerAux) {
4061:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4062: "      %s%d   gradA[%d]; //[1]; // $\\nabla a(x_q)$, Value of the auxiliary field gradient at $x_q$\n",
4063:                               &count, numeric_str, dim, 1);STRING_ERROR_CHECK("Message to short");
4064:   }
4065:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4066: "\n"
4067: "      for (int comp = 0; comp < N_comp; ++comp) {\n",
4068:                             &count);STRING_ERROR_CHECK("Message to short");
4069:   if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        u[comp] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4070:   if (useFieldDer) {
4071:     switch (dim) {
4072:     case 1:
4073:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        gradU[comp].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4074:     case 2:
4075:       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;
4076:     case 3:
4077:       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;
4078:     }
4079:   }
4080:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4081: "      }\n",
4082:                             &count);STRING_ERROR_CHECK("Message to short");
4083:   if (useFieldAux) {
4084:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      a[0] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");
4085:   }
4086:   if (useFieldDerAux) {
4087:     switch (dim) {
4088:     case 1:
4089:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      gradA[0].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4090:     case 2:
4091:       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;
4092:     case 3:
4093:       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;
4094:     }
4095:   }
4096:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4097: "      /* Get field and derivatives at this quadrature point */\n"
4098: "      for (int i = 0; i < N_b; ++i) {\n"
4099: "        for (int comp = 0; comp < N_comp; ++comp) {\n"
4100: "          const int b    = i*N_comp+comp;\n"
4101: "          const int pidx = qidx*N_bt + b;\n"
4102: "          const int uidx = cell*N_bt + b;\n"
4103: "          %s%d   realSpaceDer;\n\n",
4104:                             &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4105:   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");}
4106:   if (useFieldDer) {
4107:     switch (dim) {
4108:     case 2:
4109:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4110: "          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"
4111: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4112: "          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"
4113: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n",
4114:                            &count);STRING_ERROR_CHECK("Message to short");break;
4115:     case 3:
4116:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4117: "          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"
4118: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4119: "          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"
4120: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n"
4121: "          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"
4122: "          gradU[comp].z += u_i[uidx]*realSpaceDer.z;\n",
4123:                            &count);STRING_ERROR_CHECK("Message to short");break;
4124:     }
4125:   }
4126:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4127: "        }\n"
4128: "      }\n",
4129:                             &count);STRING_ERROR_CHECK("Message to short");
4130:   if (useFieldAux) {
4131:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"          a[0] += a_i[cell];\n", &count);STRING_ERROR_CHECK("Message to short");
4132:   }
4133:   /* Calculate residual at quadrature points: Should be generated by an weak form egine */
4134:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4135: "      /* Process values at quadrature points */\n",
4136:                             &count);STRING_ERROR_CHECK("Message to short");
4137:   switch (op) {
4138:   case LAPLACIAN:
4139:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4140:     if (useF1) {
4141:       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");}
4142:       else        {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_1[fidx] = gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4143:     }
4144:     break;
4145:   case ELASTICITY:
4146:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4147:     if (useF1) {
4148:     switch (dim) {
4149:     case 2:
4150:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4151: "      switch (cidx) {\n"
4152: "      case 0:\n"
4153: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].x + gradU[0].x);\n"
4154: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].y + gradU[1].x);\n"
4155: "        break;\n"
4156: "      case 1:\n"
4157: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].x + gradU[0].y);\n"
4158: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].y + gradU[1].y);\n"
4159: "      }\n",
4160:                            &count);STRING_ERROR_CHECK("Message to short");break;
4161:     case 3:
4162:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4163: "      switch (cidx) {\n"
4164: "      case 0:\n"
4165: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].x + gradU[0].x);\n"
4166: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].y + gradU[1].x);\n"
4167: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].z + gradU[2].x);\n"
4168: "        break;\n"
4169: "      case 1:\n"
4170: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].x + gradU[0].y);\n"
4171: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[1].y);\n"
4172: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[2].y);\n"
4173: "        break;\n"
4174: "      case 2:\n"
4175: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].x + gradU[0].z);\n"
4176: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[1].z);\n"
4177: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[2].z);\n"
4178: "      }\n",
4179:                            &count);STRING_ERROR_CHECK("Message to short");break;
4180:     }}
4181:     break;
4182:   default:
4183:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PDE operator %d is not supported", op);
4184:   }
4185:   if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_0[fidx] *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");}
4186:   if (useF1) {
4187:     switch (dim) {
4188:     case 1:
4189:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_1[fidx].x *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4190:     case 2:
4191:       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;
4192:     case 3:
4193:       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;
4194:     }
4195:   }
4196:   /* Thread transpose */
4197:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4198: "    }\n\n"
4199: "    /* ==== TRANSPOSE THREADS ==== */\n"
4200: "    barrier(CLK_LOCAL_MEM_FENCE);\n\n",
4201:                        &count);STRING_ERROR_CHECK("Message to short");
4202:   /* Basis phase */
4203:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4204: "    /* Map values at quadrature points to coefficients */\n"
4205: "    for (int c = 0; c < N_sbc; ++c) {\n"
4206: "      const int cell = c*N_bl*N_q + blbidx; /* Cell number in batch */\n"
4207: "\n"
4208: "      e_i = 0.0;\n"
4209: "      for (int q = 0; q < N_q; ++q) {\n"
4210: "        const int pidx = q*N_bt + bidx;\n"
4211: "        const int fidx = (cell*N_q + q)*N_comp + cidx;\n"
4212: "        %s%d   realSpaceDer;\n\n",
4213:                        &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");

4215:   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");}
4216:   if (useF1) {
4217:     switch (dim) {
4218:     case 2:
4219:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4220: "        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"
4221: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4222: "        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"
4223: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n",
4224:                            &count);STRING_ERROR_CHECK("Message to short");break;
4225:     case 3:
4226:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4227: "        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"
4228: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4229: "        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"
4230: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n"
4231: "        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"
4232: "        e_i           += realSpaceDer.z*f_1[fidx].z;\n",
4233:                            &count);STRING_ERROR_CHECK("Message to short");break;
4234:     }
4235:   }
4236:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4237: "      }\n"
4238: "      /* Write element vector for N_{cbc} cells at a time */\n"
4239: "      elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
4240: "    }\n"
4241: "    /* ==== Could do one write per batch ==== */\n"
4242: "  }\n"
4243: "  return;\n"
4244: "}\n",
4245:                        &count);STRING_ERROR_CHECK("Message to short");
4246:   return(0);
4247: }

4251: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
4252: {
4253:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4254:   PetscInt        dim, N_bl;
4255:   PetscBool       flg;
4256:   char           *buffer;
4257:   size_t          len;
4258:   char            errMsg[8192];
4259:   cl_int          ierr2;
4260:   PetscErrorCode  ierr;

4263:   PetscFEGetSpatialDimension(fem, &dim);
4264:   PetscMalloc1(8192, &buffer);
4265:   PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
4266:   PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
4267:   PetscOptionsHasName(fem->hdr.prefix, "-petscfe_opencl_kernel_print", &flg);
4268:   if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
4269:   len  = strlen(buffer);
4270:   *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
4271:   clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
4272:   if (ierr != CL_SUCCESS) {
4273:     clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
4274:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
4275:   }
4276:   PetscFree(buffer);
4277:   *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
4278:   return(0);
4279: }

4283: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
4284: {
4285:   const PetscInt Nblocks = N/blockSize;

4288:   if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
4289:   *z = 1;
4290:   for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
4291:     *y = Nblocks / *x;
4292:     if (*x * *y == Nblocks) break;
4293:   }
4294:   if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
4295:   return(0);
4296: }

4300: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
4301: {
4302:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4303:   PetscStageLog     stageLog;
4304:   PetscEventPerfLog eventLog = NULL;
4305:   PetscInt          stage;
4306:   PetscErrorCode    ierr;

4309:   PetscLogGetStageLog(&stageLog);
4310:   PetscStageLogGetCurrent(stageLog, &stage);
4311:   PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
4312:     /* Log performance info */
4313:   eventLog->eventInfo[ocl->residualEvent].count++;
4314:   eventLog->eventInfo[ocl->residualEvent].time  += time;
4315:   eventLog->eventInfo[ocl->residualEvent].flops += flops;
4316:   return(0);
4317: }

4321: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
4322:                                                const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4323: {
4324:   /* Nbc = batchSize */
4325:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4326:   void            (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
4327:   void            (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
4328:   PetscQuadrature   q;
4329:   PetscInt          dim;
4330:   PetscInt          N_b;    /* The number of basis functions */
4331:   PetscInt          N_comp; /* The number of basis function components */
4332:   PetscInt          N_bt;   /* The total number of scalar basis functions */
4333:   PetscInt          N_q;    /* The number of quadrature points */
4334:   PetscInt          N_bst;  /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
4335:   PetscInt          N_t;    /* The number of threads, N_bst * N_bl */
4336:   PetscInt          N_bl;   /* The number of blocks */
4337:   PetscInt          N_bc;   /* The batch size, N_bl*N_q*N_b */
4338:   PetscInt          N_cb;   /* The number of batches */
4339:   PetscInt          numFlops, f0Flops = 0, f1Flops = 0;
4340:   PetscBool         useAux      = coefficientsAux ? PETSC_TRUE : PETSC_FALSE;
4341:   PetscBool         useField    = PETSC_FALSE;
4342:   PetscBool         useFieldDer = PETSC_TRUE;
4343:   PetscBool         useF0       = PETSC_TRUE;
4344:   PetscBool         useF1       = PETSC_TRUE;
4345:   /* OpenCL variables */
4346:   cl_program        ocl_prog;
4347:   cl_kernel         ocl_kernel;
4348:   cl_event          ocl_ev;         /* The event for tracking kernel execution */
4349:   cl_ulong          ns_start;       /* Nanoseconds counter on GPU at kernel start */
4350:   cl_ulong          ns_end;         /* Nanoseconds counter on GPU at kernel stop */
4351:   cl_mem            o_jacobianInverses, o_jacobianDeterminants;
4352:   cl_mem            o_coefficients, o_coefficientsAux, o_elemVec;
4353:   float            *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
4354:   double           *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
4355:   void             *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
4356:   size_t            local_work_size[3], global_work_size[3];
4357:   size_t            realSize, x, y, z;
4358:   PetscErrorCode    ierr;

4361:   if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
4362:   PetscFEGetSpatialDimension(fem, &dim);
4363:   PetscFEGetQuadrature(fem, &q);
4364:   PetscFEGetDimension(fem, &N_b);
4365:   PetscFEGetNumComponents(fem, &N_comp);
4366:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4367:   PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
4368:   N_bt  = N_b*N_comp;
4369:   N_q   = q->numPoints;
4370:   N_bst = N_bt*N_q;
4371:   N_t   = N_bst*N_bl;
4372:   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);
4373:   /* Calculate layout */
4374:   if (Ne % (N_cb*N_bc)) { /* Remainder cells */
4375:     PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);
4376:     return(0);
4377:   }
4378:   PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
4379:   local_work_size[0]  = N_bc*N_comp;
4380:   local_work_size[1]  = 1;
4381:   local_work_size[2]  = 1;
4382:   global_work_size[0] = x * local_work_size[0];
4383:   global_work_size[1] = y * local_work_size[1];
4384:   global_work_size[2] = z * local_work_size[2];
4385:   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);
4386:   PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
4387:   /* Generate code */
4388:   if (probAux) {
4389:     PetscSpace P;
4390:     PetscInt   NfAux, order, f;

4392:     PetscDSGetNumFields(probAux, &NfAux);
4393:     for (f = 0; f < NfAux; ++f) {
4394:       PetscFE feAux;

4396:       PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
4397:       PetscFEGetBasisSpace(feAux, &P);
4398:       PetscSpaceGetOrder(P, &order);
4399:       if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
4400:     }
4401:   }
4402:   PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
4403:   /* Create buffers on the device and send data over */
4404:   PetscDataTypeGetSize(ocl->realType, &realSize);
4405:   if (sizeof(PetscReal) != realSize) {
4406:     switch (ocl->realType) {
4407:     case PETSC_FLOAT:
4408:     {
4409:       PetscInt c, b, d;

4411:       PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
4412:       for (c = 0; c < Ne; ++c) {
4413:         f_detJ[c] = (float) geom.detJ[c];
4414:         for (d = 0; d < dim*dim; ++d) {
4415:           f_invJ[c*dim*dim+d] = (float) geom.invJ[c*dim*dim+d];
4416:         }
4417:         for (b = 0; b < N_bt; ++b) {
4418:           f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
4419:         }
4420:       }
4421:       if (coefficientsAux) { /* Assume P0 */
4422:         for (c = 0; c < Ne; ++c) {
4423:           f_coeffAux[c] = (float) coefficientsAux[c];
4424:         }
4425:       }
4426:       oclCoeff      = (void *) f_coeff;
4427:       if (coefficientsAux) {
4428:         oclCoeffAux = (void *) f_coeffAux;
4429:       } else {
4430:         oclCoeffAux = NULL;
4431:       }
4432:       oclInvJ       = (void *) f_invJ;
4433:       oclDetJ       = (void *) f_detJ;
4434:     }
4435:     break;
4436:     case PETSC_DOUBLE:
4437:     {
4438:       PetscInt c, b, d;

4440:       PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
4441:       for (c = 0; c < Ne; ++c) {
4442:         d_detJ[c] = (double) geom.detJ[c];
4443:         for (d = 0; d < dim*dim; ++d) {
4444:           d_invJ[c*dim*dim+d] = (double) geom.invJ[c*dim*dim+d];
4445:         }
4446:         for (b = 0; b < N_bt; ++b) {
4447:           d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
4448:         }
4449:       }
4450:       if (coefficientsAux) { /* Assume P0 */
4451:         for (c = 0; c < Ne; ++c) {
4452:           d_coeffAux[c] = (double) coefficientsAux[c];
4453:         }
4454:       }
4455:       oclCoeff      = (void *) d_coeff;
4456:       if (coefficientsAux) {
4457:         oclCoeffAux = (void *) d_coeffAux;
4458:       } else {
4459:         oclCoeffAux = NULL;
4460:       }
4461:       oclInvJ       = (void *) d_invJ;
4462:       oclDetJ       = (void *) d_detJ;
4463:     }
4464:     break;
4465:     default:
4466:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4467:     }
4468:   } else {
4469:     oclCoeff    = (void *) coefficients;
4470:     oclCoeffAux = (void *) coefficientsAux;
4471:     oclInvJ     = (void *) geom.invJ;
4472:     oclDetJ     = (void *) geom.detJ;
4473:   }
4474:   o_coefficients         = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt    * realSize, oclCoeff,    &ierr);
4475:   if (coefficientsAux) {
4476:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclCoeffAux, &ierr);
4477:   } else {
4478:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY,                        Ne         * realSize, oclCoeffAux, &ierr);
4479:   }
4480:   o_jacobianInverses     = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ,     &ierr);
4481:   o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclDetJ,     &ierr);
4482:   o_elemVec              = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY,                       Ne*N_bt    * realSize, NULL,        &ierr);
4483:   /* Kernel launch */
4484:   clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
4485:   clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
4486:   clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
4487:   clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
4488:   clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
4489:   clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
4490:   clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
4491:   /* Read data back from device */
4492:   if (sizeof(PetscReal) != realSize) {
4493:     switch (ocl->realType) {
4494:     case PETSC_FLOAT:
4495:     {
4496:       float   *elem;
4497:       PetscInt c, b;

4499:       PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
4500:       PetscMalloc1(Ne*N_bt, &elem);
4501:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4502:       for (c = 0; c < Ne; ++c) {
4503:         for (b = 0; b < N_bt; ++b) {
4504:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4505:         }
4506:       }
4507:       PetscFree(elem);
4508:     }
4509:     break;
4510:     case PETSC_DOUBLE:
4511:     {
4512:       double  *elem;
4513:       PetscInt c, b;

4515:       PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
4516:       PetscMalloc1(Ne*N_bt, &elem);
4517:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4518:       for (c = 0; c < Ne; ++c) {
4519:         for (b = 0; b < N_bt; ++b) {
4520:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4521:         }
4522:       }
4523:       PetscFree(elem);
4524:     }
4525:     break;
4526:     default:
4527:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4528:     }
4529:   } else {
4530:     clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
4531:   }
4532:   /* Log performance */
4533:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
4534:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END,   sizeof(cl_ulong), &ns_end,   NULL);
4535:   f0Flops = 0;
4536:   switch (ocl->op) {
4537:   case LAPLACIAN:
4538:     f1Flops = useAux ? dim : 0;break;
4539:   case ELASTICITY:
4540:     f1Flops = 2*dim*dim;break;
4541:   }
4542:   numFlops = Ne*(
4543:     N_q*(
4544:       N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
4545:       /*+
4546:        N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
4547:       +
4548:       N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
4549:     +
4550:     N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
4551:   PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
4552:   /* Cleanup */
4553:   clReleaseMemObject(o_coefficients);
4554:   clReleaseMemObject(o_coefficientsAux);
4555:   clReleaseMemObject(o_jacobianInverses);
4556:   clReleaseMemObject(o_jacobianDeterminants);
4557:   clReleaseMemObject(o_elemVec);
4558:   clReleaseKernel(ocl_kernel);
4559:   clReleaseProgram(ocl_prog);
4560:   return(0);
4561: }

4565: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
4566: {
4568:   fem->ops->setfromoptions          = NULL;
4569:   fem->ops->setup                   = PetscFESetUp_Basic;
4570:   fem->ops->view                    = NULL;
4571:   fem->ops->destroy                 = PetscFEDestroy_OpenCL;
4572:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4573:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
4574:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_OpenCL;
4575:   fem->ops->integratebdresidual     = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
4576:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
4577:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
4578:   return(0);
4579: }

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

4584:   Level: intermediate

4586: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4587: M*/

4591: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
4592: {
4593:   PetscFE_OpenCL *ocl;
4594:   cl_uint         num_platforms;
4595:   cl_platform_id  platform_ids[42];
4596:   cl_uint         num_devices;
4597:   cl_device_id    device_ids[42];
4598:   cl_int          ierr2;
4599:   PetscErrorCode  ierr;

4603:   PetscNewLog(fem,&ocl);
4604:   fem->data = ocl;

4606:   /* Init Platform */
4607:   clGetPlatformIDs(42, platform_ids, &num_platforms);
4608:   if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
4609:   ocl->pf_id = platform_ids[0];
4610:   /* Init Device */
4611:   clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
4612:   if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
4613:   ocl->dev_id = device_ids[0];
4614:   /* Create context with one command queue */
4615:   ocl->ctx_id   = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
4616:   ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
4617:   /* Types */
4618:   ocl->realType = PETSC_FLOAT;
4619:   /* Register events */
4620:   PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
4621:   /* Equation handling */
4622:   ocl->op = LAPLACIAN;

4624:   PetscFEInitialize_OpenCL(fem);
4625:   return(0);
4626: }

4630: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
4631: {
4632:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

4636:   ocl->realType = realType;
4637:   return(0);
4638: }

4642: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
4643: {
4644:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

4649:   *realType = ocl->realType;
4650:   return(0);
4651: }

4653: #endif /* PETSC_HAVE_OPENCL */

4657: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
4658: {
4659:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4660:   PetscErrorCode     ierr;

4663:   CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
4664:   PetscFree(cmp->embedding);
4665:   PetscFree(cmp);
4666:   return(0);
4667: }

4671: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
4672: {
4673:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4674:   DM                 K;
4675:   PetscReal         *work, *subpoint;
4676:   PetscBLASInt      *pivots;
4677: #ifndef PETSC_USE_COMPLEX
4678:   PetscBLASInt       n, info;
4679: #endif
4680:   PetscInt           dim, pdim, spdim, j, s;
4681:   PetscErrorCode     ierr;

4684:   /* Get affine mapping from reference cell to each subcell */
4685:   PetscDualSpaceGetDM(fem->dualSpace, &K);
4686:   DMGetDimension(K, &dim);
4687:   DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
4688:   CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
4689:   /* Determine dof embedding into subelements */
4690:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
4691:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
4692:   PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
4693:   DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
4694:   for (s = 0; s < cmp->numSubelements; ++s) {
4695:     PetscInt sd = 0;

4697:     for (j = 0; j < pdim; ++j) {
4698:       PetscBool       inside;
4699:       PetscQuadrature f;
4700:       PetscInt        d, e;

4702:       PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
4703:       /* Apply transform to first point, and check that point is inside subcell */
4704:       for (d = 0; d < dim; ++d) {
4705:         subpoint[d] = -1.0;
4706:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
4707:       }
4708:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
4709:       if (inside) {cmp->embedding[s*spdim+sd++] = j;}
4710:     }
4711:     if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
4712:   }
4713:   DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
4714:   /* Construct the change of basis from prime basis to nodal basis for each subelement */
4715:   PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
4716:   PetscMalloc2(spdim,&pivots,spdim,&work);
4717:   for (s = 0; s < cmp->numSubelements; ++s) {
4718:     for (j = 0; j < spdim; ++j) {
4719:       PetscReal      *Bf;
4720:       PetscQuadrature f;
4721:       PetscInt        q, k;

4723:       PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
4724:       PetscMalloc1(f->numPoints*spdim,&Bf);
4725:       PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
4726:       for (k = 0; k < spdim; ++k) {
4727:         /* n_j \cdot \phi_k */
4728:         fem->invV[(s*spdim + j)*spdim+k] = 0.0;
4729:         for (q = 0; q < f->numPoints; ++q) {
4730:           fem->invV[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
4731:         }
4732:       }
4733:       PetscFree(Bf);
4734:     }
4735: #ifndef PETSC_USE_COMPLEX
4736:     n = spdim;
4737:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &fem->invV[s*spdim*spdim], &n, pivots, &info));
4738:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &fem->invV[s*spdim*spdim], &n, pivots, work, &n, &info));
4739: #endif
4740:   }
4741:   PetscFree2(pivots,work);
4742:   return(0);
4743: }

4747: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
4748: {
4749:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4750:   DM                 dm;
4751:   PetscInt           pdim;  /* Dimension of FE space P */
4752:   PetscInt           spdim; /* Dimension of subelement FE space P */
4753:   PetscInt           dim;   /* Spatial dimension */
4754:   PetscInt           comp;  /* Field components */
4755:   PetscInt          *subpoints;
4756:   PetscReal         *tmpB, *tmpD, *tmpH, *subpoint;
4757:   PetscInt           p, s, d, e, j, k;
4758:   PetscErrorCode     ierr;

4761:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
4762:   DMGetDimension(dm, &dim);
4763:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
4764:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
4765:   PetscFEGetNumComponents(fem, &comp);
4766:   /* Divide points into subelements */
4767:   DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
4768:   DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
4769:   for (p = 0; p < npoints; ++p) {
4770:     for (s = 0; s < cmp->numSubelements; ++s) {
4771:       PetscBool inside;

4773:       /* Apply transform, and check that point is inside cell */
4774:       for (d = 0; d < dim; ++d) {
4775:         subpoint[d] = -1.0;
4776:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
4777:       }
4778:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
4779:       if (inside) {subpoints[p] = s; break;}
4780:     }
4781:     if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
4782:   }
4783:   DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
4784:   /* Evaluate the prime basis functions at all points */
4785:   if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
4786:   if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
4787:   if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
4788:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
4789:   /* Translate to the nodal basis */
4790:   if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
4791:   if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
4792:   if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
4793:   for (p = 0; p < npoints; ++p) {
4794:     const PetscInt s = subpoints[p];

4796:     if (B) {
4797:       /* Multiply by V^{-1} (spdim x spdim) */
4798:       for (j = 0; j < spdim; ++j) {
4799:         const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
4800:         PetscInt       c;

4802:         B[i] = 0.0;
4803:         for (k = 0; k < spdim; ++k) {
4804:           B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
4805:         }
4806:         for (c = 1; c < comp; ++c) {
4807:           B[i+c] = B[i];
4808:         }
4809:       }
4810:     }
4811:     if (D) {
4812:       /* Multiply by V^{-1} (spdim x spdim) */
4813:       for (j = 0; j < spdim; ++j) {
4814:         for (d = 0; d < dim; ++d) {
4815:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
4816:           PetscInt       c;

4818:           D[i] = 0.0;
4819:           for (k = 0; k < spdim; ++k) {
4820:             D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
4821:           }
4822:           for (c = 1; c < comp; ++c) {
4823:             D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
4824:           }
4825:         }
4826:       }
4827:     }
4828:     if (H) {
4829:       /* Multiply by V^{-1} (pdim x pdim) */
4830:       for (j = 0; j < spdim; ++j) {
4831:         for (d = 0; d < dim*dim; ++d) {
4832:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
4833:           PetscInt       c;

4835:           H[i] = 0.0;
4836:           for (k = 0; k < spdim; ++k) {
4837:             H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
4838:           }
4839:           for (c = 1; c < comp; ++c) {
4840:             H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
4841:           }
4842:         }
4843:       }
4844:     }
4845:   }
4846:   DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
4847:   if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
4848:   if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
4849:   if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
4850:   return(0);
4851: }

4855: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
4856: {
4858:   fem->ops->setfromoptions          = NULL;
4859:   fem->ops->setup                   = PetscFESetUp_Composite;
4860:   fem->ops->view                    = NULL;
4861:   fem->ops->destroy                 = PetscFEDestroy_Composite;
4862:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4863:   fem->ops->gettabulation           = PetscFEGetTabulation_Composite;
4864:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
4865:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
4866:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
4867:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
4868:   return(0);
4869: }

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

4874:   Level: intermediate

4876: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4877: M*/

4881: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
4882: {
4883:   PetscFE_Composite *cmp;
4884:   PetscErrorCode     ierr;

4888:   PetscNewLog(fem, &cmp);
4889:   fem->data = cmp;

4891:   cmp->cellRefiner    = 0;
4892:   cmp->numSubelements = -1;
4893:   cmp->v0             = NULL;
4894:   cmp->jac            = NULL;

4896:   PetscFEInitialize_Composite(fem);
4897:   return(0);
4898: }

4902: PetscErrorCode PetscFECompositeExpandQuadrature(PetscFE fem, PetscQuadrature q, PetscQuadrature *qref)
4903: {
4904:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4905:   const PetscReal   *points,    *weights;
4906:   PetscReal         *pointsRef, *weightsRef;
4907:   PetscInt           dim, order, npoints, npointsRef, c, p, d, e;
4908:   PetscErrorCode     ierr;

4914:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
4915:   PetscQuadratureGetOrder(q, &order);
4916:   PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
4917:   npointsRef = npoints*cmp->numSubelements;
4918:   PetscMalloc1(npointsRef*dim,&pointsRef);
4919:   PetscMalloc1(npointsRef,&weightsRef);
4920:   for (c = 0; c < cmp->numSubelements; ++c) {
4921:     for (p = 0; p < npoints; ++p) {
4922:       for (d = 0; d < dim; ++d) {
4923:         pointsRef[(c*npoints + p)*dim+d] = cmp->v0[c*dim+d];
4924:         for (e = 0; e < dim; ++e) {
4925:           pointsRef[(c*npoints + p)*dim+d] += cmp->jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
4926:         }
4927:       }
4928:       /* Could also use detJ here */
4929:       weightsRef[c*npoints+p] = weights[p]/cmp->numSubelements;
4930:     }
4931:   }
4932:   PetscQuadratureSetOrder(*qref, order);
4933:   PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
4934:   return(0);
4935: }

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

4942:   Not collective

4944:   Input Parameter:
4945: . fe - The PetscFE

4947:   Output Parameter:
4948: . dim - The dimension

4950:   Level: intermediate

4952: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
4953: @*/
4954: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
4955: {

4961:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
4962:   return(0);
4963: }

4965: /*
4966: Purpose: Compute element vector for chunk of elements

4968: Input:
4969:   Sizes:
4970:      Ne:  number of elements
4971:      Nf:  number of fields
4972:      PetscFE
4973:        dim: spatial dimension
4974:        Nb:  number of basis functions
4975:        Nc:  number of field components
4976:        PetscQuadrature
4977:          Nq:  number of quadrature points

4979:   Geometry:
4980:      PetscCellGeometry
4981:        PetscReal v0s[Ne*dim]
4982:        PetscReal jacobians[Ne*dim*dim]        possibly *Nq
4983:        PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
4984:        PetscReal jacobianDeterminants[Ne]     possibly *Nq
4985:   FEM:
4986:      PetscFE
4987:        PetscQuadrature
4988:          PetscReal   quadPoints[Nq*dim]
4989:          PetscReal   quadWeights[Nq]
4990:        PetscReal   basis[Nq*Nb*Nc]
4991:        PetscReal   basisDer[Nq*Nb*Nc*dim]
4992:      PetscScalar coefficients[Ne*Nb*Nc]
4993:      PetscScalar elemVec[Ne*Nb*Nc]

4995:   Problem:
4996:      PetscInt f: the active field
4997:      f0, f1

4999:   Work Space:
5000:      PetscFE
5001:        PetscScalar f0[Nq*dim];
5002:        PetscScalar f1[Nq*dim*dim];
5003:        PetscScalar u[Nc];
5004:        PetscScalar gradU[Nc*dim];
5005:        PetscReal   x[dim];
5006:        PetscScalar realSpaceDer[dim];

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

5010: Input:
5011:   Sizes:
5012:      N_cb: Number of serial cell batches

5014:   Geometry:
5015:      PetscReal v0s[Ne*dim]
5016:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
5017:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5018:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
5019:   FEM:
5020:      static PetscReal   quadPoints[Nq*dim]
5021:      static PetscReal   quadWeights[Nq]
5022:      static PetscReal   basis[Nq*Nb*Nc]
5023:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
5024:      PetscScalar coefficients[Ne*Nb*Nc]
5025:      PetscScalar elemVec[Ne*Nb*Nc]

5027: ex62.c:
5028:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
5029:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
5030:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
5031:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

5033: ex52.c:
5034:   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)
5035:   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)

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

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

5044: ex52_integrateElementOpenCL.c:
5045: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
5046:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5047:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

5057:   Not collective

5059:   Input Parameters:
5060: + fem          - The PetscFE object for the field being integrated
5061: . prob         - The PetscDS specifing the discretizations and continuum functions
5062: . field        - The field being integrated
5063: . Ne           - The number of elements in the chunk
5064: . geom         - The cell geometry for each cell in the chunk
5065: . coefficients - The array of FEM basis coefficients for the elements
5066: . probAux      - The PetscDS specifing the auxiliary discretizations
5067: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5069:   Output Parameter
5070: . integral     - the integral for this field

5072:   Level: developer

5074: .seealso: PetscFEIntegrateResidual()
5075: @*/
5076: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5077:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
5078: {

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

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

5093:   Not collective

5095:   Input Parameters:
5096: + fem          - The PetscFE object for the field being integrated
5097: . prob         - The PetscDS specifing the discretizations and continuum functions
5098: . field        - The field being integrated
5099: . Ne           - The number of elements in the chunk
5100: . geom         - The cell geometry for each cell in the chunk
5101: . coefficients - The array of FEM basis coefficients for the elements
5102: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5103: . probAux      - The PetscDS specifing the auxiliary discretizations
5104: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5106:   Output Parameter
5107: . elemVec      - the element residual vectors from each element

5109:   Note:
5110: $ Loop over batch of elements (e):
5111: $   Loop over quadrature points (q):
5112: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
5113: $     Call f_0 and f_1
5114: $   Loop over element vector entries (f,fc --> i):
5115: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

5117:   Level: developer

5119: .seealso: PetscFEIntegrateResidual()
5120: @*/
5121: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5122:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5123: {

5129:   if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5130:   return(0);
5131: }

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

5138:   Not collective

5140:   Input Parameters:
5141: + fem          - The PetscFE object for the field being integrated
5142: . prob         - The PetscDS specifing the discretizations and continuum functions
5143: . field        - The field being integrated
5144: . Ne           - The number of elements in the chunk
5145: . geom         - The cell geometry for each cell in the chunk
5146: . coefficients - The array of FEM basis coefficients for the elements
5147: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5148: . probAux      - The PetscDS specifing the auxiliary discretizations
5149: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5151:   Output Parameter
5152: . elemVec      - the element residual vectors from each element

5154:   Level: developer

5156: .seealso: PetscFEIntegrateResidual()
5157: @*/
5158: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5159:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5160: {

5165:   if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5166:   return(0);
5167: }

5171: /*C
5172:   PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration

5174:   Not collective

5176:   Input Parameters:
5177: + fem          = The PetscFE object for the field being integrated
5178: . prob         - The PetscDS specifing the discretizations and continuum functions
5179: . fieldI       - The test field being integrated
5180: . fieldJ       - The basis field being integrated
5181: . Ne           - The number of elements in the chunk
5182: . geom         - The cell geometry for each cell in the chunk
5183: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5184: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5185: . probAux      - The PetscDS specifing the auxiliary discretizations
5186: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5188:   Output Parameter
5189: . elemMat      - the element matrices for the Jacobian from each element

5191:   Note:
5192: $ Loop over batch of elements (e):
5193: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5194: $     Loop over quadrature points (q):
5195: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5196: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5197: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5198: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5199: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5200: */
5201: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
5202:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5203: {

5208:   if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5209:   return(0);
5210: }

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

5217:   Not collective

5219:   Input Parameters:
5220: + fem          = The PetscFE object for the field being integrated
5221: . prob         - The PetscDS specifing the discretizations and continuum functions
5222: . fieldI       - The test field being integrated
5223: . fieldJ       - The basis field being integrated
5224: . Ne           - The number of elements in the chunk
5225: . geom         - The cell geometry for each cell in the chunk
5226: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5227: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5228: . probAux      - The PetscDS specifing the auxiliary discretizations
5229: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5231:   Output Parameter
5232: . elemMat              - the element matrices for the Jacobian from each element

5234:   Note:
5235: $ Loop over batch of elements (e):
5236: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5237: $     Loop over quadrature points (q):
5238: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5239: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5240: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5241: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5242: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5243: */
5244: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
5245:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5246: {

5251:   if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5252:   return(0);
5253: }

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

5262:   Input Parameter:
5263: . fe - The initial PetscFE

5265:   Output Parameter:
5266: . feRef - The refined PetscFE

5268:   Level: developer

5270: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5271: @*/
5272: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
5273: {
5274:   PetscSpace       P, Pref;
5275:   PetscDualSpace   Q, Qref;
5276:   DM               K, Kref;
5277:   PetscQuadrature  q, qref;
5278:   PetscInt         numComp;
5279:   PetscErrorCode   ierr;

5282:   PetscFEGetBasisSpace(fe, &P);
5283:   PetscFEGetDualSpace(fe, &Q);
5284:   PetscFEGetQuadrature(fe, &q);
5285:   PetscDualSpaceGetDM(Q, &K);
5286:   /* Create space */
5287:   PetscObjectReference((PetscObject) P);
5288:   Pref = P;
5289:   /* Create dual space */
5290:   PetscDualSpaceDuplicate(Q, &Qref);
5291:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
5292:   PetscDualSpaceSetDM(Qref, Kref);
5293:   DMDestroy(&Kref);
5294:   PetscDualSpaceSetUp(Qref);
5295:   /* Create element */
5296:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
5297:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
5298:   PetscFESetBasisSpace(*feRef, Pref);
5299:   PetscFESetDualSpace(*feRef, Qref);
5300:   PetscFEGetNumComponents(fe,    &numComp);
5301:   PetscFESetNumComponents(*feRef, numComp);
5302:   PetscFESetUp(*feRef);
5303:   PetscSpaceDestroy(&Pref);
5304:   PetscDualSpaceDestroy(&Qref);
5305:   /* Create quadrature */
5306:   PetscFECompositeExpandQuadrature(*feRef, q, &qref);
5307:   PetscFESetQuadrature(*feRef, qref);
5308:   PetscQuadratureDestroy(&qref);
5309:   return(0);
5310: }

5314: /*@
5315:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

5317:   Collective on DM

5319:   Input Parameters:
5320: + dm         - The underlying DM for the domain
5321: . dim        - The spatial dimension
5322: . numComp    - The number of components
5323: . isSimplex  - Flag for simplex reference cell, otherwise its a tensor product
5324: . prefix     - The options prefix, or NULL
5325: - qorder     - The quadrature order

5327:   Output Parameter:
5328: . fem - The PetscFE object

5330:   Level: beginner

5332: .keywords: PetscFE, finite element
5333: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
5334: @*/
5335: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
5336: {
5337:   PetscQuadrature q;
5338:   DM              K;
5339:   PetscSpace      P;
5340:   PetscDualSpace  Q;
5341:   PetscInt        order;
5342:   PetscErrorCode  ierr;

5345:   /* Create space */
5346:   PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
5347:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
5348:   PetscSpaceSetFromOptions(P);
5349:   PetscSpacePolynomialSetNumVariables(P, dim);
5350:   PetscSpaceSetUp(P);
5351:   PetscSpaceGetOrder(P, &order);
5352:   /* Create dual space */
5353:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
5354:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
5355:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
5356:   PetscDualSpaceSetDM(Q, K);
5357:   DMDestroy(&K);
5358:   PetscDualSpaceSetOrder(Q, order);
5359:   PetscDualSpaceSetFromOptions(Q);
5360:   PetscDualSpaceSetUp(Q);
5361:   /* Create element */
5362:   PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
5363:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
5364:   PetscFESetFromOptions(*fem);
5365:   PetscFESetBasisSpace(*fem, P);
5366:   PetscFESetDualSpace(*fem, Q);
5367:   PetscFESetNumComponents(*fem, numComp);
5368:   PetscFESetUp(*fem);
5369:   PetscSpaceDestroy(&P);
5370:   PetscDualSpaceDestroy(&Q);
5371:   /* Create quadrature (with specified order if given) */
5372:   if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5373:   else           {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order+1, 1), -1.0, 1.0, &q);}
5374:   PetscFESetQuadrature(*fem, q);
5375:   PetscQuadratureDestroy(&q);
5376:   return(0);
5377: }