Actual source code: dtfe.c

petsc-master 2014-12-28
Report Typos and Errors
  1: /* Basis Jet Tabulation

  3: We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We
  4: follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can
  5: be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis
  6: as a prime basis.

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

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

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

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

 16:   n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k
 17:        \delta_{ji} = \sum_k \alpha_{ki} V_{jk}
 18:                  I = V \alpha

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

 22:    \alpha = V^{-1}

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

 26: Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials
 27: (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can
 28: be implemented exactly as in FIAT using functionals $L_j$.

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

 32: We will have three objects:
 33:  - Space, P: this just need point evaluation I think
 34:  - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q
 35:  - FEM: This keeps {P, P', Q}
 36: */
 37: #include <petsc-private/petscfeimpl.h> /*I "petscfe.h" I*/
 38: #include <petsc-private/dtimpl.h>
 39: #include <petsc-private/dmpleximpl.h>
 40: #include <petscdmshell.h>
 41: #include <petscdmplex.h>
 42: #include <petscblaslapack.h>

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

 55: PetscClassId PETSCSPACE_CLASSID = 0;

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

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

 65:   Not Collective

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

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

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

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

 89:   Level: advanced

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

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

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

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

109:   Collective on PetscSpace

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

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

118:   Level: intermediate

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

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

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

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

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

152:   Not Collective

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

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

160:   Level: intermediate

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

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

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

184:   Collective on PetscSpace

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

190:   Level: developer

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

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

211: /*
212:   PetscSpaceViewFromOptions - Processes command line options to determine if/how a PetscSpace is to be viewed.

214:   Collective on PetscSpace

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

221:   Level: intermediate

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

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

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

253:   Collective on PetscSpace

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

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

261:   Level: developer

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

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

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

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

304:   Collective on PetscSpace

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

309:   Level: developer

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

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

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

328:   Collective on PetscSpace

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

333:   Level: developer

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

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

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

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

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

359:   Collective on MPI_Comm

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

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

367:   Level: beginner

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

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

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

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

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

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

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

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

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

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

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

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

443:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

576:   Level: developer

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

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

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

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

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

616:   Level: developer

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

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

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

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

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

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

783:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

980:   if (D || H) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Cannot calculate derivatives for a DG space");
981:   if (npoints != dg->quad->numPoints) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot evaluate DG space on %d points != %d size", npoints, dg->quad->numPoints);
982:   PetscMemzero(B, npoints*npoints * sizeof(PetscReal));
983:   for (p = 0; p < npoints; ++p) {
984:     for (d = 0; d < dim; ++d) {
985:       if (PetscAbsReal(points[p*dim+d] - dg->quad->points[p*dim+d]) > 1.0e-10) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot evaluate DG point (%d, %d) %g != %g", p, d, points[p*dim+d], dg->quad->points[p*dim+d]);
986:     }
987:     B[p*npoints+p] = 1.0;
988:   }
989:   return(0);
990: }

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

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

1009:   Level: intermediate

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

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

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

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

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


1037: PetscClassId PETSCDUALSPACE_CLASSID = 0;

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

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

1047:   Not Collective

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

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

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

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

1071:   Level: advanced

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

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

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

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

1091:   Collective on PetscDualSpace

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

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

1100:   Level: intermediate

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

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

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

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

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

1134:   Not Collective

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

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

1142:   Level: intermediate

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

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

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

1166:   Collective on PetscDualSpace

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

1172:   Level: developer

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

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

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

1196:   Collective on PetscDualSpace

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

1203:   Level: intermediate

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

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

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

1235:   Collective on PetscDualSpace

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

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

1243:   Level: developer

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

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

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

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

1286:   Collective on PetscDualSpace

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

1291:   Level: developer

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

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

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

1310:   Collective on PetscDualSpace

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

1315:   Level: developer

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

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

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

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

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

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

1348:   Collective on MPI_Comm

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

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

1356:   Level: beginner

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

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

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

1374:   s->order = 0;

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

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

1385:   Collective on PetscDualSpace

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

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

1393:   Level: beginner

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

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

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

1413:   Not collective

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

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

1421:   Level: intermediate

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

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

1439:   Not collective

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

1445:   Level: intermediate

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

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

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

1467:   Not collective

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

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

1475:   Level: intermediate

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

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

1493:   Not collective

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

1499:   Level: intermediate

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

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

1516:   Not collective

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

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

1525:   Level: intermediate

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

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

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

1548:   Not collective

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

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

1556:   Level: intermediate

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

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

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

1577:   Not collective

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

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

1585:   Level: intermediate

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

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

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

1606:   Collective on PetscDualSpace

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

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

1616:   Level: advanced

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

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

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

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

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

1646:   Level: developer

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

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

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

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

1687:   Input Parameters:
1688: + sp - the PetscDualSpace object
1689: - height - the height of the mesh point for which the subspace is desired

1691:   Output Parameters:
1692:   bdsp - the subspace: must be destroyed by the user

1694:   Level: advanced

1696: .seealso: PetscDualSpace
1697: @*/
1698: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
1699: {

1705:   *bdsp = NULL;
1706:   if (sp->ops->getheightsubspace) {
1707:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1708:   }
1709:   return(0);
1710: }

1714: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
1715: {
1716:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1717:   PetscReal           D   = 1.0;
1718:   PetscInt            n, i;
1719:   PetscErrorCode      ierr;

1722:   DMGetDimension(sp->dm, &n);
1723:   if (lag->simplex || !lag->continuous) {
1724:     for (i = 1; i <= n; ++i) {
1725:       D *= ((PetscReal) (order+i))/i;
1726:     }
1727:     *dim = (PetscInt) (D + 0.5);
1728:   } else {
1729:     *dim = 1;
1730:     for (i = 0; i < n; ++i) *dim *= (order+1);
1731:   }
1732:   return(0);
1733: }

1737: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1738: {
1739:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1740:   DM                  dm    = sp->dm;
1741:   PetscInt            order = sp->order;
1742:   PetscBool           disc  = lag->continuous ? PETSC_FALSE : PETSC_TRUE;
1743:   PetscSection        csection;
1744:   Vec                 coordinates;
1745:   PetscReal          *qpoints, *qweights;
1746:   PetscInt           *closure = NULL, closureSize, c;
1747:   PetscInt            depth, dim, pdimMax, pMax = 0, *pStart, *pEnd, cell, coneSize, d, n, f = 0;
1748:   PetscBool           simplex;
1749:   PetscErrorCode      ierr;

1752:   /* Classify element type */
1753:   DMGetDimension(dm, &dim);
1754:   DMPlexGetDepth(dm, &depth);
1755:   PetscCalloc1(dim+1, &lag->numDof);
1756:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
1757:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1758:   DMPlexGetConeSize(dm, pStart[depth], &coneSize);
1759:   DMGetCoordinateSection(dm, &csection);
1760:   DMGetCoordinatesLocal(dm, &coordinates);
1761:   if (depth == 1) {
1762:     if      (coneSize == dim+1)    simplex = PETSC_TRUE;
1763:     else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
1764:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1765:   }
1766:   else if (depth == dim) {
1767:     if      (coneSize == dim+1)   simplex = PETSC_TRUE;
1768:     else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
1769:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1770:   }
1771:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
1772:   lag->simplex = simplex;
1773:   PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
1774:   pdimMax *= (pEnd[dim] - pStart[dim]);
1775:   PetscMalloc1(pdimMax, &sp->functional);
1776:   for (d = 0; d <= depth; d++) {
1777:     pMax = PetscMax(pMax,pEnd[d]);
1778:   }
1779:   if (!dim) {
1780:     PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1781:     PetscMalloc1(1, &qpoints);
1782:     PetscMalloc1(1, &qweights);
1783:     PetscQuadratureSetOrder(sp->functional[f], 0);
1784:     PetscQuadratureSetData(sp->functional[f], PETSC_DETERMINE, 1, qpoints, qweights);
1785:     qpoints[0]  = 0.0;
1786:     qweights[0] = 1.0;
1787:     ++f;
1788:     lag->numDof[0] = 1;
1789:   } else {
1790:     PetscBT seen;

1792:     PetscBTCreate(pMax, &seen);
1793:     PetscBTMemzero(pMax, seen);
1794:     for (cell = pStart[depth]; cell < pEnd[depth]; ++cell) {
1795:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1796:       for (c = 0; c < closureSize*2; c += 2) {
1797:         const PetscInt p = closure[c];

1799:         if (PetscBTLookup(seen, p)) continue;
1800:         PetscBTSet(seen, p);
1801:         if ((p >= pStart[0]) && (p < pEnd[0])) {
1802:           /* Vertices */
1803:           const PetscScalar *coords;
1804:           PetscInt           dof, off, d;

1806:           if (order < 1 || disc) continue;
1807:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1808:           PetscMalloc1(dim, &qpoints);
1809:           PetscMalloc1(1, &qweights);
1810:           PetscQuadratureSetOrder(sp->functional[f], 0);
1811:           PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1812:           VecGetArrayRead(coordinates, &coords);
1813:           PetscSectionGetDof(csection, p, &dof);
1814:           PetscSectionGetOffset(csection, p, &off);
1815:           if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim);
1816:           for (d = 0; d < dof; ++d) {qpoints[d] = PetscRealPart(coords[off+d]);}
1817:           qweights[0] = 1.0;
1818:           ++f;
1819:           VecRestoreArrayRead(coordinates, &coords);
1820:           lag->numDof[0] = 1;
1821:         } else if ((p >= pStart[1]) && (p < pEnd[1])) {
1822:           /* Edges */
1823:           PetscScalar *coords;
1824:           PetscInt     num = order-1, k;

1826:           if (order < 2 || disc) continue;
1827:           coords = NULL;
1828:           DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);
1829:           if (n != dim*2) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d has %d coordinate values instead of %d", p, n, dim*2);
1830:           for (k = 1; k <= num; ++k) {
1831:             PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1832:             PetscMalloc1(dim, &qpoints);
1833:             PetscMalloc1(1, &qweights);
1834:             PetscQuadratureSetOrder(sp->functional[f], 0);
1835:             PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1836:             for (d = 0; d < dim; ++d) {qpoints[d] = k*PetscRealPart(coords[1*dim+d] - coords[0*dim+d])/order + PetscRealPart(coords[0*dim+d]);}
1837:             qweights[0] = 1.0;
1838:             ++f;
1839:           }
1840:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);
1841:           lag->numDof[1] = num;
1842:         } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) {
1843:           /* Faces */

1845:           if (disc) continue;
1846:           if ( simplex && (order < 3)) continue;
1847:           if (!simplex && (order < 2)) continue;
1848:           lag->numDof[depth-1] = 0;
1849:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces");
1850:         } else if ((p >= pStart[depth]) && (p < pEnd[depth])) {
1851:           /* Cells */
1852:           PetscInt     orderEff = lag->continuous && order ? (simplex ? order-3 : order-2) : order;
1853:           PetscReal    denom    = order ? (lag->continuous ? order : (simplex ? order+3 : order+2)) : (simplex ? 3 : 2);
1854:           PetscScalar *coords   = NULL;
1855:           PetscReal    dx = 2.0/denom, *v0, *J, *invJ, detJ;
1856:           PetscInt    *ind, *tup;
1857:           PetscInt     cdim, csize, d, d2, o;

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

1865:           PetscCalloc5(dim,&ind,dim,&tup,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1866:           DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
1867:           if (simplex || !lag->continuous) {
1868:             for (o = 0; o <= orderEff; ++o) {
1869:               PetscMemzero(ind, dim*sizeof(PetscInt));
1870:               while (ind[0] >= 0) {
1871:                 PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1872:                 PetscMalloc1(dim, &qpoints);
1873:                 PetscMalloc1(1,   &qweights);
1874:                 PetscQuadratureSetOrder(sp->functional[f], 0);
1875:                 PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1876:                 LatticePoint_Internal(dim, o, ind, tup);
1877:                 for (d = 0; d < dim; ++d) {
1878:                   qpoints[d] = v0[d];
1879:                   for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
1880:                 }
1881:                 qweights[0] = 1.0;
1882:                 ++f;
1883:               }
1884:             }
1885:           } else {
1886:             while (ind[0] >= 0) {
1887:               PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1888:               PetscMalloc1(dim, &qpoints);
1889:               PetscMalloc1(1,   &qweights);
1890:               PetscQuadratureSetOrder(sp->functional[f], 0);
1891:               PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1892:               TensorPoint_Internal(dim, orderEff+1, ind, tup);
1893:               for (d = 0; d < dim; ++d) {
1894:                 qpoints[d] = v0[d];
1895:                 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d]+1)*dx);
1896:               }
1897:               qweights[0] = 1.0;
1898:               ++f;
1899:             }
1900:           }
1901:           PetscFree5(ind,tup,v0,J,invJ);
1902:           DMPlexVecRestoreClosure(dm, csection, coordinates, p, &csize, &coords);
1903:           lag->numDof[depth] = cdim;
1904:         }
1905:       }
1906:       DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);
1907:     }
1908:     PetscBTDestroy(&seen);
1909:   }
1910:   if (pEnd[dim] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d not equal to dimension %d", f, pdimMax);
1911:   PetscFree2(pStart,pEnd);
1912:   if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
1913:   return(0);
1914: }

1918: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
1919: {
1920:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1921:   PetscErrorCode      ierr;

1924:   PetscFree(lag->numDof);
1925:   PetscFree(lag);
1926:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
1927:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
1928:   return(0);
1929: }

1933: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
1934: {
1935:   PetscInt       order;
1936:   PetscBool      cont;

1940:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
1941:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
1942:   PetscDualSpaceGetOrder(sp, &order);
1943:   PetscDualSpaceSetOrder(*spNew, order);
1944:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
1945:   PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
1946:   return(0);
1947: }

1951: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscDualSpace sp)
1952: {
1953:   PetscBool      continuous, flg;

1957:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
1958:   PetscOptionsHead("PetscDualSpace Lagrange Options");
1959:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
1960:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
1961:   PetscOptionsTail();
1962:   return(0);
1963: }

1967: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
1968: {
1969:   DM              K;
1970:   const PetscInt *numDof;
1971:   PetscInt        spatialDim, Nc, size = 0, d;
1972:   PetscErrorCode  ierr;

1975:   PetscDualSpaceGetDM(sp, &K);
1976:   PetscDualSpaceGetNumDof(sp, &numDof);
1977:   DMGetDimension(K, &spatialDim);
1978:   DMPlexGetHeightStratum(K, 0, NULL, &Nc);
1979:   if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
1980:   for (d = 0; d <= spatialDim; ++d) {
1981:     PetscInt pStart, pEnd;

1983:     DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
1984:     size += (pEnd-pStart)*numDof[d];
1985:   }
1986:   *dim = size;
1987:   return(0);
1988: }

1992: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
1993: {
1994:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

1997:   *numDof = lag->numDof;
1998:   return(0);
1999: }

2003: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2004: {
2005:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2010:   *continuous = lag->continuous;
2011:   return(0);
2012: }

2016: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2017: {
2018:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2022:   lag->continuous = continuous;
2023:   return(0);
2024: }

2028: /*@
2029:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2031:   Not Collective

2033:   Input Parameter:
2034: . sp         - the PetscDualSpace

2036:   Output Parameter:
2037: . continuous - flag for element continuity

2039:   Level: intermediate

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

2051:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2052:   return(0);
2053: }

2057: /*@
2058:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2060:   Logically Collective on PetscDualSpace

2062:   Input Parameters:
2063: + sp         - the PetscDualSpace
2064: - continuous - flag for element continuity

2066:   Options Database:
2067: . -petscdualspace_lagrange_continuity <bool>

2069:   Level: intermediate

2071: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2072: .seealso: PetscDualSpaceLagrangeGetContinuity()
2073: @*/
2074: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2075: {

2081:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2082:   return(0);
2083: }

2087: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
2088: {
2089:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2090:   PetscBool          continuous;
2091:   PetscInt           order;
2092:   PetscErrorCode     ierr;

2097:   PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
2098:   PetscDualSpaceGetOrder(sp,&order);
2099:   if (height == 0) {
2100:     PetscObjectReference((PetscObject)sp);
2101:     *bdsp = sp;
2102:   }
2103:   else if (continuous == PETSC_FALSE || !order) {
2104:     *bdsp = NULL;
2105:   }
2106:   else {
2107:     DM dm, K;
2108:     PetscInt dim;

2110:     PetscDualSpaceGetDM(sp,&dm);
2111:     DMGetDimension(dm,&dim);
2112:     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);}
2113:     PetscDualSpaceDuplicate(sp,bdsp);
2114:     PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplex, &K);
2115:     PetscDualSpaceSetDM(*bdsp, K);
2116:     DMDestroy(&K);
2117:     PetscDualSpaceSetUp(*bdsp);
2118:   }
2119:   return(0);
2120: }

2124: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2125: {
2127:   sp->ops->setfromoptions    = PetscDualSpaceSetFromOptions_Lagrange;
2128:   sp->ops->setup             = PetscDualSpaceSetUp_Lagrange;
2129:   sp->ops->view              = NULL;
2130:   sp->ops->destroy           = PetscDualSpaceDestroy_Lagrange;
2131:   sp->ops->duplicate         = PetscDualSpaceDuplicate_Lagrange;
2132:   sp->ops->getdimension      = PetscDualSpaceGetDimension_Lagrange;
2133:   sp->ops->getnumdof         = PetscDualSpaceGetNumDof_Lagrange;
2134:   sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
2135:   return(0);
2136: }

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

2141:   Level: intermediate

2143: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2144: M*/

2148: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2149: {
2150:   PetscDualSpace_Lag *lag;
2151:   PetscErrorCode      ierr;

2155:   PetscNewLog(sp,&lag);
2156:   sp->data = lag;

2158:   lag->numDof     = NULL;
2159:   lag->simplex    = PETSC_TRUE;
2160:   lag->continuous = PETSC_TRUE;

2162:   PetscDualSpaceInitialize_Lagrange(sp);
2163:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2164:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2165:   return(0);
2166: }

2170: PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
2171: {
2173:   return(0);
2174: }

2178: PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
2179: {
2180:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2181:   PetscErrorCode         ierr;

2184:   PetscFree(s);
2185:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
2186:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
2187:   return(0);
2188: }

2192: PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace *spNew)
2193: {
2194:   PetscInt       dim, d;

2198:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2199:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACESIMPLE);
2200:   PetscDualSpaceGetDimension(sp, &dim);
2201:   PetscDualSpaceSimpleSetDimension(*spNew, dim);
2202:   for (d = 0; d < dim; ++d) {
2203:     PetscQuadrature q;

2205:     PetscDualSpaceGetFunctional(sp, d, &q);
2206:     PetscDualSpaceSimpleSetFunctional(*spNew, d, q);
2207:   }
2208:   return(0);
2209: }

2213: PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscDualSpace sp)
2214: {
2216:   return(0);
2217: }

2221: PetscErrorCode PetscDualSpaceGetDimension_Simple(PetscDualSpace sp, PetscInt *dim)
2222: {
2223:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;

2226:   *dim = s->dim;
2227:   return(0);
2228: }

2232: PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
2233: {
2234:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2235:   PetscInt               f;
2236:   PetscErrorCode         ierr;

2239:   for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
2240:   PetscFree(sp->functional);
2241:   s->dim = dim;
2242:   PetscCalloc1(s->dim, &sp->functional);
2243:   return(0);
2244: }

2248: PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
2249: {
2250:   PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2251:   PetscErrorCode         ierr;

2254:   if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim);
2255:   PetscObjectReference((PetscObject) q);
2256:   sp->functional[f] = q;
2257:   return(0);
2258: }

2262: /*@
2263:   PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis

2265:   Logically Collective on PetscDualSpace

2267:   Input Parameters:
2268: + sp  - the PetscDualSpace
2269: - dim - the basis dimension

2271:   Level: intermediate

2273: .keywords: PetscDualSpace, dimension
2274: .seealso: PetscDualSpaceSimpleSetFunctional()
2275: @*/
2276: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
2277: {

2283:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
2284:   return(0);
2285: }

2289: /*@
2290:   PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space

2292:   Not Collective

2294:   Input Parameters:
2295: + sp  - the PetscDualSpace
2296: . f - the basis index
2297: - q - the basis functional

2299:   Level: intermediate

2301: .keywords: PetscDualSpace, functional
2302: .seealso: PetscDualSpaceSimpleSetDimension()
2303: @*/
2304: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
2305: {

2310:   PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
2311:   return(0);
2312: }

2316: PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
2317: {
2319:   sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple;
2320:   sp->ops->setup          = PetscDualSpaceSetUp_Simple;
2321:   sp->ops->view           = NULL;
2322:   sp->ops->destroy        = PetscDualSpaceDestroy_Simple;
2323:   sp->ops->duplicate      = PetscDualSpaceDuplicate_Simple;
2324:   sp->ops->getdimension   = PetscDualSpaceGetDimension_Simple;
2325:   sp->ops->getnumdof      = NULL;
2326:   return(0);
2327: }

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

2332:   Level: intermediate

2334: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2335: M*/

2339: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
2340: {
2341:   PetscDualSpace_Simple *s;
2342:   PetscErrorCode         ierr;

2346:   PetscNewLog(sp,&s);
2347:   sp->data = s;

2349:   s->dim = 0;

2351:   PetscDualSpaceInitialize_Simple(sp);
2352:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
2353:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
2354:   return(0);
2355: }


2358: PetscClassId PETSCFE_CLASSID = 0;

2360: PetscFunctionList PetscFEList              = NULL;
2361: PetscBool         PetscFERegisterAllCalled = PETSC_FALSE;

2365: /*@C
2366:   PetscFERegister - Adds a new PetscFE implementation

2368:   Not Collective

2370:   Input Parameters:
2371: + name        - The name of a new user-defined creation routine
2372: - create_func - The creation routine itself

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

2377:   Sample usage:
2378: .vb
2379:     PetscFERegister("my_fe", MyPetscFECreate);
2380: .ve

2382:   Then, your PetscFE type can be chosen with the procedural interface via
2383: .vb
2384:     PetscFECreate(MPI_Comm, PetscFE *);
2385:     PetscFESetType(PetscFE, "my_fe");
2386: .ve
2387:    or at runtime via the option
2388: .vb
2389:     -petscfe_type my_fe
2390: .ve

2392:   Level: advanced

2394: .keywords: PetscFE, register
2395: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()

2397: @*/
2398: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2399: {

2403:   PetscFunctionListAdd(&PetscFEList, sname, function);
2404:   return(0);
2405: }

2409: /*@C
2410:   PetscFESetType - Builds a particular PetscFE

2412:   Collective on PetscFE

2414:   Input Parameters:
2415: + fem  - The PetscFE object
2416: - name - The kind of FEM space

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

2421:   Level: intermediate

2423: .keywords: PetscFE, set, type
2424: .seealso: PetscFEGetType(), PetscFECreate()
2425: @*/
2426: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2427: {
2428:   PetscErrorCode (*r)(PetscFE);
2429:   PetscBool      match;

2434:   PetscObjectTypeCompare((PetscObject) fem, name, &match);
2435:   if (match) return(0);

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

2441:   if (fem->ops->destroy) {
2442:     (*fem->ops->destroy)(fem);
2443:     fem->ops->destroy = NULL;
2444:   }
2445:   (*r)(fem);
2446:   PetscObjectChangeTypeName((PetscObject) fem, name);
2447:   return(0);
2448: }

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

2455:   Not Collective

2457:   Input Parameter:
2458: . fem  - The PetscFE

2460:   Output Parameter:
2461: . name - The PetscFE type name

2463:   Level: intermediate

2465: .keywords: PetscFE, get, type, name
2466: .seealso: PetscFESetType(), PetscFECreate()
2467: @*/
2468: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2469: {

2475:   if (!PetscFERegisterAllCalled) {
2476:     PetscFERegisterAll();
2477:   }
2478:   *name = ((PetscObject) fem)->type_name;
2479:   return(0);
2480: }

2484: /*@C
2485:   PetscFEView - Views a PetscFE

2487:   Collective on PetscFE

2489:   Input Parameter:
2490: + fem - the PetscFE object to view
2491: - v   - the viewer

2493:   Level: developer

2495: .seealso PetscFEDestroy()
2496: @*/
2497: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
2498: {

2503:   if (!v) {
2504:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
2505:   }
2506:   if (fem->ops->view) {
2507:     (*fem->ops->view)(fem, v);
2508:   }
2509:   return(0);
2510: }

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

2517:   Collective on PetscFE

2519:   Input Parameters:
2520: + fem    - the PetscFE
2521: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
2522: - optionname - option to activate viewing

2524:   Level: intermediate

2526: .keywords: PetscFE, view, options, database
2527: .seealso: VecViewFromOptions(), MatViewFromOptions()
2528: */
2529: PetscErrorCode PetscFEViewFromOptions(PetscFE fem, const char prefix[], const char optionname[])
2530: {
2531:   PetscViewer       viewer;
2532:   PetscViewerFormat format;
2533:   PetscBool         flg;
2534:   PetscErrorCode    ierr;

2537:   if (prefix) {
2538:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), prefix, optionname, &viewer, &format, &flg);
2539:   } else {
2540:     PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), ((PetscObject) fem)->prefix, optionname, &viewer, &format, &flg);
2541:   }
2542:   if (flg) {
2543:     PetscViewerPushFormat(viewer, format);
2544:     PetscFEView(fem, viewer);
2545:     PetscViewerPopFormat(viewer);
2546:     PetscViewerDestroy(&viewer);
2547:   }
2548:   return(0);
2549: }

2553: /*@
2554:   PetscFESetFromOptions - sets parameters in a PetscFE from the options database

2556:   Collective on PetscFE

2558:   Input Parameter:
2559: . fem - the PetscFE object to set options for

2561:   Options Database:
2562: . -petscfe_num_blocks  the number of cell blocks to integrate concurrently
2563: . -petscfe_num_batches the number of cell batches to integrate serially

2565:   Level: developer

2567: .seealso PetscFEView()
2568: @*/
2569: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
2570: {
2571:   const char    *defaultType;
2572:   char           name[256];
2573:   PetscBool      flg;

2578:   if (!((PetscObject) fem)->type_name) {
2579:     defaultType = PETSCFEBASIC;
2580:   } else {
2581:     defaultType = ((PetscObject) fem)->type_name;
2582:   }
2583:   if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}

2585:   PetscObjectOptionsBegin((PetscObject) fem);
2586:   PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
2587:   if (flg) {
2588:     PetscFESetType(fem, name);
2589:   } else if (!((PetscObject) fem)->type_name) {
2590:     PetscFESetType(fem, defaultType);
2591:   }
2592:   PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
2593:   PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
2594:   if (fem->ops->setfromoptions) {
2595:     (*fem->ops->setfromoptions)(fem);
2596:   }
2597:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
2598:   PetscObjectProcessOptionsHandlers((PetscObject) fem);
2599:   PetscOptionsEnd();
2600:   PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
2601:   return(0);
2602: }

2606: /*@C
2607:   PetscFESetUp - Construct data structures for the PetscFE

2609:   Collective on PetscFE

2611:   Input Parameter:
2612: . fem - the PetscFE object to setup

2614:   Level: developer

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

2624:   if (fem->ops->setup) {(*fem->ops->setup)(fem);}
2625:   return(0);
2626: }

2630: /*@
2631:   PetscFEDestroy - Destroys a PetscFE object

2633:   Collective on PetscFE

2635:   Input Parameter:
2636: . fem - the PetscFE object to destroy

2638:   Level: developer

2640: .seealso PetscFEView()
2641: @*/
2642: PetscErrorCode PetscFEDestroy(PetscFE *fem)
2643: {

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

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

2653:   PetscFree((*fem)->numDof);
2654:   PetscFree((*fem)->invV);
2655:   PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
2656:   PetscSpaceDestroy(&(*fem)->basisSpace);
2657:   PetscDualSpaceDestroy(&(*fem)->dualSpace);
2658:   PetscQuadratureDestroy(&(*fem)->quadrature);

2660:   if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
2661:   PetscHeaderDestroy(fem);
2662:   return(0);
2663: }

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

2670:   Collective on MPI_Comm

2672:   Input Parameter:
2673: . comm - The communicator for the PetscFE object

2675:   Output Parameter:
2676: . fem - The PetscFE object

2678:   Level: beginner

2680: .seealso: PetscFESetType(), PETSCFEGALERKIN
2681: @*/
2682: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
2683: {
2684:   PetscFE        f;

2689:   PetscCitationsRegister(FECitation,&FEcite);
2690:   *fem = NULL;
2691:   PetscFEInitializePackage();

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

2696:   f->basisSpace    = NULL;
2697:   f->dualSpace     = NULL;
2698:   f->numComponents = 1;
2699:   f->numDof        = NULL;
2700:   f->invV          = NULL;
2701:   f->B             = NULL;
2702:   f->D             = NULL;
2703:   f->H             = NULL;
2704:   PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
2705:   f->blockSize     = 0;
2706:   f->numBlocks     = 1;
2707:   f->batchSize     = 0;
2708:   f->numBatches    = 1;

2710:   *fem = f;
2711:   return(0);
2712: }

2716: /*@
2717:   PetscFEGetSpatialDimension - Returns the spatial dimension of the element

2719:   Not collective

2721:   Input Parameter:
2722: . fem - The PetscFE object

2724:   Output Parameter:
2725: . dim - The spatial dimension

2727:   Level: intermediate

2729: .seealso: PetscFECreate()
2730: @*/
2731: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
2732: {
2733:   DM             dm;

2739:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
2740:   DMGetDimension(dm, dim);
2741:   return(0);
2742: }

2746: /*@
2747:   PetscFESetNumComponents - Sets the number of components in the element

2749:   Not collective

2751:   Input Parameters:
2752: + fem - The PetscFE object
2753: - comp - The number of field components

2755:   Level: intermediate

2757: .seealso: PetscFECreate()
2758: @*/
2759: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
2760: {
2763:   fem->numComponents = comp;
2764:   return(0);
2765: }

2769: /*@
2770:   PetscFEGetNumComponents - Returns the number of components in the element

2772:   Not collective

2774:   Input Parameter:
2775: . fem - The PetscFE object

2777:   Output Parameter:
2778: . comp - The number of field components

2780:   Level: intermediate

2782: .seealso: PetscFECreate()
2783: @*/
2784: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
2785: {
2789:   *comp = fem->numComponents;
2790:   return(0);
2791: }

2795: /*@
2796:   PetscFESetTileSizes - Sets the tile sizes for evaluation

2798:   Not collective

2800:   Input Parameters:
2801: + fem - The PetscFE object
2802: . blockSize - The number of elements in a block
2803: . numBlocks - The number of blocks in a batch
2804: . batchSize - The number of elements in a batch
2805: - numBatches - The number of batches in a chunk

2807:   Level: intermediate

2809: .seealso: PetscFECreate()
2810: @*/
2811: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
2812: {
2815:   fem->blockSize  = blockSize;
2816:   fem->numBlocks  = numBlocks;
2817:   fem->batchSize  = batchSize;
2818:   fem->numBatches = numBatches;
2819:   return(0);
2820: }

2824: /*@
2825:   PetscFEGetTileSizes - Returns the tile sizes for evaluation

2827:   Not collective

2829:   Input Parameter:
2830: . fem - The PetscFE object

2832:   Output Parameters:
2833: + blockSize - The number of elements in a block
2834: . numBlocks - The number of blocks in a batch
2835: . batchSize - The number of elements in a batch
2836: - numBatches - The number of batches in a chunk

2838:   Level: intermediate

2840: .seealso: PetscFECreate()
2841: @*/
2842: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
2843: {
2850:   if (blockSize)  *blockSize  = fem->blockSize;
2851:   if (numBlocks)  *numBlocks  = fem->numBlocks;
2852:   if (batchSize)  *batchSize  = fem->batchSize;
2853:   if (numBatches) *numBatches = fem->numBatches;
2854:   return(0);
2855: }

2859: /*@
2860:   PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution

2862:   Not collective

2864:   Input Parameter:
2865: . fem - The PetscFE object

2867:   Output Parameter:
2868: . sp - The PetscSpace object

2870:   Level: intermediate

2872: .seealso: PetscFECreate()
2873: @*/
2874: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
2875: {
2879:   *sp = fem->basisSpace;
2880:   return(0);
2881: }

2885: /*@
2886:   PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution

2888:   Not collective

2890:   Input Parameters:
2891: + fem - The PetscFE object
2892: - sp - The PetscSpace object

2894:   Level: intermediate

2896: .seealso: PetscFECreate()
2897: @*/
2898: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
2899: {

2905:   PetscSpaceDestroy(&fem->basisSpace);
2906:   fem->basisSpace = sp;
2907:   PetscObjectReference((PetscObject) fem->basisSpace);
2908:   return(0);
2909: }

2913: /*@
2914:   PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product

2916:   Not collective

2918:   Input Parameter:
2919: . fem - The PetscFE object

2921:   Output Parameter:
2922: . sp - The PetscDualSpace object

2924:   Level: intermediate

2926: .seealso: PetscFECreate()
2927: @*/
2928: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
2929: {
2933:   *sp = fem->dualSpace;
2934:   return(0);
2935: }

2939: /*@
2940:   PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product

2942:   Not collective

2944:   Input Parameters:
2945: + fem - The PetscFE object
2946: - sp - The PetscDualSpace object

2948:   Level: intermediate

2950: .seealso: PetscFECreate()
2951: @*/
2952: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
2953: {

2959:   PetscDualSpaceDestroy(&fem->dualSpace);
2960:   fem->dualSpace = sp;
2961:   PetscObjectReference((PetscObject) fem->dualSpace);
2962:   return(0);
2963: }

2967: /*@
2968:   PetscFEGetQuadrature - Returns the PetscQuadreture used to calculate inner products

2970:   Not collective

2972:   Input Parameter:
2973: . fem - The PetscFE object

2975:   Output Parameter:
2976: . q - The PetscQuadrature object

2978:   Level: intermediate

2980: .seealso: PetscFECreate()
2981: @*/
2982: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
2983: {
2987:   *q = fem->quadrature;
2988:   return(0);
2989: }

2993: /*@
2994:   PetscFESetQuadrature - Sets the PetscQuadreture used to calculate inner products

2996:   Not collective

2998:   Input Parameters:
2999: + fem - The PetscFE object
3000: - q - The PetscQuadrature object

3002:   Level: intermediate

3004: .seealso: PetscFECreate()
3005: @*/
3006: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
3007: {

3012:   PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
3013:   PetscQuadratureDestroy(&fem->quadrature);
3014:   fem->quadrature = q;
3015:   PetscObjectReference((PetscObject) q);
3016:   return(0);
3017: }

3021: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
3022: {
3023:   const PetscInt *numDofDual;
3024:   PetscErrorCode  ierr;

3029:   PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
3030:   if (!fem->numDof) {
3031:     DM       dm;
3032:     PetscInt dim, d;

3034:     PetscDualSpaceGetDM(fem->dualSpace, &dm);
3035:     DMGetDimension(dm, &dim);
3036:     PetscMalloc1(dim+1, &fem->numDof);
3037:     for (d = 0; d <= dim; ++d) {
3038:       fem->numDof[d] = fem->numComponents*numDofDual[d];
3039:     }
3040:   }
3041:   *numDof = fem->numDof;
3042:   return(0);
3043: }

3047: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
3048: {
3049:   PetscInt         npoints = fem->quadrature->numPoints;
3050:   const PetscReal *points  = fem->quadrature->points;
3051:   PetscErrorCode   ierr;

3058:   if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
3059:   if (B) *B = fem->B;
3060:   if (D) *D = fem->D;
3061:   if (H) *H = fem->H;
3062:   return(0);
3063: }

3067: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3068: {
3069:   DM               dm;
3070:   PetscInt         pdim; /* Dimension of FE space P */
3071:   PetscInt         dim;  /* Spatial dimension */
3072:   PetscInt         comp; /* Field components */
3073:   PetscErrorCode   ierr;

3081:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3082:   DMGetDimension(dm, &dim);
3083:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3084:   PetscFEGetNumComponents(fem, &comp);
3085:   if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
3086:   if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
3087:   if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
3088:   (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
3089:   return(0);
3090: }

3094: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3095: {
3096:   DM             dm;

3101:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3102:   if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
3103:   if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
3104:   if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
3105:   return(0);
3106: }

3110: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
3111: {
3112:   PetscFE_Basic *b = (PetscFE_Basic *) fem->data;

3116:   PetscFree(b);
3117:   return(0);
3118: }

3122: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
3123: {
3124:   PetscSpace        basis;
3125:   PetscDualSpace    dual;
3126:   PetscQuadrature   q;
3127:   PetscInt          dim, Nq;
3128:   PetscViewerFormat format;
3129:   PetscErrorCode    ierr;

3132:   PetscFEGetBasisSpace(fe, &basis);
3133:   PetscFEGetDualSpace(fe, &dual);
3134:   PetscFEGetQuadrature(fe, &q);
3135:   PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
3136:   PetscViewerGetFormat(viewer, &format);
3137:   PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
3138:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
3139:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
3140:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
3141:     PetscViewerASCIIPushTab(viewer);
3142:     PetscQuadratureView(q, viewer);
3143:     PetscViewerASCIIPopTab(viewer);
3144:   } else {
3145:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
3146:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
3147:   }
3148:   PetscViewerASCIIPushTab(viewer);
3149:   PetscSpaceView(basis, viewer);
3150:   PetscDualSpaceView(dual, viewer);
3151:   PetscViewerASCIIPopTab(viewer);
3152:   return(0);
3153: }

3157: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
3158: {
3159:   PetscBool      iascii;

3165:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
3166:   if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
3167:   return(0);
3168: }

3172: /* Construct the change of basis from prime basis to nodal basis */
3173: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
3174: {
3175:   PetscReal     *work;
3176:   PetscBLASInt  *pivots;
3177: #ifndef PETSC_USE_COMPLEX
3178:   PetscBLASInt   n, info;
3179: #endif
3180:   PetscInt       pdim, j;

3184:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3185:   PetscMalloc1(pdim*pdim,&fem->invV);
3186:   for (j = 0; j < pdim; ++j) {
3187:     PetscReal      *Bf;
3188:     PetscQuadrature f;
3189:     PetscInt        q, k;

3191:     PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
3192:     PetscMalloc1(f->numPoints*pdim,&Bf);
3193:     PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
3194:     for (k = 0; k < pdim; ++k) {
3195:       /* n_j \cdot \phi_k */
3196:       fem->invV[j*pdim+k] = 0.0;
3197:       for (q = 0; q < f->numPoints; ++q) {
3198:         fem->invV[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
3199:       }
3200:     }
3201:     PetscFree(Bf);
3202:   }
3203:   PetscMalloc2(pdim,&pivots,pdim,&work);
3204: #ifndef PETSC_USE_COMPLEX
3205:   n = pdim;
3206:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, fem->invV, &n, pivots, &info));
3207:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
3208: #endif
3209:   PetscFree2(pivots,work);
3210:   return(0);
3211: }

3215: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
3216: {

3220:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
3221:   return(0);
3222: }

3226: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
3227: {
3228:   DM               dm;
3229:   PetscInt         pdim; /* Dimension of FE space P */
3230:   PetscInt         dim;  /* Spatial dimension */
3231:   PetscInt         comp; /* Field components */
3232:   PetscReal       *tmpB, *tmpD, *tmpH;
3233:   PetscInt         p, d, j, k;
3234:   PetscErrorCode   ierr;

3237:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
3238:   DMGetDimension(dm, &dim);
3239:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3240:   PetscFEGetNumComponents(fem, &comp);
3241:   /* Evaluate the prime basis functions at all points */
3242:   if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3243:   if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3244:   if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3245:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
3246:   /* Translate to the nodal basis */
3247:   for (p = 0; p < npoints; ++p) {
3248:     if (B) {
3249:       /* Multiply by V^{-1} (pdim x pdim) */
3250:       for (j = 0; j < pdim; ++j) {
3251:         const PetscInt i = (p*pdim + j)*comp;
3252:         PetscInt       c;

3254:         B[i] = 0.0;
3255:         for (k = 0; k < pdim; ++k) {
3256:           B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
3257:         }
3258:         for (c = 1; c < comp; ++c) {
3259:           B[i+c] = B[i];
3260:         }
3261:       }
3262:     }
3263:     if (D) {
3264:       /* Multiply by V^{-1} (pdim x pdim) */
3265:       for (j = 0; j < pdim; ++j) {
3266:         for (d = 0; d < dim; ++d) {
3267:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
3268:           PetscInt       c;

3270:           D[i] = 0.0;
3271:           for (k = 0; k < pdim; ++k) {
3272:             D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
3273:           }
3274:           for (c = 1; c < comp; ++c) {
3275:             D[((p*pdim + j)*comp + c)*dim + d] = D[i];
3276:           }
3277:         }
3278:       }
3279:     }
3280:     if (H) {
3281:       /* Multiply by V^{-1} (pdim x pdim) */
3282:       for (j = 0; j < pdim; ++j) {
3283:         for (d = 0; d < dim*dim; ++d) {
3284:           const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
3285:           PetscInt       c;

3287:           H[i] = 0.0;
3288:           for (k = 0; k < pdim; ++k) {
3289:             H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
3290:           }
3291:           for (c = 1; c < comp; ++c) {
3292:             H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
3293:           }
3294:         }
3295:       }
3296:     }
3297:   }
3298:   if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3299:   if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3300:   if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3301:   return(0);
3302: }

3306: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3307:                                       const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
3308: {
3309:   const PetscInt  debug = 0;
3310:   void          (*obj_func)(const PetscScalar u[], const PetscScalar u_x[], const PetscScalar u_t[], const PetscScalar a[], const PetscScalar a_x[], const PetscScalar a_t[], const PetscReal x[], PetscScalar obj[]);
3311:   PetscQuadrature quad;
3312:   PetscScalar    *u, *u_x, *a, *a_x;
3313:   PetscReal      *x;
3314:   PetscInt        dim, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, e;
3315:   PetscErrorCode  ierr;

3318:   PetscDSGetObjective(prob, field, &obj_func);
3319:   if (!obj_func) return(0);
3320:   PetscFEGetSpatialDimension(fem, &dim);
3321:   PetscFEGetQuadrature(fem, &quad);
3322:   PetscDSGetTotalDimension(prob, &totDim);
3323:   PetscDSGetTotalDimension(probAux, &totDimAux);
3324:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
3325:   PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3326:   PetscDSGetRefCoordArrays(prob, &x, NULL);
3327:   for (e = 0; e < Ne; ++e) {
3328:     const PetscReal *v0   = geom[e].v0;
3329:     const PetscReal *J    = geom[e].J;
3330:     const PetscReal *invJ = geom[e].invJ;
3331:     const PetscReal  detJ = geom[e].detJ;
3332:     const PetscReal *quadPoints, *quadWeights;
3333:     PetscInt         Nq, q;

3335:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3336:     if (debug > 1) {
3337:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3338: #ifndef PETSC_USE_COMPLEX
3339:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3340: #endif
3341:     }
3342:     for (q = 0; q < Nq; ++q) {
3343:       PetscScalar integrand;

3345:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3346:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3347:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       NULL, u, u_x, NULL);
3348:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3349:       obj_func(u, NULL, u_x, a, NULL, a_x, x, &integrand);
3350:       integrand *= detJ*quadWeights[q];
3351:       integral[field] += PetscRealPart(integrand);
3352:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", PetscRealPart(integrand), integral[field]);}
3353:     }
3354:     cOffset    += totDim;
3355:     cOffsetAux += totDimAux;
3356:   }
3357:   return(0);
3358: }

3362: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3363:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3364: {
3365:   const PetscInt  debug = 0;
3366:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
3367:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
3368:   PetscQuadrature quad;
3369:   PetscReal     **basisField, **basisFieldDer;
3370:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3371:   PetscReal      *x;
3372:   PetscInt        dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3373:   PetscErrorCode  ierr;

3376:   PetscFEGetSpatialDimension(fem, &dim);
3377:   PetscFEGetQuadrature(fem, &quad);
3378:   PetscFEGetDimension(fem, &Nb);
3379:   PetscFEGetNumComponents(fem, &Nc);
3380:   PetscDSGetTotalDimension(prob, &totDim);
3381:   PetscDSGetFieldOffset(prob, field, &fOffset);
3382:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3383:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3384:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3385:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3386:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3387:   if (probAux) {
3388:     PetscDSGetTotalDimension(probAux, &totDimAux);
3389:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3390:   }
3391:   for (e = 0; e < Ne; ++e) {
3392:     const PetscReal *v0   = geom[e].v0;
3393:     const PetscReal *J    = geom[e].J;
3394:     const PetscReal *invJ = geom[e].invJ;
3395:     const PetscReal  detJ = geom[e].detJ;
3396:     const PetscReal *quadPoints, *quadWeights;
3397:     PetscInt         Nq, q;

3399:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3400:     if (debug > 1) {
3401:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3402: #ifndef PETSC_USE_COMPLEX
3403:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3404: #endif
3405:     }
3406:     for (q = 0; q < Nq; ++q) {
3407:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3408:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3409:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3410:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3411:       f0_func(u, u_t, u_x, a, NULL, a_x, x, &f0[q*Nc]);
3412:       f1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3413:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3414:     }
3415:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3416:     cOffset    += totDim;
3417:     cOffsetAux += totDimAux;
3418:   }
3419:   return(0);
3420: }

