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