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