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!