Actual source code: plexfem.c

petsc-3.4.4 2014-03-13
  1: #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/

  5: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
  6: {
  7:   DM_Plex *mesh = (DM_Plex*) dm->data;

 12:   *scale = mesh->scale[unit];
 13:   return(0);
 14: }

 18: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
 19: {
 20:   DM_Plex *mesh = (DM_Plex*) dm->data;

 24:   mesh->scale[unit] = scale;
 25:   return(0);
 26: }

 28: PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
 29: {
 30:   switch (i) {
 31:   case 0:
 32:     switch (j) {
 33:     case 0: return 0;
 34:     case 1:
 35:       switch (k) {
 36:       case 0: return 0;
 37:       case 1: return 0;
 38:       case 2: return 1;
 39:       }
 40:     case 2:
 41:       switch (k) {
 42:       case 0: return 0;
 43:       case 1: return -1;
 44:       case 2: return 0;
 45:       }
 46:     }
 47:   case 1:
 48:     switch (j) {
 49:     case 0:
 50:       switch (k) {
 51:       case 0: return 0;
 52:       case 1: return 0;
 53:       case 2: return -1;
 54:       }
 55:     case 1: return 0;
 56:     case 2:
 57:       switch (k) {
 58:       case 0: return 1;
 59:       case 1: return 0;
 60:       case 2: return 0;
 61:       }
 62:     }
 63:   case 2:
 64:     switch (j) {
 65:     case 0:
 66:       switch (k) {
 67:       case 0: return 0;
 68:       case 1: return 1;
 69:       case 2: return 0;
 70:       }
 71:     case 1:
 72:       switch (k) {
 73:       case 0: return -1;
 74:       case 1: return 0;
 75:       case 2: return 0;
 76:       }
 77:     case 2: return 0;
 78:     }
 79:   }
 80:   return 0;
 81: }

 85: /*@C
 86:   DMPlexCreateRigidBody - create rigid body modes from coordinates

 88:   Collective on DM

 90:   Input Arguments:
 91: + dm - the DM
 92: . section - the local section associated with the rigid field, or NULL for the default section
 93: - globalSection - the global section associated with the rigid field, or NULL for the default section

 95:   Output Argument:
 96: . sp - the null space

 98:   Note: This is necessary to take account of Dirichlet conditions on the displacements

100:   Level: advanced

102: .seealso: MatNullSpaceCreate()
103: @*/
104: PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
105: {
106:   MPI_Comm       comm;
107:   Vec            coordinates, localMode, mode[6];
108:   PetscSection   coordSection;
109:   PetscScalar   *coords;
110:   PetscInt       dim, vStart, vEnd, v, n, m, d, i, j;

114:   PetscObjectGetComm((PetscObject)dm,&comm);
115:   DMPlexGetDimension(dm, &dim);
116:   if (dim == 1) {
117:     MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
118:     return(0);
119:   }
120:   if (!section)       {DMGetDefaultSection(dm, &section);}
121:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
122:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
123:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
124:   DMPlexGetCoordinateSection(dm, &coordSection);
125:   DMGetCoordinatesLocal(dm, &coordinates);
126:   m    = (dim*(dim+1))/2;
127:   VecCreate(comm, &mode[0]);
128:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
129:   VecSetUp(mode[0]);
130:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
131:   /* Assume P1 */
132:   DMGetLocalVector(dm, &localMode);
133:   for (d = 0; d < dim; ++d) {
134:     PetscScalar values[3] = {0.0, 0.0, 0.0};

136:     values[d] = 1.0;
137:     VecSet(localMode, 0.0);
138:     for (v = vStart; v < vEnd; ++v) {
139:       DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
140:     }
141:     DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
142:     DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
143:   }
144:   VecGetArray(coordinates, &coords);
145:   for (d = dim; d < dim*(dim+1)/2; ++d) {
146:     PetscInt i, j, k = dim > 2 ? d - dim : d;

148:     VecSet(localMode, 0.0);
149:     for (v = vStart; v < vEnd; ++v) {
150:       PetscScalar values[3] = {0.0, 0.0, 0.0};
151:       PetscInt    off;

153:       PetscSectionGetOffset(coordSection, v, &off);
154:       for (i = 0; i < dim; ++i) {
155:         for (j = 0; j < dim; ++j) {
156:           values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
157:         }
158:       }
159:       DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
160:     }
161:     DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
162:     DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
163:   }
164:   VecRestoreArray(coordinates, &coords);
165:   DMRestoreLocalVector(dm, &localMode);
166:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
167:   /* Orthonormalize system */
168:   for (i = dim; i < m; ++i) {
169:     PetscScalar dots[6];

171:     VecMDot(mode[i], i, mode, dots);
172:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
173:     VecMAXPY(mode[i], i, dots, mode);
174:     VecNormalize(mode[i], NULL);
175:   }
176:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
177:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
178:   return(0);
179: }
180: /*******************************************************************************
181: This should be in a separate Discretization object, but I am not sure how to lay
182: it out yet, so I am stuffing things here while I experiment.
183: *******************************************************************************/
186: PetscErrorCode DMPlexSetFEMIntegration(DM dm,
187:                                           PetscErrorCode (*integrateResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
188:                                                                                  const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
189:                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
190:                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
191:                                           PetscErrorCode (*integrateBdResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
192:                                                                                    const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
193:                                                                                    void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]),
194:                                                                                    void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]), PetscScalar[]),
195:                                           PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[],
196:                                                                                        const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
197:                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
198:                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
199:                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
200:                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
201:                                           PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
202:                                                                                  const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
203:                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
204:                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
205:                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
206:                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]))
207: {
208:   DM_Plex *mesh = (DM_Plex*) dm->data;

212:   mesh->integrateResidualFEM       = integrateResidualFEM;
213:   mesh->integrateBdResidualFEM     = integrateBdResidualFEM;
214:   mesh->integrateJacobianActionFEM = integrateJacobianActionFEM;
215:   mesh->integrateJacobianFEM       = integrateJacobianFEM;
216:   return(0);
217: }

