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