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