3424: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3425:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3426: {
3427:   const PetscInt  debug = 0;
3428:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]);
3429:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]);
3430:   PetscQuadrature quad;
3431:   PetscReal     **basisField, **basisFieldDer;
3432:   PetscScalar    *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3433:   PetscReal      *x;
3434:   PetscInt        dim, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3435:   PetscErrorCode  ierr;

3438:   PetscFEGetSpatialDimension(fem, &dim);
3439:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3440:   PetscFEGetQuadrature(fem, &quad);
3441:   PetscFEGetDimension(fem, &Nb);
3442:   PetscFEGetNumComponents(fem, &Nc);
3443:   PetscDSGetTotalBdDimension(prob, &totDim);
3444:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3445:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3446:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3447:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3448:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3449:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3450:   if (probAux) {
3451:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3452:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3453:   }
3454:   for (e = 0; e < Ne; ++e) {
3455:     const PetscReal *v0          = geom[e].v0;
3456:     const PetscReal *n           = geom[e].n;
3457:     const PetscReal *J           = geom[e].J;
3458:     const PetscReal *invJ        = geom[e].invJ;
3459:     const PetscReal  detJ        = geom[e].detJ;
3460:     const PetscReal *quadPoints, *quadWeights;
3461:     PetscInt         Nq, q;

3463:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3464:     if (debug > 1) {
3465:       PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
3466: #ifndef PETSC_USE_COMPLEX
3467:       DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3468: #endif
3469:     }
3470:     for (q = 0; q < Nq; ++q) {
3471:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3472:       CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3473:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3474:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3475:       f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3476:       f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3477:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3478:     }
3479:     UpdateElementVec(dim-1, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3480:     cOffset    += totDim;
3481:     cOffsetAux += totDimAux;
3482:   }
3483:   return(0);
3484: }

3488: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3489:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3490: {
3491:   const PetscInt  debug      = 0;
3492:   void          (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g0[]);
3493:   void          (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[]);
3494:   void          (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[]);
3495:   void          (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[]);
3496:   PetscFE         fe;
3497:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3498:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3499:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3500:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3501:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3502:   PetscQuadrature quad;
3503:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3504:   PetscReal      *x;
3505:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3506:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3507:   PetscInt        dim, totDim, totDimAux, e;
3508:   PetscErrorCode  ierr;

3511:   PetscFEGetSpatialDimension(fem, &dim);
3512:   PetscFEGetQuadrature(fem, &quad);
3513:   PetscDSGetTotalDimension(prob, &totDim);
3514:   PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3515:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3516:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3517:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3518:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3519:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3520:   PetscFEGetDimension(fe, &NbI);
3521:   PetscFEGetNumComponents(fe, &NcI);
3522:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3523:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3524:   PetscFEGetDimension(fe, &NbJ);
3525:   PetscFEGetNumComponents(fe, &NcJ);
3526:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3527:   if (probAux) {
3528:     PetscDSGetTotalDimension(probAux, &totDimAux);
3529:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3530:   }
3531:   basisI    = basisField[fieldI];
3532:   basisJ    = basisField[fieldJ];
3533:   basisDerI = basisFieldDer[fieldI];
3534:   basisDerJ = basisFieldDer[fieldJ];
3535:   /* Initialize here in case the function is not defined */
3536:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3537:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3538:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3539:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3540:   for (e = 0; e < Ne; ++e) {
3541:     const PetscReal *v0   = geom[e].v0;
3542:     const PetscReal *J    = geom[e].J;
3543:     const PetscReal *invJ = geom[e].invJ;
3544:     const PetscReal  detJ = geom[e].detJ;
3545:     const PetscReal *quadPoints, *quadWeights;
3546:     PetscInt         Nq, q;

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

3552:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3553:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3554:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3555:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3556:       if (g0_func) {
3557:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3558:         g0_func(u, u_t, u_x, a, NULL, a_x, x, g0);
3559:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3560:       }
3561:       if (g1_func) {
3562:         PetscInt d, d2;
3563:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3564:         g1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3565:         for (fc = 0; fc < NcI; ++fc) {
3566:           for (gc = 0; gc < NcJ; ++gc) {
3567:             for (d = 0; d < dim; ++d) {
3568:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3569:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3570:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3571:             }
3572:           }
3573:         }
3574:       }
3575:       if (g2_func) {
3576:         PetscInt d, d2;
3577:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3578:         g2_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3579:         for (fc = 0; fc < NcI; ++fc) {
3580:           for (gc = 0; gc < NcJ; ++gc) {
3581:             for (d = 0; d < dim; ++d) {
3582:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3583:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3584:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3585:             }
3586:           }
3587:         }
3588:       }
3589:       if (g3_func) {
3590:         PetscInt d, d2, dp, d3;
3591:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3592:         g3_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3593:         for (fc = 0; fc < NcI; ++fc) {
3594:           for (gc = 0; gc < NcJ; ++gc) {
3595:             for (d = 0; d < dim; ++d) {
3596:               for (dp = 0; dp < dim; ++dp) {
3597:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3598:                 for (d2 = 0; d2 < dim; ++d2) {
3599:                   for (d3 = 0; d3 < dim; ++d3) {
3600:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3601:                   }
3602:                 }
3603:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3604:               }
3605:             }
3606:           }
3607:         }
3608:       }

3610:       for (f = 0; f < NbI; ++f) {
3611:         for (fc = 0; fc < NcI; ++fc) {
3612:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3613:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3614:           for (g = 0; g < NbJ; ++g) {
3615:             for (gc = 0; gc < NcJ; ++gc) {
3616:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3617:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3618:               const PetscInt fOff = eOffset+i*totDim+j;
3619:               PetscInt       d, d2;

3621:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3622:               for (d = 0; d < dim; ++d) {
3623:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3624:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3625:                 for (d2 = 0; d2 < dim; ++d2) {
3626:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3627:                 }
3628:               }
3629:             }
3630:           }
3631:         }
3632:       }
3633:     }
3634:     if (debug > 1) {
3635:       PetscInt fc, f, gc, g;

3637:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3638:       for (fc = 0; fc < NcI; ++fc) {
3639:         for (f = 0; f < NbI; ++f) {
3640:           const PetscInt i = offsetI + f*NcI+fc;
3641:           for (gc = 0; gc < NcJ; ++gc) {
3642:             for (g = 0; g < NbJ; ++g) {
3643:               const PetscInt j = offsetJ + g*NcJ+gc;
3644:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3645:             }
3646:           }
3647:           PetscPrintf(PETSC_COMM_SELF, "\n");
3648:         }
3649:       }
3650:     }
3651:     cOffset    += totDim;
3652:     cOffsetAux += totDimAux;
3653:     eOffset    += PetscSqr(totDim);
3654:   }
3655:   return(0);
3656: }

3660: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3661:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3662: {
3663:   const PetscInt  debug      = 0;
3664:   void          (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g0[]);
3665:   void          (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g1[]);
3666:   void          (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g2[]);
3667:   void          (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g3[]);
3668:   PetscFE         fe;
3669:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
3670:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3671:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
3672:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
3673:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
3674:   PetscQuadrature quad;
3675:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3676:   PetscReal      *x;
3677:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3678:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3679:   PetscInt        dim, bdim, totDim, totDimAux, e;
3680:   PetscErrorCode  ierr;

3683:   PetscFEGetQuadrature(fem, &quad);
3684:   PetscFEGetSpatialDimension(fem, &dim);
3685:   dim += 1; /* Spatial dimension is one higher than topological dimension */
3686:   bdim = dim-1;
3687:   PetscDSGetTotalBdDimension(prob, &totDim);
3688:   PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3689:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3690:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3691:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3692:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3693:   PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
3694:   PetscFEGetDimension(fe, &NbI);
3695:   PetscFEGetNumComponents(fe, &NcI);
3696:   PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
3697:   PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
3698:   PetscFEGetDimension(fe, &NbJ);
3699:   PetscFEGetNumComponents(fe, &NcJ);
3700:   PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
3701:   if (probAux) {
3702:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3703:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3704:   }
3705:   basisI    = basisField[fieldI];
3706:   basisJ    = basisField[fieldJ];
3707:   basisDerI = basisFieldDer[fieldI];
3708:   basisDerJ = basisFieldDer[fieldJ];
3709:   /* Initialize here in case the function is not defined */
3710:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3711:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3712:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3713:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3714:   for (e = 0; e < Ne; ++e) {
3715:     const PetscReal *v0   = geom[e].v0;
3716:     const PetscReal *n    = geom[e].n;
3717:     const PetscReal *J    = geom[e].J;
3718:     const PetscReal *invJ = geom[e].invJ;
3719:     const PetscReal  detJ = geom[e].detJ;
3720:     const PetscReal *quadPoints, *quadWeights;
3721:     PetscInt         Nq, q;

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

3727:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3728:       CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3729:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3730:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3731:       /* TODO: I think I have a mistake in the dim loops here */
3732:       if (g0_func) {
3733:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3734:         g0_func(u, u_t, u_x, a, NULL, a_x, x, n, g0);
3735:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3736:       }
3737:       if (g1_func) {
3738:         PetscInt d, d2;
3739:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3740:         g1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3741:         for (fc = 0; fc < NcI; ++fc) {
3742:           for (gc = 0; gc < NcJ; ++gc) {
3743:             for (d = 0; d < dim; ++d) {
3744:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
3745:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3746:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3747:             }
3748:           }
3749:         }
3750:       }
3751:       if (g2_func) {
3752:         PetscInt d, d2;
3753:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3754:         g2_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3755:         for (fc = 0; fc < NcI; ++fc) {
3756:           for (gc = 0; gc < NcJ; ++gc) {
3757:             for (d = 0; d < dim; ++d) {
3758:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
3759:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3760:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3761:             }
3762:           }
3763:         }
3764:       }
3765:       if (g3_func) {
3766:         PetscInt d, d2, dp, d3;
3767:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3768:         g3_func(u, u_t, u_x, a, NULL, a_x, x, n, g3);
3769:         for (fc = 0; fc < NcI; ++fc) {
3770:           for (gc = 0; gc < NcJ; ++gc) {
3771:             for (d = 0; d < dim; ++d) {
3772:               for (dp = 0; dp < dim; ++dp) {
3773:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3774:                 for (d2 = 0; d2 < dim; ++d2) {
3775:                   for (d3 = 0; d3 < dim; ++d3) {
3776:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3777:                   }
3778:                 }
3779:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3780:               }
3781:             }
3782:           }
3783:         }
3784:       }

3786:       for (f = 0; f < NbI; ++f) {
3787:         for (fc = 0; fc < NcI; ++fc) {
3788:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3789:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
3790:           for (g = 0; g < NbJ; ++g) {
3791:             for (gc = 0; gc < NcJ; ++gc) {
3792:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3793:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
3794:               const PetscInt fOff = eOffset+i*totDim+j;
3795:               PetscInt       d, d2;

3797:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3798:               for (d = 0; d < bdim; ++d) {
3799:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
3800:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
3801:                 for (d2 = 0; d2 < bdim; ++d2) {
3802:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
3803:                 }
3804:               }
3805:             }
3806:           }
3807:         }
3808:       }
3809:     }
3810:     if (debug > 1) {
3811:       PetscInt fc, f, gc, g;

3813:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3814:       for (fc = 0; fc < NcI; ++fc) {
3815:         for (f = 0; f < NbI; ++f) {
3816:           const PetscInt i = offsetI + f*NcI+fc;
3817:           for (gc = 0; gc < NcJ; ++gc) {
3818:             for (g = 0; g < NbJ; ++g) {
3819:               const PetscInt j = offsetJ + g*NcJ+gc;
3820:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3821:             }
3822:           }
3823:           PetscPrintf(PETSC_COMM_SELF, "\n");
3824:         }
3825:       }
3826:     }
3827:     cOffset    += totDim;
3828:     cOffsetAux += totDimAux;
3829:     eOffset    += PetscSqr(totDim);
3830:   }
3831:   return(0);
3832: }

3836: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
3837: {
3839:   fem->ops->setfromoptions          = NULL;
3840:   fem->ops->setup                   = PetscFESetUp_Basic;
3841:   fem->ops->view                    = PetscFEView_Basic;
3842:   fem->ops->destroy                 = PetscFEDestroy_Basic;
3843:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
3844:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
3845:   fem->ops->integrate               = PetscFEIntegrate_Basic;
3846:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
3847:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
3848:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
3849:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
3850:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
3851:   return(0);
3852: }

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

3857:   Level: intermediate

3859: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3860: M*/

3864: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
3865: {
3866:   PetscFE_Basic *b;

3871:   PetscNewLog(fem,&b);
3872:   fem->data = b;

3874:   PetscFEInitialize_Basic(fem);
3875:   return(0);
3876: }

3880: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
3881: {
3882:   PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;

3886:   PetscFree(na);
3887:   return(0);
3888: }

3892: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3893:                                                   const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3894: {
3895:   const PetscInt  debug = 0;
3896:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
3897:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
3898:   PetscQuadrature quad;
3899:   PetscReal     **basisField, **basisFieldDer;
3900:   PetscScalar    *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3901:   PetscReal      *x;
3902:   PetscInt        dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3903:   PetscErrorCode  ierr;

3906:   PetscFEGetSpatialDimension(fem, &dim);
3907:   PetscFEGetQuadrature(fem, &quad);
3908:   PetscFEGetDimension(fem, &Nb);
3909:   PetscFEGetNumComponents(fem, &Nc);
3910:   PetscDSGetTotalDimension(prob, &totDim);
3911:   PetscDSGetFieldOffset(prob, field, &fOffset);
3912:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3913:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3914:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3915:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3916:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3917:   if (probAux) {
3918:     PetscDSGetTotalDimension(probAux, &totDimAux);
3919:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3920:   }
3921:   for (e = 0; e < Ne; ++e) {
3922:     const PetscReal *quadPoints, *quadWeights;
3923:     PetscInt         Nq, q;

3925:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3926:     for (q = 0; q < Nq; ++q) {
3927:       const PetscReal *v0   = geom[e*Nq+q].v0;
3928:       const PetscReal *J    = geom[e*Nq+q].J;
3929:       const PetscReal *invJ = geom[e*Nq+q].invJ;
3930:       const PetscReal  detJ = geom[e*Nq+q].detJ;

3932:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3933:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3934:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3935:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3936:       f0_func(u, u_t, u_x, a, NULL, a_x, x, &f0[q*Nc]);
3937:       f1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3938:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3939:     }
3940:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3941:     cOffset    += totDim;
3942:     cOffsetAux += totDimAux;
3943:   }
3944:   return(0);
3945: }

