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: /*@C
 10:   DMPlexSNESComputeResidualCEED - Assemble the local residual for a `SNES` on a `DMPLEX` using the libCEED operator attached to the `DM`

 12:   Collective

 14:   Input Parameters:
 15: + dm   - the `DMPLEX` for which libCEED operators have been created by `DMCeedCreate()`
 16: . locX - local solution vector including ghost values
 17: - ctx  - application context (unused)

 19:   Output Parameter:
 20: . locF - local residual vector to be assembled

 22:   Level: developer

 24:   Note:
 25:   This is normally installed as the `SNES` local residual callback via `DMSNESSetFunctionLocal()` when using libCEED for the finite-element evaluation.

 27: .seealso: [](ch_snes), `SNES`, `DMPLEX`, `DMCeedCreate()`, `DMSNESSetFunctionLocal()`, `DMPlexTSComputeRHSFunctionFVMCEED()`
 28: @*/
 29: PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, PetscCtx ctx)
 30: {
 31:   Ceed       ceed;
 32:   DMCeed     sd = dm->dmceed;
 33:   CeedVector clocX, clocF;

 35:   PetscFunctionBegin;
 36:   PetscCall(DMGetCeed(dm, &ceed));
 37:   PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
 38:   PetscCall(DMCeedComputeGeometry(dm, sd));

 40:   PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
 41:   PetscCall(VecGetCeedVector(locF, ceed, &clocF));
 42:   PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
 43:   PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
 44:   PetscCall(VecRestoreCeedVector(locF, &clocF));

 46:   {
 47:     DM_Plex *mesh = (DM_Plex *)dm->data;

 49:     if (mesh->printFEM) {
 50:       PetscSection section;
 51:       Vec          locFbc;
 52:       PetscInt     pStart, pEnd, p, maxDof;
 53:       PetscScalar *zeroes;

 55:       PetscCall(DMGetLocalSection(dm, &section));
 56:       PetscCall(VecDuplicate(locF, &locFbc));
 57:       PetscCall(VecCopy(locF, locFbc));
 58:       PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
 59:       PetscCall(PetscSectionGetMaxDof(section, &maxDof));
 60:       PetscCall(PetscCalloc1(maxDof, &zeroes));
 61:       for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES));
 62:       PetscCall(PetscFree(zeroes));
 63:       PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc));
 64:       PetscCall(VecDestroy(&locFbc));
 65:     }
 66:   }
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }