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 306 PetscErrorCode ierr; 307 ierr = PetscNewLog(sp,&subsp);CHKERRQ(ierr); 308 sp->data = (void *) subsp; 309 PetscFunctionReturn(0); 310 } 311 312 static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim) 313 { 314 PetscSpace_Subspace *subsp; 315 316 PetscFunctionBegin; 317 subsp = (PetscSpace_Subspace *) sp->data; 318 *dim = subsp->Nb; 319 PetscFunctionReturn(0); 320 } 321 322 static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp) 323 { 324 const PetscReal *x; 325 const PetscReal *Jx; 326 const PetscReal *u; 327 const PetscReal *Ju; 328 PetscDualSpace dualSubspace; 329 PetscSpace origSpace; 330 PetscInt origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset; 331 PetscReal *allPoints, *allWeights, *B, *V; 332 DM dm; 333 PetscSpace_Subspace *subsp; 334 PetscErrorCode ierr; 335 336 PetscFunctionBegin; 337 subsp = (PetscSpace_Subspace *) sp->data; 338 x = subsp->x; 339 Jx = subsp->Jx; 340 u = subsp->u; 341 Ju = subsp->Ju; 342 origSpace = subsp->origSpace; 343 dualSubspace = subsp->dualSubspace; 344 ierr = PetscSpaceGetNumComponents(origSpace,&origNc);CHKERRQ(ierr); 345 ierr = PetscSpaceGetNumVariables(origSpace,&origDim);CHKERRQ(ierr); 346 ierr = PetscDualSpaceGetDM(dualSubspace,&dm);CHKERRQ(ierr); 347 ierr = DMGetDimension(dm,&subDim);CHKERRQ(ierr); 348 ierr = PetscSpaceGetDimension(origSpace,&origNb);CHKERRQ(ierr); 349 ierr = PetscDualSpaceGetDimension(dualSubspace,&subNb);CHKERRQ(ierr); 350 ierr = PetscDualSpaceGetNumComponents(dualSubspace,&subNc);CHKERRQ(ierr); 351 352 for (f = 0, numPoints = 0; f < subNb; f++) { 353 PetscQuadrature q; 354 PetscInt qNp; 355 356 ierr = PetscDualSpaceGetFunctional(dualSubspace,f,&q);CHKERRQ(ierr); 357 ierr = PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,NULL);CHKERRQ(ierr); 358 numPoints += qNp; 359 } 360 ierr = PetscMalloc1(subNb*origNb,&V);CHKERRQ(ierr); 361 ierr = PetscMalloc3(numPoints*origDim,&allPoints,numPoints*origNc,&allWeights,numPoints*origNb*origNc,&B);CHKERRQ(ierr); 362 for (f = 0, offset = 0; f < subNb; f++) { 363 PetscQuadrature q; 364 PetscInt qNp, p; 365 const PetscReal *qp; 366 const PetscReal *qw; 367 368 ierr = PetscDualSpaceGetFunctional(dualSubspace,f,&q);CHKERRQ(ierr); 369 ierr = PetscQuadratureGetData(q,NULL,NULL,&qNp,&qp,&qw);CHKERRQ(ierr); 370 for (p = 0; p < qNp; p++, offset++) { 371 if (x) { 372 for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i]; 373 } else { 374 for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0; 375 } 376 if (Jx) { 377 for (i = 0; i < origDim; i++) { 378 for (j = 0; j < subDim; j++) { 379 allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j]; 380 } 381 } 382 } else { 383 for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i]; 384 } 385 for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0; 386 if (Ju) { 387 for (i = 0; i < origNc; i++) { 388 for (j = 0; j < subNc; j++) { 389 allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i]; 390 } 391 } 392 } else { 393 for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i]; 394 } 395 } 396 } 397 ierr = PetscSpaceEvaluate(origSpace,numPoints,allPoints,B,NULL,NULL);CHKERRQ(ierr); 398 for (f = 0, offset = 0; f < subNb; f++) { 399 PetscInt b, p, s, qNp; 400 PetscQuadrature q; 401 const PetscReal *qw; 402 403 ierr = PetscDualSpaceGetFunctional(dualSubspace,f,&q);CHKERRQ(ierr); 404 ierr = PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,&qw);CHKERRQ(ierr); 405 if (u) { 406 for (b = 0; b < origNb; b++) { 407 for (s = 0; s < subNc; s++) { 408 V[f * origNb + b] += qw[s] * u[s]; 409 } 410 } 411 } else { 412 for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0; 413 } 414 for (p = 0; p < qNp; p++, offset++) { 415 for (b = 0; b < origNb; b++) { 416 for (s = 0; s < origNc; s++) { 417 V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s]; 418 } 419 } 420 } 421 } 422 /* orthnormalize rows of V */ 423 for (f = 0; f < subNb; f++) { 424 PetscReal rho = 0.0, scal; 425 426 for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]); 427 428 scal = 1. / PetscSqrtReal(rho); 429 430 for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal; 431 for (j = f + 1; j < subNb; j++) { 432 for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i]; 433 for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal; 434 } 435 } 436 ierr = PetscFree3(allPoints,allWeights,B);CHKERRQ(ierr); 437 subsp->Q = V; 438 PetscFunctionReturn(0); 439 } 440 441 static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly) 442 { 443 PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data; 444 PetscErrorCode ierr; 445 446 PetscFunctionBegin; 447 *poly = PETSC_FALSE; 448 ierr = PetscSpacePolynomialGetTensor(subsp->origSpace,poly);CHKERRQ(ierr); 449 if (*poly) { 450 if (subsp->Jx) { 451 PetscInt subDim, origDim, i, j; 452 PetscInt maxnnz; 453 454 ierr = PetscSpaceGetNumVariables(subsp->origSpace,&origDim);CHKERRQ(ierr); 455 ierr = PetscSpaceGetNumVariables(sp,&subDim);CHKERRQ(ierr); 456 maxnnz = 0; 457 for (i = 0; i < origDim; i++) { 458 PetscInt nnz = 0; 459 460 for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.); 461 maxnnz = PetscMax(maxnnz,nnz); 462 } 463 for (j = 0; j < subDim; j++) { 464 PetscInt nnz = 0; 465 466 for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.); 467 maxnnz = PetscMax(maxnnz,nnz); 468 } 469 if (maxnnz > 1) *poly = PETSC_FALSE; 470 } 471 } 472 PetscFunctionReturn(0); 473 } 474 475 static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp) 476 { 477 PetscErrorCode ierr; 478 479 PetscFunctionBegin; 480 sp->ops->setup = PetscSpaceSetUp_Subspace; 481 sp->ops->view = PetscSpaceView_Subspace; 482 sp->ops->destroy = PetscSpaceDestroy_Subspace; 483 sp->ops->getdimension = PetscSpaceGetDimension_Subspace; 484 sp->ops->evaluate = PetscSpaceEvaluate_Subspace; 485 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 489 PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace) 490 { 491 PetscSpace_Subspace *subsp; 492 PetscInt origDim, subDim, origNc, subNc, subNb; 493 PetscInt order; 494 DM dm; 495 PetscErrorCode ierr; 496 497 PetscFunctionBegin; 498 PetscValidHeaderSpecific(origSpace,PETSCSPACE_CLASSID,1); 499 PetscValidHeaderSpecific(dualSubspace,PETSCDUALSPACE_CLASSID,2); 500 if (x) PetscValidRealPointer(x,3); 501 if (Jx) PetscValidRealPointer(Jx,4); 502 if (u) PetscValidRealPointer(u,5); 503 if (Ju) PetscValidRealPointer(Ju,6); 504 PetscValidPointer(subspace,7); 505 ierr = PetscSpaceGetNumComponents(origSpace,&origNc);CHKERRQ(ierr); 506 ierr = PetscSpaceGetNumVariables(origSpace,&origDim);CHKERRQ(ierr); 507 ierr = PetscDualSpaceGetDM(dualSubspace,&dm);CHKERRQ(ierr); 508 ierr = DMGetDimension(dm,&subDim);CHKERRQ(ierr); 509 ierr = PetscDualSpaceGetDimension(dualSubspace,&subNb);CHKERRQ(ierr); 510 ierr = PetscDualSpaceGetNumComponents(dualSubspace,&subNc);CHKERRQ(ierr); 511 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace),subspace);CHKERRQ(ierr); 512 ierr = PetscSpaceSetType(*subspace,PETSCSPACESUBSPACE);CHKERRQ(ierr); 513 ierr = PetscSpaceSetNumVariables(*subspace,subDim);CHKERRQ(ierr); 514 ierr = PetscSpaceSetNumComponents(*subspace,subNc);CHKERRQ(ierr); 515 ierr = PetscSpaceGetDegree(origSpace,&order,NULL);CHKERRQ(ierr); 516 ierr = PetscSpaceSetDegree(*subspace,order,PETSC_DETERMINE);CHKERRQ(ierr); 517 subsp = (PetscSpace_Subspace *) (*subspace)->data; 518 subsp->Nb = subNb; 519 switch (copymode) { 520 case PETSC_OWN_POINTER: 521 if (x) subsp->x_alloc = x; 522 if (Jx) subsp->Jx_alloc = Jx; 523 if (u) subsp->u_alloc = u; 524 if (Ju) subsp->Ju_alloc = Ju; 525 case PETSC_USE_POINTER: 526 if (x) subsp->x = x; 527 if (Jx) subsp->Jx = Jx; 528 if (u) subsp->u = u; 529 if (Ju) subsp->Ju = Ju; 530 break; 531 case PETSC_COPY_VALUES: 532 if (x) { 533 ierr = PetscMalloc1(origDim,&subsp->x_alloc);CHKERRQ(ierr); 534 ierr = PetscMemcpy(subsp->x_alloc,x,origDim*sizeof(*subsp->x_alloc));CHKERRQ(ierr); 535 subsp->x = subsp->x_alloc; 536 } 537 if (Jx) { 538 ierr = PetscMalloc1(origDim * subDim,&subsp->Jx_alloc);CHKERRQ(ierr); 539 ierr = PetscMemcpy(subsp->Jx_alloc,Jx,origDim * subDim*sizeof(*subsp->Jx_alloc));CHKERRQ(ierr); 540 subsp->Jx = subsp->Jx_alloc; 541 } 542 if (u) { 543 ierr = PetscMalloc1(subNc,&subsp->u_alloc);CHKERRQ(ierr); 544 ierr = PetscMemcpy(subsp->u_alloc,u,subNc*sizeof(*subsp->u_alloc));CHKERRQ(ierr); 545 subsp->u = subsp->u_alloc; 546 } 547 if (Ju) { 548 ierr = PetscMalloc1(origNc * subNc,&subsp->Ju_alloc);CHKERRQ(ierr); 549 ierr = PetscMemcpy(subsp->Ju_alloc,Ju,origNc * subNc*sizeof(*subsp->Ju_alloc));CHKERRQ(ierr); 550 subsp->Ju = subsp->Ju_alloc; 551 } 552 break; 553 default: 554 SETERRQ(PetscObjectComm((PetscObject)origSpace),PETSC_ERR_ARG_OUTOFRANGE,"Unknown copy mode"); 555 } 556 ierr = PetscObjectReference((PetscObject)origSpace);CHKERRQ(ierr); 557 subsp->origSpace = origSpace; 558 ierr = PetscObjectReference((PetscObject)dualSubspace);CHKERRQ(ierr); 559 subsp->dualSubspace = dualSubspace; 560 ierr = PetscSpaceInitialize_Subspace(*subspace);CHKERRQ(ierr); 561 PetscFunctionReturn(0); 562 } 563 564