Actual source code: plexceed.c

  1: #include <petsc/private/dmpleximpl.h>

  3: #ifdef PETSC_HAVE_LIBCEED
  4: #include <petscdmplexceed.h>

  6: /* Define the map from the local vector (Lvector) to the cells (Evector) */
  7: PetscErrorCode DMPlexGetCeedRestriction(DM dm, CeedElemRestriction *ERestrict)
  8: {

 14:   if (!dm->ceedERestrict) {
 15:     PetscDS      ds;
 16:     PetscFE      fe;
 17:     PetscSpace   sp;
 18:     PetscSection s;
 19:     PetscInt     Nf, *Nc, c, P, cStart, cEnd, Ncell, cell, lsize, *erestrict, eoffset;
 20:     PetscBool    simplex;
 21:     Ceed         ceed;

 23:     DMPlexIsSimplex(dm, &simplex);
 24:     if (simplex) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "LibCEED does not work with simplices");
 25:     DMGetDS(dm, &ds);
 26:     PetscDSGetNumFields(ds, &Nf);
 27:     if (Nf != 1) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "LibCEED only works with one field right now");
 28:     PetscDSGetComponents(ds, &Nc);
 29:     PetscDSGetDiscretization(ds, 0, (PetscObject *) &fe);
 30:     PetscFEGetBasisSpace(fe, &sp);
 31:     PetscSpaceGetDegree(sp, &P, NULL);
 32:     ++P;
 33:     DMGetLocalSection(dm, &s);
 34:     DMPlexGetHeightStratum(dm, 0, &cStart,& cEnd);
 35:     Ncell = cEnd - cStart;

 37:     PetscMalloc1(Ncell*P*P, &erestrict);
 38:     for (cell = cStart, eoffset = 0; cell < cEnd; ++cell) {
 39:       PetscInt Ni, *ind, i;

 41:       DMPlexGetClosureIndices(dm, s, s, cell, PETSC_TRUE, &Ni, &ind, NULL, NULL);
 42:       for (i = 0; i < Ni; i += Nc[0]) {
 43:         for (c = 0; c < Nc[0]; ++c) {
 44:           if (ind[i+c] != ind[i] + c) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %D closure indices not interlaced", cell);
 45:         }
 46:         erestrict[eoffset++] = ind[i];
 47:       }
 48:       DMPlexRestoreClosureIndices(dm, s, s, cell, PETSC_TRUE, &Ni, &ind, NULL, NULL);
 49:     }
 50:     if (eoffset != Ncell*P*P) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Size of array %D != %D restricted dofs", Ncell*P*P, eoffset);

 52:     DMGetCeed(dm, &ceed);
 53:     PetscSectionGetStorageSize(s, &lsize);
 54:     CeedElemRestrictionCreate(ceed, Ncell, P*P, Nc[0], 1, lsize, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, &dm->ceedERestrict);
 55:     PetscFree(erestrict);
 56:   }
 57:   *ERestrict = dm->ceedERestrict;
 58:   return(0);
 59: }

 61: #endif