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