1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 4 typedef struct { 5 PetscInt dummy; 6 } PetscDualSpace_Refined; 7 8 /*@ 9 PetscDualSpaceRefinedSetCellSpaces - Set the dual spaces for the closures of each of the cells 10 in the multicell DM of a PetscDualSpace 11 12 Collective on PetscDualSpace 13 14 Input Parameters: 15 + sp - a PetscDualSpace 16 - cellSpaces - one PetscDualSpace for each of the cells. The reference count of each cell space will be incremented, 17 so the user is still responsible for these spaces afterwards 18 19 Level: intermediate 20 21 .seealso: `PetscFERefine()` 22 @*/ 23 PetscErrorCode PetscDualSpaceRefinedSetCellSpaces(PetscDualSpace sp, const PetscDualSpace cellSpaces[]) { 24 PetscFunctionBegin; 25 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 26 PetscValidPointer(cellSpaces, 2); 27 PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change cell spaces after setup is called"); 28 PetscTryMethod(sp, "PetscDualSpaceRefinedSetCellSpaces_C", (PetscDualSpace, const PetscDualSpace[]), (sp, cellSpaces)); 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode PetscDualSpaceRefinedSetCellSpaces_Refined(PetscDualSpace sp, const PetscDualSpace cellSpaces[]) { 33 DM dm; 34 PetscInt pStart, pEnd; 35 PetscInt cStart, cEnd, c; 36 37 PetscFunctionBegin; 38 dm = sp->dm; 39 PetscCheck(dm, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "PetscDualSpace must have a DM (PetscDualSpaceSetDM()) before calling PetscDualSpaceRefinedSetCellSpaces"); 40 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 41 if (!sp->pointSpaces) { PetscCall(PetscCalloc1(pEnd - pStart, &(sp->pointSpaces))); } 42 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 43 for (c = 0; c < cEnd - cStart; c++) { 44 PetscCall(PetscObjectReference((PetscObject)cellSpaces[c])); 45 PetscCall(PetscDualSpaceDestroy(&(sp->pointSpaces[c + cStart - pStart]))); 46 sp->pointSpaces[c + cStart - pStart] = cellSpaces[c]; 47 } 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode PetscDualSpaceDestroy_Refined(PetscDualSpace sp) { 52 PetscDualSpace_Refined *ref = (PetscDualSpace_Refined *)sp->data; 53 54 PetscFunctionBegin; 55 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceRefinedSetCellSpaces_C", NULL)); 56 PetscCall(PetscFree(ref)); 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode PetscDualSpaceSetUp_Refined(PetscDualSpace sp) { 61 PetscInt pStart, pEnd, depth; 62 PetscInt cStart, cEnd, c, spdim; 63 PetscInt h; 64 DM dm; 65 PetscSection section; 66 67 PetscFunctionBegin; 68 PetscCall(PetscDualSpaceGetDM(sp, &dm)); 69 PetscCall(DMPlexGetDepth(dm, &depth)); 70 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 71 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 72 for (c = cStart; c < cEnd; c++) { 73 if (sp->pointSpaces[c - pStart]) { 74 PetscInt ccStart, ccEnd; 75 PetscCheck(sp->pointSpaces[c - pStart]->k == sp->k, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have the same form degree as the refined dual space"); 76 PetscCheck(sp->pointSpaces[c - pStart]->Nc == sp->Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have the same number of components as the refined dual space"); 77 PetscCall(DMPlexGetHeightStratum(sp->pointSpaces[c - pStart]->dm, 0, &ccStart, &ccEnd)); 78 PetscCheck(ccEnd - ccStart == 1, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have a single cell themselves"); 79 } 80 } 81 for (c = cStart; c < cEnd; c++) { 82 if (sp->pointSpaces[c - pStart]) { 83 PetscBool cUniform; 84 85 PetscCall(PetscDualSpaceGetUniform(sp->pointSpaces[c - pStart], &cUniform)); 86 if (!cUniform) break; 87 } 88 if ((c > cStart) && sp->pointSpaces[c - pStart] != sp->pointSpaces[c - 1 - pStart]) break; 89 } 90 if (c < cEnd) sp->uniform = PETSC_FALSE; 91 for (h = 0; h < depth; h++) { 92 PetscInt hStart, hEnd; 93 94 PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd)); 95 for (c = hStart; c < hEnd; c++) { 96 PetscInt coneSize, e; 97 PetscDualSpace cspace = sp->pointSpaces[c - pStart]; 98 const PetscInt *cone; 99 const PetscInt *refCone; 100 101 if (!cspace) continue; 102 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 103 PetscCall(DMPlexGetCone(dm, c, &cone)); 104 PetscCall(DMPlexGetCone(cspace->dm, 0, &refCone)); 105 for (e = 0; e < coneSize; e++) { 106 PetscInt point = cone[e]; 107 PetscInt refpoint = refCone[e]; 108 PetscDualSpace espace; 109 110 PetscCall(PetscDualSpaceGetPointSubspace(cspace, refpoint, &espace)); 111 if (sp->pointSpaces[point - pStart] == NULL) { 112 PetscCall(PetscObjectReference((PetscObject)espace)); 113 sp->pointSpaces[point - pStart] = espace; 114 } 115 } 116 } 117 } 118 PetscCall(PetscDualSpaceGetSection(sp, §ion)); 119 PetscCall(PetscDualSpaceGetDimension(sp, &spdim)); 120 PetscCall(PetscMalloc1(spdim, &(sp->functional))); 121 PetscCall(PetscDualSpacePushForwardSubspaces_Internal(sp, pStart, pEnd)); 122 PetscFunctionReturn(0); 123 } 124 125 static PetscErrorCode PetscDualSpaceRefinedView_Ascii(PetscDualSpace sp, PetscViewer viewer) { 126 PetscFunctionBegin; 127 if (sp->dm && sp->pointSpaces) { 128 PetscInt pStart, pEnd; 129 PetscInt cStart, cEnd, c; 130 131 PetscCall(DMPlexGetChart(sp->dm, &pStart, &pEnd)); 132 PetscCall(DMPlexGetHeightStratum(sp->dm, 0, &cStart, &cEnd)); 133 PetscCall(PetscViewerASCIIPrintf(viewer, "Refined dual space:\n")); 134 PetscCall(PetscViewerASCIIPushTab(viewer)); 135 for (c = cStart; c < cEnd; c++) { 136 if (!sp->pointSpaces[c - pStart]) { 137 PetscCall(PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT " not set yet\n", c)); 138 } else { 139 PetscCall(PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT ":ot set yet\n", c)); 140 PetscCall(PetscViewerASCIIPushTab(viewer)); 141 PetscCall(PetscDualSpaceView(sp->pointSpaces[c - pStart], viewer)); 142 PetscCall(PetscViewerASCIIPopTab(viewer)); 143 } 144 } 145 PetscCall(PetscViewerASCIIPopTab(viewer)); 146 } else PetscCall(PetscViewerASCIIPrintf(viewer, "Refined dual space: (cell spaces not set yet)\n")); 147 PetscFunctionReturn(0); 148 } 149 150 static PetscErrorCode PetscDualSpaceView_Refined(PetscDualSpace sp, PetscViewer viewer) { 151 PetscBool iascii; 152 153 PetscFunctionBegin; 154 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 155 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 156 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 157 if (iascii) PetscCall(PetscDualSpaceRefinedView_Ascii(sp, viewer)); 158 PetscFunctionReturn(0); 159 } 160 161 static PetscErrorCode PetscDualSpaceInitialize_Refined(PetscDualSpace sp) { 162 PetscFunctionBegin; 163 sp->ops->destroy = PetscDualSpaceDestroy_Refined; 164 sp->ops->view = PetscDualSpaceView_Refined; 165 sp->ops->setfromoptions = NULL; 166 sp->ops->duplicate = NULL; 167 sp->ops->setup = PetscDualSpaceSetUp_Refined; 168 sp->ops->createheightsubspace = NULL; 169 sp->ops->createpointsubspace = NULL; 170 sp->ops->getsymmetries = NULL; 171 sp->ops->apply = PetscDualSpaceApplyDefault; 172 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 173 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault; 174 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault; 175 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault; 176 PetscFunctionReturn(0); 177 } 178 179 /*MC 180 PETSCDUALSPACEREFINED = "refined" - A PetscDualSpace object that defines the joint dual space of a group of cells, usually refined from one larger cell 181 182 Level: intermediate 183 184 .seealso: `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()` 185 M*/ 186 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Refined(PetscDualSpace sp) { 187 PetscDualSpace_Refined *ref; 188 189 PetscFunctionBegin; 190 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 191 PetscCall(PetscNewLog(sp, &ref)); 192 sp->data = ref; 193 194 PetscCall(PetscDualSpaceInitialize_Refined(sp)); 195 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceRefinedSetCellSpaces_C", PetscDualSpaceRefinedSetCellSpaces_Refined)); 196 PetscFunctionReturn(0); 197 } 198