Actual source code: dtfe.c

petsc-master 2016-08-30
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>
 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:   PetscSpaceRegisterAll();
135:   PetscFunctionListFind(PetscSpaceList, name, &r);
136:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSpace type: %s", name);

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

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

152:   Not Collective

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

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

160:   Level: intermediate

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

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

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

184:   Collective on PetscSpace

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

190:   Level: developer

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

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

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

214:   Collective on PetscSpace

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

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

222:   Level: developer

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

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

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

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

265:   Collective on PetscSpace

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

270:   Level: developer

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

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

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

289:   Collective on PetscSpace

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

294:   Level: developer

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

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

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

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

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

320:   Collective on MPI_Comm

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

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

328:   Level: beginner

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

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

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

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

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

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

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

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

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

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

378:   Level: intermediate

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

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

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

400:   Level: intermediate

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

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

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

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

427:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

519:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);
520:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", NULL);
521:   PetscFree(poly->degrees);
522:   PetscFree(poly);
523:   return(0);
524: }

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

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

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

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

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

562:   Level: developer

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

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

590: /*
591:   LatticePointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to at most 'max'.
592:                                        Ordering is lexicographic with lowest index as least significant in ordering.
593:                                        e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,0}.

595:   Input Parameters:
596: + len - The length of the tuple
597: . max - The maximum sum
598: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition

600:   Output Parameter:
601: . tup - A tuple of len integers whos sum is at most 'max'
602: */
603: static PetscErrorCode LatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
604: {
606:   while (len--) {
607:     max -= tup[len];
608:     if (!max) {
609:       tup[len] = 0;
610:       break;
611:     }
612:   }
613:   tup[++len]++;
614:   return(0);
615: }

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

622:   Input Parameters:
623: + len - The length of the tuple
624: . max - The max for all entries in the tuple
625: - ind - The current multi-index of the tuple, initialized to the 0 tuple

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

631:   Level: developer

633: .seealso: 
634: */
635: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
636: {
637:   PetscInt       i;

641:   if (len == 1) {
642:     tup[0] = ind[0]++;
643:     ind[0] = ind[0] >= max ? -1 : ind[0];
644:   } else if (max == 0) {
645:     for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
646:   } else {
647:     tup[0] = ind[0];
648:     TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
649:     if (ind[1] < 0) {
650:       ind[1] = 0;
651:       if (ind[0] == max-1) {ind[0] = -1;}
652:       else                 {++ind[0];}
653:     }
654:   }
655:   return(0);
656: }

660: /*
661:   TensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
662:                                       Ordering is lexicographic with lowest index as least significant in ordering.
663:                                       e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,1}, {0,2}, {1,2}, {2,2}.

665:   Input Parameters:
666: + len - The length of the tuple
667: . max - The maximum value
668: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition

670:   Output Parameter:
671: . tup - A tuple of len integers whos sum is at most 'max'
672: */
673: static PetscErrorCode TensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
674: {
675:   PetscInt       i;

678:   for (i = 0; i < len; i++) {
679:     if (tup[i] < max) {
680:       break;
681:     } else {
682:       tup[i] = 0;
683:     }
684:   }
685:   tup[i]++;
686:   return(0);
687: }

691: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
692: {
693:   PetscSpace_Poly *poly    = (PetscSpace_Poly *) sp->data;
694:   DM               dm      = sp->dm;
695:   PetscInt         ndegree = sp->order+1;
696:   PetscInt        *degrees = poly->degrees;
697:   PetscInt         dim     = poly->numVariables;
698:   PetscReal       *lpoints, *tmp, *LB, *LD, *LH;
699:   PetscInt        *ind, *tup;
700:   PetscInt         pdim, d, der, i, p, deg, o;
701:   PetscErrorCode   ierr;

704:   PetscSpaceGetDimension(sp, &pdim);
705:   DMGetWorkArray(dm, npoints, PETSC_REAL, &lpoints);
706:   DMGetWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
707:   if (B) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
708:   if (D) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
709:   if (H) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
710:   for (d = 0; d < dim; ++d) {
711:     for (p = 0; p < npoints; ++p) {
712:       lpoints[p] = points[p*dim+d];
713:     }
714:     PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
715:     /* LB, LD, LH (ndegree * dim x npoints) */
716:     for (deg = 0; deg < ndegree; ++deg) {
717:       for (p = 0; p < npoints; ++p) {
718:         if (B) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
719:         if (D) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
720:         if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
721:       }
722:     }
723:   }
724:   /* Multiply by A (pdim x ndegree * dim) */
725:   PetscMalloc2(dim,&ind,dim,&tup);
726:   if (B) {
727:     /* B (npoints x pdim) */
728:     if (poly->tensor) {
729:       i = 0;
730:       PetscMemzero(ind, dim * sizeof(PetscInt));
731:       while (ind[0] >= 0) {
732:         TensorPoint_Internal(dim, sp->order+1, ind, tup);
733:         for (p = 0; p < npoints; ++p) {
734:           B[p*pdim + i] = 1.0;
735:           for (d = 0; d < dim; ++d) {
736:             B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
737:           }
738:         }
739:         ++i;
740:       }
741:     } else {
742:       i = 0;
743:       for (o = 0; o <= sp->order; ++o) {
744:         PetscMemzero(ind, dim * sizeof(PetscInt));
745:         while (ind[0] >= 0) {
746:           LatticePoint_Internal(dim, o, ind, tup);
747:           for (p = 0; p < npoints; ++p) {
748:             B[p*pdim + i] = 1.0;
749:             for (d = 0; d < dim; ++d) {
750:               B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
751:             }
752:           }
753:           ++i;
754:         }
755:       }
756:     }
757:   }
758:   if (D) {
759:     /* D (npoints x pdim x dim) */
760:     if (poly->tensor) {
761:       i = 0;
762:       PetscMemzero(ind, dim * sizeof(PetscInt));
763:       while (ind[0] >= 0) {
764:         TensorPoint_Internal(dim, sp->order+1, ind, tup);
765:         for (p = 0; p < npoints; ++p) {
766:           for (der = 0; der < dim; ++der) {
767:             D[(p*pdim + i)*dim + der] = 1.0;
768:             for (d = 0; d < dim; ++d) {
769:               if (d == der) {
770:                 D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
771:               } else {
772:                 D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
773:               }
774:             }
775:           }
776:         }
777:         ++i;
778:       }
779:     } else {
780:       i = 0;
781:       for (o = 0; o <= sp->order; ++o) {
782:         PetscMemzero(ind, dim * sizeof(PetscInt));
783:         while (ind[0] >= 0) {
784:           LatticePoint_Internal(dim, o, ind, tup);
785:           for (p = 0; p < npoints; ++p) {
786:             for (der = 0; der < dim; ++der) {
787:               D[(p*pdim + i)*dim + der] = 1.0;
788:               for (d = 0; d < dim; ++d) {
789:                 if (d == der) {
790:                   D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
791:                 } else {
792:                   D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
793:                 }
794:               }
795:             }
796:           }
797:           ++i;
798:         }
799:       }
800:     }
801:   }
802:   if (H) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to code second derivatives");
803:   PetscFree2(ind,tup);
804:   if (B) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
805:   if (D) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
806:   if (H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
807:   DMRestoreWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
808:   DMRestoreWorkArray(dm, npoints, PETSC_REAL, &lpoints);
809:   return(0);
810: }

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

819:   Input Parameters:
820: + sp     - the function space object
821: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

823:   Level: beginner

825: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetOrder(), PetscSpacePolynomialSetNumVariables()
826: @*/
827: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
828: {

833:   PetscTryMethod(sp,"PetscSpacePolynomialSetTensor_C",(PetscSpace,PetscBool),(sp,tensor));
834:   return(0);
835: }

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

844:   Input Parameters:
845: . sp     - the function space object

847:   Output Parameters:
848: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space

850:   Level: beginner

852: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetOrder(), PetscSpacePolynomialSetNumVariables()
853: @*/
854: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
855: {

861:   PetscTryMethod(sp,"PetscSpacePolynomialGetTensor_C",(PetscSpace,PetscBool*),(sp,tensor));
862:   return(0);
863: }

867: static PetscErrorCode PetscSpacePolynomialSetTensor_Polynomial(PetscSpace sp, PetscBool tensor)
868: {
869:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

872:   poly->tensor = tensor;
873:   return(0);
874: }

878: static PetscErrorCode PetscSpacePolynomialGetTensor_Polynomial(PetscSpace sp, PetscBool *tensor)
879: {
880:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

885:   *tensor = poly->tensor;
886:   return(0);
887: }

891: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
892: {

896:   sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
897:   sp->ops->setup          = PetscSpaceSetUp_Polynomial;
898:   sp->ops->view           = PetscSpaceView_Polynomial;
899:   sp->ops->destroy        = PetscSpaceDestroy_Polynomial;
900:   sp->ops->getdimension   = PetscSpaceGetDimension_Polynomial;
901:   sp->ops->evaluate       = PetscSpaceEvaluate_Polynomial;
902:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Polynomial);
903:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialSetTensor_C", PetscSpacePolynomialSetTensor_Polynomial);
904:   return(0);
905: }

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

910:   Level: intermediate

912: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
913: M*/

917: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
918: {
919:   PetscSpace_Poly *poly;
920:   PetscErrorCode   ierr;

924:   PetscNewLog(sp,&poly);
925:   sp->data = poly;

927:   poly->numVariables = 0;
928:   poly->symmetric    = PETSC_FALSE;
929:   poly->tensor       = PETSC_FALSE;
930:   poly->degrees      = NULL;

932:   PetscSpaceInitialize_Polynomial(sp);
933:   return(0);
934: }

938: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
939: {
940:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

944:   poly->symmetric = sym;
945:   return(0);
946: }

950: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
951: {
952:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

957:   *sym = poly->symmetric;
958:   return(0);
959: }

963: PetscErrorCode PetscSpacePolynomialSetNumVariables(PetscSpace sp, PetscInt n)
964: {
965:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

969:   poly->numVariables = n;
970:   return(0);
971: }

975: PetscErrorCode PetscSpacePolynomialGetNumVariables(PetscSpace sp, PetscInt *n)
976: {
977:   PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;

982:   *n = poly->numVariables;
983:   return(0);
984: }

988: PetscErrorCode PetscSpaceSetFromOptions_DG(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
989: {
990:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

994:   PetscOptionsHead(PetscOptionsObject,"PetscSpace DG options");
995:   PetscOptionsInt("-petscspace_dg_num_variables", "The number of different variables, e.g. x and y", "PetscSpaceDGSetNumVariables", dg->numVariables, &dg->numVariables, NULL);
996:   PetscOptionsTail();
997:   return(0);
998: }

1002: PetscErrorCode PetscSpaceDGView_Ascii(PetscSpace sp, PetscViewer viewer)
1003: {
1004:   PetscSpace_DG    *dg = (PetscSpace_DG *) sp->data;
1005:   PetscViewerFormat format;
1006:   PetscErrorCode    ierr;

1009:   PetscViewerGetFormat(viewer, &format);
1010:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1011:     PetscViewerASCIIPrintf(viewer, "DG space in dimension %d:\n", dg->numVariables);
1012:     PetscViewerASCIIPushTab(viewer);
1013:     PetscQuadratureView(dg->quad, viewer);
1014:     PetscViewerASCIIPopTab(viewer);
1015:   } else {
1016:     PetscViewerASCIIPrintf(viewer, "DG space in dimension %d on %d points\n", dg->numVariables, dg->quad->numPoints);
1017:   }
1018:   return(0);
1019: }

1023: PetscErrorCode PetscSpaceView_DG(PetscSpace sp, PetscViewer viewer)
1024: {
1025:   PetscBool      iascii;

1031:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1032:   if (iascii) {PetscSpaceDGView_Ascii(sp, viewer);}
1033:   return(0);
1034: }

1038: PetscErrorCode PetscSpaceSetUp_DG(PetscSpace sp)
1039: {
1040:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

1044:   if (!dg->quad->points && sp->order) {
1045:     PetscDTGaussJacobiQuadrature(dg->numVariables, sp->order, -1.0, 1.0, &dg->quad);
1046:   }
1047:   return(0);
1048: }

1052: PetscErrorCode PetscSpaceDestroy_DG(PetscSpace sp)
1053: {
1054:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

1058:   PetscQuadratureDestroy(&dg->quad);
1059:   return(0);
1060: }

1064: PetscErrorCode PetscSpaceGetDimension_DG(PetscSpace sp, PetscInt *dim)
1065: {
1066:   PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;

1069:   *dim = dg->quad->numPoints;
1070:   return(0);
1071: }

1075: PetscErrorCode PetscSpaceEvaluate_DG(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
1076: {
1077:   PetscSpace_DG *dg  = (PetscSpace_DG *) sp->data;
1078:   PetscInt       dim = dg->numVariables, d, p;

1082:   if (D || H) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Cannot calculate derivatives for a DG space");
1083:   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);
1084:   PetscMemzero(B, npoints*npoints * sizeof(PetscReal));
1085:   for (p = 0; p < npoints; ++p) {
1086:     for (d = 0; d < dim; ++d) {
1087:       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]);
1088:     }
1089:     B[p*npoints+p] = 1.0;
1090:   }
1091:   return(0);
1092: }

1096: PetscErrorCode PetscSpaceInitialize_DG(PetscSpace sp)
1097: {
1099:   sp->ops->setfromoptions = PetscSpaceSetFromOptions_DG;
1100:   sp->ops->setup          = PetscSpaceSetUp_DG;
1101:   sp->ops->view           = PetscSpaceView_DG;
1102:   sp->ops->destroy        = PetscSpaceDestroy_DG;
1103:   sp->ops->getdimension   = PetscSpaceGetDimension_DG;
1104:   sp->ops->evaluate       = PetscSpaceEvaluate_DG;
1105:   return(0);
1106: }

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

1111:   Level: intermediate

1113: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
1114: M*/

1118: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_DG(PetscSpace sp)
1119: {
1120:   PetscSpace_DG *dg;

1125:   PetscNewLog(sp,&dg);
1126:   sp->data = dg;

1128:   dg->numVariables    = 0;
1129:   dg->quad->dim       = 0;
1130:   dg->quad->numPoints = 0;
1131:   dg->quad->points    = NULL;
1132:   dg->quad->weights   = NULL;

1134:   PetscSpaceInitialize_DG(sp);
1135:   return(0);
1136: }


1139: PetscClassId PETSCDUALSPACE_CLASSID = 0;

1141: PetscFunctionList PetscDualSpaceList              = NULL;
1142: PetscBool         PetscDualSpaceRegisterAllCalled = PETSC_FALSE;

1146: /*@C
1147:   PetscDualSpaceRegister - Adds a new PetscDualSpace implementation

1149:   Not Collective

1151:   Input Parameters:
1152: + name        - The name of a new user-defined creation routine
1153: - create_func - The creation routine itself

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

1158:   Sample usage:
1159: .vb
1160:     PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
1161: .ve

1163:   Then, your PetscDualSpace type can be chosen with the procedural interface via
1164: .vb
1165:     PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
1166:     PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
1167: .ve
1168:    or at runtime via the option
1169: .vb
1170:     -petscdualspace_type my_dual_space
1171: .ve

1173:   Level: advanced

1175: .keywords: PetscDualSpace, register
1176: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()

1178: @*/
1179: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
1180: {

1184:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
1185:   return(0);
1186: }

1190: /*@C
1191:   PetscDualSpaceSetType - Builds a particular PetscDualSpace

1193:   Collective on PetscDualSpace

1195:   Input Parameters:
1196: + sp   - The PetscDualSpace object
1197: - name - The kind of space

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

1202:   Level: intermediate

1204: .keywords: PetscDualSpace, set, type
1205: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
1206: @*/
1207: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
1208: {
1209:   PetscErrorCode (*r)(PetscDualSpace);
1210:   PetscBool      match;

1215:   PetscObjectTypeCompare((PetscObject) sp, name, &match);
1216:   if (match) return(0);

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

1222:   if (sp->ops->destroy) {
1223:     (*sp->ops->destroy)(sp);
1224:     sp->ops->destroy = NULL;
1225:   }
1226:   (*r)(sp);
1227:   PetscObjectChangeTypeName((PetscObject) sp, name);
1228:   return(0);
1229: }

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

1236:   Not Collective

1238:   Input Parameter:
1239: . sp  - The PetscDualSpace

1241:   Output Parameter:
1242: . name - The PetscDualSpace type name

1244:   Level: intermediate

1246: .keywords: PetscDualSpace, get, type, name
1247: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
1248: @*/
1249: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
1250: {

1256:   if (!PetscDualSpaceRegisterAllCalled) {
1257:     PetscDualSpaceRegisterAll();
1258:   }
1259:   *name = ((PetscObject) sp)->type_name;
1260:   return(0);
1261: }

1265: /*@
1266:   PetscDualSpaceView - Views a PetscDualSpace

1268:   Collective on PetscDualSpace

1270:   Input Parameter:
1271: + sp - the PetscDualSpace object to view
1272: - v  - the viewer

1274:   Level: developer

1276: .seealso PetscDualSpaceDestroy()
1277: @*/
1278: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
1279: {

1284:   if (!v) {
1285:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
1286:   }
1287:   if (sp->ops->view) {
1288:     (*sp->ops->view)(sp, v);
1289:   }
1290:   return(0);
1291: }

1295: /*@
1296:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

1298:   Collective on PetscDualSpace

1300:   Input Parameter:
1301: . sp - the PetscDualSpace object to set options for

1303:   Options Database:
1304: . -petscspace_order the approximation order of the space

1306:   Level: developer

1308: .seealso PetscDualSpaceView()
1309: @*/
1310: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
1311: {
1312:   const char    *defaultType;
1313:   char           name[256];
1314:   PetscBool      flg;

1319:   if (!((PetscObject) sp)->type_name) {
1320:     defaultType = PETSCDUALSPACELAGRANGE;
1321:   } else {
1322:     defaultType = ((PetscObject) sp)->type_name;
1323:   }
1324:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

1326:   PetscObjectOptionsBegin((PetscObject) sp);
1327:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
1328:   if (flg) {
1329:     PetscDualSpaceSetType(sp, name);
1330:   } else if (!((PetscObject) sp)->type_name) {
1331:     PetscDualSpaceSetType(sp, defaultType);
1332:   }
1333:   PetscOptionsInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
1334:   if (sp->ops->setfromoptions) {
1335:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
1336:   }
1337:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1338:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
1339:   PetscOptionsEnd();
1340:   PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
1341:   return(0);
1342: }

1346: /*@
1347:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

1349:   Collective on PetscDualSpace

1351:   Input Parameter:
1352: . sp - the PetscDualSpace object to setup

1354:   Level: developer

1356: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
1357: @*/
1358: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
1359: {

1364:   if (sp->setupcalled) return(0);
1365:   sp->setupcalled = PETSC_TRUE;
1366:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
1367:   return(0);
1368: }

1372: /*@
1373:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

1375:   Collective on PetscDualSpace

1377:   Input Parameter:
1378: . sp - the PetscDualSpace object to destroy

1380:   Level: developer

1382: .seealso PetscDualSpaceView()
1383: @*/
1384: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
1385: {
1386:   PetscInt       dim, f;

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

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

1396:   PetscDualSpaceGetDimension(*sp, &dim);
1397:   for (f = 0; f < dim; ++f) {
1398:     PetscQuadratureDestroy(&(*sp)->functional[f]);
1399:   }
1400:   PetscFree((*sp)->functional);
1401:   DMDestroy(&(*sp)->dm);

1403:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
1404:   PetscHeaderDestroy(sp);
1405:   return(0);
1406: }

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

