1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
3 #include <petscds.h>
4 #include <petsc/private/petscimpl.h>
5 #include <petsc/private/petscfeimpl.h>
6 #include <petscdmceed.h>
7 #include <petscdmplexceed.h>
8
DMPlexSNESComputeResidualCEED(DM dm,Vec locX,Vec locF,PetscCtx ctx)9 PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, PetscCtx ctx)
10 {
11 Ceed ceed;
12 DMCeed sd = dm->dmceed;
13 CeedVector clocX, clocF;
14
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));
19
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));
25
26 {
27 DM_Plex *mesh = (DM_Plex *)dm->data;
28
29 if (mesh->printFEM) {
30 PetscSection section;
31 Vec locFbc;
32 PetscInt pStart, pEnd, p, maxDof;
33 PetscScalar *zeroes;
34
35 PetscCall(DMGetLocalSection(dm, §ion));
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 }
49