Unpaired faces in 2D mesh

Errors
1

Hi there,

I am starting using pyFR for 2d incompressible flow around obstacles. The code used for generating meshes looks good, as I can visualize a body-fitted triangular mesh by opening the .msh file with gmsh. But the .msh file fails to be imported to .pyfrm due to unpaired faces. I do not understand where the problem comes from. Please check the meshing function I used. I cannot upload the .msh file by the way.

Thanks in advance for any suggestions.

    def OMesh(self, extruded=False, xmin = -10.0, xmax = 10.0, ymin = -10.0, ymax = 10.0, obstacle_element_scale = 0.2, edge_method='polygon', laplace_smoothing_iter = 200, output_dir='./'):
        def add_polygon(geom, vertices_tags):
            lines = []
            for (line_start, line_end) in zip(vertices_tags[:-1], vertices_tags[1:]):
                lines.append(geom.add_line(line_start, line_end))
            lines.append(geom.add_line(vertices_tags[-1],vertices_tags[0]))
            return lines
        
        obstacle_vertices_coords = self.curve_pts
        obstacle_top_coord = np.max(obstacle_vertices_coords,axis=0)[1]
        obstacle_bottom_coord = np.min(obstacle_vertices_coords,axis=0)[1]
        centroid = np.mean(obstacle_vertices_coords,axis=0)

        gmsh.initialize()
        mod = gmsh.model
        geom = mod.geo()

        # curvatures, is_convex = curvature(self.curve_pts[:,:2])
        obstacle_vertices = []
        for pt in obstacle_vertices_coords:
            pt_tag = geom.addPoint(*pt,obstacle_element_scale)
            obstacle_vertices.append(pt_tag)
        if edge_method == 'polygon':
            obstacle_edge = add_polygon(geom, obstacle_vertices)
        elif edge_method == 'polyline':
            obstacle_edge = [geom.addPolyline(obstacle_vertices + [obstacle_vertices[0]])]
        elif edge_method == 'spline':
            obstacle_edge = [geom.add_spline(obstacle_vertices + [obstacle_vertices[0]])]
        else:
            raise(ValueError('Invalid edge_method {edge_method}; choose one of polygon polyline or spline'))
        obstacle_loop = geom.add_curve_loop(obstacle_edge)

        outer_vertices = []
        outer_vertices.append(geom.add_point(xmax,ymax,0.0))
        outer_vertices.append(geom.add_point(xmax,ymin,0.0))
        outer_vertices.append(geom.add_point(xmin,ymin,0.0))
        outer_vertices.append(geom.add_point(xmin,ymax,0.0))
        outer_lines = add_polygon(geom, outer_vertices)
        outer_loop = geom.add_curve_loop(outer_lines)
        surf = geom.addPlaneSurface([outer_loop, obstacle_loop])

        minaniso_fields = []
#        if boundary_layer_parameters is not None:
        boundary_layer_default_parameters = {
            "Thickness": np.mean(np.linalg.norm(self.curve_pts-centroid,axis=1)),
            "Quads": 1,
            "Size": 0.01,
            "SizeFar": 0.1,
            "IntersectMetrics": 1
        }
#        boundary_layer_parameters = boundary_layer_parameters or {}
        boundary_layer_parameters = {**boundary_layer_default_parameters}
        blayer = mod.mesh.field.add("BoundaryLayer")
        mod.mesh.field.setAsBoundaryLayer(blayer)
        for param in boundary_layer_parameters:
            mod.mesh.field.setNumber(blayer, param, boundary_layer_parameters[param])
        mod.mesh.field.setNumbers(blayer, "CurvesList", [obstacle_loop])
        minaniso_fields.append(blayer)
    
        try:
            bl_offset = 1.25*boundary_layer_parameters["Thickness"]
        except:
            bl_offset = 0
        wake_refinement_default_parameters = {
            "VIn": 0.2,
            "VOut": 2.0,
            "XMax": xmax,
            "XMin": 0.0,
            "YMax": obstacle_top_coord + bl_offset,
            "YMin": obstacle_bottom_coord - bl_offset
        }
        wake_refinement_parameters = {**wake_refinement_default_parameters}
        wrbox = mod.mesh.field.add("Box")
        for param in wake_refinement_parameters:
            mod.mesh.field.setNumber(wrbox, param, wake_refinement_parameters[param])
        minaniso_fields.append(wrbox)

        if len(minaniso_fields)>0:
            minaniso = mod.mesh.field.add("MinAniso")
            mod.mesh.field.setNumbers(minaniso, "FieldsList", minaniso_fields)
            mod.mesh.field.setAsBackgroundMesh(minaniso)
            
        geom.synchronize()

        pg_out = mod.addPhysicalGroup(1,[outer_lines[0]])
        mod.setPhysicalName(1,pg_out,"out")
        
        pg_sym1 = mod.addPhysicalGroup(1,[outer_lines[1]])
        mod.setPhysicalName(1,pg_sym1,"sym1")
        
        pg_in = mod.addPhysicalGroup(1,[outer_lines[2]])
        mod.setPhysicalName(1,pg_in,"in")
        
        pg_sym2 = mod.addPhysicalGroup(1,[outer_lines[3]])
        mod.setPhysicalName(1,pg_sym2,"sym2")
        
        pg_obstacle = mod.addPhysicalGroup(1,[obstacle_loop])
        mod.setPhysicalName(1,pg_obstacle,"obstacle")
        
        pg_fluid = mod.addPhysicalGroup(2,[surf])
        mod.setPhysicalName(2,pg_fluid,"fluid")
        mod.mesh.generate(2)
        mod.mesh.optimize("Laplace2D", niter=laplace_smoothing_iter)

            
        filename = f'{output_dir}/{self.name}_{self.index}.msh'
        gmsh.option.setNumber("Mesh.MshFileVersion",2.2)
        gmsh.write(filename)
        gmsh.write(filename[:-4]+ '.vtk')

        gmsh.finalize()

        return filename,0
2

This means that pyfr was unable to pair the faces. Do you have any periodic faces in your mesh? Currently PyFR can only pair periodic faces if they are tagged as such and where the coordinates of either side differ by a translation in one cardinal direction.

If that is not the issue, I would suggest making the mesh as small as possible and having a look at the connectivity to make sure it’s doing what you think.

3

I am creating a 2d rectangle with the obstacle cut out as a hole. No periodic faces should be there.

Thank you for the advice. I will make a smaller mesh to check what is happening.