// CLUSTERS.C
// Cluster trace functions
// © 2021 Peter J. Meyer
#include "iss.h"
#define STACK_PTR_CHECK true
// These are global.
// int ct_spin_value, ct_new_spin_value;
// short int ct_i, ct_j, ct_k, ct_l, ct_r;
// short int nn_i, nn_j, nn_k, nn_l;
// int ct_size_of_cluster;
// short unsigned int ct_cluster_number;
// int ct_cluster_number_overflow;
// Variables local to this module.
static int ct_cluster_type;
static int ct_number_clusters;
static int ct_stop_on_spanning_cluster;
static int ct_get_cluster_size;
static int ct_spanning_cluster_found;
static int ct_num_spanning_clusters;
static short int best_r;
static short unsigned int ct_spin_distance;
static void trace_single_cluster(void);
static void clear_cluster_number_array(void);
static void clear_path_status(void);
static void update_path_status(void);
static int is_spanning_cluster(int total_spanning_cluster);
// Call with input data set up in *actd.
/*---------------------------------------------------------*/
void trace_all_clusters(struct _all_cluster_trace_data *actd)
{
int i, j, k, l;
ct_cluster_type = actd->cluster_type;
ct_number_clusters = actd->number_clusters;
ct_stop_on_spanning_cluster = actd->stop_on_spanning_cluster;
ct_get_cluster_size = actd->get_cluster_size;
if ( ct_stop_on_spanning_cluster )
{
ct_number_clusters = false;
ct_get_cluster_size = false;
}
ct_cluster_number = ct_num_spanning_clusters = ct_cluster_number_overflow = 0;
ct_spanning_cluster_found = false;
if ( ct_number_clusters )
clear_cluster_number_array();
mark_all_sites_unvisited_and_clear_path(); // S&B.C
ct_i = ct_j = ct_k = ct_l = -1;
// For 2d lattice ct_k and ct_l remain -1.
// For 3d lattice ct_l remains -1.
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( ct_stop_on_spanning_cluster )
if ( !IS_BOUNDARY_SITE_2(i,j) )
continue;
// A spanning cluster must include a boundary site
// so to find a spanning cluster it is sufficient
// to start traces only from boundary sites.
if ( !( ct_spin_value = spin_at_site(i,j) ) )
// Returns 0 if site is not occupied
// otherwise it returns the spin value.
continue; // because site not occupied.
if ( SITE_IS_VISITED_2(i,j) ) // ISM_DEF.C
continue;
// About to trace out a new cluster.
if ( ++ct_cluster_number == 0 )
{
ct_cluster_number_overflow++;
ct_cluster_number = 1; // Start again from 1.
}
if ( ct_get_cluster_size )
ct_size_of_cluster = 0;
// Now do cluster trace.
clear_stack(); // STACK.C
clear_path_status();
ct_i = i;
ct_j = j;
ct_spin_distance = 0;
trace_single_cluster(); // recursive function
// If counting cluster sizes (and not stopping on spanning cluster)
// then size of cluster is now in size_of_cluster.
// At this point stack_ptr should = stack and ct_spin_distance should = 0.
#if STACK_POINTER_CHECK
if ( !stack_empty() )
{
printf("\nStack not empty (CLUSTERS.C)!\n");
exit(1);
}
if ( ct_spin_distance )
{
printf("\nct_spin_distance not zero (CLUSTERS.C)!\n");
exit(1);
}
#endif
if ( ct_spanning_cluster_found && ct_stop_on_spanning_cluster )
goto done;
}
}
break;
case 3:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
ADD_CODE // 3d
}
}
}
break;
case 4:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
for ( l=0; l<size; l++ )
{
ADD_CODE // 4d
}
}
}
}
break;
}
done:
actd->num_clusters_traced = ct_cluster_number
+ (unsigned int)65535*ct_cluster_number_overflow;
actd->spanning_cluster_found = ct_spanning_cluster_found;
actd->num_spanning_clusters = ct_num_spanning_clusters;
actd->cluster_number_overflow = ct_cluster_number_overflow;
spanning_cluster_list[ct_num_spanning_clusters] = 0;
}
// This function is called (recursively) once for each site or spin.
// Only the return address is pushed onto the stack in an attempt
// to minimize the danger of system stack overflow.
// On entry the site is specified by ct_i, ct_j, ...
//
// Count number of bonds in cluster?
/*----------------------------------*/
static void trace_single_cluster(void)
{
if ( ct_get_cluster_size )
ct_size_of_cluster++;
switch ( dimensionality )
{
case 2:
MARK_SITE_AS_VISITED_2(ct_i,ct_j);
// We may wish to know if this is a spanning cluster
// even if we don't want to stop on a spanning cluster.
update_path_status();
if ( is_spanning_cluster(false) ) // false = only 1 dimension to be covered
{
ct_spanning_cluster_found = true;
if ( ct_number_clusters )
{
if ( ct_num_spanning_clusters == MAX_NUM_SPANNING_CLUSTERS )
{
printf("\nMore than %d spanning clusters found.\n",
MAX_NUM_SPANNING_CLUSTERS);
exit(1);
}
spanning_cluster_list[ct_num_spanning_clusters++] = ct_cluster_number;
}
if ( ct_stop_on_spanning_cluster )
{
MARK_SITE_IN_SPANNING_CLUSTER_2(ct_i,ct_j);
return;
}
}
if ( ct_number_clusters )
{
cluster_numbers2[ct_i][ct_j] = ct_cluster_number;
#if RECORD_SPIN_DISTANCES
spin_distances2[ct_i][ct_j] = ct_spin_distance;
#endif
}
// Now visit each of the unvisited nn spins.
loop
{
get_nn_sites(ct_i,ct_j); // Q&A.C
// nn_sites[][] now holds the indices
// of the nearest neighbour sites in all directions
// (no. of directions = coord_num).
// Not all of these are necessarily occupied.
best_r = NON_EXISTENT;
for ( ct_r=0; ct_r<coord_num; ct_r++ )
{
if ( ( nn_i = nn_sites[ct_r][0] ) == NON_EXISTENT )
continue;
// because this nn site does not exist (only with free boundary conditions)
nn_j = nn_sites[ct_r][1];
if ( !SITE_IS_OCCUPIED_2(nn_i,nn_j) )
continue; // because site is not occupied
if ( ct_cluster_type != CT_NN )
{
// ct_cluster_type is one of CT_OPEN_BOND or CT_SPIN.
if ( ! ( bonds2[ct_i][ct_j] & bit[ct_r] ) )
continue;
// because no open bond exists in direction ct_r
if ( ct_cluster_type == CT_SPIN )
{
// ct_spin_value is set in trace_all_clusters().
if ( ct_spin_value != spin_at_site(nn_i,nn_j) )
continue; // because spin value is not same
}
}
// Site exists and is occupied (or bond exists in direction ct_r).
if ( !SITE_IS_VISITED_2(nn_i,nn_j) )
{
// Found an unvisited nn site or spin.
if ( ct_stop_on_spanning_cluster )
{
// Find nn in best direction for spanning path.
if ( best_r == NON_EXISTENT )
best_r = ct_r;
else if ( directionality[ct_r] > directionality[best_r] )
best_r = ct_r;
}
else
{
best_r = ct_r;
break;
}
}
}
if ( best_r == NON_EXISTENT )
return; // because site has no unvisited nn spins.
// Have an unvisited occupied nn spin.
if ( !push_onto_stack(ct_i,ct_j) )
return; // otherwise would loop forever
ct_i = nn_sites[best_r][0];
ct_j = nn_sites[best_r][1];
ct_spin_distance++;
trace_single_cluster();
// A recursive call, which marks this site as visited.
if ( !pop_from_stack(&ct_i,&ct_j) )
{
printf("\nFatal error: Stack underflow in trace_single_cluster()!\n");
// should never occur
exit(1);
}
ct_spin_distance--;
if ( ct_spanning_cluster_found && ct_stop_on_spanning_cluster )
{
MARK_SITE_IN_SPANNING_CLUSTER_2(ct_i,ct_j);
// This usually does not mark all the sites in the spanning cluster,
// just those on the stack at the time that the spanning cluster is found.
// The marked sites usually, but not always, form a spanning path.
// In any case they show which cluster contains the spanning path.
return; // back up through the recursive calls.
}
}
break;
case 3:
ADD_CODE // 3d
break;
case 4:
ADD_CODE // 4d
break;
}
}
/*----------------------------------------*/
static void clear_cluster_number_array(void)
{
switch ( dimensionality )
{
case 2:
memset(&cluster_numbers2[0][0],0,num_sites);
break;
case 3:
memset(&cluster_numbers3[0][0][0],0,num_sites);
break;
case 4:
memset(&cluster_numbers4[0][0][0][0],0,num_sites);
break;
}
}
/*-------------------------------*/
static void clear_path_status(void)
{
int d, i;
for ( d=0; d<dimensionality; d++ )
for ( i=0; i<size; i++ )
edge[d][i] = false;
}
/*--------------------------------*/
static void update_path_status(void)
{
switch ( dimensionality )
{
case 2:
edge[0][ct_i] = edge[1][ct_j] = true;
break;
case 3:
edge[0][ct_i] = edge[1][ct_j] = edge[2][ct_k] = true;
break;
case 4:
edge[0][ct_i] = edge[1][ct_j] = edge[2][ct_k] = edge[3][ct_l]= true;
break;
}
}
/*------------------------------------------------------*/
static int is_spanning_cluster(int total_spanning_cluster)
{
int d, i;
if ( total_spanning_cluster )
{
for ( d=0; d<dimensionality; d++ )
{
for ( i=0; i<size; i++ )
{
if ( !edge[d][i] )
return ( false );
}
}
return ( true ); // because all dimensions are filled
}
else
{
for ( d=0; d<dimensionality; d++ )
{
for ( i=0; i<size; i++ )
{
if ( !edge[d][i] )
break;
}
if ( i== size )
return ( true ); // because a dimension is filled
}
return ( false );
}
}
// Assumes end of spanning cluster list marked by 0.
// Returns true if cluster_num is in the spanning cluster list, false otherwise.
/*--------------------------------------------------------*/
int in_spanning_cluster_list(unsigned short int cluster_num)
{
int i = 0;
loop
{
if ( spanning_cluster_list[i] == 0 )
return ( false );
if ( spanning_cluster_list[i++] == cluster_num )
return ( true );
}
}