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