xref: /petsc/src/dm/impls/plex/tests/ex33.c (revision 39290ff683c9b991eb3b8ed2ff3d692480fddb9f)
18b438e21SMatthew G. Knepley static char help[] = "Tests for high order geometry\n\n";
28b438e21SMatthew G. Knepley 
38b438e21SMatthew G. Knepley #include <petscdmplex.h>
48b438e21SMatthew G. Knepley #include <petscds.h>
58b438e21SMatthew G. Knepley 
68b438e21SMatthew G. Knepley typedef enum {TRANSFORM_NONE, TRANSFORM_SHEAR, TRANSFORM_ANNULUS, TRANSFORM_SHELL} Transform;
78b438e21SMatthew G. Knepley const char * const TransformTypes[] = {"none", "shear", "annulus", "shell", "Mesh Transform", "TRANSFORM_", NULL};
88b438e21SMatthew G. Knepley 
98b438e21SMatthew G. Knepley typedef struct {
108b438e21SMatthew G. Knepley   char        filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
118b438e21SMatthew G. Knepley   Transform   meshTransform;      /* Transform for initial box mesh */
128b438e21SMatthew G. Knepley   PetscReal   *transformDataReal; /* Parameters for mesh transform */
138b438e21SMatthew G. Knepley   PetscScalar *transformData;     /* Parameters for mesh transform */
148b438e21SMatthew G. Knepley   PetscReal   volume;             /* Analytical volume of the mesh */
158b438e21SMatthew G. Knepley   PetscReal   tol;                /* Tolerance for volume check */
168b438e21SMatthew G. Knepley } AppCtx;
178b438e21SMatthew G. Knepley 
188b438e21SMatthew G. Knepley PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
198b438e21SMatthew G. Knepley {
208b438e21SMatthew G. Knepley   PetscInt       n = 0, i;
218b438e21SMatthew G. Knepley   PetscErrorCode ierr;
228b438e21SMatthew G. Knepley 
238b438e21SMatthew G. Knepley   PetscFunctionBegin;
248b438e21SMatthew G. Knepley   options->filename[0]       = '\0';
258b438e21SMatthew G. Knepley   options->meshTransform     = TRANSFORM_NONE;
268b438e21SMatthew G. Knepley   options->transformDataReal = NULL;
278b438e21SMatthew G. Knepley   options->transformData     = NULL;
288b438e21SMatthew G. Knepley   options->volume            = -1.0;
298b438e21SMatthew G. Knepley   options->tol               = PETSC_SMALL;
308b438e21SMatthew G. Knepley 
318b438e21SMatthew G. Knepley   ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr);
328b438e21SMatthew G. Knepley   ierr = PetscOptionsString("-filename", "The mesh file", "ex33.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr);
338b438e21SMatthew G. Knepley   ierr = PetscOptionsEnum("-mesh_transform", "Method to transform initial box mesh <none, shear, annulus, shell>", "ex33.c", TransformTypes, (PetscEnum) options->meshTransform, (PetscEnum *) &options->meshTransform, NULL);CHKERRQ(ierr);
348b438e21SMatthew G. Knepley   switch (options->meshTransform) {
358b438e21SMatthew G. Knepley     case TRANSFORM_NONE: break;
368b438e21SMatthew G. Knepley     case TRANSFORM_SHEAR:
37*39290ff6SMatthew G. Knepley       n = 2;
388b438e21SMatthew G. Knepley       ierr = PetscMalloc1(n, &options->transformDataReal);CHKERRQ(ierr);
398b438e21SMatthew G. Knepley       for (i = 0; i < n; ++i) options->transformDataReal[i] = 1.0;
408b438e21SMatthew G. Knepley       ierr = PetscOptionsRealArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformDataReal, &n, NULL);CHKERRQ(ierr);
418b438e21SMatthew G. Knepley       break;
428b438e21SMatthew G. Knepley     case TRANSFORM_ANNULUS:
438b438e21SMatthew G. Knepley       n = 2;
448b438e21SMatthew G. Knepley       ierr = PetscMalloc1(n, &options->transformData);CHKERRQ(ierr);
458b438e21SMatthew G. Knepley       options->transformData[0] = 1.0;
468b438e21SMatthew G. Knepley       options->transformData[1] = 2.0;
478b438e21SMatthew G. Knepley       ierr = PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL);CHKERRQ(ierr);
488b438e21SMatthew G. Knepley       break;
498b438e21SMatthew G. Knepley     case TRANSFORM_SHELL:
508b438e21SMatthew G. Knepley       n = 2;
518b438e21SMatthew G. Knepley       ierr = PetscMalloc1(n, &options->transformData);CHKERRQ(ierr);
528b438e21SMatthew G. Knepley       options->transformData[0] = 1.0;
538b438e21SMatthew G. Knepley       options->transformData[1] = 2.0;
548b438e21SMatthew G. Knepley       ierr = PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL);CHKERRQ(ierr);
558b438e21SMatthew G. Knepley       break;
568b438e21SMatthew G. Knepley     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %D", options->meshTransform);
578b438e21SMatthew G. Knepley   }
588b438e21SMatthew G. Knepley   ierr = PetscOptionsReal("-volume", "The analytical volume of the mesh", "ex33.c", options->volume, &options->volume, NULL);CHKERRQ(ierr);
598b438e21SMatthew G. Knepley   ierr = PetscOptionsReal("-tol", "The tolerance for the volume check", "ex33.c", options->tol, &options->tol, NULL);CHKERRQ(ierr);
608b438e21SMatthew G. Knepley   ierr = PetscOptionsEnd();
618b438e21SMatthew G. Knepley   PetscFunctionReturn(0);
628b438e21SMatthew G. Knepley }
638b438e21SMatthew G. Knepley 
648b438e21SMatthew G. Knepley static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
658b438e21SMatthew G. Knepley                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
668b438e21SMatthew G. Knepley                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
678b438e21SMatthew G. Knepley                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
688b438e21SMatthew G. Knepley {
698b438e21SMatthew G. Knepley   const PetscInt Nc = uOff[1] - uOff[0];
708b438e21SMatthew G. Knepley   PetscInt       c;
718b438e21SMatthew G. Knepley 
728b438e21SMatthew G. Knepley   for (c = 0; c < Nc; ++c) f0[c] = u[c];
738b438e21SMatthew G. Knepley }
748b438e21SMatthew G. Knepley 
758b438e21SMatthew G. Knepley /*
768b438e21SMatthew G. Knepley   We would like to map the unit square to a quarter of the annulus between circles of radius 1 and 2. We start by mapping the straight sections, which
778b438e21SMatthew G. Knepley   will correspond to the top and bottom of our square. So
788b438e21SMatthew G. Knepley 
798b438e21SMatthew G. Knepley     (0,0)--(1,0)  ==>  (1,0)--(2,0)      Just a shift of (1,0)
808b438e21SMatthew G. Knepley     (0,1)--(1,1)  ==>  (0,1)--(0,2)      Switch x and y
818b438e21SMatthew G. Knepley 
828b438e21SMatthew G. Knepley   So it looks like we want to map each layer in y to a ray, so x is the radius and y is the angle:
838b438e21SMatthew G. Knepley 
848b438e21SMatthew G. Knepley     (x, y)  ==>  (x+1, \pi/2 y)                           in (r', \theta') space
858b438e21SMatthew G. Knepley             ==>  ((x+1) cos(\pi/2 y), (x+1) sin(\pi/2 y)) in (x', y') space
868b438e21SMatthew G. Knepley */
878b438e21SMatthew G. Knepley static void f0_annulus(PetscInt dim, PetscInt Nf, PetscInt NfAux,
888b438e21SMatthew G. Knepley                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
898b438e21SMatthew G. Knepley                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
908b438e21SMatthew G. Knepley                        PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
918b438e21SMatthew G. Knepley {
928b438e21SMatthew G. Knepley   const PetscReal ri = PetscRealPart(constants[0]);
938b438e21SMatthew G. Knepley   const PetscReal ro = PetscRealPart(constants[1]);
948b438e21SMatthew G. Knepley 
958b438e21SMatthew G. Knepley   xp[0] = (x[0] * (ro-ri) + ri) * PetscCosReal(0.5*PETSC_PI*x[1]);
968b438e21SMatthew G. Knepley   xp[1] = (x[0] * (ro-ri) + ri) * PetscSinReal(0.5*PETSC_PI*x[1]);
978b438e21SMatthew G. Knepley }
988b438e21SMatthew G. Knepley 
998b438e21SMatthew G. Knepley /*
1008b438e21SMatthew G. Knepley   We would like to map the unit cube to a hemisphere of the spherical shell between balls of radius 1 and 2. We want to map the bottom surface onto the
1018b438e21SMatthew G. Knepley   lower hemisphere and the upper surface onto the top, letting z be the radius.
1028b438e21SMatthew G. Knepley 
1038b438e21SMatthew G. Knepley     (x, y)  ==>  ((z+3)/2, \pi/2 (|x| or |y|), arctan y/x)                                                  in (r', \theta', \phi') space
1048b438e21SMatthew G. Knepley             ==>  ((z+3)/2 \cos(\theta') cos(\phi'), (z+3)/2 \cos(\theta') sin(\phi'), (z+3)/2 sin(\theta')) in (x', y', z') space
1058b438e21SMatthew G. Knepley */
1068b438e21SMatthew G. Knepley static void f0_shell(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1078b438e21SMatthew G. Knepley                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1088b438e21SMatthew G. Knepley                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1098b438e21SMatthew G. Knepley                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
1108b438e21SMatthew G. Knepley {
1118b438e21SMatthew G. Knepley   const PetscReal pi4    = PETSC_PI/4.0;
1128b438e21SMatthew G. Knepley   const PetscReal ri     = PetscRealPart(constants[0]);
1138b438e21SMatthew G. Knepley   const PetscReal ro     = PetscRealPart(constants[1]);
1148b438e21SMatthew G. Knepley   const PetscReal rp     = (x[2]+1) * 0.5*(ro-ri) + ri;
1158b438e21SMatthew G. Knepley   const PetscReal phip   = PetscAtan2Real(x[1], x[0]);
1168b438e21SMatthew G. Knepley   const PetscReal thetap = 0.5*PETSC_PI * (1.0 - ((((phip <= pi4) && (phip >= -pi4)) || ((phip >= 3.0*pi4) || (phip <= -3.0*pi4))) ? PetscAbsReal(x[0]) : PetscAbsReal(x[1])));
1178b438e21SMatthew G. Knepley 
1188b438e21SMatthew G. Knepley   xp[0] = rp * PetscCosReal(thetap) * PetscCosReal(phip);
1198b438e21SMatthew G. Knepley   xp[1] = rp * PetscCosReal(thetap) * PetscSinReal(phip);
1208b438e21SMatthew G. Knepley   xp[2] = rp * PetscSinReal(thetap);
1218b438e21SMatthew G. Knepley }
1228b438e21SMatthew G. Knepley 
123*39290ff6SMatthew G. Knepley static PetscErrorCode DMCreateCoordinateDisc(DM dm)
1248b438e21SMatthew G. Knepley {
1258b438e21SMatthew G. Knepley   DM             cdm;
1268b438e21SMatthew G. Knepley   PetscFE        fe;
127*39290ff6SMatthew G. Knepley   DMPolytopeType ct;
128*39290ff6SMatthew G. Knepley   PetscInt       dim, dE, cStart;
129*39290ff6SMatthew G. Knepley   PetscBool      simplex;
1308b438e21SMatthew G. Knepley   PetscErrorCode ierr;
1318b438e21SMatthew G. Knepley 
1328b438e21SMatthew G. Knepley   PetscFunctionBegin;
1338b438e21SMatthew G. Knepley   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
1348b438e21SMatthew G. Knepley   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1358b438e21SMatthew G. Knepley   ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr);
136*39290ff6SMatthew G. Knepley   ierr = DMPlexGetHeightStratum(cdm, 0, &cStart, NULL);CHKERRQ(ierr);
137*39290ff6SMatthew G. Knepley   ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr);
138*39290ff6SMatthew G. Knepley   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
139*39290ff6SMatthew G. Knepley   ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, dE, simplex, "geom_", -1, &fe);CHKERRQ(ierr);
140*39290ff6SMatthew G. Knepley   ierr = DMProjectCoordinates(dm, fe);CHKERRQ(ierr);
1418b438e21SMatthew G. Knepley   ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
1428b438e21SMatthew G. Knepley   PetscFunctionReturn(0);
1438b438e21SMatthew G. Knepley }
1448b438e21SMatthew G. Knepley 
1458b438e21SMatthew G. Knepley PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
1468b438e21SMatthew G. Knepley {
1478b438e21SMatthew G. Knepley   const char    *filename = ctx->filename;
1488b438e21SMatthew G. Knepley   DM             cdm;
1498b438e21SMatthew G. Knepley   PetscDS        cds;
1508b438e21SMatthew G. Knepley   size_t         len;
1518b438e21SMatthew G. Knepley   PetscMPIInt    rank;
1528b438e21SMatthew G. Knepley   PetscErrorCode ierr;
1538b438e21SMatthew G. Knepley 
1548b438e21SMatthew G. Knepley   PetscFunctionBegin;
1558b438e21SMatthew G. Knepley   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1568b438e21SMatthew G. Knepley   ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
157*39290ff6SMatthew G. Knepley   if (len) {ierr = DMPlexCreateFromFile(comm, filename, PETSC_TRUE, dm);CHKERRQ(ierr);}
158*39290ff6SMatthew G. Knepley   else     {ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);}
1598b438e21SMatthew G. Knepley   ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr);
1608b438e21SMatthew G. Knepley   ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
1618b438e21SMatthew G. Knepley   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
162*39290ff6SMatthew G. Knepley   if (!len) {ierr = DMCreateCoordinateDisc(*dm);CHKERRQ(ierr);}
1638b438e21SMatthew G. Knepley   switch (ctx->meshTransform) {
1648b438e21SMatthew G. Knepley     case TRANSFORM_NONE:
1658b438e21SMatthew G. Knepley       ierr = DMPlexRemapGeometry(*dm, 0.0, identity);CHKERRQ(ierr);
1668b438e21SMatthew G. Knepley       break;
1678b438e21SMatthew G. Knepley     case TRANSFORM_SHEAR:
1688b438e21SMatthew G. Knepley       ierr = DMPlexShearGeometry(*dm, DM_X, ctx->transformDataReal);CHKERRQ(ierr);
1698b438e21SMatthew G. Knepley       break;
1708b438e21SMatthew G. Knepley     case TRANSFORM_ANNULUS:
1718b438e21SMatthew G. Knepley       ierr = DMGetCoordinateDM(*dm, &cdm);CHKERRQ(ierr);
1728b438e21SMatthew G. Knepley       ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr);
1738b438e21SMatthew G. Knepley       ierr = PetscDSSetConstants(cds, 2, ctx->transformData);CHKERRQ(ierr);
1748b438e21SMatthew G. Knepley       ierr = DMPlexRemapGeometry(*dm, 0.0, f0_annulus);CHKERRQ(ierr);
1758b438e21SMatthew G. Knepley       break;
1768b438e21SMatthew G. Knepley     case TRANSFORM_SHELL:
1778b438e21SMatthew G. Knepley       ierr = DMGetCoordinateDM(*dm, &cdm);CHKERRQ(ierr);
1788b438e21SMatthew G. Knepley       ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr);
1798b438e21SMatthew G. Knepley       ierr = PetscDSSetConstants(cds, 2, ctx->transformData);CHKERRQ(ierr);
1808b438e21SMatthew G. Knepley       ierr = DMPlexRemapGeometry(*dm, 0.0, f0_shell);CHKERRQ(ierr);
1818b438e21SMatthew G. Knepley       break;
1828b438e21SMatthew G. Knepley     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %D", ctx->meshTransform);
1838b438e21SMatthew G. Knepley   }
1848b438e21SMatthew G. Knepley   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
1858b438e21SMatthew G. Knepley   PetscFunctionReturn(0);
1868b438e21SMatthew G. Knepley }
1878b438e21SMatthew G. Knepley 
1888b438e21SMatthew G. Knepley static void volume(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1898b438e21SMatthew G. Knepley                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1908b438e21SMatthew G. Knepley                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1918b438e21SMatthew G. Knepley                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar vol[])
1928b438e21SMatthew G. Knepley {
1938b438e21SMatthew G. Knepley   vol[0] = 1.;
1948b438e21SMatthew G. Knepley }
1958b438e21SMatthew G. Knepley 
1968b438e21SMatthew G. Knepley static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx)
1978b438e21SMatthew G. Knepley {
1988b438e21SMatthew G. Knepley   PetscDS        ds;
1998b438e21SMatthew G. Knepley   PetscFE        fe;
200*39290ff6SMatthew G. Knepley   DMPolytopeType ct;
201*39290ff6SMatthew G. Knepley   PetscInt       dim, cStart;
202*39290ff6SMatthew G. Knepley   PetscBool      simplex;
2038b438e21SMatthew G. Knepley   PetscErrorCode ierr;
2048b438e21SMatthew G. Knepley 
2058b438e21SMatthew G. Knepley   PetscFunctionBeginUser;
206*39290ff6SMatthew G. Knepley   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
207*39290ff6SMatthew G. Knepley   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);
208*39290ff6SMatthew G. Knepley   ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr);
209*39290ff6SMatthew G. Knepley   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
210*39290ff6SMatthew G. Knepley   ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fe);CHKERRQ(ierr);
2118b438e21SMatthew G. Knepley   ierr = PetscFESetName(fe, "scalar");CHKERRQ(ierr);
2128b438e21SMatthew G. Knepley   ierr = DMAddField(dm, NULL, (PetscObject) fe);
2138b438e21SMatthew G. Knepley   ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
2148b438e21SMatthew G. Knepley   ierr = DMCreateDS(dm);CHKERRQ(ierr);
2158b438e21SMatthew G. Knepley   ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
2168b438e21SMatthew G. Knepley   ierr = PetscDSSetObjective(ds, 0, volume);CHKERRQ(ierr);
2178b438e21SMatthew G. Knepley   PetscFunctionReturn(0);
2188b438e21SMatthew G. Knepley }
2198b438e21SMatthew G. Knepley 
2208b438e21SMatthew G. Knepley static PetscErrorCode CheckVolume(DM dm, AppCtx *ctx)
2218b438e21SMatthew G. Knepley {
2228b438e21SMatthew G. Knepley   Vec            u;
2238b438e21SMatthew G. Knepley   PetscScalar    result;
2248b438e21SMatthew G. Knepley   PetscReal      vol, tol = ctx->tol;
2258b438e21SMatthew G. Knepley   PetscErrorCode ierr;
2268b438e21SMatthew G. Knepley 
2278b438e21SMatthew G. Knepley   PetscFunctionBeginUser;
2288b438e21SMatthew G. Knepley   ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr);
2298b438e21SMatthew G. Knepley   ierr = DMPlexComputeIntegralFEM(dm, u, &result, ctx);CHKERRQ(ierr);
2308b438e21SMatthew G. Knepley   vol  = PetscRealPart(result);
2318b438e21SMatthew G. Knepley   ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr);
2328b438e21SMatthew G. Knepley   ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Volume: %g\n", (double) vol);CHKERRQ(ierr);
2338b438e21SMatthew G. Knepley   if (ctx->volume > 0.0 && PetscAbsReal(ctx->volume - vol) > tol) {
2348b438e21SMatthew G. Knepley     SETERRQ4(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Calculated volume %g != %g actual volume (error %g > %g tol)", (double) vol, (double) ctx->volume, (double) PetscAbsReal(ctx->volume - vol), (double) tol);
2358b438e21SMatthew G. Knepley   }
2368b438e21SMatthew G. Knepley   PetscFunctionReturn(0);
2378b438e21SMatthew G. Knepley }
2388b438e21SMatthew G. Knepley 
2398b438e21SMatthew G. Knepley int main(int argc, char **argv)
2408b438e21SMatthew G. Knepley {
241*39290ff6SMatthew G. Knepley   DM             dm;
2428b438e21SMatthew G. Knepley   AppCtx         user;
2438b438e21SMatthew G. Knepley   PetscErrorCode ierr;
2448b438e21SMatthew G. Knepley 
2458b438e21SMatthew G. Knepley   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
2468b438e21SMatthew G. Knepley   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
247*39290ff6SMatthew G. Knepley   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
248*39290ff6SMatthew G. Knepley   ierr = CreateDiscretization(dm, &user);CHKERRQ(ierr);
249*39290ff6SMatthew G. Knepley   ierr = CheckVolume(dm, &user);CHKERRQ(ierr);
250*39290ff6SMatthew G. Knepley   ierr = DMDestroy(&dm);CHKERRQ(ierr);
2518b438e21SMatthew G. Knepley   ierr = PetscFree(user.transformDataReal);CHKERRQ(ierr);
2528b438e21SMatthew G. Knepley   ierr = PetscFree(user.transformData);CHKERRQ(ierr);
2538b438e21SMatthew G. Knepley   ierr = PetscFinalize();
2548b438e21SMatthew G. Knepley   return ierr;
2558b438e21SMatthew G. Knepley }
2568b438e21SMatthew G. Knepley 
2578b438e21SMatthew G. Knepley /*TEST
2588b438e21SMatthew G. Knepley 
2598b438e21SMatthew G. Knepley   test:
2608b438e21SMatthew G. Knepley     suffix: square_0
261*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -geom_petscspace_degree 1 -volume 4.
2628b438e21SMatthew G. Knepley 
2638b438e21SMatthew G. Knepley   test:
2648b438e21SMatthew G. Knepley     suffix: square_1
265*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -geom_petscspace_degree 2 -volume 4.
2668b438e21SMatthew G. Knepley 
2678b438e21SMatthew G. Knepley   test:
2688b438e21SMatthew G. Knepley     suffix: square_2
269*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -dm_refine 1 -geom_petscspace_degree 1 -volume 4.
2708b438e21SMatthew G. Knepley 
2718b438e21SMatthew G. Knepley   test:
2728b438e21SMatthew G. Knepley     suffix: square_3
273*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -dm_refine 1 -geom_petscspace_degree 2 -volume 4.
2748b438e21SMatthew G. Knepley 
2758b438e21SMatthew G. Knepley   test:
2768b438e21SMatthew G. Knepley     suffix: cube_0
277*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -geom_petscspace_degree 1 -volume 8.
2788b438e21SMatthew G. Knepley 
2798b438e21SMatthew G. Knepley   test:
2808b438e21SMatthew G. Knepley     suffix: cube_1
281*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -geom_petscspace_degree 2 -volume 8.
2828b438e21SMatthew G. Knepley 
2838b438e21SMatthew G. Knepley   test:
2848b438e21SMatthew G. Knepley     suffix: cube_2
285*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 1 -geom_petscspace_degree 1 -volume 8.
2868b438e21SMatthew G. Knepley 
2878b438e21SMatthew G. Knepley   test:
2888b438e21SMatthew G. Knepley     suffix: cube_3
289*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 1 -geom_petscspace_degree 2 -volume 8.
2908b438e21SMatthew G. Knepley 
2918b438e21SMatthew G. Knepley   test:
2928b438e21SMatthew G. Knepley     suffix: shear_0
293*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -geom_petscspace_degree 1 -mesh_transform shear -transform_data 3.0 -volume 4.
2948b438e21SMatthew G. Knepley 
2958b438e21SMatthew G. Knepley   test:
2968b438e21SMatthew G. Knepley     suffix: shear_1
297*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -geom_petscspace_degree 2 -mesh_transform shear -transform_data 3.0 -volume 4.
2988b438e21SMatthew G. Knepley 
2998b438e21SMatthew G. Knepley   test:
3008b438e21SMatthew G. Knepley     suffix: shear_2
301*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -dm_refine 1 -geom_petscspace_degree 1 -mesh_transform shear -transform_data 3.0 -volume 4.
3028b438e21SMatthew G. Knepley 
3038b438e21SMatthew G. Knepley   test:
3048b438e21SMatthew G. Knepley     suffix: shear_3
305*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -dm_refine 1 -geom_petscspace_degree 2 -mesh_transform shear -transform_data 3.0 -volume 4.
3068b438e21SMatthew G. Knepley 
3078b438e21SMatthew G. Knepley   test:
3088b438e21SMatthew G. Knepley     suffix: shear_4
309*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -geom_petscspace_degree 1 -mesh_transform shear -transform_data 3.0 -volume 8.
3108b438e21SMatthew G. Knepley 
3118b438e21SMatthew G. Knepley   test:
3128b438e21SMatthew G. Knepley     suffix: shear_5
313*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -geom_petscspace_degree 2 -mesh_transform shear -transform_data 3.0 -volume 8.
3148b438e21SMatthew G. Knepley 
3158b438e21SMatthew G. Knepley   test:
3168b438e21SMatthew G. Knepley     suffix: shear_6
317*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 1 -geom_petscspace_degree 1 -mesh_transform shear -transform_data 3.0,4.0 -volume 8.
3188b438e21SMatthew G. Knepley 
3198b438e21SMatthew G. Knepley   test:
3208b438e21SMatthew G. Knepley     suffix: shear_7
321*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 1 -geom_petscspace_degree 2 -mesh_transform shear -transform_data 3.0,4.0 -volume 8.
3228b438e21SMatthew G. Knepley 
3238b438e21SMatthew G. Knepley   test:
3248b438e21SMatthew G. Knepley     # Area: (a+b)/2 h = 3/\sqrt{2} (sqrt{2} - 1/\sqrt{2}) = 3/2
3258b438e21SMatthew G. Knepley     suffix: annulus_0
3268b438e21SMatthew G. Knepley     requires: double
327*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -geom_petscspace_degree 1 -mesh_transform annulus -volume 1.5
3288b438e21SMatthew G. Knepley 
3298b438e21SMatthew G. Knepley   test:
3308b438e21SMatthew G. Knepley     # Area: 3/4 \pi = 2.3562
3318b438e21SMatthew G. Knepley     suffix: annulus_1
3328b438e21SMatthew G. Knepley     requires: double
333*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_refine 3 -geom_petscspace_degree 1 -mesh_transform annulus -volume 2.35619449019235 -tol .016
3348b438e21SMatthew G. Knepley 
3358b438e21SMatthew G. Knepley   test:
3368b438e21SMatthew G. Knepley     # Area: 3/4 \pi = 2.3562
3378b438e21SMatthew G. Knepley     suffix: annulus_2
3388b438e21SMatthew G. Knepley     requires: double
339*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_refine 3 -geom_petscspace_degree 2 -mesh_transform annulus -volume 2.35619449019235 -tol .0038
3408b438e21SMatthew G. Knepley 
3418b438e21SMatthew G. Knepley   test:
3428b438e21SMatthew G. Knepley     # Area: 3/4 \pi = 2.3562
3438b438e21SMatthew G. Knepley     suffix: annulus_3
3448b438e21SMatthew G. Knepley     requires: double
345*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_refine 3 -geom_petscspace_degree 3 -mesh_transform annulus -volume 2.35619449019235 -tol 2.2e-6
3468b438e21SMatthew G. Knepley 
3478b438e21SMatthew G. Knepley   test:
3488b438e21SMatthew G. Knepley     # Area: 3/4 \pi = 2.3562
3498b438e21SMatthew G. Knepley     suffix: annulus_4
3508b438e21SMatthew G. Knepley     requires: double
351*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_refine 2 -geom_petscspace_degree 2 -petscfe_default_quadrature_order 2 -mesh_transform annulus -volume 2.35619449019235 -tol .00012
3528b438e21SMatthew G. Knepley 
3538b438e21SMatthew G. Knepley   test:
3548b438e21SMatthew G. Knepley     # Area: 3/4 \pi = 2.3562
3558b438e21SMatthew G. Knepley     suffix: annulus_5
3568b438e21SMatthew G. Knepley     requires: double
357*39290ff6SMatthew G. Knepley     args: -dm_plex_box_faces 1,1 -dm_refine 2 -geom_petscspace_degree 3 -petscfe_default_quadrature_order 3 -mesh_transform annulus -volume 2.35619449019235 -tol 1.2e-7
3588b438e21SMatthew G. Knepley 
3598b438e21SMatthew G. Knepley   test:
3608b438e21SMatthew G. Knepley     suffix: shell_0
3618b438e21SMatthew G. Knepley     requires: double
362*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 1 -geom_petscspace_degree 1 -petscfe_default_quadrature_order 1 -mesh_transform shell -volume 5.633164922 -tol 1.0e-7
3638b438e21SMatthew G. Knepley 
3648b438e21SMatthew G. Knepley   test:
3658b438e21SMatthew G. Knepley     # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
3668b438e21SMatthew G. Knepley     suffix: shell_1
3678b438e21SMatthew G. Knepley     requires: double
368*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 2 -geom_petscspace_degree 1 -petscfe_default_quadrature_order 1 -mesh_transform shell -volume 14.66076571675238 -tol 3.1
3698b438e21SMatthew G. Knepley 
3708b438e21SMatthew G. Knepley   test:
3718b438e21SMatthew G. Knepley     # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
3728b438e21SMatthew G. Knepley     suffix: shell_2
3738b438e21SMatthew G. Knepley     requires: double
374*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 2 -geom_petscspace_degree 2 -petscfe_default_quadrature_order 2 -mesh_transform shell -volume 14.66076571675238 -tol .1
3758b438e21SMatthew G. Knepley 
3768b438e21SMatthew G. Knepley   test:
3778b438e21SMatthew G. Knepley     # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
3788b438e21SMatthew G. Knepley     suffix: shell_3
3798b438e21SMatthew G. Knepley     requires: double
380*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 2 -geom_petscspace_degree 3 -petscfe_default_quadrature_order 3 -mesh_transform shell -volume 14.66076571675238 -tol .02
3818b438e21SMatthew G. Knepley 
3828b438e21SMatthew G. Knepley   test:
3838b438e21SMatthew G. Knepley     # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
3848b438e21SMatthew G. Knepley     suffix: shell_4
3858b438e21SMatthew G. Knepley     requires: double
386*39290ff6SMatthew G. Knepley     args: -dm_plex_box_dim 3 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -dm_refine 2 -geom_petscspace_degree 4 -petscfe_default_quadrature_order 4 -mesh_transform shell -volume 14.66076571675238 -tol .006
387*39290ff6SMatthew G. Knepley 
388*39290ff6SMatthew G. Knepley   test:
389*39290ff6SMatthew G. Knepley     # Volume: 1.0
390*39290ff6SMatthew G. Knepley     suffix: gmsh_q2
391*39290ff6SMatthew G. Knepley     requires: double
392*39290ff6SMatthew G. Knepley     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q2.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
393*39290ff6SMatthew G. Knepley 
394*39290ff6SMatthew G. Knepley   test:
395*39290ff6SMatthew G. Knepley     # Volume: 1.0
396*39290ff6SMatthew G. Knepley     suffix: gmsh_q3
397*39290ff6SMatthew G. Knepley     requires: double
398*39290ff6SMatthew G. Knepley     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q3.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
3998b438e21SMatthew G. Knepley 
4008b438e21SMatthew G. Knepley TEST*/
401