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