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