Using a custom user element (UEL)

[bruno.zilli]

GitHub - brunozilli/beam_shear_center_test_upn100: An academic solver for...

Vlasov

Kind Regards

Bruno Zilli

[Disla]

Hello Bruno,

Thank you very much!!. I have a more elaborated section which you could add to your tests collection. Node coordinates are given in mm.

Shear_Centre_Mesh_Test1.inp

Shear center was computed as.

X  0 mm
Y  -90.28 mm
Z  -73.61 mm

Coordinates origin is :

[bruno.zilli]

Ciao Disla, and thank you so much for the kind words and support!

Yes, absolutely – I’ll get on it as soon as I can. In the meantime, would it be possible for you to just share the raw geometry file? Any format like BREP, STEP, or IGES would be perfect.

My current plan is to build the 1D midline mesh (strictly Vlasov style!) using SALOME, since I already have a custom parser working nicely with its UNV flavour. A future step will definitely be to extend the parser to handle Gmsh .msh files as well, so stay tuned!

In the meanwhile… would you be curious to try a test run yourself? All you need is a 1D mesh of just the section’s border, centred as neatly as you can on the centroid. If you want to double-check the centroid location beforehand, you can use a quick 2D mesh and the subroutine in my Poisson repo right here:
compute_section_properties.f

I must confess, I’ve found an absolutely brilliant assistant in DeepSeek – I highly recommend letting her lend you a hand with the code, she’s been marvellous!

Grazie mille ancora, and looking forward to seeing what you come up with! A presto,
Bruno

[Disla]

@bruno.zilli

would you be curious to try a test run yourself?

Unfortunately my programing skills are very limited but I have manage to reconstruct the original Step from the mesh.

Shear_Centre_Section_contour_faceted.step

[bruno.zilli]

ciao Disla, I:m on a business trip. I’ll get it on Monday. Thanks. Kind regards. Bruno

[bruno.zilli]

ciao Disla, your link didn’t work. No file step-iges. Anyway, here below is a working path most suitable for a general beam 3D, Please, everybody is very welcome to participate in this work.

┌─────────────────────────────────────────────────────────────────────────────────┐
│ FILES STRUCTURE FOR CALCULIX UEL │
└─────────────────────────────────────────────────────────────────────────────────┘

                ┌─────────────────────────────────────┐
                │  EXISTING CALCULIX FILES            │
                │  (DO NOT MODIFY - already compiled) │
                │                                     │
                │  • ccx_2.20/src/                    │
                │    - mainsp.f                       │
                │    - results.f                      │
                │    - spooles_driver.f                │
                │    - ...                            │
                └─────────────────────────────────────┘
                                  │
                                  │  CALLS
                                  ▼

┌─────────────────────────────────────────────────────────────────────────────────┐
│ NEW FILES (YOU CREATE) │
│ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 1: e_c3d_u100.f │ │
│ │ └── Main UEL subroutine called by CalculiX during assembly │ │
│ │ Computes: amatrx(18,18) and rhs(18) │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │ │
│ │ CALLS │
│ ▼ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 2: beam3d_k.f │ │
│ │ └── Computes local stiffness matrix Kel(18,18) │ │
│ │ Uses 3-point Gauss integration │ │
│ │ Calls: build_D, build_B, shape3 │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │ │
│ │ CALLS │
│ ▼ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 3: build_D.f │ │
│ │ └── Builds 6x6 constitutive matrix D │ │
│ │ Includes shear center offset (ey, ez) │ │
│ │ WARNING: condensed torsion valid only for compact sections │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 4: build_B.f │ │
│ │ └── Builds 6x18 strain-displacement matrix B │ │
│ │ Maps DOFs to strains: ε, γy, γz, κy, κz, ω │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 5: shape3.f │ │
│ │ └── Quadratic shape functions for 3-node line element │ │
│ │ Natural coordinate ξ ∈ [-1, 1] │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 6: transform_beam18.f │ │
│ │ └── Builds 18x18 rotation matrix T (local → global) │ │
│ │ T = blkdiag(R, R, R, R, R, R) │ │
│ │ R(3,3) from ex, ey, ez axes │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 7: matmul.f │ │
│ │ └── Utility: matrix multiplication (18x18) and A^T * B │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │
│ ┌─────────────────────────────────────────────────────────────────────────┐ │
│ │ FILE 8: resultsmech_u100.f │ │
│ │ └── Post-processing UEL called by CalculiX after solution │ │
│ │ Computes stresses/strains for .dat and .frd output │ │
│ └─────────────────────────────────────────────────────────────────────────┘ │
│ │
└─────────────────────────────────────────────────────────────────────────────────┘

