Actual source code: ex1.c

      1: const char help[] = "Construct and set a Lagrange dual space from options, then view it to\n"
      2:                     "understand the effects of different parameters.";

      4: #include <petscfe.h>
      5: #include <petscdmplex.h>

      7: int main(int argc, char **argv)
      8: {
      9:   PetscInt       dim;
     10:   PetscBool      tensorCell;
     11:   DM             K;
     12:   PetscDualSpace dsp;

     14:   PetscFunctionBeginUser;
     15:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     16:   dim        = 2;
     17:   tensorCell = PETSC_FALSE;
     18:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PETSCDUALSPACELAGRANGE test", "none");
     19:   PetscCall(PetscOptionsRangeInt("-dim", "The spatial dimension", "ex1.c", dim, &dim, NULL, 0, 3));
     20:   PetscCall(PetscOptionsBool("-tensor", "Whether the cell is a tensor product cell or a simplex", "ex1.c", tensorCell, &tensorCell, NULL));
     21:   PetscOptionsEnd();

     23:   PetscCall(PetscDualSpaceCreate(PETSC_COMM_WORLD, &dsp));
     24:   PetscCall(PetscObjectSetName((PetscObject)dsp, "Lagrange dual space"));
     25:   PetscCall(PetscDualSpaceSetType(dsp, PETSCDUALSPACELAGRANGE));
     26:   /* While Lagrange nodes don't require the existence of a reference cell to
     27:    * be refined, when we construct finite element dual spaces we have to be
     28:    * careful about what kind of continuity is maintained when cells are glued
     29:    * together to make a mesh.  The PetscDualSpace object is responsible for
     30:    * conveying continuity requirements to a finite element assembly routines,
     31:    * so a PetscDualSpace needs a reference element: a single element mesh,
     32:    * whose boundary points are the interstitial points in a mesh */
     33:   PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_WORLD, DMPolytopeTypeSimpleShape(dim, (PetscBool)!tensorCell), &K));
     34:   PetscCall(PetscDualSpaceSetDM(dsp, K));
     35:   /* This gives us the opportunity to change the parameters of the dual space
     36:    * from the command line, as we do in the tests below.  When
     37:    * PetscDualSpaceSetFromOptions() is called, it also enables other optional
     38:    * behavior (see the next step) */
     39:   PetscCall(PetscDualSpaceSetFromOptions(dsp));
     40:   /* This step parses the parameters of the dual space into
     41:    * sets of functionals that are assigned to each of the mesh points in K.
     42:    *
     43:    * The functionals can be accessed individually by
     44:    * PetscDualSpaceGetFunctional(), or more efficiently all at once by
     45:    * PetscDualSpaceGetAllData(), which returns a set of quadrature points
     46:    * at which to evaluate a function, and a matrix that takes those
     47:    * evaluations and turns them into the evaluation of the dual space's
     48:    * functionals on the function.
     49:    *
     50:    * (TODO: tutorial for PetscDualSpaceGetAllData() and
     51:    * PetscDualSpaceGetInteriorData().)
     52:    *
     53:    * Because we called PetscDualSpaceSetFromOptions(), we have the opportunity
     54:    * to inspect the results of PetscDualSpaceSetUp() from the command line
     55:    * with "-petscdualspace_view", followed by an optional description of how
     56:    * we would like to see the dual space (see examples in the tests below).
     57:    * */
     58:   PetscCall(PetscDualSpaceSetUp(dsp));
     59:   PetscCall(DMDestroy(&K));
     60:   PetscCall(PetscDualSpaceDestroy(&dsp));
     61:   PetscCall(PetscFinalize());
     62:   return 0;
     63: }

     65: /*TEST

     67:   # quadratic nodes on the triangle
     68:   test:
     69:     suffix: 0
     70:     filter: sed -E "s/\(\+0, \+0\)/(+0., +0.)/g"
     71:     args: -dim 2 -tensor 0 -petscdualspace_order 2 -petscdualspace_view ascii::ascii_info_detail

     73:   # linear nodes on the quadrilateral
     74:   test:
     75:     suffix: 1
     76:     args: -dim 2 -tensor 1 -petscdualspace_order 1 -petscdualspace_lagrange_tensor 1 -petscdualspace_view ascii::ascii_info_detail

     78:   # lowest order Raviart-Thomas / Nedelec edge nodes on the hexahedron
     79:   test:
     80:     suffix: 2
     81:     args: -dim 3 -tensor 1 -petscdualspace_order 1 -petscdualspace_components 3 -petscdualspace_form_degree 1 -petscdualspace_lagrange_trimmed 1 -petscdualspace_lagrange_tensor 1 -petscdualspace_view ascii::ascii_info_detail

     83:   # first order Nedelec second type face nodes on the tetrahedron
     84:   test:
     85:     suffix: 3
     86:     args: -dim 3 -tensor 0 -petscdualspace_order 1 -petscdualspace_components 3 -petscdualspace_form_degree -2 -petscdualspace_view ascii::ascii_info_detail

     88:   ## Comparing different node types

     90:   test:
     91:     suffix: 4
     92:     args: -dim 2 -tensor 0 -petscdualspace_order 3 -petscdualspace_lagrange_continuity 0 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 0 -petscdualspace_view ascii::ascii_info_detail

     94:   test:
     95:     suffix: 5
     96:     args: -dim 2 -tensor 0 -petscdualspace_order 3 -petscdualspace_lagrange_continuity 0 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1 -petscdualspace_view ascii::ascii_info_detail

     98:   test:
     99:     suffix: 6
    100:     args: -dim 2 -tensor 0 -petscdualspace_order 3 -petscdualspace_lagrange_continuity 0 -petscdualspace_lagrange_node_type gaussjacobi -petscdualspace_lagrange_node_endpoints 0 -petscdualspace_view ascii::ascii_info_detail

    102:   test:
    103:     suffix: 7
    104:     args: -dim 2 -tensor 0 -petscdualspace_order 3 -petscdualspace_lagrange_continuity 0 -petscdualspace_lagrange_node_type gaussjacobi -petscdualspace_lagrange_node_endpoints 1 -petscdualspace_view ascii::ascii_info_detail

    106: TEST*/