Matrix Transforms, Displacements, and Resultant Forces

[mbatesole]

Hello All,

I don't have an issue with scikit-fem so much as I have a general question of use.

I'd like to take a beam, constrain one end then displace the other end, but not just in translation-- I'd like to perform a 4x4 matrix transform on that end consisting of both rotation and translation then determine the resultant force.

Every package I've investigated so far only implements the translation without the rotation: Solidworks, Fusion 360, CalculiX, etc.

Is this something that scikit-fem can do? Would you be so kind as to show me how?

Thanks!

[gdmcbain]

So the points x on some face will be transformed to x' = A x + b for matrix A and vector b independent of x, as in https://en.wikipedia.org/wiki/Affine_transformation#Augmented_matrix ? Then the displacement is u (x) = x' − x = (A − I) x + b.

I think we can do this with a combination of

• the linear elasticity example, ex11, except that that has trivially A = I,

https://github.com/kinnala/scikit-fem/blob/de2f4acba96abd6acb716197330e24b7190618b4/docs/examples/ex11.py#L13-L19

and

• a way to have the essential boundary conditions depend on the coordinates, as in the parabolic velocity profile applied at the inlet to the backward facing step example, ex24:

https://github.com/kinnala/scikit-fem/blob/de2f4acba96abd6acb716197330e24b7190618b4/docs/examples/ex24.py#L76-L81

[gdmcbain]

I had a go at this. There was a bit more mucking about with numpy.einsum than I'd like, but I may have overlooked an easier way. It's basically ex11 but with a nonzero rotation matrix A about the z-axis.

ex_affine.zip

In the course of this I did notice some simplifications in selecting the degrees of freedom for the essential boundary condition in ex24 (& ex27) which I'll push shortly (branch essential-dofs).

[gdmcbain]

For the resultant force, one can integrate the stress over the 'left' boundary, but there might be a short-cut like (K @ u).reshape((3, -1))? The two methods are demonstrated for the scalar Laplace equation in ex13?

[mbatesole]

Wow, the transform seems to be working correctly! I've never been exposed to numpy.einsum previously. I'm still not quite sure what's going on there. That will take some more digging. But that's great!

On the other hand, I'm still completely stuck on the resultant force. I've had a go at this quite literally all day and I can't seem to figure it out. I've tried many different combinations, but nothing seems to be correct. I was referencing values in CalucliX, Solid Works, and just a plain old free body diagram to see if the results were even in the ball park. Nothing seems to match.

*RESULTS***** Attempt #1: 3.60e+04, -4.09e+03, -9.30e+00 Attempt #2: -3.02e+04, 7.23e+03, 2.29e+04

SOFTWARE: CalculiX: 3.41e-15, -8.19e-04, 2.90e-15 SolidWorks (Entire Model): 3.40e-07, -1.94e-06, 1.75e-10 SolidWorks (Left Face): 3.40e-07, 1.27e-03, 1.75e-10 SolidWorks Result (Entire Model): 1.97e-06 SolidWorks Result (Left Face): 1.27e-03

FREE BODY SkyCiv: (Beam Calculator) -8.00e-04 N

Here is a link to the script. The two models that are shown are ex11.py and one with just rotation applied so I could observe the difference.

https://pastebin.com/dYqkvdhj

[kinnala]

Sorry, but it's unclear to me what is meant by "resultant force". Integral of traction over the displaced face?

[mbatesole]

I don't have the mathematics knowledge to answer your question in like terms. But here's the idea in layman's terms:

If you have a beam of a certain dimension, of a certain material, prevent one end from moving, and displace the other end some amount in 3D space (rotation and/or translation) what is the force that is generated from that displacement?

I guess it's like FEM in reverse? Rather than input the force and see the displacement, I'd like to input the displacement, and determine the force.

Does that make any sense?

[kinnala]

Alright, I think I get your point.

So you would like to know the force (field) which must be applied to the object's free end so that the object will be deformed from the reference configuration into the configuration specified by a given rotation and translation matrix?

Edit 1: Or resulting force (field) at the fixed end?

Edit 2: Probably the latter now that I carefully read the comments of @gdmcbain .

[mbatesole]