2. DETAILED SUBROUTINE SPECIFICATIONS

FILE 1: e_c3d_u100.f (MAIN UEL)

Item	Description
Purpose	Main UEL entry point called by CalculiX during global matrix assembly
Called by	CalculiX mainsp.f (automatically when element type U100 is used)
Calls	beam3d_k, transform_beam18, matmul_transpose, matmul
Inputs	xl(3,3) - nodal coordinates (global)
	elvar(18) - current nodal displacements (global)
	props(12) - element properties from *UEL PROPERTY
	iel - element ID (for debugging)
	mi(*) - model information
Outputs	amatrx(18,18) - element stiffness matrix (global)
	rhs(18) - element residual vector
Modifications needed	None - create from scratch

What this subroutine does (step by step):

fortran

SUBROUTINE e_c3d_u100(rhs, amatrx, ..., xl, elvar, props, ...)

STEP 1: Extract inputs
    • x1,y1,z1 = xl(1,1), xl(2,1), xl(3,1)  ! node 1
    • x2,y2,z2 = xl(1,2), xl(2,2), xl(3,2)  ! node 2
    • x3,y3,z3 = xl(1,3), xl(2,3), xl(3,3)  ! node 3
    • E, G, A, Iy, Iz, J, ksy, ksz, yg, zg, ys, zs = props(1..12)
    • L = sqrt((x3-x1)^2 + (y3-y1)^2 + (z3-z1)^2)
    • ey = ys - yg
    • ez = zs - zg

STEP 2: Compute local stiffness matrix
    CALL beam3d_k(E, G, A, Iy, Iz, J, ksy, ksz, ey, ez, L, Kel)

STEP 3: Build transformation matrix
    CALL transform_beam18(x1, y1, z1, x3, y3, z3, T)

STEP 4: Transform to global coordinates
    CALL matmul_transpose(T, Kel, dummy)   ! dummy = T^T * Kel
    CALL matmul(dummy, T, Kglob)           ! Kglob = T^T * Kel * T

STEP 5: Set outputs
    amatrx(1..18, 1..18) = Kglob(1..18, 1..18)
    DO i = 1, 18
        rhs(i) = 0.0
        DO j = 1, 18
            rhs(i) = rhs(i) - Kglob(i,j) * elvar(j)
        END DO
    END DO

RETURN

FILE 2: beam3d_k.f (LOCAL STIFFNESS MATRIX)

Item	Description
Purpose	Compute local stiffness matrix (18x18) for 3D Timoshenko beam with shear center offset
Called by	e_c3d_u100
Calls	build_D, shape3, build_B
Inputs	E, G, A, Iy, Iz, J, ksy, ksz, ey, ez, L
Outputs	ke(18,18) - local stiffness matrix
Modifications needed	None - create from scratch

What this subroutine does (step by step):

fortran

SUBROUTINE beam3d_k(E, G, A, Iy, Iz, J, ksy, ksz, ey, ez, L, ke)

STEP 1: Initialize ke(18,18) = 0

STEP 2: Build constitutive matrix D(6,6)
    CALL build_D(E, G, A, Iy, Iz, J, ksy, ksz, ey, ez, D)

STEP 3: Gauss integration loop (3 points)
    DO g = 1, 3
        ξ = [-0.7745966692, 0.0, 0.7745966692]
        w = [5/9, 8/9, 5/9]
        
        CALL shape3(ξ, N, dNdxsi)
        
        detJ = L / 2.0
        dNdx(1..3) = dNdxsi(1..3) / detJ
        
        CALL build_B(dNdx, N, B)
        
        Bt = transpose(B)
        
        weight = detJ * w(g)
        
        ! FULL matrix multiplication (preserves coupling!)
        DO i = 1, 18
            DO j = 1, 18
                DO k = 1, 6
                    DO l = 1, 6
                        ke(i,j) = ke(i,j) + Bt(i,k) * D(k,l) * B(l,j) * weight
                    END DO
                END DO
            END DO
        END DO
    END DO

