1. petsc
  2. PETSc
  3. petsc

Commits

Matt Knepley  committed b5f81b6 Merge

Merge branch 'knepley/fix-fem-multifield' into next

* knepley/fix-fem-multifield:
DMPlex ex3: Use PetscFECreateDefault() and correct tensor tests
PetscFE: Fix PetscFECreateDefault() - The space should default to tensor on a non-simplex - The quadrature should be tensor Gauss on a non-simplex
DT: Added PetscDTGaussTensorQuadrature()
DT: Added citation for Gloub+Welsch method
PetscSection: Now global section can also be viewed

  • Participants
  • Parent commits cf2fe8f, 1ce47e0
  • Branches knepley/fix-quadrature-order, next-oct-2014

Comments (0)

Files changed (7)

File config/builder.py

View file
  • Ignore whitespace
                                                                  {'numProcs': 1, 'args': '-dim 3 -petscspace_order 2 -num_comp 3 -qorder 2 -porder 1'},
                                                                  {'numProcs': 1, 'args': '-dim 3 -petscspace_order 2 -num_comp 3 -qorder 2 -porder 2'},
                                                                  # 12-14 2D P_1 on a quadrilaterial
-                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -num_comp 2 -qorder 1', 'requires': 'Broken'},
-                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -num_comp 2 -qorder 1 -porder 1', 'requires': 'Broken'},
-                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -num_comp 2 -qorder 1 -porder 2', 'requires': 'Broken'},
+                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -petscspace_poly_tensor 0 -num_comp 2 -qorder 1', 'requires': 'Broken'},
+                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -petscspace_poly_tensor 0 -num_comp 2 -qorder 1 -porder 1', 'requires': 'Broken'},
+                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -petscspace_poly_tensor 0 -num_comp 2 -qorder 1 -porder 2', 'requires': 'Broken'},
                                                                  # 15-17 2D Q_1 on a quadrilaterial
-                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -petscspace_poly_tensor -num_comp 2 -qorder 1'},
-                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -petscspace_poly_tensor -num_comp 2 -qorder 1 -porder 1'},
-                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -petscspace_poly_tensor -num_comp 2 -qorder 1 -porder 2'}],
+                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -num_comp 2 -qorder 1'},
+                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -num_comp 2 -qorder 1 -porder 1'},
+                                                                 {'numProcs': 1, 'args': '-simplex 0 -petscspace_order 1 -num_comp 2 -qorder 1 -porder 2'}],
                         'src/dm/impls/plex/examples/tests/ex4': [# 2D Simplex 0-3
                                                                  {'numProcs': 1, 'args': '-dim 2 -cell_hybrid 0 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 1, 'args': '-dim 2 -cell_hybrid 0 -num_refinements 1 -dm_view ::ascii_info_detail'},

File include/petscdt.h

View file
  • Ignore whitespace
 PETSC_EXTERN PetscErrorCode PetscDTLegendreEval(PetscInt,const PetscReal*,PetscInt,const PetscInt*,PetscReal*,PetscReal*,PetscReal*);
 PETSC_EXTERN PetscErrorCode PetscDTGaussQuadrature(PetscInt,PetscReal,PetscReal,PetscReal*,PetscReal*);
 PETSC_EXTERN PetscErrorCode PetscDTReconstructPoly(PetscInt,PetscInt,const PetscReal*,PetscInt,const PetscReal*,PetscReal*);
+PETSC_EXTERN PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt,PetscInt,PetscReal,PetscReal,PetscQuadrature*);
 PETSC_EXTERN PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt,PetscInt,PetscReal,PetscReal,PetscQuadrature*);
 
 #endif

File src/dm/dt/interface/dt.c

View file
  • Ignore whitespace
 #include <petscdmplex.h>
 #include <petscdmshell.h>
 
+static PetscBool GaussCite       = PETSC_FALSE;
+const char       GaussCitation[] = "@article{GolubWelsch1969,\n"
+                                   "  author  = {Golub and Welsch},\n"
+                                   "  title   = {Calculation of Quadrature Rules},\n"
+                                   "  journal = {Math. Comp.},\n"
+                                   "  volume  = {23},\n"
+                                   "  number  = {106},\n"
+                                   "  pages   = {221--230},\n"
+                                   "  year    = {1969}\n}\n";
+
 #undef __FUNCT__
 #define __FUNCT__ "PetscQuadratureCreate"
 PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
   PetscBLASInt   N,LDZ,info;
 
   PetscFunctionBegin;