Yes! That's it. Does it matter which end? Wouldn't the magnitudes be the same, with the force vectors inverted?

I think it's a two part problem where first we find the deformation field after applying the transform to the free end, and then calculate the forces based on the deformations and integrate for one end. But I'm honestly talking above my pay grade here.

Here's an example I found: https://primeraeng.com/wp-content/uploads/2016/03/An-Introduction-to-the-Finite-Element-Method-for-Young-Engineers_PART-2.pdf

Page 25 and 27-- It is only looking at the force (and moment) in one dimension to simplify things.

I really appreciate you guys taking the time to help me out... especially considering my lack of general knowledge of the inner workings of FEM.

[kinnala]

This should work in master branch:

from skfem.models.elasticity import linear_stress

sigma = linear_stress(1e3, 0.3)
force = [0] * 3

left_facets = m.facets_satisfying(lambda x: x[0] == 0.0)
left_basis = FacetBasis(m, e, facets=left_facets)

for i in range(3):

    @Functional
    def traction(w):
        from skfem.helpers import sym_grad, dot
        return dot(sigma(sym_grad(w['w']))[i], w.n)

    force[i] = traction.assemble(left_basis, left_basis.interpolate(u))

For me it gives (combined with the example of @gdmcbain)

In [1]: force                                                                   
Out[1]: [81.36988270476198, -0.05297977805562402, -3.686287386450715e-18]

Does it make sense?

Note: I'm currently working on a new release 1.0.0 and have almost forgotten how the things look in 0.4.2. If you want to use 0.4.2 you need to do some modifications to how the functional traction is defined.

[mbatesole]

Had another go at this today. I think I've made some progress. https://www.dropbox.com/s/uf5rie6938xbnhk/ex_affine_result_force.py?dl=0

1) I had to adjust the affine function to work with 4x4 Matrices to get the behavior I was looking for:

2) We need something to compare the result to. Firstly I used FreeCad/Calculix to translate the right side by 1mm in the y-axis. The resulting displacement mesh has the right and left sides non-orthogonal. Using only translation with a matrix transform will keep the right side parallel to the left. It seems the same thing happens in ex11.py if you modify line 19: 'u[dofs['right'].nodal['u^2']] = 0.3' https://github.com/kinnala/scikit-fem/blob/de2f4acba96abd6acb716197330e24b7190618b4/docs/examples/ex11.py#L13-L19 However this behavior is different in FreeCad/CalculiX. So I iteratively adjusted the transformation matrix to get very similar looking displaced mesh to FreeCad/CalculiX. Secondly, I calculated by hand what the result should be on the y-axis just to get a ball park figure. Here are the results:

It looks like CalculiX and Hand Calculations are in agreement but scikit-fem is off. I'm wondering if rather than integrating over the resultant forces, we're integrating over internal stresses of the material?

https://www.dropbox.com/s/uf5rie6938xbnhk/ex_affine_result_force.py?dl=0

[kinnala]

I tried to refine the mesh in your example twice and got something which is closer to the other results at least in y-axis.

**********************RESULTS***********************
scikit-fem:          -1.53e+04, -8.67e-04, -6.60e-06
CalculiX:             3.41e-15, -8.19e-04,  2.90e-15
By Hand (y-axis)                -8.00e-04

Not sure if this is what you're after. Remark: I don't completely understand why CalculiX would have both x and z components equal to zero based on my intuition. Is the translation applied differently? Maybe we are after all calculating a different quantity?

Note: I also tried to use ElementTriP2 to get a more accurate result but realized that one of the remaining issues #83 causes the approach of using get_dofs to not work for higher-order tetrahedral elements.

[kinnala]

Okay, I managed to fix #83 .

With twice refined tetrahedral mesh and quadratic elements I get

**********************RESULTS***********************
scikit-fem:          -1.28e+04, -2.90e-04, -4.31e-06
CalculiX:             3.41e-15, -8.19e-04,  2.90e-15
By Hand (y-axis)                -8.00e-04

The order of magnitude seems to be still correct. This is most likely more accurate than the linear solution above.

In my opinion it makes sense that

• z-component is almost zero since no twisting happens in that direction
• x-component is nonzero since there will always be some pull or push if the beam is displaced in this way, and
• y-component has the same order of magnitude than your hand calculation

