Dr. Huang's CPFEM UMAT in CalculiX

Atsushi · April 20, 2024, 9:31pm

Hello, everyone.

Does anyone know how to use UMAT with Dr. Huang’s crystal plasticity finite element method(CPFEM)?
The program and input files can be downloaded from the following website.

https://www.researchgate.net/post/Can_anyone_help_me_to_implement_properly_Huangs_UMAT_for_single_crystal_plasticity_on_Abaqus

“umatcrystal_mod.f” is the UMAT program, and “umatcryspl_mod.inp” is the input file.

The input file contains several keywords that CalculiX does not recognize. Therefore, they must be modified in order to be calculated. The umatcryspl_mod2.inp is a fixed version of them.

I recompiled the program and ran the input file. The results of the run are as follows.

atsushi@atsushi:~/tmp$ ~/CalculiX_2.21_umatHuang/src/ccx_2.21_MT -i umatcryspl_mod

************************************************************

CalculiX Version 2.21, Copyright(C) 1998-2023 Guido Dhondt
CalculiX comes with ABSOLUTELY NO WARRANTY. This is free
software, and you are welcome to redistribute it under
certain conditions, see gpl.htm

************************************************************

You are using an executable made on Sat Jul 29 10:52:01 CEST 2023

  The numbers below are estimated upper bounds

  number of:

   nodes:           20
   elements:            1
   one-dimensional elements:            0
   two-dimensional elements:            0
   integration points per element:            8
   degrees of freedom per node:            3
   layers per element:            1

   distributed facial loads:            0
   distributed volumetric loads:            0
   concentrated loads:            0
   single point constraints:           24
   multiple point constraints:            1
   terms in all multiple point constraints:            1
   tie constraints:            0
   dependent nodes tied by cyclic constraints:            0
   dependent nodes in pre-tension constraints:            0

   sets:            5
   terms in all sets:           25

   materials:            1
   constants per material and temperature:          160
   temperature points per material:            2
   plastic data points per material:            0

   orientations:            0
   amplitudes:            3
   data points in all amplitudes:            3
   print requests:            0
   transformations:            0
   property cards:            0

 *ERROR reading *USER MATERIAL: anisotropic definition
   is not complete.
 *ERROR reading *USER MATERIAL. Card image:
        *DEPVAR

Segmentation fault
atsushi@atsushi:~/tmp$

Unfortunately, an error occurs.
It seems to be an error in usermaterials.f.
Also, the same error occurs in the original “umatcrystal_mod.f”.
I don’t know why.

Does anyone know?
Please tell me, any small hints are welcome.

Thanks,
Atsushi

“umatcryspl_mod2.inp”

**
** Heading +++++++++++++++++++++++++++++++++++++++++++++++++
**
*Heading
CPFEM test
**Hash: rrJYTz6R, Date: 04/11/2024, Unit system: MM_TON_S_C
**
** Nodes +++++++++++++++++++++++++++++++++++++++++++++++++++
**
*Node
1, 0.00000000E+000, 0.00000000E+000, 0.00000000E+000
2, 0.00000000E+000, 1.00000000E+001, 0.00000000E+000
3, 0.00000000E+000, 1.00000000E+001, 1.00000000E+001
4, 0.00000000E+000, 0.00000000E+000, 1.00000000E+001
5, 1.00000000E+002, 0.00000000E+000, 0.00000000E+000
6, 1.00000000E+002, 1.00000000E+001, 0.00000000E+000
7, 1.00000000E+002, 1.00000000E+001, 1.00000000E+001
8, 1.00000000E+002, 0.00000000E+000, 1.00000000E+001
9, 0.00000000E+000, 5.00000000E+000, 0.00000000E+000
10, 0.00000000E+000, 1.00000000E+001, 5.00000000E+000
11, 0.00000000E+000, 5.00000000E+000, 1.00000000E+001
12, 0.00000000E+000, 0.00000000E+000, 5.00000000E+000
13, 1.00000000E+002, 5.00000000E+000, 0.00000000E+000
14, 1.00000000E+002, 1.00000000E+001, 5.00000000E+000
15, 1.00000000E+002, 5.00000000E+000, 1.00000000E+001
16, 1.00000000E+002, 0.00000000E+000, 5.00000000E+000
17, 5.00000000E+001, 0.00000000E+000, 0.00000000E+000
18, 5.00000000E+001, 1.00000000E+001, 0.00000000E+000
19, 5.00000000E+001, 1.00000000E+001, 1.00000000E+001
20, 5.00000000E+001, 0.00000000E+000, 1.00000000E+001
**
** Elements ++++++++++++++++++++++++++++++++++++++++++++++++
**
*Element, Type=C3D20R
1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20
**
** Node sets +++++++++++++++++++++++++++++++++++++++++++++++
**
*Nset, Nset=Fix123
1
*Nset, Nset=Fix13
2
*Nset, Nset=Fix1
3, 4, 9, 10, 11, 12
*Nset, Nset=Move
5, 6, 7, 8, 13, 14, 15, 16
**
** Element sets ++++++++++++++++++++++++++++++++++++++++++++
**
*Elset, Elset=Solid
1
**
** Surfaces ++++++++++++++++++++++++++++++++++++++++++++++++
**
**
** Physical constants ++++++++++++++++++++++++++++++++++++++
**
**
** Materials +++++++++++++++++++++++++++++++++++++++++++++++
**
*Material, Name=ABAQUSNL1
*USER MATERIAL,CONSTANTS=160
**
** All the constants below must be real numbers!
**
 168400.0, 121400.0,  75400.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    c11  ,   c12  ,   c44  , (elastic constants of copper crystal)
