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