Using a custom user element (UEL)

[bruno.zilli]

Dear CalculiX community,

Good morning! I hope this message finds you well.

I am working on a custom user element (UEL) for a 3‑node beam element with 6 degrees of freedom per node. The compilation of CalculiX (version 2.22) went fine, and the uel_ symbol is present in the executable:

nm ccx_2.22 | grep uel_
0000000000010a30 T uel_

The analysis runs without problems when I use a standard element (e.g. B31), but when I switch to my user element I encounter a segmentation fault during the stress calculation. Moreover, my UEL routine does not seem to be called at all (no debug file is written, even with a simple open/write at the very beginning of the subroutine).

My input file (simplified)

*HEADING
Test UEL beam - Simplified
*NODE
1, 0.0, 0.0, 0.0
2, 2.5, 0.0, 0.0
3, 5.0, 0.0, 0.0
*USER ELEMENT, TYPE=U1, NODES=3, MAXDOF=6, INTEGRATIONPOINTS=2
*ELEMENT, TYPE=U1, ELSET=TRAVE
1, 1, 2, 3
*MATERIAL, NAME=DUMMY
*ELASTIC
2.1e11, 0.3
*SOLID SECTION, ELSET=TRAVE, MATERIAL=DUMMY
*UEL PROPERTY, ELSET=TRAVE
1, 0.01, 8.33e-6, 8.33e-6, 0.0, 0.85, 2.1e11, 0.3, 0.0, 1.0, 0.0, 0.0, 0.0
*BOUNDARY
1, 1, 6, 0.0
*STEP
*STATIC
0.01, 1.0
*CLOAD
3, 2, -10000.0
*END STEP

Warnings and error
During the run I get these warnings:

*WARNING in calinput. Card image cannot be interpreted:
*UELPROPERTY,ELSET=TRAVE
*WARNING in calinput. Card image cannot be interpreted:
1,0.01,8.33E-6,8.33E-6,0.,0.85,2.1E11,0.3,0.,1.,0.,0.,0.

Then the analysis proceeds, the matrix structure is determined, and finally a segmentation fault occurs:

Determining the structure of the matrix:
number of equations 12
number of nonzero lower triangular matrix elements 66

Using up to 1 cpu(s) for the stress calculation.

Errore di segmentazione (core dump creato)

No debug file is created, which makes me think that the UEL routine is never entered.

What I have already tried
Using *UEL PROPERTY (with a space) instead of *UELPROPERTY – the warning persists.

Adding a simple open/write at the very beginning of uel.f – the file is never created.

Compiling with -g and running under gdb – the backtrace points to a memory access issue inside the solver, not inside my UEL.

Confirming that the uel_ symbol is indeed present in the executable.

Trying different element types: TYPE=U1, TYPE=U, TYPE=-1, TYPE=100001, with and without *USER ELEMENT.

Removing the *SOLID SECTION block – but then I get a “no material assigned” error.

My questions
Is the syntax for *UEL PROPERTY correct in my input file? The persistent warning suggests that the line is not being parsed as expected.

Why would the UEL routine not be called, even though the symbol is present and the element type seems to be recognised?

Could there be something wrong with the way the user element properties are passed? (The code in uel.f reads 13 properties; I am providing exactly 13 numbers.)

Any suggestions on how to further debug this segmentation fault?

I would be very grateful for any insights or examples of working UEL input decks for CalculiX 2.22.

Thank you very much in advance for your time and help!

Kind regards,

Bruno Zilli

Italia

File uel.f – 3‑node beam element with 6 DOF per node

subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,
 &     props,nprops,coords,mcrd,nnode,u,du,v,a,time,dtime,
 &     ttime,kstep,kinc,iel,iproc,jmatyp,
 &     iflags,iperiod,iprop,iprops,nadapt,msp,isvars,noel)

implicit none

!     ================================================================
!     INPUT / OUTPUT PARAMETERS
!     ================================================================
!     Standard CalculiX UEL interface
!
!     INPUT:
!       ndofel   : number of degrees of freedom for the element
!       nrhs     : number of right-hand sides
!       nsvars   : number of state variables
!       nprops   : number of material properties
!       mcrd     : maximum number of coordinates (3 for 3D)
!       nnode    : number of nodes for the element
!       props    : array of user-defined properties
!       coords   : nodal coordinates (mcrd, nnode)
!       u        : total displacement vector
!       du       : displacement increment
!       v        : velocity
!       a        : acceleration
!       time     : time parameters
!       dtime    : time increment
!       ttime    : total time
!       kstep    : step number
!       kinc     : increment number
!       iel      : element number
!       iproc    : processor number (for parallel runs)
!       jmatyp   : element type identifier (from \*USER ELEMENT)
!       iflags   : flags for element operations
!       iperiod  : periodicity flag
!       iprop    : property index
!       iprops   : integer properties
!       nadapt   : adaptation flag
!       msp      : maximum number of state variables
!       isvars   : integer state variables
!       noel     : element number (same as iel, kept for compatibility)
!
!     OUTPUT:
!       rhs      : residual vector (forces)
!       amatrx   : stiffness matrix (tangent matrix)
!       svars    : state variables (updated)
!       energy   : energy contributions
!     ================================================================

!     Dummy arguments
integer ndofel, nrhs, nsvars, nprops, mcrd, nnode, kstep, kinc,
&        iel, iproc, jmatyp, iperiod, iprop, nadapt, msp, noel
integer iflags(*), iprops(*), isvars(\*)

  real*8 rhs(ndofel,*)
  real*8 amatrx(ndofel,*)
  real*8 svars(*)
  real*8 energy(*)
  real*8 props(*)
  real*8 coords(mcrd,*)
  real*8 u(*), du(*), v(*), a(*), time(*), dtime, ttime

!     ================================================================
!     LOCAL VARIABLES
!     ================================================================
!     Element geometry and shape functions
integer i, j, k, ii, ipt, idx
integer nope, ndof, npoints
parameter (nope=3, ndof=6, npoints=2)   ! 3 nodes, 6 DOF/node, 2 Gauss points

  real*8 xl(3,3)          ! Nodal coordinates in local array (3 coordinates, 3 nodes)
  real*8 e1(3), e2(3), e3(3)   ! Local axes vectors
  real*8 len              ! Element length
  real*8 xi, wgt          ! Gauss point coordinate and weight
  real*8 detJ, jac        ! Jacobian and its determinant
  real*8 N(3)             ! Shape functions
  real*8 dN_dxi(3)        ! Derivative of shape functions w.r.t. natural coordinate
  real*8 dNdx(3)          ! Derivative of shape functions w.r.t. global coordinates