1413:   Collective on MPI_Comm

1415:   Input Parameter:
1416: . comm - The communicator for the PetscDualSpace object

1418:   Output Parameter:
1419: . sp - The PetscDualSpace object

1421:   Level: beginner

1423: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
1424: @*/
1425: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
1426: {
1427:   PetscDualSpace s;

1432:   PetscCitationsRegister(FECitation,&FEcite);
1433:   *sp  = NULL;
1434:   PetscFEInitializePackage();

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

1438:   s->order = 0;
1439:   s->setupcalled = PETSC_FALSE;

1441:   *sp = s;
1442:   return(0);
1443: }

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

1450:   Collective on PetscDualSpace

1452:   Input Parameter:
1453: . sp - The original PetscDualSpace

1455:   Output Parameter:
1456: . spNew - The duplicate PetscDualSpace

1458:   Level: beginner

1460: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
1461: @*/
1462: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
1463: {

1469:   (*sp->ops->duplicate)(sp, spNew);
1470:   return(0);
1471: }

1475: /*@
1476:   PetscDualSpaceGetDM - Get the DM representing the reference cell

1478:   Not collective

1480:   Input Parameter:
1481: . sp - The PetscDualSpace

1483:   Output Parameter:
1484: . dm - The reference cell

1486:   Level: intermediate

1488: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
1489: @*/
1490: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
1491: {
1495:   *dm = sp->dm;
1496:   return(0);
1497: }

1501: /*@
1502:   PetscDualSpaceSetDM - Get the DM representing the reference cell

1504:   Not collective

1506:   Input Parameters:
1507: + sp - The PetscDualSpace
1508: - dm - The reference cell

1510:   Level: intermediate

1512: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
1513: @*/
1514: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
1515: {

1521:   DMDestroy(&sp->dm);
1522:   PetscObjectReference((PetscObject) dm);
1523:   sp->dm = dm;
1524:   return(0);
1525: }

1529: /*@
1530:   PetscDualSpaceGetOrder - Get the order of the dual space

1532:   Not collective

1534:   Input Parameter:
1535: . sp - The PetscDualSpace

1537:   Output Parameter:
1538: . order - The order

1540:   Level: intermediate

1542: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
1543: @*/
1544: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
1545: {
1549:   *order = sp->order;
1550:   return(0);
1551: }

1555: /*@
1556:   PetscDualSpaceSetOrder - Set the order of the dual space

1558:   Not collective

1560:   Input Parameters:
1561: + sp - The PetscDualSpace
1562: - order - The order

1564:   Level: intermediate

1566: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
1567: @*/
1568: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
1569: {
1572:   sp->order = order;
1573:   return(0);
1574: }

1578: /*@
1579:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

1581:   Not collective

1583:   Input Parameter:
1584: . sp - The PetscDualSpace

1586:   Output Parameter:
1587: . tensor - Whether the dual space has tensor layout (vs. simplicial)

1589:   Level: intermediate

1591: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
1592: @*/
1593: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
1594: {

1600:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
1601:   return(0);
1602: }

1606: /*@
1607:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

1609:   Not collective

1611:   Input Parameters:
1612: + sp - The PetscDualSpace
1613: - tensor - Whether the dual space has tensor layout (vs. simplicial)

1615:   Level: intermediate

1617: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
1618: @*/
1619: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
1620: {

1625:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
1626:   return(0);
1627: }

1631: /*@
1632:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

1634:   Not collective

1636:   Input Parameters:
1637: + sp - The PetscDualSpace
1638: - i  - The basis number

1640:   Output Parameter:
1641: . functional - The basis functional

1643:   Level: intermediate

1645: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
1646: @*/
1647: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
1648: {
1649:   PetscInt       dim;

1655:   PetscDualSpaceGetDimension(sp, &dim);
1656:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
1657:   *functional = sp->functional[i];
1658:   return(0);
1659: }

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

1666:   Not collective

1668:   Input Parameter:
1669: . sp - The PetscDualSpace

1671:   Output Parameter:
1672: . dim - The dimension

1674:   Level: intermediate

1676: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1677: @*/
1678: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
1679: {

1685:   *dim = 0;
1686:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
1687:   return(0);
1688: }

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

1695:   Not collective

1697:   Input Parameter:
1698: . sp - The PetscDualSpace

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

1703:   Level: intermediate

1705: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1706: @*/
1707: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
1708: {

1714:   (*sp->ops->getnumdof)(sp, numDof);
1715:   if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
1716:   return(0);
1717: }

1721: /*@
1722:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

1724:   Collective on PetscDualSpace

1726:   Input Parameters:
1727: + sp      - The PetscDualSpace
1728: . dim     - The spatial dimension
1729: - simplex - Flag for simplex, otherwise use a tensor-product cell

1731:   Output Parameter:
1732: . refdm - The reference cell

1734:   Level: advanced

1736: .keywords: PetscDualSpace, reference cell
1737: .seealso: PetscDualSpaceCreate(), DMPLEX
1738: @*/
1739: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1740: {

1744:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
1745:   return(0);
1746: }

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

1753:   Input Parameters:
1754: + sp      - The PetscDualSpace object
1755: . f       - The basis functional index
1756: . time    - The time
1757: . geom    - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1758: . numComp - The number of components for the function
1759: . func    - The input function
1760: - ctx     - A context for the function

1762:   Output Parameter:
1763: . value   - numComp output values

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

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

1770:   Level: developer

1772: .seealso: PetscDualSpaceCreate()
1773: @*/
1774: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFECellGeom *geom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1775: {
1776:   DM               dm;
1777:   PetscQuadrature  quad;
1778:   PetscReal        x[3];
1779:   PetscScalar     *val;
1780:   PetscInt         dim, q, c;
1781:   PetscErrorCode   ierr;

1786:   dim  = geom->dim;
1787:   PetscDualSpaceGetDM(sp, &dm);
1788:   PetscDualSpaceGetFunctional(sp, f, &quad);
1789:   DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1790:   for (c = 0; c < numComp; ++c) value[c] = 0.0;
1791:   for (q = 0; q < quad->numPoints; ++q) {
1792:     CoordinatesRefToReal(geom->dimEmbed, dim, geom->v0, geom->J, &quad->points[q*dim], x);
1793:     (*func)(geom->dimEmbed, time, x, numComp, val, ctx);
1794:     for (c = 0; c < numComp; ++c) {
1795:       value[c] += val[c]*quad->weights[q];
1796:     }
1797:   }
1798:   DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1799:   return(0);
1800: }

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

1807:   Input Parameters:
1808: + sp      - The PetscDualSpace object
1809: . f       - The basis functional index
1810: . time    - The time
1811: . geom    - A context with geometric information for this cell, we currently just use the centroid
1812: . numComp - The number of components for the function
1813: . func    - The input function
1814: - ctx     - A context for the function

1816:   Output Parameter:
1817: . value   - numComp output values

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

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

1824:   Level: developer

1826: .seealso: PetscDualSpaceCreate()
1827: @*/
1828: PetscErrorCode PetscDualSpaceApplyFVM(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFVCellGeom *geom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1829: {
1830:   DM               dm;
1831:   PetscQuadrature  quad;
1832:   PetscScalar     *val;
1833:   PetscInt         dimEmbed, q, c;
1834:   PetscErrorCode   ierr;

1839:   PetscDualSpaceGetDM(sp, &dm);
1840:   DMGetCoordinateDim(dm, &dimEmbed);
1841:   PetscDualSpaceGetFunctional(sp, f, &quad);
1842:   DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1843:   for (c = 0; c < numComp; ++c) value[c] = 0.0;
1844:   for (q = 0; q < quad->numPoints; ++q) {
1845:     (*func)(dimEmbed, time, geom->centroid, numComp, val, ctx);
1846:     for (c = 0; c < numComp; ++c) {
1847:       value[c] += val[c]*quad->weights[q];
1848:     }
1849:   }
1850:   DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1851:   return(0);
1852: }

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

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

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

1865:   Not collective

1867:   Input Parameters:
1868: + sp - the PetscDualSpace object
1869: - height - the height of the mesh point for which the subspace is desired

1871:   Output Parameters:
1872:   bdsp - the subspace: must be destroyed by the user

1874:   Level: advanced

1876: .seealso: PetscDualSpace
1877: @*/
1878: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
1879: {

1885:   *bdsp = NULL;
1886:   if (sp->ops->getheightsubspace) {
1887:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1888:   }
1889:   return(0);
1890: }

1894: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
1895: {
1896:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

1899:   *tensor = lag->tensorSpace;
1900:   return(0);
1901: }

1905: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
1906: {
1907:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

1910:   lag->tensorSpace = tensor;
1911:   return(0);
1912: }

1914: #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c)

1916: #define CartIndex(perEdge,a,b) (perEdge*(a)+b)

1920: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1921: {

1923:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1924:   PetscInt           dim, order, p;
1925:   PetscErrorCode     ierr;

1928:   PetscDualSpaceGetOrder(sp,&order);
1929:   DMGetDimension(sp->dm,&dim);
1930:   if (!dim || !lag->continuous || order < 3) return(0);
1931:   if (dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Lagrange symmetries not implemented for dim = %D > 3",dim);
1932:   if (!lag->symmetries) { /* store symmetries */
1933:     PetscDualSpace hsp;
1934:     DM             K;
1935:     PetscInt       numPoints = 1, d;
1936:     PetscInt       numFaces;
1937:     PetscInt       ***symmetries;
1938:     const PetscInt ***hsymmetries;

1940:     if (lag->simplexCell) {
1941:       numFaces = 1 + dim;
1942:       for (d = 0; d < dim; d++) numPoints = numPoints * 2 + 1;
1943:     }
1944:     else {
1945:       numPoints = PetscPowInt(3,dim);
1946:       numFaces  = 2 * dim;
1947:     }
1948:     PetscCalloc1(numPoints,&symmetries);
1949:     if (0 < dim && dim < 3) { /* compute self symmetries */
1950:       PetscInt **cellSymmetries;

1952:       lag->numSelfSym = 2 * numFaces;
1953:       lag->selfSymOff = numFaces;
1954:       PetscCalloc1(2*numFaces,&cellSymmetries);
1955:       /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
1956:       symmetries[0] = &cellSymmetries[numFaces];
1957:       if (dim == 1) {
1958:         PetscInt dofPerEdge = order - 1;

1960:         if (dofPerEdge > 1) {
1961:           PetscInt i, *reverse;

1963:           PetscMalloc1(dofPerEdge,&reverse);
1964:           for (i = 0; i < dofPerEdge; i++) reverse[i] = (dofPerEdge - 1 - i);
1965:           symmetries[0][-2] = reverse;

1967:           /* yes, this is redundant, but it makes it easier to cleanup if I don't have to worry about what not to free */
1968:           PetscMalloc1(dofPerEdge,&reverse);
1969:           for (i = 0; i < dofPerEdge; i++) reverse[i] = (dofPerEdge - 1 - i);
1970:           symmetries[0][1] = reverse;
1971:         }
1972:       } else {
1973:         PetscInt dofPerEdge = lag->simplexCell ? (order - 2) : (order - 1), s;

1975:         if (dofPerEdge > 1) {
1976:           for (s = -numFaces; s < numFaces; s++) {
1977:             PetscInt *sym, i, j, k, l;

1979:             if (!s) continue;
1980:             if (lag->simplexCell) {
1981:               PetscMalloc1((dofPerEdge * (dofPerEdge + 1))/2,&sym);
1982:               for (j = 0, l = 0; j < dofPerEdge; j++) {
1983:                 for (k = 0; k < dofPerEdge - j; k++, l++) {
1984:                   i = dofPerEdge - 1 - j - k;
1985:                   switch (s) {
1986:                   case -3:
1987:                     sym[l] = BaryIndex(dofPerEdge,i,k,j);
1988:                     break;
1989:                   case -2:
1990:                     sym[l] = BaryIndex(dofPerEdge,j,i,k);
1991:                     break;
1992:                   case -1:
1993:                     sym[l] = BaryIndex(dofPerEdge,k,j,i);
1994:                     break;
1995:                   case 1:
1996:                     sym[l] = BaryIndex(dofPerEdge,k,i,j);
1997:                     break;
1998:                   case 2:
1999:                     sym[l] = BaryIndex(dofPerEdge,j,k,i);
2000:                     break;
2001:                   }
2002:                 }
2003:               }
2004:             } else {
2005:               PetscMalloc1(dofPerEdge * dofPerEdge,&sym);
2006:               for (j = 0, l = 0; j < dofPerEdge; j++) {
2007:                 for (k = 0; k < dofPerEdge; k++, l++) {
2008:                   switch (s) {
2009:                   case -4:
2010:                     sym[l] = CartIndex(dofPerEdge,k,j);
2011:                     break;
2012:                   case -3:
2013:                     sym[l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),k);
2014:                     break;
2015:                   case -2:
2016:                     sym[l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),(dofPerEdge - 1 - j));
2017:                     break;
2018:                   case -1:
2019:                     sym[l] = CartIndex(dofPerEdge,j,(dofPerEdge - 1 - k));
2020:                     break;
2021:                   case 1:
2022:                     sym[l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),j);
2023:                     break;
2024:                   case 2:
2025:                     sym[l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),(dofPerEdge - 1 - k));
2026:                     break;
2027:                   case 3:
2028:                     sym[l] = CartIndex(dofPerEdge,k,(dofPerEdge - 1 - j));
2029:                     break;
2030:                   }
2031:                 }
2032:               }
2033:             }
2034:             symmetries[0][s] = sym;
2035:           }
2036:         }
2037:       }
2038:     }
2039:     PetscDualSpaceGetHeightSubspace(sp,1,&hsp);
2040:     PetscDualSpaceGetSymmetries(hsp,&hsymmetries,NULL);
2041:     if (hsymmetries) {
2042:       PetscBool      *seen;
2043:       const PetscInt *cone;
2044:       PetscInt       KclosureSize, *Kclosure = NULL;

2046:       PetscDualSpaceGetDM(sp,&K);
2047:       PetscCalloc1(numPoints,&seen);
2048:       DMPlexGetCone(K,0,&cone);
2049:       DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
2050:       for (p = 0; p < numFaces; p++) {
2051:         PetscInt closureSize, *closure = NULL, q;

2053:         DMPlexGetTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);
2054:         for (q = 0; q < closureSize; q++) {
2055:           PetscInt point = closure[2*q], r;

2057:           if(!seen[point]) {
2058:             for (r = 0; r < KclosureSize; r++) {
2059:               if (Kclosure[2 * r] == point) break;
2060:             }
2061:             seen[point] = PETSC_TRUE;
2062:             symmetries[r] = (PetscInt **) hsymmetries[q];
2063:           }
2064:         }
2065:         DMPlexRestoreTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);
2066:       }
2067:       DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
2068:       PetscFree(seen);
2069:     }
2070:     lag->symmetries = symmetries;
2071:   }
2072:   if (perms) *perms = (const PetscInt ***) lag->symmetries;
2073:   return(0);
2074: }

2078: PETSC_EXTERN PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
2079: {

2084:   if (perms) {
2086:     *perms = NULL;
2087:   }
2088:   if (flips) {
2090:     *flips = NULL;
2091:   }
2092:   if (sp->ops->getsymmetries) {
2093:     (sp->ops->getsymmetries)(sp,perms,flips);
2094:   }
2095:   return(0);
2096: }

2100: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
2101: {
2102:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2103:   PetscReal           D   = 1.0;
2104:   PetscInt            n, i;
2105:   PetscErrorCode      ierr;

2108:   *dim = -1;                    /* Ensure that the compiler knows *dim is set. */
2109:   DMGetDimension(sp->dm, &n);
2110:   if (!lag->tensorSpace) {
2111:     for (i = 1; i <= n; ++i) {
2112:       D *= ((PetscReal) (order+i))/i;
2113:     }
2114:     *dim = (PetscInt) (D + 0.5);
2115:   } else {
2116:     *dim = 1;
2117:     for (i = 0; i < n; ++i) *dim *= (order+1);
2118:   }
2119:   return(0);
2120: }

2124: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
2125: {
2126:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2127:   PetscBool          continuous, tensor;
2128:   PetscInt           order;
2129:   PetscErrorCode     ierr;

2134:   PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
2135:   PetscDualSpaceGetOrder(sp,&order);
2136:   if (height == 0) {
2137:     PetscObjectReference((PetscObject)sp);
2138:     *bdsp = sp;
2139:   }
2140:   else if (continuous == PETSC_FALSE || !order) {
2141:     *bdsp = NULL;
2142:   }
2143:   else {
2144:     DM dm, K;
2145:     PetscInt dim;

2147:     PetscDualSpaceGetDM(sp,&dm);
2148:     DMGetDimension(dm,&dim);
2149:     if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim);}
2150:     PetscDualSpaceDuplicate(sp,bdsp);
2151:     PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplexCell, &K);
2152:     PetscDualSpaceSetDM(*bdsp, K);
2153:     DMDestroy(&K);
2154:     PetscDualSpaceLagrangeGetTensor(sp,&tensor);
2155:     PetscDualSpaceLagrangeSetTensor(*bdsp,tensor);
2156:     PetscDualSpaceSetUp(*bdsp);
2157:   }
2158:   return(0);
2159: }

