1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petsc/private/dtimpl.h> /*I "petscdt.h" I*/ 3 4 static PetscErrorCode PetscSpacePointView_Ascii(PetscSpace sp, PetscViewer viewer) 5 { 6 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data; 7 PetscViewerFormat format; 8 9 PetscFunctionBegin; 10 PetscCall(PetscViewerGetFormat(viewer, &format)); 11 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 12 PetscCall(PetscViewerASCIIPrintf(viewer, "Point space in dimension %" PetscInt_FMT ":\n", sp->Nv)); 13 PetscCall(PetscViewerASCIIPushTab(viewer)); 14 PetscCall(PetscQuadratureView(pt->quad, viewer)); 15 PetscCall(PetscViewerASCIIPopTab(viewer)); 16 } else PetscCall(PetscViewerASCIIPrintf(viewer, "Point space in dimension %" PetscInt_FMT " on %" PetscInt_FMT " points\n", sp->Nv, pt->quad->numPoints)); 17 PetscFunctionReturn(PETSC_SUCCESS); 18 } 19 20 static PetscErrorCode PetscSpaceView_Point(PetscSpace sp, PetscViewer viewer) 21 { 22 PetscBool isascii; 23 24 PetscFunctionBegin; 25 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 26 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 27 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 28 if (isascii) PetscCall(PetscSpacePointView_Ascii(sp, viewer)); 29 PetscFunctionReturn(PETSC_SUCCESS); 30 } 31 32 static PetscErrorCode PetscSpaceSetUp_Point(PetscSpace sp) 33 { 34 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data; 35 36 PetscFunctionBegin; 37 if (!pt->quad->points && sp->degree >= 0) { 38 PetscCall(PetscQuadratureDestroy(&pt->quad)); 39 PetscCall(PetscDTStroudConicalQuadrature(sp->Nv, sp->Nc, PetscMax(sp->degree + 1, 1), -1.0, 1.0, &pt->quad)); 40 } 41 PetscFunctionReturn(PETSC_SUCCESS); 42 } 43 44 static PetscErrorCode PetscSpaceDestroy_Point(PetscSpace sp) 45 { 46 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data; 47 48 PetscFunctionBegin; 49 PetscCall(PetscQuadratureDestroy(&pt->quad)); 50 PetscCall(PetscFree(pt)); 51 PetscFunctionReturn(PETSC_SUCCESS); 52 } 53 54 static PetscErrorCode PetscSpaceGetDimension_Point(PetscSpace sp, PetscInt *dim) 55 { 56 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data; 57 58 PetscFunctionBegin; 59 *dim = pt->quad->numPoints; 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 63 static PetscErrorCode PetscSpaceEvaluate_Point(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 64 { 65 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data; 66 PetscInt dim = sp->Nv, pdim = pt->quad->numPoints, d, p, i, c; 67 68 PetscFunctionBegin; 69 PetscCheck(npoints == pt->quad->numPoints, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot evaluate Point space on %" PetscInt_FMT " points != %" PetscInt_FMT " size", npoints, pt->quad->numPoints); 70 PetscCall(PetscArrayzero(B, npoints * pdim)); 71 for (p = 0; p < npoints; ++p) { 72 for (i = 0; i < pdim; ++i) { 73 for (d = 0; d < dim; ++d) { 74 if (PetscAbsReal(points[p * dim + d] - pt->quad->points[p * dim + d]) > 1.0e-10) break; 75 } 76 if (d >= dim) { 77 B[p * pdim + i] = 1.0; 78 break; 79 } 80 } 81 } 82 /* Replicate for other components */ 83 for (c = 1; c < sp->Nc; ++c) { 84 for (p = 0; p < npoints; ++p) { 85 for (i = 0; i < pdim; ++i) B[(c * npoints + p) * pdim + i] = B[p * pdim + i]; 86 } 87 } 88 if (D) PetscCall(PetscArrayzero(D, npoints * pdim * dim)); 89 if (H) PetscCall(PetscArrayzero(H, npoints * pdim * dim * dim)); 90 PetscFunctionReturn(PETSC_SUCCESS); 91 } 92 93 static PetscErrorCode PetscSpaceInitialize_Point(PetscSpace sp) 94 { 95 PetscFunctionBegin; 96 sp->ops->setfromoptions = NULL; 97 sp->ops->setup = PetscSpaceSetUp_Point; 98 sp->ops->view = PetscSpaceView_Point; 99 sp->ops->destroy = PetscSpaceDestroy_Point; 100 sp->ops->getdimension = PetscSpaceGetDimension_Point; 101 sp->ops->evaluate = PetscSpaceEvaluate_Point; 102 PetscFunctionReturn(PETSC_SUCCESS); 103 } 104 105 /*MC 106 PETSCSPACEPOINT = "point" - A `PetscSpace` object that encapsulates functions defined on a set of quadrature points. 107 108 Level: intermediate 109 110 .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()` 111 M*/ 112 113 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Point(PetscSpace sp) 114 { 115 PetscSpace_Point *pt; 116 117 PetscFunctionBegin; 118 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 119 PetscCall(PetscNew(&pt)); 120 sp->data = pt; 121 122 sp->Nv = 0; 123 sp->maxDegree = PETSC_INT_MAX; 124 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &pt->quad)); 125 PetscCall(PetscQuadratureSetData(pt->quad, 0, 1, 0, NULL, NULL)); 126 127 PetscCall(PetscSpaceInitialize_Point(sp)); 128 PetscFunctionReturn(PETSC_SUCCESS); 129 } 130 131 /*@ 132 PetscSpacePointSetPoints - Sets the evaluation points for the space to coincide with the points of a quadrature rule 133 134 Logically Collective 135 136 Input Parameters: 137 + sp - The `PetscSpace` 138 - q - The `PetscQuadrature` defining the points 139 140 Level: intermediate 141 142 .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()` 143 @*/ 144 PetscErrorCode PetscSpacePointSetPoints(PetscSpace sp, PetscQuadrature q) 145 { 146 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data; 147 148 PetscFunctionBegin; 149 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 150 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 2); 151 PetscCall(PetscQuadratureDestroy(&pt->quad)); 152 PetscCall(PetscQuadratureDuplicate(q, &pt->quad)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 /*@ 157 PetscSpacePointGetPoints - Gets the evaluation points for the space as the points of a quadrature rule 158 159 Logically Collective 160 161 Input Parameter: 162 . sp - The `PetscSpace` 163 164 Output Parameter: 165 . q - The `PetscQuadrature` defining the points 166 167 Level: intermediate 168 169 .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()` 170 @*/ 171 PetscErrorCode PetscSpacePointGetPoints(PetscSpace sp, PetscQuadrature *q) 172 { 173 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data; 174 175 PetscFunctionBegin; 176 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 177 PetscAssertPointer(q, 2); 178 *q = pt->quad; 179 PetscFunctionReturn(PETSC_SUCCESS); 180 } 181