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