221: PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec localX)
222: {
223:   Vec            coordinates;
224:   PetscSection   section, cSection;
225:   PetscInt       dim, vStart, vEnd, v, c, d;
226:   PetscScalar   *values, *cArray;
227:   PetscReal     *coords;

231:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
232:   DMGetDefaultSection(dm, &section);
233:   DMPlexGetCoordinateSection(dm, &cSection);
234:   DMGetCoordinatesLocal(dm, &coordinates);
235:   PetscMalloc(numComp * sizeof(PetscScalar), &values);
236:   VecGetArray(coordinates, &cArray);
237:   PetscSectionGetDof(cSection, vStart, &dim);
238:   PetscMalloc(dim * sizeof(PetscReal),&coords);
239:   for (v = vStart; v < vEnd; ++v) {
240:     PetscInt dof, off;

242:     PetscSectionGetDof(cSection, v, &dof);
243:     PetscSectionGetOffset(cSection, v, &off);
244:     if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim);
245:     for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]);
246:     for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
247:     VecSetValuesSection(localX, section, v, values, mode);
248:   }
249:   VecRestoreArray(coordinates, &cArray);
250:   /* Temporary, must be replaced by a projection on the finite element basis */
251:   {
252:     PetscInt eStart = 0, eEnd = 0, e, depth;

254:     DMPlexGetLabelSize(dm, "depth", &depth);
255:     --depth;
256:     if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
257:     for (e = eStart; e < eEnd; ++e) {
258:       const PetscInt *cone = NULL;
259:       PetscInt        coneSize, d;
260:       PetscScalar    *coordsA, *coordsB;

262:       DMPlexGetConeSize(dm, e, &coneSize);
263:       DMPlexGetCone(dm, e, &cone);
264:       if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e);
265:       VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);
266:       VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);
267:       for (d = 0; d < dim; ++d) {
268:         coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d]));
269:       }
270:       for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
271:       VecSetValuesSection(localX, section, e, values, mode);
272:     }
273:   }

