2.4. Customizing a Coupled-Field Element

User-programmable features (UPFs) are available with coupled-field elements PLANE222, PLANE223, SOLID225, SOLID226, and SOLID227. Coupled-field UPFs enable you to customize element characteristics, material models, loads, constitutive equations, and output items. The UPFs give you access to all the applicable element results at any location within the element, making it easier to program nonlinearities. The coupled-field UPFs are integrated into the element standard workflow, so you can customize the element property of interest without writing a user element or a user material model from scratch. For your convenience, there are several utility functions available to help you access element results and other element information.

The following topics are available:

2.4.1. Activating Programmable Features for Supported Analysis Types

User-programmable capabilities are available with coupled-field elements PLANE222, PLANE223, SOLID225, SOLID226, and SOLID227 for the following coupled-diffusion analyses:

Structural-diffusion (KEYOPT(1) = 100001)
Thermal-diffusion (KEYOPT(1) = 100010)
Electric-diffusion (KEYOPT(1) = 100100)
Structural-thermal-diffusion (KEYOPT(1) = 100011)
Structural-electric-diffusion (KEYOPT(1) = 100101)
Structural-thermal-electric-diffusion (KEYOPT(1) = 100111)

To activate user-programmable features (UPFs), set KEYOPT(12) = 1 on the coupled-field elements used in any of the above analyses. Activating UPFs enables you to define or modify the element material properties, loads, constitutive equations, as well as element results. Multiple functions and subroutines are available to provide you convenient access the element and solution information.

2.4.2. Files on the Distribution Media

The UserProcedures.F file contains the source code for the routines that provide an interface to the coupled-field element code and allow you to customize its behavior. This file can be found on the distribution medium along with other user routines (see Studying the Mechanical APDL User Routines).

The coupled-field UPF code is implemented as a Fortran submodule programming unit that contains two procedures (Figure 2.1: UserProcedures.F routine ):

user_char: To modify element characteristics, including number of work variables, saved variables, SMISC and NMISC items, element matrix symmetry, and element nonlinearity (see Modifying the Coupled-field Element Characteristics).
user_call: To modify element calculations, such as material properties, loads, element matrices, and load vectors using the utility functions described in Customizing Element Calculations.

Both procedures must be present for UserProcedures.F to be successfully compiled and linked.

Figure 2.1: UserProcedures.F routine

You can alter the submodule name (UserProcedures in Figure 2.1: UserProcedures.F routine ) but not the module name (UserInterface).

When compiling UserProcedures.F, the UPF tools will automatically access the ancestor module file userinterface.mod that must be present at the following path:

Linux

/ansys_inc/v251/ansys/customize/include/userinterface.mod

Windows

Program Files\ANSYS Inc\V251\ansys\customize\include\userinterface.mod

2.4.3. Coding Coupled-field UPFs

The element code makes additional calls to the two subroutines user_char and user_call when user-programmable features (UPFs) for a coupled-field analysis are enabled (KEYOPT(12) = 1). These calls enable you to customize element characteristics and element calculations by modifying the user_char and user_call procedures, respectively.

Figure 2.2: Coupled-Field Element in User-programmable Mode

The next sections describe how user_char and user_call procedures interact with the standard coupled-field element behavior. The following abbreviations are used in the description of input and output arguments:

int – integer

dp – double precision

log – logical

sc – scalar

ar – array

in – input

out – output

2.4.3.1. Modifying the Coupled-field Element Characteristics

Element characteristics are programmed as an array of integer numbers (approximately 200) that define the element and control its behavior. The user_char procedure allows you to reset some of these characteristics. All other characteristics default automatically. user_char is called once from the element routine that sets these default characteristics.

The information about element characteristics is passed to user_char in the form of three data structures: num, non, and elem. The interface information (Figure 2.3: user_char Procedure Interface) resides in the UserInterface module and is not replicated in the procedure body.

Figure 2.3: user_char Procedure Interface

A brief description of user_char arguments is given in the following table.

Table 2.1: user_char Procedure Arguments

Name	Type	Content
num	ielc_num	sizes of element data arrays
non	ielc_non	element formulations
elem	ielc_elem	element definitions

Components of the num structure can be used to customize some of the element characteristics as follows:

