1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/snesimpl.h>
  3: #include <petscds.h>
  4: #include <petsc/private/petscimpl.h>
  5: #include <petsc/private/petscfeimpl.h>
  6: #include <petscdmceed.h>
  7: #include <petscdmplexceed.h>

  9: PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, PetscCtx ctx)
 10: {
 11:   Ceed       ceed;
 12:   DMCeed     sd = dm->dmceed;
 13:   CeedVector clocX, clocF;

 15:   PetscFunctionBegin;
 16:   PetscCall(DMGetCeed(dm, &ceed));
 17:   PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
 18:   PetscCall(DMCeedComputeGeometry(dm, sd));

 20:   PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
 21:   PetscCall(VecGetCeedVector(locF, ceed, &clocF));
 22:   PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
 23:   PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
 24:   PetscCall(VecRestoreCeedVector(locF, &clocF));

 26:   {
 27:     DM_Plex *mesh = (DM_Plex *)dm->data;

 29:     if (mesh->printFEM) {
 30:       PetscSection section;
 31:       Vec          locFbc;
 32:       PetscInt     pStart, pEnd, p, maxDof;
 33:       PetscScalar *zeroes;

 35:       PetscCall(DMGetLocalSection(dm, &section));
 36:       PetscCall(VecDuplicate(locF, &locFbc));
 37:       PetscCall(VecCopy(locF, locFbc));
 38:       PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
 39:       PetscCall(PetscSectionGetMaxDof(section, &maxDof));
 40:       PetscCall(PetscCalloc1(maxDof, &zeroes));
 41:       for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES));
 42:       PetscCall(PetscFree(zeroes));
 43:       PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc));
 44:       PetscCall(VecDestroy(&locFbc));
 45:     }
 46:   }
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }