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