**    MPa  ,   MPa  ,   MPa  , 
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
** constants only used for an elastic orthotropic or anisotropic material
**    MPa  ,
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
** constants only used for an elastic anisotropic material
**    MPa  ,
**
** The elastic constants above are relative to crystal axes, where
**   1 -- [100],  2 -- [010],  3 -- [001] .  These elastic constants 
**   are arranged in the following order:  
**   eight constants each line (data card)
**
** (1) isotropic: 
**     E    ,  Nu      (Young's modulus and Poisson's ratio)
**     0.
**     0.
**
** (2) cubic:
**     c11  ,  c12  ,  c44
**     0.
**     0.
**
** (3) orthotropic:
**     D1111,  D1122,  D2222,  D1133,  D2233,  D3333,  D1212,  D1313,  
**     D2233
**     0.
**
** (4) anisotropic:
**     D1111,  D1122,  D2222,  D1133,  D2233,  D3333,  D1112,  D2212,
**     D3312,  D1212,  D1113,  D2213,  D3313,  D1213,  D1313,  D1123,
**     D2223,  D3323,  D1223,  D1323,  D2323
**
**
      1.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
** number of sets of slip systems
**    --   ,
**
      1.0,      1.0,      1.0,      1.0,      1.0,      0.0,      0.0,      0.0,
**    normal to slip plane   ,      slip direction      , of the 1st set
**    --   ,   --   ,   --   ,   --   ,   --   ,   --   ,
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    normal to slip plane   ,      slip direction      , of the 2nd set
**    --   ,   --   ,   --   ,   --   ,   --   ,   --   ,
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    normal to slip plane   ,      slip direction      , of the 3rd set
**    --   ,
**
**
     -1.0,      0.0,      1.0,      0.0,      0.0,      1.0,      0.0,      0.0,
** direction in local system ,      global system       , of the 1st vector
**    ---  ,   --   ,   --   ,   --   ,   --   ,   --   ,
** (the first vector to determine crystal orientation in global system)
**
      0.0,      1.0,      0.0,      0.0,      1.0,      0.0,      0.0,      0.0,
** direction in local system ,      global system       , of the 2nd vector
**    --   ,   --   ,   --   ,   --   ,   --   ,   --   ,
** (the second vector to determine crystal orientation in global system)
**
** constraint:  The angle between two non-parallel vectors in the local
**              and global systems should be the same.  The relative 
**              difference must be less than 0.1%. 
**
**
     10.0,    0.001,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**     n   ,  adot  , of 1st set of slip systems
**    ---  ,  1/sec ,
** (power hardening exponent and hardening coefficient)
**  gammadot = adot * ( tau / g ) ** n
**
** Users who want to use their own constitutive relation may change the
**   function subprograms F and DFDX called by the subroutine 
**   STRAINRATE and provide the necessary data (no more than 8) in the 
**   above line (data card).
**
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    n    ,  adot  , of 2nd set of slip systems
**    ---  ,  1/sec ,
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    n    ,  adot  , of 3rd set of slip systems
**    --   ,  1/sec ,
**
**
    541.5,    109.5,     60.8,      0.0,      0.0,      0.0,      0.0,      0.0,
**    h0   ,  taus  ,  tau0  , of 1st set of slip systems
**    MPa  ,   MPa  ,   MPa  ,
** (initial hardening modulus, saturation stress and initial critical 
**  resolved shear stress)
**  H = H0 * { sech [ H0 * gamma / (taus - tau0 ) ] } ** 2
**
** Users who want to use their own self-hardening law may change the 
**   function subprogram HSELF called by the subroutine LATENTHARDEN 
**   and provide the necessary data (no more than 8) in the above line 
**   (data card).
**
**
      1.0,      1.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    q    ,   q1   , Latent hardening of 1st set of slip systems
