1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 typedef struct { 4 PetscDualSpace dualSubspace; 5 PetscSpace origSpace; 6 PetscReal *x; 7 PetscReal *x_alloc; 8 PetscReal *Jx; 9 PetscReal *Jx_alloc; 10 PetscReal *u; 11 PetscReal *u_alloc; 12 PetscReal *Ju; 13 PetscReal *Ju_alloc; 14 PetscReal *Q; 15 PetscInt Nb; 16 } PetscSpace_Subspace; 17 18 static PetscErrorCode PetscSpaceDestroy_Subspace(PetscSpace sp) 19 { 20 PetscSpace_Subspace *subsp; 21 22 PetscFunctionBegin; 23 subsp = (PetscSpace_Subspace *)sp->data; 24 subsp->x = NULL; 25 PetscCall(PetscFree(subsp->x_alloc)); 26 subsp->Jx = NULL; 27 PetscCall(PetscFree(subsp->Jx_alloc)); 28 subsp->u = NULL; 29 PetscCall(PetscFree(subsp->u_alloc)); 30 subsp->Ju = NULL; 31 PetscCall(PetscFree(subsp->Ju_alloc)); 32 PetscCall(PetscFree(subsp->Q)); 33 PetscCall(PetscSpaceDestroy(&subsp->origSpace)); 34 PetscCall(PetscDualSpaceDestroy(&subsp->dualSubspace)); 35 PetscCall(PetscFree(subsp)); 36 sp->data = NULL; 37 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", NULL)); 38 PetscFunctionReturn(PETSC_SUCCESS); 39 } 40 41 static PetscErrorCode PetscSpaceView_Subspace(PetscSpace sp, PetscViewer viewer) 42 { 43 PetscBool iascii; 44 PetscSpace_Subspace *subsp; 45 46 PetscFunctionBegin; 47 subsp = (PetscSpace_Subspace *)sp->data; 48 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 49 if (iascii) { 50 PetscInt origDim, subDim, origNc, subNc, o, s; 51 52 PetscCall(PetscSpaceGetNumVariables(subsp->origSpace, &origDim)); 53 PetscCall(PetscSpaceGetNumComponents(subsp->origSpace, &origNc)); 54 PetscCall(PetscSpaceGetNumVariables(sp, &subDim)); 55 PetscCall(PetscSpaceGetNumComponents(sp, &subNc)); 56 if (subsp->x) { 57 PetscCall(PetscViewerASCIIPrintf(viewer, "Subspace-to-space domain shift:\n\n")); 58 for (o = 0; o < origDim; o++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)subsp->x[o])); 59 } 60 if (subsp->Jx) { 61 PetscCall(PetscViewerASCIIPrintf(viewer, "Subspace-to-space domain transform:\n\n")); 62 for (o = 0; o < origDim; o++) { 63 PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Jx[o * subDim + 0])); 64 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 65 for (s = 1; s < subDim; s++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Jx[o * subDim + s])); 66 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 67 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 68 } 69 } 70 if (subsp->u) { 71 PetscCall(PetscViewerASCIIPrintf(viewer, "Space-to-subspace range shift:\n\n")); 72 for (o = 0; o < origNc; o++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)subsp->u[o])); 73 } 74 if (subsp->Ju) { 75 PetscCall(PetscViewerASCIIPrintf(viewer, "Space-to-subsace domain transform:\n")); 76 for (o = 0; o < origNc; o++) { 77 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 78 for (s = 0; s < subNc; s++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Ju[o * subNc + s])); 79 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 80 } 81 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 82 } 83 PetscCall(PetscViewerASCIIPrintf(viewer, "Original space:\n")); 84 } 85 PetscCall(PetscViewerASCIIPushTab(viewer)); 86 PetscCall(PetscSpaceView(subsp->origSpace, viewer)); 87 PetscCall(PetscViewerASCIIPopTab(viewer)); 88 PetscFunctionReturn(PETSC_SUCCESS); 89 } 90 91 static PetscErrorCode PetscSpaceEvaluate_Subspace(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 92 { 93 PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data; 94 PetscSpace origsp; 95 PetscInt origDim, subDim, origNc, subNc, subNb, origNb, i, j, k, l, m, n, o; 96 PetscReal *inpoints, *inB = NULL, *inD = NULL, *inH = NULL; 97 98 PetscFunctionBegin; 99 origsp = subsp->origSpace; 100 PetscCall(PetscSpaceGetNumVariables(sp, &subDim)); 101 PetscCall(PetscSpaceGetNumVariables(origsp, &origDim)); 102 PetscCall(PetscSpaceGetNumComponents(sp, &subNc)); 103 PetscCall(PetscSpaceGetNumComponents(origsp, &origNc)); 104 PetscCall(PetscSpaceGetDimension(sp, &subNb)); 105 PetscCall(PetscSpaceGetDimension(origsp, &origNb)); 106 PetscCall(DMGetWorkArray(sp->dm, npoints * origDim, MPIU_REAL, &inpoints)); 107 for (i = 0; i < npoints; i++) { 108 if (subsp->x) { 109 for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = subsp->x[j]; 110 } else { 111 for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = 0.0; 112 } 113 if (subsp->Jx) { 114 for (j = 0; j < origDim; j++) { 115 for (k = 0; k < subDim; k++) inpoints[i * origDim + j] += subsp->Jx[j * subDim + k] * points[i * subDim + k]; 116 } 117 } else { 118 for (j = 0; j < PetscMin(subDim, origDim); j++) inpoints[i * origDim + j] += points[i * subDim + j]; 119 } 120 } 121 if (B) PetscCall(DMGetWorkArray(sp->dm, npoints * origNb * origNc, MPIU_REAL, &inB)); 122 if (D) PetscCall(DMGetWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inD)); 123 if (H) PetscCall(DMGetWorkArray(sp->dm, npoints * origNb * origNc * origDim * origDim, MPIU_REAL, &inH)); 124 PetscCall(PetscSpaceEvaluate(origsp, npoints, inpoints, inB, inD, inH)); 125 if (H) { 126 PetscReal *phi, *psi; 127 128 PetscCall(DMGetWorkArray(sp->dm, origNc * origDim * origDim, MPIU_REAL, &phi)); 129 PetscCall(DMGetWorkArray(sp->dm, origNc * subDim * subDim, MPIU_REAL, &psi)); 130 for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0; 131 for (i = 0; i < subNb; i++) { 132 const PetscReal *subq = &subsp->Q[i * origNb]; 133 134 for (j = 0; j < npoints; j++) { 135 for (k = 0; k < origNc * origDim; k++) phi[k] = 0.; 136 for (k = 0; k < origNc * subDim; k++) psi[k] = 0.; 137 for (k = 0; k < origNb; k++) { 138 for (l = 0; l < origNc * origDim * origDim; l++) phi[l] += inH[(j * origNb + k) * origNc * origDim * origDim + l] * subq[k]; 139 } 140 if (subsp->Jx) { 141 for (k = 0; k < subNc; k++) { 142 for (l = 0; l < subDim; l++) { 143 for (m = 0; m < origDim; m++) { 144 for (n = 0; n < subDim; n++) { 145 for (o = 0; o < origDim; o++) psi[(k * subDim + l) * subDim + n] += subsp->Jx[m * subDim + l] * subsp->Jx[o * subDim + n] * phi[(k * origDim + m) * origDim + o]; 146 } 147 } 148 } 149 } 150 } else { 151 for (k = 0; k < subNc; k++) { 152 for (l = 0; l < PetscMin(subDim, origDim); l++) { 153 for (m = 0; m < PetscMin(subDim, origDim); m++) psi[(k * subDim + l) * subDim + m] += phi[(k * origDim + l) * origDim + m]; 154 } 155 } 156 } 157 if (subsp->Ju) { 158 for (k = 0; k < subNc; k++) { 159 for (l = 0; l < origNc; l++) { 160 for (m = 0; m < subDim * subDim; m++) H[((j * subNb + i) * subNc + k) * subDim * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim * subDim + m]; 161 } 162 } 163 } else { 164 for (k = 0; k < PetscMin(subNc, origNc); k++) { 165 for (l = 0; l < subDim * subDim; l++) H[((j * subNb + i) * subNc + k) * subDim * subDim + l] += psi[k * subDim * subDim + l]; 166 } 167 } 168 } 169 } 170 PetscCall(DMRestoreWorkArray(sp->dm, subNc * origDim, MPIU_REAL, &psi)); 171 PetscCall(DMRestoreWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi)); 172 PetscCall(DMRestoreWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inH)); 173 } 174 if (D) { 175 PetscReal *phi, *psi; 176 177 PetscCall(DMGetWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi)); 178 PetscCall(DMGetWorkArray(sp->dm, origNc * subDim, MPIU_REAL, &psi)); 179 for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0; 180 for (i = 0; i < subNb; i++) { 181 const PetscReal *subq = &subsp->Q[i * origNb]; 182 183 for (j = 0; j < npoints; j++) { 184 for (k = 0; k < origNc * origDim; k++) phi[k] = 0.; 185 for (k = 0; k < origNc * subDim; k++) psi[k] = 0.; 186 for (k = 0; k < origNb; k++) { 187 for (l = 0; l < origNc * origDim; l++) phi[l] += inD[(j * origNb + k) * origNc * origDim + l] * subq[k]; 188 } 189 if (subsp->Jx) { 190 for (k = 0; k < subNc; k++) { 191 for (l = 0; l < subDim; l++) { 192 for (m = 0; m < origDim; m++) psi[k * subDim + l] += subsp->Jx[m * subDim + l] * phi[k * origDim + m]; 193 } 194 } 195 } else { 196 for (k = 0; k < subNc; k++) { 197 for (l = 0; l < PetscMin(subDim, origDim); l++) psi[k * subDim + l] += phi[k * origDim + l]; 198 } 199 } 200 if (subsp->Ju) { 201 for (k = 0; k < subNc; k++) { 202 for (l = 0; l < origNc; l++) { 203 for (m = 0; m < subDim; m++) D[((j * subNb + i) * subNc + k) * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim + m]; 204 } 205 } 206 } else { 207 for (k = 0; k < PetscMin(subNc, origNc); k++) { 208 for (l = 0; l < subDim; l++) D[((j * subNb + i) * subNc + k) * subDim + l] += psi[k * subDim + l]; 209 } 210 } 211 } 212 } 213 PetscCall(DMRestoreWorkArray(sp->dm, subNc * origDim, MPIU_REAL, &psi)); 214 PetscCall(DMRestoreWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi)); 215 PetscCall(DMRestoreWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inD)); 216 } 217 if (B) { 218 PetscReal *phi; 219 220 PetscCall(DMGetWorkArray(sp->dm, origNc, MPIU_REAL, &phi)); 221 if (subsp->u) { 222 for (i = 0; i < npoints * subNb; i++) { 223 for (j = 0; j < subNc; j++) B[i * subNc + j] = subsp->u[j]; 224 } 225 } else { 226 for (i = 0; i < npoints * subNb * subNc; i++) B[i] = 0.0; 227 } 228 for (i = 0; i < subNb; i++) { 229 const PetscReal *subq = &subsp->Q[i * origNb]; 230 231 for (j = 0; j < npoints; j++) { 232 for (k = 0; k < origNc; k++) phi[k] = 0.; 233 for (k = 0; k < origNb; k++) { 234 for (l = 0; l < origNc; l++) phi[l] += inB[(j * origNb + k) * origNc + l] * subq[k]; 235 } 236 if (subsp->Ju) { 237 for (k = 0; k < subNc; k++) { 238 for (l = 0; l < origNc; l++) B[(j * subNb + i) * subNc + k] += subsp->Ju[k * origNc + l] * phi[l]; 239 } 240 } else { 241 for (k = 0; k < PetscMin(subNc, origNc); k++) B[(j * subNb + i) * subNc + k] += phi[k]; 242 } 243 } 244 } 245 PetscCall(DMRestoreWorkArray(sp->dm, origNc, MPIU_REAL, &phi)); 246 PetscCall(DMRestoreWorkArray(sp->dm, npoints * origNb * origNc, MPIU_REAL, &inB)); 247 } 248 PetscCall(DMRestoreWorkArray(sp->dm, npoints * origDim, MPIU_REAL, &inpoints)); 249 PetscFunctionReturn(PETSC_SUCCESS); 250 } 251 252 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp) 253 { 254 PetscSpace_Subspace *subsp; 255 256 PetscFunctionBegin; 257 PetscCall(PetscNew(&subsp)); 258 sp->data = (void *)subsp; 259 PetscFunctionReturn(PETSC_SUCCESS); 260 } 261 262 static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim) 263 { 264 PetscSpace_Subspace *subsp; 265 266 PetscFunctionBegin; 267 subsp = (PetscSpace_Subspace *)sp->data; 268 *dim = subsp->Nb; 269 PetscFunctionReturn(PETSC_SUCCESS); 270 } 271 272 static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp) 273 { 274 const PetscReal *x; 275 const PetscReal *Jx; 276 const PetscReal *u; 277 const PetscReal *Ju; 278 PetscDualSpace dualSubspace; 279 PetscSpace origSpace; 280 PetscInt origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset; 281 PetscReal *allPoints, *allWeights, *B, *V; 282 DM dm; 283 PetscSpace_Subspace *subsp; 284 285 PetscFunctionBegin; 286 subsp = (PetscSpace_Subspace *)sp->data; 287 x = subsp->x; 288 Jx = subsp->Jx; 289 u = subsp->u; 290 Ju = subsp->Ju; 291 origSpace = subsp->origSpace; 292 dualSubspace = subsp->dualSubspace; 293 PetscCall(PetscSpaceGetNumComponents(origSpace, &origNc)); 294 PetscCall(PetscSpaceGetNumVariables(origSpace, &origDim)); 295 PetscCall(PetscDualSpaceGetDM(dualSubspace, &dm)); 296 PetscCall(DMGetDimension(dm, &subDim)); 297 PetscCall(PetscSpaceGetDimension(origSpace, &origNb)); 298 PetscCall(PetscDualSpaceGetDimension(dualSubspace, &subNb)); 299 PetscCall(PetscDualSpaceGetNumComponents(dualSubspace, &subNc)); 300 301 for (f = 0, numPoints = 0; f < subNb; f++) { 302 PetscQuadrature q; 303 PetscInt qNp; 304 305 PetscCall(PetscDualSpaceGetFunctional(dualSubspace, f, &q)); 306 PetscCall(PetscQuadratureGetData(q, NULL, NULL, &qNp, NULL, NULL)); 307 numPoints += qNp; 308 } 309 PetscCall(PetscMalloc1(subNb * origNb, &V)); 310 PetscCall(PetscMalloc3(numPoints * origDim, &allPoints, numPoints * origNc, &allWeights, numPoints * origNb * origNc, &B)); 311 for (f = 0, offset = 0; f < subNb; f++) { 312 PetscQuadrature q; 313 PetscInt qNp, p; 314 const PetscReal *qp; 315 const PetscReal *qw; 316 317 PetscCall(PetscDualSpaceGetFunctional(dualSubspace, f, &q)); 318 PetscCall(PetscQuadratureGetData(q, NULL, NULL, &qNp, &qp, &qw)); 319 for (p = 0; p < qNp; p++, offset++) { 320 if (x) { 321 for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i]; 322 } else { 323 for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0; 324 } 325 if (Jx) { 326 for (i = 0; i < origDim; i++) { 327 for (j = 0; j < subDim; j++) allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j]; 328 } 329 } else { 330 for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i]; 331 } 332 for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0; 333 if (Ju) { 334 for (i = 0; i < origNc; i++) { 335 for (j = 0; j < subNc; j++) allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i]; 336 } 337 } else { 338 for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i]; 339 } 340 } 341 } 342 PetscCall(PetscSpaceEvaluate(origSpace, numPoints, allPoints, B, NULL, NULL)); 343 for (f = 0, offset = 0; f < subNb; f++) { 344 PetscInt b, p, s, qNp; 345 PetscQuadrature q; 346 const PetscReal *qw; 347 348 PetscCall(PetscDualSpaceGetFunctional(dualSubspace, f, &q)); 349 PetscCall(PetscQuadratureGetData(q, NULL, NULL, &qNp, NULL, &qw)); 350 if (u) { 351 for (b = 0; b < origNb; b++) { 352 for (s = 0; s < subNc; s++) V[f * origNb + b] += qw[s] * u[s]; 353 } 354 } else { 355 for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0; 356 } 357 for (p = 0; p < qNp; p++, offset++) { 358 for (b = 0; b < origNb; b++) { 359 for (s = 0; s < origNc; s++) V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s]; 360 } 361 } 362 } 363 /* orthnormalize rows of V */ 364 for (f = 0; f < subNb; f++) { 365 PetscReal rho = 0.0, scal; 366 367 for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]); 368 369 scal = 1. / PetscSqrtReal(rho); 370 371 for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal; 372 for (j = f + 1; j < subNb; j++) { 373 for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i]; 374 for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal; 375 } 376 } 377 PetscCall(PetscFree3(allPoints, allWeights, B)); 378 subsp->Q = V; 379 PetscFunctionReturn(PETSC_SUCCESS); 380 } 381 382 static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly) 383 { 384 PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data; 385 386 PetscFunctionBegin; 387 *poly = PETSC_FALSE; 388 PetscCall(PetscSpacePolynomialGetTensor(subsp->origSpace, poly)); 389 if (*poly) { 390 if (subsp->Jx) { 391 PetscInt subDim, origDim, i, j; 392 PetscInt maxnnz; 393 394 PetscCall(PetscSpaceGetNumVariables(subsp->origSpace, &origDim)); 395 PetscCall(PetscSpaceGetNumVariables(sp, &subDim)); 396 maxnnz = 0; 397 for (i = 0; i < origDim; i++) { 398 PetscInt nnz = 0; 399 400 for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.); 401 maxnnz = PetscMax(maxnnz, nnz); 402 } 403 for (j = 0; j < subDim; j++) { 404 PetscInt nnz = 0; 405 406 for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.); 407 maxnnz = PetscMax(maxnnz, nnz); 408 } 409 if (maxnnz > 1) *poly = PETSC_FALSE; 410 } 411 } 412 PetscFunctionReturn(PETSC_SUCCESS); 413 } 414 415 static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp) 416 { 417 PetscFunctionBegin; 418 sp->ops->setup = PetscSpaceSetUp_Subspace; 419 sp->ops->view = PetscSpaceView_Subspace; 420 sp->ops->destroy = PetscSpaceDestroy_Subspace; 421 sp->ops->getdimension = PetscSpaceGetDimension_Subspace; 422 sp->ops->evaluate = PetscSpaceEvaluate_Subspace; 423 PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace)); 424 PetscFunctionReturn(PETSC_SUCCESS); 425 } 426 /*@ 427 PetscSpaceCreateSubspace - creates a subspace from a an `origSpace` and its dual `dualSubspace` 428 429 Input Parameters: 430 + origSpace - the original `PetscSpace` 431 . dualSubspace - no idea 432 . x - no idea 433 . Jx - no idea 434 . u - no idea 435 . Ju - no idea 436 - copymode - whether to copy, borrow, or own some of the input arrays I guess 437 438 Output Parameter: 439 . subspace - the subspace 440 441 Level: advanced 442 443 .seealso: `PetscSpace`, `PetscDualSpace`, `PetscCopyMode`, `PetscSpaceType` 444 @*/ 445 PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace) 446 { 447 PetscSpace_Subspace *subsp; 448 PetscInt origDim, subDim, origNc, subNc, subNb; 449 PetscInt order; 450 DM dm; 451 452 PetscFunctionBegin; 453 PetscValidHeaderSpecific(origSpace, PETSCSPACE_CLASSID, 1); 454 PetscValidHeaderSpecific(dualSubspace, PETSCDUALSPACE_CLASSID, 2); 455 if (x) PetscValidRealPointer(x, 3); 456 if (Jx) PetscValidRealPointer(Jx, 4); 457 if (u) PetscValidRealPointer(u, 5); 458 if (Ju) PetscValidRealPointer(Ju, 6); 459 PetscValidPointer(subspace, 8); 460 PetscCall(PetscSpaceGetNumComponents(origSpace, &origNc)); 461 PetscCall(PetscSpaceGetNumVariables(origSpace, &origDim)); 462 PetscCall(PetscDualSpaceGetDM(dualSubspace, &dm)); 463 PetscCall(DMGetDimension(dm, &subDim)); 464 PetscCall(PetscDualSpaceGetDimension(dualSubspace, &subNb)); 465 PetscCall(PetscDualSpaceGetNumComponents(dualSubspace, &subNc)); 466 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace), subspace)); 467 PetscCall(PetscSpaceSetType(*subspace, PETSCSPACESUBSPACE)); 468 PetscCall(PetscSpaceSetNumVariables(*subspace, subDim)); 469 PetscCall(PetscSpaceSetNumComponents(*subspace, subNc)); 470 PetscCall(PetscSpaceGetDegree(origSpace, &order, NULL)); 471 PetscCall(PetscSpaceSetDegree(*subspace, order, PETSC_DETERMINE)); 472 subsp = (PetscSpace_Subspace *)(*subspace)->data; 473 subsp->Nb = subNb; 474 switch (copymode) { 475 case PETSC_OWN_POINTER: 476 if (x) subsp->x_alloc = x; 477 if (Jx) subsp->Jx_alloc = Jx; 478 if (u) subsp->u_alloc = u; 479 if (Ju) subsp->Ju_alloc = Ju; 480 /* fall through */ 481 case PETSC_USE_POINTER: 482 if (x) subsp->x = x; 483 if (Jx) subsp->Jx = Jx; 484 if (u) subsp->u = u; 485 if (Ju) subsp->Ju = Ju; 486 break; 487 case PETSC_COPY_VALUES: 488 if (x) { 489 PetscCall(PetscMalloc1(origDim, &subsp->x_alloc)); 490 PetscCall(PetscArraycpy(subsp->x_alloc, x, origDim)); 491 subsp->x = subsp->x_alloc; 492 } 493 if (Jx) { 494 PetscCall(PetscMalloc1(origDim * subDim, &subsp->Jx_alloc)); 495 PetscCall(PetscArraycpy(subsp->Jx_alloc, Jx, origDim * subDim)); 496 subsp->Jx = subsp->Jx_alloc; 497 } 498 if (u) { 499 PetscCall(PetscMalloc1(subNc, &subsp->u_alloc)); 500 PetscCall(PetscArraycpy(subsp->u_alloc, u, subNc)); 501 subsp->u = subsp->u_alloc; 502 } 503 if (Ju) { 504 PetscCall(PetscMalloc1(origNc * subNc, &subsp->Ju_alloc)); 505 PetscCall(PetscArraycpy(subsp->Ju_alloc, Ju, origNc * subNc)); 506 subsp->Ju = subsp->Ju_alloc; 507 } 508 break; 509 default: 510 SETERRQ(PetscObjectComm((PetscObject)origSpace), PETSC_ERR_ARG_OUTOFRANGE, "Unknown copy mode"); 511 } 512 PetscCall(PetscObjectReference((PetscObject)origSpace)); 513 subsp->origSpace = origSpace; 514 PetscCall(PetscObjectReference((PetscObject)dualSubspace)); 515 subsp->dualSubspace = dualSubspace; 516 PetscCall(PetscSpaceInitialize_Subspace(*subspace)); 517 PetscFunctionReturn(PETSC_SUCCESS); 518 } 519