2163: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
2164: {
2165:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2166:   DM                  dm    = sp->dm;
2167:   PetscInt            order = sp->order;
2168:   PetscBool           continuous;
2169:   PetscSection        csection;
2170:   Vec                 coordinates;
2171:   PetscReal          *qpoints, *qweights;
2172:   PetscInt            depth, dim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0;
2173:   PetscBool           simplex, tensorSpace;
2174:   PetscErrorCode      ierr;

2177:   /* Classify element type */
2178:   if (!order) lag->continuous = PETSC_FALSE;
2179:   continuous = lag->continuous;
2180:   DMGetDimension(dm, &dim);
2181:   DMPlexGetDepth(dm, &depth);
2182:   DMPlexGetChart(dm, &pStart, &pEnd);
2183:   PetscCalloc1(dim+1, &lag->numDof);
2184:   PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
2185:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
2186:   DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);
2187:   DMGetCoordinateSection(dm, &csection);
2188:   DMGetCoordinatesLocal(dm, &coordinates);
2189:   if (depth == 1) {
2190:     if      (coneSize == dim+1)    simplex = PETSC_TRUE;
2191:     else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
2192:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
2193:   }
2194:   else if (depth == dim) {
2195:     if      (coneSize == dim+1)   simplex = PETSC_TRUE;
2196:     else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
2197:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
2198:   }
2199:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
2200:   lag->simplexCell = simplex;
2201:   if (dim > 1 && continuous && lag->simplexCell == lag->tensorSpace) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Mismatching simplex/tensor cells and spaces only allowed for discontinuous elements");
2202:   tensorSpace    = lag->tensorSpace;
2203:   lag->height    = 0;
2204:   lag->subspaces = NULL;
2205:   if (continuous && sp->order > 0 && dim > 0) {
2206:     PetscInt i;

2208:     lag->height = dim;
2209:     PetscMalloc1(dim,&lag->subspaces);
2210:     PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);
2211:     PetscDualSpaceSetUp(lag->subspaces[0]);
2212:     for (i = 1; i < dim; i++) {
2213:       PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);
2214:       PetscObjectReference((PetscObject)(lag->subspaces[i]));
2215:     }
2216:   }
2217:   PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
2218:   pdimMax *= (pStratEnd[depth] - pStratStart[depth]);
2219:   PetscMalloc1(pdimMax, &sp->functional);
2220:   if (!dim) {
2221:     PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
2222:     PetscMalloc1(1, &qweights);
2223:     PetscQuadratureSetOrder(sp->functional[f], 0);
2224:     PetscQuadratureSetData(sp->functional[f], 0, 1, NULL, qweights);
2225:     qweights[0] = 1.0;
2226:     ++f;
2227:     lag->numDof[0] = 1;
2228:   } else {
2229:     PetscInt     *tup;
2230:     PetscReal    *v0, *hv0, *J, *invJ, detJ, hdetJ;
2231:     PetscSection section;

2233:     PetscSectionCreate(PETSC_COMM_SELF,&section);
2234:     PetscSectionSetChart(section,pStart,pEnd);
2235:     PetscCalloc5(dim+1,&tup,dim,&v0,dim,&hv0,dim*dim,&J,dim*dim,&invJ);
2236:     for (p = pStart; p < pEnd; p++) {
2237:       PetscInt       pointDim, d, nFunc = 0;
2238:       PetscDualSpace hsp;

2240:       DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
2241:       for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;}
2242:       pointDim = (depth == 1 && d == 1) ? dim : d;
2243:       hsp = ((pointDim < dim) && lag->subspaces) ? lag->subspaces[dim - pointDim - 1] : NULL;
2244:       if (hsp) {
2245:         PetscDualSpace_Lag *hlag = (PetscDualSpace_Lag *) hsp->data;
2246:         DM                 hdm;

2248:         PetscDualSpaceGetDM(hsp,&hdm);
2249:         DMPlexComputeCellGeometryFEM(hdm, 0, NULL, hv0, NULL, NULL, &hdetJ);
2250:         lag->numDof[pointDim] = nFunc = hlag->numDof[pointDim];
2251:       }
2252:       if (pointDim == dim) {
2253:         /* Cells, create for self */
2254:         PetscInt     orderEff = continuous ? (!tensorSpace ? order-1-dim : order-2) : order;
2255:         PetscReal    denom    = continuous ? order : (!tensorSpace ? order+1+dim : order+2);
2256:         PetscReal    numer    = (!simplex || !tensorSpace) ? 2. : (2./dim);
2257:         PetscReal    dx = numer/denom;
2258:         PetscInt     cdim, d, d2;

2260:         if (orderEff < 0) continue;
2261:         PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);

2263:         PetscMemzero(tup,(dim+1)*sizeof(PetscInt));
2264:         if (!tensorSpace) {
2265:           while (!tup[dim]) {
2266:             PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
2267:             PetscMalloc1(dim, &qpoints);
2268:             PetscMalloc1(1,   &qweights);
2269:             PetscQuadratureSetOrder(sp->functional[f], 0);
2270:             PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
2271:             for (d = 0; d < dim; ++d) {
2272:               qpoints[d] = v0[d];
2273:               for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
2274:             }
2275:             qweights[0] = 1.0;
2276:             ++f;
2277:             LatticePointLexicographic_Internal(dim, orderEff, tup);
2278:           }
2279:         } else {
2280:           while (!tup[dim]) {
2281:             PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
2282:             PetscMalloc1(dim, &qpoints);
2283:             PetscMalloc1(1,   &qweights);
2284:             PetscQuadratureSetOrder(sp->functional[f], 0);
2285:             PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
2286:             for (d = 0; d < dim; ++d) {
2287:               qpoints[d] = v0[d];
2288:               for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
2289:             }
2290:             qweights[0] = 1.0;
2291:             ++f;
2292:             TensorPointLexicographic_Internal(dim, orderEff, tup);
2293:           }
2294:         }
2295:         lag->numDof[dim] = cdim;
2296:       } else { /* transform functionals from subspaces */
2297:         PetscInt q;

2299:         for (q = 0; q < nFunc; q++, f++) {
2300:           PetscQuadrature fn;
2301:           PetscInt        fdim, nPoints, i;
2302:           const PetscReal *points;
2303:           const PetscReal *weights;
2304:           PetscReal       *qpoints;
2305:           PetscReal       *qweights;

2307:           PetscDualSpaceGetFunctional(hsp, q, &fn);
2308:           PetscQuadratureGetData(fn,&fdim,&nPoints,&points,&weights);
2309:           if (fdim != pointDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected height dual space dim %D, got %D",pointDim,fdim);
2310:           PetscMalloc1(nPoints * dim, &qpoints);
2311:           PetscMalloc1(nPoints, &qweights);
2312:           for (i = 0; i < nPoints; i++) {
2313:             PetscInt  j, k;
2314:             PetscReal *qp = &qpoints[i * dim];

2316:             qweights[i] = weights[i];
2317:             for (j = 0; j < dim; j++) qp[j] = v0[j];
2318:             for (j = 0; j < dim; j++) {
2319:               for (k = 0; k < pointDim; k++) qp[j] += J[dim * j + k] * (points[pointDim * i + k] - hv0[k]);
2320:             }
2321:           }
2322:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
2323:           PetscQuadratureSetOrder(sp->functional[f],0);
2324:           PetscQuadratureSetData(sp->functional[f],dim,nPoints,qpoints,qweights);
2325:         }
2326:       }
2327:       PetscSectionSetDof(section,p,lag->numDof[pointDim]);
2328:     }
2329:     PetscFree5(tup,v0,hv0,J,invJ);
2330:     PetscSectionSetUp(section);
2331:     { /* reorder to closure order */
2332:       PetscInt *key, count;
2333:       PetscQuadrature *reorder = NULL;

2335:       PetscCalloc1(f,&key);
2336:       PetscMalloc1(f,&reorder);

2338:       for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) {
2339:         PetscInt *closure = NULL, closureSize, c;

2341:         DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
2342:         for (c = 0; c < closureSize; c++) {
2343:           PetscInt point = closure[2 * c], dof, off, i;

2345:           PetscSectionGetDof(section,point,&dof);
2346:           PetscSectionGetOffset(section,point,&off);
2347:           for (i = 0; i < dof; i++) {
2348:             PetscInt fi = i + off;
2349:             if (!key[fi]) {
2350:               key[fi] = 1;
2351:               reorder[count++] = sp->functional[fi];
2352:             }
2353:           }
2354:         }
2355:         DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
2356:       }
2357:       PetscFree(sp->functional);
2358:       sp->functional = reorder;
2359:       PetscFree(key);
2360:     }
2361:     PetscSectionDestroy(&section);
2362:   }
2363:   if (pStratEnd[depth] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d not equal to dimension %d", f, pdimMax);
2364:   PetscFree2(pStratStart,pStratEnd);
2365:   if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
2366:   return(0);
2367: }

2371: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
2372: {
2373:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2374:   PetscInt            i;
2375:   PetscErrorCode      ierr;

2378:   if (lag->symmetries) {
2379:     PetscInt **selfSyms = lag->symmetries[0];

2381:     if (selfSyms) {
2382:       PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];

2384:       for (i = 0; i < lag->numSelfSym; i++) {
2385:         PetscFree(allocated[i]);
2386:       }
2387:       PetscFree(allocated);
2388:     }
2389:     PetscFree(lag->symmetries);
2390:   }
2391:   for (i = 0; i < lag->height; i++) {
2392:     PetscDualSpaceDestroy(&lag->subspaces[i]);
2393:   }
2394:   PetscFree(lag->subspaces);
2395:   PetscFree(lag->numDof);
2396:   PetscFree(lag);
2397:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
2398:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
2399:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
2400:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
2401:   return(0);
2402: }

2406: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
2407: {
2408:   PetscInt       order;
2409:   PetscBool      cont, tensor;

2413:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2414:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
2415:   PetscDualSpaceGetOrder(sp, &order);
2416:   PetscDualSpaceSetOrder(*spNew, order);
2417:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
2418:   PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
2419:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
2420:   PetscDualSpaceLagrangeSetTensor(*spNew, tensor);
2421:   return(0);
2422: }

2426: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
2427: {
2428:   PetscBool      continuous, tensor, flg;

2432:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
2433:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
2434:   PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
2435:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
2436:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
2437:   PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetContinuity", tensor, &tensor, &flg);
2438:   if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
2439:   PetscOptionsTail();
2440:   return(0);
2441: }

2445: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
2446: {
2447:   DM              K;
2448:   const PetscInt *numDof;
2449:   PetscInt        spatialDim, Nc, size = 0, d;
2450:   PetscErrorCode  ierr;

2453:   PetscDualSpaceGetDM(sp, &K);
2454:   PetscDualSpaceGetNumDof(sp, &numDof);
2455:   DMGetDimension(K, &spatialDim);
2456:   DMPlexGetHeightStratum(K, 0, NULL, &Nc);
2457:   if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
2458:   for (d = 0; d <= spatialDim; ++d) {
2459:     PetscInt pStart, pEnd;

2461:     DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
2462:     size += (pEnd-pStart)*numDof[d];
2463:   }
2464:   *dim = size;
2465:   return(0);
2466: }

2470: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
2471: {
2472:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2475:   *numDof = lag->numDof;
2476:   return(0);
2477: }

2481: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2482: {
2483:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2488:   *continuous = lag->continuous;
2489:   return(0);
2490: }

2494: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2495: {
2496:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2500:   lag->continuous = continuous;
2501:   return(0);
2502: }

2506: /*@
2507:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2509:   Not Collective

2511:   Input Parameter:
2512: . sp         - the PetscDualSpace

2514:   Output Parameter:
2515: . continuous - flag for element continuity

2517:   Level: intermediate

2519: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2520: .seealso: PetscDualSpaceLagrangeSetContinuity()
2521: @*/
2522: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2523: {

2529:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2530:   return(0);
2531: }

2535: /*@
2536:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2538:   Logically Collective on PetscDualSpace

2540:   Input Parameters:
2541: + sp         - the PetscDualSpace
2542: - continuous - flag for element continuity

2544:   Options Database:
2545: . -petscdualspace_lagrange_continuity <bool>

2547:   Level: intermediate

2549: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2550: .seealso: PetscDualSpaceLagrangeGetContinuity()
2551: @*/
2552: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2553: {

2559:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2560:   return(0);
2561: }

2565: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
2566: {
2567:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2568:   PetscErrorCode     ierr;

2573:   if (height == 0) {
2574:     *bdsp = sp;
2575:   }
2576:   else {
2577:     DM dm;
2578:     PetscInt dim;

2580:     PetscDualSpaceGetDM(sp,&dm);
2581:     DMGetDimension(dm,&dim);
2582:     if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim);}
2583:     if (height <= lag->height) {
2584:       *bdsp = lag->subspaces[height-1];
2585:     }
2586:     else {
2587:       *bdsp = NULL;
2588:     }
2589:   }
2590:   return(0);
2591: }

2595: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2596: {
2598:   sp->ops->setfromoptions    = PetscDualSpaceSetFromOptions_Lagrange;
2599:   sp->ops->setup             = PetscDualSpaceSetUp_Lagrange;
2600:   sp->ops->view              = NULL;
2601:   sp->ops->destroy           = PetscDualSpaceDestroy_Lagrange;
2602:   sp->ops->duplicate         = PetscDualSpaceDuplicate_Lagrange;
2603:   sp->ops->getdimension      = PetscDualSpaceGetDimension_Lagrange;
2604:   sp->ops->getnumdof         = PetscDualSpaceGetNumDof_Lagrange;
2605:   sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
2606:   sp->ops->getsymmetries     = PetscDualSpaceGetSymmetries_Lagrange;
2607:   return(0);
2608: }

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

2613:   Level: intermediate

2615: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2616: M*/

2620: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2621: {
2622:   PetscDualSpace_Lag *lag;
2623:   PetscErrorCode      ierr;

2627:   PetscNewLog(sp,&lag);
2628:   sp->data = lag;

2630:   lag->numDof      = NULL;
2631:   lag->simplexCell = PETSC_TRUE;
2632:   lag->tensorSpace = PETSC_FALSE;
2633:   lag->continuous  = PETSC_TRUE;

2635:   PetscDualSpaceInitialize_Lagrange(sp);
2636:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2637:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2638:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
2639:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
2640:   return(0);
2641: }

2645: PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
2646: {
2647:   PetscDualSpace_Simple *s  = (PetscDualSpace_Simple *) sp->data;
2648:   DM                     dm = sp->dm;
2649:   PetscInt               dim;
2650:   PetscErrorCode         ierr;

2653:   DMGetDimension(dm, &dim);
2654:   PetscCalloc1(dim+1, &s->numDof);
2655:   return(0);
2656: }

2660: PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
2661: {
2662:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2663:   PetscErrorCode         ierr;

2666:   PetscFree(s->numDof);
2667:   PetscFree(s);
2668:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
2669:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
2670:   return(0);
2671: }

2675: PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace *spNew)
2676: {
2677:   PetscInt       dim, d;

2681:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2682:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACESIMPLE);
2683:   PetscDualSpaceGetDimension(sp, &dim);
2684:   PetscDualSpaceSimpleSetDimension(*spNew, dim);
2685:   for (d = 0; d < dim; ++d) {
2686:     PetscQuadrature q;

2688:     PetscDualSpaceGetFunctional(sp, d, &q);
2689:     PetscDualSpaceSimpleSetFunctional(*spNew, d, q);
2690:   }
2691:   return(0);
2692: }

2696: PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
2697: {
2699:   return(0);
2700: }

2704: PetscErrorCode PetscDualSpaceGetDimension_Simple(PetscDualSpace sp, PetscInt *dim)
2705: {
2706:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

2709:   *dim = s->dim;
2710:   return(0);
2711: }

2715: PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
2716: {
2717:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2718:   DM                     dm;
2719:   PetscInt               spatialDim, f;
2720:   PetscErrorCode         ierr;

2723:   for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
2724:   PetscFree(sp->functional);
2725:   s->dim = dim;
2726:   PetscCalloc1(s->dim, &sp->functional);
2727:   PetscFree(s->numDof);
2728:   PetscDualSpaceGetDM(sp, &dm);
2729:   DMGetCoordinateDim(dm, &spatialDim);
2730:   PetscCalloc1(spatialDim+1, &s->numDof);
2731:   s->numDof[spatialDim] = dim;
2732:   return(0);
2733: }

2737: PetscErrorCode PetscDualSpaceGetNumDof_Simple(PetscDualSpace sp, const PetscInt **numDof)
2738: {
2739:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

2742:   *numDof = s->numDof;
2743:   return(0);
2744: }

2748: PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
2749: {
2750:   PetscDualSpace_Simple *s   = (PetscDualSpace_Simple *) sp->data;
2751:   PetscReal              vol = 0.0;
2752:   PetscReal             *weights;
2753:   PetscInt               Nq, p;
2754:   PetscErrorCode         ierr;

2757:   if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim);
2758:   PetscQuadratureDuplicate(q, &sp->functional[f]);
2759:   /* Reweight so that it has unit volume */
2760:   PetscQuadratureGetData(sp->functional[f], NULL, &Nq, NULL, (const PetscReal **) &weights);
2761:   for (p = 0; p < Nq; ++p) vol += weights[p];
2762:   for (p = 0; p < Nq; ++p) weights[p] /= vol;
2763:   return(0);
2764: }

2768: /*@
2769:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

2771:   Logically Collective on PetscDualSpace

2773:   Input Parameters:
2774: + sp  - the PetscDualSpace
2775: - dim - the basis dimension

2777:   Level: intermediate

2779: .keywords: PetscDualSpace, dimension
2780: .seealso: PetscDualSpaceSimpleSetFunctional()
2781: @*/
2782: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
2783: {

2789:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
2790:   return(0);
2791: }

2795: /*@
2796:   PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space

2798:   Not Collective

2800:   Input Parameters:
2801: + sp  - the PetscDualSpace
2802: . f - the basis index
2803: - q - the basis functional

2805:   Level: intermediate

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

2809: .keywords: PetscDualSpace, functional
2810: .seealso: PetscDualSpaceSimpleSetDimension()
2811: @*/
2812: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
2813: {

2818:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
2819:   return(0);
2820: }

2824: PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
2825: {
2827:   sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple;
2828:   sp->ops->setup          = PetscDualSpaceSetUp_Simple;
2829:   sp->ops->view           = NULL;
2830:   sp->ops->destroy        = PetscDualSpaceDestroy_Simple;
2831:   sp->ops->duplicate      = PetscDualSpaceDuplicate_Simple;
2832:   sp->ops->getdimension   = PetscDualSpaceGetDimension_Simple;
2833:   sp->ops->getnumdof      = PetscDualSpaceGetNumDof_Simple;
2834:   return(0);
2835: }

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

2840:   Level: intermediate

2842: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2843: M*/

2847: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
2848: {
2849:   PetscDualSpace_Simple *s;
2850:   PetscErrorCode         ierr;

2854:   PetscNewLog(sp,&s);
2855:   sp->data = s;

2857:   s->dim    = 0;
2858:   s->numDof = NULL;

2860:   PetscDualSpaceInitialize_Simple(sp);
2861:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
2862:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
2863:   return(0);
2864: }


2867: PetscClassId PETSCFE_CLASSID = 0;

2869: PetscFunctionList PetscFEList              = NULL;
2870: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

2874: /*@C
2875:   PetscFERegister - Adds a new PetscFE implementation

2877:   Not Collective

2879:   Input Parameters:
2880: + name        - The name of a new user-defined creation routine
2881: - create_func - The creation routine itself

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

2886:   Sample usage:
2887: .vb
2888:     PetscFERegister("my_fe", MyPetscFECreate);
2889: .ve

2891:   Then, your PetscFE type can be chosen with the procedural interface via
2892: .vb
2893:     PetscFECreate(MPI_Comm, PetscFE *);
2894:     PetscFESetType(PetscFE, "my_fe");
2895: .ve
2896:    or at runtime via the option
2897: .vb
2898:     -petscfe_type my_fe
2899: .ve

2901:   Level: advanced

2903: .keywords: PetscFE, register
2904: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()

2906: @*/
2907: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2908: {

2912:   PetscFunctionListAdd(&PetscFEList, sname, function);
2913:   return(0);
2914: }

2918: /*@C
2919:   PetscFESetType - Builds a particular PetscFE

2921:   Collective on PetscFE

2923:   Input Parameters:
2924: + fem  - The PetscFE object
2925: - name - The kind of FEM space

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

2930:   Level: intermediate

2932: .keywords: PetscFE, set, type
2933: .seealso: PetscFEGetType(), PetscFECreate()
2934: @*/
2935: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2936: {
2937:   PetscErrorCode (*r)(PetscFE);
2938:   PetscBool      match;

2943:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
2944:   if (match) return(0);

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

2950:   if (fem->ops->destroy) {
2951:     (*fem->ops->destroy)(fem);
2952:     fem->ops->destroy = NULL;
2953:   }
2954:   (*r)(fem);
2955:   PetscObjectChangeTypeName((PetscObject) fem, name);
2956:   return(0);
2957: }

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

2964:   Not Collective

2966:   Input Parameter:
2967: . fem  - The PetscFE

2969:   Output Parameter:
2970: . name - The PetscFE type name

2972:   Level: intermediate

2974: .keywords: PetscFE, get, type, name
2975: .seealso: PetscFESetType(), PetscFECreate()
2976: @*/
2977: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2978: {

2984:   if (!PetscFERegisterAllCalled) {
2985:     PetscFERegisterAll();
2986:   }
2987:   *name = ((PetscObject) fem)->type_name;
2988:   return(0);
2989: }

