1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #ifdef PETSC_HAVE_LIBCEED 4 #include <petscdmplexceed.h> 5 6 /* Define the map from the local vector (Lvector) to the cells (Evector) */ 7 PetscErrorCode DMPlexGetCeedRestriction(DM dm, CeedElemRestriction *ERestrict) 8 { 9 PetscErrorCode ierr; 10 11 PetscFunctionBeginUser; 12 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 13 PetscValidPointer(ERestrict, 2); 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; 22 23 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 24 if (simplex) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "LibCEED does not work with simplices"); 25 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 26 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 27 if (Nf != 1) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "LibCEED only works with one field right now"); 28 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 29 ierr = PetscDSGetDiscretization(ds, 0, (PetscObject *) &fe);CHKERRQ(ierr); 30 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 31 ierr = PetscSpaceGetDegree(sp, &P, NULL);CHKERRQ(ierr); 32 ++P; 33 ierr = DMGetLocalSection(dm, &s);CHKERRQ(ierr); 34 ierr = DMPlexGetHeightStratum(dm, 0, &cStart,& cEnd);CHKERRQ(ierr); 35 Ncell = cEnd - cStart; 36 37 ierr = PetscMalloc1(Ncell*P*P, &erestrict);CHKERRQ(ierr); 38 for (cell = cStart, eoffset = 0; cell < cEnd; ++cell) { 39 PetscInt Ni, *ind, i; 40 41 ierr = DMPlexGetClosureIndices(dm, s, s, cell, PETSC_TRUE, &Ni, &ind, NULL, NULL);CHKERRQ(ierr); 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 ierr = DMPlexRestoreClosureIndices(dm, s, s, cell, PETSC_TRUE, &Ni, &ind, NULL, NULL);CHKERRQ(ierr); 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); 51 52 ierr = DMGetCeed(dm, &ceed);CHKERRQ(ierr); 53 ierr = PetscSectionGetStorageSize(s, &lsize);CHKERRQ(ierr); 54 ierr = CeedElemRestrictionCreate(ceed, Ncell, P*P, Nc[0], 1, lsize, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, &dm->ceedERestrict);CHKERRQ(ierr); 55 ierr = PetscFree(erestrict);CHKERRQ(ierr); 56 } 57 *ERestrict = dm->ceedERestrict; 58 PetscFunctionReturn(0); 59 } 60 61 #endif 62