+  ierr = PetscCitationsRegister(GaussCitation, &GaussCite);CHKERRQ(ierr);
   /* Set up the Golub-Welsch system */
   for (i=0; i<npoints; i++) {
     x[i] = 0;                   /* diagonal is 0 */
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "PetscDTGaussTensorQuadrature"
+/*@
+  PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
+
+  Not Collective
+
+  Input Arguments:
++ dim     - The spatial dimension
+. npoints - number of points in one dimension
+. a       - left end of interval (often-1)
+- b       - right end of interval (often +1)
+
+  Output Argument:
+. q - A PetscQuadrature object
+
+  Level: intermediate
+
+.seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
+@*/
+PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
+{
+  PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k;
+  PetscReal     *x, *w, *xw, *ww;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscMalloc1(totpoints*dim,&x);CHKERRQ(ierr);
+  ierr = PetscMalloc1(totpoints,&w);CHKERRQ(ierr);
+  /* Set up the Golub-Welsch system */
+  switch (dim) {
+  case 0:
+    ierr = PetscFree(x);CHKERRQ(ierr);
+    ierr = PetscFree(w);CHKERRQ(ierr);
+    ierr = PetscMalloc1(1, &x);CHKERRQ(ierr);
+    ierr = PetscMalloc1(1, &w);CHKERRQ(ierr);
+    x[0] = 0.0;
+    w[0] = 1.0;
+    break;
+  case 1:
+    ierr = PetscDTGaussQuadrature(npoints, a, b, x, w);CHKERRQ(ierr);
+    break;
+  case 2:
+    ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr);
+    ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr);
+    for (i = 0; i < npoints; ++i) {
+      for (j = 0; j < npoints; ++j) {
+        x[(i*npoints+j)*dim+0] = xw[i];
+        x[(i*npoints+j)*dim+1] = xw[j];
+        w[i*npoints+j]         = ww[i] * ww[j];
+      }
+    }
+    ierr = PetscFree2(xw,ww);CHKERRQ(ierr);
+    break;
+  case 3:
+    ierr = PetscMalloc2(npoints,&xw,npoints,&ww);CHKERRQ(ierr);
+    ierr = PetscDTGaussQuadrature(npoints, a, b, xw, ww);CHKERRQ(ierr);
+    for (i = 0; i < npoints; ++i) {
+      for (j = 0; j < npoints; ++j) {
+        for (k = 0; k < npoints; ++k) {
+          x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
+          x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
+          x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
+          w[(i*npoints+j)*npoints+k]         = ww[i] * ww[j] * ww[k];
+        }
+      }
+    }
+    ierr = PetscFree2(xw,ww);CHKERRQ(ierr);
+    break;
+  default:
+    SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
+  }
+  ierr = PetscQuadratureCreate(PETSC_COMM_SELF, q);CHKERRQ(ierr);
+  ierr = PetscQuadratureSetData(*q, dim, totpoints, x, w);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "PetscDTFactorial_Internal"
 /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
   Not Collective
 
   Input Arguments:
-+ dim - The simplex dimension
-. order - The quadrature order
-. a - left end of interval (often-1)
-- b - right end of interval (often +1)
++ dim   - The simplex dimension
+. order - The number of points in one dimension
+. a     - left end of interval (often-1)
+- b     - right end of interval (often +1)
 
-  Output Arguments:
+  Output Argument:
 . q - A PetscQuadrature object
 
   Level: intermediate
   Karniadakis and Sherwin.
   FIAT
 
-.seealso: PetscDTGaussQuadrature()
+.seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
 @*/
 PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt order, PetscReal a, PetscReal b, PetscQuadrature *q)
 {

File src/dm/dt/interface/dtfe.c

View file
  • Ignore whitespace
   ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr);
   ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr);
   ierr = PetscSpacePolynomialSetNumVariables(P, dim);CHKERRQ(ierr);
+  ierr = PetscSpacePolynomialSetTensor(P, !isSimplex);CHKERRQ(ierr);
   ierr = PetscSpaceSetUp(P);CHKERRQ(ierr);
   ierr = PetscSpaceGetOrder(P, &order);CHKERRQ(ierr);
   /* Create dual space */
   ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr);
   ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr);
   /* Create quadrature (with specified order if given) */