3949: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3950:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3951: {
3952:   const PetscInt  debug = 0;
3953:   void          (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]);
3954:   void          (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]);
3955:   PetscQuadrature quad;
3956:   PetscReal    **basisField, **basisFieldDer;
3957:   PetscScalar    *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3958:   PetscReal      *x;
3959:   PetscInt        dim, Nb, Nc, totDim, totDimAux, 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:   PetscFEGetQuadrature(fem, &quad);
3966:   PetscFEGetDimension(fem, &Nb);
3967:   PetscFEGetNumComponents(fem, &Nc);
3968:   PetscDSGetTotalBdDimension(prob, &totDim);
3969:   PetscDSGetBdFieldOffset(prob, field, &fOffset);
3970:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3971:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3972:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3973:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3974:   PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3975:   if (probAux) {
3976:     PetscDSGetTotalBdDimension(probAux, &totDimAux);
3977:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3978:   }
3979:   for (e = 0; e < Ne; ++e) {
3980:     const PetscReal *quadPoints, *quadWeights;
3981:     PetscInt         Nq, q;

3983:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3984:     for (q = 0; q < Nq; ++q) {
3985:       const PetscReal *v0   = geom[e*Nq+q].v0;
3986:       const PetscReal *n    = geom[e*Nq+q].n;
3987:       const PetscReal *J    = geom[e*Nq+q].J;
3988:       const PetscReal *invJ = geom[e*Nq+q].invJ;
3989:       const PetscReal  detJ = geom[e*Nq+q].detJ;

3991:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
3992:       CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3993:       EvaluateFieldJets(prob,    PETSC_TRUE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
3994:       EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
3995:       f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3996:       f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3997:       TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3998:     }
3999:     UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4000:     cOffset    += totDim;
4001:     cOffsetAux += totDimAux;
4002:   }
4003:   return(0);
4004: }

4008: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
4009:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
4010: {
4011:   const PetscInt  debug      = 0;
4012:   void          (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g0[]);
4013:   void          (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[]);
4014:   void          (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[]);
4015:   void          (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[]);
4016:   PetscFE         fe;
4017:   PetscInt        cOffset    = 0; /* Offset into coefficients[] for element e */
4018:   PetscInt        cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
4019:   PetscInt        eOffset    = 0; /* Offset into elemMat[] for element e */
4020:   PetscInt        offsetI    = 0; /* Offset into an element vector for fieldI */
4021:   PetscInt        offsetJ    = 0; /* Offset into an element vector for fieldJ */
4022:   PetscQuadrature quad;
4023:   PetscScalar    *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4024:   PetscReal      *x;
4025:   PetscReal     **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
4026:   PetscInt        NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
4027:   PetscInt        dim, totDim, totDimAux, e;
4028:   PetscErrorCode  ierr;

4031:   PetscFEGetSpatialDimension(fem, &dim);
4032:   PetscFEGetQuadrature(fem, &quad);
4033:   PetscDSGetTotalDimension(prob, &totDim);
4034:   PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
4035:   PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4036:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4037:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
4038:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4039:   PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4040:   PetscFEGetDimension(fe, &NbI);
4041:   PetscFEGetNumComponents(fe, &NcI);
4042:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
4043:   PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
4044:   PetscFEGetDimension(fe, &NbJ);
4045:   PetscFEGetNumComponents(fe, &NcJ);
4046:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
4047:   if (probAux) {
4048:     PetscDSGetTotalDimension(probAux, &totDimAux);
4049:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4050:   }
4051:   basisI    = basisField[fieldI];
4052:   basisJ    = basisField[fieldJ];
4053:   basisDerI = basisFieldDer[fieldI];
4054:   basisDerJ = basisFieldDer[fieldJ];
4055:   /* Initialize here in case the function is not defined */
4056:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4057:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
4058:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
4059:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4060:   for (e = 0; e < Ne; ++e) {
4061:     const PetscReal *quadPoints, *quadWeights;
4062:     PetscInt         Nq, q;

4064:     PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4065:     for (q = 0; q < Nq; ++q) {
4066:       const PetscReal *v0   = geom[e*Nq+q].v0;
4067:       const PetscReal *J    = geom[e*Nq+q].J;
4068:       const PetscReal *invJ = geom[e*Nq+q].invJ;
4069:       const PetscReal  detJ = geom[e*Nq+q].detJ;
4070:       PetscInt         f, g, fc, gc, c;

4072:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
4073:       CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4074:       EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, &coefficients[cOffset],       &coefficients_t[cOffset], u, u_x, u_t);
4075:       EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL,                     a, a_x, NULL);
4076:       if (g0_func) {
4077:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4078:         g0_func(u, u_t, u_x, a, NULL, a_x, x, g0);
4079:         for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
4080:       }
4081:       if (g1_func) {
4082:         PetscInt d, d2;
4083:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4084:         g1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
4085:         for (fc = 0; fc < NcI; ++fc) {
4086:           for (gc = 0; gc < NcJ; ++gc) {
4087:             for (d = 0; d < dim; ++d) {
4088:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
4089:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4090:               g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4091:             }
4092:           }
4093:         }
4094:       }
4095:       if (g2_func) {
4096:         PetscInt d, d2;
4097:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4098:         g2_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
4099:         for (fc = 0; fc < NcI; ++fc) {
4100:           for (gc = 0; gc < NcJ; ++gc) {
4101:             for (d = 0; d < dim; ++d) {
4102:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
4103:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4104:               g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4105:             }
4106:           }
4107:         }
4108:       }
4109:       if (g3_func) {
4110:         PetscInt d, d2, dp, d3;
4111:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4112:         g3_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
4113:         for (fc = 0; fc < NcI; ++fc) {
4114:           for (gc = 0; gc < NcJ; ++gc) {
4115:             for (d = 0; d < dim; ++d) {
4116:               for (dp = 0; dp < dim; ++dp) {
4117:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
4118:                 for (d2 = 0; d2 < dim; ++d2) {
4119:                   for (d3 = 0; d3 < dim; ++d3) {
4120:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
4121:                   }
4122:                 }
4123:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
4124:               }
4125:             }
4126:           }
4127:         }
4128:       }

4130:       for (f = 0; f < NbI; ++f) {
4131:         for (fc = 0; fc < NcI; ++fc) {
4132:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
4133:           const PetscInt i    = offsetI+fidx; /* Element matrix row */
4134:           for (g = 0; g < NbJ; ++g) {
4135:             for (gc = 0; gc < NcJ; ++gc) {
4136:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
4137:               const PetscInt j    = offsetJ+gidx; /* Element matrix column */
4138:               const PetscInt fOff = eOffset+i*totDim+j;
4139:               PetscInt       d, d2;

4141:               elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
4142:               for (d = 0; d < dim; ++d) {
4143:                 elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
4144:                 elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
4145:                 for (d2 = 0; d2 < dim; ++d2) {
4146:                   elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
4147:                 }
4148:               }
4149:             }
4150:           }
4151:         }
4152:       }
4153:     }
4154:     if (debug > 1) {
4155:       PetscInt fc, f, gc, g;

4157:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
4158:       for (fc = 0; fc < NcI; ++fc) {
4159:         for (f = 0; f < NbI; ++f) {
4160:           const PetscInt i = offsetI + f*NcI+fc;
4161:           for (gc = 0; gc < NcJ; ++gc) {
4162:             for (g = 0; g < NbJ; ++g) {
4163:               const PetscInt j = offsetJ + g*NcJ+gc;
4164:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
4165:             }
4166:           }
4167:           PetscPrintf(PETSC_COMM_SELF, "\n");
4168:         }
4169:       }
4170:     }
4171:     cOffset    += totDim;
4172:     cOffsetAux += totDimAux;
4173:     eOffset    += PetscSqr(totDim);
4174:   }
4175:   return(0);
4176: }

4180: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
4181: {
4183:   fem->ops->setfromoptions          = NULL;
4184:   fem->ops->setup                   = PetscFESetUp_Basic;
4185:   fem->ops->view                    = NULL;
4186:   fem->ops->destroy                 = PetscFEDestroy_Nonaffine;
4187:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
4188:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
4189:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Nonaffine;
4190:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Nonaffine;
4191:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
4192:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Nonaffine;
4193:   return(0);
4194: }

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

4199:   Level: intermediate

4201: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4202: M*/

4206: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
4207: {
4208:   PetscFE_Nonaffine *na;
4209:   PetscErrorCode     ierr;

4213:   PetscNewLog(fem, &na);
4214:   fem->data = na;

4216:   PetscFEInitialize_Nonaffine(fem);
4217:   return(0);
4218: }

4220: #ifdef PETSC_HAVE_OPENCL

4224: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
4225: {
4226:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4227:   PetscErrorCode  ierr;

4230:   clReleaseCommandQueue(ocl->queue_id);
4231:   ocl->queue_id = 0;
4232:   clReleaseContext(ocl->ctx_id);
4233:   ocl->ctx_id = 0;
4234:   PetscFree(ocl);
4235:   return(0);
4236: }

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

