What is 'move'?

[r614]

Hey! I am trying to use your library for running FEA on a 3D Mesh.

I was curious what the "move" keyword in the example means here:

fe.dof.uniaxial(displacement, move=0.2, clamped=True)

Is it the magnitude of displacement per mesh unit along the displacement field, or the force applied?

If, let's say I want to see how a 3D object deforms under gravity, do I just create a displacement field along Y, and add move = -9.8 or do I need to compute the actual displacement first?

Cheers!

[r614]

Also, as a side note, I get this error when trying to initialize my region:

ValueError                                Traceback (most recent call last)
/var/folders/56/ccm6q3vs0sn1vz8q9__yt_180000gn/T/ipykernel_32703/610656645.py in <module>
      5 quadrature = fe.GaussLegendre(order=1, dim=3)
      6 
----> 7 region = fe.RegionHexahedron(fe_mesh)

~/miniforge3/envs/capstone/lib/python3.9/site-packages/felupe/region/_templates.py in __init__(self, mesh)
    101         quadrature = GaussLegendre(order=1, dim=3)
    102 
--> 103         super().__init__(mesh, element, quadrature)
    104 
    105 

~/miniforge3/envs/capstone/lib/python3.9/site-packages/felupe/region/_region.py in __init__(self, mesh, element, quadrature, grad)
     67             #
     68             # dXdr_IJpe
---> 69             self.dXdr = np.einsum(
     70                 "caI,aJp->IJpc", self.mesh.points[self.mesh.cells], self.dhdr
     71             )
ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (290,3,3)->(3,newaxis,newaxis,290,3) (8,3,8)->(3,8,newaxis,8) 

The code:

import felupe as fe 
import meshio 

mesh = meshio.read("../data/models/crosshatch.msh")
# mesh.points.shape = (132, 3)
# mesh.cells[0][1].shape = (290, 3)
# mesh.cells[0][0] = 'triangle'


fe_mesh = fe.Mesh(mesh.points, mesh.cells[0][1])
# element = fe.Hexahedron()
# quadrature = fe.GaussLegendre(order=1, dim=3)

region = fe.RegionHexahedron(fe_mesh)

[adtzlr]

Hi,

for gravity you have to assemble the force vector by yourself.

a) the move argument of uniaxial refers to the absolute displacement on the right end-face in x-direction of a cube (as visualized in the README). For an example how to construct boundaries, see the Getting Started example.

b) You can't use RegionHexahedron for a mesh with Triangles. Use RegionTriangle instead or construct your own Region with a self-defined quadrature rule.

[r614]

Here's a link to the mesh!

RegionTriangle runs into an issue with IndexError, not sure how I can fix it tbh. 🤔 :

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/var/folders/56/ccm6q3vs0sn1vz8q9__yt_180000gn/T/ipykernel_32703/3375723294.py in <module>
      5 #quadrature = fe.GaussLegendre(order=1, dim=3)
      6 
----> 7 region = fe.RegionTriangle(fe_mesh)
      8 
      9 # Bugged, Github issue: https://github.com/adtzlr/felupe/issues/139

~/miniforge3/envs/capstone/lib/python3.9/site-packages/felupe/region/_templates.py in __init__(self, mesh)
    145         quadrature = TriangleQuadrature(order=1)
    146 
--> 147         super().__init__(mesh, element, quadrature)
    148 
    149 

~/miniforge3/envs/capstone/lib/python3.9/site-packages/felupe/region/_region.py in __init__(self, mesh, element, quadrature, grad)
     72 
     73             # inverse of dXdr
---> 74             self.drdX = inv(self.dXdr)
     75 
     76             # Determinant of geometric gradient evaluated at quadrature point "p"

~/miniforge3/envs/capstone/lib/python3.9/site-packages/felupe/math/_tensor.py in inv(A, determinant, full_output, sym)
     63 
     64     if determinant is None:
---> 65         detA = det(A)
     66     else:
     67         detA = determinant

