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 PetscInt 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 PetscInt 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 PetscInt PETSCFE_CLASSID = 0; 1924 1925 #undef __FUNCT__ 1926 #define __FUNCT__ "PetscFEView" 1927 /*@C 1928 PetscFEView - Views a PetscFE 1929 1930 Collective on PetscFE 1931 1932 Input Parameter: 1933 + sp - the PetscFE object to view 1934 - v - the viewer 1935 1936 Level: developer 1937 1938 .seealso PetscFEDestroy() 1939 @*/ 1940 PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v) 1941 { 1942 PetscErrorCode ierr; 1943 1944 PetscFunctionBegin; 1945 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1946 if (!v) { 1947 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);CHKERRQ(ierr); 1948 } 1949 if (fem->ops->view) { 1950 ierr = (*fem->ops->view)(fem, v);CHKERRQ(ierr); 1951 } 1952 PetscFunctionReturn(0); 1953 } 1954 1955 #undef __FUNCT__ 1956 #define __FUNCT__ "PetscFEDestroy" 1957 /*@ 1958 PetscFEDestroy - Destroys a PetscFE object 1959 1960 Collective on PetscFE 1961 1962 Input Parameter: 1963 . fem - the PetscFE object to destroy 1964 1965 Level: developer 1966 1967 .seealso PetscFEView() 1968 @*/ 1969 PetscErrorCode PetscFEDestroy(PetscFE *fem) 1970 { 1971 PetscErrorCode ierr; 1972 1973 PetscFunctionBegin; 1974 if (!*fem) PetscFunctionReturn(0); 1975 PetscValidHeaderSpecific((*fem), PETSCFE_CLASSID, 1); 1976 1977 if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; PetscFunctionReturn(0);} 1978 ((PetscObject) (*fem))->refct = 0; 1979 /* if memory was published with AMS then destroy it */ 1980 ierr = PetscObjectAMSViewOff((PetscObject) *fem);CHKERRQ(ierr); 1981 1982 ierr = PetscSpaceDestroy(&(*fem)->basisSpace);CHKERRQ(ierr); 1983 ierr = PetscDualSpaceDestroy(&(*fem)->dualSpace);CHKERRQ(ierr); 1984 ierr = PetscQuadratureDestroy(&(*fem)->quadrature);CHKERRQ(ierr); 1985 ierr = PetscFree((*fem)->numDof);CHKERRQ(ierr); 1986 1987 if ((*fem)->ops->destroy) {ierr = (*(*fem)->ops->destroy)(*fem);CHKERRQ(ierr);} 1988 ierr = PetscHeaderDestroy(fem);CHKERRQ(ierr); 1989 PetscFunctionReturn(0); 1990 } 1991 1992 #undef __FUNCT__ 1993 #define __FUNCT__ "PetscFECreate" 1994 /*@ 1995 PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType(). 1996 1997 Collective on MPI_Comm 1998 1999 Input Parameter: 2000 . comm - The communicator for the PetscFE object 2001 2002 Output Parameter: 2003 . fem - The PetscFE object 2004 2005 Level: beginner 2006 2007 .seealso: PetscFESetType(), PETSCFEGALERKIN 2008 @*/ 2009 PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem) 2010 { 2011 PetscFE f; 2012 PetscErrorCode ierr; 2013 2014 PetscFunctionBegin; 2015 PetscValidPointer(fem, 2); 2016 *fem = NULL; 2017 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 2018 ierr = PetscFEInitializePackage();CHKERRQ(ierr); 2019 #endif 2020 2021 ierr = PetscHeaderCreate(f, _p_PetscFE, struct _PetscFEOps, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);CHKERRQ(ierr); 2022 ierr = PetscMemzero(f->ops, sizeof(struct _PetscFEOps));CHKERRQ(ierr); 2023 2024 f->basisSpace = NULL; 2025 f->dualSpace = NULL; 2026 f->numComponents = 1; 2027 f->numDof = NULL; 2028 ierr = PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));CHKERRQ(ierr); 2029 2030 *fem = f; 2031 PetscFunctionReturn(0); 2032 } 2033 2034 #undef __FUNCT__ 2035 #define __FUNCT__ "PetscFEGetDimension" 2036 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 2037 { 2038 PetscErrorCode ierr; 2039 2040 PetscFunctionBegin; 2041 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2042 PetscValidPointer(dim, 2); 2043 ierr = PetscSpaceGetDimension(fem->basisSpace, dim);CHKERRQ(ierr); 2044 PetscFunctionReturn(0); 2045 } 2046 2047 #undef __FUNCT__ 2048 #define __FUNCT__ "PetscFEGetSpatialDimension" 2049 PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim) 2050 { 2051 DM dm; 2052 PetscErrorCode ierr; 2053 2054 PetscFunctionBegin; 2055 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2056 PetscValidPointer(dim, 2); 2057 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2058 ierr = DMPlexGetDimension(dm, dim);CHKERRQ(ierr); 2059 PetscFunctionReturn(0); 2060 } 2061 2062 #undef __FUNCT__ 2063 #define __FUNCT__ "PetscFESetNumComponents" 2064 PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp) 2065 { 2066 PetscFunctionBegin; 2067 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2068 fem->numComponents = comp; 2069 PetscFunctionReturn(0); 2070 } 2071 2072 #undef __FUNCT__ 2073 #define __FUNCT__ "PetscFEGetNumComponents" 2074 PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp) 2075 { 2076 PetscFunctionBegin; 2077 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2078 PetscValidPointer(comp, 2); 2079 *comp = fem->numComponents; 2080 PetscFunctionReturn(0); 2081 } 2082 2083 #undef __FUNCT__ 2084 #define __FUNCT__ "PetscFEGetBasisSpace" 2085 PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp) 2086 { 2087 PetscFunctionBegin; 2088 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2089 PetscValidPointer(sp, 2); 2090 *sp = fem->basisSpace; 2091 PetscFunctionReturn(0); 2092 } 2093 2094 #undef __FUNCT__ 2095 #define __FUNCT__ "PetscFESetBasisSpace" 2096 PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp) 2097 { 2098 PetscErrorCode ierr; 2099 2100 PetscFunctionBegin; 2101 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2102 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 2); 2103 ierr = PetscSpaceDestroy(&fem->basisSpace);CHKERRQ(ierr); 2104 fem->basisSpace = sp; 2105 ierr = PetscObjectReference((PetscObject) fem->basisSpace);CHKERRQ(ierr); 2106 PetscFunctionReturn(0); 2107 } 2108 2109 #undef __FUNCT__ 2110 #define __FUNCT__ "PetscFEGetDualSpace" 2111 PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp) 2112 { 2113 PetscFunctionBegin; 2114 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2115 PetscValidPointer(sp, 2); 2116 *sp = fem->dualSpace; 2117 PetscFunctionReturn(0); 2118 } 2119 2120 #undef __FUNCT__ 2121 #define __FUNCT__ "PetscFESetDualSpace" 2122 PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp) 2123 { 2124 PetscErrorCode ierr; 2125 2126 PetscFunctionBegin; 2127 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2128 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 2); 2129 ierr = PetscDualSpaceDestroy(&fem->dualSpace);CHKERRQ(ierr); 2130 fem->dualSpace = sp; 2131 ierr = PetscObjectReference((PetscObject) fem->dualSpace);CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 #undef __FUNCT__ 2136 #define __FUNCT__ "PetscFEGetQuadrature" 2137 PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q) 2138 { 2139 PetscFunctionBegin; 2140 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2141 PetscValidPointer(q, 2); 2142 *q = fem->quadrature; 2143 PetscFunctionReturn(0); 2144 } 2145 2146 #undef __FUNCT__ 2147 #define __FUNCT__ "PetscFESetQuadrature" 2148 PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q) 2149 { 2150 PetscErrorCode ierr; 2151 2152 PetscFunctionBegin; 2153 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2154 ierr = PetscQuadratureDestroy(&fem->quadrature);CHKERRQ(ierr); 2155 fem->quadrature = q; 2156 PetscFunctionReturn(0); 2157 } 2158 2159 #undef __FUNCT__ 2160 #define __FUNCT__ "PetscFEGetNumDof" 2161 PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof) 2162 { 2163 const PetscInt *numDofDual; 2164 PetscErrorCode ierr; 2165 2166 PetscFunctionBegin; 2167 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2168 PetscValidPointer(numDof, 2); 2169 ierr = PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);CHKERRQ(ierr); 2170 if (!fem->numDof) { 2171 DM dm; 2172 PetscInt dim, d; 2173 2174 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2175 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2176 ierr = PetscMalloc((dim+1) * sizeof(PetscInt), &fem->numDof);CHKERRQ(ierr); 2177 for (d = 0; d <= dim; ++d) { 2178 fem->numDof[d] = fem->numComponents*numDofDual[d]; 2179 } 2180 } 2181 *numDof = fem->numDof; 2182 PetscFunctionReturn(0); 2183 } 2184 2185 #undef __FUNCT__ 2186 #define __FUNCT__ "PetscFEGetTabulation" 2187 PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H) 2188 { 2189 DM dm; 2190 PetscInt pdim; /* Dimension of FE space P */ 2191 PetscInt dim; /* Spatial dimension */ 2192 PetscInt comp; /* Field components */ 2193 PetscInt npoints = fem->quadrature.numQuadPoints; 2194 const PetscReal *points = fem->quadrature.quadPoints; 2195 PetscReal *tmpB, *tmpD, *invV; 2196 PetscInt p, d, j, k; 2197 PetscErrorCode ierr; 2198 2199 PetscFunctionBegin; 2200 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2201 PetscValidPointer(points, 3); 2202 if (B) PetscValidPointer(B, 4); 2203 if (D) PetscValidPointer(D, 5); 2204 if (H) PetscValidPointer(H, 6); 2205 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2206 2207 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2208 ierr = PetscSpaceGetDimension(fem->basisSpace, &pdim);CHKERRQ(ierr); 2209 ierr = PetscFEGetNumComponents(fem, &comp);CHKERRQ(ierr); 2210 /* 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); */ 2211 2212 if (B) { 2213 ierr = DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);CHKERRQ(ierr); 2214 ierr = DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);CHKERRQ(ierr); 2215 } 2216 if (D) { 2217 ierr = DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);CHKERRQ(ierr); 2218 ierr = DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);CHKERRQ(ierr); 2219 } 2220 if (H) {ierr = DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, H);CHKERRQ(ierr);} 2221 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? *H : NULL);CHKERRQ(ierr); 2222 2223 ierr = DMGetWorkArray(dm, pdim*pdim, PETSC_REAL, &invV);CHKERRQ(ierr); 2224 for (j = 0; j < pdim; ++j) { 2225 PetscReal *Bf; 2226 PetscQuadrature f; 2227 PetscInt q; 2228 2229 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f); 2230 ierr = DMGetWorkArray(dm, f.numQuadPoints*pdim, PETSC_REAL, &Bf);CHKERRQ(ierr); 2231 ierr = PetscSpaceEvaluate(fem->basisSpace, f.numQuadPoints, f.quadPoints, Bf, NULL, NULL);CHKERRQ(ierr); 2232 for (k = 0; k < pdim; ++k) { 2233 /* n_j \cdot \phi_k */ 2234 invV[j*pdim+k] = 0.0; 2235 for (q = 0; q < f.numQuadPoints; ++q) { 2236 invV[j*pdim+k] += Bf[q*pdim+k]*f.quadWeights[q]; 2237 } 2238 } 2239 ierr = DMRestoreWorkArray(dm, f.numQuadPoints*pdim, PETSC_REAL, &Bf);CHKERRQ(ierr); 2240 } 2241 { 2242 PetscReal *work; 2243 PetscBLASInt *pivots; 2244 PetscBLASInt n = pdim, info; 2245 2246 ierr = DMGetWorkArray(dm, pdim, PETSC_INT, &pivots);CHKERRQ(ierr); 2247 ierr = DMGetWorkArray(dm, pdim, PETSC_REAL, &work);CHKERRQ(ierr); 2248 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invV, &n, pivots, &info)); 2249 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invV, &n, pivots, work, &n, &info)); 2250 ierr = DMRestoreWorkArray(dm, pdim, PETSC_INT, &pivots);CHKERRQ(ierr); 2251 ierr = DMRestoreWorkArray(dm, pdim, PETSC_REAL, &work);CHKERRQ(ierr); 2252 } 2253 for (p = 0; p < npoints; ++p) { 2254 if (B) { 2255 /* Multiply by V^{-1} (pdim x pdim) */ 2256 for (j = 0; j < pdim; ++j) { 2257 const PetscInt i = (p*pdim + j)*comp; 2258 PetscInt c; 2259 2260 (*B)[i] = 0.0; 2261 for (k = 0; k < pdim; ++k) { 2262 (*B)[i] += invV[k*pdim+j] * tmpB[p*pdim + k]; 2263 } 2264 for (c = 1; c < comp; ++c) { 2265 (*B)[i+c] = (*B)[i]; 2266 } 2267 } 2268 } 2269 if (D) { 2270 /* Multiply by V^{-1} (pdim x pdim) */ 2271 for (j = 0; j < pdim; ++j) { 2272 for (d = 0; d < dim; ++d) { 2273 const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d; 2274 PetscInt c; 2275 2276 (*D)[i] = 0.0; 2277 for (k = 0; k < pdim; ++k) { 2278 (*D)[i] += invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d]; 2279 } 2280 for (c = 1; c < comp; ++c) { 2281 (*D)[((p*pdim + j)*comp + c)*dim + d] = (*D)[i]; 2282 } 2283 } 2284 } 2285 } 2286 } 2287 ierr = DMRestoreWorkArray(dm, pdim*pdim, PETSC_REAL, &invV);CHKERRQ(ierr); 2288 if (B) {ierr = DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);CHKERRQ(ierr);} 2289 if (D) {ierr = DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);CHKERRQ(ierr);} 2290 PetscFunctionReturn(0); 2291 } 2292 2293 #undef __FUNCT__ 2294 #define __FUNCT__ "PetscFERestoreTabulation" 2295 PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H) 2296 { 2297 DM dm; 2298 PetscErrorCode ierr; 2299 2300 PetscFunctionBegin; 2301 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2302 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2303 if (B) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, B);CHKERRQ(ierr);} 2304 if (D) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, D);CHKERRQ(ierr);} 2305 if (H) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, H);CHKERRQ(ierr);} 2306 PetscFunctionReturn(0); 2307 } 2308 2309 /* 2310 Purpose: Compute element vector for chunk of elements 2311 2312 Input: 2313 Sizes: 2314 Ne: number of elements 2315 Nf: number of fields 2316 PetscFE 2317 dim: spatial dimension 2318 Nb: number of basis functions 2319 Nc: number of field components 2320 PetscQuadrature 2321 Nq: number of quadrature points 2322 2323 Geometry: 2324 PetscCellGeometry 2325 PetscReal v0s[Ne*dim] 2326 PetscReal jacobians[Ne*dim*dim] possibly *Nq 2327 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 2328 PetscReal jacobianDeterminants[Ne] possibly *Nq 2329 FEM: 2330 PetscFE 2331 PetscQuadrature 2332 PetscReal quadPoints[Nq*dim] 2333 PetscReal quadWeights[Nq] 2334 PetscReal basis[Nq*Nb*Nc] 2335 PetscReal basisDer[Nq*Nb*Nc*dim] 2336 PetscScalar coefficients[Ne*Nb*Nc] 2337 PetscScalar elemVec[Ne*Nb*Nc] 2338 2339 Problem: 2340 PetscInt f: the active field 2341 f0, f1 2342 2343 Work Space: 2344 PetscFE 2345 PetscScalar f0[Nq*dim]; 2346 PetscScalar f1[Nq*dim*dim]; 2347 PetscScalar u[Nc]; 2348 PetscScalar gradU[Nc*dim]; 2349 PetscReal x[dim]; 2350 PetscScalar realSpaceDer[dim]; 2351 2352 Purpose: Compute element vector for N_cb batches of elements 2353 2354 Input: 2355 Sizes: 2356 N_cb: Number of serial cell batches 2357 2358 Geometry: 2359 PetscReal v0s[Ne*dim] 2360 PetscReal jacobians[Ne*dim*dim] possibly *Nq 2361 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 2362 PetscReal jacobianDeterminants[Ne] possibly *Nq 2363 FEM: 2364 static PetscReal quadPoints[Nq*dim] 2365 static PetscReal quadWeights[Nq] 2366 static PetscReal basis[Nq*Nb*Nc] 2367 static PetscReal basisDer[Nq*Nb*Nc*dim] 2368 PetscScalar coefficients[Ne*Nb*Nc] 2369 PetscScalar elemVec[Ne*Nb*Nc] 2370 2371 ex62.c: 2372 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 2373 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 2374 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 2375 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 2376 2377 ex52.c: 2378 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) 2379 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) 2380 2381 ex52_integrateElement.cu 2382 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 2383 2384 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 2385 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 2386 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 2387 2388 ex52_integrateElementOpenCL.c: 2389 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 2390 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 2391 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 2392 2393 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 2394 */ 2395 2396 #undef __FUNCT__ 2397 #define __FUNCT__ "PetscFEIntegrateResidualChunk" 2398 /*C 2399 PetscFEIntegrateResidualChunk - Produce the element residual vector for a batch of elements by quadrature integration 2400 2401 Not collective 2402 2403 Input Parameters: 2404 + dim - The spatial dimension 2405 . Nf - The number of physical fields 2406 . Nc - The total number of fields components 2407 . field - The field being integrated 2408 . Ne - The number of elements 2409 . Nq - The number of quadrature points in this field 2410 . Nb - The number of basis functions in this field 2411 . v0s - The coordinates of the initial vertex for each element (the constant part of the transform from the reference element) 2412 . jacobians - The Jacobian for each element (the linear part of the transform from the reference element) 2413 . jacobianInverses - The Jacobian inverse for each element (the linear part of the transform to the reference element) 2414 . jacobianDeterminants - The Jacobian determinant for each element 2415 2416 . quadPoints - The quadrature points 2417 . quadWeights - The quadrature weights 2418 . coefficients - The array of FEM basis coefficients for the elements 2419 . f0_func - f_0 function from the first order FEM model 2420 - f1_func - f_1 function from the first order FEM model 2421 2422 Output Parameter 2423 . elemVec - the element residual vectors from each element 2424 2425 Calling sequence of f0_func and f1_func: 2426 $ void f0(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]) 2427 2428 Note: 2429 $ Loop over batch of elements (e): 2430 $ Loop over quadrature points (q): 2431 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q 2432 $ Call f_0 and f_1 2433 $ Loop over element vector entries (f,fc --> i): 2434 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 2435 */ 2436 PetscErrorCode PetscFEIntegrateResidualChunk(PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt field, PetscCellGeometry geom, const PetscScalar coefficients[], 2437 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[]), 2438 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f1[]), 2439 PetscScalar elemVec[]) 2440 { 2441 const PetscInt debug = 0; 2442 PetscQuadrature quad; 2443 PetscScalar *f0, *f1, *u, *gradU; 2444 PetscReal *x, *realSpaceDer; 2445 PetscInt dim, numComponents = 0, cOffset = 0, eOffset = 0, e, f; 2446 PetscErrorCode ierr; 2447 2448 PetscFunctionBegin; 2449 ierr = PetscFEGetSpatialDimension(fe[0], &dim);CHKERRQ(ierr); 2450 for (f = 0; f < Nf; ++f) { 2451 PetscInt Nc; 2452 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 2453 numComponents += Nc; 2454 } 2455 ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr); 2456 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); 2457 for (e = 0; e < Ne; ++e) { 2458 const PetscReal detJ = geom.detJ[e]; 2459 const PetscReal *v0 = &geom.v0[e*dim]; 2460 const PetscReal *J = &geom.J[e*dim*dim]; 2461 const PetscReal *invJ = &geom.invJ[e*dim*dim]; 2462 const PetscInt Nq = quad.numQuadPoints; 2463 PetscInt q, f; 2464 2465 if (debug > 1) { 2466 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 2467 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 2468 } 2469 for (q = 0; q < Nq; ++q) { 2470 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 2471 PetscInt fOffset = 0; 2472 PetscInt dOffset = cOffset; 2473 const PetscReal *quadPoints = quad.quadPoints; 2474 const PetscReal *quadWeights = quad.quadWeights; 2475 PetscInt Ncomp, d, d2, f, i; 2476 2477 ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr); 2478 for (d = 0; d < numComponents; ++d) {u[d] = 0.0;} 2479 for (d = 0; d < dim*(numComponents); ++d) {gradU[d] = 0.0;} 2480 for (d = 0; d < dim; ++d) { 2481 x[d] = v0[d]; 2482 for (d2 = 0; d2 < dim; ++d2) { 2483 x[d] += J[d*dim+d2]*(quadPoints[q*dim+d2] + 1.0); 2484 } 2485 } 2486 for (f = 0; f < Nf; ++f) { 2487 PetscReal *basis, *basisDer; 2488 PetscInt Nb, Ncomp, b, comp; 2489 2490 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2491 ierr = PetscFEGetNumComponents(fe[f], &Ncomp);CHKERRQ(ierr); 2492 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2493 ierr = PetscFEGetTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2494 for (b = 0; b < Nb; ++b) { 2495 for (comp = 0; comp < Ncomp; ++comp) { 2496 const PetscInt cidx = b*Ncomp+comp; 2497 PetscInt d, g; 2498 2499 u[fOffset+comp] += coefficients[dOffset+cidx]*basis[q*Nb*Ncomp+cidx]; 2500 for (d = 0; d < dim; ++d) { 2501 realSpaceDer[d] = 0.0; 2502 for (g = 0; g < dim; ++g) { 2503 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g]; 2504 } 2505 gradU[(fOffset+comp)*dim+d] += coefficients[dOffset+cidx]*realSpaceDer[d]; 2506 } 2507 } 2508 } 2509 ierr = PetscFERestoreTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2510 if (debug > 1) { 2511 PetscInt d; 2512 for (comp = 0; comp < Ncomp; ++comp) { 2513 ierr = PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, comp, u[fOffset+comp]);CHKERRQ(ierr); 2514 for (d = 0; d < dim; ++d) { 2515 ierr = PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, comp, 'x'+d, gradU[(fOffset+comp)*dim+d]);CHKERRQ(ierr); 2516 } 2517 } 2518 } 2519 fOffset += Ncomp; 2520 dOffset += Nb*Ncomp; 2521 } 2522 2523 f0_func(u, gradU, NULL, NULL, x, &f0[q*Ncomp]); 2524 for (i = 0; i < Ncomp; ++i) { 2525 f0[q*Ncomp+i] *= detJ*quadWeights[q]; 2526 } 2527 f1_func(u, gradU, NULL, NULL, x, &f1[q*Ncomp*dim]); 2528 for (i = 0; i < Ncomp*dim; ++i) { 2529 f1[q*Ncomp*dim+i] *= detJ*quadWeights[q]; 2530 } 2531 if (debug > 1) { 2532 PetscInt c,d; 2533 for (c = 0; c < Ncomp; ++c) { 2534 ierr = PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, f0[q*Ncomp+c]);CHKERRQ(ierr); 2535 for (d = 0; d < dim; ++d) { 2536 ierr = PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, f1[(q*Ncomp + c)*dim+d]);CHKERRQ(ierr); 2537 } 2538 } 2539 } 2540 if (q == Nq-1) {cOffset = dOffset;} 2541 } 2542 for (f = 0; f < Nf; ++f) { 2543 PetscInt Nb, Ncomp, b, comp; 2544 2545 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2546 ierr = PetscFEGetNumComponents(fe[f], &Ncomp);CHKERRQ(ierr); 2547 if (f == field) { 2548 PetscReal *basis; 2549 PetscReal *basisDer; 2550 2551 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2552 ierr = PetscFEGetTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2553 for (b = 0; b < Nb; ++b) { 2554 for (comp = 0; comp < Ncomp; ++comp) { 2555 const PetscInt cidx = b*Ncomp+comp; 2556 PetscInt q; 2557 2558 elemVec[eOffset+cidx] = 0.0; 2559 for (q = 0; q < Nq; ++q) { 2560 PetscInt d, g; 2561 2562 elemVec[eOffset+cidx] += basis[q*Nb*Ncomp+cidx]*f0[q*Ncomp+comp]; 2563 for (d = 0; d < dim; ++d) { 2564 realSpaceDer[d] = 0.0; 2565 for (g = 0; g < dim; ++g) { 2566 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g]; 2567 } 2568 elemVec[eOffset+cidx] += realSpaceDer[d]*f1[(q*Ncomp+comp)*dim+d]; 2569 } 2570 } 2571 } 2572 } 2573 ierr = PetscFERestoreTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2574 if (debug > 1) { 2575 PetscInt b, comp; 2576 2577 for (b = 0; b < Nb; ++b) { 2578 for (comp = 0; comp < Ncomp; ++comp) { 2579 ierr = PetscPrintf(PETSC_COMM_SELF, " elemVec[%d,%d]: %g\n", b, comp, elemVec[eOffset+b*Ncomp+comp]);CHKERRQ(ierr); 2580 } 2581 } 2582 } 2583 } 2584 eOffset += Nb*Ncomp; 2585 } 2586 } 2587 ierr = PetscFree6(f0,f1,u,gradU,x,realSpaceDer); 2588 PetscFunctionReturn(0); 2589 } 2590