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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 = NULL; 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 = NULL; 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 428 } 429 430 /*@ 431 PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product space 432 433 Input Parameters: 434 + sp - the function space object 435 - numTensSpaces - the number of spaces 436 437 Level: intermediate 438 439 Note: 440 The name NumSubspaces is misleading because it is actually setting the number of defining spaces of the tensor product space, not a number of Subspaces of it 441 442 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 443 @*/ 444 PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces) 445 { 446 PetscFunctionBegin; 447 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 448 PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces)); 449 PetscFunctionReturn(PETSC_SUCCESS); 450 } 451 452 /*@ 453 PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product space 454 455 Input Parameter: 456 . sp - the function space object 457 458 Output Parameter: 459 . numTensSpaces - the number of spaces 460 461 Level: intermediate 462 463 Note: 464 The name NumSubspaces is misleading because it is actually getting the number of defining spaces of the tensor product space, not a number of Subspaces of it 465 466 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 467 @*/ 468 PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces) 469 { 470 PetscFunctionBegin; 471 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 472 PetscAssertPointer(numTensSpaces, 2); 473 PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces)); 474 PetscFunctionReturn(PETSC_SUCCESS); 475 } 476 477 /*@ 478 PetscSpaceTensorSetSubspace - Set a space in the tensor product space 479 480 Input Parameters: 481 + sp - the function space object 482 . s - The space number 483 - subsp - the number of spaces 484 485 Level: intermediate 486 487 Note: 488 The name SetSubspace is misleading because it is actually setting one of the defining spaces of the tensor product space, not a Subspace of it 489 490 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 491 @*/ 492 PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp) 493 { 494 PetscFunctionBegin; 495 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 496 if (subsp) PetscValidHeaderSpecific(subsp, PETSCSPACE_CLASSID, 3); 497 PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp)); 498 PetscFunctionReturn(PETSC_SUCCESS); 499 } 500 501 /*@ 502 PetscSpaceTensorGetSubspace - Get a space in the tensor product space 503 504 Input Parameters: 505 + sp - the function space object 506 - s - The space number 507 508 Output Parameter: 509 . subsp - the `PetscSpace` 510 511 Level: intermediate 512 513 Note: 514 The name GetSubspace is misleading because it is actually getting one of the defining spaces of the tensor product space, not a Subspace of it 515 516 .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()` 517 @*/ 518 PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp) 519 { 520 PetscFunctionBegin; 521 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 522 PetscAssertPointer(subsp, 3); 523 PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp)); 524 PetscFunctionReturn(PETSC_SUCCESS); 525 } 526 527 static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces) 528 { 529 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 530 PetscInt Ns; 531 532 PetscFunctionBegin; 533 PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called"); 534 Ns = tens->numTensSpaces; 535 if (numTensSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS); 536 if (Ns >= 0) { 537 PetscInt s; 538 539 for (s = 0; s < Ns; s++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[s])); 540 PetscCall(PetscFree(tens->tensspaces)); 541 } 542 Ns = tens->numTensSpaces = numTensSpaces; 543 PetscCall(PetscCalloc1(Ns, &tens->tensspaces)); 544 PetscFunctionReturn(PETSC_SUCCESS); 545 } 546 547 static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces) 548 { 549 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 550 551 PetscFunctionBegin; 552 *numTensSpaces = tens->numTensSpaces; 553 PetscFunctionReturn(PETSC_SUCCESS); 554 } 555 556 static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace) 557 { 558 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 559 PetscInt Ns; 560 561 PetscFunctionBegin; 562 PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called"); 563 Ns = tens->numTensSpaces; 564 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first"); 565 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s); 566 PetscCall(PetscObjectReference((PetscObject)subspace)); 567 PetscCall(PetscSpaceDestroy(&tens->tensspaces[s])); 568 tens->tensspaces[s] = subspace; 569 PetscFunctionReturn(PETSC_SUCCESS); 570 } 571 572 static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp) 573 { 574 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data; 575 PetscInt Nc, dim, order, i; 576 PetscSpace bsp; 577 578 PetscFunctionBegin; 579 PetscCall(PetscSpaceGetNumVariables(sp, &dim)); 580 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."); 581 PetscCall(PetscSpaceGetNumComponents(sp, &Nc)); 582 PetscCall(PetscSpaceGetDegree(sp, &order, NULL)); 583 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); 584 if (!tens->heightsubspaces) PetscCall(PetscCalloc1(dim, &tens->heightsubspaces)); 585 if (height <= dim) { 586 if (!tens->heightsubspaces[height - 1]) { 587 PetscSpace sub; 588 const char *name; 589 590 PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &bsp)); 591 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub)); 592 PetscCall(PetscObjectGetName((PetscObject)sp, &name)); 593 PetscCall(PetscObjectSetName((PetscObject)sub, name)); 594 PetscCall(PetscSpaceSetType(sub, PETSCSPACETENSOR)); 595 PetscCall(PetscSpaceSetNumComponents(sub, Nc)); 596 PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE)); 597 PetscCall(PetscSpaceSetNumVariables(sub, dim - height)); 598 PetscCall(PetscSpaceTensorSetNumSubspaces(sub, dim - height)); 599 for (i = 0; i < dim - height; i++) PetscCall(PetscSpaceTensorSetSubspace(sub, i, bsp)); 600 PetscCall(PetscSpaceSetUp(sub)); 601 tens->heightsubspaces[height - 1] = sub; 602 } 603 *subsp = tens->heightsubspaces[height - 1]; 604 } else { 605 *subsp = NULL; 606 } 607 PetscFunctionReturn(PETSC_SUCCESS); 608 } 609 610 static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace) 611 { 612 PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data; 613 PetscInt Ns; 614 615 PetscFunctionBegin; 616 Ns = tens->numTensSpaces; 617 PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first"); 618 PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s); 619 *subspace = tens->tensspaces[s]; 620 PetscFunctionReturn(PETSC_SUCCESS); 621 } 622 623 static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp) 624 { 625 PetscFunctionBegin; 626 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor; 627 sp->ops->setup = PetscSpaceSetUp_Tensor; 628 sp->ops->view = PetscSpaceView_Tensor; 629 sp->ops->destroy = PetscSpaceDestroy_Tensor; 630 sp->ops->getdimension = PetscSpaceGetDimension_Tensor; 631 sp->ops->evaluate = PetscSpaceEvaluate_Tensor; 632 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor; 633 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor)); 634 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor)); 635 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor)); 636 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor)); 637 PetscFunctionReturn(PETSC_SUCCESS); 638 } 639 640 /*MC 641 PETSCSPACETENSOR = "tensor" - A `PetscSpace` object that encapsulates a tensor product space. 642 A tensor product is created of the components of the subspaces as well. 643 644 Level: intermediate 645 646 .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceTensorSetSubspace()`, 647 `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceTensorSetNumSubspaces()` 648 M*/ 649 650 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp) 651 { 652 PetscSpace_Tensor *tens; 653 654 PetscFunctionBegin; 655 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 656 PetscCall(PetscNew(&tens)); 657 sp->data = tens; 658 659 tens->numTensSpaces = PETSC_DEFAULT; 660 661 PetscCall(PetscSpaceInitialize_Tensor(sp)); 662 PetscFunctionReturn(PETSC_SUCCESS); 663 } 664