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