1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 static PetscErrorCode PetscSpaceSetFromOptions_Ptrimmed(PetscOptionItems *PetscOptionsObject,PetscSpace sp) 4 { 5 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 6 PetscErrorCode ierr; 7 8 PetscFunctionBegin; 9 ierr = PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");CHKERRQ(ierr); 10 ierr = PetscOptionsInt("-petscspace_ptrimmed_form_degree", "form degree of trimmed space", "PetscSpacePTrimmedSetFormDegree", pt->formDegree, &(pt->formDegree), NULL);CHKERRQ(ierr); 11 ierr = PetscOptionsTail();CHKERRQ(ierr); 12 PetscFunctionReturn(0); 13 } 14 15 static PetscErrorCode PetscSpacePTrimmedView_Ascii(PetscSpace sp, PetscViewer v) 16 { 17 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 18 PetscInt f; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 f = pt->formDegree; 23 ierr = PetscViewerASCIIPrintf(v, "Trimmed polynomials %D%s-forms of degree %D\n", PetscAbsInt(f), f < 0 ? "*" : "", sp->maxDegree);CHKERRQ(ierr); 24 PetscFunctionReturn(0); 25 } 26 27 static PetscErrorCode PetscSpaceView_Ptrimmed(PetscSpace sp, PetscViewer viewer) 28 { 29 PetscBool iascii; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 34 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 35 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 36 if (iascii) {ierr = PetscSpacePTrimmedView_Ascii(sp, viewer);CHKERRQ(ierr);} 37 PetscFunctionReturn(0); 38 } 39 40 static PetscErrorCode PetscSpaceDestroy_Ptrimmed(PetscSpace sp) 41 { 42 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 43 PetscErrorCode ierr; 44 45 PetscFunctionBegin; 46 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePTrimmedGetFormDegree_C", NULL);CHKERRQ(ierr); 47 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePTrimmedSetFormDegree_C", NULL);CHKERRQ(ierr); 48 if (pt->subspaces) { 49 PetscInt d; 50 51 for (d = 0; d < sp->Nv; ++d) { 52 ierr = PetscSpaceDestroy(&pt->subspaces[d]);CHKERRQ(ierr); 53 } 54 } 55 ierr = PetscFree(pt->subspaces);CHKERRQ(ierr); 56 ierr = PetscFree(pt);CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode PetscSpaceSetUp_Ptrimmed(PetscSpace sp) 61 { 62 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 63 PetscInt Nf; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 if (pt->setupCalled) PetscFunctionReturn(0); 68 if (pt->formDegree < -sp->Nv || pt->formDegree > sp->Nv) SETERRQ3(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Form degree %D not in valid range [%D,%D]", pt->formDegree, sp->Nv, sp->Nv); 69 ierr = PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf);CHKERRQ(ierr); 70 if (sp->Nc == PETSC_DETERMINE) { 71 sp->Nc = Nf; 72 } 73 if (sp->Nc % Nf) SETERRQ2(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "Number of components %D is not a multiple of form dimension %D\n", sp->Nc, Nf); 74 if (sp->Nc != Nf) { 75 PetscSpace subsp; 76 PetscInt nCopies = sp->Nc / Nf; 77 PetscInt Nv, deg, maxDeg; 78 PetscInt formDegree = pt->formDegree; 79 const char *prefix; 80 const char *name; 81 char subname[PETSC_MAX_PATH_LEN]; 82 83 ierr = PetscSpaceSetType(sp, PETSCSPACESUM);CHKERRQ(ierr); 84 ierr = PetscSpaceSumSetConcatenate(sp, PETSC_TRUE);CHKERRQ(ierr); 85 ierr = PetscSpaceSumSetNumSubspaces(sp, nCopies);CHKERRQ(ierr); 86 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);CHKERRQ(ierr); 87 ierr = PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);CHKERRQ(ierr); 88 ierr = PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix);CHKERRQ(ierr); 89 ierr = PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_");CHKERRQ(ierr); 90 if (((PetscObject)sp)->name) { 91 ierr = PetscObjectGetName((PetscObject)sp, &name);CHKERRQ(ierr); 92 ierr = PetscSNPrintf(subname, PETSC_MAX_PATH_LEN-1, "%s sum component", name);CHKERRQ(ierr); 93 ierr = PetscObjectSetName((PetscObject)subsp, subname);CHKERRQ(ierr); 94 } else { 95 ierr = PetscObjectSetName((PetscObject)subsp, "sum component");CHKERRQ(ierr); 96 } 97 ierr = PetscSpaceSetType(subsp, PETSCSPACEPTRIMMED);CHKERRQ(ierr); 98 ierr = PetscSpaceGetNumVariables(sp, &Nv);CHKERRQ(ierr); 99 ierr = PetscSpaceSetNumVariables(subsp, Nv);CHKERRQ(ierr); 100 ierr = PetscSpaceSetNumComponents(subsp, Nf);CHKERRQ(ierr); 101 ierr = PetscSpaceGetDegree(sp, °, &maxDeg);CHKERRQ(ierr); 102 ierr = PetscSpaceSetDegree(subsp, deg, maxDeg);CHKERRQ(ierr); 103 ierr = PetscSpacePTrimmedSetFormDegree(subsp, formDegree);CHKERRQ(ierr); 104 ierr = PetscSpaceSetUp(subsp);CHKERRQ(ierr); 105 for (PetscInt i = 0; i < nCopies; i++) { 106 ierr = PetscSpaceSumSetSubspace(sp, i, subsp);CHKERRQ(ierr); 107 } 108 ierr = PetscSpaceDestroy(&subsp);CHKERRQ(ierr); 109 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 110 PetscFunctionReturn(0); 111 } 112 if (sp->degree == PETSC_DEFAULT) { 113 sp->degree = 0; 114 } else if (sp->degree < 0) { 115 SETERRQ1(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Invalid negative degree %D\n", sp->degree); 116 } 117 if (pt->formDegree == 0 || PetscAbsInt(pt->formDegree) == sp->Nv) { 118 // Convert to regular polynomial space 119 ierr = PetscSpaceSetType(sp, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 120 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 pt->setupCalled = PETSC_TRUE; 124 PetscFunctionReturn(0); 125 } 126 127 static PetscErrorCode PetscSpaceGetDimension_Ptrimmed(PetscSpace sp, PetscInt *dim) 128 { 129 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 130 PetscInt f; 131 PetscInt Nf; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 f = pt->formDegree; 136 // For PetscSpace, degree refers to the largest complete polynomial degree contained in the space which 137 // is equal to the index of a P trimmed space only for 0-forms: otherwise, the index is degree + 1 138 ierr = PetscDTPTrimmedSize(sp->Nv, f == 0 ? sp->degree : sp->degree + 1, pt->formDegree, dim);CHKERRQ(ierr); 139 ierr = PetscDTBinomialInt(sp->Nv, PetscAbsInt(pt->formDegree), &Nf);CHKERRQ(ierr); 140 *dim *= (sp->Nc / Nf); 141 PetscFunctionReturn(0); 142 } 143 144 /* 145 p in [0, npoints), i in [0, pdim), c in [0, Nc) 146 147 B[p][i][c] = B[p][i_scalar][c][c] 148 */ 149 static PetscErrorCode PetscSpaceEvaluate_Ptrimmed(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 150 { 151 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 152 DM dm = sp->dm; 153 PetscInt jet, degree, Nf, Ncopies, Njet; 154 PetscInt Nc = sp->Nc; 155 PetscInt f; 156 PetscInt dim = sp->Nv; 157 PetscReal *eval; 158 PetscInt Nb; 159 PetscErrorCode ierr; 160 161 PetscFunctionBegin; 162 if (!pt->setupCalled) { 163 ierr = PetscSpaceSetUp(sp);CHKERRQ(ierr); 164 ierr = PetscSpaceEvaluate(sp, npoints, points, B, D, H);CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 if (H) { 168 jet = 2; 169 } else if (D) { 170 jet = 1; 171 } else { 172 jet = 0; 173 } 174 f = pt->formDegree; 175 degree = f == 0 ? sp->degree : sp->degree + 1; 176 ierr = PetscDTBinomialInt(dim, PetscAbsInt(f), &Nf);CHKERRQ(ierr); 177 Ncopies = Nc / Nf; 178 if (Ncopies != 1) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_PLIB, "Multicopy spaces should have been converted to PETSCSPACESUM"); 179 ierr = PetscDTBinomialInt(dim + jet, dim, &Njet);CHKERRQ(ierr); 180 ierr = PetscDTPTrimmedSize(dim, degree, f, &Nb);CHKERRQ(ierr); 181 ierr = DMGetWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval);CHKERRQ(ierr); 182 ierr = PetscDTPTrimmedEvalJet(dim, npoints, points, degree, f, jet, eval);CHKERRQ(ierr); 183 if (B) { 184 PetscInt p_strl = Nf*Nb; 185 PetscInt b_strl = Nf; 186 PetscInt v_strl = 1; 187 188 PetscInt b_strr = Nf*Njet*npoints; 189 PetscInt v_strr = Njet*npoints; 190 PetscInt p_strr = 1; 191 192 for (PetscInt v = 0; v < Nf; v++) { 193 for (PetscInt b = 0; b < Nb; b++) { 194 for (PetscInt p = 0; p < npoints; p++) { 195 B[p*p_strl + b*b_strl + v*v_strl] = eval[b*b_strr + v*v_strr + p*p_strr]; 196 } 197 } 198 } 199 } 200 if (D) { 201 PetscInt p_strl = dim*Nb*Ncopies*Nc; 202 PetscInt b_strl = dim*Nb*Nc; 203 PetscInt c_strl = dim*(Nc + Nf); 204 PetscInt v_strl = dim; 205 PetscInt d_strl = 1; 206 207 PetscInt b_strr = Nf*Njet*npoints; 208 PetscInt v_strr = Njet*npoints; 209 PetscInt d_strr = npoints; 210 PetscInt p_strr = 1; 211 212 for (PetscInt c = 0; c < Ncopies; c++) { 213 for (PetscInt v = 0; v < Nf; v++) { 214 for (PetscInt d = 0; d < dim; d++) { 215 for (PetscInt b = 0; b < Nb; b++) { 216 for (PetscInt p = 0; p < npoints; p++) { 217 D[p*p_strl + b*b_strl + c*c_strl + v*v_strl + d*d_strl] = eval[b*b_strr + v*v_strr + (1+d)*d_strr + p*p_strr]; 218 } 219 } 220 } 221 } 222 } 223 } 224 if (H) { 225 PetscInt p_strl = dim*dim*Nb*Ncopies*Nc; 226 PetscInt b_strl = dim*dim*Nb*Nc; 227 PetscInt c_strl = dim*dim*(Nc + Nf); 228 PetscInt v_strl = dim*dim*1; 229 PetscInt d1_strl = dim; 230 PetscInt d2_strl = 1; 231 232 PetscInt b_strr = Nf*Njet*npoints; 233 PetscInt v_strr = Njet*npoints; 234 PetscInt j_strr = npoints; 235 PetscInt p_strr = 1; 236 237 PetscInt *derivs; 238 ierr = PetscCalloc1(dim, &derivs);CHKERRQ(ierr); 239 for (PetscInt d1 = 0; d1 < dim; d1++) { 240 for (PetscInt d2 = 0; d2 < dim; d2++) { 241 PetscInt j; 242 derivs[d1]++; 243 derivs[d2]++; 244 ierr = PetscDTGradedOrderToIndex(dim, derivs, &j);CHKERRQ(ierr); 245 derivs[d1]--; 246 derivs[d2]--; 247 for (PetscInt c = 0; c < Ncopies; c++) { 248 for (PetscInt v = 0; v < Nf; v++) { 249 for (PetscInt b = 0; b < Nb; b++) { 250 for (PetscInt p = 0; p < npoints; p++) { 251 H[p*p_strl + b*b_strl + c*c_strl + v*v_strl + d1*d1_strl + d2*d2_strl] = eval[b*b_strr + v*v_strr + j*j_strr + p*p_strr]; 252 } 253 } 254 } 255 } 256 } 257 } 258 ierr = PetscFree(derivs);CHKERRQ(ierr); 259 } 260 ierr = DMRestoreWorkArray(dm, Nb * Nf * Njet * npoints, MPIU_REAL, &eval);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 /*@ 265 PetscSpacePTrimmedSetFormDegree - Set the form degree of the trimmed polynomials. 266 267 Input Parameters: 268 + sp - the function space object 269 - formDegree - the form degree 270 271 Options Database: 272 . -petscspace_ptrimmed_form_degree <int> - The trimmed polynomial form degree 273 274 Level: intermediate 275 276 .seealso: PetscDTAltV, PetscDTPTrimmedEvalJet(), PetscSpacePTrimmedGetFormDegree() 277 @*/ 278 PetscErrorCode PetscSpacePTrimmedSetFormDegree(PetscSpace sp, PetscInt formDegree) 279 { 280 PetscErrorCode ierr; 281 282 PetscFunctionBegin; 283 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 284 ierr = PetscTryMethod(sp,"PetscSpacePTrimmedSetFormDegree_C",(PetscSpace,PetscInt),(sp,formDegree));CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 /*@ 289 PetscSpacePTrimmedGetFormDegree - Get the form degree of the trimmed polynomials. 290 291 Input Parameters: 292 . sp - the function space object 293 294 Output Parameters: 295 . formDegee - the form degree 296 297 Level: intermediate 298 299 .seealso: PetscDTAltV, PetscDTPTrimmedEvalJet(), PetscSpacePTrimmedSetFormDegree() 300 @*/ 301 PetscErrorCode PetscSpacePTrimmedGetFormDegree(PetscSpace sp, PetscInt *formDegree) 302 { 303 PetscErrorCode ierr; 304 305 PetscFunctionBegin; 306 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 307 PetscValidPointer(formDegree, 2); 308 ierr = PetscTryMethod(sp,"PetscSpacePTrimmedGetFormDegree_C",(PetscSpace,PetscInt*),(sp,formDegree));CHKERRQ(ierr); 309 PetscFunctionReturn(0); 310 } 311 312 static PetscErrorCode PetscSpacePTrimmedSetFormDegree_Ptrimmed(PetscSpace sp, PetscInt formDegree) 313 { 314 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 315 316 PetscFunctionBegin; 317 pt->formDegree = formDegree; 318 PetscFunctionReturn(0); 319 } 320 321 static PetscErrorCode PetscSpacePTrimmedGetFormDegree_Ptrimmed(PetscSpace sp, PetscInt *formDegree) 322 { 323 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 324 325 PetscFunctionBegin; 326 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 327 PetscValidPointer(formDegree, 2); 328 *formDegree = pt->formDegree; 329 PetscFunctionReturn(0); 330 } 331 332 static PetscErrorCode PetscSpaceGetHeightSubspace_Ptrimmed(PetscSpace sp, PetscInt height, PetscSpace *subsp) 333 { 334 PetscSpace_Ptrimmed *pt = (PetscSpace_Ptrimmed *) sp->data; 335 PetscInt dim; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 ierr = PetscSpaceGetNumVariables(sp, &dim);CHKERRQ(ierr); 340 if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim); 341 if (!pt->subspaces) {ierr = PetscCalloc1(dim, &(pt->subspaces));CHKERRQ(ierr);} 342 if ((dim - height) <= PetscAbsInt(pt->formDegree)) { 343 if (!pt->subspaces[height-1]) { 344 PetscInt Nc, degree, Nf, Ncopies, Nfsub; 345 PetscSpace sub; 346 const char *name; 347 348 ierr = PetscSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 349 ierr = PetscDTBinomialInt(dim, PetscAbsInt(pt->formDegree), &Nf);CHKERRQ(ierr); 350 ierr = PetscDTBinomialInt((dim-height), PetscAbsInt(pt->formDegree), &Nfsub);CHKERRQ(ierr); 351 Ncopies = Nf / Nc; 352 ierr = PetscSpaceGetDegree(sp, °ree, NULL);CHKERRQ(ierr); 353 354 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);CHKERRQ(ierr); 355 ierr = PetscObjectGetName((PetscObject) sp, &name);CHKERRQ(ierr); 356 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 357 ierr = PetscSpaceSetType(sub, PETSCSPACEPTRIMMED);CHKERRQ(ierr); 358 ierr = PetscSpaceSetNumComponents(sub, Nfsub * Ncopies);CHKERRQ(ierr); 359 ierr = PetscSpaceSetDegree(sub, degree, PETSC_DETERMINE);CHKERRQ(ierr); 360 ierr = PetscSpaceSetNumVariables(sub, dim-height);CHKERRQ(ierr); 361 ierr = PetscSpacePTrimmedSetFormDegree(sub, pt->formDegree);CHKERRQ(ierr); 362 ierr = PetscSpaceSetUp(sub);CHKERRQ(ierr); 363 pt->subspaces[height-1] = sub; 364 } 365 *subsp = pt->subspaces[height-1]; 366 } else { 367 *subsp = NULL; 368 } 369 PetscFunctionReturn(0); 370 } 371 372 static PetscErrorCode PetscSpaceInitialize_Ptrimmed(PetscSpace sp) 373 { 374 PetscErrorCode ierr; 375 376 PetscFunctionBegin; 377 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePTrimmedGetFormDegree_C", PetscSpacePTrimmedGetFormDegree_Ptrimmed);CHKERRQ(ierr); 378 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePTrimmedSetFormDegree_C", PetscSpacePTrimmedSetFormDegree_Ptrimmed);CHKERRQ(ierr); 379 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Ptrimmed; 380 sp->ops->setup = PetscSpaceSetUp_Ptrimmed; 381 sp->ops->view = PetscSpaceView_Ptrimmed; 382 sp->ops->destroy = PetscSpaceDestroy_Ptrimmed; 383 sp->ops->getdimension = PetscSpaceGetDimension_Ptrimmed; 384 sp->ops->evaluate = PetscSpaceEvaluate_Ptrimmed; 385 sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Ptrimmed; 386 PetscFunctionReturn(0); 387 } 388 389 /*MC 390 PETSCSPACEPTRIMMED = "ptrimmed" - A PetscSpace object that encapsulates a trimmed polynomial space. 391 392 Level: intermediate 393 394 .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType(), PetscDTPTrimmedEvalJet() 395 M*/ 396 397 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Ptrimmed(PetscSpace sp) 398 { 399 PetscSpace_Ptrimmed *pt; 400 PetscErrorCode ierr; 401 402 PetscFunctionBegin; 403 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 404 ierr = PetscNewLog(sp,&pt);CHKERRQ(ierr); 405 sp->data = pt; 406 407 pt->subspaces = NULL; 408 sp->Nc = PETSC_DETERMINE; 409 410 ierr = PetscSpaceInitialize_Ptrimmed(sp);CHKERRQ(ierr); 411 PetscFunctionReturn(0); 412 } 413 414