1 /* Discretization tools */ 2 3 #include <petscconf.h> 4 #if defined(PETSC_HAVE_MATHIMF_H) 5 #include <mathimf.h> /* this needs to be included before math.h */ 6 #endif 7 8 #include <petscdt.h> /*I "petscdt.h" I*/ /*I "petscfe.h" I*/ 9 #include <petscblaslapack.h> 10 #include <petsc-private/petscimpl.h> 11 #include <petscviewer.h> 12 #include <petscdmplex.h> 13 #include <petscdmshell.h> 14 15 #undef __FUNCT__ 16 #define __FUNCT__ "PetscQuadratureDestroy" 17 PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q) 18 { 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 ierr = PetscFree(q->quadPoints);CHKERRQ(ierr); 23 ierr = PetscFree(q->quadWeights);CHKERRQ(ierr); 24 PetscFunctionReturn(0); 25 } 26 27 #undef __FUNCT__ 28 #define __FUNCT__ "PetscDTLegendreEval" 29 /*@ 30 PetscDTLegendreEval - evaluate Legendre polynomial at points 31 32 Not Collective 33 34 Input Arguments: 35 + npoints - number of spatial points to evaluate at 36 . points - array of locations to evaluate at 37 . ndegree - number of basis degrees to evaluate 38 - degrees - sorted array of degrees to evaluate 39 40 Output Arguments: 41 + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL) 42 . D - row-oriented derivative evaluation matrix (or NULL) 43 - D2 - row-oriented second derivative evaluation matrix (or NULL) 44 45 Level: intermediate 46 47 .seealso: PetscDTGaussQuadrature() 48 @*/ 49 PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2) 50 { 51 PetscInt i,maxdegree; 52 53 PetscFunctionBegin; 54 if (!npoints || !ndegree) PetscFunctionReturn(0); 55 maxdegree = degrees[ndegree-1]; 56 for (i=0; i<npoints; i++) { 57 PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x; 58 PetscInt j,k; 59 x = points[i]; 60 pm2 = 0; 61 pm1 = 1; 62 pd2 = 0; 63 pd1 = 0; 64 pdd2 = 0; 65 pdd1 = 0; 66 k = 0; 67 if (degrees[k] == 0) { 68 if (B) B[i*ndegree+k] = pm1; 69 if (D) D[i*ndegree+k] = pd1; 70 if (D2) D2[i*ndegree+k] = pdd1; 71 k++; 72 } 73 for (j=1; j<=maxdegree; j++,k++) { 74 PetscReal p,d,dd; 75 p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j; 76 d = pd2 + (2*j-1)*pm1; 77 dd = pdd2 + (2*j-1)*pd1; 78 pm2 = pm1; 79 pm1 = p; 80 pd2 = pd1; 81 pd1 = d; 82 pdd2 = pdd1; 83 pdd1 = dd; 84 if (degrees[k] == j) { 85 if (B) B[i*ndegree+k] = p; 86 if (D) D[i*ndegree+k] = d; 87 if (D2) D2[i*ndegree+k] = dd; 88 } 89 } 90 } 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "PetscDTGaussQuadrature" 96 /*@ 97 PetscDTGaussQuadrature - create Gauss quadrature 98 99 Not Collective 100 101 Input Arguments: 102 + npoints - number of points 103 . a - left end of interval (often-1) 104 - b - right end of interval (often +1) 105 106 Output Arguments: 107 + x - quadrature points 108 - w - quadrature weights 109 110 Level: intermediate 111 112 References: 113 Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 221--230, 1969. 114 115 .seealso: PetscDTLegendreEval() 116 @*/ 117 PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w) 118 { 119 PetscErrorCode ierr; 120 PetscInt i; 121 PetscReal *work; 122 PetscScalar *Z; 123 PetscBLASInt N,LDZ,info; 124 125 PetscFunctionBegin; 126 /* Set up the Golub-Welsch system */ 127 for (i=0; i<npoints; i++) { 128 x[i] = 0; /* diagonal is 0 */ 129 if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i)); 130 } 131 ierr = PetscRealView(npoints-1,w,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 132 ierr = PetscMalloc2(npoints*npoints,PetscScalar,&Z,PetscMax(1,2*npoints-2),PetscReal,&work);CHKERRQ(ierr); 133 ierr = PetscBLASIntCast(npoints,&N);CHKERRQ(ierr); 134 LDZ = N; 135 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 136 PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info)); 137 ierr = PetscFPTrapPop();CHKERRQ(ierr); 138 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error"); 139 140 for (i=0; i<(npoints+1)/2; i++) { 141 PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */ 142 x[i] = (a+b)/2 - y*(b-a)/2; 143 x[npoints-i-1] = (a+b)/2 + y*(b-a)/2; 144 145 w[i] = w[npoints-1-i] = (b-a)*PetscSqr(0.5*PetscAbsScalar(Z[i*npoints] + Z[(npoints-i-1)*npoints])); 146 } 147 ierr = PetscFree2(Z,work);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "PetscDTFactorial_Internal" 153 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 154 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 155 PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial) 156 { 157 PetscReal f = 1.0; 158 PetscInt i; 159 160 PetscFunctionBegin; 161 for (i = 1; i < n+1; ++i) f *= i; 162 *factorial = f; 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "PetscDTComputeJacobi" 168 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x. 169 Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */ 170 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 171 { 172 PetscReal apb, pn1, pn2; 173 PetscInt k; 174 175 PetscFunctionBegin; 176 if (!n) {*P = 1.0; PetscFunctionReturn(0);} 177 if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); PetscFunctionReturn(0);} 178 apb = a + b; 179 pn2 = 1.0; 180 pn1 = 0.5 * (a - b + (apb + 2.0) * x); 181 *P = 0.0; 182 for (k = 2; k < n+1; ++k) { 183 PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0); 184 PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b); 185 PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb); 186 PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb); 187 188 a2 = a2 / a1; 189 a3 = a3 / a1; 190 a4 = a4 / a1; 191 *P = (a2 + a3 * x) * pn1 - a4 * pn2; 192 pn2 = pn1; 193 pn1 = *P; 194 } 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "PetscDTComputeJacobiDerivative" 200 /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */ 201 PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P) 202 { 203 PetscReal nP; 204 PetscErrorCode ierr; 205 206 PetscFunctionBegin; 207 if (!n) {*P = 0.0; PetscFunctionReturn(0);} 208 ierr = PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);CHKERRQ(ierr); 209 *P = 0.5 * (a + b + n + 1) * nP; 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "PetscDTMapSquareToTriangle_Internal" 215 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 216 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta) 217 { 218 PetscFunctionBegin; 219 *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0; 220 *eta = y; 221 PetscFunctionReturn(0); 222 } 223 224 #undef __FUNCT__ 225 #define __FUNCT__ "PetscDTMapCubeToTetrahedron_Internal" 226 /* Maps from [-1,1]^2 to the (-1,1) reference triangle */ 227 PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta) 228 { 229 PetscFunctionBegin; 230 *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0; 231 *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0; 232 *zeta = z; 233 PetscFunctionReturn(0); 234 } 235 236 #undef __FUNCT__ 237 #define __FUNCT__ "PetscDTGaussJacobiQuadrature1D_Internal" 238 static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w) 239 { 240 PetscInt maxIter = 100; 241 PetscReal eps = 1.0e-8; 242 PetscReal a1, a2, a3, a4, a5, a6; 243 PetscInt k; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 248 a1 = pow(2, a+b+1); 249 #if defined(PETSC_HAVE_TGAMMA) 250 a2 = tgamma(a + npoints + 1); 251 a3 = tgamma(b + npoints + 1); 252 a4 = tgamma(a + b + npoints + 1); 253 #else 254 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable."); 255 #endif 256 257 ierr = PetscDTFactorial_Internal(npoints, &a5);CHKERRQ(ierr); 258 a6 = a1 * a2 * a3 / a4 / a5; 259 /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses. 260 Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */ 261 for (k = 0; k < npoints; ++k) { 262 PetscReal r = -cos((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP; 263 PetscInt j; 264 265 if (k > 0) r = 0.5 * (r + x[k-1]); 266 for (j = 0; j < maxIter; ++j) { 267 PetscReal s = 0.0, delta, f, fp; 268 PetscInt i; 269 270 for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]); 271 ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr); 272 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);CHKERRQ(ierr); 273 delta = f / (fp - f * s); 274 r = r - delta; 275 if (fabs(delta) < eps) break; 276 } 277 x[k] = r; 278 ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);CHKERRQ(ierr); 279 w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP); 280 } 281 PetscFunctionReturn(0); 282 } 283 284 #undef __FUNCT__ 285 #define __FUNCT__ "PetscDTGaussJacobiQuadrature" 286 /*@C 287 PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex 288 289 Not Collective 290 291 Input Arguments: 292 + dim - The simplex dimension 293 . order - The quadrature order 294 . a - left end of interval (often-1) 295 - b - right end of interval (often +1) 296 297 Output Arguments: 298 . q - A PetscQuadrature object 299 300 Level: intermediate 301 302 References: 303 Karniadakis and Sherwin. 304 FIAT 305 306 .seealso: PetscDTGaussQuadrature() 307 @*/ 308 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q) 309 { 310 PetscInt npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order; 311 PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w; 312 PetscInt i, j, k; 313 PetscErrorCode ierr; 314 315 PetscFunctionBegin; 316 if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now"); 317 ierr = PetscMalloc(npoints*dim * sizeof(PetscReal), &x);CHKERRQ(ierr); 318 ierr = PetscMalloc(npoints * sizeof(PetscReal), &w);CHKERRQ(ierr); 319 switch (dim) { 320 case 1: 321 ierr = PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, x, w);CHKERRQ(ierr); 322 break; 323 case 2: 324 ierr = PetscMalloc4(order,PetscReal,&px,order,PetscReal,&wx,order,PetscReal,&py,order,PetscReal,&wy);CHKERRQ(ierr); 325 ierr = PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);CHKERRQ(ierr); 326 ierr = PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);CHKERRQ(ierr); 327 for (i = 0; i < order; ++i) { 328 for (j = 0; j < order; ++j) { 329 ierr = PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*order+j)*2+0], &x[(i*order+j)*2+1]);CHKERRQ(ierr); 330 w[i*order+j] = 0.5 * wx[i] * wy[j]; 331 } 332 } 333 ierr = PetscFree4(px,wx,py,wy);CHKERRQ(ierr); 334 break; 335 case 3: 336 ierr = PetscMalloc6(order,PetscReal,&px,order,PetscReal,&wx,order,PetscReal,&py,order,PetscReal,&wy,order,PetscReal,&pz,order,PetscReal,&wz);CHKERRQ(ierr); 337 ierr = PetscDTGaussJacobiQuadrature1D_Internal(order, 0.0, 0.0, px, wx);CHKERRQ(ierr); 338 ierr = PetscDTGaussJacobiQuadrature1D_Internal(order, 1.0, 0.0, py, wy);CHKERRQ(ierr); 339 ierr = PetscDTGaussJacobiQuadrature1D_Internal(order, 2.0, 0.0, pz, wz);CHKERRQ(ierr); 340 for (i = 0; i < order; ++i) { 341 for (j = 0; j < order; ++j) { 342 for (k = 0; k < order; ++k) { 343 ierr = PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*order+j)*order+k)*3+0], &x[((i*order+j)*order+k)*3+1], &x[((i*order+j)*order+k)*3+2]);CHKERRQ(ierr); 344 w[(i*order+j)*order+k] = 0.125 * wx[i] * wy[j] * wz[k]; 345 } 346 } 347 } 348 ierr = PetscFree6(px,wx,py,wy,pz,wz);CHKERRQ(ierr); 349 break; 350 default: 351 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim); 352 } 353 q->numQuadPoints = npoints; 354 q->quadPoints = x; 355 q->quadWeights = w; 356 PetscFunctionReturn(0); 357 } 358 359 #undef __FUNCT__ 360 #define __FUNCT__ "PetscDTPseudoInverseQR" 361 /* Overwrites A. Can only handle full-rank problems with m>=n 362 * A in column-major format 363 * Ainv in row-major format 364 * tau has length m 365 * worksize must be >= max(1,n) 366 */ 367 static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work) 368 { 369 PetscErrorCode ierr; 370 PetscBLASInt M,N,K,lda,ldb,ldwork,info; 371 PetscScalar *A,*Ainv,*R,*Q,Alpha; 372 373 PetscFunctionBegin; 374 #if defined(PETSC_USE_COMPLEX) 375 { 376 PetscInt i,j; 377 ierr = PetscMalloc2(m*n,PetscScalar,&A,m*n,PetscScalar,&Ainv);CHKERRQ(ierr); 378 for (j=0; j<n; j++) { 379 for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j]; 380 } 381 mstride = m; 382 } 383 #else 384 A = A_in; 385 Ainv = Ainv_out; 386 #endif 387 388 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr); 389 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr); 390 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr); 391 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr); 392 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 393 LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info); 394 ierr = PetscFPTrapPop();CHKERRQ(ierr); 395 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error"); 396 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 397 398 /* Extract an explicit representation of Q */ 399 Q = Ainv; 400 ierr = PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));CHKERRQ(ierr); 401 K = N; /* full rank */ 402 LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info); 403 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error"); 404 405 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 406 Alpha = 1.0; 407 ldb = lda; 408 BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb); 409 /* Ainv is Q, overwritten with inverse */ 410 411 #if defined(PETSC_USE_COMPLEX) 412 { 413 PetscInt i; 414 for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]); 415 ierr = PetscFree2(A,Ainv);CHKERRQ(ierr); 416 } 417 #endif 418 PetscFunctionReturn(0); 419 } 420 421 #undef __FUNCT__ 422 #define __FUNCT__ "PetscDTLegendreIntegrate" 423 /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */ 424 static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B) 425 { 426 PetscErrorCode ierr; 427 PetscReal *Bv; 428 PetscInt i,j; 429 430 PetscFunctionBegin; 431 ierr = PetscMalloc((ninterval+1)*ndegree*sizeof(PetscReal),&Bv);CHKERRQ(ierr); 432 /* Point evaluation of L_p on all the source vertices */ 433 ierr = PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);CHKERRQ(ierr); 434 /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */ 435 for (i=0; i<ninterval; i++) { 436 for (j=0; j<ndegree; j++) { 437 if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 438 else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j]; 439 } 440 } 441 ierr = PetscFree(Bv);CHKERRQ(ierr); 442 PetscFunctionReturn(0); 443 } 444 445 #undef __FUNCT__ 446 #define __FUNCT__ "PetscDTReconstructPoly" 447 /*@ 448 PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals 449 450 Not Collective 451 452 Input Arguments: 453 + degree - degree of reconstruction polynomial 454 . nsource - number of source intervals 455 . sourcex - sorted coordinates of source cell boundaries (length nsource+1) 456 . ntarget - number of target intervals 457 - targetx - sorted coordinates of target cell boundaries (length ntarget+1) 458 459 Output Arguments: 460 . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s] 461 462 Level: advanced 463 464 .seealso: PetscDTLegendreEval() 465 @*/ 466 PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R) 467 { 468 PetscErrorCode ierr; 469 PetscInt i,j,k,*bdegrees,worksize; 470 PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget; 471 PetscScalar *tau,*work; 472 473 PetscFunctionBegin; 474 PetscValidRealPointer(sourcex,3); 475 PetscValidRealPointer(targetx,5); 476 PetscValidRealPointer(R,6); 477 if (degree >= nsource) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Reconstruction degree %D must be less than number of source intervals %D",degree,nsource); 478 #if defined(PETSC_USE_DEBUG) 479 for (i=0; i<nsource; i++) { 480 if (sourcex[i] >= sourcex[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Source interval %D has negative orientation (%G,%G)",i,sourcex[i],sourcex[i+1]); 481 } 482 for (i=0; i<ntarget; i++) { 483 if (targetx[i] >= targetx[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Target interval %D has negative orientation (%G,%G)",i,targetx[i],targetx[i+1]); 484 } 485 #endif 486 xmin = PetscMin(sourcex[0],targetx[0]); 487 xmax = PetscMax(sourcex[nsource],targetx[ntarget]); 488 center = (xmin + xmax)/2; 489 hscale = (xmax - xmin)/2; 490 worksize = nsource; 491 ierr = PetscMalloc4(degree+1,PetscInt,&bdegrees,nsource+1,PetscReal,&sourcey,nsource*(degree+1),PetscReal,&Bsource,worksize,PetscScalar,&work);CHKERRQ(ierr); 492 ierr = PetscMalloc4(nsource,PetscScalar,&tau,nsource*(degree+1),PetscReal,&Bsinv,ntarget+1,PetscReal,&targety,ntarget*(degree+1),PetscReal,&Btarget);CHKERRQ(ierr); 493 for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale; 494 for (i=0; i<=degree; i++) bdegrees[i] = i+1; 495 ierr = PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);CHKERRQ(ierr); 496 ierr = PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);CHKERRQ(ierr); 497 for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale; 498 ierr = PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);CHKERRQ(ierr); 499 for (i=0; i<ntarget; i++) { 500 PetscReal rowsum = 0; 501 for (j=0; j<nsource; j++) { 502 PetscReal sum = 0; 503 for (k=0; k<degree+1; k++) { 504 sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j]; 505 } 506 R[i*nsource+j] = sum; 507 rowsum += sum; 508 } 509 for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */ 510 } 511 ierr = PetscFree4(bdegrees,sourcey,Bsource,work);CHKERRQ(ierr); 512 ierr = PetscFree4(tau,Bsinv,targety,Btarget);CHKERRQ(ierr); 513 PetscFunctionReturn(0); 514 } 515 516 /* Basis Jet Tabulation 517 518 We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We 519 follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can 520 be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis 521 as a prime basis. 522 523 \psi_i = \sum_k \alpha_{ki} \phi_k 524 525 Our nodal basis is defined in terms of the dual basis $n_j$ 526 527 n_j \cdot \psi_i = \delta_{ji} 528 529 and we may act on the first equation to obtain 530 531 n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k 532 \delta_{ji} = \sum_k \alpha_{ki} V_{jk} 533 I = V \alpha 534 535 so the coefficients of the nodal basis in the prime basis are 536 537 \alpha = V^{-1} 538 539 We will define the dual basis vectors $n_j$ using a quadrature rule. 540 541 Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials 542 (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can 543 be implemented exactly as in FIAT using functionals $L_j$. 544 545 I will have to count the degrees correctly for the Legendre product when we are on simplices. 546 547 We will have three objects: 548 - Space, P: this just need point evaluation I think 549 - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q 550 - FEM: This keeps {P, P', Q} 551 */ 552 #include <petsc-private/petscfeimpl.h> 553 554 PetscClassId PETSCSPACE_CLASSID = 0; 555 556 PetscFunctionList PetscSpaceList = NULL; 557 PetscBool PetscSpaceRegisterAllCalled = PETSC_FALSE; 558 559 #undef __FUNCT__ 560 #define __FUNCT__ "PetscSpaceRegister" 561 /*@C 562 PetscSpaceRegister - Adds a new PetscSpace implementation 563 564 Not Collective 565 566 Input Parameters: 567 + name - The name of a new user-defined creation routine 568 - create_func - The creation routine itself 569 570 Notes: 571 PetscSpaceRegister() may be called multiple times to add several user-defined PetscSpaces 572 573 Sample usage: 574 .vb 575 PetscSpaceRegister("my_space", MyPetscSpaceCreate); 576 .ve 577 578 Then, your PetscSpace type can be chosen with the procedural interface via 579 .vb 580 PetscSpaceCreate(MPI_Comm, PetscSpace *); 581 PetscSpaceSetType(PetscSpace, "my_space"); 582 .ve 583 or at runtime via the option 584 .vb 585 -petscspace_type my_space 586 .ve 587 588 Level: advanced 589 590 .keywords: PetscSpace, register 591 .seealso: PetscSpaceRegisterAll(), PetscSpaceRegisterDestroy() 592 593 @*/ 594 PetscErrorCode PetscSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscSpace)) 595 { 596 PetscErrorCode ierr; 597 598 PetscFunctionBegin; 599 ierr = PetscFunctionListAdd(&PetscSpaceList, sname, function);CHKERRQ(ierr); 600 PetscFunctionReturn(0); 601 } 602 603 #undef __FUNCT__ 604 #define __FUNCT__ "PetscSpaceSetType" 605 /*@C 606 PetscSpaceSetType - Builds a particular PetscSpace 607 608 Collective on PetscSpace 609 610 Input Parameters: 611 + sp - The PetscSpace object 612 - name - The kind of space 613 614 Options Database Key: 615 . -petscspace_type <type> - Sets the PetscSpace type; use -help for a list of available types 616 617 Level: intermediate 618 619 .keywords: PetscSpace, set, type 620 .seealso: PetscSpaceGetType(), PetscSpaceCreate() 621 @*/ 622 PetscErrorCode PetscSpaceSetType(PetscSpace sp, PetscSpaceType name) 623 { 624 PetscErrorCode (*r)(PetscSpace); 625 PetscBool match; 626 PetscErrorCode ierr; 627 628 PetscFunctionBegin; 629 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 630 ierr = PetscObjectTypeCompare((PetscObject) sp, name, &match);CHKERRQ(ierr); 631 if (match) PetscFunctionReturn(0); 632 633 if (!PetscSpaceRegisterAllCalled) {ierr = PetscSpaceRegisterAll();CHKERRQ(ierr);} 634 ierr = PetscFunctionListFind(PetscSpaceList, name, &r);CHKERRQ(ierr); 635 if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSpace type: %s", name); 636 637 if (sp->ops->destroy) { 638 ierr = (*sp->ops->destroy)(sp);CHKERRQ(ierr); 639 sp->ops->destroy = NULL; 640 } 641 ierr = (*r)(sp);CHKERRQ(ierr); 642 ierr = PetscObjectChangeTypeName((PetscObject) sp, name);CHKERRQ(ierr); 643 PetscFunctionReturn(0); 644 } 645 646 #undef __FUNCT__ 647 #define __FUNCT__ "PetscSpaceGetType" 648 /*@C 649 PetscSpaceGetType - Gets the PetscSpace type name (as a string) from the object. 650 651 Not Collective 652 653 Input Parameter: 654 . dm - The PetscSpace 655 656 Output Parameter: 657 . name - The PetscSpace type name 658 659 Level: intermediate 660 661 .keywords: PetscSpace, get, type, name 662 .seealso: PetscSpaceSetType(), PetscSpaceCreate() 663 @*/ 664 PetscErrorCode PetscSpaceGetType(PetscSpace sp, PetscSpaceType *name) 665 { 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 670 PetscValidCharPointer(name, 2); 671 if (!PetscSpaceRegisterAllCalled) { 672 ierr = PetscSpaceRegisterAll();CHKERRQ(ierr); 673 } 674 *name = ((PetscObject) sp)->type_name; 675 PetscFunctionReturn(0); 676 } 677 678 #undef __FUNCT__ 679 #define __FUNCT__ "PetscSpaceView" 680 /*@C 681 PetscSpaceView - Views a PetscSpace 682 683 Collective on PetscSpace 684 685 Input Parameter: 686 + sp - the PetscSpace object to view 687 - v - the viewer 688 689 Level: developer 690 691 .seealso PetscSpaceDestroy() 692 @*/ 693 PetscErrorCode PetscSpaceView(PetscSpace sp, PetscViewer v) 694 { 695 PetscErrorCode ierr; 696 697 PetscFunctionBegin; 698 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 699 if (!v) { 700 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);CHKERRQ(ierr); 701 } 702 if (sp->ops->view) { 703 ierr = (*sp->ops->view)(sp, v);CHKERRQ(ierr); 704 } 705 PetscFunctionReturn(0); 706 } 707 708 #undef __FUNCT__ 709 #define __FUNCT__ "PetscSpaceViewFromOptions" 710 /* 711 PetscSpaceViewFromOptions - Processes command line options to determine if/how a PetscSpace is to be viewed. 712 713 Collective on PetscSpace 714 715 Input Parameters: 716 + sp - the PetscSpace 717 . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd' 718 - optionname - option to activate viewing 719 720 Level: intermediate 721 722 .keywords: PetscSpace, view, options, database 723 .seealso: VecViewFromOptions(), MatViewFromOptions() 724 */ 725 PetscErrorCode PetscSpaceViewFromOptions(PetscSpace sp, const char prefix[], const char optionname[]) 726 { 727 PetscViewer viewer; 728 PetscViewerFormat format; 729 PetscBool flg; 730 PetscErrorCode ierr; 731 732 PetscFunctionBegin; 733 if (prefix) { 734 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), prefix, optionname, &viewer, &format, &flg);CHKERRQ(ierr); 735 } else { 736 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), ((PetscObject) sp)->prefix, optionname, &viewer, &format, &flg);CHKERRQ(ierr); 737 } 738 if (flg) { 739 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 740 ierr = PetscSpaceView(sp, viewer);CHKERRQ(ierr); 741 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 742 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 743 } 744 PetscFunctionReturn(0); 745 } 746 747 #undef __FUNCT__ 748 #define __FUNCT__ "PetscSpaceSetFromOptions" 749 /*@ 750 PetscSpaceSetFromOptions - sets parameters in a PetscSpace from the options database 751 752 Collective on PetscSpace 753 754 Input Parameter: 755 . sp - the PetscSpace object to set options for 756 757 Options Database: 758 . -petscspace_order the approximation order of the space 759 760 Level: developer 761 762 .seealso PetscSpaceView() 763 @*/ 764 PetscErrorCode PetscSpaceSetFromOptions(PetscSpace sp) 765 { 766 const char *defaultType; 767 char name[256]; 768 PetscBool flg; 769 PetscErrorCode ierr; 770 771 PetscFunctionBegin; 772 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 773 if (!((PetscObject) sp)->type_name) { 774 defaultType = PETSCSPACEPOLYNOMIAL; 775 } else { 776 defaultType = ((PetscObject) sp)->type_name; 777 } 778 if (!PetscSpaceRegisterAllCalled) {ierr = PetscSpaceRegisterAll();CHKERRQ(ierr);} 779 780 ierr = PetscObjectOptionsBegin((PetscObject) sp);CHKERRQ(ierr); 781 ierr = PetscOptionsList("-petscspace_type", "Linear space", "PetscSpaceSetType", PetscSpaceList, defaultType, name, 256, &flg);CHKERRQ(ierr); 782 if (flg) { 783 ierr = PetscSpaceSetType(sp, name);CHKERRQ(ierr); 784 } else if (!((PetscObject) sp)->type_name) { 785 ierr = PetscSpaceSetType(sp, defaultType);CHKERRQ(ierr); 786 } 787 ierr = PetscOptionsInt("-petscspace_order", "The approximation order", "PetscSpaceSetOrder", sp->order, &sp->order, NULL);CHKERRQ(ierr); 788 if (sp->ops->setfromoptions) { 789 ierr = (*sp->ops->setfromoptions)(sp);CHKERRQ(ierr); 790 } 791 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 792 ierr = PetscObjectProcessOptionsHandlers((PetscObject) sp);CHKERRQ(ierr); 793 ierr = PetscOptionsEnd();CHKERRQ(ierr); 794 ierr = PetscSpaceViewFromOptions(sp, NULL, "-petscspace_view");CHKERRQ(ierr); 795 PetscFunctionReturn(0); 796 } 797 798 #undef __FUNCT__ 799 #define __FUNCT__ "PetscSpaceSetUp" 800 /*@C 801 PetscSpaceSetUp - Construct data structures for the PetscSpace 802 803 Collective on PetscSpace 804 805 Input Parameter: 806 . sp - the PetscSpace object to setup 807 808 Level: developer 809 810 .seealso PetscSpaceView(), PetscSpaceDestroy() 811 @*/ 812 PetscErrorCode PetscSpaceSetUp(PetscSpace sp) 813 { 814 PetscErrorCode ierr; 815 816 PetscFunctionBegin; 817 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 818 if (sp->ops->setup) {ierr = (*sp->ops->setup)(sp);CHKERRQ(ierr);} 819 PetscFunctionReturn(0); 820 } 821 822 #undef __FUNCT__ 823 #define __FUNCT__ "PetscSpaceDestroy" 824 /*@ 825 PetscSpaceDestroy - Destroys a PetscSpace object 826 827 Collective on PetscSpace 828 829 Input Parameter: 830 . sp - the PetscSpace object to destroy 831 832 Level: developer 833 834 .seealso PetscSpaceView() 835 @*/ 836 PetscErrorCode PetscSpaceDestroy(PetscSpace *sp) 837 { 838 PetscErrorCode ierr; 839 840 PetscFunctionBegin; 841 if (!*sp) PetscFunctionReturn(0); 842 PetscValidHeaderSpecific((*sp), PETSCSPACE_CLASSID, 1); 843 844 if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; PetscFunctionReturn(0);} 845 ((PetscObject) (*sp))->refct = 0; 846 /* if memory was published with AMS then destroy it */ 847 ierr = PetscObjectAMSViewOff((PetscObject) *sp);CHKERRQ(ierr); 848 849 ierr = DMDestroy(&(*sp)->dm);CHKERRQ(ierr); 850 851 ierr = (*(*sp)->ops->destroy)(*sp);CHKERRQ(ierr); 852 ierr = PetscHeaderDestroy(sp);CHKERRQ(ierr); 853 PetscFunctionReturn(0); 854 } 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "PetscSpaceCreate" 858 /*@ 859 PetscSpaceCreate - Creates an empty PetscSpace object. The type can then be set with PetscSpaceSetType(). 860 861 Collective on MPI_Comm 862 863 Input Parameter: 864 . comm - The communicator for the PetscSpace object 865 866 Output Parameter: 867 . sp - The PetscSpace object 868 869 Level: beginner 870 871 .seealso: PetscSpaceSetType(), PETSCSPACEPOLYNOMIAL 872 @*/ 873 PetscErrorCode PetscSpaceCreate(MPI_Comm comm, PetscSpace *sp) 874 { 875 PetscSpace s; 876 PetscErrorCode ierr; 877 878 PetscFunctionBegin; 879 PetscValidPointer(sp, 2); 880 *sp = NULL; 881 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 882 ierr = PetscFEInitializePackage();CHKERRQ(ierr); 883 #endif 884 885 ierr = PetscHeaderCreate(s, _p_PetscSpace, struct _PetscSpaceOps, PETSCSPACE_CLASSID, "PetscSpace", "Linear Space", "PetscSpace", comm, PetscSpaceDestroy, PetscSpaceView);CHKERRQ(ierr); 886 ierr = PetscMemzero(s->ops, sizeof(struct _PetscSpaceOps));CHKERRQ(ierr); 887 888 s->order = 0; 889 ierr = DMShellCreate(comm, &s->dm);CHKERRQ(ierr); 890 891 *sp = s; 892 PetscFunctionReturn(0); 893 } 894 895 #undef __FUNCT__ 896 #define __FUNCT__ "PetscSpaceGetDimension" 897 /* Dimension of the space, i.e. number of basis vectors */ 898 PetscErrorCode PetscSpaceGetDimension(PetscSpace sp, PetscInt *dim) 899 { 900 PetscErrorCode ierr; 901 902 PetscFunctionBegin; 903 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 904 PetscValidPointer(dim, 2); 905 *dim = 0; 906 if (sp->ops->getdimension) {ierr = (*sp->ops->getdimension)(sp, dim);CHKERRQ(ierr);} 907 PetscFunctionReturn(0); 908 } 909 910 #undef __FUNCT__ 911 #define __FUNCT__ "PetscSpaceGetOrder" 912 PetscErrorCode PetscSpaceGetOrder(PetscSpace sp, PetscInt *order) 913 { 914 PetscFunctionBegin; 915 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 916 PetscValidPointer(order, 2); 917 *order = sp->order; 918 PetscFunctionReturn(0); 919 } 920 921 #undef __FUNCT__ 922 #define __FUNCT__ "PetscSpaceSetOrder" 923 PetscErrorCode PetscSpaceSetOrder(PetscSpace sp, PetscInt order) 924 { 925 PetscFunctionBegin; 926 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 927 sp->order = order; 928 PetscFunctionReturn(0); 929 } 930 931 #undef __FUNCT__ 932 #define __FUNCT__ "PetscSpaceEvaluate" 933 PetscErrorCode PetscSpaceEvaluate(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 934 { 935 PetscErrorCode ierr; 936 937 PetscFunctionBegin; 938 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 939 PetscValidPointer(points, 3); 940 if (B) PetscValidPointer(B, 4); 941 if (D) PetscValidPointer(D, 5); 942 if (H) PetscValidPointer(H, 6); 943 if (sp->ops->evaluate) {ierr = (*sp->ops->evaluate)(sp, npoints, points, B, D, H);CHKERRQ(ierr);} 944 PetscFunctionReturn(0); 945 } 946 947 #undef __FUNCT__ 948 #define __FUNCT__ "PetscSpaceSetFromOptions_Polynomial" 949 PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscSpace sp) 950 { 951 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 952 PetscErrorCode ierr; 953 954 PetscFunctionBegin; 955 ierr = PetscObjectOptionsBegin((PetscObject) sp);CHKERRQ(ierr); 956 ierr = PetscOptionsInt("-petscspace_num_variables", "The number of different variables, e.g. x and y", "PetscSpacePolynomialSetNumVariables", poly->numVariables, &poly->numVariables, NULL);CHKERRQ(ierr); 957 ierr = PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);CHKERRQ(ierr); 958 ierr = PetscOptionsEnd();CHKERRQ(ierr); 959 PetscFunctionReturn(0); 960 } 961 962 #undef __FUNCT__ 963 #define __FUNCT__ "PetscSpacePolynomialView_Ascii" 964 PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer viewer) 965 { 966 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 967 PetscViewerFormat format; 968 PetscErrorCode ierr; 969 970 PetscFunctionBegin; 971 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 972 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 973 ierr = PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d", poly->numVariables, sp->order);CHKERRQ(ierr); 974 } else { 975 ierr = PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d", poly->numVariables, sp->order);CHKERRQ(ierr); 976 } 977 PetscFunctionReturn(0); 978 } 979 980 #undef __FUNCT__ 981 #define __FUNCT__ "PetscSpaceView_Polynomial" 982 PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer) 983 { 984 PetscBool iascii; 985 PetscErrorCode ierr; 986 987 PetscFunctionBegin; 988 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 989 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 990 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 991 if (iascii) {ierr = PetscSpacePolynomialView_Ascii(sp, viewer);CHKERRQ(ierr);} 992 PetscFunctionReturn(0); 993 } 994 995 #undef __FUNCT__ 996 #define __FUNCT__ "PetscSpaceSetUp_Polynomial" 997 PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp) 998 { 999 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1000 PetscInt ndegree = sp->order+1; 1001 PetscInt deg; 1002 PetscErrorCode ierr; 1003 1004 PetscFunctionBegin; 1005 ierr = PetscMalloc(ndegree * sizeof(PetscInt), &poly->degrees);CHKERRQ(ierr); 1006 for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg; 1007 PetscFunctionReturn(0); 1008 } 1009 1010 #undef __FUNCT__ 1011 #define __FUNCT__ "PetscSpaceDestroy_Polynomial" 1012 PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp) 1013 { 1014 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1015 PetscErrorCode ierr; 1016 1017 PetscFunctionBegin; 1018 ierr = PetscFree(poly->degrees);CHKERRQ(ierr); 1019 ierr = PetscFree(poly);CHKERRQ(ierr); 1020 PetscFunctionReturn(0); 1021 } 1022 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "PetscSpaceGetDimension_Polynomial" 1025 PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim) 1026 { 1027 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1028 PetscInt deg = sp->order; 1029 PetscInt n = poly->numVariables, i; 1030 PetscReal D = 1.0; 1031 1032 PetscFunctionBegin; 1033 for (i = 1; i <= n; ++i) { 1034 D *= ((PetscReal) (deg+i))/i; 1035 } 1036 *dim = (PetscInt) (D + 0.5); 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "LatticePoint_Internal" 1042 /* 1043 LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'. 1044 1045 Input Parameters: 1046 + len - The length of the tuple 1047 . sum - The sum of all entries in the tuple 1048 - ind - The current multi-index of the tuple, initialized to the 0 tuple 1049 1050 Output Parameter: 1051 + ind - The multi-index of the tuple, -1 indicates the iteration has terminated 1052 . tup - A tuple of len integers addig to sum 1053 1054 Level: developer 1055 1056 .seealso: 1057 */ 1058 static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[]) 1059 { 1060 PetscInt i; 1061 PetscErrorCode ierr; 1062 1063 PetscFunctionBegin; 1064 if (len == 1) { 1065 ind[0] = -1; 1066 tup[0] = sum; 1067 } else if (sum == 0) { 1068 for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;} 1069 } else { 1070 tup[0] = sum - ind[0]; 1071 ierr = LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);CHKERRQ(ierr); 1072 if (ind[1] < 0) { 1073 if (ind[0] == sum) {ind[0] = -1;} 1074 else {ind[1] = 0; ++ind[0];} 1075 } 1076 } 1077 PetscFunctionReturn(0); 1078 } 1079 1080 #undef __FUNCT__ 1081 #define __FUNCT__ "PetscSpaceEvaluate_Polynomial" 1082 PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) 1083 { 1084 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1085 DM dm = sp->dm; 1086 PetscInt ndegree = sp->order+1; 1087 PetscInt *degrees = poly->degrees; 1088 PetscInt dim = poly->numVariables; 1089 PetscReal *lpoints, *tmp, *LB, *LD, *LH; 1090 PetscInt *ind, *tup; 1091 PetscInt pdim, d, der, i, p, deg, o; 1092 PetscErrorCode ierr; 1093 1094 PetscFunctionBegin; 1095 ierr = PetscSpaceGetDimension(sp, &pdim);CHKERRQ(ierr); 1096 ierr = DMGetWorkArray(dm, npoints, PETSC_REAL, &lpoints);CHKERRQ(ierr); 1097 ierr = DMGetWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);CHKERRQ(ierr); 1098 if (B) {ierr = DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);CHKERRQ(ierr);} 1099 if (D) {ierr = DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);CHKERRQ(ierr);} 1100 if (H) {ierr = DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);CHKERRQ(ierr);} 1101 for (d = 0; d < dim; ++d) { 1102 for (p = 0; p < npoints; ++p) { 1103 lpoints[p] = points[p*dim+d]; 1104 } 1105 ierr = PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);CHKERRQ(ierr); 1106 /* LB, LD, LH (ndegree * dim x npoints) */ 1107 for (deg = 0; deg < ndegree; ++deg) { 1108 for (p = 0; p < npoints; ++p) { 1109 if (B) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg]; 1110 if (D) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg]; 1111 if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg]; 1112 } 1113 } 1114 } 1115 /* Multiply by A (pdim x ndegree * dim) */ 1116 ierr = PetscMalloc2(dim,PetscInt,&ind,dim,PetscInt,&tup);CHKERRQ(ierr); 1117 if (B) { 1118 /* B (npoints x pdim) */ 1119 i = 0; 1120 for (o = 0; o <= sp->order; ++o) { 1121 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 1122 while (ind[0] >= 0) { 1123 ierr = LatticePoint_Internal(dim, o, ind, tup);CHKERRQ(ierr); 1124 for (p = 0; p < npoints; ++p) { 1125 B[p*pdim + i] = 1.0; 1126 for (d = 0; d < dim; ++d) { 1127 B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p]; 1128 } 1129 } 1130 ++i; 1131 } 1132 } 1133 } 1134 if (D) { 1135 /* D (npoints x pdim x dim) */ 1136 i = 0; 1137 for (o = 0; o <= sp->order; ++o) { 1138 ierr = PetscMemzero(ind, dim * sizeof(PetscInt));CHKERRQ(ierr); 1139 while (ind[0] >= 0) { 1140 ierr = LatticePoint_Internal(dim, o, ind, tup);CHKERRQ(ierr); 1141 for (p = 0; p < npoints; ++p) { 1142 for (der = 0; der < dim; ++der) { 1143 D[(p*pdim + i)*dim + der] = 1.0; 1144 for (d = 0; d < dim; ++d) { 1145 if (d == der) { 1146 D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p]; 1147 } else { 1148 D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p]; 1149 } 1150 } 1151 } 1152 } 1153 ++i; 1154 } 1155 } 1156 } 1157 if (H) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to code second derivatives"); 1158 ierr = PetscFree2(ind,tup);CHKERRQ(ierr); 1159 if (B) {ierr = DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);CHKERRQ(ierr);} 1160 if (D) {ierr = DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);CHKERRQ(ierr);} 1161 if (H) {ierr = DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);CHKERRQ(ierr);} 1162 ierr = DMRestoreWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);CHKERRQ(ierr); 1163 ierr = DMRestoreWorkArray(dm, npoints, PETSC_REAL, &lpoints);CHKERRQ(ierr); 1164 PetscFunctionReturn(0); 1165 } 1166 1167 #undef __FUNCT__ 1168 #define __FUNCT__ "PetscSpaceInitialize_Polynomial" 1169 PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp) 1170 { 1171 PetscFunctionBegin; 1172 sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial; 1173 sp->ops->setup = PetscSpaceSetUp_Polynomial; 1174 sp->ops->view = PetscSpaceView_Polynomial; 1175 sp->ops->destroy = PetscSpaceDestroy_Polynomial; 1176 sp->ops->getdimension = PetscSpaceGetDimension_Polynomial; 1177 sp->ops->evaluate = PetscSpaceEvaluate_Polynomial; 1178 PetscFunctionReturn(0); 1179 } 1180 1181 /*MC 1182 PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of linear polynomials. 1183 1184 Level: intermediate 1185 1186 .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType() 1187 M*/ 1188 1189 #undef __FUNCT__ 1190 #define __FUNCT__ "PetscSpaceCreate_Polynomial" 1191 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp) 1192 { 1193 PetscSpace_Poly *poly; 1194 PetscErrorCode ierr; 1195 1196 PetscFunctionBegin; 1197 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 1198 ierr = PetscNewLog(sp, PetscSpace_Poly, &poly);CHKERRQ(ierr); 1199 sp->data = poly; 1200 1201 poly->numVariables = 0; 1202 poly->symmetric = PETSC_FALSE; 1203 poly->degrees = NULL; 1204 1205 ierr = PetscSpaceInitialize_Polynomial(sp);CHKERRQ(ierr); 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "PetscSpacePolynomialSetSymmetric" 1211 PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym) 1212 { 1213 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1214 1215 PetscFunctionBegin; 1216 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 1217 poly->symmetric = sym; 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "PetscSpacePolynomialGetSymmetric" 1223 PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym) 1224 { 1225 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1226 1227 PetscFunctionBegin; 1228 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 1229 PetscValidPointer(sym, 2); 1230 *sym = poly->symmetric; 1231 PetscFunctionReturn(0); 1232 } 1233 1234 #undef __FUNCT__ 1235 #define __FUNCT__ "PetscSpacePolynomialSetNumVariables" 1236 PetscErrorCode PetscSpacePolynomialSetNumVariables(PetscSpace sp, PetscInt n) 1237 { 1238 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1239 1240 PetscFunctionBegin; 1241 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 1242 poly->numVariables = n; 1243 PetscFunctionReturn(0); 1244 } 1245 1246 #undef __FUNCT__ 1247 #define __FUNCT__ "PetscSpacePolynomialGetNumVariables" 1248 PetscErrorCode PetscSpacePolynomialGetNumVariables(PetscSpace sp, PetscInt *n) 1249 { 1250 PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data; 1251 1252 PetscFunctionBegin; 1253 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1); 1254 PetscValidPointer(n, 2); 1255 *n = poly->numVariables; 1256 PetscFunctionReturn(0); 1257 } 1258 1259 1260 PetscClassId PETSCDUALSPACE_CLASSID = 0; 1261 1262 PetscFunctionList PetscDualSpaceList = NULL; 1263 PetscBool PetscDualSpaceRegisterAllCalled = PETSC_FALSE; 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "PetscDualSpaceRegister" 1267 /*@C 1268 PetscDualSpaceRegister - Adds a new PetscDualSpace implementation 1269 1270 Not Collective 1271 1272 Input Parameters: 1273 + name - The name of a new user-defined creation routine 1274 - create_func - The creation routine itself 1275 1276 Notes: 1277 PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces 1278 1279 Sample usage: 1280 .vb 1281 PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate); 1282 .ve 1283 1284 Then, your PetscDualSpace type can be chosen with the procedural interface via 1285 .vb 1286 PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *); 1287 PetscDualSpaceSetType(PetscDualSpace, "my_dual_space"); 1288 .ve 1289 or at runtime via the option 1290 .vb 1291 -petscdualspace_type my_dual_space 1292 .ve 1293 1294 Level: advanced 1295 1296 .keywords: PetscDualSpace, register 1297 .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy() 1298 1299 @*/ 1300 PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace)) 1301 { 1302 PetscErrorCode ierr; 1303 1304 PetscFunctionBegin; 1305 ierr = PetscFunctionListAdd(&PetscDualSpaceList, sname, function);CHKERRQ(ierr); 1306 PetscFunctionReturn(0); 1307 } 1308 1309 #undef __FUNCT__ 1310 #define __FUNCT__ "PetscDualSpaceSetType" 1311 /*@C 1312 PetscDualSpaceSetType - Builds a particular PetscDualSpace 1313 1314 Collective on PetscDualSpace 1315 1316 Input Parameters: 1317 + sp - The PetscDualSpace object 1318 - name - The kind of space 1319 1320 Options Database Key: 1321 . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types 1322 1323 Level: intermediate 1324 1325 .keywords: PetscDualSpace, set, type 1326 .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate() 1327 @*/ 1328 PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name) 1329 { 1330 PetscErrorCode (*r)(PetscDualSpace); 1331 PetscBool match; 1332 PetscErrorCode ierr; 1333 1334 PetscFunctionBegin; 1335 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1336 ierr = PetscObjectTypeCompare((PetscObject) sp, name, &match);CHKERRQ(ierr); 1337 if (match) PetscFunctionReturn(0); 1338 1339 if (!PetscDualSpaceRegisterAllCalled) {ierr = PetscDualSpaceRegisterAll();CHKERRQ(ierr);} 1340 ierr = PetscFunctionListFind(PetscDualSpaceList, name, &r);CHKERRQ(ierr); 1341 if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name); 1342 1343 if (sp->ops->destroy) { 1344 ierr = (*sp->ops->destroy)(sp);CHKERRQ(ierr); 1345 sp->ops->destroy = NULL; 1346 } 1347 ierr = (*r)(sp);CHKERRQ(ierr); 1348 ierr = PetscObjectChangeTypeName((PetscObject) sp, name);CHKERRQ(ierr); 1349 PetscFunctionReturn(0); 1350 } 1351 1352 #undef __FUNCT__ 1353 #define __FUNCT__ "PetscDualSpaceGetType" 1354 /*@C 1355 PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object. 1356 1357 Not Collective 1358 1359 Input Parameter: 1360 . dm - The PetscDualSpace 1361 1362 Output Parameter: 1363 . name - The PetscDualSpace type name 1364 1365 Level: intermediate 1366 1367 .keywords: PetscDualSpace, get, type, name 1368 .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate() 1369 @*/ 1370 PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name) 1371 { 1372 PetscErrorCode ierr; 1373 1374 PetscFunctionBegin; 1375 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1376 PetscValidCharPointer(name, 2); 1377 if (!PetscDualSpaceRegisterAllCalled) { 1378 ierr = PetscDualSpaceRegisterAll();CHKERRQ(ierr); 1379 } 1380 *name = ((PetscObject) sp)->type_name; 1381 PetscFunctionReturn(0); 1382 } 1383 1384 #undef __FUNCT__ 1385 #define __FUNCT__ "PetscDualSpaceView" 1386 /*@C 1387 PetscDualSpaceView - Views a PetscDualSpace 1388 1389 Collective on PetscDualSpace 1390 1391 Input Parameter: 1392 + sp - the PetscDualSpace object to view 1393 - v - the viewer 1394 1395 Level: developer 1396 1397 .seealso PetscDualSpaceDestroy() 1398 @*/ 1399 PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v) 1400 { 1401 PetscErrorCode ierr; 1402 1403 PetscFunctionBegin; 1404 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1405 if (!v) { 1406 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);CHKERRQ(ierr); 1407 } 1408 if (sp->ops->view) { 1409 ierr = (*sp->ops->view)(sp, v);CHKERRQ(ierr); 1410 } 1411 PetscFunctionReturn(0); 1412 } 1413 1414 #undef __FUNCT__ 1415 #define __FUNCT__ "PetscDualSpaceViewFromOptions" 1416 /* 1417 PetscDualSpaceViewFromOptions - Processes command line options to determine if/how a PetscDualSpace is to be viewed. 1418 1419 Collective on PetscDualSpace 1420 1421 Input Parameters: 1422 + sp - the PetscDualSpace 1423 . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd' 1424 - optionname - option to activate viewing 1425 1426 Level: intermediate 1427 1428 .keywords: PetscDualSpace, view, options, database 1429 .seealso: VecViewFromOptions(), MatViewFromOptions() 1430 */ 1431 PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace sp, const char prefix[], const char optionname[]) 1432 { 1433 PetscViewer viewer; 1434 PetscViewerFormat format; 1435 PetscBool flg; 1436 PetscErrorCode ierr; 1437 1438 PetscFunctionBegin; 1439 if (prefix) { 1440 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), prefix, optionname, &viewer, &format, &flg);CHKERRQ(ierr); 1441 } else { 1442 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), ((PetscObject) sp)->prefix, optionname, &viewer, &format, &flg);CHKERRQ(ierr); 1443 } 1444 if (flg) { 1445 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 1446 ierr = PetscDualSpaceView(sp, viewer);CHKERRQ(ierr); 1447 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 1448 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1449 } 1450 PetscFunctionReturn(0); 1451 } 1452 1453 #undef __FUNCT__ 1454 #define __FUNCT__ "PetscDualSpaceSetFromOptions" 1455 /*@ 1456 PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database 1457 1458 Collective on PetscDualSpace 1459 1460 Input Parameter: 1461 . sp - the PetscDualSpace object to set options for 1462 1463 Options Database: 1464 . -petscspace_order the approximation order of the space 1465 1466 Level: developer 1467 1468 .seealso PetscDualSpaceView() 1469 @*/ 1470 PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp) 1471 { 1472 const char *defaultType; 1473 char name[256]; 1474 PetscBool flg; 1475 PetscErrorCode ierr; 1476 1477 PetscFunctionBegin; 1478 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1479 if (!((PetscObject) sp)->type_name) { 1480 defaultType = PETSCDUALSPACELAGRANGE; 1481 } else { 1482 defaultType = ((PetscObject) sp)->type_name; 1483 } 1484 if (!PetscSpaceRegisterAllCalled) {ierr = PetscSpaceRegisterAll();CHKERRQ(ierr);} 1485 1486 ierr = PetscObjectOptionsBegin((PetscObject) sp);CHKERRQ(ierr); 1487 ierr = PetscOptionsList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);CHKERRQ(ierr); 1488 if (flg) { 1489 ierr = PetscDualSpaceSetType(sp, name);CHKERRQ(ierr); 1490 } else if (!((PetscObject) sp)->type_name) { 1491 ierr = PetscDualSpaceSetType(sp, defaultType);CHKERRQ(ierr); 1492 } 1493 ierr = PetscOptionsInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);CHKERRQ(ierr); 1494 if (sp->ops->setfromoptions) { 1495 ierr = (*sp->ops->setfromoptions)(sp);CHKERRQ(ierr); 1496 } 1497 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 1498 ierr = PetscObjectProcessOptionsHandlers((PetscObject) sp);CHKERRQ(ierr); 1499 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1500 ierr = PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");CHKERRQ(ierr); 1501 PetscFunctionReturn(0); 1502 } 1503 1504 #undef __FUNCT__ 1505 #define __FUNCT__ "PetscDualSpaceSetUp" 1506 /*@C 1507 PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace 1508 1509 Collective on PetscDualSpace 1510 1511 Input Parameter: 1512 . sp - the PetscDualSpace object to setup 1513 1514 Level: developer 1515 1516 .seealso PetscDualSpaceView(), PetscDualSpaceDestroy() 1517 @*/ 1518 PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp) 1519 { 1520 PetscErrorCode ierr; 1521 1522 PetscFunctionBegin; 1523 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1524 if (sp->ops->setup) {ierr = (*sp->ops->setup)(sp);CHKERRQ(ierr);} 1525 PetscFunctionReturn(0); 1526 } 1527 1528 #undef __FUNCT__ 1529 #define __FUNCT__ "PetscDualSpaceDestroy" 1530 /*@ 1531 PetscDualSpaceDestroy - Destroys a PetscDualSpace object 1532 1533 Collective on PetscDualSpace 1534 1535 Input Parameter: 1536 . sp - the PetscDualSpace object to destroy 1537 1538 Level: developer 1539 1540 .seealso PetscDualSpaceView() 1541 @*/ 1542 PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp) 1543 { 1544 PetscInt dim, f; 1545 PetscErrorCode ierr; 1546 1547 PetscFunctionBegin; 1548 if (!*sp) PetscFunctionReturn(0); 1549 PetscValidHeaderSpecific((*sp), PETSCDUALSPACE_CLASSID, 1); 1550 1551 if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; PetscFunctionReturn(0);} 1552 ((PetscObject) (*sp))->refct = 0; 1553 /* if memory was published with AMS then destroy it */ 1554 ierr = PetscObjectAMSViewOff((PetscObject) *sp);CHKERRQ(ierr); 1555 1556 ierr = PetscDualSpaceGetDimension(*sp, &dim);CHKERRQ(ierr); 1557 for (f = 0; f < dim; ++f) { 1558 ierr = PetscQuadratureDestroy(&(*sp)->functional[f]);CHKERRQ(ierr); 1559 } 1560 ierr = PetscFree((*sp)->functional);CHKERRQ(ierr); 1561 ierr = DMDestroy(&(*sp)->dm);CHKERRQ(ierr); 1562 1563 if ((*sp)->ops->destroy) {ierr = (*(*sp)->ops->destroy)(*sp);CHKERRQ(ierr);} 1564 ierr = PetscHeaderDestroy(sp);CHKERRQ(ierr); 1565 PetscFunctionReturn(0); 1566 } 1567 1568 #undef __FUNCT__ 1569 #define __FUNCT__ "PetscDualSpaceCreate" 1570 /*@ 1571 PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType(). 1572 1573 Collective on MPI_Comm 1574 1575 Input Parameter: 1576 . comm - The communicator for the PetscDualSpace object 1577 1578 Output Parameter: 1579 . sp - The PetscDualSpace object 1580 1581 Level: beginner 1582 1583 .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE 1584 @*/ 1585 PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp) 1586 { 1587 PetscDualSpace s; 1588 PetscErrorCode ierr; 1589 1590 PetscFunctionBegin; 1591 PetscValidPointer(sp, 2); 1592 *sp = NULL; 1593 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1594 ierr = PetscFEInitializePackage();CHKERRQ(ierr); 1595 #endif 1596 1597 ierr = PetscHeaderCreate(s, _p_PetscDualSpace, struct _PetscDualSpaceOps, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);CHKERRQ(ierr); 1598 ierr = PetscMemzero(s->ops, sizeof(struct _PetscDualSpaceOps));CHKERRQ(ierr); 1599 1600 s->order = 0; 1601 1602 *sp = s; 1603 PetscFunctionReturn(0); 1604 } 1605 1606 #undef __FUNCT__ 1607 #define __FUNCT__ "PetscDualSpaceGetDM" 1608 PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm) 1609 { 1610 PetscFunctionBegin; 1611 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1612 PetscValidPointer(dm, 2); 1613 *dm = sp->dm; 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "PetscDualSpaceSetDM" 1619 PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm) 1620 { 1621 PetscErrorCode ierr; 1622 1623 PetscFunctionBegin; 1624 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1625 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1626 ierr = DMDestroy(&sp->dm);CHKERRQ(ierr); 1627 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1628 sp->dm = dm; 1629 PetscFunctionReturn(0); 1630 } 1631 1632 #undef __FUNCT__ 1633 #define __FUNCT__ "PetscDualSpaceGetOrder" 1634 PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order) 1635 { 1636 PetscFunctionBegin; 1637 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1638 PetscValidPointer(order, 2); 1639 *order = sp->order; 1640 PetscFunctionReturn(0); 1641 } 1642 1643 #undef __FUNCT__ 1644 #define __FUNCT__ "PetscDualSpaceSetOrder" 1645 PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order) 1646 { 1647 PetscFunctionBegin; 1648 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1649 sp->order = order; 1650 PetscFunctionReturn(0); 1651 } 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "PetscDualSpaceGetFunctional" 1655 PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional) 1656 { 1657 PetscInt dim; 1658 PetscErrorCode ierr; 1659 1660 PetscFunctionBegin; 1661 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1662 PetscValidPointer(functional, 3); 1663 ierr = PetscDualSpaceGetDimension(sp, &dim);CHKERRQ(ierr); 1664 if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim); 1665 *functional = sp->functional[i]; 1666 PetscFunctionReturn(0); 1667 } 1668 1669 #undef __FUNCT__ 1670 #define __FUNCT__ "PetscDualSpaceGetDimension" 1671 /* Dimension of the space, i.e. number of basis vectors */ 1672 PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim) 1673 { 1674 PetscErrorCode ierr; 1675 1676 PetscFunctionBegin; 1677 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1678 PetscValidPointer(dim, 2); 1679 *dim = 0; 1680 if (sp->ops->getdimension) {ierr = (*sp->ops->getdimension)(sp, dim);CHKERRQ(ierr);} 1681 PetscFunctionReturn(0); 1682 } 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "PetscDualSpaceGetNumDof" 1686 PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof) 1687 { 1688 PetscErrorCode ierr; 1689 1690 PetscFunctionBegin; 1691 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1692 PetscValidPointer(numDof, 2); 1693 *numDof = NULL; 1694 if (sp->ops->getnumdof) {ierr = (*sp->ops->getnumdof)(sp, numDof);CHKERRQ(ierr);} 1695 PetscFunctionReturn(0); 1696 } 1697 1698 #undef __FUNCT__ 1699 #define __FUNCT__ "PetscDualSpaceCreateReferenceCell" 1700 PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm) 1701 { 1702 DM rdm; 1703 PetscErrorCode ierr; 1704 1705 PetscFunctionBeginUser; 1706 ierr = DMCreate(PetscObjectComm((PetscObject) sp), &rdm);CHKERRQ(ierr); 1707 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 1708 ierr = DMPlexSetDimension(rdm, dim);CHKERRQ(ierr); 1709 switch (dim) { 1710 case 2: 1711 { 1712 PetscInt numPoints[2] = {3, 1}; 1713 PetscInt coneSize[4] = {3, 0, 0, 0}; 1714 PetscInt cones[3] = {1, 2, 3}; 1715 PetscInt coneOrientations[3] = {0, 0, 0}; 1716 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 1717 1718 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1719 } 1720 break; 1721 case 3: 1722 { 1723 PetscInt numPoints[2] = {4, 1}; 1724 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 1725 PetscInt cones[4] = {1, 2, 3, 4}; 1726 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 1727 PetscScalar vertexCoords[12] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0}; 1728 1729 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1730 } 1731 break; 1732 default: 1733 SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 1734 } 1735 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 1736 ierr = DMPlexCopyCoordinates(rdm, *refdm);CHKERRQ(ierr); 1737 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 1738 PetscFunctionReturn(0); 1739 } 1740 1741 #undef __FUNCT__ 1742 #define __FUNCT__ "PetscDualSpaceSetUp_Lagrange" 1743 PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp) 1744 { 1745 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 1746 DM dm = sp->dm; 1747 PetscInt order = sp->order; 1748 PetscSection csection; 1749 Vec coordinates; 1750 PetscReal *qpoints, *qweights; 1751 PetscInt depth, dim, pdim, *pStart, *pEnd, coneSize, d, n, f = 0; 1752 PetscErrorCode ierr; 1753 1754 PetscFunctionBegin; 1755 ierr = PetscDualSpaceGetDimension(sp, &pdim);CHKERRQ(ierr); 1756 ierr = PetscMalloc(pdim * sizeof(PetscQuadrature), &sp->functional);CHKERRQ(ierr); 1757 /* Classify element type */ 1758 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1759 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1760 ierr = PetscMalloc((dim+1) * sizeof(PetscInt), &lag->numDof);CHKERRQ(ierr); 1761 ierr = PetscMemzero(lag->numDof, (dim+1) * sizeof(PetscInt));CHKERRQ(ierr); 1762 ierr = PetscMalloc2(depth+1,PetscInt,&pStart,depth+1,PetscInt,&pEnd);CHKERRQ(ierr); 1763 for (d = 0; d <= depth; ++d) { 1764 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1765 } 1766 ierr = DMPlexGetConeSize(dm, pStart[depth], &coneSize);CHKERRQ(ierr); 1767 ierr = DMPlexGetCoordinateSection(dm, &csection);CHKERRQ(ierr); 1768 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1769 if (coneSize == dim+1) { 1770 PetscInt *closure = NULL, closureSize, c; 1771 1772 /* Simplex */ 1773 ierr = DMPlexGetTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1774 for (c = 0; c < closureSize*2; c += 2) { 1775 const PetscInt p = closure[c]; 1776 1777 if ((p >= pStart[0]) && (p < pEnd[0])) { 1778 /* Vertices */ 1779 const PetscScalar *coords; 1780 PetscInt dof, off, d; 1781 1782 if (order < 1) continue; 1783 sp->functional[f].numQuadPoints = 1; 1784 ierr = PetscMalloc(sp->functional[f].numQuadPoints*dim * sizeof(PetscReal), &qpoints);CHKERRQ(ierr); 1785 ierr = PetscMalloc(sp->functional[f].numQuadPoints * sizeof(PetscReal), &qweights);CHKERRQ(ierr); 1786 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 1787 ierr = PetscSectionGetDof(csection, p, &dof);CHKERRQ(ierr); 1788 ierr = PetscSectionGetOffset(csection, p, &off);CHKERRQ(ierr); 1789 if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim); 1790 for (d = 0; d < dof; ++d) {qpoints[d] = coords[off+d];} 1791 qweights[0] = 1.0; 1792 sp->functional[f].quadPoints = qpoints; 1793 sp->functional[f].quadWeights = qweights; 1794 ++f; 1795 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 1796 lag->numDof[0] = 1; 1797 } else if ((p >= pStart[1]) && (p < pEnd[1])) { 1798 /* Edges */ 1799 PetscScalar *coords; 1800 PetscInt k; 1801 1802 if (order < 2) continue; 1803 coords = NULL; 1804 ierr = DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);CHKERRQ(ierr); 1805 if (n != dim*2) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d has %d coordinate values instead of %d", p, n, dim*2); 1806 for (k = 1; k < order; ++k) { 1807 sp->functional[f].numQuadPoints = 1; 1808 ierr = PetscMalloc(sp->functional[f].numQuadPoints*dim * sizeof(PetscReal), &qpoints);CHKERRQ(ierr); 1809 ierr = PetscMalloc(sp->functional[f].numQuadPoints * sizeof(PetscReal), &qweights);CHKERRQ(ierr); 1810 for (d = 0; d < dim; ++d) {qpoints[d] = k*(coords[1*dim+d] - coords[0*dim+d])/order + coords[0*dim+d];} 1811 qweights[0] = 1.0; 1812 sp->functional[f].quadPoints = qpoints; 1813 sp->functional[f].quadWeights = qweights; 1814 ++f; 1815 } 1816 ierr = DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);CHKERRQ(ierr); 1817 lag->numDof[1] = order-1; 1818 } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) { 1819 /* Faces */ 1820 1821 if (order < 3) continue; 1822 lag->numDof[depth-1] = 0; 1823 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces"); 1824 } else if ((p >= pStart[depth]) && (p < pEnd[depth])) { 1825 /* Cells */ 1826 1827 if ((order > 0) && (order < 3)) continue; 1828 lag->numDof[depth] = 0; 1829 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement cells"); 1830 } 1831 } 1832 ierr = DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1833 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot handle cells with cone size %d", coneSize); 1834 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 1835 if (f != pdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vector %d not equal to dimension %d", f, pdim); 1836 PetscFunctionReturn(0); 1837 } 1838 1839 #undef __FUNCT__ 1840 #define __FUNCT__ "PetscDualSpaceDestroy_Lagrange" 1841 PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp) 1842 { 1843 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 1844 PetscErrorCode ierr; 1845 1846 PetscFunctionBegin; 1847 ierr = PetscFree(lag->numDof);CHKERRQ(ierr); 1848 ierr = PetscFree(lag);CHKERRQ(ierr); 1849 PetscFunctionReturn(0); 1850 } 1851 1852 #undef __FUNCT__ 1853 #define __FUNCT__ "PetscDualSpaceGetDimension_Lagrange" 1854 PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim) 1855 { 1856 PetscInt deg = sp->order; 1857 PetscReal D = 1.0; 1858 PetscInt n, i; 1859 PetscErrorCode ierr; 1860 1861 PetscFunctionBegin; 1862 /* TODO: Assumes simplices */ 1863 ierr = DMPlexGetDimension(sp->dm, &n);CHKERRQ(ierr); 1864 for (i = 1; i <= n; ++i) { 1865 D *= ((PetscReal) (deg+i))/i; 1866 } 1867 *dim = (PetscInt) (D + 0.5); 1868 PetscFunctionReturn(0); 1869 } 1870 1871 #undef __FUNCT__ 1872 #define __FUNCT__ "PetscDualSpaceGetNumDof_Lagrange" 1873 PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof) 1874 { 1875 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 1876 1877 PetscFunctionBegin; 1878 *numDof = lag->numDof; 1879 PetscFunctionReturn(0); 1880 } 1881 1882 #undef __FUNCT__ 1883 #define __FUNCT__ "PetscDualSpaceInitialize_Lagrange" 1884 PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp) 1885 { 1886 PetscFunctionBegin; 1887 sp->ops->setfromoptions = NULL; 1888 sp->ops->setup = PetscDualSpaceSetUp_Lagrange; 1889 sp->ops->view = NULL; 1890 sp->ops->destroy = PetscDualSpaceDestroy_Lagrange; 1891 sp->ops->getdimension = PetscDualSpaceGetDimension_Lagrange; 1892 sp->ops->getnumdof = PetscDualSpaceGetNumDof_Lagrange; 1893 PetscFunctionReturn(0); 1894 } 1895 1896 /*MC 1897 PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals 1898 1899 Level: intermediate 1900 1901 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 1902 M*/ 1903 1904 #undef __FUNCT__ 1905 #define __FUNCT__ "PetscDualSpaceCreate_Lagrange" 1906 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp) 1907 { 1908 PetscDualSpace_Lag *lag; 1909 PetscErrorCode ierr; 1910 1911 PetscFunctionBegin; 1912 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 1913 ierr = PetscNewLog(sp, PetscDualSpace_Lag, &lag);CHKERRQ(ierr); 1914 sp->data = lag; 1915 1916 lag->numDof = NULL; 1917 1918 ierr = PetscDualSpaceInitialize_Lagrange(sp);CHKERRQ(ierr); 1919 PetscFunctionReturn(0); 1920 } 1921 1922 1923 PetscClassId PETSCFE_CLASSID = 0; 1924 1925 PetscFunctionList PetscFEList = NULL; 1926 PetscBool PetscFERegisterAllCalled = PETSC_FALSE; 1927 1928 #undef __FUNCT__ 1929 #define __FUNCT__ "PetscFERegister" 1930 /*@C 1931 PetscFERegister - Adds a new PetscFE implementation 1932 1933 Not Collective 1934 1935 Input Parameters: 1936 + name - The name of a new user-defined creation routine 1937 - create_func - The creation routine itself 1938 1939 Notes: 1940 PetscFERegister() may be called multiple times to add several user-defined PetscFEs 1941 1942 Sample usage: 1943 .vb 1944 PetscFERegister("my_fe", MyPetscFECreate); 1945 .ve 1946 1947 Then, your PetscFE type can be chosen with the procedural interface via 1948 .vb 1949 PetscFECreate(MPI_Comm, PetscFE *); 1950 PetscFESetType(PetscFE, "my_fe"); 1951 .ve 1952 or at runtime via the option 1953 .vb 1954 -petscfe_type my_fe 1955 .ve 1956 1957 Level: advanced 1958 1959 .keywords: PetscFE, register 1960 .seealso: PetscFERegisterAll(), PetscFERegisterDestroy() 1961 1962 @*/ 1963 PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE)) 1964 { 1965 PetscErrorCode ierr; 1966 1967 PetscFunctionBegin; 1968 ierr = PetscFunctionListAdd(&PetscFEList, sname, function);CHKERRQ(ierr); 1969 PetscFunctionReturn(0); 1970 } 1971 1972 #undef __FUNCT__ 1973 #define __FUNCT__ "PetscFESetType" 1974 /*@C 1975 PetscFESetType - Builds a particular PetscFE 1976 1977 Collective on PetscFE 1978 1979 Input Parameters: 1980 + fem - The PetscFE object 1981 - name - The kind of FEM space 1982 1983 Options Database Key: 1984 . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types 1985 1986 Level: intermediate 1987 1988 .keywords: PetscFE, set, type 1989 .seealso: PetscFEGetType(), PetscFECreate() 1990 @*/ 1991 PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name) 1992 { 1993 PetscErrorCode (*r)(PetscFE); 1994 PetscBool match; 1995 PetscErrorCode ierr; 1996 1997 PetscFunctionBegin; 1998 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1999 ierr = PetscObjectTypeCompare((PetscObject) fem, name, &match);CHKERRQ(ierr); 2000 if (match) PetscFunctionReturn(0); 2001 2002 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} 2003 ierr = PetscFunctionListFind(PetscFEList, name, &r);CHKERRQ(ierr); 2004 if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name); 2005 2006 if (fem->ops->destroy) { 2007 ierr = (*fem->ops->destroy)(fem);CHKERRQ(ierr); 2008 fem->ops->destroy = NULL; 2009 } 2010 ierr = (*r)(fem);CHKERRQ(ierr); 2011 ierr = PetscObjectChangeTypeName((PetscObject) fem, name);CHKERRQ(ierr); 2012 PetscFunctionReturn(0); 2013 } 2014 2015 #undef __FUNCT__ 2016 #define __FUNCT__ "PetscFEGetType" 2017 /*@C 2018 PetscFEGetType - Gets the PetscFE type name (as a string) from the object. 2019 2020 Not Collective 2021 2022 Input Parameter: 2023 . dm - The PetscFE 2024 2025 Output Parameter: 2026 . name - The PetscFE type name 2027 2028 Level: intermediate 2029 2030 .keywords: PetscFE, get, type, name 2031 .seealso: PetscFESetType(), PetscFECreate() 2032 @*/ 2033 PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name) 2034 { 2035 PetscErrorCode ierr; 2036 2037 PetscFunctionBegin; 2038 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2039 PetscValidCharPointer(name, 2); 2040 if (!PetscFERegisterAllCalled) { 2041 ierr = PetscFERegisterAll();CHKERRQ(ierr); 2042 } 2043 *name = ((PetscObject) fem)->type_name; 2044 PetscFunctionReturn(0); 2045 } 2046 2047 #undef __FUNCT__ 2048 #define __FUNCT__ "PetscFEView" 2049 /*@C 2050 PetscFEView - Views a PetscFE 2051 2052 Collective on PetscFE 2053 2054 Input Parameter: 2055 + fem - the PetscFE object to view 2056 - v - the viewer 2057 2058 Level: developer 2059 2060 .seealso PetscFEDestroy() 2061 @*/ 2062 PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v) 2063 { 2064 PetscErrorCode ierr; 2065 2066 PetscFunctionBegin; 2067 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2068 if (!v) { 2069 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);CHKERRQ(ierr); 2070 } 2071 if (fem->ops->view) { 2072 ierr = (*fem->ops->view)(fem, v);CHKERRQ(ierr); 2073 } 2074 PetscFunctionReturn(0); 2075 } 2076 2077 #undef __FUNCT__ 2078 #define __FUNCT__ "PetscFEViewFromOptions" 2079 /* 2080 PetscFEViewFromOptions - Processes command line options to determine if/how a PetscFE is to be viewed. 2081 2082 Collective on PetscFE 2083 2084 Input Parameters: 2085 + fem - the PetscFE 2086 . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd' 2087 - optionname - option to activate viewing 2088 2089 Level: intermediate 2090 2091 .keywords: PetscFE, view, options, database 2092 .seealso: VecViewFromOptions(), MatViewFromOptions() 2093 */ 2094 PetscErrorCode PetscFEViewFromOptions(PetscFE fem, const char prefix[], const char optionname[]) 2095 { 2096 PetscViewer viewer; 2097 PetscViewerFormat format; 2098 PetscBool flg; 2099 PetscErrorCode ierr; 2100 2101 PetscFunctionBegin; 2102 if (prefix) { 2103 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), prefix, optionname, &viewer, &format, &flg);CHKERRQ(ierr); 2104 } else { 2105 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), ((PetscObject) fem)->prefix, optionname, &viewer, &format, &flg);CHKERRQ(ierr); 2106 } 2107 if (flg) { 2108 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 2109 ierr = PetscFEView(fem, viewer);CHKERRQ(ierr); 2110 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 2111 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2112 } 2113 PetscFunctionReturn(0); 2114 } 2115 2116 #undef __FUNCT__ 2117 #define __FUNCT__ "PetscFESetFromOptions" 2118 /*@ 2119 PetscFESetFromOptions - sets parameters in a PetscFE from the options database 2120 2121 Collective on PetscFE 2122 2123 Input Parameter: 2124 . fem - the PetscFE object to set options for 2125 2126 Level: developer 2127 2128 .seealso PetscFEView() 2129 @*/ 2130 PetscErrorCode PetscFESetFromOptions(PetscFE fem) 2131 { 2132 const char *defaultType; 2133 char name[256]; 2134 PetscBool flg; 2135 PetscErrorCode ierr; 2136 2137 PetscFunctionBegin; 2138 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2139 if (!((PetscObject) fem)->type_name) { 2140 defaultType = PETSCFEBASIC; 2141 } else { 2142 defaultType = ((PetscObject) fem)->type_name; 2143 } 2144 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} 2145 2146 ierr = PetscObjectOptionsBegin((PetscObject) fem);CHKERRQ(ierr); 2147 ierr = PetscOptionsList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);CHKERRQ(ierr); 2148 if (flg) { 2149 ierr = PetscFESetType(fem, name);CHKERRQ(ierr); 2150 } else if (!((PetscObject) fem)->type_name) { 2151 ierr = PetscFESetType(fem, defaultType);CHKERRQ(ierr); 2152 } 2153 if (fem->ops->setfromoptions) { 2154 ierr = (*fem->ops->setfromoptions)(fem);CHKERRQ(ierr); 2155 } 2156 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 2157 ierr = PetscObjectProcessOptionsHandlers((PetscObject) fem);CHKERRQ(ierr); 2158 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2159 ierr = PetscFEViewFromOptions(fem, NULL, "-petscfe_view");CHKERRQ(ierr); 2160 PetscFunctionReturn(0); 2161 } 2162 2163 #undef __FUNCT__ 2164 #define __FUNCT__ "PetscFESetUp" 2165 /*@C 2166 PetscFESetUp - Construct data structures for the PetscFE 2167 2168 Collective on PetscFE 2169 2170 Input Parameter: 2171 . fem - the PetscFE object to setup 2172 2173 Level: developer 2174 2175 .seealso PetscFEView(), PetscFEDestroy() 2176 @*/ 2177 PetscErrorCode PetscFESetUp(PetscFE fem) 2178 { 2179 PetscErrorCode ierr; 2180 2181 PetscFunctionBegin; 2182 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2183 if (fem->ops->setup) {ierr = (*fem->ops->setup)(fem);CHKERRQ(ierr);} 2184 PetscFunctionReturn(0); 2185 } 2186 2187 #undef __FUNCT__ 2188 #define __FUNCT__ "PetscFEDestroy" 2189 /*@ 2190 PetscFEDestroy - Destroys a PetscFE object 2191 2192 Collective on PetscFE 2193 2194 Input Parameter: 2195 . fem - the PetscFE object to destroy 2196 2197 Level: developer 2198 2199 .seealso PetscFEView() 2200 @*/ 2201 PetscErrorCode PetscFEDestroy(PetscFE *fem) 2202 { 2203 PetscErrorCode ierr; 2204 2205 PetscFunctionBegin; 2206 if (!*fem) PetscFunctionReturn(0); 2207 PetscValidHeaderSpecific((*fem), PETSCFE_CLASSID, 1); 2208 2209 if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; PetscFunctionReturn(0);} 2210 ((PetscObject) (*fem))->refct = 0; 2211 /* if memory was published with AMS then destroy it */ 2212 ierr = PetscObjectAMSViewOff((PetscObject) *fem);CHKERRQ(ierr); 2213 2214 ierr = PetscFree((*fem)->numDof);CHKERRQ(ierr); 2215 ierr = PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);CHKERRQ(ierr); 2216 ierr = PetscSpaceDestroy(&(*fem)->basisSpace);CHKERRQ(ierr); 2217 ierr = PetscDualSpaceDestroy(&(*fem)->dualSpace);CHKERRQ(ierr); 2218 ierr = PetscQuadratureDestroy(&(*fem)->quadrature);CHKERRQ(ierr); 2219 2220 if ((*fem)->ops->destroy) {ierr = (*(*fem)->ops->destroy)(*fem);CHKERRQ(ierr);} 2221 ierr = PetscHeaderDestroy(fem);CHKERRQ(ierr); 2222 PetscFunctionReturn(0); 2223 } 2224 2225 #undef __FUNCT__ 2226 #define __FUNCT__ "PetscFECreate" 2227 /*@ 2228 PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType(). 2229 2230 Collective on MPI_Comm 2231 2232 Input Parameter: 2233 . comm - The communicator for the PetscFE object 2234 2235 Output Parameter: 2236 . fem - The PetscFE object 2237 2238 Level: beginner 2239 2240 .seealso: PetscFESetType(), PETSCFEGALERKIN 2241 @*/ 2242 PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem) 2243 { 2244 PetscFE f; 2245 PetscErrorCode ierr; 2246 2247 PetscFunctionBegin; 2248 PetscValidPointer(fem, 2); 2249 *fem = NULL; 2250 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 2251 ierr = PetscFEInitializePackage();CHKERRQ(ierr); 2252 #endif 2253 2254 ierr = PetscHeaderCreate(f, _p_PetscFE, struct _PetscFEOps, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);CHKERRQ(ierr); 2255 ierr = PetscMemzero(f->ops, sizeof(struct _PetscFEOps));CHKERRQ(ierr); 2256 2257 f->basisSpace = NULL; 2258 f->dualSpace = NULL; 2259 f->numComponents = 1; 2260 f->numDof = NULL; 2261 f->B = NULL; 2262 f->D = NULL; 2263 f->H = NULL; 2264 ierr = PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));CHKERRQ(ierr); 2265 2266 *fem = f; 2267 PetscFunctionReturn(0); 2268 } 2269 2270 #undef __FUNCT__ 2271 #define __FUNCT__ "PetscFEGetDimension" 2272 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 2273 { 2274 PetscErrorCode ierr; 2275 2276 PetscFunctionBegin; 2277 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2278 PetscValidPointer(dim, 2); 2279 ierr = PetscSpaceGetDimension(fem->basisSpace, dim);CHKERRQ(ierr); 2280 PetscFunctionReturn(0); 2281 } 2282 2283 #undef __FUNCT__ 2284 #define __FUNCT__ "PetscFEGetSpatialDimension" 2285 PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim) 2286 { 2287 DM dm; 2288 PetscErrorCode ierr; 2289 2290 PetscFunctionBegin; 2291 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2292 PetscValidPointer(dim, 2); 2293 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2294 ierr = DMPlexGetDimension(dm, dim);CHKERRQ(ierr); 2295 PetscFunctionReturn(0); 2296 } 2297 2298 #undef __FUNCT__ 2299 #define __FUNCT__ "PetscFESetNumComponents" 2300 PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp) 2301 { 2302 PetscFunctionBegin; 2303 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2304 fem->numComponents = comp; 2305 PetscFunctionReturn(0); 2306 } 2307 2308 #undef __FUNCT__ 2309 #define __FUNCT__ "PetscFEGetNumComponents" 2310 PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp) 2311 { 2312 PetscFunctionBegin; 2313 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2314 PetscValidPointer(comp, 2); 2315 *comp = fem->numComponents; 2316 PetscFunctionReturn(0); 2317 } 2318 2319 #undef __FUNCT__ 2320 #define __FUNCT__ "PetscFEGetBasisSpace" 2321 PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp) 2322 { 2323 PetscFunctionBegin; 2324 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2325 PetscValidPointer(sp, 2); 2326 *sp = fem->basisSpace; 2327 PetscFunctionReturn(0); 2328 } 2329 2330 #undef __FUNCT__ 2331 #define __FUNCT__ "PetscFESetBasisSpace" 2332 PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp) 2333 { 2334 PetscErrorCode ierr; 2335 2336 PetscFunctionBegin; 2337 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2338 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 2); 2339 ierr = PetscSpaceDestroy(&fem->basisSpace);CHKERRQ(ierr); 2340 fem->basisSpace = sp; 2341 ierr = PetscObjectReference((PetscObject) fem->basisSpace);CHKERRQ(ierr); 2342 PetscFunctionReturn(0); 2343 } 2344 2345 #undef __FUNCT__ 2346 #define __FUNCT__ "PetscFEGetDualSpace" 2347 PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp) 2348 { 2349 PetscFunctionBegin; 2350 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2351 PetscValidPointer(sp, 2); 2352 *sp = fem->dualSpace; 2353 PetscFunctionReturn(0); 2354 } 2355 2356 #undef __FUNCT__ 2357 #define __FUNCT__ "PetscFESetDualSpace" 2358 PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp) 2359 { 2360 PetscErrorCode ierr; 2361 2362 PetscFunctionBegin; 2363 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2364 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 2); 2365 ierr = PetscDualSpaceDestroy(&fem->dualSpace);CHKERRQ(ierr); 2366 fem->dualSpace = sp; 2367 ierr = PetscObjectReference((PetscObject) fem->dualSpace);CHKERRQ(ierr); 2368 PetscFunctionReturn(0); 2369 } 2370 2371 #undef __FUNCT__ 2372 #define __FUNCT__ "PetscFEGetQuadrature" 2373 PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q) 2374 { 2375 PetscFunctionBegin; 2376 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2377 PetscValidPointer(q, 2); 2378 *q = fem->quadrature; 2379 PetscFunctionReturn(0); 2380 } 2381 2382 #undef __FUNCT__ 2383 #define __FUNCT__ "PetscFESetQuadrature" 2384 PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q) 2385 { 2386 PetscErrorCode ierr; 2387 2388 PetscFunctionBegin; 2389 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2390 ierr = PetscQuadratureDestroy(&fem->quadrature);CHKERRQ(ierr); 2391 fem->quadrature = q; 2392 PetscFunctionReturn(0); 2393 } 2394 2395 #undef __FUNCT__ 2396 #define __FUNCT__ "PetscFEGetNumDof" 2397 PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof) 2398 { 2399 const PetscInt *numDofDual; 2400 PetscErrorCode ierr; 2401 2402 PetscFunctionBegin; 2403 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2404 PetscValidPointer(numDof, 2); 2405 ierr = PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);CHKERRQ(ierr); 2406 if (!fem->numDof) { 2407 DM dm; 2408 PetscInt dim, d; 2409 2410 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2411 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2412 ierr = PetscMalloc((dim+1) * sizeof(PetscInt), &fem->numDof);CHKERRQ(ierr); 2413 for (d = 0; d <= dim; ++d) { 2414 fem->numDof[d] = fem->numComponents*numDofDual[d]; 2415 } 2416 } 2417 *numDof = fem->numDof; 2418 PetscFunctionReturn(0); 2419 } 2420 2421 #undef __FUNCT__ 2422 #define __FUNCT__ "PetscFEGetDefaultTabulation" 2423 PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H) 2424 { 2425 PetscInt npoints = fem->quadrature.numQuadPoints; 2426 const PetscReal *points = fem->quadrature.quadPoints; 2427 PetscErrorCode ierr; 2428 2429 PetscFunctionBegin; 2430 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2431 if (B) PetscValidPointer(B, 2); 2432 if (D) PetscValidPointer(D, 3); 2433 if (H) PetscValidPointer(H, 4); 2434 if (!fem->B) {ierr = PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);CHKERRQ(ierr);} 2435 if (B) *B = fem->B; 2436 if (D) *D = fem->D; 2437 if (H) *H = fem->H; 2438 PetscFunctionReturn(0); 2439 } 2440 2441 #undef __FUNCT__ 2442 #define __FUNCT__ "PetscFEGetTabulation" 2443 PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H) 2444 { 2445 DM dm; 2446 PetscInt pdim; /* Dimension of FE space P */ 2447 PetscInt dim; /* Spatial dimension */ 2448 PetscInt comp; /* Field components */ 2449 PetscReal *tmpB, *tmpD, *invV; 2450 PetscInt p, d, j, k; 2451 PetscErrorCode ierr; 2452 2453 PetscFunctionBegin; 2454 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2455 PetscValidPointer(points, 3); 2456 if (B) PetscValidPointer(B, 4); 2457 if (D) PetscValidPointer(D, 5); 2458 if (H) PetscValidPointer(H, 6); 2459 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2460 2461 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2462 ierr = PetscSpaceGetDimension(fem->basisSpace, &pdim);CHKERRQ(ierr); 2463 ierr = PetscFEGetNumComponents(fem, &comp);CHKERRQ(ierr); 2464 /* if (nvalues%dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate values %d must be divisible by the spatial dimension %d", nvalues, dim); */ 2465 2466 if (B) { 2467 ierr = DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);CHKERRQ(ierr); 2468 ierr = DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);CHKERRQ(ierr); 2469 } 2470 if (D) { 2471 ierr = DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);CHKERRQ(ierr); 2472 ierr = DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);CHKERRQ(ierr); 2473 } 2474 if (H) {ierr = DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, H);CHKERRQ(ierr);} 2475 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? *H : NULL);CHKERRQ(ierr); 2476 2477 ierr = DMGetWorkArray(dm, pdim*pdim, PETSC_REAL, &invV);CHKERRQ(ierr); 2478 for (j = 0; j < pdim; ++j) { 2479 PetscReal *Bf; 2480 PetscQuadrature f; 2481 PetscInt q; 2482 2483 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f); 2484 ierr = DMGetWorkArray(dm, f.numQuadPoints*pdim, PETSC_REAL, &Bf);CHKERRQ(ierr); 2485 ierr = PetscSpaceEvaluate(fem->basisSpace, f.numQuadPoints, f.quadPoints, Bf, NULL, NULL);CHKERRQ(ierr); 2486 for (k = 0; k < pdim; ++k) { 2487 /* n_j \cdot \phi_k */ 2488 invV[j*pdim+k] = 0.0; 2489 for (q = 0; q < f.numQuadPoints; ++q) { 2490 invV[j*pdim+k] += Bf[q*pdim+k]*f.quadWeights[q]; 2491 } 2492 } 2493 ierr = DMRestoreWorkArray(dm, f.numQuadPoints*pdim, PETSC_REAL, &Bf);CHKERRQ(ierr); 2494 } 2495 { 2496 PetscReal *work; 2497 PetscBLASInt *pivots; 2498 PetscBLASInt n = pdim, info; 2499 2500 ierr = DMGetWorkArray(dm, pdim, PETSC_INT, &pivots);CHKERRQ(ierr); 2501 ierr = DMGetWorkArray(dm, pdim, PETSC_REAL, &work);CHKERRQ(ierr); 2502 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invV, &n, pivots, &info)); 2503 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invV, &n, pivots, work, &n, &info)); 2504 ierr = DMRestoreWorkArray(dm, pdim, PETSC_INT, &pivots);CHKERRQ(ierr); 2505 ierr = DMRestoreWorkArray(dm, pdim, PETSC_REAL, &work);CHKERRQ(ierr); 2506 } 2507 for (p = 0; p < npoints; ++p) { 2508 if (B) { 2509 /* Multiply by V^{-1} (pdim x pdim) */ 2510 for (j = 0; j < pdim; ++j) { 2511 const PetscInt i = (p*pdim + j)*comp; 2512 PetscInt c; 2513 2514 (*B)[i] = 0.0; 2515 for (k = 0; k < pdim; ++k) { 2516 (*B)[i] += invV[k*pdim+j] * tmpB[p*pdim + k]; 2517 } 2518 for (c = 1; c < comp; ++c) { 2519 (*B)[i+c] = (*B)[i]; 2520 } 2521 } 2522 } 2523 if (D) { 2524 /* Multiply by V^{-1} (pdim x pdim) */ 2525 for (j = 0; j < pdim; ++j) { 2526 for (d = 0; d < dim; ++d) { 2527 const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d; 2528 PetscInt c; 2529 2530 (*D)[i] = 0.0; 2531 for (k = 0; k < pdim; ++k) { 2532 (*D)[i] += invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d]; 2533 } 2534 for (c = 1; c < comp; ++c) { 2535 (*D)[((p*pdim + j)*comp + c)*dim + d] = (*D)[i]; 2536 } 2537 } 2538 } 2539 } 2540 } 2541 ierr = DMRestoreWorkArray(dm, pdim*pdim, PETSC_REAL, &invV);CHKERRQ(ierr); 2542 if (B) {ierr = DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);CHKERRQ(ierr);} 2543 if (D) {ierr = DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);CHKERRQ(ierr);} 2544 PetscFunctionReturn(0); 2545 } 2546 2547 #undef __FUNCT__ 2548 #define __FUNCT__ "PetscFERestoreTabulation" 2549 PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H) 2550 { 2551 DM dm; 2552 PetscErrorCode ierr; 2553 2554 PetscFunctionBegin; 2555 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2556 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2557 if (B) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, B);CHKERRQ(ierr);} 2558 if (D) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, D);CHKERRQ(ierr);} 2559 if (H) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, H);CHKERRQ(ierr);} 2560 PetscFunctionReturn(0); 2561 } 2562 2563 #undef __FUNCT__ 2564 #define __FUNCT__ "PetscFEDestroy_Basic" 2565 PetscErrorCode PetscFEDestroy_Basic(PetscFE fem) 2566 { 2567 PetscFE_Basic *b = (PetscFE_Basic *) fem->data; 2568 PetscErrorCode ierr; 2569 2570 PetscFunctionBegin; 2571 ierr = PetscFree(b);CHKERRQ(ierr); 2572 PetscFunctionReturn(0); 2573 } 2574 2575 #undef __FUNCT__ 2576 #define __FUNCT__ "PetscFEIntegrateResidual_Basic" 2577 PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt field, PetscCellGeometry geom, const PetscScalar coefficients[], 2578 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[]), 2579 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f1[]), 2580 PetscScalar elemVec[]) 2581 { 2582 const PetscInt debug = 0; 2583 PetscQuadrature quad; 2584 PetscScalar *f0, *f1, *u, *gradU; 2585 PetscReal *x, *realSpaceDer; 2586 PetscInt dim, numComponents = 0, cOffset = 0, eOffset = 0, e, f; 2587 PetscErrorCode ierr; 2588 2589 PetscFunctionBegin; 2590 ierr = PetscFEGetSpatialDimension(fe[0], &dim);CHKERRQ(ierr); 2591 for (f = 0; f < Nf; ++f) { 2592 PetscInt Nc; 2593 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 2594 numComponents += Nc; 2595 } 2596 ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr); 2597 ierr = PetscMalloc6(quad.numQuadPoints*dim,PetscScalar,&f0,quad.numQuadPoints*dim*dim,PetscScalar,&f1,numComponents,PetscScalar,&u,numComponents*dim,PetscScalar,&gradU,dim,PetscReal,&x,dim,PetscReal,&realSpaceDer); 2598 for (e = 0; e < Ne; ++e) { 2599 const PetscReal detJ = geom.detJ[e]; 2600 const PetscReal *v0 = &geom.v0[e*dim]; 2601 const PetscReal *J = &geom.J[e*dim*dim]; 2602 const PetscReal *invJ = &geom.invJ[e*dim*dim]; 2603 const PetscInt Nq = quad.numQuadPoints; 2604 const PetscReal *quadPoints = quad.quadPoints; 2605 const PetscReal *quadWeights = quad.quadWeights; 2606 PetscInt q, f; 2607 2608 if (debug > 1) { 2609 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 2610 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 2611 } 2612 for (q = 0; q < Nq; ++q) { 2613 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 2614 PetscInt fOffset = 0; 2615 PetscInt dOffset = cOffset; 2616 PetscInt Ncomp, d, d2, f, i; 2617 2618 ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr); 2619 for (d = 0; d < numComponents; ++d) {u[d] = 0.0;} 2620 for (d = 0; d < dim*(numComponents); ++d) {gradU[d] = 0.0;} 2621 for (d = 0; d < dim; ++d) { 2622 x[d] = v0[d]; 2623 for (d2 = 0; d2 < dim; ++d2) { 2624 x[d] += J[d*dim+d2]*(quadPoints[q*dim+d2] + 1.0); 2625 } 2626 } 2627 for (f = 0; f < Nf; ++f) { 2628 PetscReal *basis, *basisDer; 2629 PetscInt Nb, Ncomp, b, comp; 2630 2631 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2632 ierr = PetscFEGetNumComponents(fe[f], &Ncomp);CHKERRQ(ierr); 2633 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2634 ierr = PetscFEGetDefaultTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2635 for (b = 0; b < Nb; ++b) { 2636 for (comp = 0; comp < Ncomp; ++comp) { 2637 const PetscInt cidx = b*Ncomp+comp; 2638 PetscInt d, g; 2639 2640 u[fOffset+comp] += coefficients[dOffset+cidx]*basis[q*Nb*Ncomp+cidx]; 2641 for (d = 0; d < dim; ++d) { 2642 realSpaceDer[d] = 0.0; 2643 for (g = 0; g < dim; ++g) { 2644 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g]; 2645 } 2646 gradU[(fOffset+comp)*dim+d] += coefficients[dOffset+cidx]*realSpaceDer[d]; 2647 } 2648 } 2649 } 2650 if (debug > 1) { 2651 PetscInt d; 2652 for (comp = 0; comp < Ncomp; ++comp) { 2653 ierr = PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, comp, u[fOffset+comp]);CHKERRQ(ierr); 2654 for (d = 0; d < dim; ++d) { 2655 ierr = PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, comp, 'x'+d, gradU[(fOffset+comp)*dim+d]);CHKERRQ(ierr); 2656 } 2657 } 2658 } 2659 fOffset += Ncomp; 2660 dOffset += Nb*Ncomp; 2661 } 2662 2663 f0_func(u, gradU, NULL, NULL, x, &f0[q*Ncomp]); 2664 for (i = 0; i < Ncomp; ++i) { 2665 f0[q*Ncomp+i] *= detJ*quadWeights[q]; 2666 } 2667 f1_func(u, gradU, NULL, NULL, x, &f1[q*Ncomp*dim]); 2668 for (i = 0; i < Ncomp*dim; ++i) { 2669 f1[q*Ncomp*dim+i] *= detJ*quadWeights[q]; 2670 } 2671 if (debug > 1) { 2672 PetscInt c,d; 2673 for (c = 0; c < Ncomp; ++c) { 2674 ierr = PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, f0[q*Ncomp+c]);CHKERRQ(ierr); 2675 for (d = 0; d < dim; ++d) { 2676 ierr = PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, f1[(q*Ncomp + c)*dim+d]);CHKERRQ(ierr); 2677 } 2678 } 2679 } 2680 if (q == Nq-1) {cOffset = dOffset;} 2681 } 2682 for (f = 0; f < Nf; ++f) { 2683 PetscInt Nb, Ncomp, b, comp; 2684 2685 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2686 ierr = PetscFEGetNumComponents(fe[f], &Ncomp);CHKERRQ(ierr); 2687 if (f == field) { 2688 PetscReal *basis; 2689 PetscReal *basisDer; 2690 2691 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2692 ierr = PetscFEGetDefaultTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2693 for (b = 0; b < Nb; ++b) { 2694 for (comp = 0; comp < Ncomp; ++comp) { 2695 const PetscInt cidx = b*Ncomp+comp; 2696 PetscInt q; 2697 2698 elemVec[eOffset+cidx] = 0.0; 2699 for (q = 0; q < Nq; ++q) { 2700 PetscInt d, g; 2701 2702 elemVec[eOffset+cidx] += basis[q*Nb*Ncomp+cidx]*f0[q*Ncomp+comp]; 2703 for (d = 0; d < dim; ++d) { 2704 realSpaceDer[d] = 0.0; 2705 for (g = 0; g < dim; ++g) { 2706 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g]; 2707 } 2708 elemVec[eOffset+cidx] += realSpaceDer[d]*f1[(q*Ncomp+comp)*dim+d]; 2709 } 2710 } 2711 } 2712 } 2713 if (debug > 1) { 2714 PetscInt b, comp; 2715 2716 for (b = 0; b < Nb; ++b) { 2717 for (comp = 0; comp < Ncomp; ++comp) { 2718 ierr = PetscPrintf(PETSC_COMM_SELF, " elemVec[%d,%d]: %g\n", b, comp, elemVec[eOffset+b*Ncomp+comp]);CHKERRQ(ierr); 2719 } 2720 } 2721 } 2722 } 2723 eOffset += Nb*Ncomp; 2724 } 2725 } 2726 ierr = PetscFree6(f0,f1,u,gradU,x,realSpaceDer); 2727 PetscFunctionReturn(0); 2728 } 2729 2730 #undef __FUNCT__ 2731 #define __FUNCT__ "PetscFEIntegrateJacobian_Basic" 2732 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt fieldI, PetscInt fieldJ, PetscCellGeometry geom, const PetscScalar coefficients[], 2733 void (*g0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g0[]), 2734 void (*g1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g1[]), 2735 void (*g2_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g2[]), 2736 void (*g3_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]), 2737 PetscScalar elemMat[]) 2738 { 2739 const PetscInt debug = 0; 2740 PetscInt cellDof = 0; /* Total number of dof on a cell */ 2741 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 2742 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 2743 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 2744 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 2745 PetscQuadrature quad; 2746 PetscScalar *g0, *g1, *g2, *g3, *u, *gradU; 2747 PetscReal *x, *realSpaceDerI, *realSpaceDerJ; 2748 PetscReal *basisI, *basisDerI, *basisJ, *basisDerJ; 2749 PetscInt NbI, NcI, NbJ, NcJ, numComponents = 0; 2750 PetscInt dim, f, e; 2751 PetscErrorCode ierr; 2752 2753 PetscFunctionBegin; 2754 ierr = PetscFEGetSpatialDimension(fe[fieldI], &dim);CHKERRQ(ierr); 2755 ierr = PetscFEGetDefaultTabulation(fe[fieldI], &basisI, &basisDerI, NULL);CHKERRQ(ierr); 2756 ierr = PetscFEGetDefaultTabulation(fe[fieldJ], &basisJ, &basisDerJ, NULL);CHKERRQ(ierr); 2757 for (f = 0; f < Nf; ++f) { 2758 PetscInt Nb, Nc; 2759 2760 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2761 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 2762 if (f == fieldI) {offsetI = cellDof; NbI = Nb; NcI = Nc;} 2763 if (f == fieldJ) {offsetJ = cellDof; NbJ = Nb; NcJ = Nc;} 2764 numComponents += Nc; 2765 cellDof += Nb*Nc; 2766 } 2767 ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr); 2768 ierr = PetscMalloc4(NcI*NcJ,PetscScalar,&g0,NcI*NcJ*dim,PetscScalar,&g1,NcI*NcJ*dim,PetscScalar,&g2,NcI*NcJ*dim*dim,PetscScalar,&g3); 2769 ierr = PetscMalloc5(numComponents,PetscScalar,&u,numComponents*dim,PetscScalar,&gradU,dim,PetscReal,&x,dim,PetscReal,&realSpaceDerI,dim,PetscReal,&realSpaceDerJ); 2770 for (e = 0; e < Ne; ++e) { 2771 const PetscReal detJ = geom.detJ[e]; 2772 const PetscReal *v0 = &geom.v0[e*dim]; 2773 const PetscReal *J = &geom.J[e*dim*dim]; 2774 const PetscReal *invJ = &geom.invJ[e*dim*dim]; 2775 const PetscInt Nq = quad.numQuadPoints; 2776 const PetscReal *quadPoints = quad.quadPoints; 2777 const PetscReal *quadWeights = quad.quadWeights; 2778 PetscInt f, g, q; 2779 2780 for (f = 0; f < NbI; ++f) { 2781 for (g = 0; g < NbJ; ++g) { 2782 for (q = 0; q < Nq; ++q) { 2783 PetscInt fOffset = 0; /* Offset into u[] for field_q (like offsetI) */ 2784 PetscInt dOffset = cOffset; /* Offset into coefficients[] for field_q */ 2785 PetscInt field_q, d, d2; 2786 PetscInt fc, gc, c; 2787 2788 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 2789 for (d = 0; d < numComponents; ++d) {u[d] = 0.0;} 2790 for (d = 0; d < dim*numComponents; ++d) {gradU[d] = 0.0;} 2791 for (d = 0; d < dim; ++d) { 2792 x[d] = v0[d]; 2793 for (d2 = 0; d2 < dim; ++d2) { 2794 x[d] += J[d*dim+d2]*(quadPoints[q*dim+d2] + 1.0); 2795 } 2796 } 2797 for (field_q = 0; field_q < Nf; ++field_q) { 2798 PetscReal *basis, *basisDer; 2799 PetscInt Nb, Ncomp, b, comp; 2800 2801 ierr = PetscFEGetDimension(fe[field_q], &Nb);CHKERRQ(ierr); 2802 ierr = PetscFEGetNumComponents(fe[field_q], &Ncomp);CHKERRQ(ierr); 2803 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2804 ierr = PetscFEGetDefaultTabulation(fe[field_q], &basis, &basisDer, NULL);CHKERRQ(ierr); 2805 for (b = 0; b < Nb; ++b) { 2806 for (comp = 0; comp < Ncomp; ++comp) { 2807 const PetscInt cidx = b*Ncomp+comp; 2808 PetscInt d1, d2; 2809 2810 u[fOffset+comp] += coefficients[dOffset+cidx]*basis[q*Nb*Ncomp+cidx]; 2811 for (d1 = 0; d1 < dim; ++d1) { 2812 realSpaceDerI[d1] = 0.0; 2813 for (d2 = 0; d2 < dim; ++d2) { 2814 realSpaceDerI[d1] += invJ[d2*dim+d1]*basisDer[(q*Nb*Ncomp+cidx)*dim+d2]; 2815 } 2816 gradU[(fOffset+comp)*dim+d1] += coefficients[dOffset+cidx]*realSpaceDerI[d1]; 2817 } 2818 } 2819 } 2820 if (debug > 1) { 2821 for (comp = 0; comp < Ncomp; ++comp) { 2822 ierr = PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, comp, u[fOffset+comp]);CHKERRQ(ierr); 2823 for (d = 0; d < dim; ++d) { 2824 ierr = PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, comp, 'x'+d, gradU[(fOffset+comp)*dim+d]);CHKERRQ(ierr); 2825 } 2826 } 2827 } 2828 fOffset += Ncomp; 2829 dOffset += Nb*Ncomp; 2830 } 2831 2832 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 2833 ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 2834 ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 2835 ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 2836 if (g0_func) { 2837 g0_func(u, gradU, NULL, NULL, x, g0); 2838 for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];} 2839 } 2840 if (g1_func) { 2841 g1_func(u, gradU, NULL, NULL, x, g1); 2842 for (c = 0; c < NcI*NcJ*dim; ++c) {g1[c] *= detJ*quadWeights[q];} 2843 } 2844 if (g2_func) { 2845 g2_func(u, gradU, NULL, NULL, x, g2); 2846 for (c = 0; c < NcI*NcJ*dim; ++c) {g2[c] *= detJ*quadWeights[q];} 2847 } 2848 if (g3_func) { 2849 g3_func(u, gradU, NULL, NULL, x, g3); 2850 for (c = 0; c < NcI*NcJ*dim*dim; ++c) {g3[c] *= detJ*quadWeights[q];} 2851 } 2852 2853 for (fc = 0; fc < NcI; ++fc) { 2854 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2855 const PetscInt i = offsetI+fidx; /* Element matrix row */ 2856 for (gc = 0; gc < NcJ; ++gc) { 2857 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2858 const PetscInt j = offsetJ+gidx; /* Element matrix column */ 2859 PetscInt d, d2; 2860 2861 for (d = 0; d < dim; ++d) { 2862 realSpaceDerI[d] = 0.0; 2863 realSpaceDerJ[d] = 0.0; 2864 for (d2 = 0; d2 < dim; ++d2) { 2865 realSpaceDerI[d] += invJ[d2*dim+d]*basisDerI[(q*NbI*NcI+fidx)*dim+d2]; 2866 realSpaceDerJ[d] += invJ[d2*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2]; 2867 } 2868 } 2869 elemMat[eOffset+i*cellDof+j] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx]; 2870 for (d = 0; d < dim; ++d) { 2871 elemMat[eOffset+i*cellDof+j] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*realSpaceDerJ[d]; 2872 elemMat[eOffset+i*cellDof+j] += realSpaceDerI[d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx]; 2873 for (d2 = 0; d2 < dim; ++d2) { 2874 elemMat[eOffset+i*cellDof+j] += realSpaceDerI[d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*realSpaceDerJ[d2]; 2875 } 2876 } 2877 } 2878 } 2879 } 2880 } 2881 } 2882 if (debug > 1) { 2883 PetscInt fc, f, gc, g; 2884 2885 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 2886 for (fc = 0; fc < NcI; ++fc) { 2887 for (f = 0; f < NbI; ++f) { 2888 const PetscInt i = offsetI + f*NcI+fc; 2889 for (gc = 0; gc < NcJ; ++gc) { 2890 for (g = 0; g < NbJ; ++g) { 2891 const PetscInt j = offsetJ + g*NcJ+gc; 2892 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, elemMat[eOffset+i*cellDof+j]);CHKERRQ(ierr); 2893 } 2894 } 2895 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 2896 } 2897 } 2898 } 2899 cOffset += cellDof; 2900 eOffset += cellDof*cellDof; 2901 } 2902 ierr = PetscFree4(g0,g1,g2,g3); 2903 ierr = PetscFree5(u,gradU,x,realSpaceDerI,realSpaceDerJ); 2904 PetscFunctionReturn(0); 2905 } 2906 2907 #undef __FUNCT__ 2908 #define __FUNCT__ "PetscFEInitialize_Basic" 2909 PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 2910 { 2911 PetscFunctionBegin; 2912 fem->ops->setfromoptions = NULL; 2913 fem->ops->setup = NULL; 2914 fem->ops->view = NULL; 2915 fem->ops->destroy = PetscFEDestroy_Basic; 2916 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 2917 fem->ops->integratebdresidual = NULL/* PetscFEIntegrateBdResidual_Basic */; 2918 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 2919 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 2920 PetscFunctionReturn(0); 2921 } 2922 2923 /*MC 2924 PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization 2925 2926 Level: intermediate 2927 2928 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 2929 M*/ 2930 2931 #undef __FUNCT__ 2932 #define __FUNCT__ "PetscFECreate_Basic" 2933 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 2934 { 2935 PetscFE_Basic *b; 2936 PetscErrorCode ierr; 2937 2938 PetscFunctionBegin; 2939 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2940 ierr = PetscNewLog(fem, PetscFE_Basic, &b);CHKERRQ(ierr); 2941 fem->data = b; 2942 2943 ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr); 2944 PetscFunctionReturn(0); 2945 } 2946 2947 /* 2948 Purpose: Compute element vector for chunk of elements 2949 2950 Input: 2951 Sizes: 2952 Ne: number of elements 2953 Nf: number of fields 2954 PetscFE 2955 dim: spatial dimension 2956 Nb: number of basis functions 2957 Nc: number of field components 2958 PetscQuadrature 2959 Nq: number of quadrature points 2960 2961 Geometry: 2962 PetscCellGeometry 2963 PetscReal v0s[Ne*dim] 2964 PetscReal jacobians[Ne*dim*dim] possibly *Nq 2965 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 2966 PetscReal jacobianDeterminants[Ne] possibly *Nq 2967 FEM: 2968 PetscFE 2969 PetscQuadrature 2970 PetscReal quadPoints[Nq*dim] 2971 PetscReal quadWeights[Nq] 2972 PetscReal basis[Nq*Nb*Nc] 2973 PetscReal basisDer[Nq*Nb*Nc*dim] 2974 PetscScalar coefficients[Ne*Nb*Nc] 2975 PetscScalar elemVec[Ne*Nb*Nc] 2976 2977 Problem: 2978 PetscInt f: the active field 2979 f0, f1 2980 2981 Work Space: 2982 PetscFE 2983 PetscScalar f0[Nq*dim]; 2984 PetscScalar f1[Nq*dim*dim]; 2985 PetscScalar u[Nc]; 2986 PetscScalar gradU[Nc*dim]; 2987 PetscReal x[dim]; 2988 PetscScalar realSpaceDer[dim]; 2989 2990 Purpose: Compute element vector for N_cb batches of elements 2991 2992 Input: 2993 Sizes: 2994 N_cb: Number of serial cell batches 2995 2996 Geometry: 2997 PetscReal v0s[Ne*dim] 2998 PetscReal jacobians[Ne*dim*dim] possibly *Nq 2999 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 3000 PetscReal jacobianDeterminants[Ne] possibly *Nq 3001 FEM: 3002 static PetscReal quadPoints[Nq*dim] 3003 static PetscReal quadWeights[Nq] 3004 static PetscReal basis[Nq*Nb*Nc] 3005 static PetscReal basisDer[Nq*Nb*Nc*dim] 3006 PetscScalar coefficients[Ne*Nb*Nc] 3007 PetscScalar elemVec[Ne*Nb*Nc] 3008 3009 ex62.c: 3010 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 3011 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 3012 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 3013 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 3014 3015 ex52.c: 3016 PetscErrorCode IntegrateLaplacianBatchCPU(PetscInt Ne, PetscInt Nb, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 3017 PetscErrorCode IntegrateElasticityBatchCPU(PetscInt Ne, PetscInt Nb, PetscInt Ncomp, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 3018 3019 ex52_integrateElement.cu 3020 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 3021 3022 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 3023 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 3024 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 3025 3026 ex52_integrateElementOpenCL.c: 3027 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 3028 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 3029 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 3030 3031 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 3032 */ 3033 3034 #undef __FUNCT__ 3035 #define __FUNCT__ "PetscFEIntegrateResidual" 3036 /*C 3037 PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration 3038 3039 Not collective 3040 3041 Input Parameters: 3042 + fem - The PetscFE object for the field being integrated 3043 . Ne - The number of elements in the chunk 3044 . Nf - The number of physical fields 3045 . fe - The PetscFE objects for each field 3046 . field - The field being integrated 3047 . geom - The cell geometry for each cell in the chunk 3048 . coefficients - The array of FEM basis coefficients for the elements 3049 . f0_func - f_0 function from the first order FEM model 3050 - f1_func - f_1 function from the first order FEM model 3051 3052 Output Parameter 3053 . elemVec - the element residual vectors from each element 3054 3055 Calling sequence of f0_func and f1_func: 3056 $ void f0(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[]) 3057 3058 Note: 3059 $ Loop over batch of elements (e): 3060 $ Loop over quadrature points (q): 3061 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q 3062 $ Call f_0 and f_1 3063 $ Loop over element vector entries (f,fc --> i): 3064 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 3065 */ 3066 PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt field, PetscCellGeometry geom, const PetscScalar coefficients[], 3067 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[]), 3068 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f1[]), 3069 PetscScalar elemVec[]) 3070 { 3071 PetscErrorCode ierr; 3072 3073 PetscFunctionBegin; 3074 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 3075 if (fem->ops->integrateresidual) {ierr = (*fem->ops->integrateresidual)(fem, Ne, Nf, fe, field, geom, coefficients, f0_func, f1_func, elemVec);CHKERRQ(ierr);} 3076 PetscFunctionReturn(0); 3077 } 3078 3079 #undef __FUNCT__ 3080 #define __FUNCT__ "PetscFEIntegrateJacobianAction" 3081 /*C 3082 PetscFEIntegrateJacobianAction - Produce the action of the element Jacobian on an element vector for a chunk of elements by quadrature integration 3083 3084 Not collective 3085 3086 Input Parameters: 3087 + fem = The PetscFE object for the field being integrated 3088 . Ne - The number of elements in the chunk 3089 . Nf - The number of physical fields 3090 . fe - The PetscFE objects for each field 3091 . field - The test field being integrated 3092 . geom - The cell geometry for each cell in the chunk 3093 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 3094 . input - The array of FEM basis coefficients for the elements for the input vector 3095 . g0_func - g_0 function from the first order FEM model 3096 . g1_func - g_1 function from the first order FEM model 3097 . g2_func - g_2 function from the first order FEM model 3098 - g3_func - g_3 function from the first order FEM model 3099 3100 Output Parameter 3101 . elemVec - the element vector for the action from each element 3102 3103 Calling sequence of g0_func, g1_func, g2_func and g3_func: 3104 $ void g0(PetscScalar u[], const PetscScalar gradU[], PetscScalar a[], const PetscScalar gradA[], PetscScalar x[], PetscScalar g0[]) 3105 3106 Note: 3107 $ Loop over batch of elements (e): 3108 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 3109 $ Loop over quadrature points (q): 3110 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 3111 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 3112 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 3113 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 3114 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 3115 */ 3116 PetscErrorCode PetscFEIntegrateJacobianAction(PetscFE fem, PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt field, PetscCellGeometry geom, const PetscScalar coefficients[], const PetscScalar input[], 3117 void (**g0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g0[]), 3118 void (**g1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g1[]), 3119 void (**g2_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g2[]), 3120 void (**g3_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]), 3121 PetscScalar elemVec[]) 3122 { 3123 PetscErrorCode ierr; 3124 3125 PetscFunctionBegin; 3126 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 3127 if (fem->ops->integratejacobianaction) {ierr = (*fem->ops->integratejacobianaction)(fem, Ne, Nf, fe, field, geom, coefficients, input, g0_func, g1_func, g2_func, g3_func, elemVec);CHKERRQ(ierr);} 3128 PetscFunctionReturn(0); 3129 } 3130 3131 #undef __FUNCT__ 3132 #define __FUNCT__ "PetscFEIntegrateJacobian" 3133 /*C 3134 PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration 3135 3136 Not collective 3137 3138 Input Parameters: 3139 + fem = The PetscFE object for the field being integrated 3140 . Ne - The number of elements in the chunk 3141 . Nf - The number of physical fields 3142 . fe - The PetscFE objects for each field 3143 . fieldI - The test field being integrated 3144 . fieldJ - The basis field being integrated 3145 . geom - The cell geometry for each cell in the chunk 3146 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 3147 . g0_func - g_0 function from the first order FEM model 3148 . g1_func - g_1 function from the first order FEM model 3149 . g2_func - g_2 function from the first order FEM model 3150 - g3_func - g_3 function from the first order FEM model 3151 3152 Output Parameter 3153 . elemMat - the element matrices for the Jacobian from each element 3154 3155 Calling sequence of g0_func, g1_func, g2_func and g3_func: 3156 $ void g0(PetscScalar u[], const PetscScalar gradU[], PetscScalar a[], const PetscScalar gradA[], PetscScalar x[], PetscScalar g0[]) 3157 3158 Note: 3159 $ Loop over batch of elements (e): 3160 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 3161 $ Loop over quadrature points (q): 3162 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 3163 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 3164 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 3165 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 3166 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 3167 */ 3168 PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt fieldI, PetscInt fieldJ, PetscCellGeometry geom, const PetscScalar coefficients[], 3169 void (*g0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g0[]), 3170 void (*g1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g1[]), 3171 void (*g2_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g2[]), 3172 void (*g3_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]), 3173 PetscScalar elemMat[]) 3174 { 3175 PetscErrorCode ierr; 3176 3177 PetscFunctionBegin; 3178 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 3179 if (fem->ops->integratejacobian) {ierr = (*fem->ops->integratejacobian)(fem, Ne, Nf, fe, fieldI, fieldJ, geom, coefficients, g0_func, g1_func, g2_func, g3_func, elemMat);CHKERRQ(ierr);} 3180 PetscFunctionReturn(0); 3181 } 3182