!     Matrices
real*8 B(6,18)          ! Strain-displacement matrix (6 strain components, 18 DOF)
real*8 D(6,6)           ! Constitutive matrix (6x6)
real*8 Ke(18,18)        ! Element stiffness matrix in global coordinates
real*8 Ke_loc(18,18)    ! Element stiffness matrix in local coordinates
real*8 tm(3,3)          ! Transformation matrix (3x3)
real*8 tmg(18,18)       ! Transformation matrix expanded to 18x18
real\*8 sg(18,18)        ! Temporary matrix for transformation

!     Section properties
real*8 area, Iyy, Izz, Iyz, shearK, E, nu, G, J_t
real*8 offset1, offset2
real\*8 dN, Ni

!     Gauss integration points
real\*8 gauss_xi(2), gauss_w(2)
data gauss_xi / -0.5773502691896257d0, 0.5773502691896257d0 /
data gauss_w / 1.0d0, 1.0d0 /

!     ================================================================
!     DEBUG OUTPUT – to verify if UEL is called
!     ================================================================
open(88, file=‘uel_called.txt’, status=‘unknown’)
write(88,*) ‘UEL called! IEL=’, iel, ’ JMATYP=', jmatyp
write(88,*) ‘NPROPS=’, nprops
write(88,\*) ‘NNODE=’, nnode, ’ NDOFEL=', ndofel
close(88)

!     ================================================================
!     READ USER-DEFINED PROPERTIES
!     ================================================================
!     Properties are passed through the \*UEL PROPERTY card
!     Expected order:
!       1: area
!       2: Iyy (moment of inertia about local y-axis)
!       3: Izz (moment of inertia about local z-axis)
!       4: Iyz (product of inertia)
!       5: shearK (shear correction factor)
!       6: E (Young’s modulus)
!       7: nu (Poisson’s ratio)
!       8-10: e2(1:3) – local y-direction vector (global coordinates)
!       11: offset1
!       12: offset2
!     ================================================================
area = props(1)
Iyy = props(2)
Izz = props(3)
Iyz = props(4)
shearK = props(5)
E = props(6)
nu = props(7)
e2(1) = props(8)
e2(2) = props(9)
e2(3) = props(10)
offset1 = props(11)
offset2 = props(12)

!     Shear modulus and torsional constant
G = E / (2.d0 \* (1.d0 + nu))
J_t = Iyy + Izz   ! Torsional constant for simple cross-sections

!     ================================================================
!     INITIALISE STIFFNESS MATRIX
!     ================================================================
do i = 1, 18
do j = 1, 18
Ke(i,j) = 0.d0
enddo
enddo

!     ================================================================
!     GAUSS INTEGRATION LOOP
!     ================================================================
do ipt = 1, npoints
xi = gauss_xi(ipt)
wgt = gauss_w(ipt)

!       Shape functions and their derivatives (quadratic 1D element)
N(1) = -0.5d0 \* xi \* (1.d0 - xi)
N(2) = (1.d0 - xi) \* (1.d0 + xi)
N(3) = 0.5d0 \* xi \* (1.d0 + xi)

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

!       Store nodal coordinates in local array
do j = 1, 3
xl(1,j) = coords(1,j)
xl(2,j) = coords(2,j)
xl(3,j) = coords(3,j)
enddo

!       Compute Jacobian
jac = 0.d0
do j = 1, 3
jac = jac + dN_dxi(1) \* xl(1,j)
jac = jac + dN_dxi(2) \* xl(2,j)
jac = jac + dN_dxi(3) \* xl(3,j)
enddo
detJ = jac / 2.d0

!       Compute element local x-axis (e1) – from node 1 to node 3
do j = 1, 3
e1(j) = (xl(j,3) - xl(j,1)) / 2.d0
enddo
len = dsqrt(e1(1)\*\*2 + e1(2)\*\*2 + e1(3)\*\*2)
do j = 1, 3
e1(j) = e1(j) / len
enddo

!       Compute local z-axis (e3) as cross product e1 × e2
e3(1) = e1(2)\*e2(3) - e1(3)\*e2(2)
e3(2) = e1(3)\*e2(1) - e1(1)\*e2(3)
e3(3) = e1(1)\*e2(2) - e1(2)\*e2(1)

!       Build transformation matrix tm (3x3) from global to local
do i = 1, 3
do j = 1, 3
tm(i,j) = 0.d0
enddo
enddo
tm(1,1) = e1(1)
tm(1,2) = e1(2)
tm(1,3) = e1(3)
tm(2,1) = e2(1)
tm(2,2) = e2(2)
tm(2,3) = e2(3)
tm(3,1) = e3(1)
tm(3,2) = e3(2)
tm(3,3) = e3(3)

!       Compute derivatives of shape functions w.r.t. global coordinates
dNdx(1) = dN_dxi(1) / detJ
dNdx(2) = dN_dxi(2) / detJ
dNdx(3) = dN_dxi(3) / detJ

!       Initialise B matrix
do i = 1, 6
do j = 1, 18
B(i,j) = 0.d0
enddo
enddo

!       Assemble B matrix (6 strain components × 18 DOF)
!       DOF ordering per node: \[u1, u2, u3, phi1, phi2, phi3\]
do ii = 1, 3
idx = (ii-1) \* 6
dN = dNdx(ii)
Ni = N(ii)

!         Axial strain
B(1, idx+1) = dN
!         Shear strains (using beam theory)
B(2, idx+2) = dN
B(2, idx+6) = -Ni
B(3, idx+3) = dN
B(3, idx+5) = Ni
!         Bending and torsional strains
B(4, idx+5) = dN
B(5, idx+6) = dN
B(6, idx+4) = dN
enddo

!       Build constitutive matrix D (6x6)
do i = 1, 6
do j = 1, 6
D(i,j) = 0.d0
enddo
enddo
D(1,1) = E \* area
D(2,2) = shearK \* G \* area
D(3,3) = shearK \* G \* area
D(4,4) = E \* Iyy
D(5,5) = E \* Izz
D(6,6) = G \* J_t

