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