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(0); 18 } 19 20 static PetscErrorCode PetscSpaceView_Point(PetscSpace sp, PetscViewer viewer) 21 { 22 PetscBool iascii; 23 24 PetscFunctionBegin; 25 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 26 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 27 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii)); 28 if (iascii) PetscCall(PetscSpacePointView_Ascii(sp, viewer)); 29 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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) {B[p*pdim+i] = 1.0; break;} 77 } 78 } 79 /* Replicate for other components */ 80 for (c = 1; c < sp->Nc; ++c) { 81 for (p = 0; p < npoints; ++p) { 82 for (i = 0; i < pdim; ++i) { 83 B[(c*npoints + p)*pdim + i] = B[p*pdim + i]; 84 } 85 } 86 } 87 if (D) PetscCall(PetscArrayzero(D, npoints*pdim*dim)); 88 if (H) PetscCall(PetscArrayzero(H, npoints*pdim*dim*dim)); 89 PetscFunctionReturn(0); 90 } 91 92 static PetscErrorCode PetscSpaceInitialize_Point(PetscSpace sp) 93 { 94 PetscFunctionBegin; 95 sp->ops->setfromoptions = NULL; 96 sp->ops->setup = PetscSpaceSetUp_Point; 97 sp->ops->view = PetscSpaceView_Point; 98 sp->ops->destroy = PetscSpaceDestroy_Point; 99 sp->ops->getdimension = PetscSpaceGetDimension_Point; 100 sp->ops->evaluate = PetscSpaceEvaluate_Point; 101 PetscFunctionReturn(0); 102 } 103 104 /*MC 105 PETSCSPACEPOINT = "point" - A PetscSpace object that encapsulates functions defined on a set of quadrature points. 106 107 Level: intermediate 108 109 .seealso: `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()` 110 M*/ 111 112 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Point(PetscSpace sp) 113 { 114 PetscSpace_Point *pt; 115 116 PetscFunctionBegin; 117 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 118 PetscCall(PetscNewLog(sp,&pt)); 119 sp->data = pt; 120 121 sp->Nv = 0; 122 sp->maxDegree = PETSC_MAX_INT; 123 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &pt->quad)); 124 PetscCall(PetscQuadratureSetData(pt->quad, 0, 1, 0, NULL, NULL)); 125 126 PetscCall(PetscSpaceInitialize_Point(sp)); 127 PetscFunctionReturn(0); 128 } 129 130 /*@ 131 PetscSpacePointSetPoints - Sets the evaluation points for the space to coincide with the points of a quadrature rule 132 133 Logically collective 134 135 Input Parameters: 136 + sp - The PetscSpace 137 - q - The PetscQuadrature defining the points 138 139 Level: intermediate 140 141 .seealso: `PetscSpaceCreate()`, `PetscSpaceSetType()` 142 @*/ 143 PetscErrorCode PetscSpacePointSetPoints(PetscSpace sp, PetscQuadrature q) 144 { 145 PetscSpace_Point *pt = (PetscSpace_Point *) sp->data; 146 147 PetscFunctionBegin; 148 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 149 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 2); 150 PetscCall(PetscQuadratureDestroy(&pt->quad)); 151 PetscCall(PetscQuadratureDuplicate(q, &pt->quad)); 152 PetscFunctionReturn(0); 153 } 154 155 /*@ 156 PetscSpacePointGetPoints - Gets the evaluation points for the space as the points of a quadrature rule 157 158 Logically collective 159 160 Input Parameter: 161 . sp - The PetscSpace 162 163 Output Parameter: 164 . q - The PetscQuadrature defining the points 165 166 Level: intermediate 167 168 .seealso: `PetscSpaceCreate()`, `PetscSpaceSetType()` 169 @*/ 170 PetscErrorCode PetscSpacePointGetPoints(PetscSpace sp, PetscQuadrature *q) 171 { 172 PetscSpace_Point *pt = (PetscSpace_Point *) sp->data; 173 174 PetscFunctionBegin; 175 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 176 PetscValidPointer(q, 2); 177 *q = pt->quad; 178 PetscFunctionReturn(0); 179 } 180