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 = PetscFree((*fem)->numDof);CHKERRQ(ierr); 1983 ierr = PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);CHKERRQ(ierr); 1984 ierr = PetscSpaceDestroy(&(*fem)->basisSpace);CHKERRQ(ierr); 1985 ierr = PetscDualSpaceDestroy(&(*fem)->dualSpace);CHKERRQ(ierr); 1986 ierr = PetscQuadratureDestroy(&(*fem)->quadrature);CHKERRQ(ierr); 1987 1988 if ((*fem)->ops->destroy) {ierr = (*(*fem)->ops->destroy)(*fem);CHKERRQ(ierr);} 1989 ierr = PetscHeaderDestroy(fem);CHKERRQ(ierr); 1990 PetscFunctionReturn(0); 1991 } 1992 1993 #undef __FUNCT__ 1994 #define __FUNCT__ "PetscFECreate" 1995 /*@ 1996 PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType(). 1997 1998 Collective on MPI_Comm 1999 2000 Input Parameter: 2001 . comm - The communicator for the PetscFE object 2002 2003 Output Parameter: 2004 . fem - The PetscFE object 2005 2006 Level: beginner 2007 2008 .seealso: PetscFESetType(), PETSCFEGALERKIN 2009 @*/ 2010 PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem) 2011 { 2012 PetscFE f; 2013 PetscErrorCode ierr; 2014 2015 PetscFunctionBegin; 2016 PetscValidPointer(fem, 2); 2017 *fem = NULL; 2018 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 2019 ierr = PetscFEInitializePackage();CHKERRQ(ierr); 2020 #endif 2021 2022 ierr = PetscHeaderCreate(f, _p_PetscFE, struct _PetscFEOps, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);CHKERRQ(ierr); 2023 ierr = PetscMemzero(f->ops, sizeof(struct _PetscFEOps));CHKERRQ(ierr); 2024 2025 f->basisSpace = NULL; 2026 f->dualSpace = NULL; 2027 f->numComponents = 1; 2028 f->numDof = NULL; 2029 f->B = NULL; 2030 f->D = NULL; 2031 f->H = NULL; 2032 ierr = PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));CHKERRQ(ierr); 2033 2034 *fem = f; 2035 PetscFunctionReturn(0); 2036 } 2037 2038 #undef __FUNCT__ 2039 #define __FUNCT__ "PetscFEGetDimension" 2040 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 2041 { 2042 PetscErrorCode ierr; 2043 2044 PetscFunctionBegin; 2045 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2046 PetscValidPointer(dim, 2); 2047 ierr = PetscSpaceGetDimension(fem->basisSpace, dim);CHKERRQ(ierr); 2048 PetscFunctionReturn(0); 2049 } 2050 2051 #undef __FUNCT__ 2052 #define __FUNCT__ "PetscFEGetSpatialDimension" 2053 PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim) 2054 { 2055 DM dm; 2056 PetscErrorCode ierr; 2057 2058 PetscFunctionBegin; 2059 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2060 PetscValidPointer(dim, 2); 2061 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2062 ierr = DMPlexGetDimension(dm, dim);CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 #undef __FUNCT__ 2067 #define __FUNCT__ "PetscFESetNumComponents" 2068 PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp) 2069 { 2070 PetscFunctionBegin; 2071 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2072 fem->numComponents = comp; 2073 PetscFunctionReturn(0); 2074 } 2075 2076 #undef __FUNCT__ 2077 #define __FUNCT__ "PetscFEGetNumComponents" 2078 PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp) 2079 { 2080 PetscFunctionBegin; 2081 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2082 PetscValidPointer(comp, 2); 2083 *comp = fem->numComponents; 2084 PetscFunctionReturn(0); 2085 } 2086 2087 #undef __FUNCT__ 2088 #define __FUNCT__ "PetscFEGetBasisSpace" 2089 PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp) 2090 { 2091 PetscFunctionBegin; 2092 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2093 PetscValidPointer(sp, 2); 2094 *sp = fem->basisSpace; 2095 PetscFunctionReturn(0); 2096 } 2097 2098 #undef __FUNCT__ 2099 #define __FUNCT__ "PetscFESetBasisSpace" 2100 PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp) 2101 { 2102 PetscErrorCode ierr; 2103 2104 PetscFunctionBegin; 2105 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2106 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 2); 2107 ierr = PetscSpaceDestroy(&fem->basisSpace);CHKERRQ(ierr); 2108 fem->basisSpace = sp; 2109 ierr = PetscObjectReference((PetscObject) fem->basisSpace);CHKERRQ(ierr); 2110 PetscFunctionReturn(0); 2111 } 2112 2113 #undef __FUNCT__ 2114 #define __FUNCT__ "PetscFEGetDualSpace" 2115 PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp) 2116 { 2117 PetscFunctionBegin; 2118 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2119 PetscValidPointer(sp, 2); 2120 *sp = fem->dualSpace; 2121 PetscFunctionReturn(0); 2122 } 2123 2124 #undef __FUNCT__ 2125 #define __FUNCT__ "PetscFESetDualSpace" 2126 PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp) 2127 { 2128 PetscErrorCode ierr; 2129 2130 PetscFunctionBegin; 2131 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2132 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 2); 2133 ierr = PetscDualSpaceDestroy(&fem->dualSpace);CHKERRQ(ierr); 2134 fem->dualSpace = sp; 2135 ierr = PetscObjectReference((PetscObject) fem->dualSpace);CHKERRQ(ierr); 2136 PetscFunctionReturn(0); 2137 } 2138 2139 #undef __FUNCT__ 2140 #define __FUNCT__ "PetscFEGetQuadrature" 2141 PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q) 2142 { 2143 PetscFunctionBegin; 2144 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2145 PetscValidPointer(q, 2); 2146 *q = fem->quadrature; 2147 PetscFunctionReturn(0); 2148 } 2149 2150 #undef __FUNCT__ 2151 #define __FUNCT__ "PetscFESetQuadrature" 2152 PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q) 2153 { 2154 PetscErrorCode ierr; 2155 2156 PetscFunctionBegin; 2157 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2158 ierr = PetscQuadratureDestroy(&fem->quadrature);CHKERRQ(ierr); 2159 fem->quadrature = q; 2160 PetscFunctionReturn(0); 2161 } 2162 2163 #undef __FUNCT__ 2164 #define __FUNCT__ "PetscFEGetNumDof" 2165 PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof) 2166 { 2167 const PetscInt *numDofDual; 2168 PetscErrorCode ierr; 2169 2170 PetscFunctionBegin; 2171 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2172 PetscValidPointer(numDof, 2); 2173 ierr = PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);CHKERRQ(ierr); 2174 if (!fem->numDof) { 2175 DM dm; 2176 PetscInt dim, d; 2177 2178 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2179 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2180 ierr = PetscMalloc((dim+1) * sizeof(PetscInt), &fem->numDof);CHKERRQ(ierr); 2181 for (d = 0; d <= dim; ++d) { 2182 fem->numDof[d] = fem->numComponents*numDofDual[d]; 2183 } 2184 } 2185 *numDof = fem->numDof; 2186 PetscFunctionReturn(0); 2187 } 2188 2189 #undef __FUNCT__ 2190 #define __FUNCT__ "PetscFEGetDefaultTabulation" 2191 PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H) 2192 { 2193 PetscInt npoints = fem->quadrature.numQuadPoints; 2194 const PetscReal *points = fem->quadrature.quadPoints; 2195 PetscErrorCode ierr; 2196 2197 PetscFunctionBegin; 2198 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2199 if (B) PetscValidPointer(B, 2); 2200 if (D) PetscValidPointer(D, 3); 2201 if (H) PetscValidPointer(H, 4); 2202 if (!fem->B) {ierr = PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);CHKERRQ(ierr);} 2203 if (B) *B = fem->B; 2204 if (D) *D = fem->D; 2205 if (H) *H = fem->H; 2206 PetscFunctionReturn(0); 2207 } 2208 2209 #undef __FUNCT__ 2210 #define __FUNCT__ "PetscFEGetTabulation" 2211 PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H) 2212 { 2213 DM dm; 2214 PetscInt pdim; /* Dimension of FE space P */ 2215 PetscInt dim; /* Spatial dimension */ 2216 PetscInt comp; /* Field components */ 2217 PetscReal *tmpB, *tmpD, *invV; 2218 PetscInt p, d, j, k; 2219 PetscErrorCode ierr; 2220 2221 PetscFunctionBegin; 2222 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2223 PetscValidPointer(points, 3); 2224 if (B) PetscValidPointer(B, 4); 2225 if (D) PetscValidPointer(D, 5); 2226 if (H) PetscValidPointer(H, 6); 2227 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2228 2229 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2230 ierr = PetscSpaceGetDimension(fem->basisSpace, &pdim);CHKERRQ(ierr); 2231 ierr = PetscFEGetNumComponents(fem, &comp);CHKERRQ(ierr); 2232 /* 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); */ 2233 2234 if (B) { 2235 ierr = DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);CHKERRQ(ierr); 2236 ierr = DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);CHKERRQ(ierr); 2237 } 2238 if (D) { 2239 ierr = DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);CHKERRQ(ierr); 2240 ierr = DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);CHKERRQ(ierr); 2241 } 2242 if (H) {ierr = DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, H);CHKERRQ(ierr);} 2243 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? *H : NULL);CHKERRQ(ierr); 2244 2245 ierr = DMGetWorkArray(dm, pdim*pdim, PETSC_REAL, &invV);CHKERRQ(ierr); 2246 for (j = 0; j < pdim; ++j) { 2247 PetscReal *Bf; 2248 PetscQuadrature f; 2249 PetscInt q; 2250 2251 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f); 2252 ierr = DMGetWorkArray(dm, f.numQuadPoints*pdim, PETSC_REAL, &Bf);CHKERRQ(ierr); 2253 ierr = PetscSpaceEvaluate(fem->basisSpace, f.numQuadPoints, f.quadPoints, Bf, NULL, NULL);CHKERRQ(ierr); 2254 for (k = 0; k < pdim; ++k) { 2255 /* n_j \cdot \phi_k */ 2256 invV[j*pdim+k] = 0.0; 2257 for (q = 0; q < f.numQuadPoints; ++q) { 2258 invV[j*pdim+k] += Bf[q*pdim+k]*f.quadWeights[q]; 2259 } 2260 } 2261 ierr = DMRestoreWorkArray(dm, f.numQuadPoints*pdim, PETSC_REAL, &Bf);CHKERRQ(ierr); 2262 } 2263 { 2264 PetscReal *work; 2265 PetscBLASInt *pivots; 2266 PetscBLASInt n = pdim, info; 2267 2268 ierr = DMGetWorkArray(dm, pdim, PETSC_INT, &pivots);CHKERRQ(ierr); 2269 ierr = DMGetWorkArray(dm, pdim, PETSC_REAL, &work);CHKERRQ(ierr); 2270 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invV, &n, pivots, &info)); 2271 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invV, &n, pivots, work, &n, &info)); 2272 ierr = DMRestoreWorkArray(dm, pdim, PETSC_INT, &pivots);CHKERRQ(ierr); 2273 ierr = DMRestoreWorkArray(dm, pdim, PETSC_REAL, &work);CHKERRQ(ierr); 2274 } 2275 for (p = 0; p < npoints; ++p) { 2276 if (B) { 2277 /* Multiply by V^{-1} (pdim x pdim) */ 2278 for (j = 0; j < pdim; ++j) { 2279 const PetscInt i = (p*pdim + j)*comp; 2280 PetscInt c; 2281 2282 (*B)[i] = 0.0; 2283 for (k = 0; k < pdim; ++k) { 2284 (*B)[i] += invV[k*pdim+j] * tmpB[p*pdim + k]; 2285 } 2286 for (c = 1; c < comp; ++c) { 2287 (*B)[i+c] = (*B)[i]; 2288 } 2289 } 2290 } 2291 if (D) { 2292 /* Multiply by V^{-1} (pdim x pdim) */ 2293 for (j = 0; j < pdim; ++j) { 2294 for (d = 0; d < dim; ++d) { 2295 const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d; 2296 PetscInt c; 2297 2298 (*D)[i] = 0.0; 2299 for (k = 0; k < pdim; ++k) { 2300 (*D)[i] += invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d]; 2301 } 2302 for (c = 1; c < comp; ++c) { 2303 (*D)[((p*pdim + j)*comp + c)*dim + d] = (*D)[i]; 2304 } 2305 } 2306 } 2307 } 2308 } 2309 ierr = DMRestoreWorkArray(dm, pdim*pdim, PETSC_REAL, &invV);CHKERRQ(ierr); 2310 if (B) {ierr = DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);CHKERRQ(ierr);} 2311 if (D) {ierr = DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);CHKERRQ(ierr);} 2312 PetscFunctionReturn(0); 2313 } 2314 2315 #undef __FUNCT__ 2316 #define __FUNCT__ "PetscFERestoreTabulation" 2317 PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H) 2318 { 2319 DM dm; 2320 PetscErrorCode ierr; 2321 2322 PetscFunctionBegin; 2323 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 2324 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 2325 if (B) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, B);CHKERRQ(ierr);} 2326 if (D) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, D);CHKERRQ(ierr);} 2327 if (H) {ierr = DMRestoreWorkArray(dm, 0, PETSC_REAL, H);CHKERRQ(ierr);} 2328 PetscFunctionReturn(0); 2329 } 2330 2331 /* 2332 Purpose: Compute element vector for chunk of elements 2333 2334 Input: 2335 Sizes: 2336 Ne: number of elements 2337 Nf: number of fields 2338 PetscFE 2339 dim: spatial dimension 2340 Nb: number of basis functions 2341 Nc: number of field components 2342 PetscQuadrature 2343 Nq: number of quadrature points 2344 2345 Geometry: 2346 PetscCellGeometry 2347 PetscReal v0s[Ne*dim] 2348 PetscReal jacobians[Ne*dim*dim] possibly *Nq 2349 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 2350 PetscReal jacobianDeterminants[Ne] possibly *Nq 2351 FEM: 2352 PetscFE 2353 PetscQuadrature 2354 PetscReal quadPoints[Nq*dim] 2355 PetscReal quadWeights[Nq] 2356 PetscReal basis[Nq*Nb*Nc] 2357 PetscReal basisDer[Nq*Nb*Nc*dim] 2358 PetscScalar coefficients[Ne*Nb*Nc] 2359 PetscScalar elemVec[Ne*Nb*Nc] 2360 2361 Problem: 2362 PetscInt f: the active field 2363 f0, f1 2364 2365 Work Space: 2366 PetscFE 2367 PetscScalar f0[Nq*dim]; 2368 PetscScalar f1[Nq*dim*dim]; 2369 PetscScalar u[Nc]; 2370 PetscScalar gradU[Nc*dim]; 2371 PetscReal x[dim]; 2372 PetscScalar realSpaceDer[dim]; 2373 2374 Purpose: Compute element vector for N_cb batches of elements 2375 2376 Input: 2377 Sizes: 2378 N_cb: Number of serial cell batches 2379 2380 Geometry: 2381 PetscReal v0s[Ne*dim] 2382 PetscReal jacobians[Ne*dim*dim] possibly *Nq 2383 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 2384 PetscReal jacobianDeterminants[Ne] possibly *Nq 2385 FEM: 2386 static PetscReal quadPoints[Nq*dim] 2387 static PetscReal quadWeights[Nq] 2388 static PetscReal basis[Nq*Nb*Nc] 2389 static PetscReal basisDer[Nq*Nb*Nc*dim] 2390 PetscScalar coefficients[Ne*Nb*Nc] 2391 PetscScalar elemVec[Ne*Nb*Nc] 2392 2393 ex62.c: 2394 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 2395 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 2396 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 2397 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 2398 2399 ex52.c: 2400 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) 2401 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) 2402 2403 ex52_integrateElement.cu 2404 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 2405 2406 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 2407 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 2408 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 2409 2410 ex52_integrateElementOpenCL.c: 2411 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 2412 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 2413 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 2414 2415 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 2416 */ 2417 2418 #undef __FUNCT__ 2419 #define __FUNCT__ "PetscFEIntegrateResidualChunk" 2420 /*C 2421 PetscFEIntegrateResidualChunk - Produce the element residual vector for a chunk of elements by quadrature integration 2422 2423 Not collective 2424 2425 Input Parameters: 2426 + Ne - The number of elements in the chunk 2427 . Nf - The number of physical fields 2428 . fe - The PetscFE objects for each field 2429 . field - The field being integrated 2430 . geom - The cell geometry for each cell in the chunk 2431 . coefficients - The array of FEM basis coefficients for the elements 2432 . f0_func - f_0 function from the first order FEM model 2433 - f1_func - f_1 function from the first order FEM model 2434 2435 Output Parameter 2436 . elemVec - the element residual vectors from each element 2437 2438 Calling sequence of f0_func and f1_func: 2439 $ void f0(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[]) 2440 2441 Note: 2442 $ Loop over batch of elements (e): 2443 $ Loop over quadrature points (q): 2444 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q 2445 $ Call f_0 and f_1 2446 $ Loop over element vector entries (f,fc --> i): 2447 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 2448 */ 2449 PetscErrorCode PetscFEIntegrateResidualChunk(PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt field, PetscCellGeometry geom, const PetscScalar coefficients[], 2450 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[]), 2451 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f1[]), 2452 PetscScalar elemVec[]) 2453 { 2454 const PetscInt debug = 0; 2455 PetscQuadrature quad; 2456 PetscScalar *f0, *f1, *u, *gradU; 2457 PetscReal *x, *realSpaceDer; 2458 PetscInt dim, numComponents = 0, cOffset = 0, eOffset = 0, e, f; 2459 PetscErrorCode ierr; 2460 2461 PetscFunctionBegin; 2462 ierr = PetscFEGetSpatialDimension(fe[0], &dim);CHKERRQ(ierr); 2463 for (f = 0; f < Nf; ++f) { 2464 PetscInt Nc; 2465 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 2466 numComponents += Nc; 2467 } 2468 ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr); 2469 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); 2470 for (e = 0; e < Ne; ++e) { 2471 const PetscReal detJ = geom.detJ[e]; 2472 const PetscReal *v0 = &geom.v0[e*dim]; 2473 const PetscReal *J = &geom.J[e*dim*dim]; 2474 const PetscReal *invJ = &geom.invJ[e*dim*dim]; 2475 const PetscInt Nq = quad.numQuadPoints; 2476 const PetscReal *quadPoints = quad.quadPoints; 2477 const PetscReal *quadWeights = quad.quadWeights; 2478 PetscInt q, f; 2479 2480 if (debug > 1) { 2481 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 2482 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 2483 } 2484 for (q = 0; q < Nq; ++q) { 2485 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 2486 PetscInt fOffset = 0; 2487 PetscInt dOffset = cOffset; 2488 PetscInt Ncomp, d, d2, f, i; 2489 2490 ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr); 2491 for (d = 0; d < numComponents; ++d) {u[d] = 0.0;} 2492 for (d = 0; d < dim*(numComponents); ++d) {gradU[d] = 0.0;} 2493 for (d = 0; d < dim; ++d) { 2494 x[d] = v0[d]; 2495 for (d2 = 0; d2 < dim; ++d2) { 2496 x[d] += J[d*dim+d2]*(quadPoints[q*dim+d2] + 1.0); 2497 } 2498 } 2499 for (f = 0; f < Nf; ++f) { 2500 PetscReal *basis, *basisDer; 2501 PetscInt Nb, Ncomp, b, comp; 2502 2503 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2504 ierr = PetscFEGetNumComponents(fe[f], &Ncomp);CHKERRQ(ierr); 2505 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2506 ierr = PetscFEGetDefaultTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2507 for (b = 0; b < Nb; ++b) { 2508 for (comp = 0; comp < Ncomp; ++comp) { 2509 const PetscInt cidx = b*Ncomp+comp; 2510 PetscInt d, g; 2511 2512 u[fOffset+comp] += coefficients[dOffset+cidx]*basis[q*Nb*Ncomp+cidx]; 2513 for (d = 0; d < dim; ++d) { 2514 realSpaceDer[d] = 0.0; 2515 for (g = 0; g < dim; ++g) { 2516 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g]; 2517 } 2518 gradU[(fOffset+comp)*dim+d] += coefficients[dOffset+cidx]*realSpaceDer[d]; 2519 } 2520 } 2521 } 2522 if (debug > 1) { 2523 PetscInt d; 2524 for (comp = 0; comp < Ncomp; ++comp) { 2525 ierr = PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, comp, u[fOffset+comp]);CHKERRQ(ierr); 2526 for (d = 0; d < dim; ++d) { 2527 ierr = PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, comp, 'x'+d, gradU[(fOffset+comp)*dim+d]);CHKERRQ(ierr); 2528 } 2529 } 2530 } 2531 fOffset += Ncomp; 2532 dOffset += Nb*Ncomp; 2533 } 2534 2535 f0_func(u, gradU, NULL, NULL, x, &f0[q*Ncomp]); 2536 for (i = 0; i < Ncomp; ++i) { 2537 f0[q*Ncomp+i] *= detJ*quadWeights[q]; 2538 } 2539 f1_func(u, gradU, NULL, NULL, x, &f1[q*Ncomp*dim]); 2540 for (i = 0; i < Ncomp*dim; ++i) { 2541 f1[q*Ncomp*dim+i] *= detJ*quadWeights[q]; 2542 } 2543 if (debug > 1) { 2544 PetscInt c,d; 2545 for (c = 0; c < Ncomp; ++c) { 2546 ierr = PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, f0[q*Ncomp+c]);CHKERRQ(ierr); 2547 for (d = 0; d < dim; ++d) { 2548 ierr = PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, f1[(q*Ncomp + c)*dim+d]);CHKERRQ(ierr); 2549 } 2550 } 2551 } 2552 if (q == Nq-1) {cOffset = dOffset;} 2553 } 2554 for (f = 0; f < Nf; ++f) { 2555 PetscInt Nb, Ncomp, b, comp; 2556 2557 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2558 ierr = PetscFEGetNumComponents(fe[f], &Ncomp);CHKERRQ(ierr); 2559 if (f == field) { 2560 PetscReal *basis; 2561 PetscReal *basisDer; 2562 2563 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2564 ierr = PetscFEGetDefaultTabulation(fe[f], &basis, &basisDer, NULL);CHKERRQ(ierr); 2565 for (b = 0; b < Nb; ++b) { 2566 for (comp = 0; comp < Ncomp; ++comp) { 2567 const PetscInt cidx = b*Ncomp+comp; 2568 PetscInt q; 2569 2570 elemVec[eOffset+cidx] = 0.0; 2571 for (q = 0; q < Nq; ++q) { 2572 PetscInt d, g; 2573 2574 elemVec[eOffset+cidx] += basis[q*Nb*Ncomp+cidx]*f0[q*Ncomp+comp]; 2575 for (d = 0; d < dim; ++d) { 2576 realSpaceDer[d] = 0.0; 2577 for (g = 0; g < dim; ++g) { 2578 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g]; 2579 } 2580 elemVec[eOffset+cidx] += realSpaceDer[d]*f1[(q*Ncomp+comp)*dim+d]; 2581 } 2582 } 2583 } 2584 } 2585 if (debug > 1) { 2586 PetscInt b, comp; 2587 2588 for (b = 0; b < Nb; ++b) { 2589 for (comp = 0; comp < Ncomp; ++comp) { 2590 ierr = PetscPrintf(PETSC_COMM_SELF, " elemVec[%d,%d]: %g\n", b, comp, elemVec[eOffset+b*Ncomp+comp]);CHKERRQ(ierr); 2591 } 2592 } 2593 } 2594 } 2595 eOffset += Nb*Ncomp; 2596 } 2597 } 2598 ierr = PetscFree6(f0,f1,u,gradU,x,realSpaceDer); 2599 PetscFunctionReturn(0); 2600 } 2601 2602 #undef __FUNCT__ 2603 #define __FUNCT__ "PetscFEIntegrateJacobianChunk" 2604 /*C 2605 PetscFEIntegrateJacobianChunk - Produce the element Jacobian for a chunk of elements by quadrature integration 2606 2607 Not collective 2608 2609 Input Parameters: 2610 + Ne - The number of elements in the chunk 2611 . Nf - The number of physical fields 2612 . fe - The PetscFE objects for each field 2613 . fieldI - The test field being integrated 2614 . fieldJ - The basis field being integrated 2615 . geom - The cell geometry for each cell in the chunk 2616 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 2617 . g0_func - g_0 function from the first order FEM model 2618 . g1_func - g_1 function from the first order FEM model 2619 . g2_func - g_2 function from the first order FEM model 2620 - g3_func - g_3 function from the first order FEM model 2621 2622 Output Parameter 2623 . elemMat - the element matrices for the Jacobian from each element 2624 2625 Calling sequence of g0_func, g1_func, g2_func and g3_func: 2626 $ void g0(PetscScalar u[], const PetscScalar gradU[], PetscScalar a[], const PetscScalar gradA[], PetscScalar x[], PetscScalar g0[]) 2627 2628 Note: 2629 $ Loop over batch of elements (e): 2630 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 2631 $ Loop over quadrature points (q): 2632 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 2633 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 2634 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 2635 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 2636 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 2637 */ 2638 PetscErrorCode PetscFEIntegrateJacobianChunk(PetscInt Ne, PetscInt Nf, PetscFE fe[], PetscInt fieldI, PetscInt fieldJ, PetscCellGeometry geom, const PetscScalar coefficients[], 2639 void (*g0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g0[]), 2640 void (*g1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g1[]), 2641 void (*g2_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g2[]), 2642 void (*g3_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]), 2643 PetscScalar elemMat[]) { 2644 const PetscInt debug = 0; 2645 PetscInt cellDof = 0; /* Total number of dof on a cell */ 2646 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 2647 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 2648 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 2649 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 2650 PetscQuadrature quad; 2651 PetscScalar *g0, *g1, *g2, *g3, *u, *gradU; 2652 PetscReal *x, *realSpaceDerI, *realSpaceDerJ; 2653 PetscReal *basisI, *basisDerI, *basisJ, *basisDerJ; 2654 PetscInt NbI, NcI, NbJ, NcJ, numComponents = 0; 2655 PetscInt dim, f, e; 2656 PetscErrorCode ierr; 2657 2658 PetscFunctionBegin; 2659 ierr = PetscFEGetSpatialDimension(fe[fieldI], &dim);CHKERRQ(ierr); 2660 ierr = PetscFEGetDefaultTabulation(fe[fieldI], &basisI, &basisDerI, NULL);CHKERRQ(ierr); 2661 ierr = PetscFEGetDefaultTabulation(fe[fieldJ], &basisJ, &basisDerJ, NULL);CHKERRQ(ierr); 2662 for (f = 0; f < Nf; ++f) { 2663 PetscInt Nb, Nc; 2664 2665 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 2666 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 2667 if (f == fieldI) {offsetI = cellDof; NbI = Nb; NcI = Nc;} 2668 if (f == fieldJ) {offsetJ = cellDof; NbJ = Nb; NcJ = Nc;} 2669 numComponents += Nc; 2670 cellDof += Nb*Nc; 2671 } 2672 ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr); 2673 ierr = PetscMalloc4(NcI*NcJ,PetscScalar,&g0,NcI*NcJ*dim,PetscScalar,&g1,NcI*NcJ*dim,PetscScalar,&g2,NcI*NcJ*dim*dim,PetscScalar,&g3); 2674 ierr = PetscMalloc5(numComponents,PetscScalar,&u,numComponents*dim,PetscScalar,&gradU,dim,PetscReal,&x,dim,PetscReal,&realSpaceDerI,dim,PetscReal,&realSpaceDerJ); 2675 for (e = 0; e < Ne; ++e) { 2676 const PetscReal detJ = geom.detJ[e]; 2677 const PetscReal *v0 = &geom.v0[e*dim]; 2678 const PetscReal *J = &geom.J[e*dim*dim]; 2679 const PetscReal *invJ = &geom.invJ[e*dim*dim]; 2680 const PetscInt Nq = quad.numQuadPoints; 2681 const PetscReal *quadPoints = quad.quadPoints; 2682 const PetscReal *quadWeights = quad.quadWeights; 2683 PetscInt f, g, q; 2684 2685 for (f = 0; f < NbI; ++f) { 2686 for (g = 0; g < NbJ; ++g) { 2687 for (q = 0; q < Nq; ++q) { 2688 PetscInt fOffset = 0; /* Offset into u[] for field_q (like offsetI) */ 2689 PetscInt dOffset = cOffset; /* Offset into coefficients[] for field_q */ 2690 PetscInt field_q, d, d2; 2691 PetscInt fc, gc, c; 2692 2693 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 2694 for (d = 0; d < numComponents; ++d) {u[d] = 0.0;} 2695 for (d = 0; d < dim*numComponents; ++d) {gradU[d] = 0.0;} 2696 for (d = 0; d < dim; ++d) { 2697 x[d] = v0[d]; 2698 for (d2 = 0; d2 < dim; ++d2) { 2699 x[d] += J[d*dim+d2]*(quadPoints[q*dim+d2] + 1.0); 2700 } 2701 } 2702 for (field_q = 0; field_q < Nf; ++field_q) { 2703 PetscReal *basis, *basisDer; 2704 PetscInt Nb, Ncomp, b, comp; 2705 2706 ierr = PetscFEGetDimension(fe[field_q], &Nb);CHKERRQ(ierr); 2707 ierr = PetscFEGetNumComponents(fe[field_q], &Ncomp);CHKERRQ(ierr); 2708 /* TODO: Hoist this tabulation out of the loops, maybe by memoizing */ 2709 ierr = PetscFEGetDefaultTabulation(fe[field_q], &basis, &basisDer, NULL);CHKERRQ(ierr); 2710 for (b = 0; b < Nb; ++b) { 2711 for (comp = 0; comp < Ncomp; ++comp) { 2712 const PetscInt cidx = b*Ncomp+comp; 2713 PetscInt d1, d2; 2714 2715 u[fOffset+comp] += coefficients[dOffset+cidx]*basis[q*Nb*Ncomp+cidx]; 2716 for (d1 = 0; d1 < dim; ++d1) { 2717 realSpaceDerI[d1] = 0.0; 2718 for (d2 = 0; d2 < dim; ++d2) { 2719 realSpaceDerI[d1] += invJ[d2*dim+d1]*basisDer[(q*Nb*Ncomp+cidx)*dim+d2]; 2720 } 2721 gradU[(fOffset+comp)*dim+d1] += coefficients[dOffset+cidx]*realSpaceDerI[d1]; 2722 } 2723 } 2724 } 2725 if (debug > 1) { 2726 for (comp = 0; comp < Ncomp; ++comp) { 2727 ierr = PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, comp, u[fOffset+comp]);CHKERRQ(ierr); 2728 for (d = 0; d < dim; ++d) { 2729 ierr = PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, comp, 'x'+d, gradU[(fOffset+comp)*dim+d]);CHKERRQ(ierr); 2730 } 2731 } 2732 } 2733 fOffset += Ncomp; 2734 dOffset += Nb*Ncomp; 2735 } 2736 2737 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 2738 ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 2739 ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 2740 ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 2741 if (g0_func) { 2742 g0_func(u, gradU, NULL, NULL, x, g0); 2743 for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];} 2744 } 2745 if (g1_func) { 2746 g1_func(u, gradU, NULL, NULL, x, g1); 2747 for (c = 0; c < NcI*NcJ*dim; ++c) {g1[c] *= detJ*quadWeights[q];} 2748 } 2749 if (g2_func) { 2750 g2_func(u, gradU, NULL, NULL, x, g2); 2751 for (c = 0; c < NcI*NcJ*dim; ++c) {g2[c] *= detJ*quadWeights[q];} 2752 } 2753 if (g3_func) { 2754 g3_func(u, gradU, NULL, NULL, x, g3); 2755 for (c = 0; c < NcI*NcJ*dim*dim; ++c) {g3[c] *= detJ*quadWeights[q];} 2756 } 2757 2758 for (fc = 0; fc < NcI; ++fc) { 2759 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2760 const PetscInt i = offsetI+fidx; /* Element matrix row */ 2761 for (gc = 0; gc < NcJ; ++gc) { 2762 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2763 const PetscInt j = offsetJ+gidx; /* Element matrix column */ 2764 PetscInt d, d2; 2765 2766 for (d = 0; d < dim; ++d) { 2767 realSpaceDerI[d] = 0.0; 2768 realSpaceDerJ[d] = 0.0; 2769 for (d2 = 0; d2 < dim; ++d2) { 2770 realSpaceDerI[d] += invJ[d2*dim+d]*basisDerI[(q*NbI*NcI+fidx)*dim+d2]; 2771 realSpaceDerJ[d] += invJ[d2*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2]; 2772 } 2773 } 2774 elemMat[eOffset+i*cellDof+j] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx]; 2775 for (d = 0; d < dim; ++d) { 2776 elemMat[eOffset+i*cellDof+j] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*realSpaceDerJ[d]; 2777 elemMat[eOffset+i*cellDof+j] += realSpaceDerI[d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx]; 2778 for (d2 = 0; d2 < dim; ++d2) { 2779 elemMat[eOffset+i*cellDof+j] += realSpaceDerI[d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*realSpaceDerJ[d2]; 2780 } 2781 } 2782 } 2783 } 2784 } 2785 } 2786 } 2787 if (debug > 1) { 2788 PetscInt fc, f, gc, g; 2789 2790 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 2791 for (fc = 0; fc < NcI; ++fc) { 2792 for (f = 0; f < NbI; ++f) { 2793 const PetscInt i = offsetI + f*NcI+fc; 2794 for (gc = 0; gc < NcJ; ++gc) { 2795 for (g = 0; g < NbJ; ++g) { 2796 const PetscInt j = offsetJ + g*NcJ+gc; 2797 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, elemMat[eOffset+i*cellDof+j]);CHKERRQ(ierr); 2798 } 2799 } 2800 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 2801 } 2802 } 2803 } 2804 cOffset += cellDof; 2805 eOffset += cellDof*cellDof; 2806 } 2807 ierr = PetscFree4(g0,g1,g2,g3); 2808 ierr = PetscFree5(u,gradU,x,realSpaceDerI,realSpaceDerJ); 2809 PetscFunctionReturn(0); 2810 } 2811