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