Actual source code: dtfe.c

petsc-3.5.2 2014-09-08
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:       if (ind[0] == max-1) {ind[0] = -1;}
636:       else                 {ind[1] = 0; ++ind[0];}
637:     }
638:   }
639:   return(0);
640: }

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

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

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

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

782:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1008:   Level: intermediate

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

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

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

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

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


1036: PetscClassId PETSCDUALSPACE_CLASSID = 0;

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

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

1046:   Not Collective

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

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

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

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

1070:   Level: advanced

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

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

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

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

1090:   Collective on PetscDualSpace

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

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

1099:   Level: intermediate

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

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

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

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

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

1133:   Not Collective

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

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

1141:   Level: intermediate

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

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

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

1165:   Collective on PetscDualSpace

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

1171:   Level: developer

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

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

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

1195:   Collective on PetscDualSpace

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

1202:   Level: intermediate

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

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

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

1234:   Collective on PetscDualSpace

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

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

1242:   Level: developer

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

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

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

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

1285:   Collective on PetscDualSpace

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

1290:   Level: developer

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

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

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

1309:   Collective on PetscDualSpace

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

1314:   Level: developer

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

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

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

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

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

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

1347:   Collective on MPI_Comm

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

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

1355:   Level: beginner

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

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

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

1373:   s->order = 0;

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

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

1384:   Collective on PetscDualSpace

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

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

1392:   Level: beginner

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

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

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

1412:   Not collective

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

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

1420:   Level: intermediate

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

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

1438:   Not collective

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

1444:   Level: intermediate

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

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

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

1466:   Not collective

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

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

1474:   Level: intermediate

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

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

1492:   Not collective

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

1498:   Level: intermediate

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

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

1515:   Not collective

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

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

1524:   Level: intermediate

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

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

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

1547:   Not collective

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

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

1555:   Level: intermediate

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

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

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

1576:   Not collective

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

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

1584:   Level: intermediate

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

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

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

1605:   Collective on PetscDualSpace

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

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

1615:   Level: advanced

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

1626:   DMCreate(PetscObjectComm((PetscObject) sp), &rdm);
1627:   DMSetType(rdm, DMPLEX);
1628:   DMPlexSetDimension(rdm, dim);
1629:   switch (dim) {
1630:   case 0:
1631:   {
1632:     PetscInt    numPoints[1]        = {1};
1633:     PetscInt    coneSize[1]         = {0};
1634:     PetscInt    cones[1]            = {0};
1635:     PetscInt    coneOrientations[1] = {0};
1636:     PetscScalar vertexCoords[1]     = {0.0};

1638:     DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1639:   }
1640:   break;
1641:   case 1:
1642:   {
1643:     PetscInt    numPoints[2]        = {2, 1};
1644:     PetscInt    coneSize[3]         = {2, 0, 0};
1645:     PetscInt    cones[2]            = {1, 2};
1646:     PetscInt    coneOrientations[2] = {0, 0};
1647:     PetscScalar vertexCoords[2]     = {-1.0,  1.0};

1649:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1650:   }
1651:   break;
1652:   case 2:
1653:     if (simplex) {
1654:       PetscInt    numPoints[2]        = {3, 1};
1655:       PetscInt    coneSize[4]         = {3, 0, 0, 0};
1656:       PetscInt    cones[3]            = {1, 2, 3};
1657:       PetscInt    coneOrientations[3] = {0, 0, 0};
1658:       PetscScalar vertexCoords[6]     = {-1.0, -1.0,  1.0, -1.0,  -1.0, 1.0};

1660:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1661:     } else {
1662:       PetscInt    numPoints[2]        = {4, 1};
1663:       PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
1664:       PetscInt    cones[4]            = {1, 2, 3, 4};
1665:       PetscInt    coneOrientations[4] = {0, 0, 0, 0};
1666:       PetscScalar vertexCoords[8]     = {-1.0, -1.0,  1.0, -1.0,  1.0, 1.0,  -1.0, 1.0};

1668:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1669:     }
1670:   break;
1671:   case 3:
1672:     if (simplex) {
1673:       PetscInt    numPoints[2]        = {4, 1};
1674:       PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
1675:       PetscInt    cones[4]            = {1, 3, 2, 4};
1676:       PetscInt    coneOrientations[4] = {0, 0, 0, 0};
1677:       PetscScalar vertexCoords[12]    = {-1.0, -1.0, -1.0,  1.0, -1.0, -1.0,  -1.0, 1.0, -1.0,  -1.0, -1.0, 1.0};

1679:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1680:     } else {
1681:       PetscInt    numPoints[2]        = {8, 1};
1682:       PetscInt    coneSize[9]         = {8, 0, 0, 0, 0, 0, 0, 0, 0};
1683:       PetscInt    cones[8]            = {1, 4, 3, 2, 5, 6, 7, 8};
1684:       PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
1685:       PetscScalar vertexCoords[24]    = {-1.0, -1.0, -1.0,  1.0, -1.0, -1.0,  1.0, 1.0, -1.0,  -1.0, 1.0, -1.0,
1686:                                          -1.0, -1.0,  1.0,  1.0, -1.0,  1.0,  1.0, 1.0,  1.0,  -1.0, 1.0,  1.0};

1688:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1689:     }
1690:   break;
1691:   default:
1692:     SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim);
1693:   }
1694:   *refdm = NULL;
1695:   DMPlexInterpolate(rdm, refdm);
1696:   DMPlexCopyCoordinates(rdm, *refdm);
1697:   DMDestroy(&rdm);
1698:   return(0);
1699: }

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

1706:   Input Parameters:
1707: + sp      - The PetscDualSpace object
1708: . f       - The basis functional index
1709: . geom    - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1710: . numComp - The number of components for the function
1711: . func    - The input function
1712: - ctx     - A context for the function

1714:   Output Parameter:
1715: . value   - numComp output values

1717:   Level: developer

1719: .seealso: PetscDualSpaceCreate()
1720: @*/
1721: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscCellGeometry geom, PetscInt numComp, void (*func)(const PetscReal [], PetscScalar *, void *), void *ctx, PetscScalar *value)
1722: {
1723:   DM               dm;
1724:   PetscQuadrature  quad;
1725:   const PetscReal *v0 = geom.v0;
1726:   const PetscReal *J  = geom.J;
1727:   PetscReal        x[3];
1728:   PetscScalar     *val;
1729:   PetscInt         dim, q, c;
1730:   PetscErrorCode   ierr;

1735:   PetscDualSpaceGetDM(sp, &dm);
1736:   DMPlexGetDimension(dm, &dim);
1737:   PetscDualSpaceGetFunctional(sp, f, &quad);
1738:   DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1739:   for (c = 0; c < numComp; ++c) value[c] = 0.0;
1740:   for (q = 0; q < quad->numPoints; ++q) {
1741:     CoordinatesRefToReal(dim, dim, v0, J, &quad->points[q*dim], x);
1742:     (*func)(x, val, ctx);
1743:     for (c = 0; c < numComp; ++c) {
1744:       value[c] += val[c]*quad->weights[q];
1745:     }
1746:   }
1747:   DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1748:   return(0);
1749: }

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

1761:   DMPlexGetDimension(sp->dm, &n);
1762:   if (lag->simplex || !lag->continuous) {
1763:     for (i = 1; i <= n; ++i) {
1764:       D *= ((PetscReal) (order+i))/i;
1765:     }
1766:     *dim = (PetscInt) (D + 0.5);
1767:   } else {
1768:     *dim = 1;
1769:     for (i = 0; i < n; ++i) *dim *= (order+1);
1770:   }
1771:   return(0);
1772: }

1776: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1777: {
1778:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1779:   DM                  dm    = sp->dm;
1780:   PetscInt            order = sp->order;
1781:   PetscBool           disc  = lag->continuous ? PETSC_FALSE : PETSC_TRUE;
1782:   PetscSection        csection;
1783:   Vec                 coordinates;
1784:   PetscReal          *qpoints, *qweights;
1785:   PetscInt           *closure = NULL, closureSize, c;
1786:   PetscInt            depth, dim, pdimMax, *pStart, *pEnd, cell, coneSize, d, n, f = 0;
1787:   PetscBool           simplex;
1788:   PetscErrorCode      ierr;

1791:   /* Classify element type */
1792:   DMPlexGetDimension(dm, &dim);
1793:   DMPlexGetDepth(dm, &depth);
1794:   PetscCalloc1(dim+1, &lag->numDof);
1795:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
1796:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1797:   DMPlexGetConeSize(dm, pStart[depth], &coneSize);
1798:   DMGetCoordinateSection(dm, &csection);
1799:   DMGetCoordinatesLocal(dm, &coordinates);
1800:   if      (coneSize == dim+1)    simplex = PETSC_TRUE;
1801:   else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
1802:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1803:   lag->simplex = simplex;
1804:   PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
1805:   pdimMax *= (pEnd[dim] - pStart[dim]);
1806:   PetscMalloc1(pdimMax, &sp->functional);
1807:   if (!dim) {
1808:     PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1809:     PetscMalloc1(1, &qpoints);
1810:     PetscMalloc1(1, &qweights);
1811:     PetscQuadratureSetData(sp->functional[f], PETSC_DETERMINE, 1, qpoints, qweights);
1812:     qpoints[0]  = 0.0;
1813:     qweights[0] = 1.0;
1814:     ++f;
1815:     lag->numDof[0] = 1;
1816:   } else {
1817:     PetscBT seen;

1819:     PetscBTCreate(pEnd[dim-1], &seen);
1820:     PetscBTMemzero(pEnd[dim-1], seen);
1821:     for (cell = pStart[depth]; cell < pEnd[depth]; ++cell) {
1822:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1823:       for (c = 0; c < closureSize*2; c += 2) {
1824:         const PetscInt p = closure[c];

1826:         if (PetscBTLookup(seen, p)) continue;
1827:         PetscBTSet(seen, p);
1828:         if ((p >= pStart[0]) && (p < pEnd[0])) {
1829:           /* Vertices */
1830:           const PetscScalar *coords;
1831:           PetscInt           dof, off, d;

1833:           if (order < 1 || disc) continue;
1834:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1835:           PetscMalloc1(dim, &qpoints);
1836:           PetscMalloc1(1, &qweights);
1837:           PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1838:           VecGetArrayRead(coordinates, &coords);
1839:           PetscSectionGetDof(csection, p, &dof);
1840:           PetscSectionGetOffset(csection, p, &off);
1841:           if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim);
1842:           for (d = 0; d < dof; ++d) {qpoints[d] = PetscRealPart(coords[off+d]);}
1843:           qweights[0] = 1.0;
1844:           ++f;
1845:           VecRestoreArrayRead(coordinates, &coords);
1846:           lag->numDof[0] = 1;
1847:         } else if ((p >= pStart[1]) && (p < pEnd[1])) {
1848:           /* Edges */
1849:           PetscScalar *coords;
1850:           PetscInt     num = order-1, k;

1852:           if (order < 2 || disc) continue;
1853:           coords = NULL;
1854:           DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);
1855:           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);
1856:           for (k = 1; k <= num; ++k) {
1857:             PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1858:             PetscMalloc1(dim, &qpoints);
1859:             PetscMalloc1(1, &qweights);
1860:             PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1861:             for (d = 0; d < dim; ++d) {qpoints[d] = k*PetscRealPart(coords[1*dim+d] - coords[0*dim+d])/order + PetscRealPart(coords[0*dim+d]);}
1862:             qweights[0] = 1.0;
1863:             ++f;
1864:           }
1865:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);
1866:           lag->numDof[1] = num;
1867:         } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) {
1868:           /* Faces */

1870:           if (disc) continue;
1871:           if ( simplex && (order < 3)) continue;
1872:           if (!simplex && (order < 2)) continue;
1873:           lag->numDof[depth-1] = 0;
1874:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces");
1875:         } else if ((p >= pStart[depth]) && (p < pEnd[depth])) {
1876:           /* Cells */
1877:           PetscInt     orderEff = lag->continuous && order ? (simplex ? order-3 : order-2) : order;
1878:           PetscReal    denom    = order ? (lag->continuous ? order : (simplex ? order+3 : order+2)) : (simplex ? 3 : 2);
1879:           PetscScalar *coords   = NULL;
1880:           PetscReal    dx = 2.0/denom, *v0, *J, *invJ, detJ;
1881:           PetscInt    *ind, *tup;
1882:           PetscInt     cdim, csize, d, d2, o;

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

1890:           PetscCalloc5(dim,&ind,dim,&tup,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1891:           DMPlexComputeCellGeometry(dm, p, v0, J, invJ, &detJ);
1892:           if (simplex || !lag->continuous) {
1893:             for (o = 0; o <= orderEff; ++o) {
1894:               PetscMemzero(ind, dim*sizeof(PetscInt));
1895:               while (ind[0] >= 0) {
1896:                 PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1897:                 PetscMalloc1(dim, &qpoints);
1898:                 PetscMalloc1(1,   &qweights);
1899:                 PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1900:                 LatticePoint_Internal(dim, o, ind, tup);
1901:                 for (d = 0; d < dim; ++d) {
1902:                   qpoints[d] = v0[d];
1903:                   for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
1904:                 }
1905:                 qweights[0] = 1.0;
1906:                 ++f;
1907:               }
1908:             }
1909:           } else {
1910:             while (ind[0] >= 0) {
1911:               PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1912:               PetscMalloc1(dim, &qpoints);
1913:               PetscMalloc1(1,   &qweights);
1914:               PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1915:               TensorPoint_Internal(dim, orderEff+1, ind, tup);
1916:               for (d = 0; d < dim; ++d) {
1917:                 qpoints[d] = v0[d];
1918:                 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d]+1)*dx);
1919:               }
1920:               qweights[0] = 1.0;
1921:               ++f;
1922:             }
1923:           }
1924:           PetscFree5(ind,tup,v0,J,invJ);
1925:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &csize, &coords);
1926:           lag->numDof[depth] = cdim;
1927:         }
1928:       }
1929:       DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);
1930:     }
1931:     PetscBTDestroy(&seen);
1932:   }
1933:   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);
1934:   PetscFree2(pStart,pEnd);
1935:   if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
1936:   return(0);
1937: }

