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);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 = PetscOptionsInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL);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->numSpaces) {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 viewer) 88 { 89 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data; 90 PetscBool uniform = PETSC_TRUE; 91 PetscInt Ns = tens->numSpaces, i, n; 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 for (i = 1; i < Ns; i++) { 96 if (tens->spaces[i] != tens->spaces[0]) {uniform = PETSC_FALSE; break;} 97 } 98 if (uniform) {ierr = PetscViewerASCIIPrintf(viewer, "Tensor space of %D subspaces (all identical)\n", Ns);CHKERRQ(ierr); 99 } else {ierr = PetscViewerASCIIPrintf(viewer, "Tensor space of %D subspaces\n", Ns);CHKERRQ(ierr);} 100 n = uniform ? 1 : Ns; 101 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 102 for (i = 0; i < n; i++) { 103 ierr = PetscSpaceView(tens->spaces[i], viewer);CHKERRQ(ierr); 104 } 105 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 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 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 116 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 117 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 118 if (iascii) {ierr = PetscSpaceTensorView_Ascii(sp, viewer);CHKERRQ(ierr);} 119 PetscFunctionReturn(0); 120 } 121 122 static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp) 123 { 124 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data; 125 PetscInt Nv, Ns, i; 126 PetscBool uniform = PETSC_TRUE; 127 PetscInt deg, maxDeg; 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 if (tens->setupCalled) PetscFunctionReturn(0); 132 ierr = PetscSpaceGetNumVariables(sp, &Nv);CHKERRQ(ierr); 133 ierr = PetscSpaceTensorGetNumSubspaces(sp, &Ns);CHKERRQ(ierr); 134 if (Ns == PETSC_DEFAULT) { 135 Ns = Nv; 136 ierr = PetscSpaceTensorSetNumSubspaces(sp, Ns);CHKERRQ(ierr); 137 } 138 if (!Ns) { 139 if (Nv) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces"); 140 } else { 141 PetscSpace s0; 142 143 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); 144 ierr = PetscSpaceTensorGetSubspace(sp, 0, &s0);CHKERRQ(ierr); 145 for (i = 1; i < Ns; i++) { 146 PetscSpace si; 147 148 ierr = PetscSpaceTensorGetSubspace(sp, i, &si);CHKERRQ(ierr); 149 if (si != s0) {uniform = PETSC_FALSE; break;} 150 } 151 if (uniform) { 152 PetscInt Nvs = Nv / Ns; 153 154 if (Nv % Ns) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONG,"Cannot use %D uniform subspaces for %D variable space\n", Ns, Nv); 155 if (!s0) {ierr = PetscSpaceTensorCreateSubspace(sp, Nvs, &s0);CHKERRQ(ierr);} 156 else {ierr = PetscObjectReference((PetscObject) s0);CHKERRQ(ierr);} 157 ierr = PetscSpaceSetUp(s0);CHKERRQ(ierr); 158 for (i = 0; i < Ns; i++) {ierr = PetscSpaceTensorSetSubspace(sp, i, s0);CHKERRQ(ierr);} 159 ierr = PetscSpaceDestroy(&s0);CHKERRQ(ierr); 160 } else { 161 for (i = 0 ; i < Ns; i++) { 162 PetscSpace si; 163 164 ierr = PetscSpaceTensorGetSubspace(sp, i, &si);CHKERRQ(ierr); 165 if (!si) {ierr = PetscSpaceTensorCreateSubspace(sp, 1, &si);CHKERRQ(ierr);} 166 else {ierr = PetscObjectReference((PetscObject) si);CHKERRQ(ierr);} 167 ierr = PetscSpaceSetUp(si);CHKERRQ(ierr); 168 ierr = PetscSpaceTensorSetSubspace(sp, i, si);CHKERRQ(ierr); 169 ierr = PetscSpaceDestroy(&si);CHKERRQ(ierr); 170 } 171 } 172 } 173 deg = PETSC_MAX_INT; 174 maxDeg = 0; 175 for (i = 0; i < Ns; i++) { 176 PetscSpace si; 177 PetscInt iDeg, iMaxDeg; 178 179 ierr = PetscSpaceTensorGetSubspace(sp, i, &si);CHKERRQ(ierr); 180 ierr = PetscSpaceGetDegree(si, &iDeg, &iMaxDeg);CHKERRQ(ierr); 181 deg = PetscMin(deg, iDeg); 182 maxDeg += iMaxDeg; 183 } 184 sp->degree = deg; 185 sp->maxDegree = maxDeg; 186 tens->uniform = uniform; 187 tens->setupCalled = PETSC_TRUE; 188 PetscFunctionReturn(0); 189 } 190 191 static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp) 192 { 193 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data; 194 PetscInt Ns, i; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 ierr = PetscSpaceTensorGetNumSubspaces(sp, &Ns);CHKERRQ(ierr); 199 for (i = 0; i < Ns; i++) {ierr = PetscSpaceDestroy(&tens->spaces[i]);CHKERRQ(ierr);} 200 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", NULL);CHKERRQ(ierr); 201 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", NULL);CHKERRQ(ierr); 202 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", NULL);CHKERRQ(ierr); 203 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", NULL);CHKERRQ(ierr); 204 ierr = PetscFree(tens->spaces);CHKERRQ(ierr); 205 ierr = PetscFree(tens);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim) 210 { 211 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data; 212 PetscInt i, Ns, Nc, d; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 if (tens->dim == PETSC_DEFAULT) { 217 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 218 Ns = tens->numSpaces; 219 Nc = sp->Nc; 220 d = 1; 221 for (i = 0; i < Ns; i++) { 222 PetscInt id; 223 224 ierr = PetscSpaceGetDimension(tens->spaces[i], &id);CHKERRQ(ierr); 225 d *= id; 226 } 227 d *= Nc; 228 tens->dim = d; 229 } 230 *dim = tens->dim; 231 PetscFunctionReturn(0); 232 } 233 234 static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 235 { 236 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data; 237 DM dm = sp->dm; 238 PetscInt Nc = sp->Nc; 239 PetscInt Nv = sp->Nv; 240 PetscInt Ns; 241 PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL; 242 PetscInt c, pdim, d, e, der, der2, i, l, si, p, s, step; 243 PetscErrorCode ierr; 244 245 PetscFunctionBegin; 246 if (!tens->setupCalled) {ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr);} 247 Ns = tens->numSpaces; 248 ierr = PetscSpaceGetDimension(sp,&pdim);CHKERRQ(ierr); 249 pdim /= Nc; 250 ierr = DMGetWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);CHKERRQ(ierr); 251 if (B || D || H) {ierr = DMGetWorkArray(dm, npoints*pdim, MPIU_REAL, &sB);CHKERRQ(ierr);} 252 if (D || H) {ierr = DMGetWorkArray(dm, npoints*pdim*Nv, MPIU_REAL, &sD);CHKERRQ(ierr);} 253 if (H) {ierr = DMGetWorkArray(dm, npoints*pdim*Nv*Nv, MPIU_REAL, &sH);CHKERRQ(ierr);} 254 if (B) { 255 for (i = 0; i < npoints*pdim; i++) B[i * Nc*Nc] = 1.; 256 } 257 if (D) { 258 for (i = 0; i < npoints*pdim; i++) { 259 for (l = 0; l < Nv; l++) { 260 D[i * Nc*Nc*Nv + l] = 1.; 261 } 262 } 263 } 264 if (H) { 265 for (i = 0; i < npoints*pdim; i++) { 266 for (l = 0; l < Nv*Nv; l++) { 267 H[i * Nc*Nc*Nv*Nv + l] = 1.; 268 } 269 } 270 } 271 for (s = 0, d = 0, step = 1; s < Ns; s++) { 272 PetscInt sNv, spdim; 273 PetscInt skip, j, k; 274 275 ierr = PetscSpaceGetNumVariables(tens->spaces[s], &sNv);CHKERRQ(ierr); 276 ierr = PetscSpaceGetDimension(tens->spaces[s], &spdim);CHKERRQ(ierr); 277 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); 278 skip = pdim / (step * spdim); 279 for (p = 0; p < npoints; ++p) { 280 for (i = 0; i < sNv; i++) { 281 lpoints[p * sNv + i] = points[p*Nv + d + i]; 282 } 283 } 284 ierr = PetscSpaceEvaluate(tens->spaces[s], npoints, lpoints, sB, sD, sH);CHKERRQ(ierr); 285 if (B) { 286 for (p = 0; p < npoints; p++) { 287 for (k = 0; k < skip; k++) { 288 for (si = 0; si < spdim; si++) { 289 for (j = 0; j < step; j++) { 290 i = (k * spdim + si) * step + j; 291 B[(pdim * p + i) * Nc * Nc] *= sB[spdim * p + si]; 292 } 293 } 294 } 295 } 296 } 297 if (D) { 298 for (p = 0; p < npoints; p++) { 299 for (k = 0; k < skip; k++) { 300 for (si = 0; si < spdim; si++) { 301 for (j = 0; j < step; j++) { 302 i = (k * spdim + si) * step + j; 303 for (der = 0; der < Nv; der++) { 304 if (der >= d && der < d + sNv) { 305 D[(pdim * p + i) * Nc*Nc*Nv + der] *= sD[(spdim * p + si) * sNv + der - d]; 306 } else { 307 D[(pdim * p + i) * Nc*Nc*Nv + der] *= sB[spdim * p + si]; 308 } 309 } 310 } 311 } 312 } 313 } 314 } 315 if (H) { 316 for (p = 0; p < npoints; p++) { 317 for (k = 0; k < skip; k++) { 318 for (si = 0; si < spdim; si++) { 319 for (j = 0; j < step; j++) { 320 i = (k * spdim + si) * step + j; 321 for (der = 0; der < Nv; der++) { 322 for (der2 = 0; der2 < Nv; der2++) { 323 if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) { 324 H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sH[((spdim * p + si) * sNv + der - d) * sNv + der2 - d]; 325 } else if (der >= d && der < d + sNv) { 326 H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sD[(spdim * p + si) * sNv + der - d]; 327 } else if (der2 >= d && der2 < d + sNv) { 328 H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sD[(spdim * p + si) * sNv + der2 - d]; 329 } else { 330 H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sB[spdim * p + si]; 331 } 332 } 333 } 334 } 335 } 336 } 337 } 338 } 339 d += sNv; 340 step *= spdim; 341 } 342 if (B && Nc > 1) { 343 /* Make direct sum basis for multicomponent space */ 344 for (p = 0; p < npoints; ++p) { 345 for (i = 0; i < pdim; ++i) { 346 for (c = 1; c < Nc; ++c) { 347 B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc]; 348 } 349 } 350 } 351 } 352 if (D && Nc > 1) { 353 /* Make direct sum basis for multicomponent space */ 354 for (p = 0; p < npoints; ++p) { 355 for (i = 0; i < pdim; ++i) { 356 for (c = 1; c < Nc; ++c) { 357 for (d = 0; d < Nv; ++d) { 358 D[((p*pdim*Nc + i*Nc + c)*Nc + c)*Nv + d] = D[(p*pdim + i)*Nc*Nc*Nv + d]; 359 } 360 } 361 } 362 } 363 } 364 if (H && Nc > 1) { 365 /* Make direct sum basis for multicomponent space */ 366 for (p = 0; p < npoints; ++p) { 367 for (i = 0; i < pdim; ++i) { 368 for (c = 1; c < Nc; ++c) { 369 for (d = 0; d < Nv; ++d) { 370 for (e = 0; e < Nv; ++e) { 371 H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*Nv + d)*Nv + e] = H[((p*pdim + i)*Nc*Nc*Nv + d)*Nv + e]; 372 } 373 } 374 } 375 } 376 } 377 } 378 if (H) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nv*Nv, MPIU_REAL, &sH);CHKERRQ(ierr);} 379 if (D || H) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nv, MPIU_REAL, &sD);CHKERRQ(ierr);} 380 if (B || D || H) {ierr = DMRestoreWorkArray(dm, npoints*pdim, MPIU_REAL, &sB);CHKERRQ(ierr);} 381 ierr = DMRestoreWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);CHKERRQ(ierr); 382 PetscFunctionReturn(0); 383 } 384 385 PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numSpaces) 386 { 387 PetscErrorCode ierr; 388 389 PetscFunctionBegin; 390 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 391 ierr = PetscTryMethod(sp,"PetscSpaceTensorSetNumSubspaces_C",(PetscSpace,PetscInt),(sp,numSpaces));CHKERRQ(ierr); 392 PetscFunctionReturn(0); 393 } 394 395 PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numSpaces) 396 { 397 PetscErrorCode ierr; 398 399 PetscFunctionBegin; 400 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 401 PetscValidIntPointer(numSpaces, 2); 402 ierr = PetscTryMethod(sp,"PetscSpaceTensorGetNumSubspaces_C",(PetscSpace,PetscInt*),(sp,numSpaces));CHKERRQ(ierr); 403 PetscFunctionReturn(0); 404 } 405 406 PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp) 407 { 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 412 if (subsp) PetscValidHeaderSpecific(subsp, PETSCSPACE_CLASSID, 3); 413 ierr = PetscTryMethod(sp,"PetscSpaceTensorSetSubspace_C",(PetscSpace,PetscInt,PetscSpace),(sp,s,subsp));CHKERRQ(ierr); 414 PetscFunctionReturn(0); 415 } 416 417 PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp) 418 { 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 423 PetscValidPointer(subsp, 3); 424 ierr = PetscTryMethod(sp,"PetscSpaceTensorGetSubspace_C",(PetscSpace,PetscInt,PetscSpace*),(sp,s,subsp));CHKERRQ(ierr); 425 PetscFunctionReturn(0); 426 } 427 428 static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numSpaces) 429 { 430 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data; 431 PetscInt Ns; 432 PetscErrorCode ierr; 433 434 PetscFunctionBegin; 435 if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change number of subspaces after setup called\n"); 436 Ns = tens->numSpaces; 437 if (numSpaces == Ns) PetscFunctionReturn(0); 438 if (Ns >= 0) { 439 PetscInt s; 440 441 for (s = 0; s < Ns; s++) {ierr = PetscSpaceDestroy(&tens->spaces[s]);CHKERRQ(ierr);} 442 ierr = PetscFree(tens->spaces);CHKERRQ(ierr); 443 } 444 Ns = tens->numSpaces = numSpaces; 445 ierr = PetscCalloc1(Ns, &tens->spaces);CHKERRQ(ierr); 446 PetscFunctionReturn(0); 447 } 448 449 static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numSpaces) 450 { 451 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data; 452 453 PetscFunctionBegin; 454 *numSpaces = tens->numSpaces; 455 PetscFunctionReturn(0); 456 } 457 458 static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace) 459 { 460 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data; 461 PetscInt Ns; 462 PetscErrorCode ierr; 463 464 PetscFunctionBegin; 465 if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change subspace after setup called\n"); 466 Ns = tens->numSpaces; 467 if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n"); 468 if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace); 469 ierr = PetscObjectReference((PetscObject)subspace);CHKERRQ(ierr); 470 ierr = PetscSpaceDestroy(&tens->spaces[s]);CHKERRQ(ierr); 471 tens->spaces[s] = subspace; 472 PetscFunctionReturn(0); 473 } 474 475 static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace) 476 { 477 PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data; 478 PetscInt Ns; 479 480 PetscFunctionBegin; 481 Ns = tens->numSpaces; 482 if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n"); 483 if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace); 484 *subspace = tens->spaces[s]; 485 PetscFunctionReturn(0); 486 } 487 488 static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp) 489 { 490 PetscErrorCode ierr; 491 492 PetscFunctionBegin; 493 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor; 494 sp->ops->setup = PetscSpaceSetUp_Tensor; 495 sp->ops->view = PetscSpaceView_Tensor; 496 sp->ops->destroy = PetscSpaceDestroy_Tensor; 497 sp->ops->getdimension = PetscSpaceGetDimension_Tensor; 498 sp->ops->evaluate = PetscSpaceEvaluate_Tensor; 499 sp->ops->getheightsubspace = NULL; 500 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor);CHKERRQ(ierr); 501 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor);CHKERRQ(ierr); 502 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor);CHKERRQ(ierr); 503 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor);CHKERRQ(ierr); 504 PetscFunctionReturn(0); 505 } 506 507 /*MC 508 PETSCSPACETENSOR = "tensor" - A PetscSpace object that encapsulates a tensor product space. 509 Subspaces are scalar spaces (num of componenents = 1), vector components 510 are assumed to be identical. 511 512 Level: intermediate 513 514 .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType() 515 M*/ 516 517 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp) 518 { 519 PetscSpace_Tensor *tens; 520 PetscErrorCode ierr; 521 522 PetscFunctionBegin; 523 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 524 ierr = PetscNewLog(sp,&tens);CHKERRQ(ierr); 525 sp->data = tens; 526 527 tens->numSpaces = PETSC_DEFAULT; 528 tens->dim = PETSC_DEFAULT; 529 530 ierr = PetscSpaceInitialize_Tensor(sp);CHKERRQ(ierr); 531 PetscFunctionReturn(0); 532 } 533 534