!       Compute local stiffness matrix: Ke_loc = B^T \* D \* B
do i = 1, 18
do j = 1, 18
Ke_loc(i,j) = 0.d0
do k = 1, 6
Ke_loc(i,j) = Ke_loc(i,j) + B(k,i) \* D(k,k) \* B(k,j)
enddo
enddo
enddo

!       Integrate: add contribution of current Gauss point
do i = 1, 18
do j = 1, 18
Ke(i,j) = Ke(i,j) + wgt \* Ke_loc(i,j) \* detJ
enddo
enddo

enddo   ! end of Gauss integration loop

!     ================================================================
!     TRANSFORM STIFFNESS MATRIX FROM LOCAL TO GLOBAL COORDINATES
!     ================================================================
!     Expand transformation matrix tm to 18x18
do i = 1, 18
do j = 1, 18
tmg(i,j) = 0.d0
enddo
enddo

  do i = 1, 3
    do j = 1, 3
      tmg(i,j)     = tm(i,j)
      tmg(i+3,j+3) = tm(i,j)
      tmg(i+6,j+6) = tm(i,j)
      tmg(i+9,j+9) = tm(i,j)
      tmg(i+12,j+12) = tm(i,j)
      tmg(i+15,j+15) = tm(i,j)
    enddo
  enddo

!     Transform: Ke_global = tmg^T \* Ke_local \* tmg
!     First multiply: sg = Ke_local \* tmg
do i = 1, 18
do j = 1, 18
sg(i,j) = 0.d0
do k = 1, 18
sg(i,j) = sg(i,j) + Ke(i,k) \* tmg(k,j)
enddo
enddo
enddo

!     Then multiply: Ke_global = tmg^T \* sg
do i = 1, 18
do j = 1, 18
Ke(i,j) = 0.d0
do k = 1, 18
Ke(i,j) = Ke(i,j) + tmg(k,i) \* sg(k,j)
enddo
enddo
enddo

!     ================================================================
!     ASSEMBLE OUTPUT MATRICES
!     ================================================================
!     Stiffness matrix
do i = 1, 18
do j = 1, 18
amatrx(i,j) = Ke(i,j)
enddo
enddo

!     Residual vector (set to zero for linear static analysis)
do i = 1, 18
rhs(i,1) = 0.d0
enddo

return
end

File beam_uel_test.inp – Input file

*HEADING
Test for 3-node beam UEL with 6 DOF per node
*NODE
1, 0.0, 0.0, 0.0
2, 2.5, 0.0, 0.0
3, 5.0, 0.0, 0.0
*USER ELEMENT, TYPE=U1, NODES=3, MAXDOF=6, INTEGRATIONPOINTS=2
*ELEMENT, TYPE=U1, ELSET=TRAVE
1, 1, 2, 3
*MATERIAL, NAME=DUMMY
*ELASTIC
2.1e11, 0.3
*SOLID SECTION, ELSET=TRAVE, MATERIAL=DUMMY
*UEL PROPERTY, ELSET=TRAVE
1, 0.01, 8.33e-6, 8.33e-6, 0.0, 0.85, 2.1e11, 0.3, 0.0, 1.0, 0.0, 0.0, 0.0
*BOUNDARY
1, 1, 6, 0.0
*STEP
*STATIC
0.01, 1.0
*CLOAD
3, 2, -10000.0
*END STEP

CalculiX version: 2.22

Compiler: gfortran

System: Caelinux / Ubuntu

[johanngil]

some initial comments:

• Which solver did you try? As you assume a memory problem inside the solver did you try other solvers as well?
• I would recommend to set the environment variable CCX_LOG_ALLOC=1 which writes all memory allocation and deallocation details including the function/subroutine to the console. This may show you the subroutine/function in which the segmentation fault may happen.
• As you anyway compile some code I would recommend to use the latest ccx version 2.23 as people in this forum may prefer to base further discussions on the newest release.

[Durbul]

@bruno.zilli

*WARNING in calinput. Card image cannot be interpreted:
*UELPROPERTY,ELSET=TRAVE
*WARNING in calinput. Card image cannot be interpreted:
1,0.01,8.33E-6,8.33E-6,0.,0.85,2.1E11,0.3,0.,1.,0.,0.,0.

The *UEL PROPERTY keyword is not supported by CalculiX, hence the warnings.
This means, however, that the two lines are simply ignored and it shouldn’t be the cause for a segfault.

@bruno.zilli

Why would the UEL routine not be called, even though the symbol is present and the element type seems to be recognised?

I couldn’t find any call to a subroutine uel in the code, so I believe this is not supported. I’ll need to check if there’s another logic similar to the Abaqus one in CalculiX. It seems, however, that the user element calls are hardcoded, so you would have to add them manually by modifying the source code.

@bruno.zilli

Any suggestions on how to further debug this segmentation fault?

My first try would be, as already mentioned, to test with CalculiX 2.23 (there has been a change in the U1 element from 2.22 to 2.23). Where exactly does the segfault occur? gdb should be able to provide a filename and line.

@bruno.zilli

I would be very grateful for any insights or examples of working UEL input decks for CalculiX 2.22.

As a temporary workaround, it could be worth checking the existing user subroutine e_c3d_u1.f. It has a different interface, but maybe it is possible to get your user element working in the e_c3d_u1 subroutine (by replacing the with yours / possibly some modifications required). But note that there are also postprocessing routines for u1, so these need to be adapted as well, if you want to check postprocessed quantities.

[bruno.zilli]

Hello

thanks everybody.

At the moment I’m in France. I’m going to the EDF code_aster user day in Saclay.

As soon as I’ll be back home, I’m going on with some tests.

About the second way, yes Lukas I agree. It would be a second chance.

I’ll keep you informed.

Kind regards

Bruno

[bruno.zilli]

Dear Lukas and everyone,

Thank you very much for your valuable suggestions and for taking the time to help me. I followed your advice and have been working on implementing a custom beam element in CalculiX 2.23, as you recommended.

I would like to share my progress and kindly ask for your further guidance.

What I have done so far
Compiled CalculiX 2.23 from source (successfully, with SPOOLES and ARPACK installed via apt).

Following your suggestion, I tried both approaches:

Modified e_c3d_u1.f – replacing the original Timoshenko beam formulation (2 nodes) with my quadratic beam element (3 nodes, 6 DOF per node)

