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