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 PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace) 428 { 429 PetscSpace_Subspace *subsp; 430 PetscInt origDim, subDim, origNc, subNc, subNb; 431 PetscInt order; 432 DM dm; 433 434 PetscFunctionBegin; 435 PetscValidHeaderSpecific(origSpace, PETSCSPACE_CLASSID, 1); 436 PetscValidHeaderSpecific(dualSubspace, PETSCDUALSPACE_CLASSID, 2); 437 if (x) PetscValidRealPointer(x, 3); 438 if (Jx) PetscValidRealPointer(Jx, 4); 439 if (u) PetscValidRealPointer(u, 5); 440 if (Ju) PetscValidRealPointer(Ju, 6); 441 PetscValidPointer(subspace, 8); 442 PetscCall(PetscSpaceGetNumComponents(origSpace, &origNc)); 443 PetscCall(PetscSpaceGetNumVariables(origSpace, &origDim)); 444 PetscCall(PetscDualSpaceGetDM(dualSubspace, &dm)); 445 PetscCall(DMGetDimension(dm, &subDim)); 446 PetscCall(PetscDualSpaceGetDimension(dualSubspace, &subNb)); 447 PetscCall(PetscDualSpaceGetNumComponents(dualSubspace, &subNc)); 448 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace), subspace)); 449 PetscCall(PetscSpaceSetType(*subspace, PETSCSPACESUBSPACE)); 450 PetscCall(PetscSpaceSetNumVariables(*subspace, subDim)); 451 PetscCall(PetscSpaceSetNumComponents(*subspace, subNc)); 452 PetscCall(PetscSpaceGetDegree(origSpace, &order, NULL)); 453 PetscCall(PetscSpaceSetDegree(*subspace, order, PETSC_DETERMINE)); 454 subsp = (PetscSpace_Subspace *)(*subspace)->data; 455 subsp->Nb = subNb; 456 switch (copymode) { 457 case PETSC_OWN_POINTER: 458 if (x) subsp->x_alloc = x; 459 if (Jx) subsp->Jx_alloc = Jx; 460 if (u) subsp->u_alloc = u; 461 if (Ju) subsp->Ju_alloc = Ju; 462 /* fall through */ 463 case PETSC_USE_POINTER: 464 if (x) subsp->x = x; 465 if (Jx) subsp->Jx = Jx; 466 if (u) subsp->u = u; 467 if (Ju) subsp->Ju = Ju; 468 break; 469 case PETSC_COPY_VALUES: 470 if (x) { 471 PetscCall(PetscMalloc1(origDim, &subsp->x_alloc)); 472 PetscCall(PetscArraycpy(subsp->x_alloc, x, origDim)); 473 subsp->x = subsp->x_alloc; 474 } 475 if (Jx) { 476 PetscCall(PetscMalloc1(origDim * subDim, &subsp->Jx_alloc)); 477 PetscCall(PetscArraycpy(subsp->Jx_alloc, Jx, origDim * subDim)); 478 subsp->Jx = subsp->Jx_alloc; 479 } 480 if (u) { 481 PetscCall(PetscMalloc1(subNc, &subsp->u_alloc)); 482 PetscCall(PetscArraycpy(subsp->u_alloc, u, subNc)); 483 subsp->u = subsp->u_alloc; 484 } 485 if (Ju) { 486 PetscCall(PetscMalloc1(origNc * subNc, &subsp->Ju_alloc)); 487 PetscCall(PetscArraycpy(subsp->Ju_alloc, Ju, origNc * subNc)); 488 subsp->Ju = subsp->Ju_alloc; 489 } 490 break; 491 default: 492 SETERRQ(PetscObjectComm((PetscObject)origSpace), PETSC_ERR_ARG_OUTOFRANGE, "Unknown copy mode"); 493 } 494 PetscCall(PetscObjectReference((PetscObject)origSpace)); 495 subsp->origSpace = origSpace; 496 PetscCall(PetscObjectReference((PetscObject)dualSubspace)); 497 subsp->dualSubspace = dualSubspace; 498 PetscCall(PetscSpaceInitialize_Subspace(*subspace)); 499 PetscFunctionReturn(PETSC_SUCCESS); 500 } 501