Created a new e_c3d_u2.f – with the same quadratic beam code, to avoid potential conflicts with the existing U1 element

My quadratic beam element uses:

3 nodes, 6 DOF per node (UX, UY, UZ, ROTX, ROTY, ROTZ)

Euler-Bernoulli formulation

General cross-section properties: E, G, A, Iy, Iz, Jt

Modified elements.f to recognize both element types:

Added U1 and U2 to the list of solid elements (after DASHPOTA)

Added nope=3 for U1 and U2 in the nope section

Added the corresponding elseif branches for U1 and U2

Created a test input file (for U2, as you suggested):

inp
*HEADING
Test beam quadratic
*NODE
1, 0.0, 0.0, 0.0
2, 100.0, 0.0, 0.0
3, 200.0, 0.0, 0.0
*USER ELEMENT, TYPE=U2, NODES=3, PROPERTIES=6, COORDINATES=3
1,2,3,4,5,6
1
*ELEMENT, TYPE=U2, ELSET=BEAM
101, 1, 2, 3
*USER PROPERTIES
210000.0, 80769.2, 100.0, 833.33, 833.33, 200.0
*BOUNDARY
1, 1, 6, 0.0
*STEP
*STATIC
*CLOAD
3, 2, -1000.0
*END STEP
The issue
When I run the test (with either U1 or U2), I get the following error:

text
*ERROR reading *USER ELEMENT
number of integration points 32700 exceeds 255
*ERROR reading *ELEMENT
nonexistent element type:
U2
It seems that the element type is not being recognized at all. The “number of integration points” error suggests that the program is misreading the *USER ELEMENT card, perhaps because it expects a different format or because the element type is unknown.

I checked the executable to verify which subroutines are present:

bash

Check if there is a generic ‘uel’ subroutine (as in Abaqus)

nm ccx_2.23 | grep uel_

No output – no ‘uel’ subroutine exists

Check if my custom ‘e_c3d_u2’ subroutine is present

nm ccx_2.23 | grep e_c3d_u2

Output: 0000000000xxxxxx T e_c3d_u2_ (the subroutine is present and linked)

This confirms that:

My custom subroutine e_c3d_u2 is correctly compiled and linked into the executable

There is no generic uel subroutine being called for user elements

This suggests that CalculiX does not use a generic uel interface like Abaqus, and that user elements must be integrated differently.

My questions
Is my approach correct? Following your suggestion, I created a new e_c3d_u2.f and tried to register it as a new element type. Is there any additional step required to make CalculiX recognize a new user element? I have read that user elements in CalculiX are not as straightforward as in Abaqus, and I may be missing something essential.

What is the correct syntax for *USER ELEMENT in CalculiX 2.23? I am using:

TYPE=U2

NODES=3

PROPERTIES=6

COORDINATES=3

Followed by a line with the DOFs: 1,2,3,4,5,6

Followed by a line with the number of integration points: 1

Is this correct? Or is there any missing parameter?

Do I need to modify additional files (such as userelements.f, calinput.f, or calculix.c) to register a new user element type?

Could the problem be related to the number of properties? The original e_c3d_u1 expects 10 properties, while I am providing only 6. Would changing PROPERTIES=6 to PROPERTIES=10 in the *USER ELEMENT card help, even if I only use the first 6?

Is there any working example of a custom user element in CalculiX 2.23 that I could use as a reference? I have searched the documentation and the forum but could not find a complete example.

System details
CalculiX version: 2.23 (compiled on 2026-04-01)

Compiler: gcc and gfortran

Dependencies: SPOOLES, ARPACK, LAPACK, BLAS installed via apt

I would be extremely grateful for any insights or suggestions. Thank you in advance for your time and expertise.

Kind regards,
Bruno

[bruno.zilli]

Dear calculiX friends, hello everybody. I post here a structure about the software. At the end I’m writing down a roadmap to modify it. The target is to implement the U2 quadratic general section element. Every idea or contribution is very welcome. Diagram 1: Simplified Structure of CalculiX (Element Flow)┌─────────────────────────────────────────────────────────────────────────────┐
│ INPUT DECK (.inp) │
│ *NODE, *ELEMENT, *MATERIAL, *STEP, *USER ELEMENT, *USER PROPERTIES, … │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ calinput.f │
│ Reads the input file and calls appropriate subroutines for each keyword │
│ - *USER ELEMENT → reads the card but does NOT register anything │
│ (only calls getnewline) │
│ - *ELEMENT → calls elements.f │
│ - USER PROPERTIES → reads properties, stores them in prop() │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ elements.f │
│ Reads *ELEMENT, determines the element type (label) and: │
│ 1. Sets solid=.true. for solid elements (including U1, U2 if added) │
│ 2. Sets nope and nopeexp (number of nodes) │
│ 3. Stores the element label in lakon(i) │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ main (ccx_2.23.c) │
│ - Reads data, initializes data structures │
│ - Calls the appropriate subroutine based on nmethod (static, freq, etc.) │
│ For nmethod=1 → calls linstatic or nonlingeo │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ linstatic / nonlingeo │
│ - Assembles the global stiffness matrix │
│ - For each element, calls the subroutine corresponding to its type │
│ The dispatcher is in e_c3d_u.f │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ e_c3d_u.f (DISPATCHER) │
│ │
│ lakonl(2:2) = ‘1’ ──────────────────────────────────────► e_c3d_u1 │
│ lakonl(2:2) = ‘2’ ──────────────────────────────────────► e_c3d_u2 │
│ lakonl(2:4) = ‘S45’ ─────────────────────────────────────► e_c3d_us45 │
│ lakonl(2:3) = ‘S3’ ─────────────────────────────────────► e_c3d_us3 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘Diagram 2: How to Add a New User Element (U2) to CalculiX ┌─────────────────────────────────────────────────────────────────────────────┐
│ STEPS TO ADD A NEW USER ELEMENT (U2) │
└─────────────────────────────────────────────────────────────────────────────┘