RETURN

FILE 3: build_D.f (CONSTITUTIVE MATRIX)

Item	Description
Purpose	Build 6x6 constitutive matrix with shear center offset
Called by	beam3d_k
Calls	None
Inputs	E, G, A, Iy, Iz, J, ksy, ksz, ey, ez
Outputs	D(6,6) - constitutive matrix
Modifications needed	None - create from scratch

Matrix structure:

text

D(6,6) = [  EA      0        0        0        0        0     ]
         [   0     GAy       0        0        0     GAy·ez  ]
         [   0      0      GAz       0        0    -GAz·ey  ]
         [   0      0        0      EIz       0     EIz·ey  ]
         [   0      0        0        0      EIy   -EIy·ez  ]
         [   0    GAy·ez  -GAz·ey   EIz·ey  -EIy·ez  GJ+term]

where:
    GAy = G * A * ksy
    GAz = G * A * ksz
    EIy = E * Iy
    EIz = E * Iz
    GJ  = G * J
    term = GAy·ez² + GAz·ey² + EIz·ey² + EIy·ez²

IMPORTANT WARNING: The term added to D(6,6) is valid ONLY for compact sections. For open thin-walled sections (UPN, L, C, I-beams), this OVERESTIMATES torsional stiffness. Use Vlasov beam (7 DOF) instead.

FILE 4: build_B.f (STRAIN-DISPLACEMENT MATRIX)

Item	Description
Purpose	Build 6x18 strain-displacement matrix B
Called by	beam3d_k
Calls	None
Inputs	dNdx(3) - derivatives of shape functions w.r.t. x
	N(3) - shape functions
Outputs	B(6,18) - strain-displacement matrix
Modifications needed	None - create from scratch

Matrix structure (non-zero entries):

text

Row 1 (ε = du/dx):
    B(1,1)  = dN1    ! node 1, DOF u
    B(1,7)  = dN2    ! node 2, DOF u
    B(1,13) = dN3    ! node 3, DOF u

Row 2 (γy = dv/dx - θz):
    B(2,2)  = dN1    ! node 1, DOF v
    B(2,8)  = dN2    ! node 2, DOF v
    B(2,14) = dN3    ! node 3, DOF v
    B(2,6)  = -N1    ! node 1, DOF θz
    B(2,12) = -N2    ! node 2, DOF θz
    B(2,18) = -N3    ! node 3, DOF θz

Row 3 (γz = dw/dx + θy):
    B(3,3)  = dN1    ! node 1, DOF w
    B(3,9)  = dN2    ! node 2, DOF w
    B(3,15) = dN3    ! node 3, DOF w
    B(3,5)  = N1     ! node 1, DOF θy
    B(3,11) = N2     ! node 2, DOF θy
    B(3,17) = N3     ! node 3, DOF θy

Row 4 (κy = dθy/dx):
    B(4,5)  = dN1    ! node 1, DOF θy
    B(4,11) = dN2    ! node 2, DOF θy
    B(4,17) = dN3    ! node 3, DOF θy

Row 5 (κz = dθz/dx):
    B(5,6)  = dN1    ! node 1, DOF θz
    B(5,12) = dN2    ! node 2, DOF θz
    B(5,18) = dN3    ! node 3, DOF θz

Row 6 (ω = dθx/dx):
    B(6,4)  = dN1    ! node 1, DOF θx
    B(6,10) = dN2    ! node 2, DOF θx
    B(6,16) = dN3    ! node 3, DOF θx

FILE 5: shape3.f (SHAPE FUNCTIONS)

Item	Description
Purpose	Compute quadratic shape functions and derivatives for 3-node line element
Called by	beam3d_k
Calls	None
Inputs	xi - natural coordinate (-1 to +1)
Outputs	N(3) - shape functions
	dNdxsi(3) - derivatives w.r.t. natural coordinate ξ
Modifications needed	None - create from scratch

Formulas:

text

N(1)   = -0.5 * xi * (1.0 - xi)
N(2)   = 1.0 - xi**2
N(3)   = 0.5 * xi * (1.0 + xi)

