Working on a Import Mesh Workaround (FerriteGmsh.jl Issue 15)

[jamesmichaelgorman]

This topic may technically be better off in FerriteGmsh, but I am hoping more people will see this and have a suggestion.

It has been known for a little bit (Ferrite-FEM/FerriteGmsh.jl#15), that when FerriteGmsh.jl imports a mesh via Gmsh that nodes within elements may be clockwise instead of counter-clockwise (for an unknown reason, as far as I can tell). I want to use Ferrite to solve problems that require non-trivial geometries, and that means importing a mesh.

My workaround is slow because I intend to cycle through every element. However, right now I am slowed down because I don't appear to be able to easily extract out nodes number from grid.cells[1], but I have not yet found an easy way to get node numbers out from the output. While I can easily output the node values, such as Triangle((1, 3, 2)) or Quadrilateral((1, 5, 9, 7)), I cannot appear to easily extract them out individually. Does someone know a function that would allow me to get that done? Currently, the only way I see to brute force this correction is to use getcoordinates() and then match the coordinates to nodes in a list (inefficient, but should work). Does anyone else have any ideas?

James

[lijas]

While I can easily output the node values, such as Triangle((1, 3, 2)) or Quadrilateral((1, 5, 9, 7)), I cannot appear to easily extract them out individually.

You can get the nodesids of a cell with grid.cells[cellid].nodes . This is a tuple which you can index as normal.

But I think the it would be better to adress the original issue: Can you orient your Gmsh surfaces correctly such that counter-clock-wise elements are generated?

[jamesmichaelgorman]

Thanks. I hadn't yet found it in the documentation.

I don't know why some .msh files are imported incorrectly by FerriteGmsh.jl. I was able to import the example .msh file into Gridap and it worked just fine, but not in Ferrite. Therefore, as mentioned in the link, this appears to be a FerriteGmsh issue (for some imported meshes) and not Gmsh output issue.

The FerriteGmsh mesh import issue hasn't been resolved yet and I don't know which cases where an imported mesh causes problems and which don't. Unfortunately, I don't (yet) have the time to provide a complete solution. I was, however, able to generate a partial solution for 1st order 2D elements: I confirmed the quadrilateral mesh worked after being sorted, and it is expected triangle mesh would work as well (not confirmed, since I typically only use hypercubes for meshing). It is provided below.

using Ferrite, FerriteGmsh, SparseArrays, Gmsh, DataFrames

function sort_grid(grid)
    sorted_el_list = Vector{typeof(grid.cells[1])}()
    sizehint!(sorted_el_list, length(grid.cells))

    for a=1:length(grid.cells)
        # first, check to see if element is oriented correctly already
        el_nodes = DataFrame(num=grid.cells[a].nodes[1],x=getcoordinates(grid,a)[1][1],y=getcoordinates(grid,a)[1][2])
        for b=2:length(grid.cells[a].nodes)
            append!(el_nodes, DataFrame(num=grid.cells[a].nodes[b],x=getcoordinates(grid,a)[b][1],y=getcoordinates(grid,a)[b][2]))
        end
        # println(el_nodes)
        centroid_x = mean(el_nodes.x)
        centroid_y = mean(el_nodes.y)
        # calculate angle from the centroid
        coord_angle = atan.(el_nodes.y .- centroid_y, el_nodes.x .- centroid_x)
        # check to see if nodes are sorted based off angles
        # angles need to go from negative to positive, therefore diff(x) > 0 
        if all(>(0), diff(coord_angle))==true
            push!(sorted_el_list, copy(grid.cells[a]))
            continue
        end
        # otherwise, we need to reorder elements
        # **this only works for 1st order 2D**
        if typeof(grid.cells[a])==Ferrite.Quadrilateral || typeof(grid.cells[a])==Ferrite.Triangle
            el_nodes[!, :coord_angle] = coord_angle
            sort!(el_nodes, [:coord_angle])
            new_order = el_nodes.num
            if typeof(grid.cells[a])==Ferrite.Quadrilateral
                push!(sorted_el_list, Quadrilateral((new_order[1], new_order[2], new_order[3], new_order[4])))
            elseif typeof(grid.cells[a])==Ferrite.Triangle
                push!(sorted_el_list, Triangle((new_order[1], new_order[2], new_order[3])))
            else
                error("Ferrite element type specified other than Quadrilateral or Triangle")
            end
        else
            error("This function currently only works for 1st order 2D elements")
        end
    end
    new_grid = Ferrite.Grid(sorted_el_list, grid.nodes, grid.cellsets, grid.nodesets, grid.facetsets, grid.vertexsets)
    return new_grid
end```

[lijas]

Do you still have access to the .geo file that was used to generate the .msh file? I am guessing there is some surface definition in the .geo that is oriented the "wrong" way.

You can also try to reverse the elements using the snippet below (but it requires the .geo file i think):

if dim == 2
        surfacetags = gmsh.model.get_entities(2) #get all surfaces
        for (surfacedim, surfacetag) in surfacetags
            n = gmsh.model.get_normal(surfacetag, [0.0,0.0])
            val = dot(n, [0.0,0.0,1.0]) 
            if val ≈ -1.0 #wrong orientation
                _, elementTags, _ = gmsh.model.mesh.getElements(surfacedim, surfacetag)
                gmsh.model.mesh.reverse_elements(elementTags[1])
            end
        end
    end