4243: /* dim     Number of spatial dimensions:          2                   */
4244: /* N_b     Number of basis functions:             generated           */
4245: /* N_{bt}  Number of total basis functions:       N_b * N_{comp}      */
4246: /* N_q     Number of quadrature points:           generated           */
4247: /* N_{bs}  Number of block cells                  LCM(N_b, N_q)       */
4248: /* N_{bst} Number of block cell components        LCM(N_{bt}, N_q)    */
4249: /* N_{bl}  Number of concurrent blocks            generated           */
4250: /* N_t     Number of threads:                     N_{bl} * N_{bs}     */
4251: /* N_{cbc} Number of concurrent basis      cells: N_{bl} * N_q        */
4252: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b        */
4253: /* N_{sbc} Number of serial     basis      cells: N_{bs} / N_q        */
4254: /* N_{sqc} Number of serial     quadrature cells: N_{bs} / N_b        */
4255: /* N_{cb}  Number of serial cell batches:         input               */
4256: /* N_c     Number of total cells:                 N_{cb}*N_{t}/N_{comp} */
4257: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
4258: {
4259:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4260:   PetscQuadrature q;
4261:   char           *string_tail   = *string_buffer;
4262:   char           *end_of_buffer = *string_buffer + buffer_length;
4263:   char            float_str[]   = "float", double_str[]  = "double";
4264:   char           *numeric_str   = &(float_str[0]);
4265:   PetscInt        op            = ocl->op;
4266:   PetscBool       useField      = PETSC_FALSE;
4267:   PetscBool       useFieldDer   = PETSC_TRUE;
4268:   PetscBool       useFieldAux   = useAux;
4269:   PetscBool       useFieldDerAux= PETSC_FALSE;
4270:   PetscBool       useF0         = PETSC_TRUE;
4271:   PetscBool       useF1         = PETSC_TRUE;
4272:   PetscReal      *basis, *basisDer;
4273:   PetscInt        dim, N_b, N_c, N_q, N_t, p, d, b, c;
4274:   size_t          count;
4275:   PetscErrorCode  ierr;

4278:   PetscFEGetSpatialDimension(fem, &dim);
4279:   PetscFEGetDimension(fem, &N_b);
4280:   PetscFEGetNumComponents(fem, &N_c);
4281:   PetscFEGetQuadrature(fem, &q);
4282:   N_q  = q->numPoints;
4283:   N_t  = N_b * N_c * N_q * N_bl;
4284:   /* Enable device extension for double precision */
4285:   if (ocl->realType == PETSC_DOUBLE) {
4286:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4287: "#if defined(cl_khr_fp64)\n"
4288: "#  pragma OPENCL EXTENSION cl_khr_fp64: enable\n"
4289: "#elif defined(cl_amd_fp64)\n"
4290: "#  pragma OPENCL EXTENSION cl_amd_fp64: enable\n"
4291: "#endif\n",
4292:                               &count);STRING_ERROR_CHECK("Message to short");
4293:     numeric_str  = &(double_str[0]);
4294:   }
4295:   /* Kernel API */
4296:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4297: "\n"
4298: "__kernel void integrateElementQuadrature(int N_cb, __global %s *coefficients, __global %s *coefficientsAux, __global %s *jacobianInverses, __global %s *jacobianDeterminants, __global %s *elemVec)\n"
4299: "{\n",
4300:                        &count, numeric_str, numeric_str, numeric_str, numeric_str, numeric_str);STRING_ERROR_CHECK("Message to short");
4301:   /* Quadrature */
4302:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4303: "  /* Quadrature points\n"
4304: "   - (x1,y1,x2,y2,...) */\n"
4305: "  const %s points[%d] = {\n",
4306:                        &count, numeric_str, N_q*dim);STRING_ERROR_CHECK("Message to short");
4307:   for (p = 0; p < N_q; ++p) {
4308:     for (d = 0; d < dim; ++d) {
4309:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->points[p*dim+d]);STRING_ERROR_CHECK("Message to short");
4310:     }
4311:   }
4312:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4313:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4314: "  /* Quadrature weights\n"
4315: "   - (v1,v2,...) */\n"
4316: "  const %s weights[%d] = {\n",
4317:                        &count, numeric_str, N_q);STRING_ERROR_CHECK("Message to short");
4318:   for (p = 0; p < N_q; ++p) {
4319:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->weights[p]);STRING_ERROR_CHECK("Message to short");
4320:   }
4321:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4322:   /* Basis Functions */
4323:   PetscFEGetDefaultTabulation(fem, &basis, &basisDer, NULL);
4324:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4325: "  /* Nodal basis function evaluations\n"
4326: "    - basis component is fastest varying, the basis function, then point */\n"
4327: "  const %s Basis[%d] = {\n",
4328:                        &count, numeric_str, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4329:   for (p = 0; p < N_q; ++p) {
4330:     for (b = 0; b < N_b; ++b) {
4331:       for (c = 0; c < N_c; ++c) {
4332:         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");
4333:       }
4334:     }
4335:   }
4336:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4337:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4338: "\n"
4339: "  /* Nodal basis function derivative evaluations,\n"
4340: "      - derivative direction is fastest varying, then basis component, then basis function, then point */\n"
4341: "  const %s%d BasisDerivatives[%d] = {\n",
4342:                        &count, numeric_str, dim, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4343:   for (p = 0; p < N_q; ++p) {
4344:     for (b = 0; b < N_b; ++b) {
4345:       for (c = 0; c < N_c; ++c) {
4346:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "(%s%d)(", &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4347:         for (d = 0; d < dim; ++d) {
4348:           if (d > 0) {
4349:             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");
4350:           } else {
4351:             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");
4352:           }
4353:         }
4354:         PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "),\n", &count);STRING_ERROR_CHECK("Message to short");
4355:       }
4356:     }
4357:   }
4358:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4359:   /* Sizes */
4360:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4361: "  const int dim    = %d;                           // The spatial dimension\n"
4362: "  const int N_bl   = %d;                           // The number of concurrent blocks\n"
4363: "  const int N_b    = %d;                           // The number of basis functions\n"
4364: "  const int N_comp = %d;                           // The number of basis function components\n"
4365: "  const int N_bt   = N_b*N_comp;                    // The total number of scalar basis functions\n"
4366: "  const int N_q    = %d;                           // The number of quadrature points\n"
4367: "  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"
4368: "  const int N_t    = N_bst*N_bl;                    // The number of threads, N_bst * N_bl\n"
4369: "  const int N_bc   = N_t/N_comp;                    // The number of cells per batch (N_b*N_q*N_bl)\n"
4370: "  const int N_sbc  = N_bst / (N_q * N_comp);\n"
4371: "  const int N_sqc  = N_bst / N_bt;\n"
4372: "  /*const int N_c    = N_cb * N_bc;*/\n"
4373: "\n"
4374: "  /* Calculated indices */\n"
4375: "  /*const int tidx    = get_local_id(0) + get_local_size(0)*get_local_id(1);*/\n"
4376: "  const int tidx    = get_local_id(0);\n"
4377: "  const int blidx   = tidx / N_bst;                  // Block number for this thread\n"
4378: "  const int bidx    = tidx %% N_bt;                   // Basis function mapped to this thread\n"
4379: "  const int cidx    = tidx %% N_comp;                 // Basis component mapped to this thread\n"
4380: "  const int qidx    = tidx %% N_q;                    // Quadrature point mapped to this thread\n"
4381: "  const int blbidx  = tidx %% N_q + blidx*N_q;        // Cell mapped to this thread in the basis phase\n"
4382: "  const int blqidx  = tidx %% N_b + blidx*N_b;        // Cell mapped to this thread in the quadrature phase\n"
4383: "  const int gidx    = get_group_id(1)*get_num_groups(0) + get_group_id(0);\n"
4384: "  const int Goffset = gidx*N_cb*N_bc;\n",
4385:                             &count, dim, N_bl, N_b, N_c, N_q);STRING_ERROR_CHECK("Message to short");
4386:   /* Local memory */
4387:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4388: "\n"
4389: "  /* Quadrature data */\n"
4390: "  %s                w;                   // $w_q$, Quadrature weight at $x_q$\n"
4391: "  __local %s         phi_i[%d];    //[N_bt*N_q];  // $\\phi_i(x_q)$, Value of the basis function $i$ at $x_q$\n"
4392: "  __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"
4393: "  /* Geometric data */\n"
4394: "  __local %s        detJ[%d]; //[N_t];           // $|J(x_q)|$, Jacobian determinant at $x_q$\n"
4395: "  __local %s        invJ[%d];//[N_t*dim*dim];   // $J^{-1}(x_q)$, Jacobian inverse at $x_q$\n",
4396:                             &count, numeric_str, numeric_str, N_b*N_c*N_q, numeric_str, dim, N_b*N_c*N_q, numeric_str, N_t,
4397:                             numeric_str, N_t*dim*dim, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4398:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4399: "  /* FEM data */\n"
4400: "  __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",
4401:                             &count, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4402:   if (useAux) {
4403:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4404: "  __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",
4405:                             &count, numeric_str, N_t);STRING_ERROR_CHECK("Message to short");
4406:   }
4407:   if (useF0) {
4408:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4409: "  /* Intermediate calculations */\n"
4410: "  __local %s         f_0[%d]; //[N_t*N_sqc];      // $f_0(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4411:                               &count, numeric_str, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4412:   }
4413:   if (useF1) {
4414:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4415: "  __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",
4416:                               &count, numeric_str, dim, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4417:   }
4418:   /* TODO: If using elasticity, put in mu/lambda coefficients */
4419:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4420: "  /* Output data */\n"
4421: "  %s                e_i;                 // Coefficient $e_i$ of the residual\n\n",
4422:                             &count, numeric_str);STRING_ERROR_CHECK("Message to short");
4423:   /* One-time loads */
4424:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4425: "  /* These should be generated inline */\n"
4426: "  /* Load quadrature weights */\n"
4427: "  w = weights[qidx];\n"
4428: "  /* Load basis tabulation \\phi_i for this cell */\n"
4429: "  if (tidx < N_bt*N_q) {\n"
4430: "    phi_i[tidx]    = Basis[tidx];\n"
4431: "    phiDer_i[tidx] = BasisDerivatives[tidx];\n"
4432: "  }\n\n",
4433:                        &count);STRING_ERROR_CHECK("Message to short");
4434:   /* Batch loads */
4435:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4436: "  for (int batch = 0; batch < N_cb; ++batch) {\n"
4437: "    /* Load geometry */\n"
4438: "    detJ[tidx] = jacobianDeterminants[Goffset+batch*N_bc+tidx];\n"
4439: "    for (int n = 0; n < dim*dim; ++n) {\n"
4440: "      const int offset = n*N_t;\n"
4441: "      invJ[offset+tidx] = jacobianInverses[(Goffset+batch*N_bc)*dim*dim+offset+tidx];\n"
4442: "    }\n"
4443: "    /* Load coefficients u_i for this cell */\n"
4444: "    for (int n = 0; n < N_bt; ++n) {\n"
4445: "      const int offset = n*N_t;\n"
4446: "      u_i[offset+tidx] = coefficients[(Goffset*N_bt)+batch*N_t*N_b+offset+tidx];\n"
4447: "    }\n",
4448:                        &count);STRING_ERROR_CHECK("Message to short");
4449:   if (useAux) {
4450:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4451: "    /* Load coefficients a_i for this cell */\n"
4452: "    /* TODO: This should not be N_t here, it should be N_bc*N_comp_aux */\n"
4453: "    a_i[tidx] = coefficientsAux[Goffset+batch*N_t+tidx];\n",
4454:                             &count);STRING_ERROR_CHECK("Message to short");
4455:   }
4456:   /* Quadrature phase */
4457:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4458: "    barrier(CLK_LOCAL_MEM_FENCE);\n"
4459: "\n"
4460: "    /* Map coefficients to values at quadrature points */\n"
4461: "    for (int c = 0; c < N_sqc; ++c) {\n"
4462: "      const int cell          = c*N_bl*N_b + blqidx;\n"
4463: "      const int fidx          = (cell*N_q + qidx)*N_comp + cidx;\n",
4464:                        &count);STRING_ERROR_CHECK("Message to short");
4465:   if (useField) {
4466:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4467: "      %s  u[%d]; //[N_comp];     // $u(x_q)$, Value of the field at $x_q$\n",
4468:                               &count, numeric_str, N_c);STRING_ERROR_CHECK("Message to short");
4469:   }
4470:   if (useFieldDer) {
4471:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4472: "      %s%d   gradU[%d]; //[N_comp]; // $\\nabla u(x_q)$, Value of the field gradient at $x_q$\n",
4473:                               &count, numeric_str, dim, N_c);STRING_ERROR_CHECK("Message to short");
4474:   }
4475:   if (useFieldAux) {
4476:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4477: "      %s  a[%d]; //[1];     // $a(x_q)$, Value of the auxiliary fields at $x_q$\n",
4478:                               &count, numeric_str, 1);STRING_ERROR_CHECK("Message to short");
4479:   }
4480:   if (useFieldDerAux) {
4481:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4482: "      %s%d   gradA[%d]; //[1]; // $\\nabla a(x_q)$, Value of the auxiliary field gradient at $x_q$\n",
4483:                               &count, numeric_str, dim, 1);STRING_ERROR_CHECK("Message to short");
4484:   }
4485:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4486: "\n"
4487: "      for (int comp = 0; comp < N_comp; ++comp) {\n",
4488:                             &count);STRING_ERROR_CHECK("Message to short");
4489:   if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        u[comp] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4490:   if (useFieldDer) {
4491:     switch (dim) {
4492:     case 1:
4493:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "        gradU[comp].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4494:     case 2:
4495:       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;
4496:     case 3:
4497:       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;
4498:     }
4499:   }
4500:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4501: "      }\n",
4502:                             &count);STRING_ERROR_CHECK("Message to short");
4503:   if (useFieldAux) {
4504:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      a[0] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");
4505:   }
4506:   if (useFieldDerAux) {
4507:     switch (dim) {
4508:     case 1:
4509:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      gradA[0].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4510:     case 2:
4511:       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;
4512:     case 3:
4513:       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;
4514:     }
4515:   }
4516:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4517: "      /* Get field and derivatives at this quadrature point */\n"
4518: "      for (int i = 0; i < N_b; ++i) {\n"
4519: "        for (int comp = 0; comp < N_comp; ++comp) {\n"
4520: "          const int b    = i*N_comp+comp;\n"
4521: "          const int pidx = qidx*N_bt + b;\n"
4522: "          const int uidx = cell*N_bt + b;\n"
4523: "          %s%d   realSpaceDer;\n\n",
4524:                             &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4525:   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");}
4526:   if (useFieldDer) {
4527:     switch (dim) {
4528:     case 2:
4529:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4530: "          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"
4531: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4532: "          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"
4533: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n",
4534:                            &count);STRING_ERROR_CHECK("Message to short");break;
4535:     case 3:
4536:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4537: "          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"
4538: "          gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4539: "          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"
4540: "          gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n"
4541: "          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"
4542: "          gradU[comp].z += u_i[uidx]*realSpaceDer.z;\n",
4543:                            &count);STRING_ERROR_CHECK("Message to short");break;
4544:     }
4545:   }
4546:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4547: "        }\n"
4548: "      }\n",
4549:                             &count);STRING_ERROR_CHECK("Message to short");
4550:   if (useFieldAux) {
4551:     PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"          a[0] += a_i[cell];\n", &count);STRING_ERROR_CHECK("Message to short");
4552:   }
4553:   /* Calculate residual at quadrature points: Should be generated by an weak form egine */
4554:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4555: "      /* Process values at quadrature points */\n",
4556:                             &count);STRING_ERROR_CHECK("Message to short");
4557:   switch (op) {
4558:   case LAPLACIAN:
4559:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4560:     if (useF1) {
4561:       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");}
4562:       else        {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_1[fidx] = gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4563:     }
4564:     break;
4565:   case ELASTICITY:
4566:     if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "      f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4567:     if (useF1) {
4568:     switch (dim) {
4569:     case 2:
4570:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4571: "      switch (cidx) {\n"
4572: "      case 0:\n"
4573: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].x + gradU[0].x);\n"
4574: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].y + gradU[1].x);\n"
4575: "        break;\n"
4576: "      case 1:\n"
4577: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].x + gradU[0].y);\n"
4578: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].y + gradU[1].y);\n"
4579: "      }\n",
4580:                            &count);STRING_ERROR_CHECK("Message to short");break;
4581:     case 3:
4582:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4583: "      switch (cidx) {\n"
4584: "      case 0:\n"
4585: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].x + gradU[0].x);\n"
4586: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].y + gradU[1].x);\n"
4587: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].z + gradU[2].x);\n"
4588: "        break;\n"
4589: "      case 1:\n"
4590: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].x + gradU[0].y);\n"
4591: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[1].y);\n"
4592: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[2].y);\n"
4593: "        break;\n"
4594: "      case 2:\n"
4595: "        f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].x + gradU[0].z);\n"
4596: "        f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[1].z);\n"
4597: "        f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[2].z);\n"
4598: "      }\n",
4599:                            &count);STRING_ERROR_CHECK("Message to short");break;
4600:     }}
4601:     break;
4602:   default:
4603:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PDE operator %d is not supported", op);
4604:   }
4605:   if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_0[fidx] *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");}
4606:   if (useF1) {
4607:     switch (dim) {
4608:     case 1:
4609:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,"      f_1[fidx].x *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4610:     case 2:
4611:       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;
4612:     case 3:
4613:       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;
4614:     }
4615:   }
4616:   /* Thread transpose */
4617:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4618: "    }\n\n"
4619: "    /* ==== TRANSPOSE THREADS ==== */\n"
4620: "    barrier(CLK_LOCAL_MEM_FENCE);\n\n",
4621:                        &count);STRING_ERROR_CHECK("Message to short");
4622:   /* Basis phase */
4623:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4624: "    /* Map values at quadrature points to coefficients */\n"
4625: "    for (int c = 0; c < N_sbc; ++c) {\n"
4626: "      const int cell = c*N_bl*N_q + blbidx; /* Cell number in batch */\n"
4627: "\n"
4628: "      e_i = 0.0;\n"
4629: "      for (int q = 0; q < N_q; ++q) {\n"
4630: "        const int pidx = q*N_bt + bidx;\n"
4631: "        const int fidx = (cell*N_q + q)*N_comp + cidx;\n"
4632: "        %s%d   realSpaceDer;\n\n",
4633:                        &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");

4635:   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");}
4636:   if (useF1) {
4637:     switch (dim) {
4638:     case 2:
4639:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4640: "        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"
4641: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4642: "        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"
4643: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n",
4644:                            &count);STRING_ERROR_CHECK("Message to short");break;
4645:     case 3:
4646:       PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4647: "        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"
4648: "        e_i           += realSpaceDer.x*f_1[fidx].x;\n"
4649: "        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"
4650: "        e_i           += realSpaceDer.y*f_1[fidx].y;\n"
4651: "        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"
4652: "        e_i           += realSpaceDer.z*f_1[fidx].z;\n",
4653:                            &count);STRING_ERROR_CHECK("Message to short");break;
4654:     }
4655:   }
4656:   PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4657: "      }\n"
4658: "      /* Write element vector for N_{cbc} cells at a time */\n"
4659: "      elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
4660: "    }\n"
4661: "    /* ==== Could do one write per batch ==== */\n"
4662: "  }\n"
4663: "  return;\n"
4664: "}\n",
4665:                        &count);STRING_ERROR_CHECK("Message to short");
4666:   return(0);
4667: }

4671: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
4672: {
4673:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4674:   PetscInt        dim, N_bl;
4675:   PetscBool       flg;
4676:   char           *buffer;
4677:   size_t          len;
4678:   char            errMsg[8192];
4679:   cl_int          ierr2;
4680:   PetscErrorCode  ierr;

4683:   PetscFEGetSpatialDimension(fem, &dim);
4684:   PetscMalloc1(8192, &buffer);
4685:   PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
4686:   PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
4687:   PetscOptionsHasName(fem->hdr.prefix, "-petscfe_opencl_kernel_print", &flg);
4688:   if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
4689:   len  = strlen(buffer);
4690:   *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
4691:   clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
4692:   if (ierr != CL_SUCCESS) {
4693:     clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
4694:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
4695:   }
4696:   PetscFree(buffer);
4697:   *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
4698:   return(0);
4699: }

4703: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
4704: {
4705:   const PetscInt Nblocks = N/blockSize;

4708:   if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
4709:   *z = 1;
4710:   for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
4711:     *y = Nblocks / *x;
4712:     if (*x * *y == Nblocks) break;
4713:   }
4714:   if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
4715:   return(0);
4716: }

4720: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
4721: {
4722:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4723:   PetscStageLog     stageLog;
4724:   PetscEventPerfLog eventLog = NULL;
4725:   PetscInt          stage;
4726:   PetscErrorCode    ierr;

4729:   PetscLogGetStageLog(&stageLog);
4730:   PetscStageLogGetCurrent(stageLog, &stage);
4731:   PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
4732:     /* Log performance info */
4733:   eventLog->eventInfo[ocl->residualEvent].count++;
4734:   eventLog->eventInfo[ocl->residualEvent].time  += time;
4735:   eventLog->eventInfo[ocl->residualEvent].flops += flops;
4736:   return(0);
4737: }