1941: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
1942: {
1943:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1944:   PetscErrorCode      ierr;

1947:   PetscFree(lag->numDof);
1948:   PetscFree(lag);
1949:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
1950:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
1951:   return(0);
1952: }

1956: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
1957: {
1958:   PetscInt       order;
1959:   PetscBool      cont;

1963:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
1964:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
1965:   PetscDualSpaceGetOrder(sp, &order);
1966:   PetscDualSpaceSetOrder(*spNew, order);
1967:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
1968:   PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
1969:   return(0);
1970: }

1974: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscDualSpace sp)
1975: {
1976:   PetscBool      continuous, flg;

1980:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
1981:   PetscOptionsHead("PetscDualSpace Lagrange Options");
1982:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
1983:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
1984:   PetscOptionsTail();
1985:   return(0);
1986: }

1990: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
1991: {
1992:   DM              K;
1993:   const PetscInt *numDof;
1994:   PetscInt        spatialDim, Nc, size = 0, d;
1995:   PetscErrorCode  ierr;

1998:   PetscDualSpaceGetDM(sp, &K);
1999:   PetscDualSpaceGetNumDof(sp, &numDof);
2000:   DMPlexGetDimension(K, &spatialDim);
2001:   DMPlexGetHeightStratum(K, 0, NULL, &Nc);
2002:   if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
2003:   for (d = 0; d <= spatialDim; ++d) {
2004:     PetscInt pStart, pEnd;

2006:     DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
2007:     size += (pEnd-pStart)*numDof[d];
2008:   }
2009:   *dim = size;
2010:   return(0);
2011: }

2015: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
2016: {
2017:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2020:   *numDof = lag->numDof;
2021:   return(0);
2022: }

2026: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2027: {
2028:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2033:   *continuous = lag->continuous;
2034:   return(0);
2035: }

2039: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2040: {
2041:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2045:   lag->continuous = continuous;
2046:   return(0);
2047: }

2051: /*@
2052:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2054:   Not Collective

2056:   Input Parameter:
2057: . sp         - the PetscDualSpace

2059:   Output Parameter:
2060: . continuous - flag for element continuity

2062:   Level: intermediate

2064: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2065: .seealso: PetscDualSpaceLagrangeSetContinuity()
2066: @*/
2067: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2068: {

2074:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2075:   return(0);
2076: }

2080: /*@
2081:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2083:   Logically Collective on PetscDualSpace

2085:   Input Parameters:
2086: + sp         - the PetscDualSpace
2087: - continuous - flag for element continuity

2089:   Options Database:
2090: . -petscdualspace_lagrange_continuity <bool>

2092:   Level: intermediate

2094: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2095: .seealso: PetscDualSpaceLagrangeGetContinuity()
2096: @*/
2097: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2098: {

2104:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2105:   return(0);
2106: }

2110: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2111: {
2113:   sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
2114:   sp->ops->setup          = PetscDualSpaceSetUp_Lagrange;
2115:   sp->ops->view           = NULL;
2116:   sp->ops->destroy        = PetscDualSpaceDestroy_Lagrange;
2117:   sp->ops->duplicate      = PetscDualSpaceDuplicate_Lagrange;
2118:   sp->ops->getdimension   = PetscDualSpaceGetDimension_Lagrange;
2119:   sp->ops->getnumdof      = PetscDualSpaceGetNumDof_Lagrange;
2120:   return(0);
2121: }

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

2126:   Level: intermediate

2128: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2129: M*/

2133: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2134: {
2135:   PetscDualSpace_Lag *lag;
2136:   PetscErrorCode      ierr;

2140:   PetscNewLog(sp,&lag);
2141:   sp->data = lag;

2143:   lag->numDof     = NULL;
2144:   lag->simplex    = PETSC_TRUE;
2145:   lag->continuous = PETSC_TRUE;

2147:   PetscDualSpaceInitialize_Lagrange(sp);
2148:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2149:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2150:   return(0);
2151: }


2154: PetscClassId PETSCFE_CLASSID = 0;

2156: PetscFunctionList PetscFEList              = NULL;
2157: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

2161: /*@C
2162:   PetscFERegister - Adds a new PetscFE implementation

2164:   Not Collective

2166:   Input Parameters:
2167: + name        - The name of a new user-defined creation routine
2168: - create_func - The creation routine itself

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

2173:   Sample usage:
2174: .vb
2175:     PetscFERegister("my_fe", MyPetscFECreate);
2176: .ve

2178:   Then, your PetscFE type can be chosen with the procedural interface via
2179: .vb
2180:     PetscFECreate(MPI_Comm, PetscFE *);
2181:     PetscFESetType(PetscFE, "my_fe");
2182: .ve
2183:    or at runtime via the option
2184: .vb
2185:     -petscfe_type my_fe
2186: .ve

2188:   Level: advanced

2190: .keywords: PetscFE, register
2191: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()

2193: @*/
2194: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2195: {

2199:   PetscFunctionListAdd(&PetscFEList, sname, function);
2200:   return(0);
2201: }

2205: /*@C
2206:   PetscFESetType - Builds a particular PetscFE

2208:   Collective on PetscFE

2210:   Input Parameters:
2211: + fem  - The PetscFE object
2212: - name - The kind of FEM space

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

2217:   Level: intermediate

2219: .keywords: PetscFE, set, type
2220: .seealso: PetscFEGetType(), PetscFECreate()
2221: @*/
2222: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2223: {
2224:   PetscErrorCode (*r)(PetscFE);
2225:   PetscBool      match;

2230:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
2231:   if (match) return(0);

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

2237:   if (fem->ops->destroy) {
2238:     (*fem->ops->destroy)(fem);
2239:     fem->ops->destroy = NULL;
2240:   }
2241:   (*r)(fem);
2242:   PetscObjectChangeTypeName((PetscObject) fem, name);
2243:   return(0);
2244: }

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

2251:   Not Collective

2253:   Input Parameter:
2254: . fem  - The PetscFE

2256:   Output Parameter:
2257: . name - The PetscFE type name

2259:   Level: intermediate

2261: .keywords: PetscFE, get, type, name
2262: .seealso: PetscFESetType(), PetscFECreate()
2263: @*/
2264: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2265: {

2271:   if (!PetscFERegisterAllCalled) {
2272:     PetscFERegisterAll();
2273:   }
2274:   *name = ((PetscObject) fem)->type_name;
2275:   return(0);
2276: }

2280: /*@C
2281:   PetscFEView - Views a PetscFE

2283:   Collective on PetscFE

2285:   Input Parameter:
2286: + fem - the PetscFE object to view
2287: - v   - the viewer

2289:   Level: developer

2291: .seealso PetscFEDestroy()
2292: @*/
2293: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
2294: {

2299:   if (!v) {
2300:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
2301:   }
2302:   if (fem->ops->view) {
2303:     (*fem->ops->view)(fem, v);
2304:   }
2305:   return(0);
2306: }

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

2313:   Collective on PetscFE

2315:   Input Parameters:
2316: + fem    - the PetscFE
2317: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
2318: - optionname - option to activate viewing

2320:   Level: intermediate

2322: .keywords: PetscFE, view, options, database
2323: .seealso: VecViewFromOptions(), MatViewFromOptions()
2324: */
2325: PetscErrorCode PetscFEViewFromOptions(PetscFE fem, const char prefix[], const char optionname[])
2326: {
2327:   PetscViewer       viewer;
2328:   PetscViewerFormat format;
2329:   PetscBool         flg;
2330:   PetscErrorCode    ierr;

2333:   if (prefix) {
2334:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), prefix, optionname, &viewer, &format, &flg);
2335:   } else {
2336:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), ((PetscObject) fem)->prefix, optionname, &viewer, &format, &flg);
2337:   }
2338:   if (flg) {
2339:     PetscViewerPushFormat(viewer, format);
2340:     PetscFEView(fem, viewer);
2341:     PetscViewerPopFormat(viewer);
2342:     PetscViewerDestroy(&viewer);
2343:   }
2344:   return(0);
2345: }

2349: /*@
2350:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

2352:   Collective on PetscFE

2354:   Input Parameter:
2355: . fem - the PetscFE object to set options for

2357:   Options Database:
2358: . -petscfe_num_blocks  the number of cell blocks to integrate concurrently
2359: . -petscfe_num_batches the number of cell batches to integrate serially

2361:   Level: developer

2363: .seealso PetscFEView()
2364: @*/
2365: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
2366: {
2367:   const char    *defaultType;
2368:   char           name[256];
2369:   PetscBool      flg;

2374:   if (!((PetscObject) fem)->type_name) {
2375:     defaultType = PETSCFEBASIC;
2376:   } else {
2377:     defaultType = ((PetscObject) fem)->type_name;
2378:   }
2379:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

2381:   PetscObjectOptionsBegin((PetscObject) fem);
2382:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
2383:   if (flg) {
2384:     PetscFESetType(fem, name);
2385:   } else if (!((PetscObject) fem)->type_name) {
2386:     PetscFESetType(fem, defaultType);
2387:   }
2388:   PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
2389:   PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
2390:   if (fem->ops->setfromoptions) {
2391:     (*fem->ops->setfromoptions)(fem);
2392:   }
2393:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
2394:   PetscObjectProcessOptionsHandlers((PetscObject) fem);
2395:   PetscOptionsEnd();
2396:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
2397:   return(0);
2398: }

2402: /*@C
2403:   PetscFESetUp - Construct data structures for the PetscFE

2405:   Collective on PetscFE

2407:   Input Parameter:
2408: . fem - the PetscFE object to setup

2410:   Level: developer

2412: .seealso PetscFEView(), PetscFEDestroy()
2413: @*/
2414: PetscErrorCode PetscFESetUp(PetscFE fem)
2415: {

2420:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
2421:   return(0);
2422: }

2426: /*@
2427:   PetscFEDestroy - Destroys a PetscFE object

2429:   Collective on PetscFE

2431:   Input Parameter:
2432: . fem - the PetscFE object to destroy

2434:   Level: developer

2436: .seealso PetscFEView()
2437: @*/
2438: PetscErrorCode PetscFEDestroy(PetscFE *fem)
2439: {

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

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

2449:   PetscFree((*fem)->numDof);
2450:   PetscFree((*fem)->invV);
2451:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
2452:   PetscSpaceDestroy(&(*fem)->basisSpace);
2453:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
2454:   PetscQuadratureDestroy(&(*fem)->quadrature);

2456:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
2457:   PetscHeaderDestroy(fem);
2458:   return(0);
2459: }

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

2466:   Collective on MPI_Comm

2468:   Input Parameter:
2469: . comm - The communicator for the PetscFE object

2471:   Output Parameter:
2472: . fem - The PetscFE object

2474:   Level: beginner

2476: .seealso: PetscFESetType(), PETSCFEGALERKIN
2477: @*/
2478: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
2479: {
2480:   PetscFE        f;

2485:   PetscCitationsRegister(FECitation,&FEcite);
2486:   *fem = NULL;
2487:   PetscFEInitializePackage();

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

2492:   f->basisSpace    = NULL;
2493:   f->dualSpace     = NULL;
2494:   f->numComponents = 1;
2495:   f->numDof        = NULL;
2496:   f->invV          = NULL;
2497:   f->B             = NULL;
2498:   f->D             = NULL;
2499:   f->H             = NULL;
2500:   PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
2501:   f->blockSize     = 0;
2502:   f->numBlocks     = 1;
2503:   f->batchSize     = 0;
2504:   f->numBatches    = 1;

2506:   *fem = f;
2507:   return(0);
2508: }

2512: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
2513: {
2514:   DM             dm;

2520:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2521:   DMPlexGetDimension(dm, dim);
2522:   return(0);
2523: }

2527: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
2528: {
2531:   fem->numComponents = comp;
2532:   return(0);
2533: }

2537: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
2538: {
2542:   *comp = fem->numComponents;
2543:   return(0);
2544: }

2548: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
2549: {
2552:   fem->blockSize  = blockSize;
2553:   fem->numBlocks  = numBlocks;
2554:   fem->batchSize  = batchSize;
2555:   fem->numBatches = numBatches;
2556:   return(0);
2557: }

2561: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
2562: {
2569:   if (blockSize)  *blockSize  = fem->blockSize;
2570:   if (numBlocks)  *numBlocks  = fem->numBlocks;
2571:   if (batchSize)  *batchSize  = fem->batchSize;
2572:   if (numBatches) *numBatches = fem->numBatches;
2573:   return(0);
2574: }

2578: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
2579: {
2583:   *sp = fem->basisSpace;
2584:   return(0);
2585: }

2589: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
2590: {

2596:   PetscSpaceDestroy(&fem->basisSpace);
2597:   fem->basisSpace = sp;
2598:   PetscObjectReference((PetscObject) fem->basisSpace);
2599:   return(0);
2600: }

2604: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
2605: {
2609:   *sp = fem->dualSpace;
2610:   return(0);
2611: }

2615: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
2616: {

2622:   PetscDualSpaceDestroy(&fem->dualSpace);
2623:   fem->dualSpace = sp;
2624:   PetscObjectReference((PetscObject) fem->dualSpace);
2625:   return(0);
2626: }

2630: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
2631: {
2635:   *q = fem->quadrature;
2636:   return(0);
2637: }

2641: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
2642: {

2647:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
2648:   PetscQuadratureDestroy(&fem->quadrature);
2649:   fem->quadrature = q;
2650:   PetscObjectReference((PetscObject) q);
2651:   return(0);
2652: }

2656: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
2657: {
2658:   const PetscInt *numDofDual;
2659:   PetscErrorCode  ierr;

2664:   PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
2665:   if (!fem->numDof) {
2666:     DM       dm;
2667:     PetscInt dim, d;

2669:     PetscDualSpaceGetDM(fem->dualSpace, &dm);
2670:     DMPlexGetDimension(dm, &dim);
2671:     PetscMalloc1((dim+1), &fem->numDof);
2672:     for (d = 0; d <= dim; ++d) {
2673:       fem->numDof[d] = fem->numComponents*numDofDual[d];
2674:     }
2675:   }
2676:   *numDof = fem->numDof;
2677:   return(0);
2678: }

2682: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
2683: {
2684:   PetscInt         npoints = fem->quadrature->numPoints;
2685:   const PetscReal *points  = fem->quadrature->points;
2686:   PetscErrorCode   ierr;

2693:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
2694:   if (B) *B = fem->B;
2695:   if (D) *D = fem->D;
2696:   if (H) *H = fem->H;
2697:   return(0);
2698: }

2702: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
2703: {
2704:   DM               dm;
2705:   PetscInt         pdim; /* Dimension of FE space P */
2706:   PetscInt         dim;  /* Spatial dimension */
2707:   PetscInt         comp; /* Field components */
2708:   PetscErrorCode   ierr;

2716:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2717:   DMPlexGetDimension(dm, &dim);
2718:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2719:   PetscFEGetNumComponents(fem, &comp);
2720:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
2721:   if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
2722:   if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
2723:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
2724:   return(0);
2725: }

2729: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
2730: {
2731:   DM             dm;

2736:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2737:   if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
2738:   if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
2739:   if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
2740:   return(0);
2741: }

2745: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
2746: {
2747:   PetscFE_Basic *b = (PetscFE_Basic *) fem->data;

2751:   PetscFree(b);
2752:   return(0);
2753: }

2757: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
2758: {
2759:   PetscSpace        basis;
2760:   PetscDualSpace    dual;
2761:   PetscQuadrature   q;
2762:   PetscInt          dim, Nq;
2763:   PetscViewerFormat format;
2764:   PetscErrorCode    ierr;

2767:   PetscFEGetBasisSpace(fe, &basis);
2768:   PetscFEGetDualSpace(fe, &dual);
2769:   PetscFEGetQuadrature(fe, &q);
2770:   PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
2771:   PetscViewerGetFormat(viewer, &format);
2772:   PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
2773:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2774:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
2775:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
2776:     PetscViewerASCIIPushTab(viewer);
2777:     PetscQuadratureView(q, viewer);
2778:     PetscViewerASCIIPopTab(viewer);
2779:   } else {
2780:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
2781:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
2782:   }
2783:   PetscViewerASCIIPushTab(viewer);
2784:   PetscSpaceView(basis, viewer);
2785:   PetscDualSpaceView(dual, viewer);
2786:   PetscViewerASCIIPopTab(viewer);
2787:   return(0);
2788: }

2792: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
2793: {
2794:   PetscBool      iascii;

2800:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2801:   if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
2802:   return(0);
2803: }

2807: /* Construct the change of basis from prime basis to nodal basis */
2808: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
2809: {
2810:   PetscReal     *work;
2811:   PetscBLASInt  *pivots;
2812: #ifndef PETSC_USE_COMPLEX
2813:   PetscBLASInt   n, info;
2814: #endif
2815:   PetscInt       pdim, j;

2819:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2820:   PetscMalloc1(pdim*pdim,&fem->invV);
2821:   for (j = 0; j < pdim; ++j) {
2822:     PetscReal      *Bf;
2823:     PetscQuadrature f;
2824:     PetscInt        q, k;

2826:     PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
2827:     PetscMalloc1(f->numPoints*pdim,&Bf);
2828:     PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
2829:     for (k = 0; k < pdim; ++k) {
2830:       /* n_j \cdot \phi_k */
2831:       fem->invV[j*pdim+k] = 0.0;
2832:       for (q = 0; q < f->numPoints; ++q) {
2833:         fem->invV[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
2834:       }
2835:     }
2836:     PetscFree(Bf);
2837:   }
2838:   PetscMalloc2(pdim,&pivots,pdim,&work);
2839: #ifndef PETSC_USE_COMPLEX
2840:   n = pdim;
2841:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, fem->invV, &n, pivots, &info));
2842:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
2843: #endif
2844:   PetscFree2(pivots,work);
2845:   return(0);
2846: }

2850: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
2851: {

2855:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
2856:   return(0);
2857: }

2861: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
2862: {
2863:   DM               dm;
2864:   PetscInt         pdim; /* Dimension of FE space P */
2865:   PetscInt         dim;  /* Spatial dimension */
2866:   PetscInt         comp; /* Field components */
2867:   PetscReal       *tmpB, *tmpD, *tmpH;
2868:   PetscInt         p, d, j, k;
2869:   PetscErrorCode   ierr;

2872:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2873:   DMPlexGetDimension(dm, &dim);
2874:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2875:   PetscFEGetNumComponents(fem, &comp);
2876:   /* Evaluate the prime basis functions at all points */
2877:   if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
2878:   if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
2879:   if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
2880:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
2881:   /* Translate to the nodal basis */
2882:   for (p = 0; p < npoints; ++p) {
2883:     if (B) {
2884:       /* Multiply by V^{-1} (pdim x pdim) */
2885:       for (j = 0; j < pdim; ++j) {
2886:         const PetscInt i = (p*pdim + j)*comp;
2887:         PetscInt       c;

2889:         B[i] = 0.0;
2890:         for (k = 0; k < pdim; ++k) {
2891:           B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
2892:         }
2893:         for (c = 1; c < comp; ++c) {
2894:           B[i+c] = B[i];
2895:         }
2896:       }
2897:     }
2898:     if (D) {
2899:       /* Multiply by V^{-1} (pdim x pdim) */
2900:       for (j = 0; j < pdim; ++j) {
2901:         for (d = 0; d < dim; ++d) {
2902:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
2903:           PetscInt       c;

2905:           D[i] = 0.0;
2906:           for (k = 0; k < pdim; ++k) {
2907:             D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
2908:           }
2909:           for (c = 1; c < comp; ++c) {
2910:             D[((p*pdim + j)*comp + c)*dim + d] = D[i];
2911:           }
2912:         }
2913:       }
2914:     }
2915:     if (H) {
2916:       /* Multiply by V^{-1} (pdim x pdim) */
2917:       for (j = 0; j < pdim; ++j) {
2918:         for (d = 0; d < dim*dim; ++d) {
2919:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
2920:           PetscInt       c;

2922:           H[i] = 0.0;
2923:           for (k = 0; k < pdim; ++k) {
2924:             H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
2925:           }
2926:           for (c = 1; c < comp; ++c) {
2927:             H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
2928:           }
2929:         }
2930:       }
2931:     }
2932:   }
2933:   if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
2934:   if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
2935:   if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
2936:   return(0);
2937: }

