discretelement_mod.f

Go to the documentation of this file.

! -*- f90 -*-
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!   Module name: DISCRETELEMENT                                        C
!   Purpose: DES mod file                                              C
!            Common Block containing DEM conditions                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      MODULE discretelement

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param, only: dim_m
      IMPLICIT NONE
!-----------------------------------------------

! Define interface - needed when passing arrays of assumed size
! This function is used to identify fluid grid index (i,j or k) of a
! particle in a given x-, y-, z- line of space
      INTERFACE
         INTEGER FUNCTION des_getindexfrompos(LIM1,LIM2,PART_POS,&
            grid_pos,axis,axis_index)
            INTEGER :: lim1, lim2
            DOUBLE PRECISION :: part_pos
            DOUBLE PRECISION, DIMENSION(:) :: grid_pos
            CHARACTER(LEN=1) :: axis,axis_index
         END FUNCTION des_getindexfrompos
      END INTERFACE

      INTEGER, PARAMETER :: dim_m_tri = 55 ! 10th triangular number

! Total number of particles in simulation: read from input or generated
      INTEGER particles

! Start particle tracking quantities
!----------------------------------------------------------------->>>
! Generally for inlet/outlet related routines but also employed in
! tracking for parallelization

! Dynamic particle states:

! NONEXISTENT: This state is used with the inlet/outlet to skip
! indices that do not represent particles in the system or indices
! that represent particles that have exited the system.

! ENTERING: This state identifies a particle as 'new' if true.
! Particles with a classification of 'new' do not react when in contact
! with a wall or another particle, however existing particles do collide
! and interact with 'new' particles. The classification allows new
! particles to push particles already in the system out of the way when
! entering to prevent overlap.  This flag is also used when the center
! of a particle crosses a dem outlet (i.e. an exiting particle; see
! EXITING) so that the particle will maintain its present trajectory
! until it has fully exited the system

! EXITING: This state identifies a particle as 'exiting' if true.
! If a particle initiates contact with a wall surface designated as a
! des outlet, this flag is set to true. With this classification the
! location of the particle is checked to assess if the particle has
! fully exited the system.  At this point, the particle is removed
! from the list.

! GHOST, ENTERING_GHOST, EXITING_GHOST: for ghost particles

      INTEGER(KIND=1), DIMENSION(:), ALLOCATABLE :: particle_state ! (PARTICLES)

      INTEGER, PARAMETER :: nonexistent=0
      INTEGER, PARAMETER :: normal_particle=1
      INTEGER, PARAMETER :: entering_particle=2
      INTEGER, PARAMETER :: exiting_particle=3
      INTEGER, PARAMETER :: normal_ghost=4
      INTEGER, PARAMETER :: entering_ghost=5
      INTEGER, PARAMETER :: exiting_ghost=6

! PARALLEL PROCESSING: explanation of variables in parallel architecture
! pip - particles in each processor (includes the ghost particles)
! max_pip - maximum allocated particles in processor

! Number of particles in the system (current)
      INTEGER pip
! Global sum of particles (excluding ghost) in the system
      INTEGER tot_par
! Maximum particles permitted in the system at once
      INTEGER max_pip

! End particle tracking quantities
!-----------------------------------------------------------------<<<

! For parallel processing: global id of particles
      INTEGER, DIMENSION(:), ALLOCATABLE :: iglobal_id
! Ghost count of particles on each processor
      INTEGER :: ighost_cnt
! Maximum global id, new particles global id will be assigned based
! on this value
      Integer :: imax_global_id


! If gener_part_config is true, then the particle_input.dat file does
! not need to be supplied nor does the total number of particles as
! these are determined based on the specified volume fraction (vol_frac)
! in the specified domain (des_eps_xyzstart)
      LOGICAL :: gener_part_config
      DOUBLE PRECISION ::  vol_frac(dim_m)
      DOUBLE PRECISION :: des_eps_xstart, &
      des_eps_ystart, des_eps_zstart
! volume of the IC region for computing number of particles to be seeded
      DOUBLE PRECISION, dimension(:), allocatable :: vol_ic_region!(DIMENSION_IC)
! number of particles for each phase corresponding to the IC number. This will
! be real particles for DEM but parcels or computational particles for PIC model
      INTEGER, dimension(:,:), allocatable :: part_mphase_byic!(DIMENSION_IC, DIM_M)