Step 1: Create the element subroutine
┌─────────────────────────────────────────────────────────────────────────────┐
│ e_c3d_u2.f │
│ Contains the stiffness matrix calculation for your element │
│ - 3 nodes, 6 DOF per node (UX, UY, UZ, ROTX, ROTY, ROTZ) │
│ - Euler-Bernoulli beam formulation │
│ - General cross-section properties: E, G, A, Iy, Iz, Jt │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
Step 2: Register the element in elements.f
┌─────────────────────────────────────────────────────────────────────────────┐
│ elements.f │
│ Add U2 to two sections: │
│ 1. List of solid elements (solid=.true.) │
│ 2. nope section (nope=3, nopeexp=3) │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
Step 3: Add a branch in the dispatcher
┌─────────────────────────────────────────────────────────────────────────────┐
│ e_c3d_u.f │
│ Add: │
│ elseif(lakonl(2:2).eq.‘2’) then │
│ call e_c3d_u2(…) │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
Step 4: Recompile and link
┌─────────────────────────────────────────────────────────────────────────────┐
│ $ gfortran -c -O2 e_c3d_u2.f │
│ $ gfortran -c -O2 elements.f │
│ $ gfortran -c -O2 e_c3d_u.f │
│ $ gfortran -o ccx_2.23 *.o -lspooles … │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
Step 5: Input file syntax
┌─────────────────────────────────────────────────────────────────────────────┐
│ *USER ELEMENT, TYPE=U2, NODES=3, PROPERTIES=6, COORDINATES=3 │
│ 1,2,3,4,5,6 │
│ 1 │
│ *ELEMENT, TYPE=U2, ELSET=BEAM │
│ 101, 1, 2, 3 │
│ *USER PROPERTIES │
│ 210000.0, 80769.2, 100.0, 833.33, 833.33, 200.0 │
└─────────────────────────────────────────────────────────────────────────────┘Explanatory Notes
Why U2 is not recognized out of the box: in CalculiX, user elements are not handled through a generic uel subroutine (as in Abaqus). Instead, element types are hardcoded. The dispatcher e_c3d_u.f only recognizes:

U1 → calls e_c3d_u1

US45 → calls e_c3d_us45

US3 → calls e_c3d_us3

There is no branch for U2. Therefore, even if elements.f correctly reads the element, the dispatcher does not know which subroutine to call.

[Durbul]

Thanks for following up and sharing your findings!
I can confirm that there is no generic functionality to call user elements and that the source code needs to be modified for new implementations such as a U2 element. In my initial proposal I suggested (maybe not clear enough) to replace the U1 routines such that no further modifications are necessary. But it looks like you figured out how to add an additional one yourself.

Note that in addition to the element routine e_c3d_u1 there are postprocessing routines extrapolate_u1 and resultsmech_u1 (called in extrapolate_u and resultsmech_u respectively).
If you have a working and validated implementation, it would be great if you could share your code in a PR on GitHub

[bruno.zilli]

ciao Lukas, thank you for your kind reply; this will be for sure: a gitub with all the stuff. I’ll begin to ride the tide on tuesday.

Ciao!

Bruno

[bruno.zilli]

Hello Lukas and everybody

the idea is to bypass the calculiX standard workflow. This is a general scheme with a very little modification of the original code.

┌─────────────────────────────────────────────────────────────────────────────┐
│ INPUT DECK (.inp) │
│ *USER ELEMENT, TYPE=UZ, NODES=3, PROPERTIES=6, COORDINATES=3 │
│ 1,2,3,4,5,6 │
│ 1 │
│ *ELEMENT, TYPE=UZ, ELSET=BEAM │
│ 101, 1, 2, 3 │
│ *USER PROPERTIES │
│ 210000.0, 80769.2, 100.0, 833.33, 833.33, 200.0 │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ calinput.f (MODIFIED) │
│ Add a new branch for *USER ELEMENT, TYPE=UZ: │
│ │
│ elseif(textpart(1)(1:12).eq.‘*USERELEMENT’) then │
│ … existing code … │
│ if(textpart(i)(1:6).eq.‘TYPE=UZ’) then │
│ call userelements_z(…) ← NEW external file │
│ │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ userelements_z.f (NEW FILE) │
│ Registers UZ element in a custom array iuel_z │
│ Stores: │
│ - number of integration points (1) │
│ - degrees of freedom per node (6) │
│ - number of nodes (3) │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ elements_z.f (NEW FILE) │
│ Reads *ELEMENT, TYPE=UZ and: │
│ 1. Sets solid=.true. │
│ 2. Sets nope=3, nopeexp=3 │
│ 3. Stores label ‘UZ’ in lakon_z(i) │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ main (ccx_2.23.c) - NO CHANGE │
│ Unmodified - calls standard solvers │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ linstatic / nonlingeo - NO CHANGE │
│ Unmodified │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ e_c3d_z.f (NEW DISPATCHER) │
│ Custom dispatcher for UZ elements: │
│ │
│ if(lakonl_z(2:2).eq.‘Z’) then │
│ call e_c3d_uz(…) ← YOUR beam element │
│ endif │
└─────────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────────┐
│ e_c3d_uz.f (YOUR ELEMENT) │
│ YOUR ELEMENT SUBROUTINE │
│ - 3 nodes, 6 DOF per node (UX, UY, UZ, ROTX, ROTY, ROTZ) │
│ - Euler-Bernoulli beam formulation │
│ - General cross-section: E, G, A, Iy, Iz, Jt │
│ - Analytical stiffness matrix (no Gauss integration needed) │
│ - Returns: stiffness matrix s(18,18), mass matrix sm(18,18), │
│ load vector ff(18) │
└─────────────────────────────────────────────────────────────────────────────┘

[bruno.zilli]

Hello friends, here below is the roadmap for an abacus-like beam element. Every idea is highly appreciated. Kind Regards

COMPLETE ARCHITECTURE OF THE BEAM ELEMENT WITH GENERAL SECTION

Workflow Diagram

text

┌─────────────────────────────────────────────────────────────────────────────────────┐
│                           SECTION MESH (2D TRIANGULAR)                              │
│                         (nodes coordinates, connectivity)                           │
└─────────────────────────────────────────────────────────────────────────────────────┘
                                          │
                                          ▼