2941: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
2942:                                       const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
2943: {
2944:   const PetscInt  debug = 0;
2945:   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[]);
2946:   PetscQuadrature quad;
2947:   PetscScalar    *u, *u_x, *a, *a_x;
2948:   PetscReal      *x;
2949:   PetscInt        dim, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, e;
2950:   PetscErrorCode  ierr;

2953:   PetscDSGetObjective(prob, field, &obj_func);
2954:   if (!obj_func) return(0);
2955:   PetscFEGetSpatialDimension(fem, &dim);
2956:   PetscFEGetQuadrature(fem, &quad);
2957:   PetscDSGetTotalDimension(prob, &totDim);
2958:   PetscDSGetTotalDimension(probAux, &totDimAux);
2959:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
2960:   PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
2961:   PetscDSGetRefCoordArrays(prob, &x, NULL);
2962:   for (e = 0; e < Ne; ++e) {
2963:     const PetscReal *v0   = &geom.v0[e*dim];
2964:     const PetscReal *J    = &geom.J[e*dim*dim];
2965:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
2966:     const PetscReal  detJ = geom.detJ[e];
2967:     const PetscReal *quadPoints, *quadWeights;
2968:     PetscInt         Nq, q;

2970:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
2971:     if (debug > 1) {
2972:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
2973: #ifndef PETSC_USE_COMPLEX
2974:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
2975: #endif
2976:     }
2977:     for (q = 0; q < Nq; ++q) {
2978:       PetscScalar integrand;

2980:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
2981:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
2982:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       NULL, u, u_x, NULL);
2983:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
2984:       obj_func(u, NULL, u_x, a, NULL, a_x, x, &integrand);
2985:       integrand *= detJ*quadWeights[q];
2986:       integral[field] += PetscRealPart(integrand);
2987:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", PetscRealPart(integrand), integral[field]);}
2988:     }
2989:     cOffset    += totDim;
2990:     cOffsetAux += totDimAux;
2991:   }
2992:   return(0);
2993: }

2997: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
2998:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
2999: {
3000:   const PetscInt  debug = 0;
3001:   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[]);
3002:   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[]);
3003:   PetscQuadrature quad;
3004:   PetscReal     **basisField, **basisFieldDer;
3005:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3006:   PetscReal      *x;
3007:   PetscInt        dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3008:   PetscErrorCode  ierr;

3011:   PetscFEGetSpatialDimension(fem, &dim);
3012:   PetscFEGetQuadrature(fem, &quad);
3013:   PetscFEGetDimension(fem, &Nb);
3014:   PetscFEGetNumComponents(fem, &Nc);
3015:   PetscDSGetTotalDimension(prob, &totDim);
3016:   PetscDSGetFieldOffset(prob, field, &fOffset);
3017:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3018:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3019:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3020:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3021:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3022:   if (probAux) {
3023:     PetscDSGetTotalDimension(probAux, &totDimAux);
3024:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3025:   }
3026:   for (e = 0; e < Ne; ++e) {
3027:     const PetscReal *v0   = &geom.v0[e*dim];
3028:     const PetscReal *J    = &geom.J[e*dim*dim];
3029:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
3030:     const PetscReal  detJ = geom.detJ[e];
3031:     const PetscReal *quadPoints, *quadWeights;
3032:     PetscInt         Nq, q;

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

3059: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3060:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3061: {
3062:   const PetscInt  debug = 0;
3063:   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[]);
3064:   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[]);
3065:   PetscQuadrature quad;
3066:   PetscReal     **basisField, **basisFieldDer;
3067:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3068:   PetscReal      *x;
3069:   PetscInt        dim, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3070:   PetscErrorCode  ierr;

3073:   PetscFEGetSpatialDimension(fem, &dim);
3074:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3075:   PetscFEGetQuadrature(fem, &quad);
3076:   PetscFEGetDimension(fem, &Nb);
3077:   PetscFEGetNumComponents(fem, &Nc);
3078:   PetscDSGetTotalBdDimension(prob, &totDim);
3079:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3080:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3081:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3082:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3083:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3084:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3085:   if (probAux) {
3086:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3087:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3088:   }
3089:   for (e = 0; e < Ne; ++e) {
3090:     const PetscReal *v0          = &geom.v0[e*dim];
3091:     const PetscReal *n           = &geom.n[e*dim];
3092:     const PetscReal *J           = &geom.J[e*dim*dim];
3093:     const PetscReal *invJ        = &geom.invJ[e*dim*dim];
3094:     const PetscReal  detJ        = geom.detJ[e];
3095:     const PetscReal *quadPoints, *quadWeights;
3096:     PetscInt         Nq, q;

3098:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3099:     if (debug > 1) {
3100:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3101: #ifndef PETSC_USE_COMPLEX
3102:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3103: #endif
3104:     }
3105:     for (q = 0; q < Nq; ++q) {
3106:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3107:       CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3108:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3109:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3110:       f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3111:       f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3112:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3113:     }
3114:     UpdateElementVec(dim-1, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3115:     cOffset    += totDim;
3116:     cOffsetAux += totDimAux;
3117:   }
3118:   return(0);
3119: }

3123: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3124:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3125: {
3126:   const PetscInt  debug      = 0;
3127:   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[]);
3128:   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[]);
3129:   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[]);
3130:   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[]);
3131:   PetscFE         fe;
3132:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3133:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3134:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3135:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3136:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3137:   PetscQuadrature quad;
3138:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3139:   PetscReal      *x;
3140:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3141:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3142:   PetscInt        dim, totDim, totDimAux, e;
3143:   PetscErrorCode  ierr;

3146:   PetscFEGetSpatialDimension(fem, &dim);
3147:   PetscFEGetQuadrature(fem, &quad);
3148:   PetscDSGetTotalDimension(prob, &totDim);
3149:   PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3150:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3151:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3152:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3153:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3154:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3155:   PetscFEGetDimension(fe, &NbI);
3156:   PetscFEGetNumComponents(fe, &NcI);
3157:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3158:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3159:   PetscFEGetDimension(fe, &NbJ);
3160:   PetscFEGetNumComponents(fe, &NcJ);
3161:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3162:   if (probAux) {
3163:     PetscDSGetTotalDimension(probAux, &totDimAux);
3164:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3165:   }
3166:   basisI    = basisField[fieldI];
3167:   basisJ    = basisField[fieldJ];
3168:   basisDerI = basisFieldDer[fieldI];
3169:   basisDerJ = basisFieldDer[fieldJ];
3170:   /* Initialize here in case the function is not defined */
3171:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3172:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3173:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3174:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3175:   for (e = 0; e < Ne; ++e) {
3176:     const PetscReal *v0   = &geom.v0[e*dim];
3177:     const PetscReal *J    = &geom.J[e*dim*dim];
3178:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
3179:     const PetscReal  detJ = geom.detJ[e];
3180:     const PetscReal *quadPoints, *quadWeights;
3181:     PetscInt         Nq, q;

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

3187:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3188:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3189:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3190:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3191:       if (g0_func) {
3192:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3193:         g0_func(u, u_t, u_x, a, NULL, a_x, x, g0);
3194:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3195:       }
3196:       if (g1_func) {
3197:         PetscInt d, d2;
3198:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3199:         g1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3200:         for (fc = 0; fc < NcI; ++fc) {
3201:           for (gc = 0; gc < NcJ; ++gc) {
3202:             for (d = 0; d < dim; ++d) {
3203:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3204:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3205:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3206:             }
3207:           }
3208:         }
3209:       }
3210:       if (g2_func) {
3211:         PetscInt d, d2;
3212:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3213:         g2_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3214:         for (fc = 0; fc < NcI; ++fc) {
3215:           for (gc = 0; gc < NcJ; ++gc) {
3216:             for (d = 0; d < dim; ++d) {
3217:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3218:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3219:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3220:             }
3221:           }
3222:         }
3223:       }
3224:       if (g3_func) {
3225:         PetscInt d, d2, dp, d3;
3226:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3227:         g3_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3228:         for (fc = 0; fc < NcI; ++fc) {
3229:           for (gc = 0; gc < NcJ; ++gc) {
3230:             for (d = 0; d < dim; ++d) {
3231:               for (dp = 0; dp < dim; ++dp) {
3232:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3233:                 for (d2 = 0; d2 < dim; ++d2) {
3234:                   for (d3 = 0; d3 < dim; ++d3) {
3235:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3236:                   }
3237:                 }
3238:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3239:               }
3240:             }
3241:           }
3242:         }
3243:       }

3245:       for (f = 0; f < NbI; ++f) {
3246:         for (fc = 0; fc < NcI; ++fc) {
3247:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3248:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3249:           for (g = 0; g < NbJ; ++g) {
3250:             for (gc = 0; gc < NcJ; ++gc) {
3251:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3252:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3253:               const PetscInt fOff = eOffset+i*totDim+j;
3254:               PetscInt       d, d2;

3256:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3257:               for (d = 0; d < dim; ++d) {
3258:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3259:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3260:                 for (d2 = 0; d2 < dim; ++d2) {
3261:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3262:                 }
3263:               }
3264:             }
3265:           }
3266:         }
3267:       }
3268:     }
3269:     if (debug > 1) {
3270:       PetscInt fc, f, gc, g;

3272:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3273:       for (fc = 0; fc < NcI; ++fc) {
3274:         for (f = 0; f < NbI; ++f) {
3275:           const PetscInt i = offsetI + f*NcI+fc;
3276:           for (gc = 0; gc < NcJ; ++gc) {
3277:             for (g = 0; g < NbJ; ++g) {
3278:               const PetscInt j = offsetJ + g*NcJ+gc;
3279:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3280:             }
3281:           }
3282:           PetscPrintf(PETSC_COMM_SELF, "\n");
3283:         }
3284:       }
3285:     }
3286:     cOffset    += totDim;
3287:     cOffsetAux += totDimAux;
3288:     eOffset    += PetscSqr(totDim);
3289:   }
3290:   return(0);
3291: }

3295: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3296:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3297: {
3298:   const PetscInt  debug      = 0;
3299:   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[]);
3300:   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[]);
3301:   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[]);
3302:   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[]);
3303:   PetscFE         fe;
3304:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3305:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3306:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3307:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3308:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3309:   PetscQuadrature quad;
3310:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3311:   PetscReal      *x;
3312:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3313:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3314:   PetscInt        dim, bdim, totDim, totDimAux, e;
3315:   PetscErrorCode  ierr;