! Number of real particles by IC and by solid phase. Only relevant for PIC model
      double precision, dimension(:,:), allocatable :: realpart_mphase_byic!(DIMENSION_IC, DIM_M)
! The number of particles that belong to solid phase M according to the
! vol_frac and particle diameter. this information is used when
! gener_part_config is invoked for initialization
! This will be removed soon as PART_MPHASE_BYIC will be used from now on
      INTEGER part_mphase(dim_m)

! The number of real particles that belong to solid phase M during the initialization.
! It is equal to Part_mphase for DEM but implies real number of particles for PIC model
      double precision realpart_mphase(dim_m)
! Assigns the initial particle velocity distribution based on user
! specified mean and standard deviation (regardless if already set
! within particle_input.dat)
      DOUBLE PRECISION pvel_mean, pvel_stdev

! Output/debug controls
!----------------------------------------------------------------->>>
! Logic that controls whether to print data dem simulations (granular or
! coupled)
      LOGICAL print_des_data
      CHARACTER(LEN=255) :: vtp_dir

! logic that controls if des run time messages are printed on screen or not
      LOGICAL print_des_screen

! This specifies the file type used for outputting DES data
! options are :
!    TECPLOT - data is written in Tecplot format
!    undefined - data is written in ParaView format (default)
      CHARACTER(LEN=64) :: des_output_type

! Used sporadically to control screen dumps (for debug purposes)
      LOGICAL :: debug_des

! Single particle no. index that is followed if debugging
      INTEGER focus_particle

! Output file count for .vtp type files and for tecplot files;
! for vtp output used to label .vtp file names in sequential order
! and is saved so restarts begin at the correct count
      INTEGER vtp_findex, tecplot_findex
! End Output/debug controls
!-----------------------------------------------------------------<<<

! DES - Continuum
      LOGICAL discrete_element
      LOGICAL des_continuum_coupled

! DES - Invoke hybrid model where both the DEM and continuum model
! are employed to describe solids
      LOGICAL des_continuum_hybrid

! DES -
! With this logic the particles see the fluid but the fluid does
! not see the particles.
      LOGICAL des_oneway_coupled

! Only used when coupled and represents the number of times a pure
! granular flow simulation is run before the actual coupled simulation
! is started (i.e. for particle settling w/o fluid forces)
      INTEGER nfactor

! Drag
      LOGICAL tsuji_drag

! Collision model, options are as follows
!   linear spring dashpot model (default/undefined)
!   'hertzian' model
      CHARACTER(LEN=64) :: des_coll_model
      INTEGER des_coll_model_enum
      INTEGER,PARAMETER ::  hertzian=0
      INTEGER,PARAMETER ::  lsd=1

! Integration method, options are as follows
!   'euler' first-order scheme (default)
!   'adams_bashforth' second-order scheme (by T.Li)
      CHARACTER(LEN=64) :: des_intg_method
      LOGICAL intg_adams_bashforth
      LOGICAL intg_euler

! Value of solids time step based on particle properties
      DOUBLE PRECISION dtsolid
! Run time value of simulation time used in dem simulation
      DOUBLE PRECISION s_time


! Neighbor search related quantities
!----------------------------------------------------------------->>>
! Neighbor search method, options are as follows
!   1= nsquare, 2=quadtree, 3=octree, 4=grid/cell based search
      INTEGER des_neighbor_search

! Quantities used to determine whether neighbor search should be called
      INTEGER neighbor_search_n
      DOUBLE PRECISION neighbor_search_rad_ratio
      LOGICAL do_nsearch

! Flag on whether to have DES_*_OLD arrays, if either Adams Bashforth or PIC is used
      LOGICAL do_old

! Factor muliplied by sum of radii in grid based neighbor search and
! nsquare search method.  increases the effective radius of a particle
! for detecting particle contacts
      DOUBLE PRECISION factor_rlm

! Stores number of neighbors based on neighbor search
      INTEGER, DIMENSION(:), ALLOCATABLE :: neighbor_index
      INTEGER, DIMENSION(:), ALLOCATABLE :: neighbor_index_old
      INTEGER, DIMENSION(:), ALLOCATABLE :: neighbors
      INTEGER, DIMENSION(:), ALLOCATABLE :: neighbors_old
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: pft_neighbor
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: pft_neighbor_old

      INTEGER :: neigh_num