4741: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4742:                                                const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4743: {
4744:   /* Nbc = batchSize */
4745:   PetscFE_OpenCL   *ocl = (PetscFE_OpenCL *) fem->data;
4746:   void            (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
4747:   void            (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
4748:   PetscQuadrature   q;
4749:   PetscInt          dim;
4750:   PetscInt          N_b;    /* The number of basis functions */
4751:   PetscInt          N_comp; /* The number of basis function components */
4752:   PetscInt          N_bt;   /* The total number of scalar basis functions */
4753:   PetscInt          N_q;    /* The number of quadrature points */
4754:   PetscInt          N_bst;  /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
4755:   PetscInt          N_t;    /* The number of threads, N_bst * N_bl */
4756:   PetscInt          N_bl;   /* The number of blocks */
4757:   PetscInt          N_bc;   /* The batch size, N_bl*N_q*N_b */
4758:   PetscInt          N_cb;   /* The number of batches */
4759:   PetscInt          numFlops, f0Flops = 0, f1Flops = 0;
4760:   PetscBool         useAux      = probAux ? PETSC_TRUE : PETSC_FALSE;
4761:   PetscBool         useField    = PETSC_FALSE;
4762:   PetscBool         useFieldDer = PETSC_TRUE;
4763:   PetscBool         useF0       = PETSC_TRUE;
4764:   PetscBool         useF1       = PETSC_TRUE;
4765:   /* OpenCL variables */
4766:   cl_program        ocl_prog;
4767:   cl_kernel         ocl_kernel;
4768:   cl_event          ocl_ev;         /* The event for tracking kernel execution */
4769:   cl_ulong          ns_start;       /* Nanoseconds counter on GPU at kernel start */
4770:   cl_ulong          ns_end;         /* Nanoseconds counter on GPU at kernel stop */
4771:   cl_mem            o_jacobianInverses, o_jacobianDeterminants;
4772:   cl_mem            o_coefficients, o_coefficientsAux, o_elemVec;
4773:   float            *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
4774:   double           *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
4775:   PetscReal        *r_invJ = NULL, *r_detJ = NULL;
4776:   void             *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
4777:   size_t            local_work_size[3], global_work_size[3];
4778:   size_t            realSize, x, y, z;
4779:   PetscErrorCode    ierr;

4782:   if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
4783:   PetscFEGetSpatialDimension(fem, &dim);
4784:   PetscFEGetQuadrature(fem, &q);
4785:   PetscFEGetDimension(fem, &N_b);
4786:   PetscFEGetNumComponents(fem, &N_comp);
4787:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4788:   PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
4789:   N_bt  = N_b*N_comp;
4790:   N_q   = q->numPoints;
4791:   N_bst = N_bt*N_q;
4792:   N_t   = N_bst*N_bl;
4793:   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);
4794:   /* Calculate layout */
4795:   if (Ne % (N_cb*N_bc)) { /* Remainder cells */
4796:     PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);
4797:     return(0);
4798:   }
4799:   PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
4800:   local_work_size[0]  = N_bc*N_comp;
4801:   local_work_size[1]  = 1;
4802:   local_work_size[2]  = 1;
4803:   global_work_size[0] = x * local_work_size[0];
4804:   global_work_size[1] = y * local_work_size[1];
4805:   global_work_size[2] = z * local_work_size[2];
4806:   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);
4807:   PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
4808:   /* Generate code */
4809:   if (probAux) {
4810:     PetscSpace P;
4811:     PetscInt   NfAux, order, f;

4813:     PetscDSGetNumFields(probAux, &NfAux);
4814:     for (f = 0; f < NfAux; ++f) {
4815:       PetscFE feAux;

4817:       PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
4818:       PetscFEGetBasisSpace(feAux, &P);
4819:       PetscSpaceGetOrder(P, &order);
4820:       if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
4821:     }
4822:   }
4823:   PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
4824:   /* Create buffers on the device and send data over */
4825:   PetscDataTypeGetSize(ocl->realType, &realSize);
4826:   if (sizeof(PetscReal) != realSize) {
4827:     switch (ocl->realType) {
4828:     case PETSC_FLOAT:
4829:     {
4830:       PetscInt c, b, d;

4832:       PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
4833:       for (c = 0; c < Ne; ++c) {
4834:         f_detJ[c] = (float) geom[c].detJ;
4835:         for (d = 0; d < dim*dim; ++d) {
4836:           f_invJ[c*dim*dim+d] = (float) geom[c].invJ[d];
4837:         }
4838:         for (b = 0; b < N_bt; ++b) {
4839:           f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
4840:         }
4841:       }
4842:       if (coefficientsAux) { /* Assume P0 */
4843:         for (c = 0; c < Ne; ++c) {
4844:           f_coeffAux[c] = (float) coefficientsAux[c];
4845:         }
4846:       }
4847:       oclCoeff      = (void *) f_coeff;
4848:       if (coefficientsAux) {
4849:         oclCoeffAux = (void *) f_coeffAux;
4850:       } else {
4851:         oclCoeffAux = NULL;
4852:       }
4853:       oclInvJ       = (void *) f_invJ;
4854:       oclDetJ       = (void *) f_detJ;
4855:     }
4856:     break;
4857:     case PETSC_DOUBLE:
4858:     {
4859:       PetscInt c, b, d;

4861:       PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
4862:       for (c = 0; c < Ne; ++c) {
4863:         d_detJ[c] = (double) geom[c].detJ;
4864:         for (d = 0; d < dim*dim; ++d) {
4865:           d_invJ[c*dim*dim+d] = (double) geom[c].invJ[d];
4866:         }
4867:         for (b = 0; b < N_bt; ++b) {
4868:           d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
4869:         }
4870:       }
4871:       if (coefficientsAux) { /* Assume P0 */
4872:         for (c = 0; c < Ne; ++c) {
4873:           d_coeffAux[c] = (double) coefficientsAux[c];
4874:         }
4875:       }
4876:       oclCoeff      = (void *) d_coeff;
4877:       if (coefficientsAux) {
4878:         oclCoeffAux = (void *) d_coeffAux;
4879:       } else {
4880:         oclCoeffAux = NULL;
4881:       }
4882:       oclInvJ       = (void *) d_invJ;
4883:       oclDetJ       = (void *) d_detJ;
4884:     }
4885:     break;
4886:     default:
4887:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4888:     }
4889:   } else {
4890:     PetscInt c, d;

4892:     PetscMalloc2(Ne*dim*dim,&r_invJ,Ne,&r_detJ);
4893:     for (c = 0; c < Ne; ++c) {
4894:       r_detJ[c] = geom[c].detJ;
4895:       for (d = 0; d < dim*dim; ++d) {
4896:         r_invJ[c*dim*dim+d] = geom[c].invJ[d];
4897:       }
4898:     }
4899:     oclCoeff    = (void *) coefficients;
4900:     oclCoeffAux = (void *) coefficientsAux;
4901:     oclInvJ     = (void *) r_invJ;
4902:     oclDetJ     = (void *) r_detJ;
4903:   }
4904:   o_coefficients         = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt    * realSize, oclCoeff,    &ierr);
4905:   if (coefficientsAux) {
4906:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclCoeffAux, &ierr);
4907:   } else {
4908:     o_coefficientsAux    = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY,                        Ne         * realSize, oclCoeffAux, &ierr);
4909:   }
4910:   o_jacobianInverses     = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ,     &ierr);
4911:   o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne         * realSize, oclDetJ,     &ierr);
4912:   o_elemVec              = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY,                       Ne*N_bt    * realSize, NULL,        &ierr);
4913:   /* Kernel launch */
4914:   clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
4915:   clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
4916:   clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
4917:   clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
4918:   clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
4919:   clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
4920:   clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
4921:   /* Read data back from device */
4922:   if (sizeof(PetscReal) != realSize) {
4923:     switch (ocl->realType) {
4924:     case PETSC_FLOAT:
4925:     {
4926:       float   *elem;
4927:       PetscInt c, b;

4929:       PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
4930:       PetscMalloc1(Ne*N_bt, &elem);
4931:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4932:       for (c = 0; c < Ne; ++c) {
4933:         for (b = 0; b < N_bt; ++b) {
4934:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4935:         }
4936:       }
4937:       PetscFree(elem);
4938:     }
4939:     break;
4940:     case PETSC_DOUBLE:
4941:     {
4942:       double  *elem;
4943:       PetscInt c, b;

4945:       PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
4946:       PetscMalloc1(Ne*N_bt, &elem);
4947:       clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4948:       for (c = 0; c < Ne; ++c) {
4949:         for (b = 0; b < N_bt; ++b) {
4950:           elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4951:         }
4952:       }
4953:       PetscFree(elem);
4954:     }
4955:     break;
4956:     default:
4957:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4958:     }
4959:   } else {
4960:     PetscFree2(r_invJ,r_detJ);
4961:     clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
4962:   }
4963:   /* Log performance */
4964:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
4965:   clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END,   sizeof(cl_ulong), &ns_end,   NULL);
4966:   f0Flops = 0;
4967:   switch (ocl->op) {
4968:   case LAPLACIAN:
4969:     f1Flops = useAux ? dim : 0;break;
4970:   case ELASTICITY:
4971:     f1Flops = 2*dim*dim;break;
4972:   }
4973:   numFlops = Ne*(
4974:     N_q*(
4975:       N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
4976:       /*+
4977:        N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
4978:       +
4979:       N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
4980:     +
4981:     N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
4982:   PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
4983:   /* Cleanup */
4984:   clReleaseMemObject(o_coefficients);
4985:   clReleaseMemObject(o_coefficientsAux);
4986:   clReleaseMemObject(o_jacobianInverses);
4987:   clReleaseMemObject(o_jacobianDeterminants);
4988:   clReleaseMemObject(o_elemVec);
4989:   clReleaseKernel(ocl_kernel);
4990:   clReleaseProgram(ocl_prog);
4991:   return(0);
4992: }

4996: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
4997: {
4999:   fem->ops->setfromoptions          = NULL;
5000:   fem->ops->setup                   = PetscFESetUp_Basic;
5001:   fem->ops->view                    = NULL;
5002:   fem->ops->destroy                 = PetscFEDestroy_OpenCL;
5003:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
5004:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
5005:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_OpenCL;
5006:   fem->ops->integratebdresidual     = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
5007:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
5008:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
5009:   return(0);
5010: }

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

5015:   Level: intermediate

5017: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5018: M*/

5022: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
5023: {
5024:   PetscFE_OpenCL *ocl;
5025:   cl_uint         num_platforms;
5026:   cl_platform_id  platform_ids[42];
5027:   cl_uint         num_devices;
5028:   cl_device_id    device_ids[42];
5029:   cl_int          ierr2;
5030:   PetscErrorCode  ierr;

5034:   PetscNewLog(fem,&ocl);
5035:   fem->data = ocl;

5037:   /* Init Platform */
5038:   clGetPlatformIDs(42, platform_ids, &num_platforms);
5039:   if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
5040:   ocl->pf_id = platform_ids[0];
5041:   /* Init Device */
5042:   clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
5043:   if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
5044:   ocl->dev_id = device_ids[0];
5045:   /* Create context with one command queue */
5046:   ocl->ctx_id   = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
5047:   ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
5048:   /* Types */
5049:   ocl->realType = PETSC_FLOAT;
5050:   /* Register events */
5051:   PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
5052:   /* Equation handling */
5053:   ocl->op = LAPLACIAN;

5055:   PetscFEInitialize_OpenCL(fem);
5056:   return(0);
5057: }

5061: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
5062: {
5063:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

5067:   ocl->realType = realType;
5068:   return(0);
5069: }

5073: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
5074: {
5075:   PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;

5080:   *realType = ocl->realType;
5081:   return(0);
5082: }

5084: #endif /* PETSC_HAVE_OPENCL */

5088: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
5089: {
5090:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5091:   PetscErrorCode     ierr;

5094:   CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5095:   PetscFree(cmp->embedding);
5096:   PetscFree(cmp);
5097:   return(0);
5098: }

5102: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
5103: {
5104:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5105:   DM                 K;
5106:   PetscReal         *work, *subpoint;
5107:   PetscBLASInt      *pivots;
5108: #ifndef PETSC_USE_COMPLEX
5109:   PetscBLASInt       n, info;
5110: #endif
5111:   PetscInt           dim, pdim, spdim, j, s;
5112:   PetscErrorCode     ierr;

5115:   /* Get affine mapping from reference cell to each subcell */
5116:   PetscDualSpaceGetDM(fem->dualSpace, &K);
5117:   DMGetDimension(K, &dim);
5118:   DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
5119:   CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5120:   /* Determine dof embedding into subelements */
5121:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5122:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
5123:   PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
5124:   DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
5125:   for (s = 0; s < cmp->numSubelements; ++s) {
5126:     PetscInt sd = 0;

5128:     for (j = 0; j < pdim; ++j) {
5129:       PetscBool       inside;
5130:       PetscQuadrature f;
5131:       PetscInt        d, e;

5133:       PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
5134:       /* Apply transform to first point, and check that point is inside subcell */
5135:       for (d = 0; d < dim; ++d) {
5136:         subpoint[d] = -1.0;
5137:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
5138:       }
5139:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5140:       if (inside) {cmp->embedding[s*spdim+sd++] = j;}
5141:     }
5142:     if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
5143:   }
5144:   DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
5145:   /* Construct the change of basis from prime basis to nodal basis for each subelement */
5146:   PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
5147:   PetscMalloc2(spdim,&pivots,spdim,&work);
5148:   for (s = 0; s < cmp->numSubelements; ++s) {
5149:     for (j = 0; j < spdim; ++j) {
5150:       PetscReal      *Bf;
5151:       PetscQuadrature f;
5152:       PetscInt        q, k;

5154:       PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
5155:       PetscMalloc1(f->numPoints*spdim,&Bf);
5156:       PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
5157:       for (k = 0; k < spdim; ++k) {
5158:         /* n_j \cdot \phi_k */
5159:         fem->invV[(s*spdim + j)*spdim+k] = 0.0;
5160:         for (q = 0; q < f->numPoints; ++q) {
5161:           fem->invV[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
5162:         }
5163:       }
5164:       PetscFree(Bf);
5165:     }
5166: #ifndef PETSC_USE_COMPLEX
5167:     n = spdim;
5168:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &fem->invV[s*spdim*spdim], &n, pivots, &info));
5169:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &fem->invV[s*spdim*spdim], &n, pivots, work, &n, &info));
5170: #endif
5171:   }
5172:   PetscFree2(pivots,work);
5173:   return(0);
5174: }