275:   PetscFree(coords);
276:   PetscFree(values);
277: #if 0
278:   const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin());
279:   PetscReal      detJ;

281:   PetscMalloc(localDof * sizeof(PetscScalar), &values);
282:   PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);
283:   ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true);

285:   for (PetscInt c = cStart; c < cEnd; ++c) {
286:     ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV);
287:     const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
288:     const int                          oSize   = pV.getSize();
289:     int                                v       = 0;

291:     DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);
292:     for (PetscInt cl = 0; cl < oSize; ++cl) {
293:       const PetscInt fDim;

295:       PetscSectionGetDof(oPoints[cl], &fDim);
296:       if (pointDim) {
297:         for (PetscInt d = 0; d < fDim; ++d, ++v) {
298:           values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
299:         }
300:       }
301:     }
302:     DMPlexVecSetClosure(dm, NULL, localX, c, values);
303:     pV.clear();
304:   }
305:   PetscFree2(v0,J);
306:   PetscFree(values);
307: #endif
308:   return(0);
309: }

313: /*@C
314:   DMPlexProjectFunction - This projects the given function into the function space provided.

316:   Input Parameters:
317: + dm      - The DM
318: . numComp - The number of components (functions)
319: . funcs   - The coordinate functions to evaluate
320: - mode    - The insertion mode for values

322:   Output Parameter:
323: . X - vector

325:   Level: developer

327:   Note:
328:   This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector.
329:   We will eventually fix it.

331: .seealso: DMPlexComputeL2Diff()
332: @*/
333: PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec X)
334: {
335:   Vec            localX;

339:   DMGetLocalVector(dm, &localX);
340:   DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);
341:   DMLocalToGlobalBegin(dm, localX, mode, X);
342:   DMLocalToGlobalEnd(dm, localX, mode, X);
343:   DMRestoreLocalVector(dm, &localX);
344:   return(0);
345: }

349: /*@C
350:   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

352:   Input Parameters:
353: + dm    - The DM
354: . quad  - The PetscQuadrature object for each field
355: . funcs - The functions to evaluate for each field component
356: - X     - The coefficient vector u_h

358:   Output Parameter:
359: . diff - The diff ||u - u_h||_2

361:   Level: developer

363: .seealso: DMPlexProjectFunction()
364: @*/
365: PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], void (**funcs)(const PetscReal [], PetscScalar *), Vec X, PetscReal *diff)
366: {
367:   const PetscInt debug = 0;
368:   PetscSection   section;
369:   Vec            localX;
370:   PetscReal     *coords, *v0, *J, *invJ, detJ;
371:   PetscReal      localDiff = 0.0;
372:   PetscInt       dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;

376:   DMPlexGetDimension(dm, &dim);
377:   DMGetDefaultSection(dm, &section);
378:   PetscSectionGetNumFields(section, &numFields);
379:   DMGetLocalVector(dm, &localX);
380:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
381:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
382:   for (field = 0; field < numFields; ++field) {
383:     numComponents += quad[field].numComponents;
384:   }
385:   DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);
386:   PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
387:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
388:   for (c = cStart; c < cEnd; ++c) {
389:     PetscScalar *x;
390:     PetscReal    elemDiff = 0.0;

392:     DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
393:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
394:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

396:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
397:       const PetscInt   numQuadPoints = quad[field].numQuadPoints;
398:       const PetscReal *quadPoints    = quad[field].quadPoints;
399:       const PetscReal *quadWeights   = quad[field].quadWeights;
400:       const PetscInt   numBasisFuncs = quad[field].numBasisFuncs;
401:       const PetscInt   numBasisComps = quad[field].numComponents;
402:       const PetscReal *basis         = quad[field].basis;
403:       PetscInt         q, d, e, fc, f;

405:       if (debug) {
406:         char title[1024];
407:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
408:         DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
409:       }
410:       for (q = 0; q < numQuadPoints; ++q) {
411:         for (d = 0; d < dim; d++) {
412:           coords[d] = v0[d];
413:           for (e = 0; e < dim; e++) {
414:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
415:           }
416:         }
417:         for (fc = 0; fc < numBasisComps; ++fc) {
418:           PetscScalar funcVal;
419:           PetscScalar interpolant = 0.0;

421:           (*funcs[comp+fc])(coords, &funcVal);
422:           for (f = 0; f < numBasisFuncs; ++f) {
423:             const PetscInt fidx = f*numBasisComps+fc;
424:             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
425:           }
426:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal))*quadWeights[q]*detJ);}
427:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal))*quadWeights[q]*detJ;
428:         }
429:       }
430:       comp        += numBasisComps;
431:       fieldOffset += numBasisFuncs*numBasisComps;
432:     }
433:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
434:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
435:     localDiff += elemDiff;
436:   }
437:   PetscFree4(coords,v0,J,invJ);
438:   DMRestoreLocalVector(dm, &localX);
439:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
440:   *diff = PetscSqrtReal(*diff);
441:   return(0);
442: }