2993: /*@C
2994:   PetscFEView - Views a PetscFE

2996:   Collective on PetscFE

2998:   Input Parameter:
2999: + fem - the PetscFE object to view
3000: - v   - the viewer

3002:   Level: developer

3004: .seealso PetscFEDestroy()
3005: @*/
3006: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
3007: {

3012:   if (!v) {
3013:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
3014:   }
3015:   if (fem->ops->view) {
3016:     (*fem->ops->view)(fem, v);
3017:   }
3018:   return(0);
3019: }

3023: /*@
3024:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

3026:   Collective on PetscFE

3028:   Input Parameter:
3029: . fem - the PetscFE object to set options for

3031:   Options Database:
3032: . -petscfe_num_blocks  the number of cell blocks to integrate concurrently
3033: . -petscfe_num_batches the number of cell batches to integrate serially

3035:   Level: developer

3037: .seealso PetscFEView()
3038: @*/
3039: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
3040: {
3041:   const char    *defaultType;
3042:   char           name[256];
3043:   PetscBool      flg;

3048:   if (!((PetscObject) fem)->type_name) {
3049:     defaultType = PETSCFEBASIC;
3050:   } else {
3051:     defaultType = ((PetscObject) fem)->type_name;
3052:   }
3053:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

3055:   PetscObjectOptionsBegin((PetscObject) fem);
3056:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
3057:   if (flg) {
3058:     PetscFESetType(fem, name);
3059:   } else if (!((PetscObject) fem)->type_name) {
3060:     PetscFESetType(fem, defaultType);
3061:   }
3062:   PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
3063:   PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
3064:   if (fem->ops->setfromoptions) {
3065:     (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
3066:   }
3067:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
3068:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
3069:   PetscOptionsEnd();
3070:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
3071:   return(0);
3072: }

3076: /*@C
3077:   PetscFESetUp - Construct data structures for the PetscFE

3079:   Collective on PetscFE

3081:   Input Parameter:
3082: . fem - the PetscFE object to setup

3084:   Level: developer

3086: .seealso PetscFEView(), PetscFEDestroy()
3087: @*/
3088: PetscErrorCode PetscFESetUp(PetscFE fem)
3089: {

3094:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
3095:   return(0);
3096: }

3100: /*@
3101:   PetscFEDestroy - Destroys a PetscFE object

3103:   Collective on PetscFE

3105:   Input Parameter:
3106: . fem - the PetscFE object to destroy

3108:   Level: developer

3110: .seealso PetscFEView()
3111: @*/
3112: PetscErrorCode PetscFEDestroy(PetscFE *fem)
3113: {

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

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

3123:   PetscFree((*fem)->numDof);
3124:   PetscFree((*fem)->invV);
3125:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
3126:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
3127:   PetscSpaceDestroy(&(*fem)->basisSpace);
3128:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
3129:   PetscQuadratureDestroy(&(*fem)->quadrature);

3131:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
3132:   PetscHeaderDestroy(fem);
3133:   return(0);
3134: }

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

3141:   Collective on MPI_Comm

3143:   Input Parameter:
3144: . comm - The communicator for the PetscFE object

3146:   Output Parameter:
3147: . fem - The PetscFE object

3149:   Level: beginner

3151: .seealso: PetscFESetType(), PETSCFEGALERKIN
3152: @*/
3153: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
3154: {
3155:   PetscFE        f;

3160:   PetscCitationsRegister(FECitation,&FEcite);
3161:   *fem = NULL;
3162:   PetscFEInitializePackage();

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

3166:   f->basisSpace    = NULL;
3167:   f->dualSpace     = NULL;
3168:   f->numComponents = 1;
3169:   f->numDof        = NULL;
3170:   f->invV          = NULL;
3171:   f->B             = NULL;
3172:   f->D             = NULL;
3173:   f->H             = NULL;
3174:   PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
3175:   f->blockSize     = 0;
3176:   f->numBlocks     = 1;
3177:   f->batchSize     = 0;
3178:   f->numBatches    = 1;

3180:   *fem = f;
3181:   return(0);
3182: }

3186: /*@
3187:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

3189:   Not collective

3191:   Input Parameter:
3192: . fem - The PetscFE object

3194:   Output Parameter:
3195: . dim - The spatial dimension

3197:   Level: intermediate

3199: .seealso: PetscFECreate()
3200: @*/
3201: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
3202: {
3203:   DM             dm;

3209:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3210:   DMGetDimension(dm, dim);
3211:   return(0);
3212: }

3216: /*@
3217:   PetscFESetNumComponents - Sets the number of components in the element

3219:   Not collective

3221:   Input Parameters:
3222: + fem - The PetscFE object
3223: - comp - The number of field components

3225:   Level: intermediate

3227: .seealso: PetscFECreate()
3228: @*/
3229: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
3230: {
3233:   fem->numComponents = comp;
3234:   return(0);
3235: }

3239: /*@
3240:   PetscFEGetNumComponents - Returns the number of components in the element

3242:   Not collective

3244:   Input Parameter:
3245: . fem - The PetscFE object

3247:   Output Parameter:
3248: . comp - The number of field components

3250:   Level: intermediate

3252: .seealso: PetscFECreate()
3253: @*/
3254: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
3255: {
3259:   *comp = fem->numComponents;
3260:   return(0);
3261: }

3265: /*@
3266:   PetscFESetTileSizes - Sets the tile sizes for evaluation

3268:   Not collective

3270:   Input Parameters:
3271: + fem - The PetscFE object
3272: . blockSize - The number of elements in a block
3273: . numBlocks - The number of blocks in a batch
3274: . batchSize - The number of elements in a batch
3275: - numBatches - The number of batches in a chunk

3277:   Level: intermediate

3279: .seealso: PetscFECreate()
3280: @*/
3281: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
3282: {
3285:   fem->blockSize  = blockSize;
3286:   fem->numBlocks  = numBlocks;
3287:   fem->batchSize  = batchSize;
3288:   fem->numBatches = numBatches;
3289:   return(0);
3290: }

3294: /*@
3295:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

3297:   Not collective

3299:   Input Parameter:
3300: . fem - The PetscFE object

3302:   Output Parameters:
3303: + blockSize - The number of elements in a block
3304: . numBlocks - The number of blocks in a batch
3305: . batchSize - The number of elements in a batch
3306: - numBatches - The number of batches in a chunk

3308:   Level: intermediate

3310: .seealso: PetscFECreate()
3311: @*/
3312: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
3313: {
3320:   if (blockSize)  *blockSize  = fem->blockSize;
3321:   if (numBlocks)  *numBlocks  = fem->numBlocks;
3322:   if (batchSize)  *batchSize  = fem->batchSize;
3323:   if (numBatches) *numBatches = fem->numBatches;
3324:   return(0);
3325: }

3329: /*@
3330:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

3332:   Not collective

3334:   Input Parameter:
3335: . fem - The PetscFE object

3337:   Output Parameter:
3338: . sp - The PetscSpace object

3340:   Level: intermediate

3342: .seealso: PetscFECreate()
3343: @*/
3344: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
3345: {
3349:   *sp = fem->basisSpace;
3350:   return(0);
3351: }

3355: /*@
3356:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

3358:   Not collective

3360:   Input Parameters:
3361: + fem - The PetscFE object
3362: - sp - The PetscSpace object

3364:   Level: intermediate

3366: .seealso: PetscFECreate()
3367: @*/
3368: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
3369: {

3375:   PetscSpaceDestroy(&fem->basisSpace);
3376:   fem->basisSpace = sp;
3377:   PetscObjectReference((PetscObject) fem->basisSpace);
3378:   return(0);
3379: }

3383: /*@
3384:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

3386:   Not collective

3388:   Input Parameter:
3389: . fem - The PetscFE object

3391:   Output Parameter:
3392: . sp - The PetscDualSpace object

3394:   Level: intermediate

3396: .seealso: PetscFECreate()
3397: @*/
3398: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
3399: {
3403:   *sp = fem->dualSpace;
3404:   return(0);
3405: }

3409: /*@
3410:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

3412:   Not collective

3414:   Input Parameters:
3415: + fem - The PetscFE object
3416: - sp - The PetscDualSpace object

3418:   Level: intermediate

3420: .seealso: PetscFECreate()
3421: @*/
3422: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
3423: {

3429:   PetscDualSpaceDestroy(&fem->dualSpace);
3430:   fem->dualSpace = sp;
3431:   PetscObjectReference((PetscObject) fem->dualSpace);
3432:   return(0);
3433: }

3437: /*@
3438:   PetscFEGetQuadrature - Returns the PetscQuadreture used to calculate inner products

3440:   Not collective

3442:   Input Parameter:
3443: . fem - The PetscFE object

3445:   Output Parameter:
3446: . q - The PetscQuadrature object

3448:   Level: intermediate

3450: .seealso: PetscFECreate()
3451: @*/
3452: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
3453: {
3457:   *q = fem->quadrature;
3458:   return(0);
3459: }

3463: /*@
3464:   PetscFESetQuadrature - Sets the PetscQuadreture used to calculate inner products

3466:   Not collective

3468:   Input Parameters:
3469: + fem - The PetscFE object
3470: - q - The PetscQuadrature object

3472:   Level: intermediate

3474: .seealso: PetscFECreate()
3475: @*/
3476: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
3477: {

3482:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
3483:   PetscQuadratureDestroy(&fem->quadrature);
3484:   fem->quadrature = q;
3485:   PetscObjectReference((PetscObject) q);
3486:   return(0);
3487: }

3491: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
3492: {
3493:   const PetscInt *numDofDual;
3494:   PetscErrorCode  ierr;

3499:   PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
3500:   if (!fem->numDof) {
3501:     DM       dm;
3502:     PetscInt dim, d;

3504:     PetscDualSpaceGetDM(fem->dualSpace, &dm);
3505:     DMGetDimension(dm, &dim);
3506:     PetscMalloc1(dim+1, &fem->numDof);
3507:     for (d = 0; d <= dim; ++d) {
3508:       fem->numDof[d] = fem->numComponents*numDofDual[d];
3509:     }
3510:   }
3511:   *numDof = fem->numDof;
3512:   return(0);
3513: }

3517: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
3518: {
3519:   PetscInt         npoints;
3520:   const PetscReal *points;
3521:   PetscErrorCode   ierr;

3528:   PetscQuadratureGetData(fem->quadrature, NULL, &npoints, &points, NULL);
3529:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
3530:   if (B) *B = fem->B;
3531:   if (D) *D = fem->D;
3532:   if (H) *H = fem->H;
3533:   return(0);
3534: }

3538: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **F)
3539: {
3540:   PetscErrorCode   ierr;

3545:   if (!fem->F) {
3546:     PetscDualSpace  sp;
3547:     DM              dm;
3548:     const PetscInt *cone;
3549:     PetscReal      *centroids;
3550:     PetscInt        dim, numFaces, f;

3552:     PetscFEGetDualSpace(fem, &sp);
3553:     PetscDualSpaceGetDM(sp, &dm);
3554:     DMGetDimension(dm, &dim);
3555:     DMPlexGetConeSize(dm, 0, &numFaces);
3556:     DMPlexGetCone(dm, 0, &cone);
3557:     PetscMalloc1(numFaces*dim, &centroids);
3558:     for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, &centroids[f*dim], NULL);}
3559:     PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
3560:     PetscFree(centroids);
3561:   }
3562:   *F = fem->F;
3563:   return(0);
3564: }

3568: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3569: {
3570:   DM               dm;
3571:   PetscInt         pdim; /* Dimension of FE space P */
3572:   PetscInt         dim;  /* Spatial dimension */
3573:   PetscInt         comp; /* Field components */
3574:   PetscErrorCode   ierr;

3582:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3583:   DMGetDimension(dm, &dim);
3584:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3585:   PetscFEGetNumComponents(fem, &comp);
3586:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
3587:   if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
3588:   if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
3589:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
3590:   return(0);
3591: }

3595: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3596: {
3597:   DM             dm;

3602:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3603:   if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
3604:   if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
3605:   if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
3606:   return(0);
3607: }

3611: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
3612: {
3613:   PetscFE_Basic *b = (PetscFE_Basic *) fem->data;

3617:   PetscFree(b);
3618:   return(0);
3619: }

3623: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
3624: {
3625:   PetscSpace        basis;
3626:   PetscDualSpace    dual;
3627:   PetscQuadrature   q = NULL;
3628:   PetscInt          dim, Nq;
3629:   PetscViewerFormat format;
3630:   PetscErrorCode    ierr;

3633:   PetscFEGetBasisSpace(fe, &basis);
3634:   PetscFEGetDualSpace(fe, &dual);
3635:   PetscFEGetQuadrature(fe, &q);
3636:   PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
3637:   PetscViewerGetFormat(viewer, &format);
3638:   PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
3639:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
3640:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
3641:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
3642:     PetscViewerASCIIPushTab(viewer);
3643:     PetscQuadratureView(q, viewer);
3644:     PetscViewerASCIIPopTab(viewer);
3645:   } else {
3646:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
3647:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
3648:   }
3649:   PetscViewerASCIIPushTab(viewer);
3650:   PetscSpaceView(basis, viewer);
3651:   PetscDualSpaceView(dual, viewer);
3652:   PetscViewerASCIIPopTab(viewer);
3653:   return(0);
3654: }

3658: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
3659: {
3660:   PetscBool      iascii;

3666:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
3667:   if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
3668:   return(0);
3669: }

3673: /* Construct the change of basis from prime basis to nodal basis */
3674: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
3675: {
3676:   PetscScalar   *work, *invVscalar;
3677:   PetscBLASInt  *pivots;
3678:   PetscBLASInt   n, info;
3679:   PetscInt       pdim, j;

3683:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3684:   PetscMalloc1(pdim*pdim,&fem->invV);
3685: #if defined(PETSC_USE_COMPLEX)
3686:   PetscMalloc1(pdim*pdim,&invVscalar);
3687: #else
3688:   invVscalar = fem->invV;
3689: #endif
3690:   for (j = 0; j < pdim; ++j) {
3691:     PetscReal      *Bf;
3692:     PetscQuadrature f;
3693:     PetscInt        q, k;

3695:     PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
3696:     PetscMalloc1(f->numPoints*pdim,&Bf);
3697:     PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
3698:     for (k = 0; k < pdim; ++k) {
3699:       /* n_j \cdot \phi_k */
3700:       invVscalar[j*pdim+k] = 0.0;
3701:       for (q = 0; q < f->numPoints; ++q) {
3702:         invVscalar[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
3703:       }
3704:     }
3705:     PetscFree(Bf);
3706:   }
3707:   PetscMalloc2(pdim,&pivots,pdim,&work);
3708:   n = pdim;
3709:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info));
3710:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info));
3711: #if defined(PETSC_USE_COMPLEX)
3712:   for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]);
3713:   PetscFree(invVscalar);
3714: #endif
3715:   PetscFree2(pivots,work);
3716:   return(0);
3717: }

3721: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
3722: {

3726:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
3727:   return(0);
3728: }

3732: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
3733: {
3734:   DM               dm;
3735:   PetscInt         pdim; /* Dimension of FE space P */
3736:   PetscInt         dim;  /* Spatial dimension */
3737:   PetscInt         comp; /* Field components */
3738:   PetscReal       *tmpB, *tmpD, *tmpH;
3739:   PetscInt         p, d, j, k;
3740:   PetscErrorCode   ierr;

3743:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3744:   DMGetDimension(dm, &dim);
3745:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3746:   PetscFEGetNumComponents(fem, &comp);
3747:   /* Evaluate the prime basis functions at all points */
3748:   if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3749:   if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3750:   if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3751:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
3752:   /* Translate to the nodal basis */
3753:   for (p = 0; p < npoints; ++p) {
3754:     if (B) {
3755:       /* Multiply by V^{-1} (pdim x pdim) */
3756:       for (j = 0; j < pdim; ++j) {
3757:         const PetscInt i = (p*pdim + j)*comp;
3758:         PetscInt       c;

3760:         B[i] = 0.0;
3761:         for (k = 0; k < pdim; ++k) {
3762:           B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
3763:         }
3764:         for (c = 1; c < comp; ++c) {
3765:           B[i+c] = B[i];
3766:         }
3767:       }
3768:     }
3769:     if (D) {
3770:       /* Multiply by V^{-1} (pdim x pdim) */
3771:       for (j = 0; j < pdim; ++j) {
3772:         for (d = 0; d < dim; ++d) {
3773:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
3774:           PetscInt       c;

3776:           D[i] = 0.0;
3777:           for (k = 0; k < pdim; ++k) {
3778:             D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
3779:           }
3780:           for (c = 1; c < comp; ++c) {
3781:             D[((p*pdim + j)*comp + c)*dim + d] = D[i];
3782:           }
3783:         }
3784:       }
3785:     }
3786:     if (H) {
3787:       /* Multiply by V^{-1} (pdim x pdim) */
3788:       for (j = 0; j < pdim; ++j) {
3789:         for (d = 0; d < dim*dim; ++d) {
3790:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
3791:           PetscInt       c;

3793:           H[i] = 0.0;
3794:           for (k = 0; k < pdim; ++k) {
3795:             H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
3796:           }
3797:           for (c = 1; c < comp; ++c) {
3798:             H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
3799:           }
3800:         }
3801:       }
3802:     }
3803:   }
3804:   if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3805:   if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3806:   if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3807:   return(0);
3808: }

3812: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3813:                                       const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
3814: {
3815:   const PetscInt  debug = 0;
3816:   PetscPointFunc  obj_func;
3817:   PetscQuadrature quad;
3818:   PetscScalar    *u, *u_x, *a, *a_x;
3819:   PetscReal      *x;
3820:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3821:   PetscInt        dim, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
3822:   PetscErrorCode  ierr;

3825:   PetscDSGetObjective(prob, field, &obj_func);
3826:   if (!obj_func) return(0);
3827:   PetscFEGetSpatialDimension(fem, &dim);
3828:   PetscFEGetQuadrature(fem, &quad);
3829:   PetscDSGetNumFields(prob, &Nf);
3830:   PetscDSGetTotalDimension(prob, &totDim);
3831:   PetscDSGetComponentOffsets(prob, &uOff);
3832:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3833:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
3834:   PetscDSGetRefCoordArrays(prob, &x, NULL);
3835:   if (probAux) {
3836:     PetscDSGetNumFields(probAux, &NfAux);
3837:     PetscDSGetTotalDimension(probAux, &totDimAux);
3838:     PetscDSGetComponentOffsets(probAux, &aOff);
3839:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3840:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3841:   }
3842:   for (e = 0; e < Ne; ++e) {
3843:     const PetscReal *v0   = geom[e].v0;
3844:     const PetscReal *J    = geom[e].J;
3845:     const PetscReal *invJ = geom[e].invJ;
3846:     const PetscReal  detJ = geom[e].detJ;
3847:     const PetscReal *quadPoints, *quadWeights;
3848:     PetscInt         Nq, q;

3850:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3851:     if (debug > 1) {
3852:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3853: #ifndef PETSC_USE_COMPLEX
3854:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3855: #endif
3856:     }
3857:     for (q = 0; q < Nq; ++q) {
3858:       PetscScalar integrand;

3860:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3861:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3862:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       NULL, u, u_x, NULL);
3863:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3864:       obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &integrand);
3865:       integrand *= detJ*quadWeights[q];
3866:       integral[field] += PetscRealPart(integrand);
3867:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", PetscRealPart(integrand), integral[field]);}
3868:     }
3869:     cOffset    += totDim;
3870:     cOffsetAux += totDimAux;
3871:   }
3872:   return(0);
3873: }

