Issue with the incompressible elasticity example of JuAFEM package

[Leooop]

Hi,

I plan to use the JuAFEM package as a basis for my phd work dealing with brittle damage modeling.
First of all, I would like to thank the authors(@kristoffer.carlsson et @fredrikekre) for this beautifully written package. I explored its interface by running some of the examples provided in the documentation. It turns out that I get weird results for simple cases.
I modified the example of the 2D incompressible elastic material to obtain a square domain with normal compressive boundary conditions on the upper and lower boundaries, without constraining displacements anywhere on the domain. The only changes made to the example are therefore :

This function is used instead of create_cook_grid :

function create_square_grid(nx, ny)
    grid = generate_grid(Triangle, (nx, ny), Vec(0.0,0.0), Vec(50.0,50.0));
    return grid
end;

I removed the Dirichlet boundary conditions in create_bc(dh), and the traction boundary conditions contributions to the force vector are replaced with:

@inbounds for face in 1:nfaces(cell)
        if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "top")
            reinit!(facevalues_u, cell, face)
            for q_point in 1:getnquadpoints(facevalues_u)
                dΓ = getdetJdV(facevalues_u, q_point)
                for i in 1:n_basefuncs_u
                    δu = shape_value(facevalues_u, q_point, i)
                    fe[i] += (δu ⋅ t) * dΓ
                end
            end
        elseif onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "bottom")
            reinit!(facevalues_u, cell, face)
            for q_point in 1:getnquadpoints(facevalues_u)
                dΓ = getdetJdV(facevalues_u, q_point)
                for i in 1:n_basefuncs_u
                    δu = shape_value(facevalues_u, q_point, i)
                    fe[i] += (δu ⋅ -t) * dΓ
                end
            end

        end
    end

with t = Vec{2}((0.0,-1e8)), Emod = 70e9 and ν = 0.4999999

With the quadratic/linear interpolation functions and a 50*50 elements domain, I’m getting a displacement field that doesn’t show the expected behavior.

Do you have any ideas as to why this is happening?
The same type of erroneous displacement field appears when I force only the vertical displacement of the upper and lower ends of the domain. The solution is only acceptable when the horizontal displacement of these same edges is set to zero.

I take this opportunity to ask you if you plan to develop a way to import meshes generated elsewhere?

Thank you!

[kristoffer.carlsson]

@Leooop

without constraining displacements anywhere on the domain.

Don’t you still need to at least prescribe enough Dirichlet boundary conditions to remove rigid body motion so the system will not be singular? This is generally the case for standard formulations of stationary continuum mechanical problems.

For your example you should divide the domain along the symmetry line and put symmetric Dirichlet boundary conditions on it.

@Leooop

I take this opportunity to ask you if you plan to develop a way to import meshes generated elsewhere?

Yeah, there is code spread around in various repos that does this for Abaqus and Gmsh files but we haven’t put any “official” way of doing it into the JuaFEM repo yet.

[mohamed82008]

@Leooop

I take this opportunity to ask you if you plan to develop a way to import meshes generated elsewhere?

I import inp files to JuAFEM in TopOpt using a regex-based parser https://github.com/mohamed82008/TopOpt.jl/tree/master/src/TopOptProblems/IO/INP/Parser. But not all the features of inp are supported. You might want to read and extend it.

[Leooop]

@kristoffer.carlsson

Don’t you still need to at least prescribe enough Dirichlet boundary conditions to remove rigid body motion so the system will not be singular? This is generally the case for standard formulations of stationary continuum mechanical problems.

Oh, makes sense indeed! It’s the first time I explore this kind of problem. I worked on a finite difference geodynamical model of the lithosphere before, but I never had to ask myself this question.

@kristoffer.carlsson

For your example you should divide the domain along the symmetry line and put symmetric Dirichlet boundary conditions on it.

Considering only the left part of my previous domain, and implementing Dirichlet boundary condition ux = 0 along the right boundary, as well as ux = uy = 0 on the bottom right corner node, with still the symmetric compressive Neumann boundary condition on the top and bottom boundaries, it now works like a charm, thanks!

Great to know that mesh import is available as well!

@mohamed82008

I import inp files to JuAFEM in TopOpt using a regex-based parser https://github.com/mohamed82008/TopOpt.jl/tree/master/src/TopOptProblems/IO/INP/Parser . But not all the features of inp are supported. You might want to read and extend it.

Thanks, i’ll take a look at it!