446: /*@
447:   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user

449:   Input Parameters:
450: + dm - The mesh
451: . X  - Local input vector
452: - user - The user context

454:   Output Parameter:
455: . F  - Local output vector

457:   Note:
458:   The second member of the user context must be an FEMContext.

460:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
461:   like a GPU, or vectorize on a multicore machine.

463:   Level: developer

465: .seealso: DMPlexComputeJacobianActionFEM()
466: @*/
467: PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
468: {
469:   DM_Plex         *mesh   = (DM_Plex*) dm->data;
470:   PetscFEM        *fem    = (PetscFEM*) &((DM*) user)[1];
471:   PetscQuadrature *quad   = fem->quad;
472:   PetscQuadrature *quadBd = fem->quadBd;
473:   PetscSection     section;
474:   PetscReal       *v0, *n, *J, *invJ, *detJ;
475:   PetscScalar     *elemVec, *u;
476:   PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
477:   PetscInt         cellDof, numComponents;
478:   PetscBool        has;
479:   PetscErrorCode   ierr;

482:   /* PetscLogEventBegin(ResidualFEMEvent,0,0,0,0); */
483:   DMPlexGetDimension(dm, &dim);
484:   DMGetDefaultSection(dm, &section);
485:   PetscSectionGetNumFields(section, &numFields);
486:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
487:   numCells = cEnd - cStart;
488:   for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
489:     cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
490:     numComponents += quad[field].numComponents;
491:   }
492:   DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);
493:   VecSet(F, 0.0);
494:   PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);
495:   for (c = cStart; c < cEnd; ++c) {
496:     PetscScalar *x;
497:     PetscInt     i;

499:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
500:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
501:     DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);