5178: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
5179: {
5180:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5181:   DM                 dm;
5182:   PetscInt           pdim;  /* Dimension of FE space P */
5183:   PetscInt           spdim; /* Dimension of subelement FE space P */
5184:   PetscInt           dim;   /* Spatial dimension */
5185:   PetscInt           comp;  /* Field components */
5186:   PetscInt          *subpoints;
5187:   PetscReal         *tmpB, *tmpD, *tmpH, *subpoint;
5188:   PetscInt           p, s, d, e, j, k;
5189:   PetscErrorCode     ierr;

5192:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
5193:   DMGetDimension(dm, &dim);
5194:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
5195:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5196:   PetscFEGetNumComponents(fem, &comp);
5197:   /* Divide points into subelements */
5198:   DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
5199:   DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
5200:   for (p = 0; p < npoints; ++p) {
5201:     for (s = 0; s < cmp->numSubelements; ++s) {
5202:       PetscBool inside;

5204:       /* Apply transform, and check that point is inside cell */
5205:       for (d = 0; d < dim; ++d) {
5206:         subpoint[d] = -1.0;
5207:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
5208:       }
5209:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5210:       if (inside) {subpoints[p] = s; break;}
5211:     }
5212:     if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
5213:   }
5214:   DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
5215:   /* Evaluate the prime basis functions at all points */
5216:   if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5217:   if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5218:   if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5219:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
5220:   /* Translate to the nodal basis */
5221:   if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
5222:   if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
5223:   if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
5224:   for (p = 0; p < npoints; ++p) {
5225:     const PetscInt s = subpoints[p];

5227:     if (B) {
5228:       /* Multiply by V^{-1} (spdim x spdim) */
5229:       for (j = 0; j < spdim; ++j) {
5230:         const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
5231:         PetscInt       c;

5233:         B[i] = 0.0;
5234:         for (k = 0; k < spdim; ++k) {
5235:           B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
5236:         }
5237:         for (c = 1; c < comp; ++c) {
5238:           B[i+c] = B[i];
5239:         }
5240:       }
5241:     }
5242:     if (D) {
5243:       /* Multiply by V^{-1} (spdim x spdim) */
5244:       for (j = 0; j < spdim; ++j) {
5245:         for (d = 0; d < dim; ++d) {
5246:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
5247:           PetscInt       c;

5249:           D[i] = 0.0;
5250:           for (k = 0; k < spdim; ++k) {
5251:             D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
5252:           }
5253:           for (c = 1; c < comp; ++c) {
5254:             D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
5255:           }
5256:         }
5257:       }
5258:     }
5259:     if (H) {
5260:       /* Multiply by V^{-1} (pdim x pdim) */
5261:       for (j = 0; j < spdim; ++j) {
5262:         for (d = 0; d < dim*dim; ++d) {
5263:           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
5264:           PetscInt       c;

5266:           H[i] = 0.0;
5267:           for (k = 0; k < spdim; ++k) {
5268:             H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
5269:           }
5270:           for (c = 1; c < comp; ++c) {
5271:             H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
5272:           }
5273:         }
5274:       }
5275:     }
5276:   }
5277:   DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
5278:   if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5279:   if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5280:   if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5281:   return(0);
5282: }

5286: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
5287: {
5289:   fem->ops->setfromoptions          = NULL;
5290:   fem->ops->setup                   = PetscFESetUp_Composite;
5291:   fem->ops->view                    = NULL;
5292:   fem->ops->destroy                 = PetscFEDestroy_Composite;
5293:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
5294:   fem->ops->gettabulation           = PetscFEGetTabulation_Composite;
5295:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
5296:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
5297:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
5298:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
5299:   return(0);
5300: }

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

5305:   Level: intermediate

5307: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5308: M*/

5312: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
5313: {
5314:   PetscFE_Composite *cmp;
5315:   PetscErrorCode     ierr;

5319:   PetscNewLog(fem, &cmp);
5320:   fem->data = cmp;

5322:   cmp->cellRefiner    = REFINER_NOOP;
5323:   cmp->numSubelements = -1;
5324:   cmp->v0             = NULL;
5325:   cmp->jac            = NULL;

5327:   PetscFEInitialize_Composite(fem);
5328:   return(0);
5329: }

5333: PetscErrorCode PetscFECompositeExpandQuadrature(PetscFE fem, PetscQuadrature q, PetscQuadrature *qref)
5334: {
5335:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5336:   const PetscReal   *points,    *weights;
5337:   PetscReal         *pointsRef, *weightsRef;
5338:   PetscInt           dim, order, npoints, npointsRef, c, p, d, e;
5339:   PetscErrorCode     ierr;

5345:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
5346:   PetscQuadratureGetOrder(q, &order);
5347:   PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
5348:   npointsRef = npoints*cmp->numSubelements;
5349:   PetscMalloc1(npointsRef*dim,&pointsRef);
5350:   PetscMalloc1(npointsRef,&weightsRef);
5351:   for (c = 0; c < cmp->numSubelements; ++c) {
5352:     for (p = 0; p < npoints; ++p) {
5353:       for (d = 0; d < dim; ++d) {
5354:         pointsRef[(c*npoints + p)*dim+d] = cmp->v0[c*dim+d];
5355:         for (e = 0; e < dim; ++e) {
5356:           pointsRef[(c*npoints + p)*dim+d] += cmp->jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
5357:         }
5358:       }
5359:       /* Could also use detJ here */
5360:       weightsRef[c*npoints+p] = weights[p]/cmp->numSubelements;
5361:     }
5362:   }
5363:   PetscQuadratureSetOrder(*qref, order);
5364:   PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
5365:   return(0);
5366: }

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

5373:   Not collective

5375:   Input Parameter:
5376: . fe - The PetscFE

5378:   Output Parameter:
5379: . dim - The dimension

5381:   Level: intermediate

5383: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
5384: @*/
5385: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
5386: {

5392:   if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
5393:   return(0);
5394: }

5396: /*
5397: Purpose: Compute element vector for chunk of elements

5399: Input:
5400:   Sizes:
5401:      Ne:  number of elements
5402:      Nf:  number of fields
5403:      PetscFE
5404:        dim: spatial dimension
5405:        Nb:  number of basis functions
5406:        Nc:  number of field components
5407:        PetscQuadrature
5408:          Nq:  number of quadrature points

5410:   Geometry:
5411:      PetscFECellGeom[Ne] possibly *Nq
5412:        PetscReal v0s[dim]
5413:        PetscReal n[dim]
5414:        PetscReal jacobians[dim*dim]
5415:        PetscReal jacobianInverses[dim*dim]
5416:        PetscReal jacobianDeterminants
5417:   FEM:
5418:      PetscFE
5419:        PetscQuadrature
5420:          PetscReal   quadPoints[Nq*dim]
5421:          PetscReal   quadWeights[Nq]
5422:        PetscReal   basis[Nq*Nb*Nc]
5423:        PetscReal   basisDer[Nq*Nb*Nc*dim]
5424:      PetscScalar coefficients[Ne*Nb*Nc]
5425:      PetscScalar elemVec[Ne*Nb*Nc]

5427:   Problem:
5428:      PetscInt f: the active field
5429:      f0, f1

5431:   Work Space:
5432:      PetscFE
5433:        PetscScalar f0[Nq*dim];
5434:        PetscScalar f1[Nq*dim*dim];
5435:        PetscScalar u[Nc];
5436:        PetscScalar gradU[Nc*dim];
5437:        PetscReal   x[dim];
5438:        PetscScalar realSpaceDer[dim];

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

5442: Input:
5443:   Sizes:
5444:      N_cb: Number of serial cell batches

5446:   Geometry:
5447:      PetscReal v0s[Ne*dim]
5448:      PetscReal jacobians[Ne*dim*dim]        possibly *Nq
5449:      PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5450:      PetscReal jacobianDeterminants[Ne]     possibly *Nq
5451:   FEM:
5452:      static PetscReal   quadPoints[Nq*dim]
5453:      static PetscReal   quadWeights[Nq]
5454:      static PetscReal   basis[Nq*Nb*Nc]
5455:      static PetscReal   basisDer[Nq*Nb*Nc*dim]
5456:      PetscScalar coefficients[Ne*Nb*Nc]
5457:      PetscScalar elemVec[Ne*Nb*Nc]

5459: ex62.c:
5460:   PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
5461:                                                const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
5462:                                                void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
5463:                                                void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])

5465: ex52.c:
5466:   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)
5467:   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)

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

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

5476: ex52_integrateElementOpenCL.c:
5477: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
5478:                                                      const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5479:                                                      PetscLogEvent event, PetscInt debug, PetscInt pde_op)

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

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

5489:   Not collective

5491:   Input Parameters:
5492: + fem          - The PetscFE object for the field being integrated
5493: . prob         - The PetscDS specifing the discretizations and continuum functions
5494: . field        - The field being integrated
5495: . Ne           - The number of elements in the chunk
5496: . geom         - The cell geometry for each cell in the chunk
5497: . coefficients - The array of FEM basis coefficients for the elements
5498: . probAux      - The PetscDS specifing the auxiliary discretizations
5499: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5501:   Output Parameter
5502: . integral     - the integral for this field

5504:   Level: developer

5506: .seealso: PetscFEIntegrateResidual()
5507: @*/
5508: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5509:                                 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
5510: {

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

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

5525:   Not collective

5527:   Input Parameters:
5528: + fem          - The PetscFE object for the field being integrated
5529: . prob         - The PetscDS specifing the discretizations and continuum functions
5530: . field        - The field being integrated
5531: . Ne           - The number of elements in the chunk
5532: . geom         - The cell geometry for each cell in the chunk
5533: . coefficients - The array of FEM basis coefficients for the elements
5534: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5535: . probAux      - The PetscDS specifing the auxiliary discretizations
5536: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5538:   Output Parameter
5539: . elemVec      - the element residual vectors from each element

5541:   Note:
5542: $ Loop over batch of elements (e):
5543: $   Loop over quadrature points (q):
5544: $     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
5545: $     Call f_0 and f_1
5546: $   Loop over element vector entries (f,fc --> i):
5547: $     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)

5549:   Level: developer

5551: .seealso: PetscFEIntegrateResidual()
5552: @*/
5553: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5554:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5555: {

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

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

5570:   Not collective

5572:   Input Parameters:
5573: + fem          - The PetscFE object for the field being integrated
5574: . prob         - The PetscDS specifing the discretizations and continuum functions
5575: . field        - The field being integrated
5576: . Ne           - The number of elements in the chunk
5577: . geom         - The cell geometry for each cell in the chunk
5578: . coefficients - The array of FEM basis coefficients for the elements
5579: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5580: . probAux      - The PetscDS specifing the auxiliary discretizations
5581: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5583:   Output Parameter
5584: . elemVec      - the element residual vectors from each element

5586:   Level: developer

5588: .seealso: PetscFEIntegrateResidual()
5589: @*/
5590: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5591:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5592: {

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

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

5606:   Not collective

5608:   Input Parameters:
5609: + fem          = The PetscFE object for the field being integrated
5610: . prob         - The PetscDS specifing the discretizations and continuum functions
5611: . fieldI       - The test field being integrated
5612: . fieldJ       - The basis field being integrated
5613: . Ne           - The number of elements in the chunk
5614: . geom         - The cell geometry for each cell in the chunk
5615: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5616: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5617: . probAux      - The PetscDS specifing the auxiliary discretizations
5618: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5620:   Output Parameter
5621: . elemMat      - the element matrices for the Jacobian from each element

5623:   Note:
5624: $ Loop over batch of elements (e):
5625: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5626: $     Loop over quadrature points (q):
5627: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5628: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5629: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5630: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5631: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5632: */
5633: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5634:                                         const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5635: {

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

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

5649:   Not collective

5651:   Input Parameters:
5652: + fem          = The PetscFE object for the field being integrated
5653: . prob         - The PetscDS specifing the discretizations and continuum functions
5654: . fieldI       - The test field being integrated
5655: . fieldJ       - The basis field being integrated
5656: . Ne           - The number of elements in the chunk
5657: . geom         - The cell geometry for each cell in the chunk
5658: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5659: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5660: . probAux      - The PetscDS specifing the auxiliary discretizations
5661: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements

5663:   Output Parameter
5664: . elemMat              - the element matrices for the Jacobian from each element

5666:   Note:
5667: $ Loop over batch of elements (e):
5668: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
5669: $     Loop over quadrature points (q):
5670: $       Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5671: $         elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5672: $                      + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5673: $                      + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5674: $                      + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5675: */
5676: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5677:                                           const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5678: {

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

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

5694:   Input Parameter:
5695: . fe - The initial PetscFE

5697:   Output Parameter:
5698: . feRef - The refined PetscFE

5700:   Level: developer

5702: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5703: @*/
5704: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
5705: {
5706:   PetscSpace       P, Pref;
5707:   PetscDualSpace   Q, Qref;
5708:   DM               K, Kref;
5709:   PetscQuadrature  q, qref;
5710:   PetscInt         numComp;
5711:   PetscErrorCode   ierr;

5714:   PetscFEGetBasisSpace(fe, &P);
5715:   PetscFEGetDualSpace(fe, &Q);
5716:   PetscFEGetQuadrature(fe, &q);
5717:   PetscDualSpaceGetDM(Q, &K);
5718:   /* Create space */
5719:   PetscObjectReference((PetscObject) P);
5720:   Pref = P;
5721:   /* Create dual space */
5722:   PetscDualSpaceDuplicate(Q, &Qref);
5723:   DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
5724:   PetscDualSpaceSetDM(Qref, Kref);
5725:   DMDestroy(&Kref);
5726:   PetscDualSpaceSetUp(Qref);
5727:   /* Create element */
5728:   PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
5729:   PetscFESetType(*feRef, PETSCFECOMPOSITE);
5730:   PetscFESetBasisSpace(*feRef, Pref);
5731:   PetscFESetDualSpace(*feRef, Qref);
5732:   PetscFEGetNumComponents(fe,    &numComp);
5733:   PetscFESetNumComponents(*feRef, numComp);
5734:   PetscFESetUp(*feRef);
5735:   PetscSpaceDestroy(&Pref);
5736:   PetscDualSpaceDestroy(&Qref);
5737:   /* Create quadrature */
5738:   PetscFECompositeExpandQuadrature(*feRef, q, &qref);
5739:   PetscFESetQuadrature(*feRef, qref);
5740:   PetscQuadratureDestroy(&qref);
5741:   return(0);
5742: }

5746: /*@
5747:   PetscFECreateDefault - Create a PetscFE for basic FEM computation

5749:   Collective on DM

5751:   Input Parameters:
5752: + dm         - The underlying DM for the domain
5753: . dim        - The spatial dimension
5754: . numComp    - The number of components
5755: . isSimplex  - Flag for simplex reference cell, otherwise its a tensor product
5756: . prefix     - The options prefix, or NULL
5757: - qorder     - The quadrature order

5759:   Output Parameter:
5760: . fem - The PetscFE object

5762:   Level: beginner

5764: .keywords: PetscFE, finite element
5765: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
5766: @*/
5767: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
5768: {
5769:   PetscQuadrature q;
5770:   DM              K;
5771:   PetscSpace      P;
5772:   PetscDualSpace  Q;
5773:   PetscInt        order;
5774:   PetscErrorCode  ierr;

5777:   /* Create space */
5778:   PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
5779:   PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
5780:   PetscSpaceSetFromOptions(P);
5781:   PetscSpacePolynomialSetNumVariables(P, dim);
5782:   PetscSpaceSetUp(P);
5783:   PetscSpaceGetOrder(P, &order);
5784:   /* Create dual space */
5785:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
5786:   PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
5787:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
5788:   PetscDualSpaceSetDM(Q, K);
5789:   DMDestroy(&K);
5790:   PetscDualSpaceSetOrder(Q, order);
5791:   PetscDualSpaceSetFromOptions(Q);
5792:   PetscDualSpaceSetUp(Q);
5793:   /* Create element */
5794:   PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
5795:   PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
5796:   PetscFESetFromOptions(*fem);
5797:   PetscFESetBasisSpace(*fem, P);
5798:   PetscFESetDualSpace(*fem, Q);
5799:   PetscFESetNumComponents(*fem, numComp);
5800:   PetscFESetUp(*fem);
5801:   PetscSpaceDestroy(&P);
5802:   PetscDualSpaceDestroy(&Q);
5803:   /* Create quadrature (with specified order if given) */
5804:   if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5805:   else           {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5806:   PetscFESetQuadrature(*fem, q);
5807:   PetscQuadratureDestroy(&q);
5808:   return(0);
5809: }