! Quantities used for reporting: max no. neighbors and max overlap
! that exists during last solid time step of dem simulation
      DOUBLE PRECISION overlap_max

! The number of i, j, k divisions in the grid used to perform the
! cell based neighbor search
      INTEGER :: desgridsearch_imax, desgridsearch_jmax, &
                 desgridsearch_kmax

! End neighbor search related quantities
!-----------------------------------------------------------------<<<

! User specified dimension of the system (by default 2D, but if 3D system is
! desired then it must be explicitly specified)
      INTEGER, PARAMETER :: dimn = 3

! Variable that is set to the number of walls in the system
      INTEGER nwalls

! Position of domain boundaries generally given as
!   (0, xlength, 0, ylength, 0, zlength)
      DOUBLE PRECISION wx1, ex2, by1, ty2, sz1, nz2

! X, Y, Z position of cell faces of computational fluid grid
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xe  !(0:DIMENSION_I)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: yn  !(0:DIMENSION_J)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: zt  !(0:DIMENSION_K)

! Wall normal vector
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: wall_normal  !(NWALLS,3)

! Gravity vector and magnitude
      DOUBLE PRECISION :: grav(3)
      DOUBLE PRECISION :: grav_mag


! Periodic wall BC
      LOGICAL des_periodic_walls
      LOGICAL des_periodic_walls_x
      LOGICAL des_periodic_walls_y
      LOGICAL des_periodic_walls_z


! Lees & Edwards wall BC (lost in current DEM)
!----------------------------------------------------------------->>>
! Logic for Lees & Edwards BC (T = turn on LE BC)
      LOGICAL des_le_bc
! Relative velocity of LE boundaries (distance/time)
      DOUBLE PRECISION des_le_rel_vel
! Shear direction
!   2D options are DUDY or DVDX
!   3D options are DUDY, DUDZ, DVDX, DVDZ, DWDX or DWDY
!   Note that all other directions are treated as periodic boundaries
      CHARACTER(LEN=4) :: des_le_shear_dir
! End LE BC
!-----------------------------------------------------------------<<<



! Particle-particle and Particle-wall collision model parameters
!----------------------------------------------------------------->>>
! Spring contants
      DOUBLE PRECISION kn, kn_w  !Normal
      DOUBLE PRECISION kt, kt_w, kt_fac, kt_w_fac  ! Tangential factors = KT/KN and KT_w/KN_w, resp.

! Damping coeffients
      DOUBLE PRECISION eta_des_n, eta_n_w  !Normal
      DOUBLE PRECISION eta_des_t, eta_t_w  !Tangential

! Flag to use van der Hoef et al. (2006) model for adjusting the rotation of the
! contact plane
      LOGICAL :: use_vdh_dem_model
! Tangential damping factors, eta_t = eta_t_factor * eta_N
      DOUBLE PRECISION des_etat_fac, des_etat_w_fac

! Damping coeffients in array form
      DOUBLE PRECISION :: des_etan(dim_m, dim_m), des_etan_wall(dim_m)
      DOUBLE PRECISION :: des_etat(dim_m, dim_m), des_etat_wall(dim_m)

! Friction coeficients
      DOUBLE PRECISION mew, mew_w

! coeff of restituion input in one D array, solid solid
! Tangential rest. coef. are used for hertzian collision model but not linear
      DOUBLE PRECISION des_en_input(dim_m_tri)
      DOUBLE PRECISION des_et_input(dim_m_tri)

! coeff of restitution input in one D array, solid wall
      DOUBLE PRECISION des_en_wall_input(dim_m)
      DOUBLE PRECISION des_et_wall_input(dim_m)

! Hertzian collision model:
      DOUBLE PRECISION :: e_young(dim_m), ew_young
      DOUBLE PRECISION :: v_poisson(dim_m), vw_poisson
      DOUBLE PRECISION :: hert_kn(dim_m, dim_m), hert_kwn(dim_m)
      DOUBLE PRECISION :: hert_kt(dim_m, dim_m), hert_kwt(dim_m)