dNdxsi(1) = -0.5 + xi
dNdxsi(2) = -2.0 * xi
dNdxsi(3) = 0.5 + xi

FILE 6: transform_beam18.f (ROTATION MATRIX)

Item	Description
Purpose	Build 18x18 transformation matrix from local to global coordinates
Called by	e_c3d_u100
Calls	None
Inputs	x1,y1,z1 - coordinates of node 1
	x3,y3,z3 - coordinates of node 3
Outputs	T(18,18) - rotation matrix
Modifications needed	None - create from scratch

Axis construction:

text

ex = (x3-x1, y3-y1, z3-z1) / L    (local x-axis, from node1 to node3)

Reference vector v:
    v = (0,0,1) normally
    if ex nearly parallel to (0,0,1): v = (0,1,0)

ez = (ex × v) / |ex × v|           (local z-axis)

ey = ez × ex                        (local y-axis)

R(3,3) = [ex; ey; ez]              (rotation matrix)

T(18,18) = blkdiag(R, R, R, R, R, R)

FILE 7: matmul.f (UTILITY)

Item	Description
Purpose	Matrix multiplication utilities for 18x18 matrices
Called by	e_c3d_u100
Calls	None
Inputs	A(18,18), B(18,18)
Outputs	C(18,18)
Modifications needed	None - create from scratch

Subroutines included:

fortran

SUBROUTINE matmul(A, B, C)
    C = A * B
    (standard matrix multiplication)

SUBROUTINE matmul_transpose(A, B, C)
    C = A^T * B
    (transpose of first argument)

FILE 8: resultsmech_u100.f (POST-PROCESSING)

Item	Description
Purpose	Compute stresses/strains for output (.dat and .frd files)
Called by	CalculiX results.f after solution
Calls	beam3d_k (optional, to reuse stiffness logic)
Inputs	v(18) - nodal displacements (global) at end of step
	xl(3,3) - nodal coordinates
	props(12) - element properties
Outputs	elout(nvar, nnode) - extrapolated nodal results for CGX
Modifications needed	None - create from scratch

What this subroutine does (step by step):

fortran

SUBROUTINE resultsmech_u100(v, xl, props, elout, ...)

STEP 1: Extract inputs (same as e_c3d_u100)
    • x1,y1,z1, x2,y2,z2, x3,y3,z3 from xl
    • E, G, A, Iy, Iz, J, ksy, ksz, yg, zg, ys, zs from props
    • L, ey, ez

STEP 2: Transform displacements from global to local
    CALL transform_beam18(x1,y1,z1, x3,y3,z3, T)
    v_local(1..18) = T(1..18, 1..18) * v(1..18)  (or T^T * v? Check convention)

STEP 3: Compute strains and stresses at Gauss points
    (Similar to beam3d_k but using B matrix and D matrix)
    
    DO g = 1, 3
        ξ = [-0.7746, 0.0, 0.7746]
        CALL shape3(ξ, N, dNdxsi)
        dNdx = dNdxsi / (L/2)
        CALL build_B(dNdx, N, B)
        
        strains(1..6) = B(1..6, 1..18) * v_local(1..18)
        stresses(1..6) = D(1..6, 1..6) * strains(1..6)
    END DO

STEP 4: Extrapolate to nodes (for contour plots)
    • Use extrapolation matrix (3x3) for quadratic element
    • elout(1,1..3) = axial stress at nodes 1,2,3
    • elout(2,1..3) = bending stress My at nodes
    • etc.

RETURN

3. SUMMARY TABLE

File	Purpose	Called by	Calls	Status
e_c3d_u100.f	Main UEL entry point	CalculiX mainsp.f	beam3d_k, transform_beam18, matmul_*	CREATE
beam3d_k.f	Local stiffness matrix (18×18)	e_c3d_u100	build_D, shape3, build_B	CREATE
build_D.f	Constitutive matrix (6×6)	beam3d_k	none	CREATE
build_B.f	Strain-displacement matrix (6×18)	beam3d_k	none	CREATE
shape3.f	Quadratic shape functions	beam3d_k	none	CREATE
transform_beam18.f	Rotation matrix (18×18)	e_c3d_u100, resultsmech_u100	none	CREATE
matmul.f	Matrix multiplication utilities	e_c3d_u100	none	CREATE
resultsmech_u100.f	Post-processing UEL	CalculiX results.f	build_D, build_B, shape3, transform_beam18	CREATE

