Requesting help with the error "det(J) is not positive" when using custom Delaunay grids

[zkk960317]

I need to use a random triangle mesh. But when I define my own mesh, and replace the original mesh in the "heat_equation.jl" example, I get the following error in function reinit!( ) :
ERROR: ArgumentError: det(J) is not positive: det(J) = -1.0736838778944306

The elements of the Delaunay grid is shown in the figure by Paraview, where green num is cellnum and yellow num is nodenum.

Following is the code:

using Ferrite, SparseArrays
using VoronoiDelaunay
function generate_delaunay2D_grid(dims, esize)
    r0=0.6 #Random factor
    ne = (Int(round(dims[1] / esize)), Int(round(dims[2] / esize)))
    ne_t1 = (ne[1] - 1) * (ne[2] - 1)
    ne_t2 = ne_t1 + 2 * (ne[1] + ne[2])
    rdx = rand(ne[1] - 1, ne[2] - 1) .* esize .* r0 .- esize * r0 * 0.5
    rdy = rand(ne[1] - 1, ne[2] - 1) .* esize .* r0 .- esize * r0 * 0.5
    for i = 1:ne[1]-1, j = 1:ne[2]-1
        rdx[i, j] += i * esize
        rdy[i, j] += j * esize
    end
    width = max_coord - min_coord
    ratio = 1 / max(dims[1], dims[2]) * width
    rdx = rdx * ratio
    rdy = rdy * ratio
    nodes = Point2D[Point(min_coord + rdx[i], min_coord + rdy[i]) for i in 1:ne_t1]
    ### add points at boundaries
    nodes = [nodes; [Point2D(min_coord, min_coord + esize * i * ratio) for i = 0:ne[2]];
        [Point2D(max_coord, min_coord + esize * i * ratio) for i = 0:ne[2]];
        [Point2D(min_coord + esize * i * ratio, min_coord) for i = 1:ne[1]-1];
        [Point2D(min_coord + esize * i * ratio, min_coord + width * dims[2] / dims[1]) for i = 1:ne[1]-1]]
    tess = DelaunayTessellation()
    sizehint!(tess, ne_t2)
    @time push!(tess, nodes)
    dic_nodes = Dict((nodes[i], i) for i = 1:lastindex(nodes))
    # i=0
    cells=Vector{Triangle}()
    sizehint!(cells,lastindex(tess._trigs))
    for DelaunayID in tess
        # i+=1
        cells =push!(cells, Triangle((get(dic_nodes, DelaunayID._a, 0),
        get(dic_nodes, DelaunayID._b, 0),
        get(dic_nodes, DelaunayID._c, 0))))
    end

    # cells = [Quadrilateral3D(cell.nodes) for cell in _grid.cells]
    nodes1 = [Node(((nodes[i]._x - min_coord) / ratio, (nodes[i]._y - min_coord) / ratio)) for i = 1:lastindex(nodes)]
    grid = Grid(cells, nodes1)

    return grid
end;

size = (20, 20)
grid = generate_delaunay2D_grid(size, 1)
# grid = generate_delaunay_grid(Quadrilateral, (20, 20));

dim = 2
ip = Lagrange{dim,RefTetrahedron,1}()
qr = QuadratureRule{dim,RefTetrahedron}(2)
cellvalues = CellScalarValues(qr, ip);

dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh);

K = create_sparsity_pattern(dh)

ch = ConstraintHandler(dh);

addfaceset!(grid, "left", (x) -> x[1] ≈ 0.0)
addfaceset!(grid, "right", (x) -> x[1] ≈ size[1])
addfaceset!(grid, "top", (x) -> x[2] ≈ 0.0)
addfaceset!(grid, "bottom", (x) -> x[2] ≈ size[2])

∂Ω = union(getfaceset.((grid, ), ["left", "right", "top", "bottom"])...);

dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
add!(ch, dbc);

close!(ch)
update!(ch, 0.0);

function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellScalarValues)
    n_basefuncs = getnbasefunctions(cellvalues)
    # Reset to 0
    fill!(Ke, 0)
    fill!(fe, 0)
    # Loop over quadrature points
    for q_point in 1:getnquadpoints(cellvalues)
        # Get the quadrature weight
        dΩ = getdetJdV(cellvalues, q_point)
        # Loop over test shape functions
        for i in 1:n_basefuncs
            δu = shape_value(cellvalues, q_point, i)
            ∇δu = shape_gradient(cellvalues, q_point, i)
            # Add contribution to fe
            fe[i] += δu * dΩ
            # Loop over trial shape functions
            for j in 1:n_basefuncs
                ∇u = shape_gradient(cellvalues, q_point, j)
                # Add contribution to Ke
                Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
            end
        end
    end
    return Ke, fe
end

function assemble_global(cellvalues::CellScalarValues, K::SparseMatrixCSC, dh::DofHandler)
    # Allocate the element stiffness matrix and element force vector
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    # Allocate global force vector f
    f = zeros(ndofs(dh))
    # Create an assembler
    assembler = start_assemble(K, f)
    # Loop over all cels
    for cell in CellIterator(dh)
        # Reinitialize cellvalues for this cell
        reinit!(cellvalues, cell)
        # Compute element contribution
        assemble_element!(Ke, fe, cellvalues)
        # Assemble Ke and fe into K and f
        assemble!(assembler, celldofs(cell), Ke, fe)
    end
    return K, f
end

K, f = assemble_global(cellvalues, K, dh);


apply!(K, f, ch)
u = K \ f;

# vtk_grid("heat_equation", dh) do vtk
#     vtk_point_data(vtk, [1:lastindex(grid.nodes);], "NodeNum")
#     vtk_cell_data(vtk, [1:lastindex(grid.cells);], "CellNum")
# end

[koehlerson]

The order of your elements is opposite to the one used in Ferrite, which produces a det(J) which is probably correct, but with negative sign. Same issue can happen when reading a mesh from gmsh Ferrite-FEM/FerriteGmsh.jl#15

I think you can do

        cells =push!(cells, Triangle((get(dic_nodes, DelaunayID._a, 0),
        get(dic_nodes, DelaunayID._c, 0),
        get(dic_nodes, DelaunayID._b, 0))))

[zkk960317]

The order of your elements is opposite to the one used in Ferrite, which produces a det(J) which is probably correct, but with negative sign. Same issue can happen when reading a mesh from gmsh Ferrite-FEM/FerriteGmsh.jl#15

I think you can do

        cells =push!(cells, Triangle((get(dic_nodes, DelaunayID._a, 0),
        get(dic_nodes, DelaunayID._c, 0),
        get(dic_nodes, DelaunayID._b, 0))))

Thank you very much! The problem has been solved.