3877: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3878:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
3879: {
3880:   const PetscInt  debug = 0;
3881:   PetscPointFunc  f0_func;
3882:   PetscPointFunc  f1_func;
3883:   PetscQuadrature quad;
3884:   PetscReal     **basisField, **basisFieldDer;
3885:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3886:   PetscReal      *x;
3887:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3888:   PetscInt        dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3889:   PetscErrorCode  ierr;

3892:   PetscFEGetSpatialDimension(fem, &dim);
3893:   PetscFEGetQuadrature(fem, &quad);
3894:   PetscFEGetDimension(fem, &Nb);
3895:   PetscFEGetNumComponents(fem, &Nc);
3896:   PetscDSGetNumFields(prob, &Nf);
3897:   PetscDSGetTotalDimension(prob, &totDim);
3898:   PetscDSGetComponentOffsets(prob, &uOff);
3899:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3900:   PetscDSGetFieldOffset(prob, field, &fOffset);
3901:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3902:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3903:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3904:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3905:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3906:   if (probAux) {
3907:     PetscDSGetNumFields(probAux, &NfAux);
3908:     PetscDSGetTotalDimension(probAux, &totDimAux);
3909:     PetscDSGetComponentOffsets(probAux, &aOff);
3910:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3911:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3912:   }
3913:   for (e = 0; e < Ne; ++e) {
3914:     const PetscReal *v0   = geom[e].v0;
3915:     const PetscReal *J    = geom[e].J;
3916:     const PetscReal *invJ = geom[e].invJ;
3917:     const PetscReal  detJ = geom[e].detJ;
3918:     const PetscReal *quadPoints, *quadWeights;
3919:     PetscInt         Nq, q;

3921:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3922:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3923:     PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
3924:     if (debug > 1) {
3925:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3926: #ifndef PETSC_USE_COMPLEX
3927:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3928: #endif
3929:     }
3930:     for (q = 0; q < Nq; ++q) {
3931:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3932:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3933:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3934:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3935:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, &f0[q*Nc]);
3936:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, refSpaceDer);
3937:       TransformF(dim, dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3938:     }
3939:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3940:     cOffset    += totDim;
3941:     cOffsetAux += totDimAux;
3942:   }
3943:   return(0);
3944: }

3948: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3949:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
3950: {
3951:   const PetscInt  debug = 0;
3952:   PetscBdPointFunc f0_func;
3953:   PetscBdPointFunc f1_func;
3954:   PetscQuadrature quad;
3955:   PetscReal     **basisField, **basisFieldDer;
3956:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3957:   PetscReal      *x;
3958:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3959:   PetscInt        dim, Nf, NfAux = 0, bdim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3960:   PetscErrorCode  ierr;

3963:   PetscFEGetSpatialDimension(fem, &dim);
3964:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3965:   bdim = dim-1;
3966:   PetscFEGetQuadrature(fem, &quad);
3967:   PetscFEGetDimension(fem, &Nb);
3968:   PetscFEGetNumComponents(fem, &Nc);
3969:   PetscDSGetNumFields(prob, &Nf);
3970:   PetscDSGetTotalBdDimension(prob, &totDim);
3971:   PetscDSGetComponentBdOffsets(prob, &uOff);
3972:   PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
3973:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3974:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3975:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3976:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3977:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3978:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3979:   if (probAux) {
3980:     PetscDSGetNumFields(probAux, &NfAux);
3981:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3982:     PetscDSGetComponentBdOffsets(probAux, &aOff);
3983:     PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
3984:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3985:   }
3986:   for (e = 0; e < Ne; ++e) {
3987:     const PetscReal *v0          = geom[e].v0;
3988:     const PetscReal *n           = geom[e].n;
3989:     const PetscReal *J           = geom[e].J;
3990:     const PetscReal *invJ        = geom[e].invJ;
3991:     const PetscReal  detJ        = geom[e].detJ;
3992:     const PetscReal *quadPoints, *quadWeights;
3993:     PetscInt         Nq, q;

3995:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3996:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3997:     PetscMemzero(f1, Nq*Nc*bdim * sizeof(PetscScalar));
3998:     if (debug > 1) {
3999:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
4000: #ifndef PETSC_USE_COMPLEX
4001:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
4002: #endif
4003:     }
4004:     for (q = 0; q < Nq; ++q) {
4005:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4006:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
4007:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4008:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4009:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, &f0[q*Nc]);
4010:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, refSpaceDer);
4011:       TransformF(dim, bdim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4012:     }
4013:     UpdateElementVec(bdim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4014:     cOffset    += totDim;
4015:     cOffsetAux += totDimAux;
4016:   }
4017:   return(0);
4018: }

4022: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
4023:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
4024: {
4025:   const PetscInt  debug      = 0;
4026:   PetscPointJac   g0_func;
4027:   PetscPointJac   g1_func;
4028:   PetscPointJac   g2_func;
4029:   PetscPointJac   g3_func;
4030:   PetscFE         fe;
4031:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
4032:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
4033:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
4034:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
4035:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
4036:   PetscQuadrature quad;
4037:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
4038:   PetscReal      *x;
4039:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
4040:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4041:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
4042:   PetscInt        dim, Nf, NfAux = 0, totDim, totDimAux = 0, e;
4043:   PetscErrorCode  ierr;

4046:   PetscFEGetSpatialDimension(fem, &dim);
4047:   PetscFEGetQuadrature(fem, &quad);
4048:   PetscDSGetNumFields(prob, &Nf);
4049:   PetscDSGetTotalDimension(prob, &totDim);
4050:   PetscDSGetComponentOffsets(prob, &uOff);
4051:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4052:   switch(jtype) {
4053:   case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4054:   case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4055:   case PETSCFE_JACOBIAN:     PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4056:   }
4057:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
4058:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
4059:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4060:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
4061:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4062:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4063:   PetscFEGetDimension(fe, &NbI);
4064:   PetscFEGetNumComponents(fe, &NcI);
4065:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
4066:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
4067:   PetscFEGetDimension(fe, &NbJ);
4068:   PetscFEGetNumComponents(fe, &NcJ);
4069:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
4070:   if (probAux) {
4071:     PetscDSGetNumFields(probAux, &NfAux);
4072:     PetscDSGetTotalDimension(probAux, &totDimAux);
4073:     PetscDSGetComponentOffsets(probAux, &aOff);
4074:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4075:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4076:   }
4077:   basisI    = basisField[fieldI];
4078:   basisJ    = basisField[fieldJ];
4079:   basisDerI = basisFieldDer[fieldI];
4080:   basisDerJ = basisFieldDer[fieldJ];
4081:   /* Initialize here in case the function is not defined */
4082:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4083:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
4084:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
4085:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4086:   for (e = 0; e < Ne; ++e) {
4087:     const PetscReal *v0   = geom[e].v0;
4088:     const PetscReal *J    = geom[e].J;
4089:     const PetscReal *invJ = geom[e].invJ;
4090:     const PetscReal  detJ = geom[e].detJ;
4091:     const PetscReal *quadPoints, *quadWeights;
4092:     PetscInt         Nq, q;

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

4098:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4099:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4100:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4101:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4102:       if (g0_func) {
4103:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4104:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, g0);
4105:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
4106:       }
4107:       if (g1_func) {
4108:         PetscInt d, d2;
4109:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4110:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4111:         for (fc = 0; fc < NcI; ++fc) {
4112:           for (gc = 0; gc < NcJ; ++gc) {
4113:             for (d = 0; d < dim; ++d) {
4114:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
4115:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4116:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4117:             }
4118:           }
4119:         }
4120:       }
4121:       if (g2_func) {
4122:         PetscInt d, d2;
4123:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4124:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4125:         for (fc = 0; fc < NcI; ++fc) {
4126:           for (gc = 0; gc < NcJ; ++gc) {
4127:             for (d = 0; d < dim; ++d) {
4128:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
4129:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4130:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4131:             }
4132:           }
4133:         }
4134:       }
4135:       if (g3_func) {
4136:         PetscInt d, d2, dp, d3;
4137:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4138:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4139:         for (fc = 0; fc < NcI; ++fc) {
4140:           for (gc = 0; gc < NcJ; ++gc) {
4141:             for (d = 0; d < dim; ++d) {
4142:               for (dp = 0; dp < dim; ++dp) {
4143:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
4144:                 for (d2 = 0; d2 < dim; ++d2) {
4145:                   for (d3 = 0; d3 < dim; ++d3) {
4146:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
4147:                   }
4148:                 }
4149:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
4150:               }
4151:             }
4152:           }
4153:         }
4154:       }

4156:       for (f = 0; f < NbI; ++f) {
4157:         for (fc = 0; fc < NcI; ++fc) {
4158:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
4159:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
4160:           for (g = 0; g < NbJ; ++g) {
4161:             for (gc = 0; gc < NcJ; ++gc) {
4162:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
4163:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
4164:               const PetscInt fOff = eOffset+i*totDim+j;
4165:               PetscInt       d, d2;

4167:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
4168:               for (d = 0; d < dim; ++d) {
4169:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
4170:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
4171:                 for (d2 = 0; d2 < dim; ++d2) {
4172:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
4173:                 }
4174:               }
4175:             }
4176:           }
4177:         }
4178:       }
4179:     }
4180:     if (debug > 1) {
4181:       PetscInt fc, f, gc, g;

4183:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
4184:       for (fc = 0; fc < NcI; ++fc) {
4185:         for (f = 0; f < NbI; ++f) {
4186:           const PetscInt i = offsetI + f*NcI+fc;
4187:           for (gc = 0; gc < NcJ; ++gc) {
4188:             for (g = 0; g < NbJ; ++g) {
4189:               const PetscInt j = offsetJ + g*NcJ+gc;
4190:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
4191:             }
4192:           }
4193:           PetscPrintf(PETSC_COMM_SELF, "\n");
4194:         }
4195:       }
4196:     }
4197:     cOffset    += totDim;
4198:     cOffsetAux += totDimAux;
4199:     eOffset    += PetscSqr(totDim);
4200:   }
4201:   return(0);
4202: }

4206: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
4207:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
4208: {
4209:   const PetscInt  debug      = 0;
4210:   PetscBdPointJac g0_func;
4211:   PetscBdPointJac g1_func;
4212:   PetscBdPointJac g2_func;
4213:   PetscBdPointJac g3_func;
4214:   PetscFE         fe;
4215:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
4216:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
4217:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
4218:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
4219:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
4220:   PetscQuadrature quad;
4221:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
4222:   PetscReal      *x;
4223:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
4224:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4225:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
4226:   PetscInt        dim, Nf, NfAux = 0, bdim, totDim, totDimAux = 0, e;
4227:   PetscErrorCode  ierr;

4230:   PetscFEGetQuadrature(fem, &quad);
4231:   PetscFEGetSpatialDimension(fem, &dim);
4232:   dim += 1; /* Spatial dimension is one higher than topological dimension */
4233:   bdim = dim-1;
4234:   PetscDSGetNumFields(prob, &Nf);
4235:   PetscDSGetTotalBdDimension(prob, &totDim);
4236:   PetscDSGetComponentBdOffsets(prob, &uOff);
4237:   PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
4238:   PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
4239:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
4240:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4241:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
4242:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
4243:   PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
4244:   PetscFEGetDimension(fe, &NbI);
4245:   PetscFEGetNumComponents(fe, &NcI);
4246:   PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
4247:   PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
4248:   PetscFEGetDimension(fe, &NbJ);
4249:   PetscFEGetNumComponents(fe, &NcJ);
4250:   PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
4251:   if (probAux) {
4252:     PetscDSGetNumFields(probAux, &NfAux);
4253:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
4254:     PetscDSGetComponentBdOffsets(probAux, &aOff);
4255:     PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
4256:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4257:   }
4258:   basisI    = basisField[fieldI];
4259:   basisJ    = basisField[fieldJ];
4260:   basisDerI = basisFieldDer[fieldI];
4261:   basisDerJ = basisFieldDer[fieldJ];
4262:   /* Initialize here in case the function is not defined */
4263:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4264:   PetscMemzero(g1, NcI*NcJ*bdim * sizeof(PetscScalar));
4265:   PetscMemzero(g2, NcI*NcJ*bdim * sizeof(PetscScalar));
4266:   PetscMemzero(g3, NcI*NcJ*bdim*bdim * sizeof(PetscScalar));
4267:   for (e = 0; e < Ne; ++e) {
4268:     const PetscReal *v0   = geom[e].v0;
4269:     const PetscReal *n    = geom[e].n;
4270:     const PetscReal *J    = geom[e].J;
4271:     const PetscReal *invJ = geom[e].invJ;
4272:     const PetscReal  detJ = geom[e].detJ;
4273:     const PetscReal *quadPoints, *quadWeights;
4274:     PetscInt         Nq, q;

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

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

4338:       for (f = 0; f < NbI; ++f) {
4339:         for (fc = 0; fc < NcI; ++fc) {
4340:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
4341:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
4342:           for (g = 0; g < NbJ; ++g) {
4343:             for (gc = 0; gc < NcJ; ++gc) {
4344:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
4345:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
4346:               const PetscInt fOff = eOffset+i*totDim+j;
4347:               PetscInt       d, d2;

4349:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
4350:               for (d = 0; d < bdim; ++d) {
4351:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
4352:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
4353:                 for (d2 = 0; d2 < bdim; ++d2) {
4354:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
4355:                 }
4356:               }
4357:             }
4358:           }
4359:         }
4360:       }
4361:     }
4362:     if (debug > 1) {
4363:       PetscInt fc, f, gc, g;

4365:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
4366:       for (fc = 0; fc < NcI; ++fc) {
4367:         for (f = 0; f < NbI; ++f) {
4368:           const PetscInt i = offsetI + f*NcI+fc;
4369:           for (gc = 0; gc < NcJ; ++gc) {
4370:             for (g = 0; g < NbJ; ++g) {
4371:               const PetscInt j = offsetJ + g*NcJ+gc;
4372:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
4373:             }
4374:           }
4375:           PetscPrintf(PETSC_COMM_SELF, "\n");
4376:         }
4377:       }
4378:     }
4379:     cOffset    += totDim;
4380:     cOffsetAux += totDimAux;
4381:     eOffset    += PetscSqr(totDim);
4382:   }
4383:   return(0);
4384: }

4388: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
4389: {
4391:   fem->ops->setfromoptions          = NULL;
4392:   fem->ops->setup                   = PetscFESetUp_Basic;
4393:   fem->ops->view                    = PetscFEView_Basic;
4394:   fem->ops->destroy                 = PetscFEDestroy_Basic;
4395:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4396:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
4397:   fem->ops->integrate               = PetscFEIntegrate_Basic;
4398:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
4399:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
4400:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
4401:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
4402:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
4403:   return(0);
4404: }

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

4409:   Level: intermediate

4411: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4412: M*/

4416: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
4417: {
4418:   PetscFE_Basic *b;

4423:   PetscNewLog(fem,&b);
4424:   fem->data = b;

4426:   PetscFEInitialize_Basic(fem);
4427:   return(0);
4428: }

4432: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
4433: {
4434:   PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;

4438:   PetscFree(na);
4439:   return(0);
4440: }

4444: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4445:                                                   const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
4446: {
4447:   const PetscInt  debug = 0;
4448:   PetscPointFunc  f0_func;
4449:   PetscPointFunc  f1_func;
4450:   PetscQuadrature quad;
4451:   PetscReal     **basisField, **basisFieldDer;
4452:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
4453:   PetscReal      *x;
4454:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4455:   PetscInt        dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
4456:   PetscErrorCode  ierr;

4459:   PetscFEGetSpatialDimension(fem, &dim);
4460:   PetscFEGetQuadrature(fem, &quad);
4461:   PetscFEGetDimension(fem, &Nb);
4462:   PetscFEGetNumComponents(fem, &Nc);
4463:   PetscDSGetNumFields(prob, &Nf);
4464:   PetscDSGetTotalDimension(prob, &totDim);
4465:   PetscDSGetComponentOffsets(prob, &uOff);
4466:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4467:   PetscDSGetFieldOffset(prob, field, &fOffset);
4468:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4469:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
4470:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4471:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
4472:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4473:   if (probAux) {
4474:     PetscDSGetNumFields(probAux, &NfAux);
4475:     PetscDSGetTotalDimension(probAux, &totDimAux);
4476:     PetscDSGetComponentOffsets(probAux, &aOff);
4477:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4478:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4479:   }
4480:   for (e = 0; e < Ne; ++e) {
4481:     const PetscReal *quadPoints, *quadWeights;
4482:     PetscInt         Nq, q;

4484:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4485:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
4486:     PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
4487:     for (q = 0; q < Nq; ++q) {
4488:       const PetscReal *v0   = geom[e*Nq+q].v0;
4489:       const PetscReal *J    = geom[e*Nq+q].J;
4490:       const PetscReal *invJ = geom[e*Nq+q].invJ;
4491:       const PetscReal  detJ = geom[e*Nq+q].detJ;

4493:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4494:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4495:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4496:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4497:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, &f0[q*Nc]);
4498:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, refSpaceDer);
4499:       TransformF(dim, dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4500:     }
4501:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4502:     cOffset    += totDim;
4503:     cOffsetAux += totDimAux;
4504:   }
4505:   return(0);
4506: }

4510: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4511:                                                     const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
4512: {
4513:   const PetscInt   debug = 0;
4514:   PetscBdPointFunc f0_func;
4515:   PetscBdPointFunc f1_func;
4516:   PetscQuadrature  quad;
4517:   PetscReal      **basisField, **basisFieldDer;
4518:   PetscScalar     *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4519:   PetscReal       *x;
4520:   PetscInt        *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4521:   PetscInt         dim, Nf, NfAux = 0, bdim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
4522:   PetscErrorCode   ierr;

4525:   PetscFEGetSpatialDimension(fem, &dim);
4526:   dim += 1; /* Spatial dimension is one higher than topological dimension */
4527:   bdim = dim-1;
4528:   PetscFEGetQuadrature(fem, &quad);
4529:   PetscFEGetDimension(fem, &Nb);
4530:   PetscFEGetNumComponents(fem, &Nc);
4531:   PetscDSGetNumFields(prob, &Nf);
4532:   PetscDSGetTotalBdDimension(prob, &totDim);
4533:   PetscDSGetComponentBdOffsets(prob, &uOff);
4534:   PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
4535:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
4536:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
4537:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4538:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4539:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
4540:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
4541:   if (probAux) {
4542:     PetscDSGetNumFields(probAux, &NfAux);
4543:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
4544:     PetscDSGetComponentBdOffsets(probAux, &aOff);
4545:     PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
4546:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4547:   }
4548:   for (e = 0; e < Ne; ++e) {
4549:     const PetscReal *quadPoints, *quadWeights;
4550:     PetscInt         Nq, q;

4552:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4553:     PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
4554:     PetscMemzero(f1, Nq*Nc*bdim * sizeof(PetscScalar));
4555:     for (q = 0; q < Nq; ++q) {
4556:       const PetscReal *v0   = geom[e*Nq+q].v0;
4557:       const PetscReal *n    = geom[e*Nq+q].n;
4558:       const PetscReal *J    = geom[e*Nq+q].J;
4559:       const PetscReal *invJ = geom[e*Nq+q].invJ;
4560:       const PetscReal  detJ = geom[e*Nq+q].detJ;

4562:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4563:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
4564:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4565:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4566:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, &f0[q*Nc]);
4567:       if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, x, n, refSpaceDer);
4568:       TransformF(dim, bdim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4569:     }
4570:     UpdateElementVec(bdim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4571:     cOffset    += totDim;
4572:     cOffsetAux += totDimAux;
4573:   }
4574:   return(0);
4575: }