-  ierr = PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);CHKERRQ(ierr);
+  if (isSimplex) {ierr = PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);CHKERRQ(ierr);}
+  else           {ierr = PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);CHKERRQ(ierr);}
   ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr);
   ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr);
   PetscFunctionReturn(0);

File src/dm/impls/plex/examples/tests/ex3.c

View file
  • Ignore whitespace
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "SetupElement"
-PetscErrorCode SetupElement(DM dm, AppCtx *user)
-{
-  PetscFE         fem;
-  DM              K;
-  PetscSpace      P;
-  PetscDualSpace  Q;
-  PetscReal      *B, *D;
-  PetscQuadrature q;
-  const PetscInt  dim   = user->dim;
-  PetscInt        order, numBasisFunc, numBasisComp;
-  PetscErrorCode  ierr;
-
-  PetscFunctionBegin;
-  /* Create space */
-  ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);CHKERRQ(ierr);
-  ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr);
-  ierr = PetscSpacePolynomialSetNumVariables(P, dim);CHKERRQ(ierr);
-  ierr = PetscSpaceSetUp(P);CHKERRQ(ierr);
-  ierr = PetscSpaceGetOrder(P, &order);CHKERRQ(ierr);
-  /* Create dual space */
-  ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);CHKERRQ(ierr);
-  ierr = PetscDualSpaceCreateReferenceCell(Q, dim, user->simplex, &K);CHKERRQ(ierr);
-  ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr);
-  ierr = DMDestroy(&K);CHKERRQ(ierr);
-  ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr);
-  ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr);
-  ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr);
-  /* Create quadrature */
-  ierr = PetscDTGaussJacobiQuadrature(dim, PetscMax(user->qorder, order), -1.0, 1.0, &q);CHKERRQ(ierr);
-  /* Create element */
-  ierr = PetscFECreate(PetscObjectComm((PetscObject) dm), &fem);CHKERRQ(ierr);
-  ierr = PetscFESetBasisSpace(fem, P);CHKERRQ(ierr);
-  ierr = PetscFESetDualSpace(fem, Q);CHKERRQ(ierr);
-  ierr = PetscFESetNumComponents(fem, user->numComponents);CHKERRQ(ierr);
-  ierr = PetscFESetQuadrature(fem, q);CHKERRQ(ierr);
-  ierr = PetscFESetFromOptions(fem);CHKERRQ(ierr);
-  ierr = PetscFESetUp(fem);CHKERRQ(ierr);
-  ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr);
-  ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr);
-  ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr);
-  ierr = PetscFEGetDimension(fem, &numBasisFunc);CHKERRQ(ierr);
-  ierr = PetscFEGetNumComponents(fem, &numBasisComp);CHKERRQ(ierr);
-  ierr = PetscFEGetDefaultTabulation(fem, &B, &D, NULL);CHKERRQ(ierr);
-  user->fe = fem;
-  user->fem.fe = &user->fe;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
 #define __FUNCT__ "SetupSection"
 PetscErrorCode SetupSection(DM dm, AppCtx *user)
 {
   ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
-  ierr = SetupElement(dm, &user);CHKERRQ(ierr);
+  ierr = PetscFECreateDefault(dm, user.dim, user.numComponents, user.simplex, NULL, user.qorder, &user.fe);CHKERRQ(ierr);
+  user.fem.fe = &user.fe;
   ierr = SetupSection(dm, &user);CHKERRQ(ierr);
   ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr);
   ierr = CheckFunctions(dm, user.porder, u, &user);CHKERRQ(ierr);

File src/dm/impls/plex/examples/tests/output/ex3_17.out

View file
  • Ignore whitespace
-Tests FAIL for order 2 at tolerance 1e-10 error 0.0174594
+Tests FAIL for order 2 at tolerance 1e-10 error 0.0277778
 Tests pass for order 2 derivatives at tolerance 1e-10

File src/vec/is/utils/vsectionis.c

View file
  • Ignore whitespace
     }
   }
   ierr = PetscFree2(neg,recv);CHKERRQ(ierr);
+  ierr = PetscSectionViewFromOptions(*gsection,NULL,"-section_global_view");CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }