Dr. Huang's CPFEM UMAT in CalculiX

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

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.

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)