1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 static PetscErrorCode PetscSpaceTensorCreateSubspace(PetscSpace space, PetscInt Nvs, PetscInt Ncs, PetscSpace *subspace) 4 { 5 PetscInt degree; 6 const char *prefix; 7 const char *name; 8 char subname[PETSC_MAX_PATH_LEN]; 9 10 PetscFunctionBegin; 11 PetscCall(PetscSpaceGetDegree(space, °ree, NULL)); 12 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)space, &prefix)); 13 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)space), subspace)); 14 PetscCall(PetscSpaceSetType(*subspace, PETSCSPACEPOLYNOMIAL)); 15 PetscCall(PetscSpaceSetNumVariables(*subspace, Nvs)); 16 PetscCall(PetscSpaceSetNumComponents(*subspace, Ncs)); 17 PetscCall(PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE)); 18 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix)); 19 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "tensorcomp_")); 20 if (((PetscObject)space)->name) { 21 PetscCall(PetscObjectGetName((PetscObject)space, &name)); 22 PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s tensor component", name)); 23 PetscCall(PetscObjectSetName((PetscObject)*subspace, subname)); 24 } else PetscCall(PetscObjectSetName((PetscObject)*subspace, "tensor component")); 25 PetscFunctionReturn(0); 26 } 27 28 static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscSpace sp, PetscOptionItems *PetscOptionsObject) 29 { 30 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 31 PetscInt Ns, Nc, i, Nv, deg; 32 PetscBool uniform = PETSC_TRUE; 33 34 PetscFunctionBegin; 35 PetscCall(PetscSpaceGetNumVariables(sp, &Nv)); 36 if (!Nv) PetscFunctionReturn(0); 37 PetscCall(PetscSpaceGetNumComponents(sp, &Nc)); 38 PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns)); 39 PetscCall(PetscSpaceGetDegree(sp, °, NULL)); 40 if (Ns > 1) { 41 PetscSpace s0; 42 43 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0)); 44 for (i = 1; i < Ns; i++) { 45 PetscSpace si; 46 47 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si)); 48 if (si != s0) { 49 uniform = PETSC_FALSE; 50 break; 51 } 52 } 53 } 54 Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv, 1) : Ns; 55 PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace tensor options"); 56 PetscCall(PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL, 0)); 57 PetscCall(PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL)); 58 PetscOptionsHeadEnd(); 59 PetscCheck(Ns >= 0 && (Nv <= 0 || Ns != 0), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space made up of %" PetscInt_FMT " spaces", Ns); 60 PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv); 61 if (Ns != tens->numTensSpaces) PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns)); 62 if (uniform) { 63 PetscInt Nvs = Nv / Ns; 64 PetscInt Ncs; 65 PetscSpace subspace; 66 67 PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv); 68 Ncs = (PetscInt)PetscPowReal((PetscReal)Nc, 1. / Ns); 69 PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc); 70 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &subspace)); 71 if (!subspace) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &subspace)); 72 else PetscCall(PetscObjectReference((PetscObject)subspace)); 73 PetscCall(PetscSpaceSetFromOptions(subspace)); 74 for (i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace)); 75 PetscCall(PetscSpaceDestroy(&subspace)); 76 } else { 77 for (i = 0; i < Ns; i++) { 78 PetscSpace subspace; 79 80 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &subspace)); 81 if (!subspace) { 82 char tprefix[128]; 83 84 PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &subspace)); 85 PetscCall(PetscSNPrintf(tprefix, 128, "%d_", (int)i)); 86 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix)); 87 } else PetscCall(PetscObjectReference((PetscObject)subspace)); 88 PetscCall(PetscSpaceSetFromOptions(subspace)); 89 PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace)); 90 PetscCall(PetscSpaceDestroy(&subspace)); 91 } 92 } 93 PetscFunctionReturn(0); 94 } 95 96 static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v) 97 { 98 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 99 PetscBool uniform = PETSC_TRUE; 100 PetscInt Ns = tens->numTensSpaces, i, n; 101 102 PetscFunctionBegin; 103 for (i = 1; i < Ns; i++) { 104 if (tens->tensspaces[i] != tens->tensspaces[0]) { 105 uniform = PETSC_FALSE; 106 break; 107 } 108 } 109 if (uniform) PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces (all identical)\n", Ns)); 110 else PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces\n", Ns)); 111 n = uniform ? 1 : Ns; 112 for (i = 0; i < n; i++) { 113 PetscCall(PetscViewerASCIIPushTab(v)); 114 PetscCall(PetscSpaceView(tens->tensspaces[i], v)); 115 PetscCall(PetscViewerASCIIPopTab(v)); 116 } 117 PetscFunctionReturn(0); 118 } 119 120 static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer) 121 { 122 PetscBool iascii; 123 124 PetscFunctionBegin; 125 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 126 if (iascii) PetscCall(PetscSpaceTensorView_Ascii(sp, viewer)); 127 PetscFunctionReturn(0); 128 } 129 130 static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp) 131 { 132 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 133 PetscInt Nc, Nv, Ns; 134 PetscBool uniform = PETSC_TRUE; 135 PetscInt deg, maxDeg; 136 PetscInt Ncprod; 137 138 PetscFunctionBegin; 139 if (tens->setupCalled) PetscFunctionReturn(0); 140 PetscCall(PetscSpaceGetNumVariables(sp, &Nv)); 141 PetscCall(PetscSpaceGetNumComponents(sp, &Nc)); 142 PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns)); 143 if (Ns == PETSC_DEFAULT) { 144 Ns = Nv; 145 PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns)); 146 } 147 if (!Ns) { 148 SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces"); 149 } else { 150 PetscSpace s0; 151 152 PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv); 153 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0)); 154 for (PetscInt i = 1; i < Ns; i++) { 155 PetscSpace si; 156 157 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si)); 158 if (si != s0) { 159 uniform = PETSC_FALSE; 160 break; 161 } 162 } 163 if (uniform) { 164 PetscInt Nvs = Nv / Ns; 165 PetscInt Ncs; 166 167 PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv); 168 Ncs = (PetscInt)(PetscPowReal((PetscReal)Nc, 1. / Ns)); 169 PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc); 170 if (!s0) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &s0)); 171 else PetscCall(PetscObjectReference((PetscObject)s0)); 172 PetscCall(PetscSpaceSetUp(s0)); 173 for (PetscInt i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, s0)); 174 PetscCall(PetscSpaceDestroy(&s0)); 175 Ncprod = PetscPowInt(Ncs, Ns); 176 } else { 177 PetscInt Nvsum = 0; 178 179 Ncprod = 1; 180 for (PetscInt i = 0; i < Ns; i++) { 181 PetscInt Nvs, Ncs; 182 PetscSpace si; 183 184 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si)); 185 if (!si) PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &si)); 186 else PetscCall(PetscObjectReference((PetscObject)si)); 187 PetscCall(PetscSpaceSetUp(si)); 188 PetscCall(PetscSpaceTensorSetSubspace(sp, i, si)); 189 PetscCall(PetscSpaceGetNumVariables(si, &Nvs)); 190 PetscCall(PetscSpaceGetNumComponents(si, &Ncs)); 191 PetscCall(PetscSpaceDestroy(&si)); 192 193 Nvsum += Nvs; 194 Ncprod *= Ncs; 195 } 196 197 PetscCheck(Nvsum == Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Sum of subspace variables %" PetscInt_FMT " does not equal the number of variables %" PetscInt_FMT, Nvsum, Nv); 198 PetscCheck(Nc % Ncprod == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Product of subspace components %" PetscInt_FMT " does not divide the number of components %" PetscInt_FMT, Ncprod, Nc); 199 } 200 if (Ncprod != Nc) { 201 PetscInt Ncopies = Nc / Ncprod; 202 PetscInt Nv = sp->Nv; 203 const char *prefix; 204 const char *name; 205 char subname[PETSC_MAX_PATH_LEN]; 206 PetscSpace subsp; 207 208 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp)); 209 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix)); 210 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix)); 211 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_")); 212 if (((PetscObject)sp)->name) { 213 PetscCall(PetscObjectGetName((PetscObject)sp, &name)); 214 PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name)); 215 PetscCall(PetscObjectSetName((PetscObject)subsp, subname)); 216 } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component")); 217 PetscCall(PetscSpaceSetType(subsp, PETSCSPACETENSOR)); 218 PetscCall(PetscSpaceSetNumVariables(subsp, Nv)); 219 PetscCall(PetscSpaceSetNumComponents(subsp, Ncprod)); 220 PetscCall(PetscSpaceTensorSetNumSubspaces(subsp, Ns)); 221 for (PetscInt i = 0; i < Ns; i++) { 222 PetscSpace ssp; 223 224 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &ssp)); 225 PetscCall(PetscSpaceTensorSetSubspace(subsp, i, ssp)); 226 } 227 PetscCall(PetscSpaceSetUp(subsp)); 228 PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM)); 229 PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ncopies)); 230 for (PetscInt i = 0; i < Ncopies; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp)); 231 PetscCall(PetscSpaceDestroy(&subsp)); 232 PetscCall(PetscSpaceSetUp(sp)); 233 PetscFunctionReturn(0); 234 } 235 } 236 deg = PETSC_MAX_INT; 237 maxDeg = 0; 238 for (PetscInt i = 0; i < Ns; i++) { 239 PetscSpace si; 240 PetscInt iDeg, iMaxDeg; 241 242 PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si)); 243 PetscCall(PetscSpaceGetDegree(si, &iDeg, &iMaxDeg)); 244 deg = PetscMin(deg, iDeg); 245 maxDeg += iMaxDeg; 246 } 247 sp->degree = deg; 248 sp->maxDegree = maxDeg; 249 tens->uniform = uniform; 250 tens->setupCalled = PETSC_TRUE; 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp) 255 { 256 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 257 PetscInt Ns, i; 258 259 PetscFunctionBegin; 260 Ns = tens->numTensSpaces; 261 if (tens->heightsubspaces) { 262 PetscInt d; 263 264 /* sp->Nv is the spatial dimension, so it is equal to the number 265 * of subspaces on higher co-dimension points */ 266 for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&tens->heightsubspaces[d])); 267 } 268 PetscCall(PetscFree(tens->heightsubspaces)); 269 for (i = 0; i < Ns; i++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[i])); 270 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", NULL)); 271 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", NULL)); 272 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", NULL)); 273 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", NULL)); 274 PetscCall(PetscFree(tens->tensspaces)); 275 PetscCall(PetscFree(tens)); 276 PetscFunctionReturn(0); 277 } 278 279 static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim) 280 { 281 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 282 PetscInt i, Ns, d; 283 284 PetscFunctionBegin; 285 PetscCall(PetscSpaceSetUp(sp)); 286 Ns = tens->numTensSpaces; 287 d = 1; 288 for (i = 0; i < Ns; i++) { 289 PetscInt id; 290 291 PetscCall(PetscSpaceGetDimension(tens->tensspaces[i], &id)); 292 d *= id; 293 } 294 *dim = d; 295 PetscFunctionReturn(0); 296 } 297 298 static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 299 { 300 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 301 DM dm = sp->dm; 302 PetscInt Nc = sp->Nc; 303 PetscInt Nv = sp->Nv; 304 PetscInt Ns; 305 PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL; 306 PetscInt pdim; 307 308 PetscFunctionBegin; 309 if (!tens->setupCalled) { 310 PetscCall(PetscSpaceSetUp(sp)); 311 PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H)); 312 PetscFunctionReturn(0); 313 } 314 Ns = tens->numTensSpaces; 315 PetscCall(PetscSpaceGetDimension(sp, &pdim)); 316 PetscCall(DMGetWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints)); 317 if (B || D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB)); 318 if (D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD)); 319 if (H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH)); 320 if (B) { 321 for (PetscInt i = 0; i < npoints * pdim * Nc; i++) B[i] = 1.; 322 } 323 if (D) { 324 for (PetscInt i = 0; i < npoints * pdim * Nc * Nv; i++) D[i] = 1.; 325 } 326 if (H) { 327 for (PetscInt i = 0; i < npoints * pdim * Nc * Nv * Nv; i++) H[i] = 1.; 328 } 329 for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) { 330 PetscInt sNv, sNc, spdim; 331 PetscInt vskip, cskip; 332 333 PetscCall(PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv)); 334 PetscCall(PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc)); 335 PetscCall(PetscSpaceGetDimension(tens->tensspaces[s], &spdim)); 336 PetscCheck(!(pdim % vstep) && !(pdim % spdim), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", pdim %" PetscInt_FMT ", s %" PetscInt_FMT ", vstep %" PetscInt_FMT ", spdim %" PetscInt_FMT, Nv, Ns, pdim, s, vstep, spdim); 337 PetscCheck(!(Nc % cstep) && !(Nc % sNc), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", Nc %" PetscInt_FMT ", s %" PetscInt_FMT ", cstep %" PetscInt_FMT ", sNc %" PetscInt_FMT, Nv, Ns, Nc, s, cstep, spdim); 338 vskip = pdim / (vstep * spdim); 339 cskip = Nc / (cstep * sNc); 340 for (PetscInt p = 0; p < npoints; ++p) { 341 for (PetscInt i = 0; i < sNv; i++) lpoints[p * sNv + i] = points[p * Nv + d + i]; 342 } 343 PetscCall(PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH)); 344 if (B) { 345 for (PetscInt k = 0; k < vskip; k++) { 346 for (PetscInt si = 0; si < spdim; si++) { 347 for (PetscInt j = 0; j < vstep; j++) { 348 PetscInt i = (k * spdim + si) * vstep + j; 349 350 for (PetscInt l = 0; l < cskip; l++) { 351 for (PetscInt sc = 0; sc < sNc; sc++) { 352 for (PetscInt m = 0; m < cstep; m++) { 353 PetscInt c = (l * sNc + sc) * cstep + m; 354 355 for (PetscInt p = 0; p < npoints; p++) B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc]; 356 } 357 } 358 } 359 } 360 } 361 } 362 } 363 if (D) { 364 for (PetscInt k = 0; k < vskip; k++) { 365 for (PetscInt si = 0; si < spdim; si++) { 366 for (PetscInt j = 0; j < vstep; j++) { 367 PetscInt i = (k * spdim + si) * vstep + j; 368 369 for (PetscInt l = 0; l < cskip; l++) { 370 for (PetscInt sc = 0; sc < sNc; sc++) { 371 for (PetscInt m = 0; m < cstep; m++) { 372 PetscInt c = (l * sNc + sc) * cstep + m; 373 374 for (PetscInt der = 0; der < Nv; der++) { 375 if (der >= d && der < d + sNv) { 376 for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d]; 377 } else { 378 for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sB[(spdim * p + si) * sNc + sc]; 379 } 380 } 381 } 382 } 383 } 384 } 385 } 386 } 387 } 388 if (H) { 389 for (PetscInt k = 0; k < vskip; k++) { 390 for (PetscInt si = 0; si < spdim; si++) { 391 for (PetscInt j = 0; j < vstep; j++) { 392 PetscInt i = (k * spdim + si) * vstep + j; 393 394 for (PetscInt l = 0; l < cskip; l++) { 395 for (PetscInt sc = 0; sc < sNc; sc++) { 396 for (PetscInt m = 0; m < cstep; m++) { 397 PetscInt c = (l * sNc + sc) * cstep + m; 398 399 for (PetscInt der = 0; der < Nv; der++) { 400 for (PetscInt der2 = 0; der2 < Nv; der2++) { 401 if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) { 402 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sH[(((spdim * p + si) * sNc + sc) * sNv + der - d) * sNv + der2 - d]; 403 } else if (der >= d && der < d + sNv) { 404 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d]; 405 } else if (der2 >= d && der2 < d + sNv) { 406 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der2 - d]; 407 } else { 408 for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc]; 409 } 410 } 411 } 412 } 413 } 414 } 415 } 416 } 417 } 418 } 419 d += sNv; 420 vstep *= spdim; 421 cstep *= sNc; 422 } 423 if (H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH)); 424 if (D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD)); 425 if (B || D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB)); 426 PetscCall(DMRestoreWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints)); 427 PetscFunctionReturn(0); 428 } 429 430 /*@ 431 PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product 432 433 Input Parameters: 434 + sp - the function space object 435 - numTensSpaces - the number of spaces 436 437 Level: intermediate 438 439 .seealso: `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 440 @*/ 441 PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces) 442 { 443 PetscFunctionBegin; 444 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 445 PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces)); 446 PetscFunctionReturn(0); 447 } 448 449 /*@ 450 PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product 451 452 Input Parameter: 453 . sp - the function space object 454 455 Output Parameter: 456 . numTensSpaces - the number of spaces 457 458 Level: intermediate 459 460 .seealso: `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 461 @*/ 462 PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces) 463 { 464 PetscFunctionBegin; 465 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 466 PetscValidIntPointer(numTensSpaces, 2); 467 PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces)); 468 PetscFunctionReturn(0); 469 } 470 471 /*@ 472 PetscSpaceTensorSetSubspace - Set a space in the tensor product 473 474 Input Parameters: 475 + sp - the function space object 476 . s - The space number 477 - subsp - the number of spaces 478 479 Level: intermediate 480 481 .seealso: `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 482 @*/ 483 PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp) 484 { 485 PetscFunctionBegin; 486 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 487 if (subsp) PetscValidHeaderSpecific(subsp, PETSCSPACE_CLASSID, 3); 488 PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp)); 489 PetscFunctionReturn(0); 490 } 491 492 /*@ 493 PetscSpaceTensorGetSubspace - Get a space in the tensor product 494 495 Input Parameters: 496 + sp - the function space object 497 - s - The space number 498 499 Output Parameter: 500 . subsp - the PetscSpace 501 502 Level: intermediate 503 504 .seealso: `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 505 @*/ 506 PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp) 507 { 508 PetscFunctionBegin; 509 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 510 PetscValidPointer(subsp, 3); 511 PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp)); 512 PetscFunctionReturn(0); 513 } 514 515 static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces) 516 { 517 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 518 PetscInt Ns; 519 520 PetscFunctionBegin; 521 PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called"); 522 Ns = tens->numTensSpaces; 523 if (numTensSpaces == Ns) PetscFunctionReturn(0); 524 if (Ns >= 0) { 525 PetscInt s; 526 527 for (s = 0; s < Ns; s++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[s])); 528 PetscCall(PetscFree(tens->tensspaces)); 529 } 530 Ns = tens->numTensSpaces = numTensSpaces; 531 PetscCall(PetscCalloc1(Ns, &tens->tensspaces)); 532 PetscFunctionReturn(0); 533 } 534 535 static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces) 536 { 537 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 538 539 PetscFunctionBegin; 540 *numTensSpaces = tens->numTensSpaces; 541 PetscFunctionReturn(0); 542 } 543 544 static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace) 545 { 546 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 547 PetscInt Ns; 548 549 PetscFunctionBegin; 550 PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called"); 551 Ns = tens->numTensSpaces; 552 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first"); 553 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s); 554 PetscCall(PetscObjectReference((PetscObject)subspace)); 555 PetscCall(PetscSpaceDestroy(&tens->tensspaces[s])); 556 tens->tensspaces[s] = subspace; 557 PetscFunctionReturn(0); 558 } 559 560 static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp) 561 { 562 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 563 PetscInt Nc, dim, order, i; 564 PetscSpace bsp; 565 566 PetscFunctionBegin; 567 PetscCall(PetscSpaceGetNumVariables(sp, &dim)); 568 PetscCheck(tens->uniform && tens->numTensSpaces == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can only get a generic height subspace of a uniform tensor space of 1d spaces."); 569 PetscCall(PetscSpaceGetNumComponents(sp, &Nc)); 570 PetscCall(PetscSpaceGetDegree(sp, &order, NULL)); 571 PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim); 572 if (!tens->heightsubspaces) PetscCall(PetscCalloc1(dim, &tens->heightsubspaces)); 573 if (height <= dim) { 574 if (!tens->heightsubspaces[height - 1]) { 575 PetscSpace sub; 576 const char *name; 577 578 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &bsp)); 579 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub)); 580 PetscCall(PetscObjectGetName((PetscObject)sp, &name)); 581 PetscCall(PetscObjectSetName((PetscObject)sub, name)); 582 PetscCall(PetscSpaceSetType(sub, PETSCSPACETENSOR)); 583 PetscCall(PetscSpaceSetNumComponents(sub, Nc)); 584 PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE)); 585 PetscCall(PetscSpaceSetNumVariables(sub, dim - height)); 586 PetscCall(PetscSpaceTensorSetNumSubspaces(sub, dim - height)); 587 for (i = 0; i < dim - height; i++) PetscCall(PetscSpaceTensorSetSubspace(sub, i, bsp)); 588 PetscCall(PetscSpaceSetUp(sub)); 589 tens->heightsubspaces[height - 1] = sub; 590 } 591 *subsp = tens->heightsubspaces[height - 1]; 592 } else { 593 *subsp = NULL; 594 } 595 PetscFunctionReturn(0); 596 } 597 598 static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace) 599 { 600 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 601 PetscInt Ns; 602 603 PetscFunctionBegin; 604 Ns = tens->numTensSpaces; 605 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first"); 606 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s); 607 *subspace = tens->tensspaces[s]; 608 PetscFunctionReturn(0); 609 } 610 611 static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp) 612 { 613 PetscFunctionBegin; 614 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor; 615 sp->ops->setup = PetscSpaceSetUp_Tensor; 616 sp->ops->view = PetscSpaceView_Tensor; 617 sp->ops->destroy = PetscSpaceDestroy_Tensor; 618 sp->ops->getdimension = PetscSpaceGetDimension_Tensor; 619 sp->ops->evaluate = PetscSpaceEvaluate_Tensor; 620 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor; 621 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor)); 622 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor)); 623 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor)); 624 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor)); 625 PetscFunctionReturn(0); 626 } 627 628 /*MC 629 PETSCSPACETENSOR = "tensor" - A PetscSpace object that encapsulates a tensor product space. 630 A tensor product is created of the components of the subspaces as well. 631 632 Level: intermediate 633 634 .seealso: `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()` 635 M*/ 636 637 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp) 638 { 639 PetscSpace_Tensor *tens; 640 641 PetscFunctionBegin; 642 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 643 PetscCall(PetscNew(&tens)); 644 sp->data = tens; 645 646 tens->numTensSpaces = PETSC_DEFAULT; 647 648 PetscCall(PetscSpaceInitialize_Tensor(sp)); 649 PetscFunctionReturn(0); 650 } 651