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