! Realistic material properties for correcting contact areas:
! Actual simulations will probably use LSD with softened spring
! or hertzian with unrealistically small E_Young
      DOUBLE PRECISION :: e_young_actual(dim_m), ew_young_actual
      DOUBLE PRECISION :: v_poisson_actual(dim_m), vw_poisson_actual
      DOUBLE PRECISION :: hert_kn_actual(dim_m, dim_m), hert_kwn_actual(dim_m)

!     Coefficients for computing binary (Hertzian) collision time.  The actual
!     collision time is TAU_C_Base_Actual * V^(-1/5).
      DOUBLE PRECISION :: tau_c_base_actual(dim_m, dim_m), tauw_c_base_actual(dim_m)
      DOUBLE PRECISION :: tau_c_base_sim(dim_m, dim_m), tauw_c_base_sim(dim_m)

! End particle-particle and particle-wall collision model parameters
!-----------------------------------------------------------------<<<


! Particle attributes: radius, density, mass, moment of inertia
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: des_radius !(PARTICLES)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: ro_sol     !(PARTICLES)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: pvol       !(PARTICLES)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: pmass      !(PARTICLES)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: omoi       !(PARTICLES)

! Additional quantities
      DOUBLE PRECISION :: min_radius, max_radius

! number of discrete 'solids phases'
      INTEGER des_mmax

! Old and new particle positions, velocities (translational and
! rotational)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: des_pos_old  !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: des_pos_new  !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: des_vel_old  !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: des_vel_new  !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: omega_old    !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: omega_new    !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: ppos         !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: des_acc_old  !(PARTICLES,3)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: rot_acc_old  !(PARTICLES,3)

! Particle orientation
      LOGICAL :: particle_orientation = .false.
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: orientation  !(3,PARTICLES)
      DOUBLE PRECISION, DIMENSION(3) :: init_orientation = (/0.0,1.0,0.0/)

! Defining user defined allocatable array
      INTEGER :: des_usr_var_size
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: des_usr_var  !(PARTICLES,3)

! Total force and torque on each particle
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: fc    !(3,PARTICLES)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: tow   !(3,PARTICLES)

!     particle can collide with at most COLLISION_ARRAY_MAX facets simultaneously
      INTEGER :: collision_array_max = 8

! (COLLISION_ARRAY_MAX,PARTICLES)
!     -1 value indicates no collision
      INTEGER, ALLOCATABLE :: wall_collision_facet_id(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: wall_collision_pft(:,:,:)

! Store the number of particles in a computational fluid cell
      INTEGER, DIMENSION(:), ALLOCATABLE :: pinc  ! (DIMENSION_3)

! For each particle track its i, j, k & ijk location on the fluid grid
! and solids phase no.:
      INTEGER, DIMENSION(:,:), ALLOCATABLE :: pijk ! (PARTICLES,5)=>I,J,K,IJK,M
!-----------------------------------------------------------------<<<


! note that thse variables are needed since the existing variables (i.e.
! f_gs, f_ss, etc) are also being used to store the information between
! the gas and continuous solids phases.
! drag coefficient between gas phase and discrete particle 'phases'
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: f_gds
! drag coefficient between continuous solids phases and discrete
! particle 'phases'
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: f_sds

! the following should probably be local to the subroutine
! solve_vel_star they are only needed when invoking the non-interpolated
! version of drag wherein they are used to determine part of the source
! in the gas-phase momentum balances, in the routine to determine
! coefficients for the pressure correction equation and in the partial
! elimination algorithm
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: vxf_gds
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: vxf_sds

! the contribution of solids-particle drag to the mth phase continuum
! solids momentum A matrix (center coefficient)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: sdrag_am
! the contribution of solids-particle drag to the to mth phase continuum
! solids momentum B vector
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: sdrag_bm


! the coefficient add to gas momentum A matrix
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: drag_am
! the coefficient add to gas momentum B matrix
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: drag_bm

! Explicitly calculated fluid-particle drag force.
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: drag_fc

! An intermediate array used in calculation of mean solids velocity
! by backward interpolation, i.e., when INTERP_DES_MEAN_FIELDS is true.
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE ::des_vel_node
                        !(DIMENSION_3,3,DIMENSION_M)

! An intermediate array used in calculation of solids volume fraction
! by backward interpolation, i.e., when INTERP_DES_MEAN_FIELDS is true.
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE ::  des_rops_node
                        !(DIMENSION_3,DIMENSION_M)

      DOUBLE PRECISION, DIMENSION(:,:,:), POINTER :: weightp

! the gas-particle drag coefficient
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: f_gp
                        !(PARTICLES)
      DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: wtderivp
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: sstencil
      DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: gstencil, vstencil, pgradstencil

! stencil for interpolation of solids pressure
      DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE ::  psgradstencil

! stencil for interpolation of solids velocity
      DOUBLE PRECISION, DIMENSION(:,:,:,:,:), ALLOCATABLE::  vel_sol_stencil


! quantities are set in subroutine set_interpolation_scheme
! order = order of the interpolation method, ob2l = (order+1)/2,
! ob2r = order/2
      CHARACTER(LEN=7):: scheme, interp_scheme
      INTEGER:: order, ob2l, ob2r

! END of interpolation related data
!-----------------------------------------------------------------<<<

      ! Volume of each node. Used to obtain Eulerian fields
      double precision, allocatable, dimension(:) :: des_vol_node
                        !(DIMENSION_3,dimn,dimension_m)

! Variable to track pressure force in computational fluid cell (ijk)
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: p_force
                        !(DIMN, DIMENSION_3)

! Granular temperature in a fluid cell
! Global average velocity: obtained by averaging over all the particles
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: des_vel_avg
                        !(3)

! Global granular energy & temp: obtained by averaging over all the particles
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: global_gran_energy
                        !(3)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: global_gran_temp
                        !(3)

! Kinetic and potential energy of the system: obtained by averaging
! over all particles
! Added rotational kinetic energy (DES_ROTE)
      DOUBLE PRECISION des_ke, des_pe, des_rote

! Logic for bed height calculations (T = turn on bed height
! calculations)
      LOGICAL des_calc_bedheight
! Used to track bed height of solids phase M
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: bed_height
                       !(dimension_m)

! MAX velocity of particles in each direction
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: des_vel_max


! Flag to turn on/off optimizing the list of facets at each des grid cell

      LOGICAL :: minimize_des_facet_list
!-----------------------------------------------------------------<<<


! Start Cohesion
!----------------------------------------------------------------->>>
! Includes square-well type model and a van der waals type model

! Switch to turn cohesion on and off (set in mfix.dat)
      LOGICAL use_cohesion


! Van der Waals constants (set in mfix.dat)
      LOGICAL square_well, van_der_waals
      DOUBLE PRECISION hamaker_constant
      DOUBLE PRECISION vdw_inner_cutoff ! (in cm)
      DOUBLE PRECISION vdw_outer_cutoff
      DOUBLE PRECISION wall_hamaker_constant
      DOUBLE PRECISION wall_vdw_inner_cutoff
      DOUBLE PRECISION wall_vdw_outer_cutoff
      DOUBLE PRECISION asperities ! average radius of asperities (default zero)

! Store postcohesive
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: postcohesive
                        !(PARTICLES)
! Store cluster information array for postprocessing
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: postcluster

! Variables for van der waals cohesion calculations:
! Surface energy used to calculate cohesive force for low separation distances
! in Van der Waals model (this variable is calculated at the beginning of each
! simulation to ensure the van der Waals force is continuous at the inner cutoff)
      DOUBLE PRECISION surface_energy
      DOUBLE PRECISION wall_surface_energy

! END Cohesion
!-----------------------------------------------------------------<<<

      LOGICAL :: des_explicitly_coupled

      integer, dimension(:),allocatable :: dg_pijk,dg_pijkprv

! variable to clean the ghost cells
      logical,dimension(:),allocatable :: ighost_updated
      integer :: max_isize

      CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
!                                                                     !
!  Subroutine: DES_CROSSPRDCT                                         !
!  Purpose: Calculate the cross product of two 3D vectors.            !
!                                                                     !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      FUNCTION des_crossprdct(XX,YY)

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Input vectors
      DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: xx, yy
! Result: cross product of vectors
      DOUBLE PRECISION, DIMENSION(3) :: des_crossprdct
!......................................................................!

      des_crossprdct(1) = xx(2)*yy(3) - xx(3)*yy(2)
      des_crossprdct(2) = xx(3)*yy(1) - xx(1)*yy(3)
      des_crossprdct(3) = xx(1)*yy(2) - xx(2)*yy(1)

      END FUNCTION des_crossprdct

      END MODULE discretelement