num%item
where:
item - name of the characteristics.

Table 2.2: num Structure Components

*`Item`*	Description	Type
rcon	number of real constants	int sc
smisc	number of summable miscellaneous items
nmisc	number of non-summable miscellaneous items
usvr	number of user saved variables
uwrk	number of user variables

For example,

num%rcon = 5

sets the number of real constants to 5, and

num%smisc = num%smisc + 2

adds 2 to the default number of summable miscellaneous. Note that you can only increment the number of real constants (num%rcon), the number of summable (num%smisc), or non-summable (num%nmisc) variables. In other words, the custom number of real constants or miscellaneous variables cannot not be less than the number of standard real constants miscellaneous variables listed in element descriptions (Element Reference).

Components of the non structure can be used to set the characteristics that control the element formulation:

non%item
where:
item - name of the element property.

Table 2.3: non Structure Components

Item Description Type Value

sym

element matrix

int sc

0 – symmetric

1 – unsymmetric

2 – symmetric, but may become unsymmetric due to certain material models

lin

element linearity

0 – linear element

1 – nonlinear element

For example,

non%sym = 1

sets the element matrix to nonsymmetric.

When modifying the above element characteristics, you can access (use in read-only mode) the following components of the elem structure:

elem%item
where:
item - name of the characteristics.

Table 2.4: elem Structure Components

*`Item`*	Decription	Type
type	element type	int sc
name	element name
keyop(kyop)	element keyoption (kyop = 1 from 18)

The following example illustrates customization of the element characteristics using the user_call procedure.

Example 2.1: Setting Element Characteristics Using user_char

2.4.3.2. Customizing Element Calculations

Typical element calculations have three distinct phases: data preparation, integration point calculations, and the output of element results (Figure 2.4: Element Calculation Workflow).

Data preparation — Material properties and element loads are first evaluated during the data preparation stage.

Integration point calculations — Material properties, loads, and element results are then evaluated at each integration point using available solution information, for example, element results from the previous iteration. Element matrices and load vectors are formed in the same integration loop using numerical integration over the element volume. Note that data preparation and numerical integration are performed for each physics and coupled-field effect.

Output of element results — During the output stage, all calculated element quantities are processed, and element results are written to the result file.

The user_call procedure is called multiple times during the three stages of element calculations. You can use user_call to customize the evaluation of material properties and loads at both the data preparation and integration point stage, modify the constitutive relations at the integration point stage, and modify the element miscellaneous items in the output phase of element calculations.

Figure 2.4: Element Calculation Workflow

The information about the coupled-field element that can be passed to and used in user_call is in the form of four structures: num, elem, solu, and scope. The interface information (Figure 2.5: user_call Procedure Interface) resides in the UserInterface module and is not replicated in the procedure body.

Figure 2.5: user_call Procedure Interface

A brief description of user_call arguments is given in the following table.

Table 2.5: user_call Procedure Arguments

Name	Type	Content
num	`user_num`	Numbers that characterize element entities.
elem	`user_elem`	Components and procedures to access element characteristics, nodes, coordinates, shape functions, and results. Procedures to perform operations on element results and print these results. Procedures to define or modify material properties, loads, constitutive equations, and output quantities
solu	`user_solu`	Solution settings and data to control element calculations.
scope	`user_scope`	Parameters to define the scope of element calculations.

The components and type-bound procedures of these structures are described in detail below.

Element numbers (num)

The num structure consists of frequently used element numerical property data. All the components of this structure are integer scalar parameters. You can use these parameters to define the sizes of variables in the user_call procedure or control do-loop and if-constructs:

num%item
where:
item – name of element property

Table 2.6: Element Numerical Property Data

*`Item`*	Description	Type
dim	element dimensionality (2 or 3)	int sc
nodes	number of element nodes
cnodes	number of element corner nodes
edges	number of element edges
faces	number of faces
ptsf	number of surface load values per face
fluxes	number of element fluxes (= faces*ptsf)
compn	number of stress or strain vector components
comp	number of field or flux vector components (= dim)
intp	number of integration points
dofs	number of degrees of freedom (DOFs) per node
rcon	number of real constants^[a]
uwrk	number of user work variables^[a]
usvr	number of user saved variables ^[a]
smisc	number of SMISC items ^[a]
nmisc	number of NMISC items ^[a]

^[a]Default or set in user_char

Example 2.2: Using num Components in user_call

Element Data (elem)

The elem structure consists of components and type-bound procedures that let you access the following data and use it to customize element behavior:

Element characteristics and attributes, data arrays.
Nodal coordinates, element shape functions, and shape function derivatives.
Nodal solution results, such as displacements, temperature, electric potential, and concentration.
Element results, such as stress and strains, thermal gradient and flux, electric field and current density, concentration gradient and diffusion flux, as well as miscellaneous output items.

The elem type-bound procedures also allow you to:

Perform some operations on the element nodal data, such as interpolation for missing midside nodes and gradient calculation.
Print nodal and element results in a formatted fashion.

The procedures that allow you to customize material properties, loads, constitutive relations, and element output (the so-called set_ and add_ subroutines) are also implemented as type-bound procedures of the elem structure.

Element data can be accessed as elem components

elem%item

or type-bound procedures

elem%item()
where:
item - name of the element entity

Table 2.7: Element Characteristics and Attributes

Item Description Output Type

id element number int sc

ip() current integration point number

type element type

name element name

mat material number

real real constant set number

csys element coordinate system number

keyop() element keyoptions int array (18)

nodes() element nodes int array (num%nodes)

shape

element form shape:

0 - quadrilateral or brick

1 - triangle or tetrahedron

int sc

ashape

actual element shape:

0 - quadrilateral or brick

1 - triangle or wedge

3 - pyramid

midn

element midside nodes key:

.true. if element has midsize nodes

log sc

vol() element volume dp sc

For example,

elem%id

returns the element number, and

elem%ip()

returns the integration point number for the current call to user_call. The use of some typical element properties is illustrated in the following example.

Example 2.3: Example Using elem Components in user_call

The elem structure gives you access to several arrays of element data, listed in the following table. You can access the whole array as

elem%item()

or a single element of this array as

elem%item(n)
where:
item – name of the data array
n – array entry number

Table 2.8: Element Data Arrays

*`Item`*	Description	Output Type
rcon	element real constants	dp ar (num%rcon)
smisc	summable miscellaneous items	dp ar (num%smisc)
nmisc	non-summable miscellaneous items	dp ar (num%nmisc)
usvr^[a]	user saved variables	dp ar (num%usvr)
uwrk^[a]	user work variables	dp ar (num%uwrk)

^[a]Can be overwritten.

The arrays are read-only, except for the usvr and uwrk arrays that you can both read from and write into. Arrays usvr and uwrk can be used for data storage. Data will be saved between multiple calls made to user_call between the Start and Finish stages of the element calculations (Figure 2.4: Element Calculation Workflow). Data stored in the usvr array is also written to the .esav file.

Example 2.4: Example Using elem Data Arrays in user_call

The following element properties (Table 2.9: Element Coordinates and Shape Data) require you to specify the location where the property is calculated:

elem%item(loc)
where:
item - name of the element property
loc - location within the element:
none = nodes
0 = centroid
ip = integration point number
(s,t,r) = natural coordinates of a point within the element

Table 2.9: Element Coordinates and Shape Data

*`Item`*	Description	Location	Output Type
xyz^[a]	X, Y, Z coordinates	nodes	dp ar(3,num%nodes)
		centroid	dp ar(3)
		integration point	dp ar(3)
		(s,t,r)	dp ar(3)
shpF	shape functions	centroid integration point (s,t,r)	dp ar (num%nodes)
shpFd^[a]	shape function derivatives	centroid integration point (s,t,r)	dp ar (num%dim,num%nodes)
vol^[a]	volume	integration point	dp sc

^[a]If large deflections are active (NLGEOM,ON) in the coupled-field analyses with structural DOFs, element coordinates, volume, and shape function derivatives with respect to X, Y, Z coordinates are updated based on the current displacements.

Example 2.5: Getting Element Coordinates and Shape Function Data

You can also perform the following operations on a set of scalar nodal values.