503:     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
504:     DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);
505:   }
506:   for (field = 0; field < numFields; ++field) {
507:     const PetscInt numQuadPoints = quad[field].numQuadPoints;
508:     const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
509:     void           (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[field];
510:     void           (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[field];
511:     /* Conforming batches */
512:     PetscInt blockSize  = numBasisFuncs*numQuadPoints;
513:     PetscInt numBlocks  = 1;
514:     PetscInt batchSize  = numBlocks * blockSize;
515:     PetscInt numBatches = numBatchesTmp;
516:     PetscInt numChunks  = numCells / (numBatches*batchSize);
517:     /* Remainder */
518:     PetscInt numRemainder = numCells % (numBatches * batchSize);
519:     PetscInt offset       = numCells - numRemainder;

521:     (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);
522:     (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
523:                                          f0, f1, &elemVec[offset*cellDof]);
524:   }
525:   for (c = cStart; c < cEnd; ++c) {
526:     if (mesh->printFEM > 1) {DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);}
527:     DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);
528:   }
529:   PetscFree6(u,v0,J,invJ,detJ,elemVec);
530:   /* Integration over the boundary:
531:      - This can probably be generalized to integration over a set of labels, however
532:        the idea here is to do integration where we need the cell normal
533:      - We can replace hardcoding with a registration process, and this is how we hook
534:        up the system to something like FEniCS
535:   */
536:   DMPlexHasLabel(dm, "boundary", &has);
537:   if (has && quadBd) {
538:     DMLabel         label;
539:     IS              pointIS;
540:     const PetscInt *points;
541:     PetscInt        numPoints, p;

543:     DMPlexGetLabel(dm, "boundary", &label);
544:     DMLabelGetStratumSize(label, 1, &numPoints);
545:     DMLabelGetStratumIS(label, 1, &pointIS);
546:     ISGetIndices(pointIS, &points);
547:     for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
548:       cellDof       += quadBd[field].numBasisFuncs*quadBd[field].numComponents;
549:       numComponents += quadBd[field].numComponents;
550:     }
551:     PetscMalloc7(numPoints*cellDof,PetscScalar,&u,numPoints*dim,PetscReal,&v0,numPoints*dim,PetscReal,&n,numPoints*dim*dim,PetscReal,&J,numPoints*dim*dim,PetscReal,&invJ,numPoints,PetscReal,&detJ,numPoints*cellDof,PetscScalar,&elemVec);
552:     for (p = 0; p < numPoints; ++p) {
553:       const PetscInt point = points[p];
554:       PetscScalar   *x;
555:       PetscInt       i;

557:       /* TODO: Add normal determination here */
558:       DMPlexComputeCellGeometry(dm, point, &v0[p*dim], &J[p*dim*dim], &invJ[p*dim*dim], &detJ[p]);
559:       if (detJ[p] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[p], point);
560:       DMPlexVecGetClosure(dm, NULL, X, point, NULL, &x);

562:       for (i = 0; i < cellDof; ++i) u[p*cellDof+i] = x[i];
563:       DMPlexVecRestoreClosure(dm, NULL, X, point, NULL, &x);
564:     }
565:     for (field = 0; field < numFields; ++field) {
566:       const PetscInt numQuadPoints = quadBd[field].numQuadPoints;
567:       const PetscInt numBasisFuncs = quadBd[field].numBasisFuncs;
568:       void           (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[field];
569:       void           (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[field];
570:       /* Conforming batches */
571:       PetscInt blockSize  = numBasisFuncs*numQuadPoints;
572:       PetscInt numBlocks  = 1;
573:       PetscInt batchSize  = numBlocks * blockSize;
574:       PetscInt numBatches = numBatchesTmp;
575:       PetscInt numChunks  = numPoints / (numBatches*batchSize);
576:       /* Remainder */
577:       PetscInt numRemainder = numPoints % (numBatches * batchSize);
578:       PetscInt offset       = numPoints - numRemainder;

580:       (*mesh->integrateBdResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quadBd, u, v0, n, J, invJ, detJ, f0, f1, elemVec);
581:       (*mesh->integrateBdResidualFEM)(numRemainder, numFields, field, quadBd, &u[offset*cellDof], &v0[offset*dim], &n[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
582:                                              f0, f1, &elemVec[offset*cellDof]);
583:     }
584:     for (p = 0; p < numPoints; ++p) {
585:       const PetscInt point = points[p];

587:       if (mesh->printFEM > 1) {DMPrintCellVector(point, "Residual", cellDof, &elemVec[p*cellDof]);}
588:       DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[p*cellDof], ADD_VALUES);
589:     }
590:     ISRestoreIndices(pointIS, &points);
591:     ISDestroy(&pointIS);
592:     PetscFree7(u,v0,n,J,invJ,detJ,elemVec);
593:   }
594:   if (mesh->printFEM) {
595:     PetscMPIInt rank, numProcs;
596:     PetscInt    p;

598:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
599:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
600:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Residual:\n");
601:     for (p = 0; p < numProcs; ++p) {
602:       if (p == rank) {
603:         Vec f;

605:         VecDuplicate(F, &f);
606:         VecCopy(F, f);
607:         VecChop(f, 1.0e-10);
608:         VecView(f, PETSC_VIEWER_STDOUT_SELF);
609:         VecDestroy(&f);
610:         PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
611:       }
612:       PetscBarrier((PetscObject) dm);
613:     }
614:   }
615:   /* PetscLogEventEnd(ResidualFEMEvent,0,0,0,0); */
616:   return(0);
617: }