I'm wondering if rather than integrating over the resultant forces, we're integrating over internal stresses of the material?

Stress tensor (= "internal stress") multiplied by a unit vector gives a force (field) called "traction". Traction on the boundary of the domain is something that can be specified as a boundary condition and is equivalent to an external force (field). Here we do not specify the traction as a boundary condition but simply integrate the resulting traction over the left boundary:

@Functional
def traction(w):
    from skfem.helpers import sym_grad, dot
    return dot(sigma(sym_grad(w['w']))[i], w.n)

force[i] = traction.assemble(left_basis, left_basis.interpolate(u))

Here sigma(...) is the stress tensor and w.n is the normal unit vector. The last line integrates over the left boundary. So we are only considering boundary stresses and not internal stresses.

[mbatesole]

I appreciate you explaining all of this to me. I'm used to just hitting the big gray button that looks like this: [SOLVE]. So I'm learning a lot and really appreciate your patience. :D

Would you mind popping over here and seeing how they determine their reactions forces? https://fenicsproject.discourse.group/t/computation-of-reaction-force/129

Are we doing it the same way?

[kinnala]

That seems to be a slightly different approach, let me try it out.

[gdmcbain]

Thanks for the link to the 'Computation of reaction force'. That looks like a domain integral rather than a boundary integral. I think this is like what I was alluding to on Friday. Often one can swap between domain integrals and boundary integrals using some version of 'integrating by parts' or the divergence theorem (a.k.a. Gauss's or Ostrogradsky's, Green's, or Stokes's theorem, ...). I'm very familiar with this in heat transfer and incompressible fluid flow, but a little less sure exactly how it goes in elasticity and I don't have my library here at the moment.

I think in elasticity the equivalence can be discussed in terms of energy and work. The increase in the elastic energy stored in the elastic domain equals the work done by the external forces. The external forces here are the traction applied at the 'right' to produce the imposed displacement and the traction applied at the 'left' to keep the displacement there zero. The work done by a traction is the integral of the scalar product of the traction vector with the displacement vector. (The traction is in turn the scalar product of the stress tensor with the normal vector of the surface on which it acts.)

The elastic energy is something like the integral over the domain of the scalar product (double contraction) of the stress and strain tensors; I think that that's what's meant by the FEniCS expression

assemble(inner(sigma(u,p),epsilon(v))*dx)

A good way to derive these relations begins with the weak formulation of the problem, replacing the test function with the solution. But this may be going into more detail that the OP is interested in at the moment. This is what is used for the 'Laplace with mixed boundary conditions' ex13 where the conductance is computed from the domain integral of Joule heating

u @ A @ u

and the current from a boundary integral of flux

for port, boundary in mesh.boundaries.items():
    fbasis = FacetBasis(mesh, elements, facets=boundary)
    current[port] = asm(port_flux, fbasis, w=fbasis.interpolate(u))

and they give numerically the same answer (or equal & opposite for the currents at the two ends):

conductance: {'skfem': 0.4416548065660517, 'exact': 0.4412712003053032}
Current in through ports: {'ground': -0.44160370313261493, 'positive': 0.4416152264597093}

The same kind of thinking applies to continuum mechanics, it's just more involved because this involves vectors and tensors instead of scalars and vectors.

Although mathematically equivalent for the continuous problem, for the discrete problem, the domain approach is often more accurate and easier to implement, as argued in 'Three ways to compute multiport inertance'.

[gdmcbain]

Another difference following from the higher tensorial order of elastic deformation versus potential flow is that the 'Resultant force' also includes a torque.

A resultant force is the single force and associated torque obtained by combining a system of forces and torques acting on a rigid body.

In lumped models, this is the difference between 'pinned' and 'clamped' joints.

I don't know whether the torque is of interest here, but if this example evolves into one for the documentation, it might be a good idea to calculate the resultant torque too.

Actually there are a couple of concepts overlapping here: 'resultant' and 'reaction'. The system of tractions appears at the 'left' in reaction to the action of the displacements imposed at 'right'. This system can then be integrated to give the 'resultant' force (& torque).

[kinnala]

Oh crap, there were some serious mistakes in my earlier calculations. I accidentally... 1) ...initialized sigma with wrong material parameters 2) ...used the deformed mesh when postprocessing Now I'm at least getting consistent results with the two approaches that I'm using:

**********************RESULTS***********************
scikit-fem 2:            -2.72e+04, -1.15e+03, 2.87e+00
scikit-fem:          -3.40e+04, -3.25e+03, -2.56e+01
CalculiX:             3.41e-15, -8.19e-04,  2.90e-15
By Hand (y-axis)                -8.00e-04

Here scikit-fem 2 is the domain integral approach and scikit-fem is the boundary integral approach from before.

Now the question remains that how can we compare to the results of @mbatesole . Now I think since you seem to use mm as a unit of length (based on the comments of the code), you should use MPa as a unit of Young's modulus. This gives much better correspondence:

**********************RESULTS***********************
scikit-fem 2:            -2.72e-02, -1.15e-03, 2.87e-06
scikit-fem:          -3.40e-02, -3.25e-03, -2.56e-05
CalculiX:             3.41e-15, -8.19e-04,  2.90e-15
By Hand (y-axis)                -8.00e-04

This is still with linear elements. To get closer to numerical convergence, I'm refining the mesh and using second order tetrahedrons.

**********************RESULTS***********************
scikit-fem 2:            -2.73e-02, -1.01e-03, -3.42e-08
scikit-fem:          -2.86e-02, -1.21e-03, -1.68e-05
CalculiX:             3.41e-15, -8.19e-04,  2.90e-15
By Hand (y-axis)                -8.00e-04

Theoretically speaking approach 2 is more accurate (as @gdmcbain mentioned above), so I think this is as close as I can get now with a laptop and this numerical approach. Here is the code:

import pygmsh
import numpy as np
from skfem import *
from skfem.models.elasticity import linear_elasticity, lame_parameters
from skfem.io.meshio import from_meshio
import transformations as tf

BEAM_WIDTH = .635
BEAM_HEIGHT = 0.4318 # .017in=0.4318mm;  .019in=0.4826mm
BEAM_LENGTH = 4

def create_beam():
    geom = pygmsh.opencascade.Geometry(characteristic_length_min=.3, characteristic_length_max=.3 )
    rectangle = geom.add_rectangle([0, 0, 0], BEAM_LENGTH, BEAM_HEIGHT)
    geom.extrude(rectangle, [0, 0, BEAM_WIDTH])

    mesh = pygmsh.generate_mesh(geom)
    return mesh

def affine(x, y, z):
    result =  np.einsum('ijk->kij', np.einsum('ij,j...', a - np.eye(a.shape[0]), np.stack([x, y, z, np.ones(x.shape)])) + b)
    return result[:3,:]

mesh = create_beam()
m = from_meshio(mesh)
m.refine(2)

# offset mesh to get origin in the middle of the left side
m.p[1] -= BEAM_HEIGHT/2.0
m.p[2] -= BEAM_WIDTH/2.0

# Stiffness Matrix
e1 = ElementTetP2()
e = ElementVectorH1(e1)
ib = InteriorBasis(m, e)
K = asm(linear_elasticity(*lame_parameters(4, 0.4)), ib)

dofs = {
    'left': ib.get_dofs(lambda x: x[0] == 0.0),
    'right': ib.get_dofs(lambda x: x[0] == BEAM_LENGTH),}

u = np.zeros(K.shape[0])
right_facets = m.facets_satisfying(lambda x: x[0] == BEAM_LENGTH)
right_basis = FacetBasis(m, e, facets=right_facets)
right_dofs = right_basis.get_dofs(right_facets).all()

# Transformation Matrix
# Here I've iteratively adjusted the transformation until the mesh displacement matches very
# closely the output of CalucliX for a displacement on the y-axis of 1.0 mm
a = tf.rotation_matrix(np.radians(20), [0, 0, 1], [BEAM_LENGTH, 0, 0])
b = [.1, 1, 0, 1]

# Solve the FEA for displacement
u[right_dofs] = L2_projection(affine, right_basis, right_dofs)

I = ib.complement_dofs(dofs)
u = solve(*condense(K, 0 * u, I=I, x=u))

