xref: /petsc/src/snes/utils/libceed/dmplexsnesceed.c (revision ce78bad369055609e946c9d2c25ea67a45873e27)
1*ce78bad3SBarry Smith #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2*ce78bad3SBarry Smith #include <petsc/private/snesimpl.h>   /*I "petscsnes.h"   I*/
3*ce78bad3SBarry Smith #include <petscds.h>
4*ce78bad3SBarry Smith #include <petsc/private/petscimpl.h>
5*ce78bad3SBarry Smith #include <petsc/private/petscfeimpl.h>
6*ce78bad3SBarry Smith #include <petscdmceed.h>
7*ce78bad3SBarry Smith #include <petscdmplexceed.h>
8*ce78bad3SBarry Smith 
9*ce78bad3SBarry Smith PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user)
10*ce78bad3SBarry Smith {
11*ce78bad3SBarry Smith   Ceed       ceed;
12*ce78bad3SBarry Smith   DMCeed     sd = dm->dmceed;
13*ce78bad3SBarry Smith   CeedVector clocX, clocF;
14*ce78bad3SBarry Smith 
15*ce78bad3SBarry Smith   PetscFunctionBegin;
16*ce78bad3SBarry Smith   PetscCall(DMGetCeed(dm, &ceed));
17*ce78bad3SBarry Smith   PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
18*ce78bad3SBarry Smith   PetscCall(DMCeedComputeGeometry(dm, sd));
19*ce78bad3SBarry Smith 
20*ce78bad3SBarry Smith   PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
21*ce78bad3SBarry Smith   PetscCall(VecGetCeedVector(locF, ceed, &clocF));
22*ce78bad3SBarry Smith   PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
23*ce78bad3SBarry Smith   PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
24*ce78bad3SBarry Smith   PetscCall(VecRestoreCeedVector(locF, &clocF));
25*ce78bad3SBarry Smith 
26*ce78bad3SBarry Smith   {
27*ce78bad3SBarry Smith     DM_Plex *mesh = (DM_Plex *)dm->data;
28*ce78bad3SBarry Smith 
29*ce78bad3SBarry Smith     if (mesh->printFEM) {
30*ce78bad3SBarry Smith       PetscSection section;
31*ce78bad3SBarry Smith       Vec          locFbc;
32*ce78bad3SBarry Smith       PetscInt     pStart, pEnd, p, maxDof;
33*ce78bad3SBarry Smith       PetscScalar *zeroes;
34*ce78bad3SBarry Smith 
35*ce78bad3SBarry Smith       PetscCall(DMGetLocalSection(dm, &section));
36*ce78bad3SBarry Smith       PetscCall(VecDuplicate(locF, &locFbc));
37*ce78bad3SBarry Smith       PetscCall(VecCopy(locF, locFbc));
38*ce78bad3SBarry Smith       PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
39*ce78bad3SBarry Smith       PetscCall(PetscSectionGetMaxDof(section, &maxDof));
40*ce78bad3SBarry Smith       PetscCall(PetscCalloc1(maxDof, &zeroes));
41*ce78bad3SBarry Smith       for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES));
42*ce78bad3SBarry Smith       PetscCall(PetscFree(zeroes));
43*ce78bad3SBarry Smith       PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc));
44*ce78bad3SBarry Smith       PetscCall(VecDestroy(&locFbc));
45*ce78bad3SBarry Smith     }
46*ce78bad3SBarry Smith   }
47*ce78bad3SBarry Smith   PetscFunctionReturn(PETSC_SUCCESS);
48*ce78bad3SBarry Smith }
49