**    --   ,   --   ,
** (ratios of latent to self-hardening in the same and different sets 
**   of slip systems)
**
** Users who want to use their own latent-hardening may change the 
**   function subprogram HLATNT called by the subroutine LATENTHARDEN 
**   and provide the additional data (beyond the self-hardening data, 
**   no more than 8) in the above line (data card).
**
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    h0   ,  taus  ,  tau0  , of 2nd set of slip systems
**    MPa  ,   MPa  ,   Mpa  ,
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    q    ,   q1   , of 2nd set of slip systems
**    --   ,   --   ,
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    h0   ,  taus  ,  tau0  , of 3rd set of slip systems
**    MPa  ,   MPa  ,   MPa  ,
**
      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**    q    ,   q1   , of 3rd set of slip systems
**    --   ,   --   ,
**
**
      0.5,      1.0,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,
**   THETA , NLGEOM ,
**    --   ,   --   ,
**
** THETA:  implicit integration parameter, between 0 and 1
**
** NLGEOM:  parameter determining whether finite deformation of single 
**   crystal is considered
**
**   NLGEOM=0. --- small deformation
**   otherwise --- finite rotation and finite strain,  Users must 
**                 declare "NLGEOM" in the input file, at the *STEP 
**                 card
**
**
      1.0,     10.0,   1.0E-5,      0.0,      0.0,      0.0,      0.0,      0.0
**  ITRATN , ITRMAX , GAMERR ,
**    --   ,   --   ,   --   ,
** ITRATN:  parameter determining whether iteration method is used to 
**   solve increments of stresses and state variables in terms of 
**   strain increments
**
**   ITRATN=0. --- no iteration
**   otherwise --- iteration
**
** ITRMAX:  maximum number of iterations
**
** GAMERR:  absolute error of shear strains in slip systems
**
**
*DEPVAR
125
**
** Sections ++++++++++++++++++++++++++++++++++++++++++++++++
**
*Solid section, Elset=Solid, Material=ABAQUSNL1
**
** Pre-tension sections ++++++++++++++++++++++++++++++++++++
**
**
** Constraints +++++++++++++++++++++++++++++++++++++++++++++
**
**
** Surface interactions ++++++++++++++++++++++++++++++++++++
**
**
** Contact pairs +++++++++++++++++++++++++++++++++++++++++++
**
**
** Amplitudes ++++++++++++++++++++++++++++++++++++++++++++++
**
**
** Initial conditions ++++++++++++++++++++++++++++++++++++++
**
**
** Steps +++++++++++++++++++++++++++++++++++++++++++++++++++
**
**
** Step-1 ++++++++++++++++++++++++++++++++++++++++++++++++++
**
*Step, Nlgeom, Inc=500
*Static, Solver=Spooles
0.0005, 1, 1E-07, 0.2
**
** Output frequency ++++++++++++++++++++++++++++++++++++++++
**
*Output, Frequency=1
**
** Boundary conditions +++++++++++++++++++++++++++++++++++++
**
*Boundary, op=New
** Name: Displacement_rotation-1
*Boundary
Fix123, 1, 1, 0
** Name: Displacement_rotation-2
*Boundary
Fix123, 2, 2, 0
** Name: Displacement_rotation-3
*Boundary
Fix123, 3, 3, 0
** Name: Displacement_rotation-4
*Boundary
Fix13, 1, 1, 0
** Name: Displacement_rotation-5
*Boundary
Fix13, 3, 3, 0
** Name: Displacement_rotation-6
*Boundary
Fix1, 1, 1, 0
** Name: Move
*Boundary
Move, 1, 1, 0.1
Move, 2, 2, 0
Move, 3, 3, 0
**
** Loads +++++++++++++++++++++++++++++++++++++++++++++++++++
**
*Cload, op=New
*Dload, op=New
**
** Defined fields ++++++++++++++++++++++++++++++++++++++++++
**
**
** History outputs +++++++++++++++++++++++++++++++++++++++++
**
**
** Field outputs +++++++++++++++++++++++++++++++++++++++++++
**
*Node file
U
*El file
S, E, NOE
*Output, Frequency=1
**
** End step ++++++++++++++++++++++++++++++++++++++++++++++++
**
*End step

jbr · April 20, 2024, 10:00pm

Did you modify all the necessary files? It looks like you need to define the depvars and potentially other initial variables.

I looked at one of the files, and this will not simply compile as given, for example:

include 'aba_param.inc' → Needs to be removed.
EXTERNAL F → It needs to be removed or fixed with the corresponding global variable. This is not recommended unless you truly know what you are doing.
Other internal subroutines → Need to check one by one and make sure they do not exist in ccx right now.

