xref: /petsc/src/dm/impls/plex/plexfem.c (revision f15a11d28b36f01875e426cb22da92fd19b081e4) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 #include <petscfe.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "DMPlexGetScale"
7 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
8 {
9   DM_Plex *mesh = (DM_Plex*) dm->data;
10 
11   PetscFunctionBegin;
12   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
13   PetscValidPointer(scale, 3);
14   *scale = mesh->scale[unit];
15   PetscFunctionReturn(0);
16 }
17 
18 #undef __FUNCT__
19 #define __FUNCT__ "DMPlexSetScale"
20 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
21 {
22   DM_Plex *mesh = (DM_Plex*) dm->data;
23 
24   PetscFunctionBegin;
25   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
26   mesh->scale[unit] = scale;
27   PetscFunctionReturn(0);
28 }
29 
30 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
31 {
32   switch (i) {
33   case 0:
34     switch (j) {
35     case 0: return 0;
36     case 1:
37       switch (k) {
38       case 0: return 0;
39       case 1: return 0;
40       case 2: return 1;
41       }
42     case 2:
43       switch (k) {
44       case 0: return 0;
45       case 1: return -1;
46       case 2: return 0;
47       }
48     }
49   case 1:
50     switch (j) {
51     case 0:
52       switch (k) {
53       case 0: return 0;
54       case 1: return 0;
55       case 2: return -1;
56       }
57     case 1: return 0;
58     case 2:
59       switch (k) {
60       case 0: return 1;
61       case 1: return 0;
62       case 2: return 0;
63       }
64     }
65   case 2:
66     switch (j) {
67     case 0:
68       switch (k) {
69       case 0: return 0;
70       case 1: return 1;
71       case 2: return 0;
72       }
73     case 1:
74       switch (k) {
75       case 0: return -1;
76       case 1: return 0;
77       case 2: return 0;
78       }
79     case 2: return 0;
80     }
81   }
82   return 0;
83 }
84 
85 #undef __FUNCT__
86 #define __FUNCT__ "DMPlexCreateRigidBody"
87 /*@C
88   DMPlexCreateRigidBody - create rigid body modes from coordinates
89 
90   Collective on DM
91 
92   Input Arguments:
93 + dm - the DM
94 . section - the local section associated with the rigid field, or NULL for the default section
95 - globalSection - the global section associated with the rigid field, or NULL for the default section
96 
97   Output Argument:
98 . sp - the null space
99 
100   Note: This is necessary to take account of Dirichlet conditions on the displacements
101 
102   Level: advanced
103 
104 .seealso: MatNullSpaceCreate()
105 @*/
106 PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
107 {
108   MPI_Comm       comm;
109   Vec            coordinates, localMode, mode[6];
110   PetscSection   coordSection;
111   PetscScalar   *coords;
112   PetscInt       dim, vStart, vEnd, v, n, m, d, i, j;
113   PetscErrorCode ierr;
114 
115   PetscFunctionBegin;
116   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
117   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
118   if (dim == 1) {
119     ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr);
120     PetscFunctionReturn(0);
121   }
122   if (!section)       {ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);}
123   if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);}
124   ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr);
125   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
126   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
127   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
128   m    = (dim*(dim+1))/2;
129   ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr);
130   ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr);
131   ierr = VecSetUp(mode[0]);CHKERRQ(ierr);
132   for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);}
133   /* Assume P1 */
134   ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr);
135   for (d = 0; d < dim; ++d) {
136     PetscScalar values[3] = {0.0, 0.0, 0.0};
137 
138     values[d] = 1.0;
139     ierr      = VecSet(localMode, 0.0);CHKERRQ(ierr);
140     for (v = vStart; v < vEnd; ++v) {
141       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
142     }
143     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
144     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
145   }
146   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
147   for (d = dim; d < dim*(dim+1)/2; ++d) {
148     PetscInt i, j, k = dim > 2 ? d - dim : d;
149 
150     ierr = VecSet(localMode, 0.0);CHKERRQ(ierr);
151     for (v = vStart; v < vEnd; ++v) {
152       PetscScalar values[3] = {0.0, 0.0, 0.0};
153       PetscInt    off;
154 
155       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
156       for (i = 0; i < dim; ++i) {
157         for (j = 0; j < dim; ++j) {
158           values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
159         }
160       }
161       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
162     }
163     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
164     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
165   }
166   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
167   ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr);
168   for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);}
169   /* Orthonormalize system */
170   for (i = dim; i < m; ++i) {
171     PetscScalar dots[6];
172 
173     ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr);
174     for (j = 0; j < i; ++j) dots[j] *= -1.0;
175     ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr);
176     ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);
177   }
178   ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr);
179   for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);}
180   PetscFunctionReturn(0);
181 }
182 
183 #undef __FUNCT__
184 #define __FUNCT__ "DMPlexProjectFunctionLocal"
185 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
186 {
187   PetscDualSpace *sp;
188   PetscSection    section;
189   PetscScalar    *values;
190   PetscReal      *v0, *J, detJ;
191   PetscInt        numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, c, f, d, v, comp;
192   PetscErrorCode  ierr;
193 
194   PetscFunctionBegin;
195   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
196   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
197   ierr = PetscMalloc1(numFields, &sp);CHKERRQ(ierr);
198   for (f = 0; f < numFields; ++f) {
199     ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr);
200     ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
201     ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
202     totDim += spDim*numComp;
203   }
204   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
205   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
206   ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr);
207   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
208   ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
209   ierr = PetscMalloc2(dim,&v0,dim*dim,&J);CHKERRQ(ierr);
210   for (c = cStart; c < cEnd; ++c) {
211     PetscCellGeometry geom;
212 
213     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
214     geom.v0   = v0;
215     geom.J    = J;
216     geom.detJ = &detJ;
217     for (f = 0, v = 0; f < numFields; ++f) {
218       void * const ctx = ctxs ? ctxs[f] : NULL;
219       ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
220       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
221       for (d = 0; d < spDim; ++d) {
222         if (funcs[f]) {
223           ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr);
224         } else {
225           for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
226         }
227         v += numComp;
228       }
229     }
230     ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr);
231   }
232   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
233   ierr = PetscFree2(v0,J);CHKERRQ(ierr);
234   ierr = PetscFree(sp);CHKERRQ(ierr);
235   PetscFunctionReturn(0);
236 }
237 
238 #undef __FUNCT__
239 #define __FUNCT__ "DMPlexProjectFunction"
240 /*@C
241   DMPlexProjectFunction - This projects the given function into the function space provided.
242 
243   Input Parameters:
244 + dm      - The DM
245 . fe      - The PetscFE associated with the field
246 . funcs   - The coordinate functions to evaluate, one per field
247 . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
248 - mode    - The insertion mode for values
249 
250   Output Parameter:
251 . X - vector
252 
253   Level: developer
254 
255 .seealso: DMPlexComputeL2Diff()
256 @*/
257 PetscErrorCode DMPlexProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
258 {
259   Vec            localX;
260   PetscErrorCode ierr;
261 
262   PetscFunctionBegin;
263   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
264   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
265   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr);
266   ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
267   ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
268   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
269   PetscFunctionReturn(0);
270 }
271 
272 #undef __FUNCT__
273 #define __FUNCT__ "DMPlexComputeL2Diff"
274 /*@C
275   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
276 
277   Input Parameters:
278 + dm    - The DM
279 . fe    - The PetscFE object for each field
280 . funcs - The functions to evaluate for each field component
281 . ctxs  - Optional array of contexts to pass to each function, or NULL.
282 - X     - The coefficient vector u_h
283 
284   Output Parameter:
285 . diff - The diff ||u - u_h||_2
286 
287   Level: developer
288 
289 .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff()
290 @*/
291 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
292 {
293   const PetscInt  debug = 0;
294   PetscSection    section;
295   PetscQuadrature quad;
296   Vec             localX;
297   PetscScalar    *funcVal;
298   PetscReal      *coords, *v0, *J, *invJ, detJ;
299   PetscReal       localDiff = 0.0;
300   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
301   PetscErrorCode  ierr;
302 
303   PetscFunctionBegin;
304   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
305   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
306   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
307   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
308   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
309   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
310   for (field = 0; field < numFields; ++field) {
311     PetscInt Nc;
312 
313     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
314     numComponents += Nc;
315   }
316   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
317   ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
318   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
319   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
320   for (c = cStart; c < cEnd; ++c) {
321     PetscScalar *x = NULL;
322     PetscReal    elemDiff = 0.0;
323 
324     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
325     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
326     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
327 
328     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
329       void * const     ctx = ctxs ? ctxs[field] : NULL;
330       const PetscReal *quadPoints, *quadWeights;
331       PetscReal       *basis;
332       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;
333 
334       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
335       ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr);
336       ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr);
337       ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr);
338       if (debug) {
339         char title[1024];
340         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
341         ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
342       }
343       for (q = 0; q < numQuadPoints; ++q) {
344         for (d = 0; d < dim; d++) {
345           coords[d] = v0[d];
346           for (e = 0; e < dim; e++) {
347             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
348           }
349         }
350         (*funcs[field])(coords, funcVal, ctx);
351         for (fc = 0; fc < numBasisComps; ++fc) {
352           PetscScalar interpolant = 0.0;
353 
354           for (f = 0; f < numBasisFuncs; ++f) {
355             const PetscInt fidx = f*numBasisComps+fc;
356             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
357           }
358           if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);}
359           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
360         }
361       }
362       comp        += numBasisComps;
363       fieldOffset += numBasisFuncs*numBasisComps;
364     }
365     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
366     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
367     localDiff += elemDiff;
368   }
369   ierr  = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr);
370   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
371   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
372   *diff = PetscSqrtReal(*diff);
373   PetscFunctionReturn(0);
374 }
375 
376 #undef __FUNCT__
377 #define __FUNCT__ "DMPlexComputeL2GradientDiff"
378 /*@C
379   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
380 
381   Input Parameters:
382 + dm    - The DM
383 . fe    - The PetscFE object for each field
384 . funcs - The gradient functions to evaluate for each field component
385 . ctxs  - Optional array of contexts to pass to each function, or NULL.
386 . X     - The coefficient vector u_h
387 - n     - The vector to project along
388 
389   Output Parameter:
390 . diff - The diff ||(grad u - grad u_h) . n||_2
391 
392   Level: developer
393 
394 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
395 @*/
396 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
397 {
398   const PetscInt  debug = 0;
399   PetscSection    section;
400   PetscQuadrature quad;
401   Vec             localX;
402   PetscScalar    *funcVal, *interpolantVec;
403   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
404   PetscReal       localDiff = 0.0;
405   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
406   PetscErrorCode  ierr;
407 
408   PetscFunctionBegin;
409   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
410   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
411   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
412   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
413   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
414   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
415   for (field = 0; field < numFields; ++field) {
416     PetscInt Nc;
417 
418     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
419     numComponents += Nc;
420   }
421   /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */
422   ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr);
423   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
424   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
425   for (c = cStart; c < cEnd; ++c) {
426     PetscScalar *x = NULL;
427     PetscReal    elemDiff = 0.0;
428 
429     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
430     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
431     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
432 
433     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
434       void * const     ctx = ctxs ? ctxs[field] : NULL;
435       const PetscReal *quadPoints, *quadWeights;
436       PetscReal       *basisDer;
437       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
438 
439       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
440       ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
441       ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr);
442       ierr = PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);CHKERRQ(ierr);
443       if (debug) {
444         char title[1024];
445         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
446         ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr);
447       }
448       for (q = 0; q < numQuadPoints; ++q) {
449         for (d = 0; d < dim; d++) {
450           coords[d] = v0[d];
451           for (e = 0; e < dim; e++) {
452             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
453           }
454         }
455         (*funcs[field])(coords, n, funcVal, ctx);
456         for (fc = 0; fc < Ncomp; ++fc) {
457           PetscScalar interpolant = 0.0;
458 
459           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
460           for (f = 0; f < Nb; ++f) {
461             const PetscInt fidx = f*Ncomp+fc;
462 
463             for (d = 0; d < dim; ++d) {
464               realSpaceDer[d] = 0.0;
465               for (g = 0; g < dim; ++g) {
466                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
467               }
468               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
469             }
470           }
471           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
472           if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);}
473           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
474         }
475       }
476       comp        += Ncomp;
477       fieldOffset += Nb*Ncomp;
478     }
479     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
480     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
481     localDiff += elemDiff;
482   }
483   ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
484   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
485   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
486   *diff = PetscSqrtReal(*diff);
487   PetscFunctionReturn(0);
488 }
489 
490 #undef __FUNCT__
491 #define __FUNCT__ "DMPlexComputeResidualFEM"
492 /*@
493   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
494 
495   Input Parameters:
496 + dm - The mesh
497 . X  - Local input vector
498 - user - The user context
499 
500   Output Parameter:
501 . F  - Local output vector
502 
503   Note:
504   The first member of the user context must be an FEMContext.
505 
506   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
507   like a GPU, or vectorize on a multicore machine.
508 
509   Level: developer
510 
511 .seealso: DMPlexComputeJacobianActionFEM()
512 @*/
513 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
514 {
515   DM_Plex          *mesh  = (DM_Plex *) dm->data;
516   PetscFEM         *fem   = (PetscFEM *) user;
517   PetscFE          *fe    = fem->fe;
518   PetscFE          *feAux = fem->feAux;
519   PetscFE          *feBd  = fem->feBd;
520   const char       *name  = "Residual";
521   DM                dmAux;
522   Vec               A;
523   PetscQuadrature   q;
524   PetscCellGeometry geom;
525   PetscSection      section, sectionAux;
526   PetscReal        *v0, *J, *invJ, *detJ;
527   PetscScalar      *elemVec, *u, *a = NULL;
528   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
529   PetscInt          cellDof = 0, numComponents = 0;
530   PetscInt          cellDofAux = 0, numComponentsAux = 0;
531   PetscErrorCode    ierr;
532 
533   PetscFunctionBegin;
534   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
535   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
536   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
537   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
538   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
539   numCells = cEnd - cStart;
540   for (f = 0; f < Nf; ++f) {
541     PetscInt Nb, Nc;
542 
543     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
544     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
545     cellDof       += Nb*Nc;
546     numComponents += Nc;
547   }
548   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
549   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
550   if (dmAux) {
551     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
552     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
553   }
554   for (f = 0; f < NfAux; ++f) {
555     PetscInt Nb, Nc;
556 
557     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
558     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
559     cellDofAux       += Nb*Nc;
560     numComponentsAux += Nc;
561   }
562   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
563   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
564   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
565   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
566   for (c = cStart; c < cEnd; ++c) {
567     PetscScalar *x = NULL;
568     PetscInt     i;
569 
570     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
571     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
572     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
573     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
574     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
575     if (dmAux) {
576       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
577       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
578       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
579     }
580   }
581   for (f = 0; f < Nf; ++f) {
582     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f];
583     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f];
584     PetscInt numQuadPoints, Nb;
585     /* Conforming batches */
586     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
587     /* Remainder */
588     PetscInt Nr, offset;
589 
590     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
591     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
592     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
593     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
594     blockSize = Nb*numQuadPoints;
595     batchSize = numBlocks * blockSize;
596     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
597     numChunks = numCells / (numBatches*batchSize);
598     Ne        = numChunks*numBatches*batchSize;
599     Nr        = numCells % (numBatches*batchSize);
600     offset    = numCells - Nr;
601     geom.v0   = v0;
602     geom.J    = J;
603     geom.invJ = invJ;
604     geom.detJ = detJ;
605     ierr = PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
606     geom.v0   = &v0[offset*dim];
607     geom.J    = &J[offset*dim*dim];
608     geom.invJ = &invJ[offset*dim*dim];
609     geom.detJ = &detJ[offset];
610     ierr = PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
611   }
612   for (c = cStart; c < cEnd; ++c) {
613     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
614     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
615   }
616   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
617   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
618   if (feBd) {
619     DMLabel         label, depth;
620     IS              pointIS;
621     const PetscInt *points;
622     PetscInt        dep, numPoints, p, numFaces;
623     PetscReal      *n;
624 
625     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
626     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
627     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
628     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
629     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
630     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
631       PetscInt Nb, Nc;
632 
633       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
634       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
635       cellDof       += Nb*Nc;
636       numComponents += Nc;
637     }
638     for (p = 0, numFaces = 0; p < numPoints; ++p) {
639       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
640       if (dep == dim-1) ++numFaces;
641     }
642     ierr = PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);CHKERRQ(ierr);
643     for (p = 0, f = 0; p < numPoints; ++p) {
644       const PetscInt point = points[p];
645       PetscScalar   *x     = NULL;
646       PetscInt       i;
647 
648       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
649       if (dep != dim-1) continue;
650       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
651       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
652       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
653       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
654       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
655       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
656       ++f;
657     }
658     for (f = 0; f < Nf; ++f) {
659       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f];
660       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f];
661       PetscInt numQuadPoints, Nb;
662       /* Conforming batches */
663       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
664       /* Remainder */
665       PetscInt Nr, offset;
666 
667       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
668       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
669       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
670       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
671       blockSize = Nb*numQuadPoints;
672       batchSize = numBlocks * blockSize;
673       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
674       numChunks = numFaces / (numBatches*batchSize);
675       Ne        = numChunks*numBatches*batchSize;
676       Nr        = numFaces % (numBatches*batchSize);
677       offset    = numFaces - Nr;
678       geom.v0   = v0;
679       geom.n    = n;
680       geom.J    = J;
681       geom.invJ = invJ;
682       geom.detJ = detJ;
683       ierr = PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
684       geom.v0   = &v0[offset*dim];
685       geom.n    = &n[offset*dim];
686       geom.J    = &J[offset*dim*dim];
687       geom.invJ = &invJ[offset*dim*dim];
688       geom.detJ = &detJ[offset];
689       ierr = PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
690     }
691     for (p = 0, f = 0; p < numPoints; ++p) {
692       const PetscInt point = points[p];
693 
694       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
695       if (dep != dim-1) continue;
696       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
697       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
698       ++f;
699     }
700     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
701     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
702     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
703   }
704   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
705   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
706   PetscFunctionReturn(0);
707 }
708 
709 #undef __FUNCT__
710 #define __FUNCT__ "DMPlexComputeIFunctionFEM"
711 /*@
712   DMPlexComputeIFunctionFEM - Form the local implicit function F from the local input X, X_t using pointwise functions specified by the user
713 
714   Input Parameters:
715 + dm - The mesh
716 . X  - Local input vector
717 . X_t  - Time derivative of the local input vector
718 - user - The user context
719 
720   Output Parameter:
721 . F  - Local output vector
722 
723   Note:
724   The first member of the user context must be an FEMContext.
725 
726   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
727   like a GPU, or vectorize on a multicore machine.
728 
729   Level: developer
730 
731 .seealso: DMPlexComputeResidualFEM()
732 @*/
733 PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, Vec X, Vec X_t, Vec F, void *user)
734 {
735   DM_Plex          *mesh  = (DM_Plex *) dm->data;
736   PetscFEM         *fem   = (PetscFEM *) user;
737   PetscFE          *fe    = fem->fe;
738   PetscFE          *feAux = fem->feAux;
739   PetscFE          *feBd  = fem->feBd;
740   const char       *name  = "Residual";
741   DM                dmAux;
742   Vec               A;
743   PetscQuadrature   q;
744   PetscCellGeometry geom;
745   PetscSection      section, sectionAux;
746   PetscReal        *v0, *J, *invJ, *detJ;
747   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
748   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
749   PetscInt          cellDof = 0, numComponents = 0;
750   PetscInt          cellDofAux = 0, numComponentsAux = 0;
751   PetscErrorCode    ierr;
752 
753   PetscFunctionBegin;
754   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
755   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
756   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
757   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
758   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
759   numCells = cEnd - cStart;
760   for (f = 0; f < Nf; ++f) {
761     PetscInt Nb, Nc;
762 
763     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
764     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
765     cellDof       += Nb*Nc;
766     numComponents += Nc;
767   }
768   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
769   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
770   if (dmAux) {
771     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
772     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
773   }
774   for (f = 0; f < NfAux; ++f) {
775     PetscInt Nb, Nc;
776 
777     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
778     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
779     cellDofAux       += Nb*Nc;
780     numComponentsAux += Nc;
781   }
782   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
783   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
784   ierr = PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
785   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
786   for (c = cStart; c < cEnd; ++c) {
787     PetscScalar *x = NULL, *x_t = NULL;
788     PetscInt     i;
789 
790     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
791     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
792     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
793     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
794     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
795     ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
796     for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i];
797     ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
798     if (dmAux) {
799       PetscScalar *x_a = NULL;
800       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
801       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i];
802       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
803     }
804   }
805   for (f = 0; f < Nf; ++f) {
806     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f];
807     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f];
808     PetscInt numQuadPoints, Nb;
809     /* Conforming batches */
810     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
811     /* Remainder */
812     PetscInt Nr, offset;
813 
814     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
815     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
816     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
817     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
818     blockSize = Nb*numQuadPoints;
819     batchSize = numBlocks * blockSize;
820     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
821     numChunks = numCells / (numBatches*batchSize);
822     Ne        = numChunks*numBatches*batchSize;
823     Nr        = numCells % (numBatches*batchSize);
824     offset    = numCells - Nr;
825     geom.v0   = v0;
826     geom.J    = J;
827     geom.invJ = invJ;
828     geom.detJ = detJ;
829     ierr = PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
830     geom.v0   = &v0[offset*dim];
831     geom.J    = &J[offset*dim*dim];
832     geom.invJ = &invJ[offset*dim*dim];
833     geom.detJ = &detJ[offset];
834     ierr = PetscFEIntegrateIFunction(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
835   }
836   for (c = cStart; c < cEnd; ++c) {
837     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
838     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
839   }
840   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
841   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
842   if (feBd) {
843     DMLabel         label, depth;
844     IS              pointIS;
845     const PetscInt *points;
846     PetscInt        dep, numPoints, p, numFaces;
847     PetscReal      *n;
848 
849     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
850     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
851     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
852     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
853     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
854     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
855       PetscInt Nb, Nc;
856 
857       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
858       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
859       cellDof       += Nb*Nc;
860       numComponents += Nc;
861     }
862     for (p = 0, numFaces = 0; p < numPoints; ++p) {
863       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
864       if (dep == dim-1) ++numFaces;
865     }
866     ierr = PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);CHKERRQ(ierr);
867     ierr = PetscMalloc1(numFaces*cellDof,&u_t);CHKERRQ(ierr);
868     for (p = 0, f = 0; p < numPoints; ++p) {
869       const PetscInt point = points[p];
870       PetscScalar   *x     = NULL;
871       PetscInt       i;
872 
873       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
874       if (dep != dim-1) continue;
875       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
876       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
877       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
878       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
879       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
880       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
881       ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
882       for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i];
883       ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
884       ++f;
885     }
886     for (f = 0; f < Nf; ++f) {
887       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f];
888       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f];
889       PetscInt numQuadPoints, Nb;
890       /* Conforming batches */
891       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
892       /* Remainder */
893       PetscInt Nr, offset;
894 
895       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
896       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
897       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
898       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
899       blockSize = Nb*numQuadPoints;
900       batchSize = numBlocks * blockSize;
901       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
902       numChunks = numFaces / (numBatches*batchSize);
903       Ne        = numChunks*numBatches*batchSize;
904       Nr        = numFaces % (numBatches*batchSize);
905       offset    = numFaces - Nr;
906       geom.v0   = v0;
907       geom.n    = n;
908       geom.J    = J;
909       geom.invJ = invJ;
910       geom.detJ = detJ;
911       ierr = PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
912       geom.v0   = &v0[offset*dim];
913       geom.n    = &n[offset*dim];
914       geom.J    = &J[offset*dim*dim];
915       geom.invJ = &invJ[offset*dim*dim];
916       geom.detJ = &detJ[offset];
917       ierr = PetscFEIntegrateBdIFunction(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
918     }
919     for (p = 0, f = 0; p < numPoints; ++p) {
920       const PetscInt point = points[p];
921 
922       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
923       if (dep != dim-1) continue;
924       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
925       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
926       ++f;
927     }
928     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
929     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
930     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
931     ierr = PetscFree(u_t);CHKERRQ(ierr);
932   }
933   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
934   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
935   PetscFunctionReturn(0);
936 }
937 
938 #undef __FUNCT__
939 #define __FUNCT__ "DMPlexComputeJacobianActionFEM"
940 /*@C
941   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
942 
943   Input Parameters:
944 + dm - The mesh
945 . J  - The Jacobian shell matrix
946 . X  - Local input vector
947 - user - The user context
948 
949   Output Parameter:
950 . F  - Local output vector
951 
952   Note:
953   The first member of the user context must be an FEMContext.
954 
955   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
956   like a GPU, or vectorize on a multicore machine.
957 
958   Level: developer
959 
960 .seealso: DMPlexComputeResidualFEM()
961 @*/
962 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
963 {
964   DM_Plex          *mesh = (DM_Plex *) dm->data;
965   PetscFEM         *fem  = (PetscFEM *) user;
966   PetscFE          *fe   = fem->fe;
967   PetscQuadrature   quad;
968   PetscCellGeometry geom;
969   PetscSection      section;
970   JacActionCtx     *jctx;
971   PetscReal        *v0, *J, *invJ, *detJ;
972   PetscScalar      *elemVec, *u, *a;
973   PetscInt          dim, numFields, field, numCells, cStart, cEnd, c;
974   PetscInt          cellDof = 0;
975   PetscErrorCode    ierr;
976 
977   PetscFunctionBegin;
978   /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
979   ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
980   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
981   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
982   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
983   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
984   numCells = cEnd - cStart;
985   for (field = 0; field < numFields; ++field) {
986     PetscInt Nb, Nc;
987 
988     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
989     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
990     cellDof += Nb*Nc;
991   }
992   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
993   ierr = PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&a,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
994   for (c = cStart; c < cEnd; ++c) {
995     PetscScalar *x = NULL;
996     PetscInt     i;
997 
998     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
999     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1000     ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1001     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1002     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1003     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1004     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
1005     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1006   }
1007   for (field = 0; field < numFields; ++field) {
1008     PetscInt numQuadPoints, Nb;
1009     /* Conforming batches */
1010     PetscInt numBlocks  = 1;
1011     PetscInt numBatches = 1;
1012     PetscInt numChunks, Ne, blockSize, batchSize;
1013     /* Remainder */
1014     PetscInt Nr, offset;
1015 
1016     ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr);
1017     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
1018     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1019     blockSize = Nb*numQuadPoints;
1020     batchSize = numBlocks * blockSize;
1021     numChunks = numCells / (numBatches*batchSize);
1022     Ne        = numChunks*numBatches*batchSize;
1023     Nr        = numCells % (numBatches*batchSize);
1024     offset    = numCells - Nr;
1025     geom.v0   = v0;
1026     geom.J    = J;
1027     geom.invJ = invJ;
1028     geom.detJ = detJ;
1029     ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
1030     geom.v0   = &v0[offset*dim];
1031     geom.J    = &J[offset*dim*dim];
1032     geom.invJ = &invJ[offset*dim*dim];
1033     geom.detJ = &detJ[offset];
1034     ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof],
1035                                           fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr);
1036   }
1037   for (c = cStart; c < cEnd; ++c) {
1038     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
1039     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
1040   }
1041   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1042   if (mesh->printFEM) {
1043     PetscMPIInt rank, numProcs;
1044     PetscInt    p;
1045 
1046     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
1047     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
1048     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr);
1049     for (p = 0; p < numProcs; ++p) {
1050       if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
1051       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
1052     }
1053   }
1054   /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
1055   PetscFunctionReturn(0);
1056 }
1057 
1058 #undef __FUNCT__
1059 #define __FUNCT__ "DMPlexComputeJacobianFEM"
1060 /*@
1061   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1062 
1063   Input Parameters:
1064 + dm - The mesh
1065 . X  - Local input vector
1066 - user - The user context
1067 
1068   Output Parameter:
1069 . Jac  - Jacobian matrix
1070 
1071   Note:
1072   The first member of the user context must be an FEMContext.
1073 
1074   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1075   like a GPU, or vectorize on a multicore machine.
1076 
1077   Level: developer
1078 
1079 .seealso: FormFunctionLocal()
1080 @*/
1081 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1082 {
1083   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1084   PetscFEM         *fem   = (PetscFEM *) user;
1085   PetscFE          *fe    = fem->fe;
1086   PetscFE          *feAux = fem->feAux;
1087   const char       *name  = "Jacobian";
1088   DM                dmAux;
1089   Vec               A;
1090   PetscQuadrature   quad;
1091   PetscCellGeometry geom;
1092   PetscSection      section, globalSection, sectionAux;
1093   PetscReal        *v0, *J, *invJ, *detJ;
1094   PetscScalar      *elemMat, *u, *a;
1095   PetscInt          dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1096   PetscInt          cellDof = 0, numComponents = 0;
1097   PetscInt          cellDofAux = 0, numComponentsAux = 0;
1098   PetscBool         isShell;
1099   PetscErrorCode    ierr;
1100 
1101   PetscFunctionBegin;
1102   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1103   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
1104   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1105   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1106   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1107   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1108   numCells = cEnd - cStart;
1109   for (f = 0; f < Nf; ++f) {
1110     PetscInt Nb, Nc;
1111 
1112     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
1113     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1114     cellDof       += Nb*Nc;
1115     numComponents += Nc;
1116   }
1117   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1118   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1119   if (dmAux) {
1120     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1121     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
1122   }
1123   for (f = 0; f < NfAux; ++f) {
1124     PetscInt Nb, Nc;
1125 
1126     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
1127     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
1128     cellDofAux       += Nb*Nc;
1129     numComponentsAux += Nc;
1130   }
1131   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
1132   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1133   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr);
1134   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
1135   for (c = cStart; c < cEnd; ++c) {
1136     PetscScalar *x = NULL;
1137     PetscInt     i;
1138 
1139     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1140     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1141     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1142     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1143     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1144     if (dmAux) {
1145       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1146       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
1147       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1148     }
1149   }
1150   ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1151   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1152     PetscInt numQuadPoints, Nb;
1153     /* Conforming batches */
1154     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1155     /* Remainder */
1156     PetscInt Nr, offset;
1157 
1158     ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr);
1159     ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr);
1160     ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1161     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1162     blockSize = Nb*numQuadPoints;
1163     batchSize = numBlocks * blockSize;
1164     ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1165     numChunks = numCells / (numBatches*batchSize);
1166     Ne        = numChunks*numBatches*batchSize;
1167     Nr        = numCells % (numBatches*batchSize);
1168     offset    = numCells - Nr;
1169     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1170       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ];
1171       void   (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ];
1172       void   (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ];
1173       void   (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ];
1174 
1175       geom.v0   = v0;
1176       geom.J    = J;
1177       geom.invJ = invJ;
1178       geom.detJ = detJ;
1179       ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
1180       geom.v0   = &v0[offset*dim];
1181       geom.J    = &J[offset*dim*dim];
1182       geom.invJ = &invJ[offset*dim*dim];
1183       geom.detJ = &detJ[offset];
1184       ierr = PetscFEIntegrateJacobian(fe[fieldI], Nr, Nf, fe, fieldI, fieldJ, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr);
1185     }
1186   }
1187   for (c = cStart; c < cEnd; ++c) {
1188     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
1189     ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
1190   }
1191   ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1192   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1193   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1194   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1195   if (mesh->printFEM) {
1196     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1197     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
1198     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1199   }
1200   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1201   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
1202   if (isShell) {
1203     JacActionCtx *jctx;
1204 
1205     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1206     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
1207   }
1208   PetscFunctionReturn(0);
1209 }
1210 
1211 #if 0
1212 
1213 static void g0_identity_1d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1214 {
1215   g3[0] = 1.0;
1216 }
1217 
1218 static void g0_identity_2d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1219 {
1220   g3[0] = g3[3] = 1.0;
1221 }
1222 
1223 static void g0_identity_3d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1224 {
1225   g3[0] = g3[4] = g3[8] = 1.0;
1226 }
1227 
1228 #undef __FUNCT__
1229 #define __FUNCT__ "DMPlexComputeInterpolatorFEMBroken"
1230 PetscErrorCode DMPlexComputeInterpolatorFEMBroken(DM dmc, DM dmf, Mat I, void *user)
1231 {
1232   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1233   PetscFEM         *fem   = (PetscFEM *) user;
1234   PetscFE          *fe    = fem->fe;
1235   const char       *name  = "Interpolator";
1236   PetscFE          *feRef;
1237   PetscQuadrature   quad, quadOld;
1238   PetscCellGeometry geom;
1239   PetscSection      fsection, fglobalSection;
1240   PetscSection      csection, cglobalSection;
1241   PetscReal        *v0, *J, *invJ, *detJ;
1242   PetscScalar      *elemMat;
1243   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1244   PetscInt          rCellDof = 0, cCellDof = 0, numComponents = 0;
1245   PetscErrorCode    ierr;
1246 
1247   PetscFunctionBegin;
1248 #if 0
1249   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1250 #endif
1251   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1252   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1253   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1254   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1255   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1256   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1257   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1258   numCells = cEnd - cStart;
1259   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1260   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1261     ierr = PetscFERefine(fe[fieldI], &feRef[fieldI]);CHKERRQ(ierr);
1262   }
1263   for (f = 0; f < Nf; ++f) {
1264     PetscInt rNb, cNb, Nc;
1265 
1266     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1267     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1268     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1269     numComponents += Nc;
1270     rCellDof += rNb*Nc;
1271     cCellDof += cNb*Nc;
1272   }
1273   ierr = MatZeroEntries(I);CHKERRQ(ierr);
1274   ierr = PetscMalloc5(numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1275   for (c = cStart; c < cEnd; ++c) {
1276     ierr = DMPlexComputeCellGeometry(dmc, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1277     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1278   }
1279   ierr = PetscMemzero(elemMat, numCells*rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1280   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1281     PetscInt numQuadPoints, Nb;
1282     /* Conforming batches */
1283     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1284     /* Remainder */
1285     PetscInt Nr, offset;
1286 
1287     /* Make new fine FE which refines the ref cell and the quadrature rule */
1288     ierr = PetscFEGetQuadrature(feRef[fieldI], &quad);CHKERRQ(ierr);
1289     ierr = PetscFEGetDimension(feRef[fieldI], &Nb);CHKERRQ(ierr);
1290     ierr = PetscFEGetTileSizes(feRef[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1291     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1292     blockSize = Nb*numQuadPoints;
1293     batchSize = numBlocks * blockSize;
1294     ierr = PetscFESetTileSizes(feRef[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1295     numChunks = numCells / (numBatches*batchSize);
1296     Ne        = numChunks*numBatches*batchSize;
1297     Nr        = numCells % (numBatches*batchSize);
1298     offset    = numCells - Nr;
1299 
1300     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1301       /* void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; */
1302       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = g0_identity_2d_static;
1303 
1304       /* Replace quadrature in coarse FE with refined quadrature */
1305       ierr = PetscFEGetQuadrature(fe[fieldJ], &quadOld);CHKERRQ(ierr);
1306       ierr = PetscObjectReference((PetscObject) quadOld);CHKERRQ(ierr);
1307       ierr = PetscFESetQuadrature(fe[fieldJ], quad);CHKERRQ(ierr);
1308       geom.v0   = v0;
1309       geom.J    = J;
1310       geom.invJ = invJ;
1311       geom.detJ = detJ;
1312       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Ne, Nf, feRef, fieldI, fe, fieldJ, geom, g0, elemMat);CHKERRQ(ierr);
1313       geom.v0   = &v0[offset*dim];
1314       geom.J    = &J[offset*dim*dim];
1315       geom.invJ = &invJ[offset*dim*dim];
1316       geom.detJ = &detJ[offset];
1317       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Nr, Nf, feRef, fieldI, fe, fieldJ, geom, g0, &elemMat[offset*rCellDof*cCellDof]);CHKERRQ(ierr);
1318       ierr = PetscFESetQuadrature(fe[fieldJ], quadOld);CHKERRQ(ierr);
1319     }
1320   }
1321   for (c = cStart; c < cEnd; ++c) {
1322     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, rCellDof, cCellDof, &elemMat[c*rCellDof*cCellDof]);CHKERRQ(ierr);}
1323     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, I, c, &elemMat[c*rCellDof*cCellDof], ADD_VALUES);CHKERRQ(ierr);
1324   }
1325   ierr = PetscFree5(v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1326   ierr = PetscFree(feRef);CHKERRQ(ierr);
1327   ierr = MatAssemblyBegin(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1328   ierr = MatAssemblyEnd(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1329   if (mesh->printFEM) {
1330     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1331     ierr = MatChop(I, 1.0e-10);CHKERRQ(ierr);
1332     ierr = MatView(I, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1333   }
1334 #if 0
1335   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1336 #endif
1337   PetscFunctionReturn(0);
1338 }
1339 #endif
1340 
1341 #undef __FUNCT__
1342 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1343 /*@
1344   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1345 
1346   Input Parameters:
1347 + dmf  - The fine mesh
1348 . dmc  - The coarse mesh
1349 - user - The user context
1350 
1351   Output Parameter:
1352 . In  - The interpolation matrix
1353 
1354   Note:
1355   The first member of the user context must be an FEMContext.
1356 
1357   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1358   like a GPU, or vectorize on a multicore machine.
1359 
1360   Level: developer
1361 
1362 .seealso: DMPlexComputeJacobianFEM()
1363 @*/
1364 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1365 {
1366   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1367   PetscFEM         *fem   = (PetscFEM *) user;
1368   PetscFE          *fe    = fem->fe;
1369   const char       *name  = "Interpolator";
1370   PetscFE          *feRef;
1371   PetscSection      fsection, fglobalSection;
1372   PetscSection      csection, cglobalSection;
1373   PetscScalar      *elemMat;
1374   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
1375   PetscInt          rCellDof = 0, cCellDof = 0;
1376   PetscErrorCode    ierr;
1377 
1378   PetscFunctionBegin;
1379 #if 0
1380   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1381 #endif
1382   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1383   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1384   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1385   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1386   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1387   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1388   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1389   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1390   for (f = 0; f < Nf; ++f) {
1391     PetscInt rNb, cNb, Nc;
1392 
1393     ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr);
1394     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1395     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1396     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1397     rCellDof += rNb*Nc;
1398     cCellDof += cNb*Nc;
1399   }
1400   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1401   ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1402   ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1403   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1404     PetscDualSpace   Qref;
1405     PetscQuadrature  f;
1406     const PetscReal *qpoints, *qweights;
1407     PetscReal       *points;
1408     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1409 
1410     /* Compose points from all dual basis functionals */
1411     ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr);
1412     ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1413     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1414     for (i = 0; i < fpdim; ++i) {
1415       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1416       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1417       npoints += Np;
1418     }
1419     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1420     for (i = 0, k = 0; i < fpdim; ++i) {
1421       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1422       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1423       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1424     }
1425 
1426     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1427       PetscReal *B;
1428       PetscInt   NcJ, cpdim, j;
1429 
1430       /* Evaluate basis at points */
1431       ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr);
1432       if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ);
1433       ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr);
1434       /* For now, fields only interpolate themselves */
1435       if (fieldI == fieldJ) {
1436         ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1437         for (i = 0, k = 0; i < fpdim; ++i) {
1438           ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1439           ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1440           for (p = 0; p < Np; ++p, ++k) {
1441             for (j = 0; j < cpdim; ++j) {
1442               for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cCellDof + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p];
1443             }
1444           }
1445         }
1446         ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1447       }
1448       offsetJ += cpdim*NcJ;
1449     }
1450     offsetI += fpdim*Nc;
1451     ierr = PetscFree(points);CHKERRQ(ierr);
1452   }
1453   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);}
1454   for (c = cStart; c < cEnd; ++c) {
1455     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1456   }
1457   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1458   ierr = PetscFree(feRef);CHKERRQ(ierr);
1459   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1460   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1461   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1462   if (mesh->printFEM) {
1463     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1464     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1465     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1466   }
1467 #if 0
1468   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1469 #endif
1470   PetscFunctionReturn(0);
1471 }
1472 
1473 #undef __FUNCT__
1474 #define __FUNCT__ "DMPlexAddBoundary"
1475 /* The ids can be overridden by the command line option -bc_<boundary name> */
1476 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
1477 {
1478   DM_Plex       *mesh = (DM_Plex *) dm->data;
1479   DMBoundary     b;
1480   PetscErrorCode ierr;
1481 
1482   PetscFunctionBegin;
1483   ierr = PetscNew(&b);CHKERRQ(ierr);
1484   ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr);
1485   ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr);
1486   ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr);
1487   b->essential   = isEssential;
1488   b->field       = field;
1489   b->func        = bcFunc;
1490   b->numids      = numids;
1491   b->ctx         = ctx;
1492   b->next        = mesh->boundary;
1493   mesh->boundary = b;
1494   PetscFunctionReturn(0);
1495 }
1496 
1497 #undef __FUNCT__
1498 #define __FUNCT__ "DMPlexGetNumBoundary"
1499 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd)
1500 {
1501   DM_Plex   *mesh = (DM_Plex *) dm->data;
1502   DMBoundary b    = mesh->boundary;
1503 
1504   PetscFunctionBegin;
1505   *numBd = 0;
1506   while (b) {++(*numBd); b = b->next;}
1507   PetscFunctionReturn(0);
1508 }
1509 
1510 #undef __FUNCT__
1511 #define __FUNCT__ "DMPlexGetBoundary"
1512 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
1513 {
1514   DM_Plex   *mesh = (DM_Plex *) dm->data;
1515   DMBoundary b    = mesh->boundary;
1516   PetscInt   n    = 0;
1517 
1518   PetscFunctionBegin;
1519   while (b) {
1520     if (n == bd) break;
1521     b = b->next;
1522     ++n;
1523   }
1524   if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
1525   if (isEssential) {
1526     PetscValidPointer(isEssential, 3);
1527     *isEssential = b->essential;
1528   }
1529   if (name) {
1530     PetscValidPointer(name, 4);
1531     *name = b->name;
1532   }
1533   if (field) {
1534     PetscValidPointer(field, 5);
1535     *field = b->field;
1536   }
1537   if (func) {
1538     PetscValidPointer(func, 6);
1539     *func = b->func;
1540   }
1541   if (numids) {
1542     PetscValidPointer(numids, 7);
1543     *numids = b->numids;
1544   }
1545   if (ids) {
1546     PetscValidPointer(ids, 8);
1547     *ids = b->ids;
1548   }
1549   if (ctx) {
1550     PetscValidPointer(ctx, 9);
1551     *ctx = b->ctx;
1552   }
1553   PetscFunctionReturn(0);
1554 }
1555