xref: /petsc/src/dm/impls/plex/tests/ex37.c (revision 6dd63270497ad23dcf16ae500a87ff2b2a0b7474)
1 static const char help[] = "Test of CAD functionality";
2 
3 #include <petscdmplex.h>
4 
5 /* TODO
6   - Fix IGES
7   - Test tessellation using -dm_plex_egads_with_tess
8 */
9 
10 typedef struct {
11   char      filename[PETSC_MAX_PATH_LEN];
12   PetscBool volumeMesh;
13 } AppCtx;
14 
15 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
16 {
17   PetscFunctionBeginUser;
18   options->filename[0] = '\0';
19   options->volumeMesh  = PETSC_TRUE;
20 
21   PetscOptionsBegin(comm, "", "EGADSPlex Problem Options", "EGADSlite");
22   PetscCall(PetscOptionsString("-filename", "The CAD file", "ex37.c", options->filename, options->filename, sizeof(options->filename), NULL));
23   PetscCall(PetscOptionsBool("-volume_mesh", "Create a volume mesh", "ex37.c", options->volumeMesh, &options->volumeMesh, NULL));
24   PetscOptionsEnd();
25   PetscFunctionReturn(PETSC_SUCCESS);
26 }
27 
28 static PetscErrorCode ComputeVolume(DM dm)
29 {
30   PetscObject obj = (PetscObject)dm;
31   DMLabel     bodyLabel, faceLabel, edgeLabel;
32   double      surface = 0., volume = 0., vol;
33   PetscInt    dim, pStart, pEnd, p, pid;
34   const char *name;
35 
36   PetscFunctionBeginUser;
37   PetscCall(DMGetDimension(dm, &dim));
38   if (dim < 2) PetscFunctionReturn(PETSC_SUCCESS);
39   PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
40   PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
41   PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
42 
43   PetscCall(DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd));
44   for (p = pStart; p < pEnd; ++p) {
45     PetscCall(DMLabelGetValue(dim == 2 ? faceLabel : bodyLabel, p, &pid));
46     if (pid >= 0) {
47       PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &vol, NULL, NULL));
48       volume += vol;
49     }
50   }
51   PetscCall(DMPlexGetHeightStratum(dm, 1, &pStart, &pEnd));
52   for (p = pStart; p < pEnd; ++p) {
53     PetscCall(DMLabelGetValue(dim == 2 ? edgeLabel : faceLabel, p, &pid));
54     if (pid >= 0) {
55       PetscCall(DMPlexComputeCellGeometryFVM(dm, p, &vol, NULL, NULL));
56       surface += vol;
57     }
58   }
59 
60   PetscCall(PetscObjectGetName(obj, &name));
61   PetscCall(PetscPrintf(PetscObjectComm(obj), "DM %s: Surface Area = %.6e Volume = %.6e\n", name ? name : "", surface, volume));
62   PetscFunctionReturn(PETSC_SUCCESS);
63 }
64 
65 int main(int argc, char *argv[])
66 {
67   DM     surface, dm;
68   AppCtx ctx;
69 
70   PetscFunctionBeginUser;
71   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
72   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
73   PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ctx.filename, "ex37_plex", PETSC_TRUE, &surface));
74   PetscCall(PetscObjectSetName((PetscObject)surface, "CAD Surface"));
75   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)surface, "sur_"));
76   PetscCall(DMSetFromOptions(surface));
77   PetscCall(DMViewFromOptions(surface, NULL, "-dm_view"));
78   PetscCall(ComputeVolume(surface));
79 
80   if (ctx.volumeMesh) {
81     PetscCall(DMPlexGenerate(surface, "tetgen", PETSC_TRUE, &dm));
82     PetscCall(PetscObjectSetName((PetscObject)dm, "CAD Mesh"));
83     PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE));
84     PetscCall(DMViewFromOptions(dm, NULL, "-pre_dm_view"));
85 
86     PetscCall(DMPlexInflateToGeomModel(dm, PETSC_TRUE));
87     PetscCall(DMViewFromOptions(dm, NULL, "-inf_dm_view"));
88 
89     PetscCall(DMSetFromOptions(dm));
90     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
91     PetscCall(ComputeVolume(dm));
92     PetscCall(DMDestroy(&dm));
93   }
94   PetscCall(DMDestroy(&surface));
95   PetscCall(PetscFinalize());
96   return 0;
97 }
98 
99 /*TEST
100 
101   build:
102     requires: egads tetgen datafilespath
103 
104   test:
105     suffix: sphere_0
106     args: -filename ${DATAFILESPATH}/meshes/cad/sphere_example.egadslite -dm_refine 1 -sur_dm_view -dm_plex_check_all -dm_plex_geom_print_model -sur_dm_plex_view_labels "EGADS Body ID","EGADS Face ID","EGADS Edge ID"
107 
108   test:
109     suffix: sphere_egads
110     args: -filename ${DATAFILESPATH}/meshes/cad/sphere_example.egads -dm_refine 1 -sur_dm_view -dm_plex_check_all -dm_plex_geom_print_model -sur_dm_plex_view_labels "EGADS Body ID","EGADS Face ID","EGADS Edge ID"
111 
112   test:
113     suffix: sphere_iges
114     TODO: broken
115     args: -filename ${DATAFILESPATH}/meshes/cad/sphere_example.igs -dm_refine 1 -sur_dm_view -dm_plex_check_all -dm_plex_geom_print_model -sur_dm_plex_view_labels "EGADS Body ID","EGADS Face ID","EGADS Edge ID"
116 
117   test:
118     suffix: sphere_step
119     args: -filename ${DATAFILESPATH}/meshes/cad/sphere_example.stp -dm_refine 1 -sur_dm_view -dm_plex_check_all -dm_plex_geom_print_model -sur_dm_plex_view_labels "EGADS Body ID","EGADS Face ID","EGADS Edge ID"
120 
121   test:
122     suffix: nozzle_0
123     args: -filename ${DATAFILESPATH}/meshes/cad/nozzle.egadslite -sur_dm_refine 1 -sur_dm_view -dm_plex_check_all
124 
125   test:
126     suffix: nozzle_egads
127     args: -filename ${DATAFILESPATH}/meshes/cad/nozzle.egads -sur_dm_refine 1 -sur_dm_view -dm_plex_check_all
128 
129   test:
130     suffix: nozzle_iges
131     args: -filename ${DATAFILESPATH}/meshes/cad/nozzle.igs -sur_dm_refine 1 -sur_dm_view -dm_plex_check_all
132 
133   test:
134     suffix: nozzle_step
135     args: -filename ${DATAFILESPATH}/meshes/cad/nozzle.stp -sur_dm_refine 1 -sur_dm_view -dm_plex_check_all
136 
137 TEST*/
138