621: /*@C
622:   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user

624:   Input Parameters:
625: + dm - The mesh
626: . J  - The Jacobian shell matrix
627: . X  - Local input vector
628: - user - The user context

630:   Output Parameter:
631: . F  - Local output vector

633:   Note:
634:   The second member of the user context must be an FEMContext.

636:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
637:   like a GPU, or vectorize on a multicore machine.

639:   Level: developer

641: .seealso: DMPlexComputeResidualFEM()
642: @*/
643: PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
644: {
645:   DM_Plex         *mesh = (DM_Plex*) dm->data;
646:   PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
647:   PetscQuadrature *quad = fem->quad;
648:   PetscSection     section;
649:   JacActionCtx    *jctx;
650:   PetscReal       *v0, *J, *invJ, *detJ;
651:   PetscScalar     *elemVec, *u, *a;
652:   PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
653:   PetscInt         cellDof = 0;
654:   PetscErrorCode   ierr;

657:   /* PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0); */
658:   MatShellGetContext(Jac, &jctx);
659:   DMPlexGetDimension(dm, &dim);
660:   DMGetDefaultSection(dm, &section);
661:   PetscSectionGetNumFields(section, &numFields);
662:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
663:   numCells = cEnd - cStart;
664:   for (field = 0; field < numFields; ++field) {
665:     cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
666:   }
667:   VecSet(F, 0.0);
668:   PetscMalloc7(numCells*cellDof,PetscScalar,&u,numCells*cellDof,PetscScalar,&a,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);
669:   for (c = cStart; c < cEnd; ++c) {
670:     PetscScalar *x;
671:     PetscInt     i;

673:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
674:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
675:     DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);
676:     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
677:     DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);
678:     DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);
679:     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
680:     DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);
681:   }
682:   for (field = 0; field < numFields; ++field) {
683:     const PetscInt numQuadPoints = quad[field].numQuadPoints;
684:     const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
685:     /* Conforming batches */
686:     PetscInt blockSize  = numBasisFuncs*numQuadPoints;
687:     PetscInt numBlocks  = 1;
688:     PetscInt batchSize  = numBlocks * blockSize;
689:     PetscInt numBatches = numBatchesTmp;
690:     PetscInt numChunks  = numCells / (numBatches*batchSize);
691:     /* Remainder */
692:     PetscInt numRemainder = numCells % (numBatches * batchSize);
693:     PetscInt offset       = numCells - numRemainder;

695:     (*mesh->integrateJacobianActionFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, a, v0, J, invJ, detJ, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);
696:     (*mesh->integrateJacobianActionFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &a[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
697:                                                fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);
698:   }
699:   for (c = cStart; c < cEnd; ++c) {
700:     if (mesh->printFEM > 1) {DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);}
701:     DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);
702:   }
703:   PetscFree7(u,a,v0,J,invJ,detJ,elemVec);
704:   if (mesh->printFEM) {
705:     PetscMPIInt rank, numProcs;
706:     PetscInt    p;

708:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
709:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
710:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");
711:     for (p = 0; p < numProcs; ++p) {
712:       if (p == rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
713:       PetscBarrier((PetscObject) dm);
714:     }
715:   }
716:   /* PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0); */
717:   return(0);
718: }

