1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 static PetscErrorCode PetscSpaceSetFromOptions_WXY(PetscOptionItems *PetscOptionsObject,PetscSpace sp) 4 { 5 PetscErrorCode ierr; 6 7 PetscFunctionBegin; 8 ierr = PetscOptionsHead(PetscOptionsObject,"PetscSpace WXY options");CHKERRQ(ierr); 9 ierr = PetscOptionsTail();CHKERRQ(ierr); 10 PetscFunctionReturn(0); 11 } 12 13 static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v) 14 { 15 PetscErrorCode ierr; 16 17 PetscFunctionBegin; 18 ierr = PetscViewerASCIIPrintf(v, "WXY space of degree %" PetscInt_FMT "\n", sp->degree);CHKERRQ(ierr); 19 PetscFunctionReturn(0); 20 } 21 22 static PetscErrorCode PetscSpaceView_WXY(PetscSpace sp, PetscViewer viewer) 23 { 24 PetscBool iascii; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 29 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 30 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 31 if (iascii) {ierr = PetscSpacePolynomialView_Ascii(sp, viewer);CHKERRQ(ierr);} 32 PetscFunctionReturn(0); 33 } 34 35 static PetscErrorCode PetscSpaceDestroy_WXY(PetscSpace sp) 36 { 37 PetscSpace_WXY *wxy = (PetscSpace_WXY *) sp->data; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = PetscFree(wxy);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 static PetscErrorCode PetscSpaceSetUp_WXY(PetscSpace sp) 46 { 47 PetscSpace_WXY *wxy = (PetscSpace_WXY *) sp->data; 48 49 PetscFunctionBegin; 50 if (wxy->setupCalled) PetscFunctionReturn(0); 51 PetscCheck(sp->degree >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT " invalid\n", sp->degree); 52 sp->maxDegree = sp->degree; 53 wxy->setupCalled = PETSC_TRUE; 54 PetscFunctionReturn(0); 55 } 56 57 static PetscErrorCode PetscSpaceGetDimension_WXY(PetscSpace sp, PetscInt *dim) 58 { 59 PetscFunctionBegin; 60 *dim = 6; 61 PetscFunctionReturn(0); 62 } 63 64 static PetscErrorCode PetscSpaceEvaluate_WXY(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 65 { 66 PetscSpace_WXY *wxy = (PetscSpace_WXY *) sp->data; 67 PetscInt dim = sp->Nv; 68 PetscInt Nb = 6; 69 PetscInt Nc = 3; 70 71 PetscFunctionBegin; 72 if (!wxy->setupCalled) { 73 PetscErrorCode ierr; 74 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 75 ierr = PetscSpaceEvaluate(sp, npoints, points, B, D, H);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 PetscCheck((sp->Nc == 3) && (sp->Nv == 3), PETSC_COMM_SELF, PETSC_ERR_PLIB, "WXY space must have 3 variables and 3 components"); 79 if (B) { 80 PetscInt p_inc = Nb*Nc; 81 PetscInt b_inc = Nc; 82 PetscInt c_inc = 1; 83 84 for (PetscInt p = 0; p < npoints; p++) { 85 const PetscReal x = points[p*dim+0]; 86 const PetscReal y = points[p*dim+1]; 87 const PetscReal z = points[p*dim+2]; 88 89 /* {2 y z, 0, 0} */ 90 B[p*p_inc + 0*b_inc + 0*c_inc] = 2.*y*z; 91 B[p*p_inc + 0*b_inc + 1*c_inc] = 0.; 92 B[p*p_inc + 0*b_inc + 2*c_inc] = 0.; 93 /* {0, 2 x z, 0} */ 94 B[p*p_inc + 1*b_inc + 0*c_inc] = 0.; 95 B[p*p_inc + 1*b_inc + 1*c_inc] = 2.*x*z; 96 B[p*p_inc + 1*b_inc + 2*c_inc] = 0.; 97 /* {0, 2 y z, -z^2} */ 98 B[p*p_inc + 2*b_inc + 0*c_inc] = 0.; 99 B[p*p_inc + 2*b_inc + 1*c_inc] = 2.*y*z; 100 B[p*p_inc + 2*b_inc + 2*c_inc] = -z*z; 101 /* {2 x z, 0, -z^2} */ 102 B[p*p_inc + 3*b_inc + 0*c_inc] = 2.*x*z; 103 B[p*p_inc + 3*b_inc + 1*c_inc] = 0.; 104 B[p*p_inc + 3*b_inc + 2*c_inc] = -z*z; 105 /* {x^2, x y, -3 x z} */ 106 B[p*p_inc + 4*b_inc + 0*c_inc] = x*x; 107 B[p*p_inc + 4*b_inc + 1*c_inc] = x*y; 108 B[p*p_inc + 4*b_inc + 2*c_inc] = -3.*x*z; 109 /* {x y, y^2, -3 y z} */ 110 B[p*p_inc + 5*b_inc + 0*c_inc] = x*y; 111 B[p*p_inc + 5*b_inc + 1*c_inc] = y*y; 112 B[p*p_inc + 5*b_inc + 2*c_inc] = -3.*y*z; 113 } 114 } 115 if (D) { 116 PetscInt p_inc = Nb*Nc*dim; 117 PetscInt b_inc = Nc*dim; 118 PetscInt c_inc = dim; 119 120 for (PetscInt p = 0; p < npoints; p++) { 121 const PetscReal x = points[p*dim+0]; 122 const PetscReal y = points[p*dim+1]; 123 const PetscReal z = points[p*dim+2]; 124 125 /* {2 y z, 0, 0} */ 126 D[p*p_inc + 0*b_inc + 0*c_inc + 0] = 0.; 127 D[p*p_inc + 0*b_inc + 0*c_inc + 1] = 2.*z; 128 D[p*p_inc + 0*b_inc + 0*c_inc + 2] = 2.*y; 129 D[p*p_inc + 0*b_inc + 1*c_inc + 0] = 0.; 130 D[p*p_inc + 0*b_inc + 1*c_inc + 1] = 0.; 131 D[p*p_inc + 0*b_inc + 1*c_inc + 2] = 0.; 132 D[p*p_inc + 0*b_inc + 2*c_inc + 0] = 0.; 133 D[p*p_inc + 0*b_inc + 2*c_inc + 1] = 0.; 134 D[p*p_inc + 0*b_inc + 2*c_inc + 2] = 0.; 135 /* {0, 2 x z, 0} */ 136 D[p*p_inc + 1*b_inc + 0*c_inc + 0] = 0.; 137 D[p*p_inc + 1*b_inc + 0*c_inc + 1] = 0.; 138 D[p*p_inc + 1*b_inc + 0*c_inc + 2] = 0.; 139 D[p*p_inc + 1*b_inc + 1*c_inc + 0] = 2.*z; 140 D[p*p_inc + 1*b_inc + 1*c_inc + 1] = 0.; 141 D[p*p_inc + 1*b_inc + 1*c_inc + 2] = 2.*x; 142 D[p*p_inc + 1*b_inc + 2*c_inc + 0] = 0.; 143 D[p*p_inc + 1*b_inc + 2*c_inc + 1] = 0.; 144 D[p*p_inc + 1*b_inc + 2*c_inc + 2] = 0.; 145 /* {0, 2 y z, -z^2} */ 146 D[p*p_inc + 2*b_inc + 0*c_inc + 0] = 0.; 147 D[p*p_inc + 2*b_inc + 0*c_inc + 1] = 0.; 148 D[p*p_inc + 2*b_inc + 0*c_inc + 2] = 0.; 149 D[p*p_inc + 2*b_inc + 1*c_inc + 0] = 0.; 150 D[p*p_inc + 2*b_inc + 1*c_inc + 1] = 2.*z; 151 D[p*p_inc + 2*b_inc + 1*c_inc + 2] = 2.*y; 152 D[p*p_inc + 2*b_inc + 2*c_inc + 0] = 0.; 153 D[p*p_inc + 2*b_inc + 2*c_inc + 1] = 0.; 154 D[p*p_inc + 2*b_inc + 2*c_inc + 2] = -2.*z; 155 /* {2 x z, 0, -z^2} */ 156 D[p*p_inc + 3*b_inc + 0*c_inc + 0] = 2.*z; 157 D[p*p_inc + 3*b_inc + 0*c_inc + 1] = 0.; 158 D[p*p_inc + 3*b_inc + 0*c_inc + 2] = 2.*x; 159 D[p*p_inc + 3*b_inc + 1*c_inc + 0] = 0.; 160 D[p*p_inc + 3*b_inc + 1*c_inc + 1] = 0.; 161 D[p*p_inc + 3*b_inc + 1*c_inc + 2] = 0.; 162 D[p*p_inc + 3*b_inc + 2*c_inc + 0] = 0.; 163 D[p*p_inc + 3*b_inc + 2*c_inc + 1] = 0.; 164 D[p*p_inc + 3*b_inc + 2*c_inc + 2] = -2.*z; 165 /* {x^2, x y, -3 x z} */ 166 D[p*p_inc + 4*b_inc + 0*c_inc + 0] = 2.*x; 167 D[p*p_inc + 4*b_inc + 0*c_inc + 1] = 0.; 168 D[p*p_inc + 4*b_inc + 0*c_inc + 2] = 0.; 169 D[p*p_inc + 4*b_inc + 1*c_inc + 0] = y; 170 D[p*p_inc + 4*b_inc + 1*c_inc + 1] = x; 171 D[p*p_inc + 4*b_inc + 1*c_inc + 2] = 0.; 172 D[p*p_inc + 4*b_inc + 2*c_inc + 0] = -3.*z; 173 D[p*p_inc + 4*b_inc + 2*c_inc + 1] = 0.; 174 D[p*p_inc + 4*b_inc + 2*c_inc + 2] = -3.*x; 175 /* {x y, y^2, -3 y z} */ 176 D[p*p_inc + 5*b_inc + 0*c_inc + 0] = y; 177 D[p*p_inc + 5*b_inc + 0*c_inc + 1] = x; 178 D[p*p_inc + 5*b_inc + 0*c_inc + 2] = 0.; 179 D[p*p_inc + 5*b_inc + 1*c_inc + 0] = 0.; 180 D[p*p_inc + 5*b_inc + 1*c_inc + 1] = 2.*y; 181 D[p*p_inc + 5*b_inc + 1*c_inc + 2] = 0.; 182 D[p*p_inc + 5*b_inc + 2*c_inc + 0] = 0.; 183 D[p*p_inc + 5*b_inc + 2*c_inc + 1] = -3.*z; 184 D[p*p_inc + 5*b_inc + 2*c_inc + 2] = -3.*y; 185 } 186 } 187 if (H) { 188 PetscInt p_inc = Nb*Nc*dim*dim; 189 PetscInt b_inc = Nc*dim*dim; 190 PetscInt c_inc = dim*dim; 191 192 for (PetscInt p = 0; p < npoints; p++) { 193 /* {2 y z, 0, 0} */ 194 D[p*p_inc + 0*b_inc + 0*c_inc + 0] = 0.; 195 D[p*p_inc + 0*b_inc + 0*c_inc + 1] = 0.; 196 D[p*p_inc + 0*b_inc + 0*c_inc + 2] = 0.; 197 D[p*p_inc + 0*b_inc + 0*c_inc + 3] = 0.; 198 D[p*p_inc + 0*b_inc + 0*c_inc + 4] = 0.; 199 D[p*p_inc + 0*b_inc + 0*c_inc + 5] = 2.; 200 D[p*p_inc + 0*b_inc + 0*c_inc + 6] = 0.; 201 D[p*p_inc + 0*b_inc + 0*c_inc + 7] = 2.; 202 D[p*p_inc + 0*b_inc + 0*c_inc + 8] = 0.; 203 D[p*p_inc + 0*b_inc + 1*c_inc + 0] = 0.; 204 D[p*p_inc + 0*b_inc + 1*c_inc + 1] = 0.; 205 D[p*p_inc + 0*b_inc + 1*c_inc + 2] = 0.; 206 D[p*p_inc + 0*b_inc + 1*c_inc + 3] = 0.; 207 D[p*p_inc + 0*b_inc + 1*c_inc + 4] = 0.; 208 D[p*p_inc + 0*b_inc + 1*c_inc + 5] = 0.; 209 D[p*p_inc + 0*b_inc + 1*c_inc + 6] = 0.; 210 D[p*p_inc + 0*b_inc + 1*c_inc + 7] = 0.; 211 D[p*p_inc + 0*b_inc + 1*c_inc + 8] = 0.; 212 D[p*p_inc + 0*b_inc + 2*c_inc + 0] = 0.; 213 D[p*p_inc + 0*b_inc + 2*c_inc + 1] = 0.; 214 D[p*p_inc + 0*b_inc + 2*c_inc + 2] = 0.; 215 D[p*p_inc + 0*b_inc + 2*c_inc + 3] = 0.; 216 D[p*p_inc + 0*b_inc + 2*c_inc + 4] = 0.; 217 D[p*p_inc + 0*b_inc + 2*c_inc + 5] = 0.; 218 D[p*p_inc + 0*b_inc + 2*c_inc + 6] = 0.; 219 D[p*p_inc + 0*b_inc + 2*c_inc + 7] = 0.; 220 D[p*p_inc + 0*b_inc + 2*c_inc + 8] = 0.; 221 /* {0, 2 x z, 0} */ 222 D[p*p_inc + 1*b_inc + 0*c_inc + 0] = 0.; 223 D[p*p_inc + 1*b_inc + 0*c_inc + 1] = 0.; 224 D[p*p_inc + 1*b_inc + 0*c_inc + 2] = 0.; 225 D[p*p_inc + 1*b_inc + 0*c_inc + 3] = 0.; 226 D[p*p_inc + 1*b_inc + 0*c_inc + 4] = 0.; 227 D[p*p_inc + 1*b_inc + 0*c_inc + 5] = 0.; 228 D[p*p_inc + 1*b_inc + 0*c_inc + 6] = 0.; 229 D[p*p_inc + 1*b_inc + 0*c_inc + 7] = 0.; 230 D[p*p_inc + 1*b_inc + 0*c_inc + 8] = 0.; 231 D[p*p_inc + 1*b_inc + 1*c_inc + 0] = 0.; 232 D[p*p_inc + 1*b_inc + 1*c_inc + 1] = 0.; 233 D[p*p_inc + 1*b_inc + 1*c_inc + 2] = 2.; 234 D[p*p_inc + 1*b_inc + 1*c_inc + 3] = 0.; 235 D[p*p_inc + 1*b_inc + 1*c_inc + 4] = 0.; 236 D[p*p_inc + 1*b_inc + 1*c_inc + 5] = 0.; 237 D[p*p_inc + 1*b_inc + 1*c_inc + 6] = 2.; 238 D[p*p_inc + 1*b_inc + 1*c_inc + 7] = 0.; 239 D[p*p_inc + 1*b_inc + 1*c_inc + 8] = 0.; 240 D[p*p_inc + 1*b_inc + 2*c_inc + 0] = 0.; 241 D[p*p_inc + 1*b_inc + 2*c_inc + 1] = 0.; 242 D[p*p_inc + 1*b_inc + 2*c_inc + 2] = 0.; 243 D[p*p_inc + 1*b_inc + 2*c_inc + 3] = 0.; 244 D[p*p_inc + 1*b_inc + 2*c_inc + 4] = 0.; 245 D[p*p_inc + 1*b_inc + 2*c_inc + 5] = 0.; 246 D[p*p_inc + 1*b_inc + 2*c_inc + 6] = 0.; 247 D[p*p_inc + 1*b_inc + 2*c_inc + 7] = 0.; 248 D[p*p_inc + 1*b_inc + 2*c_inc + 8] = 0.; 249 /* {0, 2 y z, -z^2} */ 250 D[p*p_inc + 2*b_inc + 0*c_inc + 0] = 0.; 251 D[p*p_inc + 2*b_inc + 0*c_inc + 1] = 0.; 252 D[p*p_inc + 2*b_inc + 0*c_inc + 2] = 0.; 253 D[p*p_inc + 2*b_inc + 0*c_inc + 3] = 0.; 254 D[p*p_inc + 2*b_inc + 0*c_inc + 4] = 0.; 255 D[p*p_inc + 2*b_inc + 0*c_inc + 5] = 0.; 256 D[p*p_inc + 2*b_inc + 0*c_inc + 6] = 0.; 257 D[p*p_inc + 2*b_inc + 0*c_inc + 7] = 0.; 258 D[p*p_inc + 2*b_inc + 0*c_inc + 8] = 0.; 259 D[p*p_inc + 2*b_inc + 1*c_inc + 0] = 0.; 260 D[p*p_inc + 2*b_inc + 1*c_inc + 1] = 0.; 261 D[p*p_inc + 2*b_inc + 1*c_inc + 2] = 0.; 262 D[p*p_inc + 2*b_inc + 1*c_inc + 3] = 0.; 263 D[p*p_inc + 2*b_inc + 1*c_inc + 4] = 0.; 264 D[p*p_inc + 2*b_inc + 1*c_inc + 5] = 2.; 265 D[p*p_inc + 2*b_inc + 1*c_inc + 6] = 0.; 266 D[p*p_inc + 2*b_inc + 1*c_inc + 7] = 2.; 267 D[p*p_inc + 2*b_inc + 1*c_inc + 8] = 0.; 268 D[p*p_inc + 2*b_inc + 2*c_inc + 0] = 0.; 269 D[p*p_inc + 2*b_inc + 2*c_inc + 1] = 0.; 270 D[p*p_inc + 2*b_inc + 2*c_inc + 2] = 0.; 271 D[p*p_inc + 2*b_inc + 2*c_inc + 3] = 0.; 272 D[p*p_inc + 2*b_inc + 2*c_inc + 4] = 0.; 273 D[p*p_inc + 2*b_inc + 2*c_inc + 5] = 0.; 274 D[p*p_inc + 2*b_inc + 2*c_inc + 6] = 0.; 275 D[p*p_inc + 2*b_inc + 2*c_inc + 7] = 0.; 276 D[p*p_inc + 2*b_inc + 2*c_inc + 8] = -2.; 277 /* {2 x z, 0, -z^2} */ 278 D[p*p_inc + 3*b_inc + 0*c_inc + 0] = 0.; 279 D[p*p_inc + 3*b_inc + 0*c_inc + 1] = 0.; 280 D[p*p_inc + 3*b_inc + 0*c_inc + 2] = 2.; 281 D[p*p_inc + 3*b_inc + 0*c_inc + 3] = 0.; 282 D[p*p_inc + 3*b_inc + 0*c_inc + 4] = 0.; 283 D[p*p_inc + 3*b_inc + 0*c_inc + 5] = 0.; 284 D[p*p_inc + 3*b_inc + 0*c_inc + 6] = 2.; 285 D[p*p_inc + 3*b_inc + 0*c_inc + 7] = 0.; 286 D[p*p_inc + 3*b_inc + 0*c_inc + 8] = 0.; 287 D[p*p_inc + 3*b_inc + 1*c_inc + 0] = 0.; 288 D[p*p_inc + 3*b_inc + 1*c_inc + 1] = 0.; 289 D[p*p_inc + 3*b_inc + 1*c_inc + 2] = 0.; 290 D[p*p_inc + 3*b_inc + 1*c_inc + 3] = 0.; 291 D[p*p_inc + 3*b_inc + 1*c_inc + 4] = 0.; 292 D[p*p_inc + 3*b_inc + 1*c_inc + 5] = 0.; 293 D[p*p_inc + 3*b_inc + 1*c_inc + 6] = 0.; 294 D[p*p_inc + 3*b_inc + 1*c_inc + 7] = 0.; 295 D[p*p_inc + 3*b_inc + 1*c_inc + 8] = 0.; 296 D[p*p_inc + 3*b_inc + 2*c_inc + 0] = 0.; 297 D[p*p_inc + 3*b_inc + 2*c_inc + 1] = 0.; 298 D[p*p_inc + 3*b_inc + 2*c_inc + 2] = 0.; 299 D[p*p_inc + 3*b_inc + 2*c_inc + 3] = 0.; 300 D[p*p_inc + 3*b_inc + 2*c_inc + 4] = 0.; 301 D[p*p_inc + 3*b_inc + 2*c_inc + 5] = 0.; 302 D[p*p_inc + 3*b_inc + 2*c_inc + 6] = 0.; 303 D[p*p_inc + 3*b_inc + 2*c_inc + 7] = 0.; 304 D[p*p_inc + 3*b_inc + 2*c_inc + 8] = -2.; 305 /* {x^2, x y, -3 x z} */ 306 D[p*p_inc + 4*b_inc + 0*c_inc + 0] = 2.; 307 D[p*p_inc + 4*b_inc + 0*c_inc + 1] = 0.; 308 D[p*p_inc + 4*b_inc + 0*c_inc + 2] = 0.; 309 D[p*p_inc + 4*b_inc + 0*c_inc + 3] = 0.; 310 D[p*p_inc + 4*b_inc + 0*c_inc + 4] = 0.; 311 D[p*p_inc + 4*b_inc + 0*c_inc + 5] = 0.; 312 D[p*p_inc + 4*b_inc + 0*c_inc + 6] = 0.; 313 D[p*p_inc + 4*b_inc + 0*c_inc + 7] = 0.; 314 D[p*p_inc + 4*b_inc + 0*c_inc + 8] = 0.; 315 D[p*p_inc + 4*b_inc + 1*c_inc + 0] = 0.; 316 D[p*p_inc + 4*b_inc + 1*c_inc + 1] = 1.; 317 D[p*p_inc + 4*b_inc + 1*c_inc + 2] = 0.; 318 D[p*p_inc + 4*b_inc + 1*c_inc + 3] = 1.; 319 D[p*p_inc + 4*b_inc + 1*c_inc + 4] = 0.; 320 D[p*p_inc + 4*b_inc + 1*c_inc + 5] = 0.; 321 D[p*p_inc + 4*b_inc + 1*c_inc + 6] = 0.; 322 D[p*p_inc + 4*b_inc + 1*c_inc + 7] = 0.; 323 D[p*p_inc + 4*b_inc + 1*c_inc + 8] = 0.; 324 D[p*p_inc + 4*b_inc + 2*c_inc + 0] = 0.; 325 D[p*p_inc + 4*b_inc + 2*c_inc + 1] = 0.; 326 D[p*p_inc + 4*b_inc + 2*c_inc + 2] = -3.; 327 D[p*p_inc + 4*b_inc + 2*c_inc + 3] = 0.; 328 D[p*p_inc + 4*b_inc + 2*c_inc + 4] = 0.; 329 D[p*p_inc + 4*b_inc + 2*c_inc + 5] = 0.; 330 D[p*p_inc + 4*b_inc + 2*c_inc + 6] = -3.; 331 D[p*p_inc + 4*b_inc + 2*c_inc + 7] = 0.; 332 D[p*p_inc + 4*b_inc + 2*c_inc + 8] = 0.; 333 /* {x y, y^2, -3 y z} */ 334 D[p*p_inc + 4*b_inc + 0*c_inc + 0] = 2.; 335 D[p*p_inc + 4*b_inc + 0*c_inc + 1] = 0.; 336 D[p*p_inc + 4*b_inc + 0*c_inc + 2] = 0.; 337 D[p*p_inc + 4*b_inc + 0*c_inc + 3] = 0.; 338 D[p*p_inc + 4*b_inc + 0*c_inc + 4] = 0.; 339 D[p*p_inc + 4*b_inc + 0*c_inc + 5] = 0.; 340 D[p*p_inc + 4*b_inc + 0*c_inc + 6] = 0.; 341 D[p*p_inc + 4*b_inc + 0*c_inc + 7] = 0.; 342 D[p*p_inc + 4*b_inc + 0*c_inc + 8] = 0.; 343 D[p*p_inc + 4*b_inc + 1*c_inc + 0] = 0.; 344 D[p*p_inc + 4*b_inc + 1*c_inc + 1] = 1.; 345 D[p*p_inc + 4*b_inc + 1*c_inc + 2] = 0.; 346 D[p*p_inc + 4*b_inc + 1*c_inc + 3] = 1.; 347 D[p*p_inc + 4*b_inc + 1*c_inc + 4] = 0.; 348 D[p*p_inc + 4*b_inc + 1*c_inc + 5] = 0.; 349 D[p*p_inc + 4*b_inc + 1*c_inc + 6] = 0.; 350 D[p*p_inc + 4*b_inc + 1*c_inc + 7] = 0.; 351 D[p*p_inc + 4*b_inc + 1*c_inc + 8] = 0.; 352 D[p*p_inc + 4*b_inc + 2*c_inc + 0] = 0.; 353 D[p*p_inc + 4*b_inc + 2*c_inc + 1] = 0.; 354 D[p*p_inc + 4*b_inc + 2*c_inc + 2] = -3.; 355 D[p*p_inc + 4*b_inc + 2*c_inc + 3] = 0.; 356 D[p*p_inc + 4*b_inc + 2*c_inc + 4] = 0.; 357 D[p*p_inc + 4*b_inc + 2*c_inc + 5] = 0.; 358 D[p*p_inc + 4*b_inc + 2*c_inc + 6] = -3.; 359 D[p*p_inc + 4*b_inc + 2*c_inc + 7] = 0.; 360 D[p*p_inc + 4*b_inc + 2*c_inc + 8] = 0.; 361 } 362 } 363 PetscFunctionReturn(0); 364 } 365 366 static PetscErrorCode PetscSpaceGetHeightSubspace_WXY(PetscSpace sp, PetscInt height, PetscSpace *subsp) 367 { 368 SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Do not know how to do this"); 369 } 370 371 static PetscErrorCode PetscSpaceInitialize_WXY(PetscSpace sp) 372 { 373 PetscFunctionBegin; 374 sp->ops->setfromoptions = PetscSpaceSetFromOptions_WXY; 375 sp->ops->setup = PetscSpaceSetUp_WXY; 376 sp->ops->view = PetscSpaceView_WXY; 377 sp->ops->destroy = PetscSpaceDestroy_WXY; 378 sp->ops->getdimension = PetscSpaceGetDimension_WXY; 379 sp->ops->evaluate = PetscSpaceEvaluate_WXY; 380 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_WXY; 381 PetscFunctionReturn(0); 382 } 383 384 /*MC 385 PETSCSPACEWXY = "wxy" - A PetscSpace object that encapsulates the Wheeler-Xu-Yotov enrichments. 386 $ curl {{0, 0, y^2 z}, {x z^2, 0, 0}, {y z^2, 0, 0}, {0, -x z^2, 0}, {0, -3/2 x^2 z, -1/2 x^2 y}, {3/2 y^2 z, 0, 1/2 y^2 x}} 387 $ = {{2 y z, 0, 0}, {0, 2 x z, 0}, {0, 2 y z, -z^2}, {2 x z, 0, -z^2}, {x^2, x y, -3 x z}, {x y, y^2, -3 y z}} 388 389 Level: intermediate 390 391 .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType() 392 M*/ 393 394 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_WXY(PetscSpace sp) 395 { 396 PetscSpace_WXY *wxy; 397 PetscErrorCode ierr; 398 399 PetscFunctionBegin; 400 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 401 ierr = PetscNewLog(sp,&wxy);CHKERRQ(ierr); 402 sp->data = wxy; 403 sp->degree = 2; 404 405 ierr = PetscSpaceInitialize_WXY(sp);CHKERRQ(ierr); 406 PetscFunctionReturn(0); 407 } 408