# Determine the resultant forces. Approach 1.
from skfem.models.elasticity import linear_stress

sigma = linear_stress(*lame_parameters(4, 0.4))
force1 = [0] * 3

left_facets = m.facets_satisfying(lambda x: x[0] == 0.0)
left_basis = FacetBasis(m, e, facets=left_facets)

for i in range(3):

    @Functional
    def traction(w):
        from skfem.helpers import sym_grad, dot
        return dot(sigma(sym_grad(w['w']))[i], w.n)

    force1[i] = traction.assemble(left_basis, left_basis.interpolate(u))

# Determine the resultant forces. Approach 2.
y = K @ u
force2 = [0] * 3

for i in range(3):
    ido = np.zeros(u.shape[0])
    dimdofs = ib.nodal_dofs[i] if len(ib.edge_dofs) == 0 else\
        np.union1d(ib.nodal_dofs[i], ib.edge_dofs[i])
    ido[dimdofs] = 1
    ido[dofs['right'].all()] = 0
    force2[i] = ido @ y

if __name__ == "__main__":
    print()
    print("**********************RESULTS***********************")
    print("scikit-fem 2:\t\t\t {:.2e}, {:.2e}, {:.2e}".format(force2[0], force2[1], force2[2]))
    print("scikit-fem:\t\t\t {:.2e}, {:.2e}, {:.2e}".format(force1[0], force1[1], force1[2]))
    print("CalculiX:\t\t\t  {:.2e}, {:.2e},  {:.2e}".format(3.408515e-15, -8.185536e-04, 2.900187e-15))
    print("By Hand (y-axis)\t\t\t    {:.2e}".format( -.0008))

    from os.path import splitext
    from sys import argv

    m.save(splitext(argv[0])[0] + '.vtk')

In practice I'd probably use hexahedral elements or at least a tensor product mesh because of the apparent symmetries in the problem.

[mbatesole]

Wow thank you guys so much!

@gdmcbain: Your explanations taught me a lot! Personally, I'd love to see the resultant force broken down into forces and torques.

@kinnala: Here's one thing that's super interesting to me... going back to https://github.com/kinnala/scikit-fem/issues/357#issuecomment-608976910 I listed the values I got running the simulation through SolidWorks.