Atsushi · April 21, 2024, 12:03am

Thanks for telling me know.

Sorry, I forgot to attach the file.
I used “aba_param.inc” to compile umat.f.

I have not changed about “EXTERNAL F”.

How to compile
Replace the file name “umatcrystal_mod.f” with “umat.f”,
Recompile with make command.
Now you can compile.

Are you saying that this work is not enough to implement UMAT?

“aba_param.inc”

      implicit real*8(a-h,o-z)
      parameter (nprecd=2)

k3rrig4n · January 23, 2026, 7:45am

Hello. Have you been able to solve this problem?

johanngil · January 23, 2026, 11:22am

Please replace the ccx (version 2.23) fortran files allocation.f, usermaterials.f, headings.f and umat_abaqusnl.f by the attached modified file versions (Unique Download Link | WeTransfer) and recompile ccx. Then the 2 attached inp (be aware of the user material name being used!) files with the Huang/Kysar single crystal material runs with ccx as it provides some console output during the execution. The results can be postprocessed as well.

There are a couple of aspects why we see here in the forum recurring problems with user defined materials especially when they have >8 material parameters:

The original inp file for Abaqus has the temperature only once after the last parameter. Calculix needs it as last parameter in every line (in case temperature dependency is needed). So the syntax in the inp file seems to be different in Abaqus vs. ccx to get it run with the same umat Fortran subroutine.
The ccx manual (version 2.23) states on page 630: „if CONSTANTS is a multiple of 8“. I don’t think that this is correct, „if CONSTANTS is greater than 8“ may be a more precise statement.
I have not found any ccx example standard test file with a user material having >8 constants being defined for multiple temperatures. I have tried different syntax variants for multi-temperature setups for CONSTANTS > 8, but I did not manage to get it run. So the attached Fortran files are expected to only work for temperature independent user material definitions.

k3rrig4n · January 26, 2026, 5:49am

Thanks for the answer. I tried your suggestion for solving this problem, but i faced other issue. I replaced fortran files, added aba_param.inc file, i renamed Dr. Huang’s umat into umat_HUANG_KYSAR_SINGLE_CRYSTAL.f, i also modified umat_main.f and makefile.inc files but i got this error:

Is there any advise on how to solve that error?

johanngil · January 26, 2026, 4:12pm

Why did you modify umat_main? umat_HUANG_KYSAR_SINGLE_CRYSTAL.f will be called from umat_abaqusnl(), not from umat_main.f. Keep in mind that the umat was originally developed for Abaqus why the sequence of calls should be umat_main() –> umat_abaqusnl() –> umat_HUANG_KYSAR_SINGLE_CRYSTAL().

The error in your screenshot means that the linker could not find subroutine umat_HUANG_KYSAR_SINGLE_CRYSTAL(). Please make sure that the subroutine in umat_HUANG_KYSAR_SINGLE_CRYSTAL.f has the name umat_HUANG_KYSAR_SINGLE_CRYSTAL(). You need to change its original name.

Please further check if umat_HUANG_KYSAR_SINGLE_CRYSTAL.o was successfully created by the compiler to enable the linker to find it.

k3rrig4n · January 27, 2026, 4:24am

Thank you for the reply. Turns out name of subroutine in umat is wrong. I managed to recompile it, but your example’s arent working. It say’s: ***ERROR - first vector is zero. I tried to change 1st and 2nd vectors, but it didnt work. Can you please help on that one?

UPD: Looks like Calculix is not reading PROPS.

UPD2: The input parameters are not passed from umat_abaqusnl to umat_HUANG_KYSAR_SINGLE_CRYSTAL

johanngil · January 27, 2026, 4:01pm

I successfully reran both examples on my side: The original inp model example umat_huang_kysar_cryspl04.inp for the UMAT shows during the run console output like:

After completion of the run I see the deformed brick in the postprocessor:

k3rrig4n · January 28, 2026, 8:23am

I am sorry to bother you that much, but can you please send your version of umat? Is it different from the version above?

johanngil · January 28, 2026, 3:46pm

Here is my umat version of the Huang/Kysar crystal plasticity model: Unique Download Link | WeTransfer

k3rrig4n · January 29, 2026, 9:26am

Thank you very much. It is working now.

Topic		Replies	Views
How to create a user material?	5	1364	July 21, 2022
Umatht.f for transient thermal analysis	9	1460	July 16, 2024
CalculiX not reading PROPS	12	1078	January 21, 2022
Segmentation fault while using UMAT in CalculiX	6	628	June 13, 2021
Using abaqus umat from shared library	8	484	January 14, 2025

Dr. Huang's CPFEM UMAT in CalculiX

Related topics