3318:   PetscFEGetQuadrature(fem, &quad);
3319:   PetscFEGetSpatialDimension(fem, &dim);
3320:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3321:   bdim = dim-1;
3322:   PetscDSGetTotalBdDimension(prob, &totDim);
3323:   PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3324:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3325:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3326:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3327:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3328:   PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
3329:   PetscFEGetDimension(fe, &NbI);
3330:   PetscFEGetNumComponents(fe, &NcI);
3331:   PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
3332:   PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
3333:   PetscFEGetDimension(fe, &NbJ);
3334:   PetscFEGetNumComponents(fe, &NcJ);
3335:   PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
3336:   if (probAux) {
3337:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3338:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3339:   }
3340:   basisI    = basisField[fieldI];
3341:   basisJ    = basisField[fieldJ];
3342:   basisDerI = basisFieldDer[fieldI];
3343:   basisDerJ = basisFieldDer[fieldJ];
3344:   /* Initialize here in case the function is not defined */
3345:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3346:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3347:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3348:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3349:   for (e = 0; e < Ne; ++e) {
3350:     const PetscReal *v0   = &geom.v0[e*dim];
3351:     const PetscReal *n    = &geom.n[e*dim];
3352:     const PetscReal *J    = &geom.J[e*dim*dim];
3353:     const PetscReal *invJ = &geom.invJ[e*dim*dim];
3354:     const PetscReal  detJ = geom.detJ[e];
3355:     const PetscReal *quadPoints, *quadWeights;
3356:     PetscInt         Nq, q;

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

3362:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3363:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3364:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3365:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3366:       /* TODO: I think I have a mistake in the dim loops here */
3367:       if (g0_func) {
3368:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3369:         g0_func(u, u_t, u_x, a, NULL, a_x, x, n, g0);
3370:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3371:       }
3372:       if (g1_func) {
3373:         PetscInt d, d2;
3374:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3375:         g1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3376:         for (fc = 0; fc < NcI; ++fc) {
3377:           for (gc = 0; gc < NcJ; ++gc) {
3378:             for (d = 0; d < dim; ++d) {
3379:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3380:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3381:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3382:             }
3383:           }
3384:         }
3385:       }
3386:       if (g2_func) {
3387:         PetscInt d, d2;
3388:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3389:         g2_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3390:         for (fc = 0; fc < NcI; ++fc) {
3391:           for (gc = 0; gc < NcJ; ++gc) {
3392:             for (d = 0; d < dim; ++d) {
3393:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3394:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3395:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3396:             }
3397:           }
3398:         }
3399:       }
3400:       if (g3_func) {
3401:         PetscInt d, d2, dp, d3;
3402:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3403:         g3_func(u, u_t, u_x, a, NULL, a_x, x, n, g3);
3404:         for (fc = 0; fc < NcI; ++fc) {
3405:           for (gc = 0; gc < NcJ; ++gc) {
3406:             for (d = 0; d < dim; ++d) {
3407:               for (dp = 0; dp < dim; ++dp) {
3408:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3409:                 for (d2 = 0; d2 < dim; ++d2) {
3410:                   for (d3 = 0; d3 < dim; ++d3) {
3411:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3412:                   }
3413:                 }
3414:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3415:               }
3416:             }
3417:           }
3418:         }
3419:       }

3421:       for (f = 0; f < NbI; ++f) {
3422:         for (fc = 0; fc < NcI; ++fc) {
3423:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3424:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3425:           for (g = 0; g < NbJ; ++g) {
3426:             for (gc = 0; gc < NcJ; ++gc) {
3427:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3428:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3429:               const PetscInt fOff = eOffset+i*totDim+j;
3430:               PetscInt       d, d2;

3432:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3433:               for (d = 0; d < bdim; ++d) {
3434:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
3435:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
3436:                 for (d2 = 0; d2 < bdim; ++d2) {
3437:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
3438:                 }
3439:               }
3440:             }
3441:           }
3442:         }
3443:       }
3444:     }
3445:     if (debug > 1) {
3446:       PetscInt fc, f, gc, g;

3448:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3449:       for (fc = 0; fc < NcI; ++fc) {
3450:         for (f = 0; f < NbI; ++f) {
3451:           const PetscInt i = offsetI + f*NcI+fc;
3452:           for (gc = 0; gc < NcJ; ++gc) {
3453:             for (g = 0; g < NbJ; ++g) {
3454:               const PetscInt j = offsetJ + g*NcJ+gc;
3455:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3456:             }
3457:           }
3458:           PetscPrintf(PETSC_COMM_SELF, "\n");
3459:         }
3460:       }
3461:     }
3462:     cOffset    += totDim;
3463:     cOffsetAux += totDimAux;
3464:     eOffset    += PetscSqr(totDim);
3465:   }
3466:   return(0);
3467: }

3471: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
3472: {
3474:   fem->ops->setfromoptions          = NULL;
3475:   fem->ops->setup                   = PetscFESetUp_Basic;
3476:   fem->ops->view                    = PetscFEView_Basic;
3477:   fem->ops->destroy                 = PetscFEDestroy_Basic;
3478:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
3479:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
3480:   fem->ops->integrate               = PetscFEIntegrate_Basic;
3481:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
3482:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
3483:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
3484:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
3485:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
3486:   return(0);
3487: }

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

3492:   Level: intermediate

3494: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3495: M*/

3499: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
3500: {
3501:   PetscFE_Basic *b;

3506:   PetscNewLog(fem,&b);
3507:   fem->data = b;

3509:   PetscFEInitialize_Basic(fem);
3510:   return(0);
3511: }

3515: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
3516: {
3517:   PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;

3521:   PetscFree(na);
3522:   return(0);
3523: }

3527: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3528:                                                   const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3529: {
3530:   const PetscInt  debug = 0;
3531:   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[]);
3532:   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[]);
3533:   PetscQuadrature quad;
3534:   PetscReal     **basisField, **basisFieldDer;
3535:   PetscScalar    *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3536:   PetscReal      *x;
3537:   PetscInt        dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3538:   PetscErrorCode  ierr;

3541:   PetscFEGetSpatialDimension(fem, &dim);
3542:   PetscFEGetQuadrature(fem, &quad);
3543:   PetscFEGetDimension(fem, &Nb);
3544:   PetscFEGetNumComponents(fem, &Nc);
3545:   PetscDSGetTotalDimension(prob, &totDim);
3546:   PetscDSGetFieldOffset(prob, field, &fOffset);
3547:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3548:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3549:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3550:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3551:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3552:   if (probAux) {
3553:     PetscDSGetTotalDimension(probAux, &totDimAux);
3554:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3555:   }
3556:   for (e = 0; e < Ne; ++e) {
3557:     const PetscReal *quadPoints, *quadWeights;
3558:     PetscInt         Nq, q;

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

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

3584: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3585:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3586: {
3587:   const PetscInt  debug = 0;
3588:   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[]);
3589:   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[]);
3590:   PetscQuadrature quad;
3591:   PetscReal    **basisField, **basisFieldDer;
3592:   PetscScalar    *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3593:   PetscReal      *x;
3594:   PetscInt        dim, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3595:   PetscErrorCode  ierr;

3598:   PetscFEGetSpatialDimension(fem, &dim);
3599:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3600:   PetscFEGetQuadrature(fem, &quad);
3601:   PetscFEGetDimension(fem, &Nb);
3602:   PetscFEGetNumComponents(fem, &Nc);
3603:   PetscDSGetTotalBdDimension(prob, &totDim);
3604:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3605:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3606:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3607:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3608:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3609:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3610:   if (probAux) {
3611:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3612:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3613:   }
3614:   for (e = 0; e < Ne; ++e) {
3615:     const PetscReal *quadPoints, *quadWeights;
3616:     PetscInt         Nq, q;

3618:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3619:     for (q = 0; q < Nq; ++q) {
3620:       const PetscReal *v0   = &geom.v0[(e*Nq+q)*dim];
3621:       const PetscReal *n    = &geom.n[(e*Nq+q)*dim];
3622:       const PetscReal *J    = &geom.J[(e*Nq+q)*dim*dim];
3623:       const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3624:       const PetscReal  detJ = geom.detJ[e*Nq+q];

3626:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3627:       CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3628:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3629:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3630:       f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3631:       f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3632:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3633:     }
3634:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3635:     cOffset    += totDim;
3636:     cOffsetAux += totDimAux;
3637:   }
3638:   return(0);
3639: }

3643: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3644:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3645: {
3646:   const PetscInt  debug      = 0;
3647:   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[]);
3648:   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[]);
3649:   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[]);
3650:   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[]);
3651:   PetscFE         fe;
3652:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3653:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3654:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3655:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3656:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3657:   PetscQuadrature quad;
3658:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3659:   PetscReal      *x;
3660:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3661:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3662:   PetscInt        dim, totDim, totDimAux, e;
3663:   PetscErrorCode  ierr;

3666:   PetscFEGetSpatialDimension(fem, &dim);
3667:   PetscFEGetQuadrature(fem, &quad);
3668:   PetscDSGetTotalDimension(prob, &totDim);
3669:   PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3670:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3671:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3672:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3673:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3674:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3675:   PetscFEGetDimension(fe, &NbI);
3676:   PetscFEGetNumComponents(fe, &NcI);
3677:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3678:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3679:   PetscFEGetDimension(fe, &NbJ);
3680:   PetscFEGetNumComponents(fe, &NcJ);
3681:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3682:   if (probAux) {
3683:     PetscDSGetTotalDimension(probAux, &totDimAux);
3684:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3685:   }
3686:   basisI    = basisField[fieldI];
3687:   basisJ    = basisField[fieldJ];
3688:   basisDerI = basisFieldDer[fieldI];
3689:   basisDerJ = basisFieldDer[fieldJ];
3690:   /* Initialize here in case the function is not defined */
3691:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3692:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3693:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3694:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3695:   for (e = 0; e < Ne; ++e) {
3696:     const PetscReal *quadPoints, *quadWeights;
3697:     PetscInt         Nq, q;

3699:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3700:     for (q = 0; q < Nq; ++q) {
3701:       const PetscReal *v0   = &geom.v0[(e*Nq+q)*dim];
3702:       const PetscReal *J    = &geom.J[(e*Nq+q)*dim*dim];
3703:       const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3704:       const PetscReal  detJ = geom.detJ[e*Nq+q];
3705:       PetscInt         f, g, fc, gc, c;

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

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

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

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

3815: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
3816: {
3818:   fem->ops->setfromoptions          = NULL;
3819:   fem->ops->setup                   = PetscFESetUp_Basic;
3820:   fem->ops->view                    = NULL;
3821:   fem->ops->destroy                 = PetscFEDestroy_Nonaffine;
3822:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
3823:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
3824:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Nonaffine;
3825:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Nonaffine;
3826:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
3827:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Nonaffine;
3828:   return(0);
3829: }

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

3834:   Level: intermediate

3836: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3837: M*/

3841: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
3842: {
3843:   PetscFE_Nonaffine *na;
3844:   PetscErrorCode     ierr;

3848:   PetscNewLog(fem, &na);
3849:   fem->data = na;

3851:   PetscFEInitialize_Nonaffine(fem);
3852:   return(0);
3853: }

3855: #ifdef PETSC_HAVE_OPENCL

3859: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
3860: {
3861:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
3862:   PetscErrorCode  ierr;

3865:   clReleaseCommandQueue(ocl->queue_id);
3866:   ocl->queue_id = 0;
3867:   clReleaseContext(ocl->ctx_id);
3868:   ocl->ctx_id = 0;
3869:   PetscFree(ocl);
3870:   return(0);
3871: }

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

3878: /* dim     Number of spatial dimensions:          2                   */
3879: /* N_b     Number of basis functions:             generated           */
3880: /* N_{bt}  Number of total basis functions:       N_b * N_{comp}      */
3881: /* N_q     Number of quadrature points:           generated           */
3882: /* N_{bs}  Number of block cells                  LCM(N_b, N_q)       */
3883: /* N_{bst} Number of block cell components        LCM(N_{bt}, N_q)    */
3884: /* N_{bl}  Number of concurrent blocks            generated           */
3885: /* N_t     Number of threads:                     N_{bl} * N_{bs}     */
3886: /* N_{cbc} Number of concurrent basis      cells: N_{bl} * N_q        */
3887: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b        */
3888: /* N_{sbc} Number of serial     basis      cells: N_{bs} / N_q        */
3889: /* N_{sqc} Number of serial     quadrature cells: N_{bs} / N_b        */
3890: /* N_{cb}  Number of serial cell batches:         input               */
3891: /* N_c     Number of total cells:                 N_{cb}*N_{t}/N_{comp} */
3892: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
3893: {
3894:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
3895:   PetscQuadrature q;
3896:   char           *string_tail   = *string_buffer;
3897:   char           *end_of_buffer = *string_buffer + buffer_length;
3898:   char            float_str[]   = "float", double_str[]  = "double";
3899:   char           *numeric_str   = &(float_str[0]);
3900:   PetscInt        op            = ocl->op;
3901:   PetscBool       useField      = PETSC_FALSE;
3902:   PetscBool       useFieldDer   = PETSC_TRUE;
3903:   PetscBool       useFieldAux   = useAux;
3904:   PetscBool       useFieldDerAux= PETSC_FALSE;
3905:   PetscBool       useF0         = PETSC_TRUE;
3906:   PetscBool       useF1         = PETSC_TRUE;
3907:   PetscReal      *basis, *basisDer;
3908:   PetscInt        dim, N_b, N_c, N_q, N_t, p, d, b, c;
3909:   size_t          count;
3910:   PetscErrorCode  ierr;

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

4270:   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");}
4271:   if (useF1) {
4272:     switch (dim) {
4273:     case 2:
4274:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4275: "        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"
4276: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4277: "        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"
4278: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n",
4279:                            &count);STRING_ERROR_CHECK("Message to short");break;
4280:     case 3:
4281:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4282: "        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"
4283: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4284: "        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"
4285: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n"
4286: "        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"
4287: "        e_i           += realSpaceDer.z*f_1[fidx].z;\n",
4288:                            &count);STRING_ERROR_CHECK("Message to short");break;
4289:     }
4290:   }
4291:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4292: "      }\n"
4293: "      /* Write element vector for N_{cbc} cells at a time */\n"
4294: "      elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
4295: "    }\n"
4296: "    /* ==== Could do one write per batch ==== */\n"
4297: "  }\n"
4298: "  return;\n"
4299: "}\n",
4300:                        &count);STRING_ERROR_CHECK("Message to short");
4301:   return(0);
4302: }

4306: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
4307: {
4308:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4309:   PetscInt        dim, N_bl;
4310:   PetscBool       flg;
4311:   char           *buffer;
4312:   size_t          len;
4313:   char            errMsg[8192];
4314:   cl_int          ierr2;
4315:   PetscErrorCode  ierr;

4318:   PetscFEGetSpatialDimension(fem, &dim);
4319:   PetscMalloc1(8192, &buffer);
4320:   PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
4321:   PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
4322:   PetscOptionsHasName(fem->hdr.prefix, "-petscfe_opencl_kernel_print", &flg);
4323:   if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
4324:   len  = strlen(buffer);
4325:   *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
4326:   clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
4327:   if (ierr != CL_SUCCESS) {
4328:     clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
4329:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
4330:   }
4331:   PetscFree(buffer);
4332:   *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
4333:   return(0);
4334: }

4338: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
4339: {
4340:   const PetscInt Nblocks = N/blockSize;

4343:   if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
4344:   *z = 1;
4345:   for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
4346:     *y = Nblocks / *x;
4347:     if (*x * *y == Nblocks) break;
4348:   }
4349:   if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
4350:   return(0);
4351: }

4355: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
4356: {
4357:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4358:   PetscStageLog     stageLog;
4359:   PetscEventPerfLog eventLog = NULL;
4360:   PetscInt          stage;
4361:   PetscErrorCode    ierr;

4364:   PetscLogGetStageLog(&stageLog);
4365:   PetscStageLogGetCurrent(stageLog, &stage);
4366:   PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
4367:     /* Log performance info */
4368:   eventLog->eventInfo[ocl->residualEvent].count++;
4369:   eventLog->eventInfo[ocl->residualEvent].time  += time;
4370:   eventLog->eventInfo[ocl->residualEvent].flops += flops;
4371:   return(0);
4372: }

4376: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
4377:                                                const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4378: {
4379:   /* Nbc = batchSize */
4380:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4381:   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[]);
4382:   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[]);
4383:   PetscQuadrature   q;
4384:   PetscInt          dim;
4385:   PetscInt          N_b;    /* The number of basis functions */
4386:   PetscInt          N_comp; /* The number of basis function components */
4387:   PetscInt          N_bt;   /* The total number of scalar basis functions */
4388:   PetscInt          N_q;    /* The number of quadrature points */
4389:   PetscInt          N_bst;  /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
4390:   PetscInt          N_t;    /* The number of threads, N_bst * N_bl */
4391:   PetscInt          N_bl;   /* The number of blocks */
4392:   PetscInt          N_bc;   /* The batch size, N_bl*N_q*N_b */
4393:   PetscInt          N_cb;   /* The number of batches */
4394:   PetscInt          numFlops, f0Flops = 0, f1Flops = 0;
4395:   PetscBool         useAux      = coefficientsAux ? PETSC_TRUE : PETSC_FALSE;
4396:   PetscBool         useField    = PETSC_FALSE;
4397:   PetscBool         useFieldDer = PETSC_TRUE;
4398:   PetscBool         useF0       = PETSC_TRUE;
4399:   PetscBool         useF1       = PETSC_TRUE;
4400:   /* OpenCL variables */
4401:   cl_program        ocl_prog;
4402:   cl_kernel         ocl_kernel;
4403:   cl_event          ocl_ev;         /* The event for tracking kernel execution */
4404:   cl_ulong          ns_start;       /* Nanoseconds counter on GPU at kernel start */
4405:   cl_ulong          ns_end;         /* Nanoseconds counter on GPU at kernel stop */
4406:   cl_mem            o_jacobianInverses, o_jacobianDeterminants;
4407:   cl_mem            o_coefficients, o_coefficientsAux, o_elemVec;
4408:   float            *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
4409:   double           *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
4410:   void             *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
4411:   size_t            local_work_size[3], global_work_size[3];
4412:   size_t            realSize, x, y, z;
4413:   PetscErrorCode    ierr;

4416:   if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
4417:   PetscFEGetSpatialDimension(fem, &dim);
4418:   PetscFEGetQuadrature(fem, &q);
4419:   PetscFEGetDimension(fem, &N_b);
4420:   PetscFEGetNumComponents(fem, &N_comp);
4421:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4422:   PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
4423:   N_bt  = N_b*N_comp;
4424:   N_q   = q->numPoints;
4425:   N_bst = N_bt*N_q;
4426:   N_t   = N_bst*N_bl;
4427:   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);
4428:   /* Calculate layout */
4429:   if (Ne % (N_cb*N_bc)) { /* Remainder cells */
4430:     PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);
4431:     return(0);
4432:   }
4433:   PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
4434:   local_work_size[0]  = N_bc*N_comp;
4435:   local_work_size[1]  = 1;
4436:   local_work_size[2]  = 1;
4437:   global_work_size[0] = x * local_work_size[0];
4438:   global_work_size[1] = y * local_work_size[1];
4439:   global_work_size[2] = z * local_work_size[2];
4440:   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);
4441:   PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
4442:   /* Generate code */
4443:   if (probAux) {
4444:     PetscSpace P;
4445:     PetscInt   NfAux, order, f;

4447:     PetscDSGetNumFields(probAux, &NfAux);
4448:     for (f = 0; f < NfAux; ++f) {
4449:       PetscFE feAux;

4451:       PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
4452:       PetscFEGetBasisSpace(feAux, &P);
4453:       PetscSpaceGetOrder(P, &order);
4454:       if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
4455:     }
4456:   }
4457:   PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
4458:   /* Create buffers on the device and send data over */
4459:   PetscDataTypeGetSize(ocl->realType, &realSize);
4460:   if (sizeof(PetscReal) != realSize) {
4461:     switch (ocl->realType) {
4462:     case PETSC_FLOAT:
4463:     {
4464:       PetscInt c, b, d;

4466:       PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
4467:       for (c = 0; c < Ne; ++c) {
4468:         f_detJ[c] = (float) geom.detJ[c];
4469:         for (d = 0; d < dim*dim; ++d) {
4470:           f_invJ[c*dim*dim+d] = (float) geom.invJ[c*dim*dim+d];
4471:         }
4472:         for (b = 0; b < N_bt; ++b) {
4473:           f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
4474:         }
4475:       }
4476:       if (coefficientsAux) { /* Assume P0 */
4477:         for (c = 0; c < Ne; ++c) {
4478:           f_coeffAux[c] = (float) coefficientsAux[c];
4479:         }
4480:       }
4481:       oclCoeff      = (void *) f_coeff;
4482:       if (coefficientsAux) {
4483:         oclCoeffAux = (void *) f_coeffAux;
4484:       } else {
4485:         oclCoeffAux = NULL;
4486:       }
4487:       oclInvJ       = (void *) f_invJ;
4488:       oclDetJ       = (void *) f_detJ;
4489:     }
4490:     break;
4491:     case PETSC_DOUBLE:
4492:     {
4493:       PetscInt c, b, d;

4495:       PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
4496:       for (c = 0; c < Ne; ++c) {
4497:         d_detJ[c] = (double) geom.detJ[c];
4498:         for (d = 0; d < dim*dim; ++d) {
4499:           d_invJ[c*dim*dim+d] = (double) geom.invJ[c*dim*dim+d];
4500:         }
4501:         for (b = 0; b < N_bt; ++b) {
4502:           d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
4503:         }
4504:       }
4505:       if (coefficientsAux) { /* Assume P0 */
4506:         for (c = 0; c < Ne; ++c) {
4507:           d_coeffAux[c] = (double) coefficientsAux[c];
4508:         }
4509:       }
4510:       oclCoeff      = (void *) d_coeff;
4511:       if (coefficientsAux) {
4512:         oclCoeffAux = (void *) d_coeffAux;
4513:       } else {
4514:         oclCoeffAux = NULL;
4515:       }
4516:       oclInvJ       = (void *) d_invJ;
4517:       oclDetJ       = (void *) d_detJ;
4518:     }
4519:     break;
4520:     default:
4521:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4522:     }
4523:   } else {
4524:     oclCoeff    = (void *) coefficients;
4525:     oclCoeffAux = (void *) coefficientsAux;
4526:     oclInvJ     = (void *) geom.invJ;
4527:     oclDetJ     = (void *) geom.detJ;
4528:   }
4529:   o_coefficients         = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt    * realSize, oclCoeff,    &ierr);
4530:   if (coefficientsAux) {
4531:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclCoeffAux, &ierr);
4532:   } else {
4533:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY,                        Ne         * realSize, oclCoeffAux, &ierr);
4534:   }
4535:   o_jacobianInverses     = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ,     &ierr);
4536:   o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclDetJ,     &ierr);
4537:   o_elemVec              = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY,                       Ne*N_bt    * realSize, NULL,        &ierr);
4538:   /* Kernel launch */
4539:   clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
4540:   clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
4541:   clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
4542:   clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
4543:   clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
4544:   clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
4545:   clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
4546:   /* Read data back from device */
4547:   if (sizeof(PetscReal) != realSize) {
4548:     switch (ocl->realType) {
4549:     case PETSC_FLOAT:
4550:     {
4551:       float   *elem;
4552:       PetscInt c, b;

4554:       PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
4555:       PetscMalloc1(Ne*N_bt, &elem);
4556:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4557:       for (c = 0; c < Ne; ++c) {
4558:         for (b = 0; b < N_bt; ++b) {
4559:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4560:         }
4561:       }
4562:       PetscFree(elem);
4563:     }
4564:     break;
4565:     case PETSC_DOUBLE:
4566:     {
4567:       double  *elem;
4568:       PetscInt c, b;

4570:       PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
4571:       PetscMalloc1(Ne*N_bt, &elem);
4572:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4573:       for (c = 0; c < Ne; ++c) {
4574:         for (b = 0; b < N_bt; ++b) {
4575:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4576:         }
4577:       }
4578:       PetscFree(elem);
4579:     }
4580:     break;
4581:     default:
4582:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4583:     }
4584:   } else {
4585:     clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
4586:   }
4587:   /* Log performance */
4588:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
4589:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END,   sizeof(cl_ulong), &ns_end,   NULL);
4590:   f0Flops = 0;
4591:   switch (ocl->op) {
4592:   case LAPLACIAN:
4593:     f1Flops = useAux ? dim : 0;break;
4594:   case ELASTICITY:
4595:     f1Flops = 2*dim*dim;break;
4596:   }
4597:   numFlops = Ne*(
4598:     N_q*(
4599:       N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
4600:       /*+
4601:        N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
4602:       +
4603:       N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
4604:     +
4605:     N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
4606:   PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
4607:   /* Cleanup */
4608:   clReleaseMemObject(o_coefficients);
4609:   clReleaseMemObject(o_coefficientsAux);
4610:   clReleaseMemObject(o_jacobianInverses);
4611:   clReleaseMemObject(o_jacobianDeterminants);
4612:   clReleaseMemObject(o_elemVec);
4613:   clReleaseKernel(ocl_kernel);
4614:   clReleaseProgram(ocl_prog);
4615:   return(0);
4616: }

4620: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
4621: {
4623:   fem->ops->setfromoptions          = NULL;
4624:   fem->ops->setup                   = PetscFESetUp_Basic;
4625:   fem->ops->view                    = NULL;
4626:   fem->ops->destroy                 = PetscFEDestroy_OpenCL;
4627:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4628:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
4629:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_OpenCL;
4630:   fem->ops->integratebdresidual     = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
4631:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
4632:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
4633:   return(0);
4634: }

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

4639:   Level: intermediate

4641: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4642: M*/

4646: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
4647: {
4648:   PetscFE_OpenCL *ocl;
4649:   cl_uint         num_platforms;
4650:   cl_platform_id  platform_ids[42];
4651:   cl_uint         num_devices;
4652:   cl_device_id    device_ids[42];
4653:   cl_int          ierr2;
4654:   PetscErrorCode  ierr;

4658:   PetscNewLog(fem,&ocl);
4659:   fem->data = ocl;

4661:   /* Init Platform */
4662:   clGetPlatformIDs(42, platform_ids, &num_platforms);
4663:   if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
4664:   ocl->pf_id = platform_ids[0];
4665:   /* Init Device */
4666:   clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
4667:   if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
4668:   ocl->dev_id = device_ids[0];
4669:   /* Create context with one command queue */
4670:   ocl->ctx_id   = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
4671:   ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
4672:   /* Types */
4673:   ocl->realType = PETSC_FLOAT;
4674:   /* Register events */
4675:   PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
4676:   /* Equation handling */
4677:   ocl->op = LAPLACIAN;

4679:   PetscFEInitialize_OpenCL(fem);
4680:   return(0);
4681: }

4685: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
4686: {
4687:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

4691:   ocl->realType = realType;
4692:   return(0);
4693: }

4697: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
4698: {
4699:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

4704:   *realType = ocl->realType;
4705:   return(0);
4706: }

4708: #endif /* PETSC_HAVE_OPENCL */

4712: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
4713: {
4714:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4715:   PetscErrorCode     ierr;

4718:   CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
4719:   PetscFree(cmp->embedding);
4720:   PetscFree(cmp);
4721:   return(0);
4722: }

4726: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
4727: {
4728:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4729:   DM                 K;
4730:   PetscReal         *work, *subpoint;
4731:   PetscBLASInt      *pivots;
4732: #ifndef PETSC_USE_COMPLEX
4733:   PetscBLASInt       n, info;
4734: #endif
4735:   PetscInt           dim, pdim, spdim, j, s;
4736:   PetscErrorCode     ierr;

4739:   /* Get affine mapping from reference cell to each subcell */
4740:   PetscDualSpaceGetDM(fem->dualSpace, &K);
4741:   DMPlexGetDimension(K, &dim);
4742:   DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
4743:   CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
4744:   /* Determine dof embedding into subelements */
4745:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
4746:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
4747:   PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
4748:   DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
4749:   for (s = 0; s < cmp->numSubelements; ++s) {
4750:     PetscInt sd = 0;

4752:     for (j = 0; j < pdim; ++j) {
4753:       PetscBool       inside;
4754:       PetscQuadrature f;
4755:       PetscInt        d, e;

4757:       PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
4758:       /* Apply transform to first point, and check that point is inside subcell */
4759:       for (d = 0; d < dim; ++d) {
4760:         subpoint[d] = -1.0;
4761:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
4762:       }
4763:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
4764:       if (inside) {cmp->embedding[s*spdim+sd++] = j;}
4765:     }
4766:     if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
4767:   }
4768:   DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
4769:   /* Construct the change of basis from prime basis to nodal basis for each subelement */
4770:   PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
4771:   PetscMalloc2(spdim,&pivots,spdim,&work);
4772:   for (s = 0; s < cmp->numSubelements; ++s) {
4773:     for (j = 0; j < spdim; ++j) {
4774:       PetscReal      *Bf;
4775:       PetscQuadrature f;
4776:       PetscInt        q, k;

4778:       PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
4779:       PetscMalloc1(f->numPoints*spdim,&Bf);
4780:       PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
4781:       for (k = 0; k < spdim; ++k) {
4782:         /* n_j \cdot \phi_k */
4783:         fem->invV[(s*spdim + j)*spdim+k] = 0.0;
4784:         for (q = 0; q < f->numPoints; ++q) {
4785:           fem->invV[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
4786:         }
4787:       }
4788:       PetscFree(Bf);
4789:     }
4790: #ifndef PETSC_USE_COMPLEX
4791:     n = spdim;
4792:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &fem->invV[s*spdim*spdim], &n, pivots, &info));
4793:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &fem->invV[s*spdim*spdim], &n, pivots, work, &n, &info));
4794: #endif
4795:   }
4796:   PetscFree2(pivots,work);
4797:   return(0);
4798: }

4802: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
4803: {
4804:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4805:   DM                 dm;
4806:   PetscInt           pdim;  /* Dimension of FE space P */
4807:   PetscInt           spdim; /* Dimension of subelement FE space P */
4808:   PetscInt           dim;   /* Spatial dimension */
4809:   PetscInt           comp;  /* Field components */
4810:   PetscInt          *subpoints;
4811:   PetscReal         *tmpB, *tmpD, *tmpH, *subpoint;
4812:   PetscInt           p, s, d, e, j, k;
4813:   PetscErrorCode     ierr;

4816:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
4817:   DMPlexGetDimension(dm, &dim);
4818:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
4819:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
4820:   PetscFEGetNumComponents(fem, &comp);
4821:   /* Divide points into subelements */
4822:   DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
4823:   DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
4824:   for (p = 0; p < npoints; ++p) {
4825:     for (s = 0; s < cmp->numSubelements; ++s) {
4826:       PetscBool inside;

4828:       /* Apply transform, and check that point is inside cell */
4829:       for (d = 0; d < dim; ++d) {
4830:         subpoint[d] = -1.0;
4831:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
4832:       }
4833:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
4834:       if (inside) {subpoints[p] = s; break;}
4835:     }
4836:     if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
4837:   }
4838:   DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
4839:   /* Evaluate the prime basis functions at all points */
4840:   if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
4841:   if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
4842:   if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
4843:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
4844:   /* Translate to the nodal basis */
4845:   if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
4846:   if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
4847:   if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
4848:   for (p = 0; p < npoints; ++p) {
4849:     const PetscInt s = subpoints[p];

4851:     if (B) {
4852:       /* Multiply by V^{-1} (spdim x spdim) */
4853:       for (j = 0; j < spdim; ++j) {
4854:         const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
4855:         PetscInt       c;

4857:         B[i] = 0.0;
4858:         for (k = 0; k < spdim; ++k) {
4859:           B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
4860:         }
4861:         for (c = 1; c < comp; ++c) {
4862:           B[i+c] = B[i];
4863:         }
4864:       }
4865:     }
4866:     if (D) {
4867:       /* Multiply by V^{-1} (spdim x spdim) */
4868:       for (j = 0; j < spdim; ++j) {
4869:         for (d = 0; d < dim; ++d) {
4870:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
4871:           PetscInt       c;

4873:           D[i] = 0.0;
4874:           for (k = 0; k < spdim; ++k) {
4875:             D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
4876:           }
4877:           for (c = 1; c < comp; ++c) {
4878:             D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
4879:           }
4880:         }
4881:       }
4882:     }
4883:     if (H) {
4884:       /* Multiply by V^{-1} (pdim x pdim) */
4885:       for (j = 0; j < spdim; ++j) {
4886:         for (d = 0; d < dim*dim; ++d) {
4887:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
4888:           PetscInt       c;

4890:           H[i] = 0.0;
4891:           for (k = 0; k < spdim; ++k) {
4892:             H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
4893:           }
4894:           for (c = 1; c < comp; ++c) {
4895:             H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
4896:           }
4897:         }
4898:       }
4899:     }
4900:   }
4901:   DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
4902:   if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
4903:   if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
4904:   if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
4905:   return(0);
4906: }

4910: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
4911: {
4913:   fem->ops->setfromoptions          = NULL;
4914:   fem->ops->setup                   = PetscFESetUp_Composite;
4915:   fem->ops->view                    = NULL;
4916:   fem->ops->destroy                 = PetscFEDestroy_Composite;
4917:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4918:   fem->ops->gettabulation           = PetscFEGetTabulation_Composite;
4919:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
4920:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
4921:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
4922:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
4923:   return(0);
4924: }

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

4929:   Level: intermediate

4931: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4932: M*/

4936: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
4937: {
4938:   PetscFE_Composite *cmp;
4939:   PetscErrorCode     ierr;

4943:   PetscNewLog(fem, &cmp);
4944:   fem->data = cmp;

4946:   cmp->cellRefiner    = 0;
4947:   cmp->numSubelements = -1;
4948:   cmp->v0             = NULL;
4949:   cmp->jac            = NULL;

4951:   PetscFEInitialize_Composite(fem);
4952:   return(0);
4953: }