4. FILES THAT DO NOT NEED MODIFICATION (CalculiX native)

File	Reason
mainsp.f	Already handles UEL calls - no changes needed
results.f	Already handles resultsmech_* calls - no changes needed
spooles_driver.f	Linear solver - unchanged
readinput.f	Already parses *UEL PROPERTY cards - no changes needed
readelement.f	Already recognizes U100 elements - no changes needed

5. COMPILATION COMMAND

bash

# Compile all new subroutines into a shared library
gfortran -shared -fPIC -o e_c3d_u100.so \
    e_c3d_u100.f \
    beam3d_k.f \
    build_D.f \
    build_B.f \
    shape3.f \
    transform_beam18.f \
    matmul.f \
    resultsmech_u100.f

[Disla]

The likn has probably expired. Here it is fresh again. 3 days

Shear_Centre_Section_contour_faceted.step

[bruno.zilli]

========================================
SHEAR CENTRE RESULTS

X_sc: 0.00 mm
Y_sc: 12.96 mm

These results come from a Vlasov‑type integration along the reconstructed midline chain, with automatic closed‑loop detection and drift correction.
Disclaimer:
This code is provided “as is”, for educational and research purposes only. No warranty is expressed or implied. Use at your own risk. The results shown above are also subject to this disclaimer. May contain traces of Fortran 77 and late‑night debugging. Ciao, Bruno

[Disla]

I’m obtaining the Shear center at :

-68.41 mm from the cdg

Edit: Another reference could be this one. It also shows the centre should be expected to be offset to the left side of the chanel.

[Calc_em]

Output from *BEAM SECTION GENERATE procedure in Abaqus for this profile:

section stiffness
**EA,EI11,EI12,EI22,GJ
6.98090E+08, 8.31698E+11,-3.30296E+05, 8.11200E+12, 2.23938E+09
section inertia
** rhoA, rhoI11, rhoI12, rhoI22, cmx, cmy
2.65939E-05, 3.17112E-02, 4.64524E-07, 0.30905 , 2.56025E-05,-6.58314E-07
*transverse shear stiffness
**GA11,GA22,GA12
2.68496E+08, 2.68496E+08, 0.0000
*centroid
2.56025E-05,-6.58314E-07
*shear center
-2.44622E-04, -68.409

One could also use this Python package dedicated to such calculations: sectionproperties documentation

[bruno.zilli]

Hi friends, many thanks! Now I’m gonna fix it. See you soon and thank you indeed. Bruno

[bruno.zilli]

=========================================
SHEAR CENTRE CALCULATION

File: meshes/Mesh_3.unv

Read mesh from: meshes/Mesh_3.unv
Nodes: 554
Elements: 554
Triangles: 0
Edges: 554
Nodes: 554
Elements: 554
Edges: 554
Triangles: 0

Processing 1D perimeter mesh…

Calculated properties from perimeter:
Thickness: 5.0000000000000000 mm
Area: 6704.2167726742919 mm^2
Centroid: ( -3.1921607946547903E-006 , 0.17464837638863079 ) mm
Ixx: 7994988.8810224254 mm^4
Iyy: 77657816.980200902 mm^4
Ixy: 1.7696475572355381 mm^4

Centroid (X, Y): -3.4680506477501864E-007 3.5731297815342828
Total edges: 554
Chain edges: 554
DEBUG: Spline nodes = 554
DEBUG: Total arc length = 1339.8123352271568
DEBUG: area_total = 1339.8123352271568
DEBUG: omega_mean = -13129.270156692324
DEBUG: omega end (shifted) = 20008.661368402009
DEBUG: x_sc = -0.77669954382714079
DEBUG: y_sc = -72.911157947566096

========================================
Shear Center Results:
X_sc = -0.77669989063220557 mm
Y_sc = -69.338028166031819 mm

=========================================
SHEAR CENTRE RESULTS

X_sc: -0.77669989063220557 mm
Y_sc: -69.338028166031819 mm
something like Abacus; but it’s a work in progress, don’t keep it for good (disclaimer!). Ciao, and grazie for the support. Bruno