┌─────────────────────────────────────────────────────────────────────────────────────┐
│                         compute_section_properties.f                                │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  INPUT:                                                                     │    │
│  │  • nnode_sec – number of nodes in the section mesh                         │    │
│  │  • ntri_sec  – number of triangles in the section mesh                     │    │
│  │  • conn_sec(3, ntri_sec) – connectivity                                    │    │
│  │  • y_sec(nnode_sec), z_sec(nnode_sec) – nodal coordinates                  │    │
│  │  • E, G – material properties                                              │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  PROCESS:                                                                   │    │
│  │  1. Compute area A                                                         │    │
│  │  2. Compute centroid (y_c, z_c)                                            │    │
│  │  3. Compute moments of inertia (I_y, I_z, I_yz)                            │    │
│  │  4. Compute torsion constant J (approximate: I_y + I_z)                    │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  OUTPUT:                                                                     │    │
│  │  • A, y_c, z_c                                                              │    │
│  │  • I_y, I_z, I_yz                                                          │    │
│  │  • J (preliminary)                                                          │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
└─────────────────────────────────────────────────────────────────────────────────────┘
                                          │
                                          ▼
┌─────────────────────────────────────────────────────────────────────────────────────┐
│                            compute_torsion_J.f                                     │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  INPUT: (same mesh) + E, G                                                  │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  FEM SOLVER FOR TORSION (Saint-Venant):                                     │    │
│  │                                                                             │    │
│  │  assemble_poisson_2d.f                                                      │    │
│  │      └─ triangle_stiffness.f                                                │    │
│  │          → K (stiffness matrix)                                            │    │
│  │          → f (RHS = 2)                                                     │    │
│  │                                                                             │    │
│  │  apply_dirichlet_bc.f                                                       │    │
│  │      → ψ = 0 on boundary nodes                                             │    │
│  │                                                                             │    │
│  │  solve_linear_system.f                                                      │    │
│  │      → K·ψ = f → ψ (warping function)                                      │    │
│  │                                                                             │    │
│  │  compute_J.f                                                                │    │
│  │      → J = 2 ∫ ψ dA                                                        │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  OUTPUT:                                                                     │    │
│  │  • J (exact, from Saint-Venant torsion)                                    │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
└─────────────────────────────────────────────────────────────────────────────────────┘
                                          │
                                          ▼
┌─────────────────────────────────────────────────────────────────────────────────────┐
│                           compute_shear_center.f                                   │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  INPUT: (same mesh) + E, G                                                  │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  FEM SOLVER FOR SHEAR CENTER (two cases):                                   │    │
│  │                                                                             │    │
│  │  build_shear_rhs.f                                                          │    │
│  │      → Case 1 (Vy=1): f = z                                                │    │
│  │      → Case 2 (Vz=1): f = -y                                               │    │
│  │                                                                             │    │
│  │  solve_shear_problem.f (same assembly, BC, solver)                          │    │
│  │      → φ (shear stress function)                                           │    │
│  │                                                                             │    │
│  │  compute_torque.f                                                           │    │
│  │      → M_torque for each case                                              │    │
│  │                                                                             │    │
│  │  compute_shear_center.f                                                     │    │
│  │      → y_s = M_torque_case1                                                │    │
│  │      → z_s = -M_torque_case2                                               │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  OUTPUT:                                                                     │    │
│  │  • y_s, z_s (shear centre coordinates)                                      │    │
│  │  • k_y, k_z (shear correction factors)                                     │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
└─────────────────────────────────────────────────────────────────────────────────────┘
                                          │
                                          ▼
┌─────────────────────────────────────────────────────────────────────────────────────┐
│                              build_D_full.f                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  INPUT:                                                                     │    │
│  │  • E, G, A, I_y, I_z, J, k_y, k_z, y_c, z_c, y_s, z_s                      │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  COMPUTE OFFSETS:                                                            │    │
│  │  • e_y = y_s - y_c                                                          │    │
│  │  • e_z = z_s - z_c                                                          │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  D(6,6) ASSEMBLY:                                                            │    │
│  │                                                                             │    │
│  │  D(1,1) = E·A                                                               │    │
│  │  D(2,2) = k_y·G·A                                                           │    │
│  │  D(3,3) = k_z·G·A                                                           │    │
│  │  D(4,4) = G·J                                                               │    │
│  │  D(5,5) = E·I_z                                                             │    │
│  │  D(6,6) = E·I_y                                                             │    │
│  │                                                                             │    │
│  │  D(1,5) = -E·A·e_z     D(5,1) = D(1,5)                                     │    │
│  │  D(1,6) =  E·A·e_y     D(6,1) = D(1,6)                                     │    │
│  │  D(5,6) =  E·I_yz      D(6,5) = D(5,6)                                     │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  OUTPUT:                                                                     │    │
│  │  • D(6,6) – complete constitutive matrix for the section                   │    │
│  │  • elcon(12) – packed parameters for ccX                                   │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
└─────────────────────────────────────────────────────────────────────────────────────┘
                                          │
                                          ▼
┌─────────────────────────────────────────────────────────────────────────────────────┐
│                          beam_material_D.f (IN ELEMENT)                            │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  INPUT: elcon(12)                                                           │    │
│  │  OUTPUT: D(6,6)                                                             │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
└─────────────────────────────────────────────────────────────────────────────────────┘
                                          │
                                          ▼
┌─────────────────────────────────────────────────────────────────────────────────────┐
│                          beam_general_element.f                                     │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  FOR EACH GAUSS POINT:                                                      │    │
│  │                                                                             │    │
│  │  1. beam_xnor.f → local coordinate system (e1, e2, e3)                     │    │
│  │                                                                             │    │
│  │  2. beam_shape.f → N(3), B(6,18)                                           │    │
│  │                                                                             │    │
│  │  3. beam_material_D.f → D(6,6) from elcon                                  │    │
│  │                                                                             │    │
│  │  4. beam_jacobian.f → detJ, dNdx                                           │    │
│  │                                                                             │    │
│  │  5. ke = ke + Bᵀ·D·B·w·detJ                                                │    │
│  │                                                                             │    │
│  │  6. beam_transform.f → transform ke to global coordinates                  │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  OUTPUT:                                                                     │    │
│  │  • ke(18,18) – element stiffness matrix in global coordinates              │    │
│  │  • fe(18) – element force vector                                           │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
└─────────────────────────────────────────────────────────────────────────────────────┘
                                          │
                                          ▼