4579: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
4580:                                                   const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
4581: {
4582:   const PetscInt  debug      = 0;
4583:   PetscPointJac   g0_func;
4584:   PetscPointJac   g1_func;
4585:   PetscPointJac   g2_func;
4586:   PetscPointJac   g3_func;
4587:   PetscFE         fe;
4588:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
4589:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
4590:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
4591:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
4592:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
4593:   PetscQuadrature quad;
4594:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4595:   PetscReal      *x;
4596:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
4597:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
4598:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4599:   PetscInt        dim, Nf, NfAux = 0, totDim, totDimAux = 0, e;
4600:   PetscErrorCode  ierr;

4603:   PetscFEGetSpatialDimension(fem, &dim);
4604:   PetscFEGetQuadrature(fem, &quad);
4605:   PetscDSGetNumFields(prob, &Nf);
4606:   PetscDSGetTotalDimension(prob, &totDim);
4607:   PetscDSGetComponentOffsets(prob, &uOff);
4608:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4609:   switch(jtype) {
4610:   case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4611:   case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4612:   case PETSCFE_JACOBIAN:     PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4613:   }
4614:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4615:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4616:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
4617:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4618:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4619:   PetscFEGetDimension(fe, &NbI);
4620:   PetscFEGetNumComponents(fe, &NcI);
4621:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
4622:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
4623:   PetscFEGetDimension(fe, &NbJ);
4624:   PetscFEGetNumComponents(fe, &NcJ);
4625:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
4626:   if (probAux) {
4627:     PetscDSGetNumFields(probAux, &NfAux);
4628:     PetscDSGetTotalDimension(probAux, &totDimAux);
4629:     PetscDSGetComponentOffsets(probAux, &aOff);
4630:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4631:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4632:   }
4633:   basisI    = basisField[fieldI];
4634:   basisJ    = basisField[fieldJ];
4635:   basisDerI = basisFieldDer[fieldI];
4636:   basisDerJ = basisFieldDer[fieldJ];
4637:   /* Initialize here in case the function is not defined */
4638:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4639:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
4640:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
4641:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4642:   for (e = 0; e < Ne; ++e) {
4643:     const PetscReal *quadPoints, *quadWeights;
4644:     PetscInt         Nq, q;

4646:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4647:     for (q = 0; q < Nq; ++q) {
4648:       const PetscReal *v0   = geom[e*Nq+q].v0;
4649:       const PetscReal *J    = geom[e*Nq+q].J;
4650:       const PetscReal *invJ = geom[e*Nq+q].invJ;
4651:       const PetscReal  detJ = geom[e*Nq+q].detJ;
4652:       PetscInt         f, g, fc, gc, c;

4654:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4655:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4656:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4657:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4658:       if (g0_func) {
4659:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4660:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, g0);
4661:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
4662:       }
4663:       if (g1_func) {
4664:         PetscInt d, d2;
4665:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4666:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4667:         for (fc = 0; fc < NcI; ++fc) {
4668:           for (gc = 0; gc < NcJ; ++gc) {
4669:             for (d = 0; d < dim; ++d) {
4670:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
4671:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4672:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4673:             }
4674:           }
4675:         }
4676:       }
4677:       if (g2_func) {
4678:         PetscInt d, d2;
4679:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4680:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4681:         for (fc = 0; fc < NcI; ++fc) {
4682:           for (gc = 0; gc < NcJ; ++gc) {
4683:             for (d = 0; d < dim; ++d) {
4684:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
4685:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4686:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4687:             }
4688:           }
4689:         }
4690:       }
4691:       if (g3_func) {
4692:         PetscInt d, d2, dp, d3;
4693:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4694:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, x, refSpaceDer);
4695:         for (fc = 0; fc < NcI; ++fc) {
4696:           for (gc = 0; gc < NcJ; ++gc) {
4697:             for (d = 0; d < dim; ++d) {
4698:               for (dp = 0; dp < dim; ++dp) {
4699:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
4700:                 for (d2 = 0; d2 < dim; ++d2) {
4701:                   for (d3 = 0; d3 < dim; ++d3) {
4702:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
4703:                   }
4704:                 }
4705:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
4706:               }
4707:             }
4708:           }
4709:         }
4710:       }

4712:       for (f = 0; f < NbI; ++f) {
4713:         for (fc = 0; fc < NcI; ++fc) {
4714:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
4715:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
4716:           for (g = 0; g < NbJ; ++g) {
4717:             for (gc = 0; gc < NcJ; ++gc) {
4718:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
4719:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
4720:               const PetscInt fOff = eOffset+i*totDim+j;
4721:               PetscInt       d, d2;

4723:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
4724:               for (d = 0; d < dim; ++d) {
4725:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
4726:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
4727:                 for (d2 = 0; d2 < dim; ++d2) {
4728:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
4729:                 }
4730:               }
4731:             }
4732:           }
4733:         }
4734:       }
4735:     }
4736:     if (debug > 1) {
4737:       PetscInt fc, f, gc, g;

4739:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
4740:       for (fc = 0; fc < NcI; ++fc) {
4741:         for (f = 0; f < NbI; ++f) {
4742:           const PetscInt i = offsetI + f*NcI+fc;
4743:           for (gc = 0; gc < NcJ; ++gc) {
4744:             for (g = 0; g < NbJ; ++g) {
4745:               const PetscInt j = offsetJ + g*NcJ+gc;
4746:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
4747:             }
4748:           }
4749:           PetscPrintf(PETSC_COMM_SELF, "\n");
4750:         }
4751:       }
4752:     }
4753:     cOffset    += totDim;
4754:     cOffsetAux += totDimAux;
4755:     eOffset    += PetscSqr(totDim);
4756:   }
4757:   return(0);
4758: }

4762: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
4763: {
4765:   fem->ops->setfromoptions          = NULL;
4766:   fem->ops->setup                   = PetscFESetUp_Basic;
4767:   fem->ops->view                    = NULL;
4768:   fem->ops->destroy                 = PetscFEDestroy_Nonaffine;
4769:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4770:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
4771:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Nonaffine;
4772:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Nonaffine;
4773:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
4774:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Nonaffine;
4775:   return(0);
4776: }

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

4781:   Level: intermediate

4783: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4784: M*/

4788: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
4789: {
4790:   PetscFE_Nonaffine *na;
4791:   PetscErrorCode     ierr;

4795:   PetscNewLog(fem, &na);
4796:   fem->data = na;

4798:   PetscFEInitialize_Nonaffine(fem);
4799:   return(0);
4800: }

4802: #ifdef PETSC_HAVE_OPENCL

4806: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
4807: {
4808:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4809:   PetscErrorCode  ierr;

4812:   clReleaseCommandQueue(ocl->queue_id);
4813:   ocl->queue_id = 0;
4814:   clReleaseContext(ocl->ctx_id);
4815:   ocl->ctx_id = 0;
4816:   PetscFree(ocl);
4817:   return(0);
4818: }

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

4825: /* dim     Number of spatial dimensions:          2                   */
4826: /* N_b     Number of basis functions:             generated           */
4827: /* N_{bt}  Number of total basis functions:       N_b * N_{comp}      */
4828: /* N_q     Number of quadrature points:           generated           */
4829: /* N_{bs}  Number of block cells                  LCM(N_b, N_q)       */
4830: /* N_{bst} Number of block cell components        LCM(N_{bt}, N_q)    */
4831: /* N_{bl}  Number of concurrent blocks            generated           */
4832: /* N_t     Number of threads:                     N_{bl} * N_{bs}     */
4833: /* N_{cbc} Number of concurrent basis      cells: N_{bl} * N_q        */
4834: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b        */
4835: /* N_{sbc} Number of serial     basis      cells: N_{bs} / N_q        */
4836: /* N_{sqc} Number of serial     quadrature cells: N_{bs} / N_b        */
4837: /* N_{cb}  Number of serial cell batches:         input               */
4838: /* N_c     Number of total cells:                 N_{cb}*N_{t}/N_{comp} */
4839: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
4840: {
4841:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4842:   PetscQuadrature q;
4843:   char           *string_tail   = *string_buffer;
4844:   char           *end_of_buffer = *string_buffer + buffer_length;
4845:   char            float_str[]   = "float", double_str[]  = "double";
4846:   char           *numeric_str   = &(float_str[0]);
4847:   PetscInt        op            = ocl->op;
4848:   PetscBool       useField      = PETSC_FALSE;
4849:   PetscBool       useFieldDer   = PETSC_TRUE;
4850:   PetscBool       useFieldAux   = useAux;
4851:   PetscBool       useFieldDerAux= PETSC_FALSE;
4852:   PetscBool       useF0         = PETSC_TRUE;
4853:   PetscBool       useF1         = PETSC_TRUE;
4854:   PetscReal      *basis, *basisDer;
4855:   PetscInt        dim, N_b, N_c, N_q, N_t, p, d, b, c;
4856:   size_t          count;
4857:   PetscErrorCode  ierr;

4860:   PetscFEGetSpatialDimension(fem, &dim);
4861:   PetscFEGetDimension(fem, &N_b);
4862:   PetscFEGetNumComponents(fem, &N_c);
4863:   PetscFEGetQuadrature(fem, &q);
4864:   N_q  = q->numPoints;
4865:   N_t  = N_b * N_c * N_q * N_bl;
4866:   /* Enable device extension for double precision */
4867:   if (ocl->realType == PETSC_DOUBLE) {
4868:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4869: "#if defined(cl_khr_fp64)\n"
4870: "#  pragma OPENCL EXTENSION cl_khr_fp64: enable\n"
4871: "#elif defined(cl_amd_fp64)\n"
4872: "#  pragma OPENCL EXTENSION cl_amd_fp64: enable\n"
4873: "#endif\n",
4874:                               &count);STRING_ERROR_CHECK("Message to short");
4875:     numeric_str  = &(double_str[0]);
4876:   }
4877:   /* Kernel API */
4878:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4879: "\n"
4880: "__kernel void integrateElementQuadrature(int N_cb, __global %s *coefficients, __global %s *coefficientsAux, __global %s *jacobianInverses, __global %s *jacobianDeterminants, __global %s *elemVec)\n"
4881: "{\n",
4882:                        &count, numeric_str, numeric_str, numeric_str, numeric_str, numeric_str);STRING_ERROR_CHECK("Message to short");
4883:   /* Quadrature */
4884:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4885: "  /* Quadrature points\n"
4886: "   - (x1,y1,x2,y2,...) */\n"
4887: "  const %s points[%d] = {\n",
4888:                        &count, numeric_str, N_q*dim);STRING_ERROR_CHECK("Message to short");
4889:   for (p = 0; p < N_q; ++p) {
4890:     for (d = 0; d < dim; ++d) {
4891:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->points[p*dim+d]);STRING_ERROR_CHECK("Message to short");
4892:     }
4893:   }
4894:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4895:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4896: "  /* Quadrature weights\n"
4897: "   - (v1,v2,...) */\n"
4898: "  const %s weights[%d] = {\n",
4899:                        &count, numeric_str, N_q);STRING_ERROR_CHECK("Message to short");
4900:   for (p = 0; p < N_q; ++p) {
4901:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->weights[p]);STRING_ERROR_CHECK("Message to short");
4902:   }
4903:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4904:   /* Basis Functions */
4905:   PetscFEGetDefaultTabulation(fem, &basis, &basisDer, NULL);
4906:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4907: "  /* Nodal basis function evaluations\n"
4908: "    - basis component is fastest varying, the basis function, then point */\n"
4909: "  const %s Basis[%d] = {\n",
4910:                        &count, numeric_str, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4911:   for (p = 0; p < N_q; ++p) {
4912:     for (b = 0; b < N_b; ++b) {
4913:       for (c = 0; c < N_c; ++c) {
4914:         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");
4915:       }
4916:     }
4917:   }
4918:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4919:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4920: "\n"
4921: "  /* Nodal basis function derivative evaluations,\n"
4922: "      - derivative direction is fastest varying, then basis component, then basis function, then point */\n"
4923: "  const %s%d BasisDerivatives[%d] = {\n",
4924:                        &count, numeric_str, dim, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4925:   for (p = 0; p < N_q; ++p) {
4926:     for (b = 0; b < N_b; ++b) {
4927:       for (c = 0; c < N_c; ++c) {
4928:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "(%s%d)(", &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4929:         for (d = 0; d < dim; ++d) {
4930:           if (d > 0) {
4931:             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");
4932:           } else {
4933:             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");
4934:           }
4935:         }
4936:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "),\n", &count);STRING_ERROR_CHECK("Message to short");
4937:       }
4938:     }
4939:   }
4940:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4941:   /* Sizes */
4942:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4943: "  const int dim    = %d;                           // The spatial dimension\n"
4944: "  const int N_bl   = %d;                           // The number of concurrent blocks\n"
4945: "  const int N_b    = %d;                           // The number of basis functions\n"
4946: "  const int N_comp = %d;                           // The number of basis function components\n"
4947: "  const int N_bt   = N_b*N_comp;                    // The total number of scalar basis functions\n"
4948: "  const int N_q    = %d;                           // The number of quadrature points\n"
4949: "  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"
4950: "  const int N_t    = N_bst*N_bl;                    // The number of threads, N_bst * N_bl\n"
4951: "  const int N_bc   = N_t/N_comp;                    // The number of cells per batch (N_b*N_q*N_bl)\n"
4952: "  const int N_sbc  = N_bst / (N_q * N_comp);\n"
4953: "  const int N_sqc  = N_bst / N_bt;\n"
4954: "  /*const int N_c    = N_cb * N_bc;*/\n"
4955: "\n"
4956: "  /* Calculated indices */\n"
4957: "  /*const int tidx    = get_local_id(0) + get_local_size(0)*get_local_id(1);*/\n"
4958: "  const int tidx    = get_local_id(0);\n"
4959: "  const int blidx   = tidx / N_bst;                  // Block number for this thread\n"
4960: "  const int bidx    = tidx %% N_bt;                   // Basis function mapped to this thread\n"
4961: "  const int cidx    = tidx %% N_comp;                 // Basis component mapped to this thread\n"
4962: "  const int qidx    = tidx %% N_q;                    // Quadrature point mapped to this thread\n"
4963: "  const int blbidx  = tidx %% N_q + blidx*N_q;        // Cell mapped to this thread in the basis phase\n"
4964: "  const int blqidx  = tidx %% N_b + blidx*N_b;        // Cell mapped to this thread in the quadrature phase\n"
4965: "  const int gidx    = get_group_id(1)*get_num_groups(0) + get_group_id(0);\n"
4966: "  const int Goffset = gidx*N_cb*N_bc;\n",
4967:                             &count, dim, N_bl, N_b, N_c, N_q);STRING_ERROR_CHECK("Message to short");
4968:   /* Local memory */
4969:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4970: "\n"
4971: "  /* Quadrature data */\n"
4972: "  %s                w;                   // $w_q$, Quadrature weight at $x_q$\n"
4973: "  __local %s         phi_i[%d];    //[N_bt*N_q];  // $\\phi_i(x_q)$, Value of the basis function $i$ at $x_q$\n"
4974: "  __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"
4975: "  /* Geometric data */\n"
4976: "  __local %s        detJ[%d]; //[N_t];           // $|J(x_q)|$, Jacobian determinant at $x_q$\n"
4977: "  __local %s        invJ[%d];//[N_t*dim*dim];   // $J^{-1}(x_q)$, Jacobian inverse at $x_q$\n",
4978:                             &count, numeric_str, numeric_str, N_b*N_c*N_q, numeric_str, dim, N_b*N_c*N_q, numeric_str, N_t,
4979:                             numeric_str, N_t*dim*dim, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4980:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4981: "  /* FEM data */\n"
4982: "  __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",
4983:                             &count, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4984:   if (useAux) {
4985:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4986: "  __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",
4987:                             &count, numeric_str, N_t);STRING_ERROR_CHECK("Message to short");
4988:   }
4989:   if (useF0) {
4990:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4991: "  /* Intermediate calculations */\n"
4992: "  __local %s         f_0[%d]; //[N_t*N_sqc];      // $f_0(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4993:                               &count, numeric_str, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4994:   }
4995:   if (useF1) {
4996:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4997: "  __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",
4998:                               &count, numeric_str, dim, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4999:   }
5000:   /* TODO: If using elasticity, put in mu/lambda coefficients */
5001:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5002: "  /* Output data */\n"
5003: "  %s                e_i;                 // Coefficient $e_i$ of the residual\n\n",
5004:                             &count, numeric_str);STRING_ERROR_CHECK("Message to short");
5005:   /* One-time loads */
5006:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5007: "  /* These should be generated inline */\n"
5008: "  /* Load quadrature weights */\n"
5009: "  w = weights[qidx];\n"
5010: "  /* Load basis tabulation \\phi_i for this cell */\n"
5011: "  if (tidx < N_bt*N_q) {\n"
5012: "    phi_i[tidx]    = Basis[tidx];\n"
5013: "    phiDer_i[tidx] = BasisDerivatives[tidx];\n"
5014: "  }\n\n",
5015:                        &count);STRING_ERROR_CHECK("Message to short");
5016:   /* Batch loads */
5017:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5018: "  for (int batch = 0; batch < N_cb; ++batch) {\n"
5019: "    /* Load geometry */\n"
5020: "    detJ[tidx] = jacobianDeterminants[Goffset+batch*N_bc+tidx];\n"
5021: "    for (int n = 0; n < dim*dim; ++n) {\n"
5022: "      const int offset = n*N_t;\n"
5023: "      invJ[offset+tidx] = jacobianInverses[(Goffset+batch*N_bc)*dim*dim+offset+tidx];\n"
5024: "    }\n"
5025: "    /* Load coefficients u_i for this cell */\n"
5026: "    for (int n = 0; n < N_bt; ++n) {\n"
5027: "      const int offset = n*N_t;\n"
5028: "      u_i[offset+tidx] = coefficients[(Goffset*N_bt)+batch*N_t*N_b+offset+tidx];\n"
5029: "    }\n",
5030:                        &count);STRING_ERROR_CHECK("Message to short");
5031:   if (useAux) {
5032:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5033: "    /* Load coefficients a_i for this cell */\n"
5034: "    /* TODO: This should not be N_t here, it should be N_bc*N_comp_aux */\n"
5035: "    a_i[tidx] = coefficientsAux[Goffset+batch*N_t+tidx];\n",
5036:                             &count);STRING_ERROR_CHECK("Message to short");
5037:   }
5038:   /* Quadrature phase */
5039:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5040: "    barrier(CLK_LOCAL_MEM_FENCE);\n"
5041: "\n"
5042: "    /* Map coefficients to values at quadrature points */\n"
5043: "    for (int c = 0; c < N_sqc; ++c) {\n"
5044: "      const int cell          = c*N_bl*N_b + blqidx;\n"
5045: "      const int fidx          = (cell*N_q + qidx)*N_comp + cidx;\n",
5046:                        &count);STRING_ERROR_CHECK("Message to short");
5047:   if (useField) {
5048:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5049: "      %s  u[%d]; //[N_comp];     // $u(x_q)$, Value of the field at $x_q$\n",
5050:                               &count, numeric_str, N_c);STRING_ERROR_CHECK("Message to short");
5051:   }
5052:   if (useFieldDer) {
5053:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5054: "      %s%d   gradU[%d]; //[N_comp]; // $\\nabla u(x_q)$, Value of the field gradient at $x_q$\n",
5055:                               &count, numeric_str, dim, N_c);STRING_ERROR_CHECK("Message to short");
5056:   }
5057:   if (useFieldAux) {
5058:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5059: "      %s  a[%d]; //[1];     // $a(x_q)$, Value of the auxiliary fields at $x_q$\n",
5060:                               &count, numeric_str, 1);STRING_ERROR_CHECK("Message to short");
5061:   }
5062:   if (useFieldDerAux) {
5063:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5064: "      %s%d   gradA[%d]; //[1]; // $\\nabla a(x_q)$, Value of the auxiliary field gradient at $x_q$\n",
5065:                               &count, numeric_str, dim, 1);STRING_ERROR_CHECK("Message to short");
5066:   }
5067:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5068: "\n"
5069: "      for (int comp = 0; comp < N_comp; ++comp) {\n",
5070:                             &count);STRING_ERROR_CHECK("Message to short");
5071:   if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        u[comp] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
5072:   if (useFieldDer) {
5073:     switch (dim) {
5074:     case 1:
5075:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        gradU[comp].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
5076:     case 2:
5077:       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;
5078:     case 3:
5079:       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;
5080:     }
5081:   }
5082:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5083: "      }\n",
5084:                             &count);STRING_ERROR_CHECK("Message to short");
5085:   if (useFieldAux) {
5086:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      a[0] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");
5087:   }
5088:   if (useFieldDerAux) {
5089:     switch (dim) {
5090:     case 1:
5091:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      gradA[0].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
5092:     case 2:
5093:       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;
5094:     case 3:
5095:       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;
5096:     }
5097:   }
5098:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5099: "      /* Get field and derivatives at this quadrature point */\n"
5100: "      for (int i = 0; i < N_b; ++i) {\n"
5101: "        for (int comp = 0; comp < N_comp; ++comp) {\n"
5102: "          const int b    = i*N_comp+comp;\n"
5103: "          const int pidx = qidx*N_bt + b;\n"
5104: "          const int uidx = cell*N_bt + b;\n"
5105: "          %s%d   realSpaceDer;\n\n",
5106:                             &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
5107:   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");}
5108:   if (useFieldDer) {
5109:     switch (dim) {
5110:     case 2:
5111:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5112: "          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"
5113: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
5114: "          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"
5115: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n",
5116:                            &count);STRING_ERROR_CHECK("Message to short");break;
5117:     case 3:
5118:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5119: "          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"
5120: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
5121: "          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"
5122: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n"
5123: "          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"
5124: "          gradU[comp].z += u_i[uidx]*realSpaceDer.z;\n",
5125:                            &count);STRING_ERROR_CHECK("Message to short");break;
5126:     }
5127:   }
5128:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5129: "        }\n"
5130: "      }\n",
5131:                             &count);STRING_ERROR_CHECK("Message to short");
5132:   if (useFieldAux) {
5133:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"          a[0] += a_i[cell];\n", &count);STRING_ERROR_CHECK("Message to short");
5134:   }
5135:   /* Calculate residual at quadrature points: Should be generated by an weak form egine */
5136:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5137: "      /* Process values at quadrature points */\n",
5138:                             &count);STRING_ERROR_CHECK("Message to short");
5139:   switch (op) {
5140:   case LAPLACIAN:
5141:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
5142:     if (useF1) {
5143:       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");}
5144:       else        {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_1[fidx] = gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
5145:     }
5146:     break;
5147:   case ELASTICITY:
5148:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
5149:     if (useF1) {
5150:     switch (dim) {
5151:     case 2:
5152:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5153: "      switch (cidx) {\n"
5154: "      case 0:\n"
5155: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].x + gradU[0].x);\n"
5156: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].y + gradU[1].x);\n"
5157: "        break;\n"
5158: "      case 1:\n"
5159: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].x + gradU[0].y);\n"
5160: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].y + gradU[1].y);\n"
5161: "      }\n",
5162:                            &count);STRING_ERROR_CHECK("Message to short");break;
5163:     case 3:
5164:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5165: "      switch (cidx) {\n"
5166: "      case 0:\n"
5167: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].x + gradU[0].x);\n"
5168: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].y + gradU[1].x);\n"
5169: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].z + gradU[2].x);\n"
5170: "        break;\n"
5171: "      case 1:\n"
5172: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].x + gradU[0].y);\n"
5173: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[1].y);\n"
5174: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[2].y);\n"
5175: "        break;\n"
5176: "      case 2:\n"
5177: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].x + gradU[0].z);\n"
5178: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[1].z);\n"
5179: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[2].z);\n"
5180: "      }\n",
5181:                            &count);STRING_ERROR_CHECK("Message to short");break;
5182:     }}
5183:     break;
5184:   default:
5185:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PDE operator %d is not supported", op);
5186:   }
5187:   if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_0[fidx] *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");}
5188:   if (useF1) {
5189:     switch (dim) {
5190:     case 1:
5191:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_1[fidx].x *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
5192:     case 2:
5193:       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;
5194:     case 3:
5195:       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;
5196:     }
5197:   }
5198:   /* Thread transpose */
5199:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5200: "    }\n\n"
5201: "    /* ==== TRANSPOSE THREADS ==== */\n"
5202: "    barrier(CLK_LOCAL_MEM_FENCE);\n\n",
5203:                        &count);STRING_ERROR_CHECK("Message to short");
5204:   /* Basis phase */
5205:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5206: "    /* Map values at quadrature points to coefficients */\n"
5207: "    for (int c = 0; c < N_sbc; ++c) {\n"
5208: "      const int cell = c*N_bl*N_q + blbidx; /* Cell number in batch */\n"
5209: "\n"
5210: "      e_i = 0.0;\n"
5211: "      for (int q = 0; q < N_q; ++q) {\n"
5212: "        const int pidx = q*N_bt + bidx;\n"
5213: "        const int fidx = (cell*N_q + q)*N_comp + cidx;\n"
5214: "        %s%d   realSpaceDer;\n\n",
5215:                        &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");