Midside node interpolation
call elem%midnod(nodes,var)
where:
nodes – list of element nodes with missing midside nodes (in)
var – array of nodal variables to be interpolated (in/out)
This subroutine modifies the array of nodal variables (var) by filling up the missing midside node entries with interpolated values.
Table 2.10: Element Midside Node Interpolation
Name Description Input Type Input / Output Type
midnod midside node interpolation int ar(num%nodes) dp ar(num%nodes)

The following example mimics a heat generation body load defined on the outer radius of a circular 2-D model.
Example 2.6: Interpolating Element Heat Generation Load
Gradient calculation
elem%grad(var,loc)
where:
var – array of nodal variables,
loc – location within the element:
0 = centroid
ip = integration point number
(s,t,r) = natural coordinates of a point within the element
The resulting array contains the derivatives of the nodal variable var with respect to X, Y, and Z in the element coordinate system.
Table 2.11: Element Gradient Operation
Name Description Input Type Location Output Type
grad gradient dp ar(num%nodes) centroid dp ar(num%nodes)
integration point
(s,t,r)

Example 2.7: Calculating Thermal Gradient at Different Element Locations

Name	Description	Input Type	Input / Output Type
midnod	midside node interpolation	int ar(num%nodes)	dp ar(num%nodes)

Name	Description	Input Type	Location	Output Type
grad	gradient	dp ar(num%nodes)	centroid	dp ar(num%nodes)
integration point
(s,t,r)

You can access nodal and element results using the elem type-bound procedures (see Table 2.12: Element Results).

elem%item(loc)
where:
item - name of the element result
loc - location within the element:
none = nodes
0 = centroid
ip = integration point number
(s,t,r) = natural coordinates of a point

Table 2.12: Element Results

Name Description	Location	Output
u displacements	nodes	dp ar(num%dim,num%nodes)
	centroid	dp ar(num,%dim)
	integration point	dp ar(num,%dim)
	(s,t,r)	dp ar(num,%dim)

temp temperature volt electric potential conc concentration	nodes	dp ar(num%nodes)
	centroid	dp sc
	integration point	dp sc
	(s,t,r)	dp sc

eptt total strain epel elastic strain epth thermal strain epdi diffusion strain eppl plastic strain s stress	nodes	dp ar(num%compn,num%nodes)
	centroid	dp ar(num%compn)
	integration point	dp ar(num%compn)
	(s,t,r)	dp ar(num%compn)

tg thermal gradient tf thermal flux ef electric field jc electric current cg concentration gradient df diffusion flux	nodes	dp ar(num%comp,num%nodes)
	centroid	dp ar(num%comp)
	integration point	dp ar(num%comp)
	(s,t,r)	dp ar(num%comp)

Example 2.8: Obtaining Nodal and Element Results

You can print element results in a formatted way using the elem type-bound print_ procedures (see Table 2.13: Print Element Results).

call elem%print_item(loc)
where:
item - name of the element result
loc - location within the element:
none = nodes
0 = centroid
ip = integration point number
(s,t,r) = natural coordinates of a point

Note: elem%print_ routines are not supported for simulations using distributed-memoy parallel (DMP) processing.

Table 2.13: Print Element Results

`Item`	Description	Location
u	displacement	nodes centroid integration point (s,t,r)
temp volt conc	temperature electric potential concentration
eptt epel epth epdi eppl s	total strain elastic strain thermal strain diffusion strain plastic strain stress
tg tf ef jc cg df	thermal gradient thermal flux electric field electric current concentration gradient diffusion flux
smisc	summable miscellaneous	-
nmisc	non-summable miscellaneous	-

Example 2.9: Printing Nodal and Element Results

Printout for element 20:

To define or modify materials, loads, constitutive equations, and output quantities, use the set_ and add_ elem type-bound procedures, respectively.

To assign a value var to the element property item, use:
call elem%set_item_level(var)
To modify the element property item by adding the value var to it, use:
call elem%add_item_level(var)
where:
item - name of the user-programmable element quantity
level - level of customization:
p = data prep
ip = integration point
o = output
var - variable that overwrites (set_) or is being added (add_) to the element quantity item

User-programmable element quantities that can be defined or modified by the set_ and add_ elem type-bound procedures are listed in the following table.

Table 2.14: User-programmable Element Quantities

