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