1 static char help[] = "Tests DMPlex Gmsh reader.\n\n"; 2 3 #include <petscdmplex.h> 4 5 #if !defined(PETSC_GMSH_EXE) 6 #define PETSC_GMSH_EXE "gmsh" 7 #endif 8 9 #include <petscds.h> 10 11 static void one(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar value[]) 12 { 13 value[0] = (PetscReal)1; 14 } 15 16 static PetscErrorCode CreateFE(DM dm) 17 { 18 DM cdm; 19 PetscSpace P; 20 PetscDualSpace Q; 21 DM K; 22 PetscFE fe; 23 DMPolytopeType ptype; 24 25 PetscInt dim, k; 26 PetscBool isSimplex; 27 28 PetscDS ds; 29 30 PetscFunctionBeginUser; 31 PetscCall(DMGetCoordinateDM(dm, &cdm)); 32 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe)); 33 PetscCall(PetscFEGetBasisSpace(fe, &P)); 34 PetscCall(PetscFEGetDualSpace(fe, &Q)); 35 PetscCall(PetscDualSpaceGetDM(Q, &K)); 36 PetscCall(DMGetDimension(K, &dim)); 37 PetscCall(PetscSpaceGetDegree(P, &k, NULL)); 38 PetscCall(DMPlexGetCellType(K, 0, &ptype)); 39 switch (ptype) { 40 case DM_POLYTOPE_QUADRILATERAL: 41 case DM_POLYTOPE_HEXAHEDRON: 42 isSimplex = PETSC_FALSE; 43 break; 44 default: 45 isSimplex = PETSC_TRUE; 46 break; 47 } 48 49 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, isSimplex, k, PETSC_DETERMINE, &fe)); 50 PetscCall(PetscFESetName(fe, "scalar")); 51 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 52 PetscCall(PetscFEDestroy(&fe)); 53 PetscCall(DMCreateDS(dm)); 54 55 PetscCall(DMGetDS(dm, &ds)); 56 PetscCall(PetscDSSetObjective(ds, 0, one)); 57 PetscFunctionReturn(PETSC_SUCCESS); 58 } 59 60 static PetscErrorCode CheckIntegral(DM dm, PetscReal integral, PetscReal tol) 61 { 62 Vec u; 63 PetscReal rval; 64 PetscScalar result; 65 66 PetscFunctionBeginUser; 67 PetscCall(DMGetGlobalVector(dm, &u)); 68 PetscCall(DMPlexComputeIntegralFEM(dm, u, &result, NULL)); 69 PetscCall(DMRestoreGlobalVector(dm, &u)); 70 rval = PetscRealPart(result); 71 if (integral > 0 && PetscAbsReal(integral - rval) > tol) { 72 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Calculated value %g != %g actual value (error %g > %g tol)\n", (double)rval, (double)integral, (double)PetscAbsReal(integral - rval), (double)tol)); 73 } 74 PetscFunctionReturn(PETSC_SUCCESS); 75 } 76 77 int main(int argc, char **argv) 78 { 79 DM dm; 80 const char *const mshlist[] = {"seg", "tri", "qua", "tet", "wed", "hex", "vtx", "B2tri", "B2qua", "B3tet", "B3hex"}; 81 const char *const fmtlist[] = {"msh22", "msh40", "msh41"}; 82 PetscInt msh = 5; 83 PetscInt fmt = 2; 84 PetscBool bin = PETSC_TRUE; 85 PetscInt dim = 3; 86 PetscInt order = 2; 87 88 const char cmdtemplate[] = "%s -format %s %s -%d -order %d %s -o %s"; 89 char gmsh[PETSC_MAX_PATH_LEN] = PETSC_GMSH_EXE; 90 char tag[PETSC_MAX_PATH_LEN], path[PETSC_MAX_PATH_LEN]; 91 char geo[PETSC_MAX_PATH_LEN], geodir[PETSC_MAX_PATH_LEN] = "."; 92 char out[PETSC_MAX_PATH_LEN], outdir[PETSC_MAX_PATH_LEN] = "."; 93 char cmd[PETSC_MAX_PATH_LEN * 4]; 94 PetscBool set, flg; 95 FILE *fp; 96 97 PetscFunctionBeginUser; 98 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 99 100 PetscCall(PetscStrncpy(geodir, "${PETSC_DIR}/share/petsc/datafiles/meshes", sizeof(geodir))); 101 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "GMSH", path, sizeof(path), &set)); 102 if (set) PetscCall(PetscStrncpy(gmsh, path, sizeof(gmsh))); 103 PetscCall(PetscOptionsGetString(NULL, NULL, "-gmsh", gmsh, sizeof(gmsh), NULL)); 104 PetscCall(PetscOptionsGetString(NULL, NULL, "-dir", geodir, sizeof(geodir), NULL)); 105 PetscCall(PetscOptionsGetString(NULL, NULL, "-out", outdir, sizeof(outdir), NULL)); 106 PetscCall(PetscOptionsGetEList(NULL, NULL, "-msh", mshlist, PETSC_STATIC_ARRAY_LENGTH(mshlist), &msh, NULL)); 107 PetscCall(PetscOptionsGetEList(NULL, NULL, "-fmt", fmtlist, PETSC_STATIC_ARRAY_LENGTH(fmtlist), &fmt, NULL)); 108 PetscCall(PetscOptionsGetBool(NULL, NULL, "-bin", &bin, NULL)); 109 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL)); 110 PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL)); 111 if (fmt == 1) bin = PETSC_FALSE; /* Recent Gmsh releases cannot generate msh40+binary format*/ 112 113 { /* This test requires Gmsh >= 4.2.0 */ 114 char space[PETSC_MAX_PATH_LEN]; 115 int inum = 0, major = 0, minor = 0, micro = 0; 116 PetscCall(PetscSNPrintf(cmd, sizeof(cmd), "%s -info", gmsh)); 117 PetscCall(PetscPOpen(PETSC_COMM_SELF, NULL, cmd, "r", &fp)); 118 if (fp) inum = fscanf(fp, "Version %s %d.%d.%d", space, &major, &minor, µ); 119 PetscCall(PetscPClose(PETSC_COMM_SELF, fp)); 120 if (inum != 4 || major < 4 || (major == 4 && minor < 2)) { 121 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Gmsh>=4.2.0 not available\n")); 122 goto finish; 123 } 124 } 125 126 PetscCall(PetscSNPrintf(tag, sizeof(tag), "%s-%d-%d-%s%s", mshlist[msh], (int)dim, (int)order, fmtlist[fmt], bin ? "-bin" : "")); 127 PetscCall(PetscSNPrintf(geo, sizeof(geo), "%s/gmsh-%s.geo", geodir, mshlist[msh])); 128 PetscCall(PetscSNPrintf(out, sizeof(out), "%s/mesh-%s.msh", outdir, tag)); 129 PetscCall(PetscStrreplace(PETSC_COMM_SELF, geo, path, sizeof(path))); 130 PetscCall(PetscFixFilename(path, geo)); 131 PetscCall(PetscStrreplace(PETSC_COMM_SELF, out, path, sizeof(path))); 132 PetscCall(PetscFixFilename(path, out)); 133 PetscCall(PetscTestFile(geo, 'r', &flg)); 134 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "File not found: %s", geo); 135 136 PetscCall(PetscSNPrintf(cmd, sizeof(cmd), cmdtemplate, gmsh, fmtlist[fmt], bin ? "-bin" : "", (int)dim, (int)order, geo, out)); 137 PetscCall(PetscPOpen(PETSC_COMM_SELF, NULL, cmd, "r", &fp)); 138 PetscCall(PetscPClose(PETSC_COMM_SELF, fp)); 139 140 PetscCall(DMPlexCreateFromFile(PETSC_COMM_SELF, out, "ex99_plex", PETSC_TRUE, &dm)); 141 PetscCall(PetscSNPrintf(tag, sizeof(tag), "mesh-%s", mshlist[msh])); 142 PetscCall(PetscObjectSetName((PetscObject)dm, tag)); 143 PetscCall(DMSetFromOptions(dm)); 144 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 145 { 146 PetscBool check; 147 PetscReal integral = 0, tol = (PetscReal)1.0e-4; 148 PetscCall(PetscOptionsGetReal(NULL, NULL, "-integral", &integral, &check)); 149 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 150 if (check) { 151 PetscCall(CreateFE(dm)); 152 PetscCall(CheckIntegral(dm, integral, tol)); 153 } 154 } 155 PetscCall(DMDestroy(&dm)); 156 157 finish: 158 PetscCall(PetscFinalize()); 159 return 0; 160 } 161 162 /*TEST 163 164 build: 165 requires: defined(PETSC_HAVE_POPEN) 166 167 test: 168 requires: defined(PETSC_GMSH_EXE) 169 args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes 170 args: -msh {{vtx}separate_output} 171 args: -order 1 172 args: -fmt {{msh22 msh40 msh41}} -bin {{0 1}} 173 args: -dm_view ::ascii_info_detail 174 args: -dm_plex_check_all 175 args: -dm_plex_gmsh_highorder false 176 args: -dm_plex_boundary_label marker 177 args: -dm_plex_gmsh_spacedim 3 178 179 test: 180 requires: defined(PETSC_GMSH_EXE) 181 args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes 182 args: -msh {{seg tri qua tet wed hex}separate_output} 183 args: -order {{1 2 3 7}} 184 args: -fmt {{msh22 msh40 msh41}} -bin {{0 1}} 185 args: -dm_view ::ascii_info_detail 186 args: -dm_plex_check_all 187 args: -dm_plex_gmsh_highorder false 188 args: -dm_plex_boundary_label marker 189 190 testset: 191 suffix: B2 # 2D ball 192 output_file: output/empty.out 193 requires: defined(PETSC_GMSH_EXE) 194 args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes 195 args: -msh {{B2tri B2qua}} 196 args: -dim 2 -integral 3.141592653589793 # pi 197 args: -order {{2 3 4 5 6 7 8 9}} -tol 0.05 198 199 testset: 200 suffix: B2_bnd # 2D ball boundary 201 output_file: output/empty.out 202 requires: defined(PETSC_GMSH_EXE) 203 args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes 204 args: -dm_plex_gmsh_spacedim 2 205 args: -msh {{B2tri B2qua}} 206 args: -dim 1 -integral 6.283185307179586 # 2*pi 207 args: -order {{2 3 4 5 6 7 8 9}} -tol 0.05 208 209 testset: 210 suffix: B3 # 3D ball 211 output_file: output/empty.out 212 requires: defined(PETSC_GMSH_EXE) 213 args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes 214 args: -msh {{B3tet B3hex}} 215 args: -dim 3 -integral 4.1887902047863905 # 4/3*pi 216 args: -order {{2 3 4 5}} -tol 0.20 217 218 testset: 219 suffix: B3_bnd # 3D ball boundary 220 output_file: output/empty.out 221 requires: defined(PETSC_GMSH_EXE) 222 args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes 223 args: -dm_plex_gmsh_spacedim 3 224 args: -msh {{B3tet B3hex}} 225 args: -dim 2 -integral 12.566370614359172 # 4*pi 226 args: -order {{2 3 4 5 6 7 8 9}} -tol 0.25 227 228 testset: 229 suffix: B_lin # linear discretizations 230 requires: defined(PETSC_GMSH_EXE) 231 args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes 232 args: -dm_plex_gmsh_highorder true 233 args: -dm_plex_gmsh_project true 234 args: -dm_plex_gmsh_project_petscspace_degree {{1 2 3}separate_output} 235 args: -dm_plex_gmsh_fe_view 236 args: -dm_plex_gmsh_project_fe_view 237 args: -order 1 -tol 1e-4 238 test: 239 suffix: dim-1 240 args: -dm_plex_gmsh_spacedim 2 241 args: -msh {{B2tri B2qua}separate_output} 242 args: -dim 1 -integral 5.656854249492381 # 4*sqrt(2) 243 test: 244 suffix: dim-2 245 args: -dm_plex_gmsh_spacedim 2 246 args: -msh {{B2tri B2qua}separate_output} 247 args: -dim 2 -integral 2.000000000000000 # 2 248 test: 249 suffix: dim-2_msh-B3tet 250 args: -dm_plex_gmsh_spacedim 3 251 args: -msh B3tet -dim 2 -integral 9.914478 252 test: 253 suffix: dim-2_msh-B3hex 254 args: -dm_plex_gmsh_spacedim 3 255 args: -msh B3hex -dim 2 -integral 8.000000 256 test: 257 suffix: dim-3_msh-B3tet 258 args: -dm_plex_gmsh_spacedim 3 259 args: -msh B3tet -dim 3 -integral 2.666649 260 test: 261 suffix: dim-3_msh-B3hex 262 args: -dm_plex_gmsh_spacedim 3 263 args: -msh B3hex -dim 3 -integral 1.539600 264 265 TEST*/ 266