4957: PetscErrorCode PetscFECompositeExpandQuadrature(PetscFE fem, PetscQuadrature q, PetscQuadrature *qref)
4958: {
4959:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4960:   const PetscReal   *points,    *weights;
4961:   PetscReal         *pointsRef, *weightsRef;
4962:   PetscInt           dim, npoints, npointsRef, c, p, d, e;
4963:   PetscErrorCode     ierr;

4969:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
4970:   PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
4971:   npointsRef = npoints*cmp->numSubelements;
4972:   PetscMalloc1(npointsRef*dim,&pointsRef);
4973:   PetscMalloc1(npointsRef,&weightsRef);
4974:   for (c = 0; c < cmp->numSubelements; ++c) {
4975:     for (p = 0; p < npoints; ++p) {
4976:       for (d = 0; d < dim; ++d) {
4977:         pointsRef[(c*npoints + p)*dim+d] = cmp->v0[c*dim+d];
4978:         for (e = 0; e < dim; ++e) {
4979:           pointsRef[(c*npoints + p)*dim+d] += cmp->jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
4980:         }
4981:       }
4982:       /* Could also use detJ here */
4983:       weightsRef[c*npoints+p] = weights[p]/cmp->numSubelements;
4984:     }
4985:   }
4986:   PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
4987:   return(0);
4988: }

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

4995:   Not collective

4997:   Input Parameter:
4998: . fe - The PetscFE

5000:   Output Parameter:
5001: . dim - The dimension

5003:   Level: intermediate

5005: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
5006: @*/
5007: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
5008: {

5014:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
5015:   return(0);
5016: }

5018: /*
5019: Purpose: Compute element vector for chunk of elements

5021: Input:
5022:   Sizes:
5023:      Ne:  number of elements
5024:      Nf:  number of fields
5025:      PetscFE
5026:        dim: spatial dimension
5027:        Nb:  number of basis functions
5028:        Nc:  number of field components
5029:        PetscQuadrature
5030:          Nq:  number of quadrature points

5032:   Geometry:
5033:      PetscCellGeometry
5034:        PetscReal v0s[Ne*dim]
5035:        PetscReal jacobians[Ne*dim*dim]        possibly *Nq
5036:        PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5037:        PetscReal jacobianDeterminants[Ne]     possibly *Nq
5038:   FEM:
5039:      PetscFE
5040:        PetscQuadrature
5041:          PetscReal   quadPoints[Nq*dim]
5042:          PetscReal   quadWeights[Nq]
5043:        PetscReal   basis[Nq*Nb*Nc]
5044:        PetscReal   basisDer[Nq*Nb*Nc*dim]
5045:      PetscScalar coefficients[Ne*Nb*Nc]
5046:      PetscScalar elemVec[Ne*Nb*Nc]

5048:   Problem:
5049:      PetscInt f: the active field
5050:      f0, f1

5052:   Work Space:
5053:      PetscFE
5054:        PetscScalar f0[Nq*dim];
5055:        PetscScalar f1[Nq*dim*dim];
5056:        PetscScalar u[Nc];
5057:        PetscScalar gradU[Nc*dim];
5058:        PetscReal   x[dim];
5059:        PetscScalar realSpaceDer[dim];

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

5063: Input:
5064:   Sizes:
5065:      N_cb: Number of serial cell batches

5067:   Geometry:
5068:      PetscReal v0s[Ne*dim]
5069:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
5070:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5071:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
5072:   FEM:
5073:      static PetscReal   quadPoints[Nq*dim]
5074:      static PetscReal   quadWeights[Nq]
5075:      static PetscReal   basis[Nq*Nb*Nc]
5076:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
5077:      PetscScalar coefficients[Ne*Nb*Nc]
5078:      PetscScalar elemVec[Ne*Nb*Nc]

5080: ex62.c:
5081:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
5082:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
5083:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
5084:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

5086: ex52.c:
5087:   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)
5088:   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)

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

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

5097: ex52_integrateElementOpenCL.c:
5098: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
5099:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5100:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

5110:   Not collective

5112:   Input Parameters:
5113: + fem          - The PetscFE object for the field being integrated
5114: . prob         - The PetscDS specifing the discretizations and continuum functions
5115: . field        - The field being integrated
5116: . Ne           - The number of elements in the chunk
5117: . geom         - The cell geometry for each cell in the chunk
5118: . coefficients - The array of FEM basis coefficients for the elements
5119: . probAux      - The PetscDS specifing the auxiliary discretizations
5120: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5122:   Output Parameter
5123: . integral     - the integral for this field

5125:   Level: developer

5127: .seealso: PetscFEIntegrateResidual()
5128: @*/
5129: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5130:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
5131: {

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

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

5146:   Not collective

5148:   Input Parameters:
5149: + fem          - The PetscFE object for the field being integrated
5150: . prob         - The PetscDS specifing the discretizations and continuum functions
5151: . field        - The field being integrated
5152: . Ne           - The number of elements in the chunk
5153: . geom         - The cell geometry for each cell in the chunk
5154: . coefficients - The array of FEM basis coefficients for the elements
5155: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5156: . probAux      - The PetscDS specifing the auxiliary discretizations
5157: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5159:   Output Parameter
5160: . elemVec      - the element residual vectors from each element

5162:   Note:
5163: $ Loop over batch of elements (e):
5164: $   Loop over quadrature points (q):
5165: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
5166: $     Call f_0 and f_1
5167: $   Loop over element vector entries (f,fc --> i):
5168: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

5170:   Level: developer

5172: .seealso: PetscFEIntegrateResidual()
5173: @*/
5174: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5175:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5176: {

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

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

5191:   Not collective

5193:   Input Parameters:
5194: + fem          - The PetscFE object for the field being integrated
5195: . prob         - The PetscDS specifing the discretizations and continuum functions
5196: . field        - The field being integrated
5197: . Ne           - The number of elements in the chunk
5198: . geom         - The cell geometry for each cell in the chunk
5199: . coefficients - The array of FEM basis coefficients for the elements
5200: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5201: . probAux      - The PetscDS specifing the auxiliary discretizations
5202: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5204:   Output Parameter
5205: . elemVec      - the element residual vectors from each element

5207:   Level: developer

5209: .seealso: PetscFEIntegrateResidual()
5210: @*/
5211: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5212:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5213: {

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

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

5227:   Not collective

5229:   Input Parameters:
5230: + fem          = The PetscFE object for the field being integrated
5231: . prob         - The PetscDS specifing the discretizations and continuum functions
5232: . fieldI       - The test field being integrated
5233: . fieldJ       - The basis field being integrated
5234: . Ne           - The number of elements in the chunk
5235: . geom         - The cell geometry for each cell in the chunk
5236: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5237: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5238: . probAux      - The PetscDS specifing the auxiliary discretizations
5239: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5241:   Output Parameter
5242: . elemMat      - the element matrices for the Jacobian from each element

5244:   Note:
5245: $ Loop over batch of elements (e):
5246: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5247: $     Loop over quadrature points (q):
5248: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5249: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5250: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5251: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5252: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5253: */
5254: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
5255:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5256: {

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

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

5270:   Not collective

5272:   Input Parameters:
5273: + fem          = The PetscFE object for the field being integrated
5274: . prob         - The PetscDS specifing the discretizations and continuum functions
5275: . fieldI       - The test field being integrated
5276: . fieldJ       - The basis field being integrated
5277: . Ne           - The number of elements in the chunk
5278: . geom         - The cell geometry for each cell in the chunk
5279: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5280: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5281: . probAux      - The PetscDS specifing the auxiliary discretizations
5282: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5284:   Output Parameter
5285: . elemMat              - the element matrices for the Jacobian from each element

5287:   Note:
5288: $ Loop over batch of elements (e):
5289: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5290: $     Loop over quadrature points (q):
5291: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5292: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5293: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5294: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5295: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5296: */
5297: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
5298:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5299: {

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

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

5315:   Input Parameter:
5316: . fe - The initial PetscFE

5318:   Output Parameter:
5319: . feRef - The refined PetscFE

5321:   Level: developer

5323: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5324: @*/
5325: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
5326: {
5327:   PetscSpace       P, Pref;
5328:   PetscDualSpace   Q, Qref;
5329:   DM               K, Kref;
5330:   PetscQuadrature  q, qref;
5331:   PetscInt         numComp;
5332:   PetscErrorCode   ierr;

5335:   PetscFEGetBasisSpace(fe, &P);
5336:   PetscFEGetDualSpace(fe, &Q);
5337:   PetscFEGetQuadrature(fe, &q);
5338:   PetscDualSpaceGetDM(Q, &K);
5339:   /* Create space */
5340:   PetscObjectReference((PetscObject) P);
5341:   Pref = P;
5342:   /* Create dual space */
5343:   PetscDualSpaceDuplicate(Q, &Qref);
5344:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
5345:   PetscDualSpaceSetDM(Qref, Kref);
5346:   DMDestroy(&Kref);
5347:   PetscDualSpaceSetUp(Qref);
5348:   /* Create element */
5349:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
5350:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
5351:   PetscFESetBasisSpace(*feRef, Pref);
5352:   PetscFESetDualSpace(*feRef, Qref);
5353:   PetscFEGetNumComponents(fe,    &numComp);
5354:   PetscFESetNumComponents(*feRef, numComp);
5355:   PetscFESetUp(*feRef);
5356:   PetscSpaceDestroy(&Pref);
5357:   PetscDualSpaceDestroy(&Qref);
5358:   /* Create quadrature */
5359:   PetscFECompositeExpandQuadrature(*feRef, q, &qref);
5360:   PetscFESetQuadrature(*feRef, qref);
5361:   PetscQuadratureDestroy(&qref);
5362:   return(0);
5363: }

5367: /*@
5368:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

5370:   Collective on DM

5372:   Input Parameters:
5373: + dm         - The underlying DM for the domain
5374: . dim        - The spatial dimension
5375: . numComp    - The number of components
5376: . isSimplex  - Flag for simplex reference cell, otherwise its a tensor product
5377: . prefix     - The options prefix, or NULL
5378: - qorder     - The quadrature order

5380:   Output Parameter:
5381: . fem - The PetscFE object

5383:   Level: beginner

5385: .keywords: PetscFE, finite element
5386: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
5387: @*/
5388: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
5389: {
5390:   PetscQuadrature q;
5391:   DM              K;
5392:   PetscSpace      P;
5393:   PetscDualSpace  Q;
5394:   PetscInt        order;
5395:   PetscErrorCode  ierr;

5398:   /* Create space */
5399:   PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
5400:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
5401:   PetscSpaceSetFromOptions(P);
5402:   PetscSpacePolynomialSetNumVariables(P, dim);
5403:   PetscSpaceSetUp(P);
5404:   PetscSpaceGetOrder(P, &order);
5405:   /* Create dual space */
5406:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
5407:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
5408:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
5409:   PetscDualSpaceSetDM(Q, K);
5410:   DMDestroy(&K);
5411:   PetscDualSpaceSetOrder(Q, order);
5412:   PetscDualSpaceSetFromOptions(Q);
5413:   PetscDualSpaceSetUp(Q);
5414:   /* Create element */
5415:   PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
5416:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
5417:   PetscFESetFromOptions(*fem);
5418:   PetscFESetBasisSpace(*fem, P);
5419:   PetscFESetDualSpace(*fem, Q);
5420:   PetscFESetNumComponents(*fem, numComp);
5421:   PetscFESetUp(*fem);
5422:   PetscSpaceDestroy(&P);
5423:   PetscDualSpaceDestroy(&Q);
5424:   /* Create quadrature (with specified order if given) */
5425:   if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5426:   else           {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5427:   PetscFESetQuadrature(*fem, q);
5428:   PetscQuadratureDestroy(&q);
5429:   return(0);
5430: }