722: /*@
723:   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.

725:   Input Parameters:
726: + dm - The mesh
727: . X  - Local input vector
728: - user - The user context

730:   Output Parameter:
731: . Jac  - Jacobian matrix

733:   Note:
734:   The second member of the user context must be an FEMContext.

736:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
737:   like a GPU, or vectorize on a multicore machine.

739:   Level: developer

741: .seealso: FormFunctionLocal()
742: @*/
743: PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user)
744: {
745:   DM_Plex         *mesh = (DM_Plex*) dm->data;
746:   PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
747:   PetscQuadrature *quad = fem->quad;
748:   PetscSection     section;
749:   PetscReal       *v0, *J, *invJ, *detJ;
750:   PetscScalar     *elemMat, *u;
751:   PetscInt         dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c;
752:   PetscInt         cellDof = 0, numComponents = 0;
753:   PetscBool        isShell;
754:   PetscErrorCode   ierr;

757:   /* PetscLogEventBegin(JacobianFEMEvent,0,0,0,0); */
758:   DMPlexGetDimension(dm, &dim);
759:   DMGetDefaultSection(dm, &section);
760:   PetscSectionGetNumFields(section, &numFields);
761:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
762:   numCells = cEnd - cStart;
763:   for (field = 0; field < numFields; ++field) {
764:     cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
765:     numComponents += quad[field].numComponents;
766:   }
767:   DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);
768:   MatZeroEntries(JacP);
769:   PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof*cellDof,PetscScalar,&elemMat);
770:   for (c = cStart; c < cEnd; ++c) {
771:     PetscScalar *x;
772:     PetscInt     i;

774:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
775:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
776:     DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);

778:     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
779:     DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);
780:   }
781:   PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));
782:   for (fieldI = 0; fieldI < numFields; ++fieldI) {
783:     const PetscInt numQuadPoints = quad[fieldI].numQuadPoints;
784:     const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs;
785:     PetscInt       fieldJ;

787:     for (fieldJ = 0; fieldJ < numFields; ++fieldJ) {
788:       void (*g0)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*numFields+fieldJ];
789:       void (*g1)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*numFields+fieldJ];
790:       void (*g2)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*numFields+fieldJ];
791:       void (*g3)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*numFields+fieldJ];
792:       /* Conforming batches */
793:       PetscInt blockSize  = numBasisFuncs*numQuadPoints;
794:       PetscInt numBlocks  = 1;
795:       PetscInt batchSize  = numBlocks * blockSize;
796:       PetscInt numBatches = numBatchesTmp;
797:       PetscInt numChunks  = numCells / (numBatches*batchSize);
798:       /* Remainder */
799:       PetscInt numRemainder = numCells % (numBatches * batchSize);
800:       PetscInt offset       = numCells - numRemainder;

802:       (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);
803:       (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
804:                                            g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);
805:     }
806:   }
807:   for (c = cStart; c < cEnd; ++c) {
808:     if (mesh->printFEM > 1) {DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);}
809:     DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);
810:   }
811:   PetscFree6(u,v0,J,invJ,detJ,elemMat);

813:   /* Assemble matrix, using the 2-step process:
814:        MatAssemblyBegin(), MatAssemblyEnd(). */
815:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
816:   MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);

818:   if (mesh->printFEM) {
819:     PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");
820:     MatChop(JacP, 1.0e-10);
821:     MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
822:   }
823:   /* PetscLogEventEnd(JacobianFEMEvent,0,0,0,0); */
824:   PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);
825:   if (isShell) {
826:     JacActionCtx *jctx;

828:     MatShellGetContext(Jac, &jctx);
829:     VecCopy(X, jctx->u);
830:   }
831:   *str = SAME_NONZERO_PATTERN;
832:   return(0);
833: }