5217:   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");}
5218:   if (useF1) {
5219:     switch (dim) {
5220:     case 2:
5221:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5222: "        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"
5223: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
5224: "        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"
5225: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n",
5226:                            &count);STRING_ERROR_CHECK("Message to short");break;
5227:     case 3:
5228:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5229: "        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"
5230: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
5231: "        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"
5232: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n"
5233: "        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"
5234: "        e_i           += realSpaceDer.z*f_1[fidx].z;\n",
5235:                            &count);STRING_ERROR_CHECK("Message to short");break;
5236:     }
5237:   }
5238:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
5239: "      }\n"
5240: "      /* Write element vector for N_{cbc} cells at a time */\n"
5241: "      elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
5242: "    }\n"
5243: "    /* ==== Could do one write per batch ==== */\n"
5244: "  }\n"
5245: "  return;\n"
5246: "}\n",
5247:                        &count);STRING_ERROR_CHECK("Message to short");
5248:   return(0);
5249: }

5253: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
5254: {
5255:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
5256:   PetscInt        dim, N_bl;
5257:   PetscBool       flg;
5258:   char           *buffer;
5259:   size_t          len;
5260:   char            errMsg[8192];
5261:   cl_int          ierr2;
5262:   PetscErrorCode  ierr;

5265:   PetscFEGetSpatialDimension(fem, &dim);
5266:   PetscMalloc1(8192, &buffer);
5267:   PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
5268:   PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
5269:   PetscOptionsHasName(((PetscObject)fem)->options,((PetscObject)fem)->prefix, "-petscfe_opencl_kernel_print", &flg);
5270:   if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
5271:   len  = strlen(buffer);
5272:   *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
5273:   clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
5274:   if (ierr != CL_SUCCESS) {
5275:     clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
5276:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
5277:   }
5278:   PetscFree(buffer);
5279:   *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
5280:   return(0);
5281: }

5285: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
5286: {
5287:   const PetscInt Nblocks = N/blockSize;

5290:   if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
5291:   *z = 1;
5292:   for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
5293:     *y = Nblocks / *x;
5294:     if (*x * *y == Nblocks) break;
5295:   }
5296:   if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
5297:   return(0);
5298: }

5302: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
5303: {
5304:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
5305:   PetscStageLog     stageLog;
5306:   PetscEventPerfLog eventLog = NULL;
5307:   PetscInt          stage;
5308:   PetscErrorCode    ierr;

5311:   PetscLogGetStageLog(&stageLog);
5312:   PetscStageLogGetCurrent(stageLog, &stage);
5313:   PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
5314:     /* Log performance info */
5315:   eventLog->eventInfo[ocl->residualEvent].count++;
5316:   eventLog->eventInfo[ocl->residualEvent].time  += time;
5317:   eventLog->eventInfo[ocl->residualEvent].flops += flops;
5318:   return(0);
5319: }

5323: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5324:                                                const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
5325: {
5326:   /* Nbc = batchSize */
5327:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
5328:   PetscPointFunc    f0_func;
5329:   PetscPointFunc    f1_func;
5330:   PetscQuadrature   q;
5331:   PetscInt          dim;
5332:   PetscInt          N_b;    /* The number of basis functions */
5333:   PetscInt          N_comp; /* The number of basis function components */
5334:   PetscInt          N_bt;   /* The total number of scalar basis functions */
5335:   PetscInt          N_q;    /* The number of quadrature points */
5336:   PetscInt          N_bst;  /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
5337:   PetscInt          N_t;    /* The number of threads, N_bst * N_bl */
5338:   PetscInt          N_bl;   /* The number of blocks */
5339:   PetscInt          N_bc;   /* The batch size, N_bl*N_q*N_b */
5340:   PetscInt          N_cb;   /* The number of batches */
5341:   PetscInt          numFlops, f0Flops = 0, f1Flops = 0;
5342:   PetscBool         useAux      = probAux ? PETSC_TRUE : PETSC_FALSE;
5343:   PetscBool         useField    = PETSC_FALSE;
5344:   PetscBool         useFieldDer = PETSC_TRUE;
5345:   PetscBool         useF0       = PETSC_TRUE;
5346:   PetscBool         useF1       = PETSC_TRUE;
5347:   /* OpenCL variables */
5348:   cl_program        ocl_prog;
5349:   cl_kernel         ocl_kernel;
5350:   cl_event          ocl_ev;         /* The event for tracking kernel execution */
5351:   cl_ulong          ns_start;       /* Nanoseconds counter on GPU at kernel start */
5352:   cl_ulong          ns_end;         /* Nanoseconds counter on GPU at kernel stop */
5353:   cl_mem            o_jacobianInverses, o_jacobianDeterminants;
5354:   cl_mem            o_coefficients, o_coefficientsAux, o_elemVec;
5355:   float            *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
5356:   double           *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
5357:   PetscReal        *r_invJ = NULL, *r_detJ = NULL;
5358:   void             *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
5359:   size_t            local_work_size[3], global_work_size[3];
5360:   size_t            realSize, x, y, z;
5361:   PetscErrorCode    ierr;

5364:   if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
5365:   PetscFEGetSpatialDimension(fem, &dim);
5366:   PetscFEGetQuadrature(fem, &q);
5367:   PetscFEGetDimension(fem, &N_b);
5368:   PetscFEGetNumComponents(fem, &N_comp);
5369:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
5370:   PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
5371:   N_bt  = N_b*N_comp;
5372:   N_q   = q->numPoints;
5373:   N_bst = N_bt*N_q;
5374:   N_t   = N_bst*N_bl;
5375:   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);
5376:   /* Calculate layout */
5377:   if (Ne % (N_cb*N_bc)) { /* Remainder cells */
5378:     PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);
5379:     return(0);
5380:   }
5381:   PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
5382:   local_work_size[0]  = N_bc*N_comp;
5383:   local_work_size[1]  = 1;
5384:   local_work_size[2]  = 1;
5385:   global_work_size[0] = x * local_work_size[0];
5386:   global_work_size[1] = y * local_work_size[1];
5387:   global_work_size[2] = z * local_work_size[2];
5388:   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);
5389:   PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
5390:   /* Generate code */
5391:   if (probAux) {
5392:     PetscSpace P;
5393:     PetscInt   NfAux, order, f;

5395:     PetscDSGetNumFields(probAux, &NfAux);
5396:     for (f = 0; f < NfAux; ++f) {
5397:       PetscFE feAux;

5399:       PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
5400:       PetscFEGetBasisSpace(feAux, &P);
5401:       PetscSpaceGetOrder(P, &order);
5402:       if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
5403:     }
5404:   }
5405:   PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
5406:   /* Create buffers on the device and send data over */
5407:   PetscDataTypeGetSize(ocl->realType, &realSize);
5408:   if (sizeof(PetscReal) != realSize) {
5409:     switch (ocl->realType) {
5410:     case PETSC_FLOAT:
5411:     {
5412:       PetscInt c, b, d;

5414:       PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
5415:       for (c = 0; c < Ne; ++c) {
5416:         f_detJ[c] = (float) geom[c].detJ;
5417:         for (d = 0; d < dim*dim; ++d) {
5418:           f_invJ[c*dim*dim+d] = (float) geom[c].invJ[d];
5419:         }
5420:         for (b = 0; b < N_bt; ++b) {
5421:           f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
5422:         }
5423:       }
5424:       if (coefficientsAux) { /* Assume P0 */
5425:         for (c = 0; c < Ne; ++c) {
5426:           f_coeffAux[c] = (float) coefficientsAux[c];
5427:         }
5428:       }
5429:       oclCoeff      = (void *) f_coeff;
5430:       if (coefficientsAux) {
5431:         oclCoeffAux = (void *) f_coeffAux;
5432:       } else {
5433:         oclCoeffAux = NULL;
5434:       }
5435:       oclInvJ       = (void *) f_invJ;
5436:       oclDetJ       = (void *) f_detJ;
5437:     }
5438:     break;
5439:     case PETSC_DOUBLE:
5440:     {
5441:       PetscInt c, b, d;

5443:       PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
5444:       for (c = 0; c < Ne; ++c) {
5445:         d_detJ[c] = (double) geom[c].detJ;
5446:         for (d = 0; d < dim*dim; ++d) {
5447:           d_invJ[c*dim*dim+d] = (double) geom[c].invJ[d];
5448:         }
5449:         for (b = 0; b < N_bt; ++b) {
5450:           d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
5451:         }
5452:       }
5453:       if (coefficientsAux) { /* Assume P0 */
5454:         for (c = 0; c < Ne; ++c) {
5455:           d_coeffAux[c] = (double) coefficientsAux[c];
5456:         }
5457:       }
5458:       oclCoeff      = (void *) d_coeff;
5459:       if (coefficientsAux) {
5460:         oclCoeffAux = (void *) d_coeffAux;
5461:       } else {
5462:         oclCoeffAux = NULL;
5463:       }
5464:       oclInvJ       = (void *) d_invJ;
5465:       oclDetJ       = (void *) d_detJ;
5466:     }
5467:     break;
5468:     default:
5469:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
5470:     }
5471:   } else {
5472:     PetscInt c, d;

5474:     PetscMalloc2(Ne*dim*dim,&r_invJ,Ne,&r_detJ);
5475:     for (c = 0; c < Ne; ++c) {
5476:       r_detJ[c] = geom[c].detJ;
5477:       for (d = 0; d < dim*dim; ++d) {
5478:         r_invJ[c*dim*dim+d] = geom[c].invJ[d];
5479:       }
5480:     }
5481:     oclCoeff    = (void *) coefficients;
5482:     oclCoeffAux = (void *) coefficientsAux;
5483:     oclInvJ     = (void *) r_invJ;
5484:     oclDetJ     = (void *) r_detJ;
5485:   }
5486:   o_coefficients         = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt    * realSize, oclCoeff,    &ierr);
5487:   if (coefficientsAux) {
5488:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclCoeffAux, &ierr);
5489:   } else {
5490:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY,                        Ne         * realSize, oclCoeffAux, &ierr);
5491:   }
5492:   o_jacobianInverses     = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ,     &ierr);
5493:   o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclDetJ,     &ierr);
5494:   o_elemVec              = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY,                       Ne*N_bt    * realSize, NULL,        &ierr);
5495:   /* Kernel launch */
5496:   clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
5497:   clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
5498:   clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
5499:   clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
5500:   clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
5501:   clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
5502:   clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
5503:   /* Read data back from device */
5504:   if (sizeof(PetscReal) != realSize) {
5505:     switch (ocl->realType) {
5506:     case PETSC_FLOAT:
5507:     {
5508:       float   *elem;
5509:       PetscInt c, b;

5511:       PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
5512:       PetscMalloc1(Ne*N_bt, &elem);
5513:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
5514:       for (c = 0; c < Ne; ++c) {
5515:         for (b = 0; b < N_bt; ++b) {
5516:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
5517:         }
5518:       }
5519:       PetscFree(elem);
5520:     }
5521:     break;
5522:     case PETSC_DOUBLE:
5523:     {
5524:       double  *elem;
5525:       PetscInt c, b;

5527:       PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
5528:       PetscMalloc1(Ne*N_bt, &elem);
5529:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
5530:       for (c = 0; c < Ne; ++c) {
5531:         for (b = 0; b < N_bt; ++b) {
5532:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
5533:         }
5534:       }
5535:       PetscFree(elem);
5536:     }
5537:     break;
5538:     default:
5539:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
5540:     }
5541:   } else {
5542:     PetscFree2(r_invJ,r_detJ);
5543:     clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
5544:   }
5545:   /* Log performance */
5546:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
5547:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END,   sizeof(cl_ulong), &ns_end,   NULL);
5548:   f0Flops = 0;
5549:   switch (ocl->op) {
5550:   case LAPLACIAN:
5551:     f1Flops = useAux ? dim : 0;break;
5552:   case ELASTICITY:
5553:     f1Flops = 2*dim*dim;break;
5554:   }
5555:   numFlops = Ne*(
5556:     N_q*(
5557:       N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
5558:       /*+
5559:        N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
5560:       +
5561:       N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
5562:     +
5563:     N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
5564:   PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
5565:   /* Cleanup */
5566:   clReleaseMemObject(o_coefficients);
5567:   clReleaseMemObject(o_coefficientsAux);
5568:   clReleaseMemObject(o_jacobianInverses);
5569:   clReleaseMemObject(o_jacobianDeterminants);
5570:   clReleaseMemObject(o_elemVec);
5571:   clReleaseKernel(ocl_kernel);
5572:   clReleaseProgram(ocl_prog);
5573:   return(0);
5574: }

5578: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
5579: {
5581:   fem->ops->setfromoptions          = NULL;
5582:   fem->ops->setup                   = PetscFESetUp_Basic;
5583:   fem->ops->view                    = NULL;
5584:   fem->ops->destroy                 = PetscFEDestroy_OpenCL;
5585:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
5586:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
5587:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_OpenCL;
5588:   fem->ops->integratebdresidual     = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
5589:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
5590:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
5591:   return(0);
5592: }

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

5597:   Level: intermediate

5599: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5600: M*/

5604: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
5605: {
5606:   PetscFE_OpenCL *ocl;
5607:   cl_uint         num_platforms;
5608:   cl_platform_id  platform_ids[42];
5609:   cl_uint         num_devices;
5610:   cl_device_id    device_ids[42];
5611:   cl_int          ierr2;
5612:   PetscErrorCode  ierr;

5616:   PetscNewLog(fem,&ocl);
5617:   fem->data = ocl;

5619:   /* Init Platform */
5620:   clGetPlatformIDs(42, platform_ids, &num_platforms);
5621:   if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
5622:   ocl->pf_id = platform_ids[0];
5623:   /* Init Device */
5624:   clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
5625:   if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
5626:   ocl->dev_id = device_ids[0];
5627:   /* Create context with one command queue */
5628:   ocl->ctx_id   = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
5629:   ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
5630:   /* Types */
5631:   ocl->realType = PETSC_FLOAT;
5632:   /* Register events */
5633:   PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
5634:   /* Equation handling */
5635:   ocl->op = LAPLACIAN;

5637:   PetscFEInitialize_OpenCL(fem);
5638:   return(0);
5639: }

5643: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
5644: {
5645:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

5649:   ocl->realType = realType;
5650:   return(0);
5651: }

5655: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
5656: {
5657:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

5662:   *realType = ocl->realType;
5663:   return(0);
5664: }

5666: #endif /* PETSC_HAVE_OPENCL */

5670: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
5671: {
5672:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5673:   PetscErrorCode     ierr;

5676:   CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5677:   PetscFree(cmp->embedding);
5678:   PetscFree(cmp);
5679:   return(0);
5680: }

5684: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
5685: {
5686:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5687:   DM                 K;
5688:   PetscReal         *subpoint;
5689:   PetscBLASInt      *pivots;
5690:   PetscBLASInt       n, info;
5691:   PetscScalar       *work, *invVscalar;
5692:   PetscInt           dim, pdim, spdim, j, s;
5693:   PetscErrorCode     ierr;

5696:   /* Get affine mapping from reference cell to each subcell */
5697:   PetscDualSpaceGetDM(fem->dualSpace, &K);
5698:   DMGetDimension(K, &dim);
5699:   DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
5700:   CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5701:   /* Determine dof embedding into subelements */
5702:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5703:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
5704:   PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
5705:   DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
5706:   for (s = 0; s < cmp->numSubelements; ++s) {
5707:     PetscInt sd = 0;

5709:     for (j = 0; j < pdim; ++j) {
5710:       PetscBool       inside;
5711:       PetscQuadrature f;
5712:       PetscInt        d, e;

5714:       PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
5715:       /* Apply transform to first point, and check that point is inside subcell */
5716:       for (d = 0; d < dim; ++d) {
5717:         subpoint[d] = -1.0;
5718:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
5719:       }
5720:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5721:       if (inside) {cmp->embedding[s*spdim+sd++] = j;}
5722:     }
5723:     if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
5724:   }
5725:   DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
5726:   /* Construct the change of basis from prime basis to nodal basis for each subelement */
5727:   PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
5728:   PetscMalloc2(spdim,&pivots,spdim,&work);
5729: #if defined(PETSC_USE_COMPLEX)
5730:   PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);
5731: #else
5732:   invVscalar = fem->invV;
5733: #endif
5734:   for (s = 0; s < cmp->numSubelements; ++s) {
5735:     for (j = 0; j < spdim; ++j) {
5736:       PetscReal      *Bf;
5737:       PetscQuadrature f;
5738:       PetscInt        q, k;

5740:       PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
5741:       PetscMalloc1(f->numPoints*spdim,&Bf);
5742:       PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
5743:       for (k = 0; k < spdim; ++k) {
5744:         /* n_j \cdot \phi_k */
5745:         invVscalar[(s*spdim + j)*spdim+k] = 0.0;
5746:         for (q = 0; q < f->numPoints; ++q) {
5747:           invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
5748:         }
5749:       }
5750:       PetscFree(Bf);
5751:     }
5752:     n = spdim;
5753:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info));
5754:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info));
5755:   }
5756: #if defined(PETSC_USE_COMPLEX)
5757:   for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
5758:   PetscFree(invVscalar);
5759: #endif
5760:   PetscFree2(pivots,work);
5761:   return(0);
5762: }

