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