SolidWorks (Left Face): 3.40e-07, 1.27e-03, `1.75e-10

As you can see the value for the y-axis is pretty much spot on with your latest results. So now I'm wondering what is going on with the Hand and CalculiX measurements?!

Running the example with second order tetrahedrons causes this error on my machine:

/usr/local/lib/python3.6/dist-packages/scikits/umfpack/umfpack.py:563: UmfpackWarning: Singular matrix
  warnings.warn('Singular matrix', UmfpackWarning)
/usr/local/lib/python3.6/dist-packages/scikits/umfpack/umfpack.py:712: UmfpackWarning: Zeroing nan and inf entries...
  warnings.warn('Zeroing nan and inf entries...', UmfpackWarning)
/usr/local/lib/python3.6/dist-packages/scikits/umfpack/umfpack.py:717: RuntimeWarning: divide by zero encountered in double_scalars
  econd = 1.0 / self.info[UMFPACK_RCOND]
/usr/local/lib/python3.6/dist-packages/scikits/umfpack/umfpack.py:721: UmfpackWarning: (almost) singular matrix! (estimated cond. number: inf)
  warnings.warn(msg, UmfpackWarning)

Looks like the error happens during the call to L2_projection.

For fun I attempted to run the example with a hexahedral mesh as you suggested in https://github.com/kinnala/scikit-fem/issues/357#issuecomment-609729303

**********************RESULTS***********************
scikit-fem 2:            -2.70e-02, -1.24e-03, -1.83e-15
scikit-fem:          -4.26e-02, -1.01e-02, -2.09e-15
CalculiX:             3.41e-15, -8.19e-04,  2.90e-15
By Hand (y-axis)                -8.00e-04

It agrees with the previous data.

I'm curios to know why the blue on the right end doesn't extend to the edges of the face in the y-dimension. Is that what you'd expect, or did I set something up incorrectly?

For sake of completeness, here's the same setup but with the linear tetrahedral mesh.

And, here's my code to generate the hexahedral mesh:

def create_hex_beam():
    geom = pygmsh.built_in.Geometry()
    geom.add_raw_code('Mesh.RecombineAll = 1;')
    geom.add_raw_code('Mesh.Recombine3DAll = 1;')

    #Algorithm: 1:MeshAdapt, 2:Automatic, 3:Delaunay, 6:Frontal-Delaunay,
    #7:BAMG, 8:Frontal-Delaunay for Quads, 9:Packing of Parallelograms
    #geom.add_raw_code('Mesh.Algorithm = 8;')

    # Subdivision: 0:None, 1:all quadrangles, 2:all hexadedra, 3:barycentric
    geom.add_raw_code('Mesh.SubdivisionAlgorithm = 2;')
    xmin = 0
    xmax = BEAM_LENGTH
    ymin = -BEAM_HEIGHT/2.0
    ymax = BEAM_HEIGHT/2.0

    rectangle = geom.add_rectangle(xmin, xmax, ymin, ymax,z=0)
    geom.set_transfinite_lines(rectangle.lines, size=8)#, bump=.2)
    geom.set_transfinite_surface(rectangle.surface)
    geom.extrude(rectangle,[0, 0, BEAM_WIDTH], num_layers=4, recombine=True )

    mesh = pygmsh.generate_mesh(geom)
    m = from_meshio(mesh)
    m.p[2] -= BEAM_WIDTH / 2.0
    e = ElementHex1()
    return(m, e)

[gdmcbain]

The hexahedral mesh should be easy to make with skfem.MeshHex.init_tensor. (Although it'd be really nice to have ElementHex2 #327 or a 20-node serendipity element. TODO!)

I was also thinking it should be possible to simply load the Calculix mesh directly into scikit-fem, thanks to meshio.

I did have Calculix installed on my Ubuntu machine in the office and will look into setting it up on this Windows 10 laptop later on. I find these comparisons between codes extremely helpful in general.

[gdmcbain]

The domain integral approach didn't come out as I had imagined. (Not that everything every comes out as I had imagined!) I'll work through this more closely. My idea is to repeat the analysis of 'Three ways…' for elasticity instead of potential flow. It will be good to have this written up anyway.

[gdmcbain]

Ah, I see, yes, the domain integral approach isn't so straightforward here unless the displacement imposed over the 'right' is uniform; otherwise it doesn't factor out. (The analysis might be generalized but...that remains TODO.)

Here's my sketch of a derivation.

In the absence of body forces, equilibrium is

div σ = 0

where σ is the stress tensor.

The weak formulation is that the inner product (integral over domain of double contraction) of this with any test vector field v vanishes.

(v, div σ) = 0, for all v

Integrate by parts, using symmetry of σ,

(v, n . σ)_boundary = (sym grad v, σ)

(Have I dropped a factor of a half on the right-hand side? There are lots of factors of a half in the theory of elasticity.)

The boundary integral can be split up into the sum of:

• over the free boundary where n . div σ = 0, which vanishes
• over the parts with prescribed traction n . div σ, but that's nowhere here
• over the parts with prescribed displacement, here 'right' and 'left'

• over parts with other (e.g. 'Robin') boundary conditions but none of that here

  (v, n . σ)_right + (v, n . σ)_left = (sym grad v, σ)

Put v = u, the displacement, noting:

• sym grad u = ε, the strain. (Or is it half or twice the strain?)

• u = 0 on 'left'

  (u, n . σ)_right = (ε, σ).

Now the difficulty appears: if u were uniform, it could be factored outside the inner product on the left leaving an expression for the resultant force in terms of the elastic strain energy (or half or twice it) on the right, e.g. if u = u0 i

(i, n . σ)_right = (ε, σ) / u0

and the left-hand side is i . F on 'right', but in the OP u varies over the 'right' and more care is required, or maybe it's not worthwhile.

Apologies for carelessness with regard to factors of a half and doubtless idiosyncratic notation.

[gdmcbain]

I think extending the previous to torque would begin with multiplying the partial differential equation by the position vector x relative to the point about which the resultant torque is reckoned. Is it a cross-product, x × div σ = 0?

[gdmcbain]

I'm not understanding something simple. In

# Transformation Matrix
# Here I've iteratively adjusted the transformation until the mesh displacement matches very
# closely the output of CalucliX for a displacement on the y-axis of 1.0 mm
a = tf.rotation_matrix(np.radians(20), [0, 0, 1], [BEAM_LENGTH, 0, 0])
b = [.1, 1, 0, 1]

Why isn't 'a displacement on the y-axis of 1.0 mm' just represented by A = 0 and b = [0, 1, 0]? (And why does the b above have four components? Has it been redefined from u = x'- x and x'= Ax + b?)

[gdmcbain]

Hey vtkplotter.Plotter looks great! Three-dimensional plotting in matplotlib is awful; it distorts the three axes uncontrollably #44. ParaView, Gmsh, & VisIt work but are too heavy; this kind of thing (from ex_affine_result_force.py)

 from vtkplotter import Plotter

    m.save(splitext(argv[0])[0] + '.vtk')

    vp = Plotter(bg=np.array([.223,.223,.223]),axes=2)
    beam1 = vp.load(splitext(argv[0])[0] + '.vtk')
    beam1.color('steelblue')
    vp.show([beam1])

is exactly what I'd like to be able to do in scikit-fem scripts.

[gdmcbain]

Could you share the CalculiX script? I think I've installed CalculiX 2.10 O. K. from https://github.com/GeneralElectric/CalculiX but am not all familiar with the input syntax.

[mbatesole]

vtkplotter really is great! I don't know how much you've played with it, but you can interact with it as well. If you hit the h key when the window is active it'll print() you a list of commands you can use. Even better, the developer has so many examples showing off its capabilities (including user interface elements). I was going to actually open a new issue here and suggest you guys take a look at it. I find it MUCH better than matplotlib. https://github.com/marcomusy/vtkplotter

I'd love it if you were able to use this to publish something!

Here's the CalculiX .inp generated by FreeCad/Fem: https://www.dropbox.com/s/wl85ggzerpzuo8c/FEMMeshGmsh.inp?dl=0

[mbatesole]

Oh about the matrix transforms. I don't know the mathematical term, but the order of operation is important. (My formal training is in orthodontics by the way so if you want to use big words about teeth, let me know ;) ) Multiplying in one direction give one result, in the other it gives another result. The cool thing about matrix transforms is you can combine a bunch of matrices and get a final transform. With a 3x3 matrix the rotation happens from the origin, but if you extend that to 4x4 you can rotate around an arbitrary point in space (the fourth column).

So in this example, I wanted to first rotate the right face around it's center, then translate it in space. Rotating the face around the origin, then translating it gives a completely different spatial result.

[gdmcbain]

I don't know the mathematical term, but the order of operation is important.

This could be ‘associative’ or ‘commutative’. Transforms A, B, & C are associative if (AB)Cx = A(BC)x for all x; transforms A & B commute if A(Bx) = B(Ax) for all x. Our transforms are associative, so the order or combining (AB)C or A(BC) isn't important but not commutative so AB≠BA, in general.

The biggest word I know in orthodontics is 'OW!'; feel free to drop more in.

But what I was confused with above, with why b was a 4-vector instead of a 3-vector, was when earlier I had looked up Affine_transformation#Augmented_matrix, I had found that

to perform a 4x4 matrix transform on that end consisting of both rotation and translation

it had the form

where y, x, and b are 3-vectors, A is 3x3.

I think now on looking at it again, @kinnala had just tacked the 1 that's always there in the lower-right hand corner of the 4×4 matrix onto the end of the 3-vector b to make it a 4-vector, so it's just a change of notation.

Going back to

So in this example, I wanted to first rotate the right face around it's center, then translate it in space. Rotating the face around the origin, then translating it gives a completely different spatial result.

The effect of two transformations, x1 = A0 x0 + b0, x2 = A1 x1 + b1 = A1 (A0 x0 + b0) + b1, when the first is a rotation about the origin (b0 = 0) and the second a translation (A1 = I) is x2 = A0 x0 + b1, which is just one 4×4 transformation with 3×3 rotation matrix A0 and translation 3-vector b1. Translating about the centre would be more involved; we'd need know where the centre was. So working through the interpretation of rotation about an arbitrary point.

With a 3x3 matrix the rotation happens from the origin, but if you extend that to 4x4 you can rotate around an arbitrary point in space (the fourth column).

If the centre of rotation is c, I think we need to:

1. Translate centre to origin (A=I, b= −c)
2. Rotate about origin (A=A, b=0)
3. Translate back (A=I, b=+c)

That's (N. B.: sequence of transformations is right to left)

[ I | c ] [ A | 0 ] [ I | –c ] = [ A | c – Ac ]
[ 0 | 1 ] [ 0 | 1 ] [ 0 |  1 ]   [ 0 |      1 ]

i.e. a 4×4 transformation with rotation matrix A and translation b = c − Ac. Is that right?

[gdmcbain]

Thank you for the FEMMeshGmsh.inp. It doesn't load into meshio:

meshio.read('FEMMeshGmsh.inp')

raises a ValueError. I might investigate this further over at meshio, but it's not immediately a concern for scikit-fem. It would have been nice to have been able to load the mesh, but no matter.

The file does run through ccx_2.16 which I have installed on Ubuntu 19.10.

Looking into FEMMeshGmsh.inp in Emacs, I see NSET=ConstraintFixed (line 786) and NSET=ConstraintDisplacement (line 825) which I presume correspond to 'left' and 'right' and then on line 897:

ConstraintDisplacement,2,2,1.0

I think this means that from the second to the second component of displacement is 1.0. This doesn't look like what has been applied in the pictures above. I would have thought that this would constrain NSET=ConstraintFixed and NSET=ConstraintDisplacement to remain parallel, just with the latter translated in its plane in the y-direction.

So far I've only run FEMeshGmsh.inf in ccx_2.16 which has produced a bunch of files (.12d, .cvg, .dat, .frd, .sta). (I haven't worked out how to render these graphically yet.) FEMMeshGmsh.dat says that the 'total force' for 'set CONSTRAINTFIXED' is

        3.408515E-15 -8.185536E-04  2.900187E-15

so essentially all in the y-direction.

[gdmcbain]

I reported the meshio.read error upstream at nschloe/meshio#783.

On looking more closely into FEMMeshGmsh.inp, I don't think we're going to be able to read it here anyway: TYPE=C3D10 corresponds to meshio type "tetra10"; scikit-fem only reads simplicial meshes like "tetra4".

[kinnala]

It should be in principle be possible to read second-order meshes now that we have the option to use curved meshes (in ex31), but I haven't tested it and there is no simple interface for that.

[marcomusy]

Thanks guys for the nice words :) i'm happy that you found it useful. vtkplotter has now changed name to vedo but the functionality is the same. Recent developments (still abit experimental) which might be of interest to you are visualization of tet meshes, including cutting/slicing/thresholding/isosurfacing e.g.

more examples here

[kinnala]

@marcomusy Thanks for reaching out. We are looking into a way to make a simple wrapper for transforming meshes and solution fields from scikit-fem to vtkplotter (now vedo) in https://github.com/kinnala/scikit-fem/issues/361. Our goal is to add a simple module to skfem.visuals which figures out the nasty details on how to transform between the data structures of the two libraries. In scikit-fem we have classes MeshTet and MeshHex that describe tetrahedral and hexahedral meshes.

I could quite easily adapt this example to read scikit-fem mesh into vedo: https://github.com/marcomusy/vedo/blob/master/vedo/examples/tetmesh/tet_build.py Basically we have points and tets as in the above example available as MeshTet.p and MeshTet.t. The solutions are given also as arrays similar to scal.

Do you find this the correct interface to use for this kind of use case or is there something else we should use?

[marcomusy]

@kinnala I think that's exactly the way to go! As you already have got everything in class MeshTet it should be quite straightforward to connect the two libraries.

btw i was also planning to add a vedo.HexMesh class so you might exploit that too in sciki-fem

[kinnala]

@kinnala I think that's exactly the way to go! As you already have got everything in class MeshTet it should be quite straightforward to connect the two libraries.

btw i was also planning to add a vedo.HexMesh class so you might exploit that too in sciki-fem

Thanks for the answer. We continue tracking this in #361 .