*`Item`*	Description	Level	Variable Type
tcond	thermal conductivity matrix	p	dp ar(num%comp,num%comp)
tcond	thermal conductivity matrix	ip	dp ar(num%comp,num%comp)
hgen	heat generation rate	p	dp ar(num%nodes)
hgen	heat generation rate	ip	dp sc
hflux	heat flux surface load	p	dp ar(num%fluxes)
tf	thermal flux	ip	dp ar(num%comp)

econd	electric conductivity matrix	p	dp ar(num%comp,num%comp)
econd	electric conductivity matrix	ip	dp ar(num%comp,num%comp)
jc	integration point current density	ip	dp ar(num%comp)

diffu	diffusivity matrix	p	dp ar(num%comp,num%comp)
diffu	diffusivity matrix	ip	dp ar(num%comp,num%comp)
csat	saturated concentration	p	dp sc
csat	saturated concentration	ip	dp sc
dgen	heat generation rate	p	dp ar(num%nodes)
dgen	heat generation rate	ip	dp sc
dflux	diffusion flux surface load	p	dp ar(num%fluxes)
df	diffusion flux	ip	dp ar(num%comp)

smisc^[a]	summable miscellaneous	o	dp ar(num%smisc)
nmisc^[a]	non-summable miscellaneous	o	dp ar(num%nmisc)

^[a]Only the set_ procedure can be used with the smisc and nmisc items.

Example 2.10: Defining Integration Point Thermal Conductivity as a Linear Function of Temperature

When modifying the user-programmable element quantities using the set_ and add_ subroutines (Table 2.14: User-programmable Element Quantities), you may need to access the nodal DOF and element calculated results. While the nodal DOF results are always available at the start of the element calculation, the element calculated results such as strains, fields, fluxes, and miscellaneous items are calculated in a sequence shown on the left of Figure 2.6: Element Results Calculation vs. User-Programmable Quantities Workflow. Therefore, element calculated results may or may not be ready at the time when a set_ or add_ subroutine is called to update a user-programmable quantity as shown on the right of Figure 2.6: Element Results Calculation vs. User-Programmable Quantities Workflow.

Figure 2.6: Element Results Calculation vs. User-Programmable Quantities Workflow

For example, when customizing thermal quantities such as thermal conductivity (tcond) or heat generation rate (hgen) in a structural-thermal-diffusion analysis, you can access structural results such as stress and elastic strain using elem%s(loc) and elem%epel(loc) functions, respectively, because the standard structural results are calculated before the thermal quantities are customized. On the other hand, standard diffusion results such as concentration gradient and diffusion flux are calculated after the customization of thermal properties, and results reported by elem%cg(loc) and elem%df(loc) during thermal calculations may not be up to date.

When a needed standard element result has not yet been calculated, you can independently derive it from the DOF solution using the element shape functions (Table 2.9: Element Coordinates and Shape Data) or the gradient operation (Table 2.11: Element Gradient Operation).

In the following example, both the thermal flux (tf) and the electric current density (jc) are customized at the integration point level to mimic the Seebeck-Peltier effects without specifying MP,SBKX that would invoke the standard coupled-field effect (Thermoelectrics in the Theory Reference). The integration point (i) thermal gradient needed for the calculation of the Seebeck electric current density (jc_usr) is already available as elem%tg(i) when elem%add_jc_ip(jc_usr) is called and does not need to be recalculated. On the other hand, the Peltier thermal flux calculation needs the total electric current density that will become available only in the next (electric) stage of the element workflow (Figure 2.6: Element Results Calculation vs. User-Programmable Quantities Workflow). To overcome this limitation, you can independently calculate the total current density (jc_usr) needed for the Peltier heat contribution to thermal flux as shown in the example below.

Example 2.11: Customizing Thermoelectric Constitutive Equations

Note that the set_ and add_ procedures are designed so that the respective element variable is updated only when the program is ready to accept it. Therefore, the sequence in which these subroutines are called within user_call is irrelevant.

Solution Settings and Controls (solu)

The structure solu contains frequently used solution settings that you can use to control the calculations in user_call. They can be accessed as:

solu%item
where:
item – name of the solution parameter

Table 2.15: Solution Settings

Item Description Output Type Output Value

antype analysis type int sc ANTYPE command

nlgeom

geometric nonlinearities

int sc

0 – off

1 – on

nropt

Newton-Raphson option

int sc

NROPT command:

0 – auto

1 – full or unsym

2 – modi

3 – init

nrunsym

Newton-Raphson option set to UNSYM

int sc

0 – default

1 – if NROPT,UNSYM

xtrap

extrapolation of integration point results to element corner nodes (ERESX command)

int sc

0 – copy

1 – extrapolate except in elements with active plasticity

2 – extrapolate

timint

transient effects (TIMINT command)

int sc

0 – off

1 – on

ldstep current load step number int sc

isubst current substep number int sc

nsubst maximum number of substeps int sc NSUBST command setting

ieqitr current iteration number int sc

neqitr maximum number of iterations int sc NEQIT command setting

output

element output pass (result calculation)

int sc

0 – off

1 – on

svrupd

history (state) variables update

int sc

0 – not updated

1 – updated

converged

converged solution

int sc

0 – non converged

1 – converged

time%beg time at the beginning of the current step dp sc

time%end time at the end of the current load step dp sc TIME command setting

time%cur current time value dp sc

time%inc current time increment dp sc

time%prv previous time value dp sc

time%ino previous time increment dp sc

toffst temperature offset from absolute zero to zero dp sc TOFFST command setting

Example 2.12: Using Solution Information

Scope of calculations (scope)

The user_call procedure is called multiple times from different parts of the element code, including every integration point of each physics (Figure 2.4: Element Calculation Workflow). It may be inefficient to repeat the same calculations every time the user_call procedure is called. The components of the scope data structure allow you to narrow down the calculations to a specific location in the element workflow.

Using the scope variables, you can limit your code to the physics, data preparation or integration loop of that physics, to a specific element quantity, or to the element output calculations:

scope%calc
scope%calc_level
scope%calc_level_item
where:
calc - area of element calculations
level - level of customization
prep = data prep
intp = integration point
item - name of the user-programmable element quantity

All variables in the scope data structure are of logical type.

Table 2.16: Parameters to Define the Scope of User Calculations

Name	Description	Type
thermal	thermal calculations	log sc
thermal_prep	thermal data prep
thermal_prep_tcond thermal_prep_hgen thermal_prep_hflux	input tcond, hgen, hflux calculation
thermal_intp	thermal integration loop
thermal_intp_hgen thermal_intp_tcond thermal_intp_tf	integration point tcond, hgen, tf calculation

electric	electric calculations	log sc
electric_prep	electric data prep
electric_prep_econd	input econd calculation
electric_intp	electric integration loop
electric_intp_econd electric_intp_jc	integration point tcond, jc

diffusion	diffusion calculations	log sc
diffusion_prep	diffusion data prep
diffusion_prep_diffu diffusion_prep_csat diffusion_prep_dgen diffusion_prep_dflux	input diffu, csat, dgen, dflux calculation
diffusion_intp	diffusion integration loop
diffusion_intp_diffu diffusion_intp_csat diffusion_intp_dgen diffusion_intp_df	integration point diffu, csat, dgen, df calculation

output	output calculations	log sc
output_smisc output_nmisc	smisc and nmisc calculation	log sc

A more efficient version of the user_call Example 2.11: Customizing Thermoelectric Constitutive Equations is shown in the following example.

Example 2.13: Using a Scope Variable in user_call

2.4.4. Example: User implementation of stress- and electro-migration analysis

This example demonstrates the coupled-field user-programmable features (UPFs) by reproducing the standard element calculations performed when the migration model (TB,MIGR) is specified The coupled-field UPFs are activated by setting KEYOPT(12) to 1 on a coupled-field element with diffusion degrees-of-freedom. The diffusion and migration parameters (D0, Ea/k, V/k, and Ze/k) are input as real constants instead of TB,MIGR that triggers the standard calculations (Figure 2.7: Standard vs. User-programmable APDL Input).

Figure 2.7: Standard vs. User-programmable APDL Input

The following calculations are performed in UserProcedures.F:

the integration point diffusivity matrix D is customized using the elem%set_diffu_ip procedure
the diffusion flux constitutive equation is customized to include stress-migration flux (df_s) and electro-migration flux (df_e), and added to the total diffusion flux using the elem%add_df_ip procedure
The summable miscellaneous items corresponding to the pure diffusion flux (DFCX, DFCY, DFCZ,DFCSUM ), stress-migration flux (DFSX, DFSY, DFSZ, DFSSUM), and electro-migration flux (DFEX, DFEY, DFEZ,DFESUM), respectively, are calculated at the element centroid and output using the elem%set_smisc_o procedure

*deck,UserProcedures       USERDISTRIB       parallel             
      submodule (UserInterface) UserMigr

#include "impcom.inc"
              
      contains

        module procedure user_char
           num%rcon = 5           ! number of real constants
           non%lin = 1            ! nonlinear analyses
        end procedure user_char

        module procedure user_call
           integer :: i
           double precision :: D0,Ea_k,V_k,Z_k,C,T,Dt,
     x       D(num%dim,num%dim),
     x       stress(num%compn,num%nodes),sigH(num%nodes),
     x       df_c(num%comp),df_s(num%comp),df_e(num%comp),
     x       df_s_m,df_e_m,df_c_m,smisc(num%smisc)

!      --- get migration parameters as real constants
           D0 = elem%rcon(2)      ! initial diffusivity coefficient
           Ea_k = elem%rcon(3)    ! activation energy/k 
           V_k = elem%rcon(4)     ! atomic volume/k
           Z_k = elem%rcon(5)     ! charge number/k

!      --- nodal stress
           stress = elem%s()
!      --- hydrostatic stress
           sigH = (stress(1,:) + stress(2,:) + stress(3,:))/3. 

!      --- integration point: 
           if (scope%diffusion_intp) then
              i = elem%ip()          ! number
              C = elem%conc(i)       ! concentration
              T = elem%temp(i) + solu%toffst   ! temperature
              Dt = D0*exp(-Ea_k/T)   ! diffusivity
              D = 0.0d0              ! diffusivity matrix
              D(1,1) = Dt; D(2,2) = Dt; D(3,3) = Dt

!         --- set integration point diffusivity
              call elem%set_diffu_ip(D)
           
!         --- hydrostatic stress gradient driven migration
              df_s = (V_k*C/T)*matmul(D,elem%grad(sigH,i))

!         --- electric field driven migration
              df_e = (Z_k*C/T)*matmul(D,elem%ef(i))

!         --- add df_s to the total diffusion flux
              call elem%add_df_ip(df_s)

!         --- add df_e to the total diffusion flux
              call elem%add_df_ip(df_e)
           end if

!      --- centroidal calculations
           if (scope%output) then
              C = elem%conc(0)       ! concentration
              T = elem%temp(0) + solu%toffst   ! temperature
              Dt = D0*exp(-Ea_k/T)   ! diffusivity

              df_c = -Dt*elem%grad(elem%conc(),0)
              df_s = Dt*V_k*C*elem%grad(sigH,0)/T
              df_e = Dt*Z_k*C*elem%ef(0)/T

              df_c_m = sqrt(df_c(1)**2 + df_c(2)**2 + df_c(3)**2)
              df_s_m = sqrt(df_s(1)**2 + df_s(2)**2 + df_s(3)**2)
              df_e_m = sqrt(df_e(1)**2 + df_e(2)**2 + df_e(3)**2)

!         --- save magnitudes of diffusion fluxes as nmiscs
              smisc = elem%smisc()    ! copy standard smisc
!         --- store diffusion flux
              smisc(2) = df_c(1)     
              smisc(3) = df_c(2)     
              if (num%comp==3) smisc(4) = df_c(3)
              smisc(5) = df_c_m
!         --- store stress migration flux
              smisc(6) = df_s(1)
              smisc(7) = df_s(2)
              if (num%comp==3) smisc(8) = df_s(3)
              smisc(9) = df_s_m
!         --- store electric migration flux
              smisc(14) = df_e(1)
              smisc(15) = df_e(2)
              if (num%comp==3) smisc(16) = df_e(3)
              smisc(17) = df_e_m

              call elem%set_smisc_o(smisc)
              if (solu%isubst==10) call elem%print_smisc()
           end if

        end procedure user_call

      end submodule UserMigr