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