~/miniforge3/envs/capstone/lib/python3.9/site-packages/felupe/math/_tensor.py in det(A)
    102     if A.shape[0] == 3:
    103         detA = (
--> 104             A[0, 0] * A[1, 1] * A[2, 2]
    105             + A[0, 1] * A[1, 2] * A[2, 0]
    106             + A[0, 2] * A[1, 0] * A[2, 1]

IndexError: index 2 is out of bounds for axis 1 with size 2

[adtzlr]

Thanks for the mesh - will have a look at what is wrong with RegionTriangle.

[r614]

Gotcha, so here's the problem. I generate the mesh using gmsh from a STL file, and it creates a triangle mesh which doesn't seem to work with either RegionTriangle or RegionHexahedron. Do you know any workarounds we might be able to try?

[adtzlr]

RegionTriangle has to be used with meshes consisting of triangular elements. RegionHexahedron: use a mesh with hex-elements, and so on. Just remove the third component of your mesh point array and you're fine. If you plan to do some shell-analysis: sorry, FElupe doesn't support this out of the box.

[adtzlr]

Also, plain strain and plain stress conditions have to be created manually (no umat templates included). FElupe is good for problems involving 3d solids (tetra or hexahedron elements) and axisymmetric ring elements (triangle or quad elements). 2d cases require some extra work. Just as a side note 📓 . I still hope you enjoy using FElupe and make some progress in the setup of your problem.

[r614]

Gotcha, that's pretty helpful. Our goal is basically to determine if a 3D structure (STL file) will collapse under gravity, so we are trying to run through a few different possibilities, haha. Do you have any recommendations for libraries/software we could take a look at to do this programmatically?

[adtzlr]

You should probably create a volume mesh (with tetrahedrons) in gmsh and import that in FElupe. What kind of material would you like to use? With the help of my other package MatADi a lot of hyperelastic materials are available. Do you plan to use something (nearly) incompressible or a compressible material? You could also have a look at scikit-fem.

[adtzlr]

Hi @r614, there is a new tutorial located at https://felupe.readthedocs.io/en/latest/tutorial/beam.html. Hope that helps!

[adtzlr]

Hey @r614, here is the solution. I re-created the mesh with pygmsh for you.

import felupe as fe
import numpy as np
import pygmsh

with pygmsh.occ.Geometry() as geom:

    geom.characteristic_length_min = 0.05
    geom.characteristic_length_max = 0.05

    rectangles = []
    holes = []
    
    rectangles.append(geom.add_rectangle([-0.9, 0.0, 0], 0.9, 0.9))
    rectangles.append(geom.add_rectangle([-1.0, 0.9, 0], 1.0, 1.0))
    rectangles.append(geom.add_rectangle([-1.1, 1.9, 0], 1.1, 1.1))
    rectangles.append(geom.add_rectangle([-1.2, 3.0, 0], 1.2, 1.2))
    
    rectangles.append(geom.add_rectangle([ 0.0, 0.0, 0], 0.8, 0.8))
    rectangles.append(geom.add_rectangle([ 0.0, 0.8, 0], 0.7, 0.7))
    rectangles.append(geom.add_rectangle([ 0.0, 1.5, 0], 0.6, 0.6))
    rectangles.append(geom.add_rectangle([ 0.0, 2.1, 0], 0.5, 0.5))
    rectangles.append(geom.add_rectangle([ 0.0, 2.6, 0], 0.4, 0.4))

    holes.append(geom.add_rectangle([-0.8, 0.1, 0], 0.7, 0.7))
    holes.append(geom.add_rectangle([-0.9, 1.0, 0], 0.8, 0.8))
    holes.append(geom.add_rectangle([-1.0, 2.0, 0], 0.9, 0.9))
    holes.append(geom.add_rectangle([-1.1, 3.1, 0], 1.0, 1.0))
    
    holes.append(geom.add_rectangle([ 0.1, 0.1, 0], 0.6, 0.6))
    holes.append(geom.add_rectangle([ 0.1, 0.9, 0], 0.5, 0.5))
    holes.append(geom.add_rectangle([ 0.1, 1.6, 0], 0.4, 0.4))
    holes.append(geom.add_rectangle([ 0.1, 2.2, 0], 0.3, 0.3))
    holes.append(geom.add_rectangle([ 0.1, 2.7, 0], 0.2, 0.2))

    keep = geom.boolean_union(rectangles)
    remove = geom.boolean_union(holes)
    
    surface = geom.boolean_difference(keep, remove)

    geom.set_recombined_surfaces(surface)
    crosshatch = geom.generate_mesh()

mesh = fe.mesh.expand(
    mesh=fe.mesh.Mesh(
        points=crosshatch.points[:, :2], 
        cells=crosshatch.cells[1][1], 
        cell_type="hexahedron"
    ),
    n=6,
    z=0.2,
)
region = fe.RegionHexahedron(mesh)
displacement = fe.Field(region, dim=3)

bounds = {"fixed": fe.dof.Boundary(displacement, fy=lambda y: y==0)}
dof0, dof1 = fe.dof.partition(displacement, bounds)

umat = fe.LinearElastic(E=210000, nu=0.3)
density = 7850 * 1e-12
gravity = np.array([0,-9.81, 0]) * 1e3

bodyforce = fe.IntegralForm(
    fun=density * gravity.reshape(-1, 1, 1),
    v=displacement,
    dV=region.dV,
).assemble()

stiffness = fe.IntegralForm(
    fun=umat.elasticity(region=region),
    v=displacement,
    dV=region.dV,
    u=displacement,
    grad_v=True,
    grad_u=True
).assemble(parallel=True)

system = fe.solve.partition(displacement, stiffness, dof1, dof0, r=-bodyforce)
displacement += fe.solve.solve(*system)

fe.save(region, displacement, filename="crosshatch_gravity.vtk")

[adtzlr]

For FElupe >=4.0.0, the code has to be adopted (introduction of field containers and list-style in- and output of forms).

import felupe as fe
import numpy as np
import pygmsh

with pygmsh.occ.Geometry() as geom:

    geom.characteristic_length_min = 0.05
    geom.characteristic_length_max = 0.05

    rectangles = []
    holes = []
    
    rectangles.append(geom.add_rectangle([-0.9, 0.0, 0], 0.9, 0.9))
    rectangles.append(geom.add_rectangle([-1.0, 0.9, 0], 1.0, 1.0))
    rectangles.append(geom.add_rectangle([-1.1, 1.9, 0], 1.1, 1.1))
    rectangles.append(geom.add_rectangle([-1.2, 3.0, 0], 1.2, 1.2))
    
    rectangles.append(geom.add_rectangle([ 0.0, 0.0, 0], 0.8, 0.8))
    rectangles.append(geom.add_rectangle([ 0.0, 0.8, 0], 0.7, 0.7))
    rectangles.append(geom.add_rectangle([ 0.0, 1.5, 0], 0.6, 0.6))
    rectangles.append(geom.add_rectangle([ 0.0, 2.1, 0], 0.5, 0.5))
    rectangles.append(geom.add_rectangle([ 0.0, 2.6, 0], 0.4, 0.4))

    holes.append(geom.add_rectangle([-0.8, 0.1, 0], 0.7, 0.7))
    holes.append(geom.add_rectangle([-0.9, 1.0, 0], 0.8, 0.8))
    holes.append(geom.add_rectangle([-1.0, 2.0, 0], 0.9, 0.9))
    holes.append(geom.add_rectangle([-1.1, 3.1, 0], 1.0, 1.0))
    
    holes.append(geom.add_rectangle([ 0.1, 0.1, 0], 0.6, 0.6))
    holes.append(geom.add_rectangle([ 0.1, 0.9, 0], 0.5, 0.5))
    holes.append(geom.add_rectangle([ 0.1, 1.6, 0], 0.4, 0.4))
    holes.append(geom.add_rectangle([ 0.1, 2.2, 0], 0.3, 0.3))
    holes.append(geom.add_rectangle([ 0.1, 2.7, 0], 0.2, 0.2))

    keep = geom.boolean_union(rectangles)
    remove = geom.boolean_union(holes)
    
    surface = geom.boolean_difference(keep, remove)

    geom.set_recombined_surfaces(surface)
    crosshatch = geom.generate_mesh()

mesh = fe.mesh.expand(
    fe.Mesh(
        points=crosshatch.points[:, :2], 
        cells=crosshatch.cells[1].data, 
        cell_type="quad"
    ),
    n=6,
    z=0.2,
)
region = fe.RegionHexahedron(mesh)
displacement = fe.Field(region, dim=3)
field = fe.FieldContainer([displacement])

bounds = {"fixed": fe.dof.Boundary(displacement, fy=lambda y: y==0)}
dof0, dof1 = fe.dof.partition(field, bounds)

umat = fe.LinearElastic(E=210000, nu=0.3)
density = 7850 * 1e-12
gravity = np.array([0,-9.81, 0]) * 1e3

bodyforce = fe.IntegralForm(
    fun=[density * gravity.reshape(-1, 1, 1)],
    v=field,
    dV=region.dV,
    grad_v=[False],
).assemble()

stiffness = fe.IntegralForm(
    fun=umat.elasticity(region=region),
    v=field,
    dV=region.dV,
    u=field,
).assemble(parallel=True)

system = fe.solve.partition(field, stiffness, dof1, dof0, r=-bodyforce)
field += fe.solve.solve(*system)

fe.save(region, field, filename="crosshatch_gravity.vtk")

[r614]

Hey! Sorry for the late reply, its finals season haha. The new tutorial and solution are amazing, thank you so much!

I'll give this another run once the next term starts, and come back with more questions (and maybe PRs, haha)