┌─────────────────────────────────────────────────────────────────────────────────────┐
│                              e_c3d_u.f (ccX CORE)                                   │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  DISPATCH:                                                                   │    │
│  │  if (lakon(iel)(1:7) == 'BEAMGEN') then                                     │    │
│  │      call beam_general_element(...)                                         │    │
│  │  end if                                                                     │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  GLOBAL ASSEMBLY (add_sm_st, add_sm_ei)                                     │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  SOLVER (SPOOLES / PARDISO / iterative)                                     │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
│                                          │                                         │
│                                          ▼                                         │
│  ┌─────────────────────────────────────────────────────────────────────────────┐    │
│  │  OUTPUT (.frd) – results for visualisation                                 │    │
│  └─────────────────────────────────────────────────────────────────────────────┘    │
└─────────────────────────────────────────────────────────────────────────────────────┘

---

COMPLETE LIST OF FILES

A. Section Analysis (2D FEM) – New Files

File Name	Purpose	Input	Output
compute_section_properties.f	Basic geometric properties	Mesh (nodes, triangles)	A, y_c, z_c, I_y, I_z, I_yz, J_approx
compute_torsion_J.f	Saint-Venant torsion (FEM)	Mesh + E, G	J_exact, ψ (warping)
compute_shear_center.f	Shear centre and correction factors	Mesh + E, G	y_s, z_s, k_y, k_z
assemble_poisson_2d.f	FEM assembly for Poisson equation	Mesh	K, f
triangle_stiffness.f	T3 element stiffness	Triangle coordinates	ke(3,3), fe(3), area
apply_dirichlet_bc.f	Apply ψ=0 on boundary	K, f	Modified K, f
solve_linear_system.f	Linear solver (Gauss – temporary)	K, f	ψ
compute_J.f	Integrate J from ψ	ψ, mesh	J
build_shear_rhs.f	Build RHS for shear problems	y, z, idir	f
compute_torque.f	Compute torque from φ	φ, mesh	M_torque
build_D_full.f	Build full D(6,6) from properties	A, I, J, offsets	D(6,6), elcon(12)

B. Beam Element (1D) – New Files

File Name	Purpose	Input	Output
beam_general_element.f	Main element routine	co, kon, ipkon, lakon	ke(18,18), fe(18)
beam_xnor.f	Local coordinate system	xl(3,3)	e1(3), e2(3), e3(3)
beam_shape.f	Shape functions and B-matrix	xi, detJ	N(3), B(6,18)
beam_jacobian.f	Jacobian and derivatives	xl, dN_dxi	detJ, dNdx(3)
beam_material_D.f	Build D from elcon	elcon(12)	D(6,6)
beam_transform.f	Transform ke to global	ke_loc(18,18), tm(3,3)	ke_glob(18,18)
beamintscheme.f	Integration points and weights	npoints	xi, w

C. Modified Existing Files

File	Modification
e_c3d_u.f	Add dispatch for BEAMGEN element type
Makefile	Add new files to compilation list
Makefile.inc	Add BEAMGEN to element types
beamgeneralsections.f	Parse *BEAM SECTION for general sections (optional)

---

VARIABLES SUMMARY

Geometric Properties (Section Mesh)

Variable	Description	Units
A	Cross-sectional area	m²
y_c, z_c	Centroid coordinates	m
I_y, I_z	Moments of inertia	m⁴
I_yz	Product of inertia	m⁴
J	Torsion constant	m⁴
y_s, z_s	Shear centre coordinates	m
k_y, k_z	Shear correction factors	-

Material Properties

Variable	Description	Units
E	Young’s modulus	Pa
G	Shear modulus (E/(2(1+ν)))	Pa
ν	Poisson’s ratio	-

Constitutive Matrix D(6,6)

Index	Component	Meaning
(1,1)	EA	Axial stiffness
(2,2)	k_y·G·A	Shear stiffness (Y direction)
(3,3)	k_z·G·A	Shear stiffness (Z direction)
(4,4)	G·J	Torsional stiffness
(5,5)	E·I_z	Bending stiffness (Y direction)
(6,6)	E·I_y	Bending stiffness (Z direction)
(1,5)	-E·A·e_z	Axial–bending coupling (from shear centre offset)
(1,6)	E·A·e_y	Axial–bending coupling (from shear centre offset)
(5,6)	E·I_yz	Bending–bending coupling (non-principal axes)

elcon Array (12 parameters for ccX)

Index	Variable	Description
1	EA	Axial stiffness
2	k_y·G·A	Shear stiffness (Y)
3	k_z·G·A	Shear stiffness (Z)
4	G·J	Torsional stiffness
5	E·I_y	Bending stiffness (Z)
6	E·I_z	Bending stiffness (Y)
7	c1	Coupling (axial–bending Y) = -E·A·e_z
8	c2	Coupling (axial–bending Z) = E·A·e_y
9	c3	Coupling (shear Y–bending Y)
10	c4	Coupling (shear Z–bending Z)
11	c5	Coupling (bending Y–bending Z) = E·I_yz
12	c6	Coupling (shear Y–shear Z)

---

ccX INTEGRATION POINTS

1. Input File (.inp)

text

*BEAM SECTION, ELSET=BEAMS, MATERIAL=STEEL, SECTION=GENERAL
12 parameters (elcon)
0.0, 1.0, 0.0   (orientation vector)

2. Element Dispatch (e_c3d_u.f)

fortran

if (lakonl(iel)(1:7) .eq. 'BEAMGEN') then
    call beam_general_element(...)
end if

3. Compilation (Makefile)

makefile

OCCXF = beam_general_element.o beam_xnor.o beam_shape.o beam_jacobian.o \
        beam_material_D.o beam_transform.o compute_section_properties.o \
        compute_torsion_J.o compute_shear_center.o assemble_poisson_2d.o \
        triangle_stiffness.o apply_dirichlet_bc.o solve_linear_system.o \
        compute_J.o build_shear_rhs.o compute_torque.o build_D_full.o \
        uel.o $(SCCXF:.f=.o)

---

DEVELOPMENT ROADMAP

Phase 0: Foundation (Testing outside ccX)

• Write compute_section_properties.f

• Test on rectangular section (compare with analytical formulas)

Phase 1: Torsion (Saint-Venant)

• Write assemble_poisson_2d.f, triangle_stiffness.f

• Write apply_dirichlet_bc.f, solve_linear_system.f

• Write compute_J.f

• Test on circular section (J = πR⁴/2)

Phase 2: Shear Centre

• Write build_shear_rhs.f, compute_torque.f

• Write compute_shear_center.f

• Test on rectangular section (y_s = y_c, z_s = z_c)

Phase 3: D Matrix

• Write build_D_full.f