5766: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
5767: {
5768:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5769:   DM                 dm;
5770:   PetscInt           pdim;  /* Dimension of FE space P */
5771:   PetscInt           spdim; /* Dimension of subelement FE space P */
5772:   PetscInt           dim;   /* Spatial dimension */
5773:   PetscInt           comp;  /* Field components */
5774:   PetscInt          *subpoints;
5775:   PetscReal         *tmpB, *tmpD, *tmpH, *subpoint;
5776:   PetscInt           p, s, d, e, j, k;
5777:   PetscErrorCode     ierr;

5780:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
5781:   DMGetDimension(dm, &dim);
5782:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
5783:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5784:   PetscFEGetNumComponents(fem, &comp);
5785:   /* Divide points into subelements */
5786:   DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
5787:   DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
5788:   for (p = 0; p < npoints; ++p) {
5789:     for (s = 0; s < cmp->numSubelements; ++s) {
5790:       PetscBool inside;

5792:       /* Apply transform, and check that point is inside cell */
5793:       for (d = 0; d < dim; ++d) {
5794:         subpoint[d] = -1.0;
5795:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
5796:       }
5797:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5798:       if (inside) {subpoints[p] = s; break;}
5799:     }
5800:     if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
5801:   }
5802:   DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
5803:   /* Evaluate the prime basis functions at all points */
5804:   if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5805:   if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5806:   if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5807:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
5808:   /* Translate to the nodal basis */
5809:   if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
5810:   if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
5811:   if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
5812:   for (p = 0; p < npoints; ++p) {
5813:     const PetscInt s = subpoints[p];

5815:     if (B) {
5816:       /* Multiply by V^{-1} (spdim x spdim) */
5817:       for (j = 0; j < spdim; ++j) {
5818:         const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
5819:         PetscInt       c;

5821:         B[i] = 0.0;
5822:         for (k = 0; k < spdim; ++k) {
5823:           B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
5824:         }
5825:         for (c = 1; c < comp; ++c) {
5826:           B[i+c] = B[i];
5827:         }
5828:       }
5829:     }
5830:     if (D) {
5831:       /* Multiply by V^{-1} (spdim x spdim) */
5832:       for (j = 0; j < spdim; ++j) {
5833:         for (d = 0; d < dim; ++d) {
5834:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
5835:           PetscInt       c;

5837:           D[i] = 0.0;
5838:           for (k = 0; k < spdim; ++k) {
5839:             D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
5840:           }
5841:           for (c = 1; c < comp; ++c) {
5842:             D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
5843:           }
5844:         }
5845:       }
5846:     }
5847:     if (H) {
5848:       /* Multiply by V^{-1} (pdim x pdim) */
5849:       for (j = 0; j < spdim; ++j) {
5850:         for (d = 0; d < dim*dim; ++d) {
5851:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
5852:           PetscInt       c;

5854:           H[i] = 0.0;
5855:           for (k = 0; k < spdim; ++k) {
5856:             H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
5857:           }
5858:           for (c = 1; c < comp; ++c) {
5859:             H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
5860:           }
5861:         }
5862:       }
5863:     }
5864:   }
5865:   DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
5866:   if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5867:   if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5868:   if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5869:   return(0);
5870: }

5874: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
5875: {
5877:   fem->ops->setfromoptions          = NULL;
5878:   fem->ops->setup                   = PetscFESetUp_Composite;
5879:   fem->ops->view                    = NULL;
5880:   fem->ops->destroy                 = PetscFEDestroy_Composite;
5881:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
5882:   fem->ops->gettabulation           = PetscFEGetTabulation_Composite;
5883:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
5884:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
5885:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
5886:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
5887:   return(0);
5888: }

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

5893:   Level: intermediate

5895: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5896: M*/

5900: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
5901: {
5902:   PetscFE_Composite *cmp;
5903:   PetscErrorCode     ierr;

5907:   PetscNewLog(fem, &cmp);
5908:   fem->data = cmp;

5910:   cmp->cellRefiner    = REFINER_NOOP;
5911:   cmp->numSubelements = -1;
5912:   cmp->v0             = NULL;
5913:   cmp->jac            = NULL;

5915:   PetscFEInitialize_Composite(fem);
5916:   return(0);
5917: }

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

5924:   Not collective

5926:   Input Parameter:
5927: . fem - The PetscFE object

5929:   Output Parameters:
5930: + blockSize - The number of elements in a block
5931: . numBlocks - The number of blocks in a batch
5932: . batchSize - The number of elements in a batch
5933: - numBatches - The number of batches in a chunk

5935:   Level: intermediate

5937: .seealso: PetscFECreate()
5938: @*/
5939: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
5940: {
5941:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;

5949:   return(0);
5950: }

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

5957:   Not collective

5959:   Input Parameter:
5960: . fe - The PetscFE

5962:   Output Parameter:
5963: . dim - The dimension

5965:   Level: intermediate

5967: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
5968: @*/
5969: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
5970: {

5976:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
5977:   return(0);
5978: }

5980: /*
5981: Purpose: Compute element vector for chunk of elements

5983: Input:
5984:   Sizes:
5985:      Ne:  number of elements
5986:      Nf:  number of fields
5987:      PetscFE
5988:        dim: spatial dimension
5989:        Nb:  number of basis functions
5990:        Nc:  number of field components
5991:        PetscQuadrature
5992:          Nq:  number of quadrature points

5994:   Geometry:
5995:      PetscFECellGeom[Ne] possibly *Nq
5996:        PetscReal v0s[dim]
5997:        PetscReal n[dim]
5998:        PetscReal jacobians[dim*dim]
5999:        PetscReal jacobianInverses[dim*dim]
6000:        PetscReal jacobianDeterminants
6001:   FEM:
6002:      PetscFE
6003:        PetscQuadrature
6004:          PetscReal   quadPoints[Nq*dim]
6005:          PetscReal   quadWeights[Nq]
6006:        PetscReal   basis[Nq*Nb*Nc]
6007:        PetscReal   basisDer[Nq*Nb*Nc*dim]
6008:      PetscScalar coefficients[Ne*Nb*Nc]
6009:      PetscScalar elemVec[Ne*Nb*Nc]

6011:   Problem:
6012:      PetscInt f: the active field
6013:      f0, f1

6015:   Work Space:
6016:      PetscFE
6017:        PetscScalar f0[Nq*dim];
6018:        PetscScalar f1[Nq*dim*dim];
6019:        PetscScalar u[Nc];
6020:        PetscScalar gradU[Nc*dim];
6021:        PetscReal   x[dim];
6022:        PetscScalar realSpaceDer[dim];

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

6026: Input:
6027:   Sizes:
6028:      N_cb: Number of serial cell batches

6030:   Geometry:
6031:      PetscReal v0s[Ne*dim]
6032:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
6033:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
6034:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
6035:   FEM:
6036:      static PetscReal   quadPoints[Nq*dim]
6037:      static PetscReal   quadWeights[Nq]
6038:      static PetscReal   basis[Nq*Nb*Nc]
6039:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
6040:      PetscScalar coefficients[Ne*Nb*Nc]
6041:      PetscScalar elemVec[Ne*Nb*Nc]

6043: ex62.c:
6044:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
6045:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
6046:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
6047:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

6049: ex52.c:
6050:   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)
6051:   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)

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

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

6060: ex52_integrateElementOpenCL.c:
6061: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
6062:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
6063:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

6073:   Not collective

6075:   Input Parameters:
6076: + fem          - The PetscFE object for the field being integrated
6077: . prob         - The PetscDS specifying the discretizations and continuum functions
6078: . field        - The field being integrated
6079: . Ne           - The number of elements in the chunk
6080: . geom         - The cell geometry for each cell in the chunk
6081: . coefficients - The array of FEM basis coefficients for the elements
6082: . probAux      - The PetscDS specifying the auxiliary discretizations
6083: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

6085:   Output Parameter
6086: . integral     - the integral for this field

6088:   Level: developer

6090: .seealso: PetscFEIntegrateResidual()
6091: @*/
6092: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
6093:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
6094: {

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

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

6109:   Not collective

6111:   Input Parameters:
6112: + fem          - The PetscFE object for the field being integrated
6113: . prob         - The PetscDS specifying the discretizations and continuum functions
6114: . field        - The field being integrated
6115: . Ne           - The number of elements in the chunk
6116: . geom         - The cell geometry for each cell in the chunk
6117: . coefficients - The array of FEM basis coefficients for the elements
6118: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
6119: . probAux      - The PetscDS specifying the auxiliary discretizations
6120: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
6121: - t            - The time

6123:   Output Parameter
6124: . elemVec      - the element residual vectors from each element

6126:   Note:
6127: $ Loop over batch of elements (e):
6128: $   Loop over quadrature points (q):
6129: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
6130: $     Call f_0 and f_1
6131: $   Loop over element vector entries (f,fc --> i):
6132: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

6134:   Level: developer

6136: .seealso: PetscFEIntegrateResidual()
6137: @*/
6138: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
6139:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
6140: {

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

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

6155:   Not collective

6157:   Input Parameters:
6158: + fem          - The PetscFE object for the field being integrated
6159: . prob         - The PetscDS specifying the discretizations and continuum functions
6160: . field        - The field being integrated
6161: . Ne           - The number of elements in the chunk
6162: . geom         - The cell geometry for each cell in the chunk
6163: . coefficients - The array of FEM basis coefficients for the elements
6164: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
6165: . probAux      - The PetscDS specifying the auxiliary discretizations
6166: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
6167: - t            - The time

6169:   Output Parameter
6170: . elemVec      - the element residual vectors from each element

6172:   Level: developer

6174: .seealso: PetscFEIntegrateResidual()
6175: @*/
6176: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
6177:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
6178: {

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

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

6192:   Not collective

6194:   Input Parameters:
6195: + fem          - The PetscFE object for the field being integrated
6196: . prob         - The PetscDS specifying the discretizations and continuum functions
6197: . jtype        - The type of matrix pointwise functions that should be used
6198: . fieldI       - The test field being integrated
6199: . fieldJ       - The basis field being integrated
6200: . Ne           - The number of elements in the chunk
6201: . geom         - The cell geometry for each cell in the chunk
6202: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
6203: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
6204: . probAux      - The PetscDS specifying the auxiliary discretizations
6205: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
6206: . t            - The time
6207: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

6209:   Output Parameter
6210: . elemMat      - the element matrices for the Jacobian from each element

6212:   Note:
6213: $ Loop over batch of elements (e):
6214: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
6215: $     Loop over quadrature points (q):
6216: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
6217: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
6218: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
6219: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
6220: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
6221: */
6222: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
6223:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
6224: {

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

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

6238:   Not collective

6240:   Input Parameters:
6241: + fem          = The PetscFE object for the field being integrated
6242: . prob         - The PetscDS specifying the discretizations and continuum functions
6243: . fieldI       - The test field being integrated
6244: . fieldJ       - The basis field being integrated
6245: . Ne           - The number of elements in the chunk
6246: . geom         - The cell geometry for each cell in the chunk
6247: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
6248: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
6249: . probAux      - The PetscDS specifying the auxiliary discretizations
6250: . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
6251: . t            - The time
6252: - u_tShift     - A multiplier for the dF/du_t term (as opposed to the dF/du term)

6254:   Output Parameter
6255: . elemMat              - the element matrices for the Jacobian from each element

6257:   Note:
6258: $ Loop over batch of elements (e):
6259: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
6260: $     Loop over quadrature points (q):
6261: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
6262: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
6263: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
6264: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
6265: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
6266: */
6267: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
6268:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
6269: {

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

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

6285:   Input Parameter:
6286: . fe - The initial PetscFE

6288:   Output Parameter:
6289: . feRef - The refined PetscFE

6291:   Level: developer

6293: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
6294: @*/
6295: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
6296: {
6297:   PetscSpace       P, Pref;
6298:   PetscDualSpace   Q, Qref;
6299:   DM               K, Kref;
6300:   PetscQuadrature  q, qref;
6301:   const PetscReal *v0, *jac;
6302:   PetscInt         numComp, numSubelements;
6303:   PetscErrorCode   ierr;

6306:   PetscFEGetBasisSpace(fe, &P);
6307:   PetscFEGetDualSpace(fe, &Q);
6308:   PetscFEGetQuadrature(fe, &q);
6309:   PetscDualSpaceGetDM(Q, &K);
6310:   /* Create space */
6311:   PetscObjectReference((PetscObject) P);
6312:   Pref = P;
6313:   /* Create dual space */
6314:   PetscDualSpaceDuplicate(Q, &Qref);
6315:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
6316:   PetscDualSpaceSetDM(Qref, Kref);
6317:   DMDestroy(&Kref);
6318:   PetscDualSpaceSetUp(Qref);
6319:   /* Create element */
6320:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
6321:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
6322:   PetscFESetBasisSpace(*feRef, Pref);
6323:   PetscFESetDualSpace(*feRef, Qref);
6324:   PetscFEGetNumComponents(fe,    &numComp);
6325:   PetscFESetNumComponents(*feRef, numComp);
6326:   PetscFESetUp(*feRef);
6327:   PetscSpaceDestroy(&Pref);
6328:   PetscDualSpaceDestroy(&Qref);
6329:   /* Create quadrature */
6330:   PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
6331:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
6332:   PetscFESetQuadrature(*feRef, qref);
6333:   PetscQuadratureDestroy(&qref);
6334:   return(0);
6335: }

6339: /*@
6340:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

6342:   Collective on DM

6344:   Input Parameters:
6345: + dm         - The underlying DM for the domain
6346: . dim        - The spatial dimension
6347: . numComp    - The number of components
6348: . isSimplex  - Flag for simplex reference cell, otherwise its a tensor product
6349: . prefix     - The options prefix, or NULL
6350: - qorder     - The quadrature order

6352:   Output Parameter:
6353: . fem - The PetscFE object

6355:   Level: beginner

6357: .keywords: PetscFE, finite element
6358: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
6359: @*/
6360: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
6361: {
6362:   PetscQuadrature q;
6363:   DM              K;
6364:   PetscSpace      P;
6365:   PetscDualSpace  Q;
6366:   PetscInt        order;
6367:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
6368:   PetscErrorCode  ierr;

6371:   /* Create space */
6372:   PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
6373:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
6374:   PetscSpaceSetFromOptions(P);
6375:   PetscSpacePolynomialSetNumVariables(P, dim);
6376:   PetscSpaceSetUp(P);
6377:   PetscSpaceGetOrder(P, &order);
6378:   PetscSpacePolynomialGetTensor(P, &tensor);
6379:   /* Create dual space */
6380:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
6381:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
6382:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
6383:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
6384:   PetscDualSpaceSetDM(Q, K);
6385:   DMDestroy(&K);
6386:   PetscDualSpaceSetOrder(Q, order);
6387:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
6388:   PetscDualSpaceSetFromOptions(Q);
6389:   PetscDualSpaceSetUp(Q);
6390:   /* Create element */
6391:   PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
6392:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
6393:   PetscFESetFromOptions(*fem);
6394:   PetscFESetBasisSpace(*fem, P);
6395:   PetscFESetDualSpace(*fem, Q);
6396:   PetscFESetNumComponents(*fem, numComp);
6397:   PetscFESetUp(*fem);
6398:   PetscSpaceDestroy(&P);
6399:   PetscDualSpaceDestroy(&Q);
6400:   /* Create quadrature (with specified order if given) */
6401:   if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
6402:   else           {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
6403:   PetscFESetQuadrature(*fem, q);
6404:   PetscQuadratureDestroy(&q);
6405:   return(0);
6406: }