• Test on rectangular section (diagonal D)

• Test on asymmetric section (coupled D)

Phase 4: Beam Element (1D)

• Write beam_shape.f (B-matrix, 3 nodes, 18 DOF)

• Write beam_jacobian.f, beam_xnor.f

• Write beam_material_D.f, beam_transform.f

• Write beam_general_element.f

Phase 5: ccX Integration

• Modify e_c3d_u.f for dispatch

• Update Makefile and Makefile.inc

• Test with simple cantilever beam

---

SUMMARY

Component	Files	Status
Section analysis (2D FEM)	11 new files	To implement
Beam element (1D)	7 new files	To implement
ccX integration	3 modified files	To implement
Total	18 new + 3 modified	Planned

[bruno.zilli]

Errata corrige

Torsion Constant J - Correction
the approximation J ≈ I_y + I_z is incorrect for non-circular sections. The Saint-Venant torsion problem gives:J = ∫ (ψ_y² + ψ_z²) dA (where ψ is warping function)
OR
J = 2∫ ψ dA (using Prandtl stress function)

[bruno.zilli]

Errata corrige:

The shear center is NOT computed directly from torque. Correct formulation: ! For unit shear force Vy = 1:
! M_torque = ∫ (σ_xy·z - σ_xz·y) dA
! y_s = -M_torque_case1 / Vy
! z_s = M_torque_case2 / Vz

[bruno.zilli]

first step

┌─────────────────────────────────────────────────────────────────┐
│ INPUT FILE (.inp for CalculiX) │
│ *BEAM GENERAL SECTION, ELSET=BEAMS, SECTION=FILES │
│ section_1.unv, section_2.unv, section_3.unv │
└─────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────┐
│ read_section_mesh_unv.f (NEW) │
│ • Reads .unv file (I-DEAS Universal File format) │
│ • Extracts nodes (coordinates) and triangles (connectivity) │
└─────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────┐
│ compute_all_sections_properties.f (NEW) │
│ • Manages array of sections (nsections) │
│ • Calls compute_section_properties for each section │
│ • Stores results in multidimensional arrays │
└─────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────┐
│ section_properties_cache (MODULE/GLOBAL COMMON BLOCK) │
│ • properties(nsections) = {A, y_c, z_c, I_y, I_z, I_yz, …} │
│ • Fast access during FEM analysis │
└─────────────────────────────────────────────────────────────────┘

[bruno.zilli]

BEAM SECTION TOOL - Compute cross-section properties from UNV (Salome) meshes

hello friends; this is a little first step: GitHub - brunozilli/beam_section_test · GitHub

every your idea or insight is very welcome. Kind Regards. Bruno

[xyont]

beam section properties based on mesh is great specifically on warping (Vlasov torsion) thank you, integration in solver with predefined common section may interesting. Also, it can be good without dependency in Salome UNV mesh, CGX capable to mesh triangle internally and saved in MSH format of Abaqus or Tetgen mesh format available in distribution.

[bruno.zilli]

Hello Dr. Dhont,

thank you very much for your kind reply and for your valuable suggestions.

Here is another little step forward, implementing the FEM calculation of the torsional constant J (Prandtl stress function) using the triangular mesh already available:

https://github.com/brunozilli/beam_section_J

The code is still experimental, but it already solves ∇²φ = -2 on the cross‑section and computes J = 2∫φ dA.

Any feedback or advice on how to improve the accuracy (especially for the HEB200 test case) would be most welcome.

Kind Regards,
Bruno Zilli

[bruno.zilli]

Dear CalculiX Community,

I hope you are well. I would like to share a small tool I have been developing, in case it might be useful to others.

beam_shear_center computes the shear centre of beam cross-sections using FEM (Poisson problems with Neumann boundary conditions). It is written in Fortran 77, reads UNV meshes, and uses LAPACK.

The code is available here:

GitHub - brunozilli/beam_shear_center

It is still under active development, but has been validated on symmetric sections (HEB200). Feedback, suggestions, or contributions would be most welcome.

Thank you very much for your time.

With kind regards,

Bruno Zilli

[bruno.zilli]

Hello calculiX friends; I did a simple test on a L 100x100x10 section. ========================================
Test: L100x100x10 45° rotated and CG centered

Read mesh from: meshes/Mesh_test_L.unv

Nodes: 650
Triangles: 1103
Mesh read: nodes = 650 elements = 1103
Section Properties:
Area: 1915.0990398593233
Centroid (y_c, z_c): 8.5161379311373582E-003 8.5161379311339755E-003
Iy: 1764444.2544620314 Iz: 1764444.2544620261 Iyz: -1035288.8769560433

DEBUG: compute_shear_center (FINAL VERSION)

DEBUG: nn= 650 ne= 1103
DEBUG: Iy= 1764444.2544620314 Iz= 1764444.2544620261 Iyz= -1035288.8769560433
DEBUG: y_bar = 8.5161379311367875E-003 z_bar = 8.5161379311339876E-003
DEBUG: total_area = 1915.0990398593233
DEBUG: assembly completed
DEBUG: RHS normalized
DEBUG: singularity fixed (node 1 pinned)
DEBUG: DGELS completed successfully
DEBUG: M_y = 68562244.724935561
DEBUG: M_z = -68563444.054936543
DEBUG: detI = 2041440468355.1592
DEBUG: y_s = 24.488241786613262
DEBUG: z_s = -24.489886607508929
Shear Center (y_s, z_s): 24.488241786613262 -24.489886607508929

Expected values (from theory):
y_s ≈ 6.8 mm (from centroid)
z_s ≈ 6.8 mm (from centroid)
(absolute from corner: ~35.7 mm)
WARNING: y_s != z_s, diff = 48.978128394122194
PASS: Shear center is NOT at centroid (expected)
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL link: GitHub - brunozilli/beam_shear_center · GitHub

[bruno.zilli]

Actually, the mesh is oriented straight with the corner in the origin, and the edges parallel to the axes.

[bruno.zilli]

Hello calculiX friends, happy Sunday! here we go!

GitHub - brunozilli/beam_shear_center_test_L: Shear center calculation with L section 100x100x10...

It works for HEB200 and a L 100x100x10 oriented with y, z positives.

Please, you’re very welcome to stretch-test-eventually eatinonego everything. Thank you in advance.

Kind regards (Ciao!)

Bruno Zilli