Logo Search packages:      
Sourcecode: gamgi version File versions  Download package

gamgi_chem_polytope.c

/******************************************
 *
 * $GAMGI/src/chem/gamgi_chem_polytope.c
 *
 * Copyright (C) 2008 Carlos Pereira
 *
 * Distributed under the terms of the GNU
 * General Public License: $GAMGI/LICENSE
 *
 */

/*****************************************************
 *                   Bibliography                    *
 *                                                   *
 * Rourke J., Computational Geometry in C,           *
 * Cambridge University Press (1998)                 *
 *                                                   *
 * Berg M., Kreveld M., Overmars M., Schwarzkopf O., *
 * Computational Geometry, Springer Verlag (2000)    *
 *****************************************************/

#include "gamgi_engine.h"
#include "gamgi_mesa.h"

#include "gamgi_math_vector.h"
#include "gamgi_engine_create.h"
#include "gamgi_engine_start.h"
#include "gamgi_engine_dlist.h"
#include "gamgi_engine_list.h"

#define GAMGI_CHEM_POLYTOPE_WIDTH 0.1

#define GAMGI_CHEM_POLYTOPE_FACES_START  20
#define GAMGI_CHEM_POLYTOPE_FACES_STEP   10

#define GAMGI_CHEM_POLYTOPE_EDGES_START  20
#define GAMGI_CHEM_POLYTOPE_EDGES_STEP   10

#define GAMGI_CHEM_POLYTOPE_FACE_START   20
#define GAMGI_CHEM_POLYTOPE_FACE_STEP    10

#define GAMGI_CHEM_POLYTOPE_VERTEX_START 20
#define GAMGI_CHEM_POLYTOPE_VERTEX_STEP  10

#define GAMGI_CHEM_POLYTOPE_LOOPS_START  20
#define GAMGI_CHEM_POLYTOPE_LOOPS_STEP   10

#define GAMGI_CHEM_POLYTOPE_REMOVE 4
#define GAMGI_CHEM_POLYTOPE_MODIFY 3
#define GAMGI_CHEM_POLYTOPE_KEEP   2
#define GAMGI_CHEM_POLYTOPE_USE    1

#define GAMGI_CHEM_POLYTOPE_TOLERANCE_COPLANAR 1.0E-4
#define GAMGI_CHEM_POLYTOPE_TOLERANCE_VISIBLE  1.0E-4

typedef struct _gamgi_polytope {

gamgi_dlist *atoms;
int n_atoms;

gamgi_slist **cells;
int n_cells;

int central;
int neighbour;
int coordination;
double radius;
gamgi_bool periodic;
double offset;

gamgi_enum style;
float *color;

/***************************************************************
 *                          arrays format                      *
 *                                                             *
 * double *vertices_xyz: vertex_1_x, vertex_1_y, vertex_1_z,   *
 * ...,  vertex_n_x, vertex_n_y, vertex_n_z                    *
 *                                                             *
 * int **faces: face 1, ..., face n                            *
 *                                                             *
 * int **vertices: vertex 1, ..., vertex n                     *
 *                                                             *
 * int **edges: edges 1, ..., edges n                          *
 *                                                             *
 * int *vertex: total, state, number of edges, number          *
 * of conflicts, face 1, offset 1, ..., face n, offset n       *
 *                                                             *
 * int *face: total, state, number of conflicts, vertex 1,     *
 * offset 1, ..., vertex n, offset n, number of edges, edge 1, * 
 * before 1 next 1, ..., edge n, before n, next n              *
 *                                                             *
 * int *edge: total, state, face 1, face 2, vertex 1, vertex 2 *
 ***************************************************************/

double *xyz;

int **vertices;
int n_vertices;

int **edges;
int n_edges;

int **faces;
int n_faces;

} gamgi_polytope;

static int *static_export_points (gamgi_group *group, gamgi_polytope *polytope)
{
double *xyz, *points;
int **vertices;
int *vertex, *index;
int n_vertices;
int n, i;

xyz = polytope->xyz;
vertices = polytope->vertices;
n_vertices = polytope->n_vertices;

/******************************************
 * build the index array to redirect the  *
 * used vertices to their final positions *
 ******************************************/

n = 0;
index = (int *) malloc (n_vertices * sizeof (int));

for (i = 0; i < n_vertices; i++)
  {
  vertex = vertices[i];
  if (vertex[2] > 0) index[i] = n++;
  }

/******************************************************
 * build the points array that will contain the       *
 * positions of all the vertices of the polyhedron    *
 * contained in the group object linked to the parent *
 ******************************************************/

group->n_points = n;
points = (double *) malloc (3 * n * sizeof (double));
group->points = points;

for (i = 0; i < n_vertices; i++)
  {
  vertex = vertices[i];
  if (vertex[2] > 0)
    {
    n = index[i];
    points[3 * n + 0] = xyz[3 * i + 0];
    points[3 * n + 1] = xyz[3 * i + 1];
    points[3 * n + 2] = xyz[3 * i + 2];
    }
  }

return index;
}

static void static_export_loops (int *index, 
gamgi_group *group, gamgi_polytope *polytope)
{
int *loops;
int **faces, **edges;
int *face, *edge;
int f, v, e, e_start;
int offset, next;
int n_faces, n_conflicts, n_loops;
int i, j;

faces = polytope->faces;
n_faces = polytope->n_faces;
edges = polytope->edges;

group->n_loops = GAMGI_CHEM_POLYTOPE_LOOPS_START;
group->loops = (int *) malloc (group->n_loops * sizeof (int));

n_loops = 1;
offset = 1;
for (i = 0; i < n_faces; i++)
  {
  /********************
   * find valid faces *
   ********************/

  face = faces[i];
  if (face[1] == GAMGI_CHEM_POLYTOPE_USE) continue;
  loops[0]++;

  /*************************
   * find first valid edge *
   *************************/

  n_conflicts = face[2];
  for (j = 4 + 2 * n_conflicts; j < face[0] - 1; j += 3)
    if (face[j] < 0) continue;

  /**********************************
   * if needed increase loops array *
   **********************************/

  if (++n_loops > group->n_loops)
    {
    group->n_loops += GAMGI_CHEM_POLYTOPE_LOOPS_STEP;
    group->loops = (int *) realloc (group->loops, group->n_loops * sizeof (int));
    }

  loops[offset] = 0;
  e_start = face[j];
  next = face[j + 2];
  do
    {
    e = face[next];
    next = face[next + 2];

    edge = edges[e];
    if (edge[2] == f) v = edge[4];
    else v = edge[5];

    /****************************************************
     * redirect valid vertices to final array positions *
     ****************************************************/

    v = index[v];

    /**********************************
     * if needed increase loops array *
     **********************************/

    if (++n_loops > group->n_loops)
      {
      group->n_loops += GAMGI_CHEM_POLYTOPE_LOOPS_STEP;
      group->loops = (int *) realloc (group->loops, group->n_loops * sizeof (int));
      }

    loops[offset]++;
    loops[offset + loops[offset]] = v;
    } while (e != e_start);

  offset = loops[offset + loops[offset]] + 1;
  }

/**********************************
 * set final size for loops array *
 **********************************/

group->n_loops = n_loops;
group->loops = (int *) realloc (group->loops, group->n_loops * sizeof (int));
}

static void static_export (gamgi_atom *atom, gamgi_polytope *polytope)
{
gamgi_group *group;
int *index;

/**************************
 * create container group *
 **************************/

group = gamgi_engine_create_group ();
gamgi_engine_start_group (group);

/*******************************
 * export points, loops arrays *
 *******************************/

index = static_export_points (group, polytope);
static_export_loops (index, group, polytope);

/**************
 * set visual *
 **************/

if (polytope->style == GAMGI_MESA_SOLID) group->faces = TRUE;

if (polytope->color == NULL)
  {
  group->red = atom->red;
  group->green = atom->green;
  group->blue = atom->blue;
  }
else
  {
  group->red = polytope->color[0];
  group->green = polytope->color[1];
  group->blue = polytope->color[2];
  }

free (index);
}

static int static_before (int offset, int f, gamgi_polytope *polytope)
{
int *face, *edge;
int e;

face = polytope->faces[f];
offset = face[offset + 1];

e = face[offset];
edge = polytope->edges[e];

/***********************************************************
 * if the previous edge has face_hidden as first face      *
 * then get its first vertex, otherwise get the second     *
 * (this is the vertex farther away from the visible edge) *
 ***********************************************************/

if (edge[2] == f) return edge[4];

return edge[5];
}

static double static_mix (int v1, int v2, int v3,
int v_conflict, gamgi_polytope *polytope)
{
double *xyz;
double u[3], v[3], w[3];

xyz = polytope->xyz;

/***************************************************
 * v1, v2, v3 must be orientated counter-clockwise *
 * when seen from the outside of the polytope!     *
 *                                                 *
 *  by construction, u,v are non-colinear vectors  *
 ***************************************************/

gamgi_math_vector_sub (xyz + 3 * v2, xyz + 3 * v1, u);
gamgi_math_vector_sub (xyz + 3 * v3, xyz + 3 * v2, v);
gamgi_math_vector_sub (xyz + 3 * v_conflict, xyz + 3 * v2, w);

return gamgi_math_vector_mix (u, v, w);
}

static gamgi_bool static_visible (int v1, int v2, int v3, int v, 
gamgi_polytope *polytope)
{
double mix;

/***************************************************
 * v1, v2, v3 must be orientated counter-clockwise *
 * when seen from the outside of the polytope!     *
 ***************************************************/

mix = static_mix (v1, v2, v3, v, polytope);
if (mix > GAMGI_CHEM_POLYTOPE_TOLERANCE_VISIBLE) return TRUE;

return FALSE;
}

static gamgi_bool static_coplanar (int v_before, int v_start, 
int v_end, int v_new, gamgi_polytope *polytope)
{
double *xyz;
double u[3], v[3], w[3];
double dot;

xyz = polytope->xyz;

/*************************************************
 * by construction, u,v are non-colinear vectors *
 *************************************************/

gamgi_math_vector_sub (xyz + 3 * v_start, xyz + 3 * v_before, u);
gamgi_math_vector_sub (xyz + 3 * v_end, xyz + 3 * v_start, v);
gamgi_math_vector_cross (u, v, w);

gamgi_math_vector_sub (xyz + 3 * v_new, xyz + 3 * v_end, u);
dot = gamgi_math_vector_dot (w, u);

if (fabs (dot) < GAMGI_CHEM_POLYTOPE_TOLERANCE_COPLANAR) return TRUE;

return FALSE;
}

static void static_vertex_start (int v, gamgi_polytope *polytope)
{
int *vertex;
int i;

vertex = polytope->vertices[v];

vertex[1] = -1;
vertex[2] = 0;
vertex[3] = 0;
for (i = 4; i < vertex[0]; i++) vertex[i] = -1;
}

static int *static_vertex_create (int v, gamgi_polytope *polytope)
{
int *vertex;
int total;

/*****************************************
 * create vertex array, with 0 conflicts *
 *****************************************/

total = GAMGI_CHEM_POLYTOPE_VERTEX_START;
vertex = (int *) malloc (total * sizeof (int));
vertex[0] = total;

polytope->vertices[v] = vertex;
static_vertex_start (v, polytope);

return vertex;
}

static int static_vertex_conflict_recreate (int v, gamgi_polytope *polytope)
{
int **vertices;
int *vertex;
int i, j;

vertices = polytope->vertices;

/******************************************
 * look for a old face conflict position  *
 *                                        *
 * the first position for conflicts is 4  *
 ******************************************/

vertex = vertices[v];

for (i = 4; i < vertex[0] - 1; i+= 2)
  if (vertex[i] < 0) return i;

/**********************************************************
 * if that is not available increase and initialize array *
 **********************************************************/

vertex[0] += GAMGI_CHEM_POLYTOPE_VERTEX_STEP;
vertex =  (int *) realloc (vertex, vertex[0] * sizeof (int));
vertices[v] = vertex;

for (j = i; j < vertex[0]; j++) vertex[i] = -1;

return i;
}

static void static_edge_start (int e, gamgi_polytope *polytope)
{
int *edge;

/*****************************************************
 * this face is no longer in use and can be recycled *
 *****************************************************/

edge = polytope->edges[e];
edge[1] = GAMGI_CHEM_POLYTOPE_USE;
}

static void static_edge_remove (int e, gamgi_polytope *polytope)
{
int *edge, *vertex;
int v;

edge = polytope->edges[e];

/*********************************************************
 * for both vertices decrease counting (number of edges) *
 *********************************************************/

v = edge[4];
vertex = polytope->vertices[v];
vertex[2]--;

v = edge[5];
vertex = polytope->vertices[v];
vertex[2]--;

/*****************************************************
 * this face is no longer in use and can be recycled *
 *****************************************************/

static_edge_start (e, polytope);
}

static void static_edges_remove (int v_new, gamgi_polytope *polytope)
{
int **faces, **edges;
int *vertex_new, *face_visible, *edge_visible;
int f_visible, e_visible;
int i, j, n_conflicts;

faces = polytope->faces;
edges = polytope->edges;

/*******************************************************************
 * use the vertex conflict data to remove all doubly visible edges *
 *******************************************************************/

vertex_new = polytope->vertices[v_new];

for (i = 4; i < vertex_new[0] - 1; i += 2)
  {
  if (vertex_new[i] < 0) continue;

  f_visible = vertex_new[i];
  face_visible = faces[f_visible];

  n_conflicts = face_visible[2];
  for (j = 4 + 2 * n_conflicts; j < face_visible[0] - 2; j += 3)
    {
    e_visible = face_visible[4 + 2*n_conflicts + j];
    if (e_visible < 0) continue;

    edge_visible = edges[e_visible];

    if (edge_visible[1] == GAMGI_CHEM_POLYTOPE_REMOVE)
      static_edge_remove (e_visible, polytope);
    else
      edge_visible[1] = GAMGI_CHEM_POLYTOPE_KEEP;
    }
  }

}

static int *static_edge_create (int e, gamgi_polytope *polytope)
{
int *edge;
int total;

/**************************************************
 * create edge array, with 2 faces and 2 vertices *
 **************************************************/

total = 6;
edge = (int *) malloc (total * sizeof (int));
edge[0] = total;

polytope->edges[e] = edge;
static_edge_start (e, polytope);

return edge;
}

static int static_edges_create (gamgi_polytope *polytope)
{
int **edges;
int *edge;
int n_edges, i, j;

edges = polytope->edges;
n_edges = polytope->n_edges;

for (i = 0; i < n_edges; i++)
  {
  /*************************
   * reuse old edges array *
   *************************/

  edge = edges[i];
  if (edge[1] == GAMGI_CHEM_POLYTOPE_USE)
    {
    edge[1] = GAMGI_CHEM_POLYTOPE_KEEP;
    return i;
    }
  }

/************************
 * increase edges array *
 ************************/

n_edges += GAMGI_CHEM_POLYTOPE_EDGES_STEP;
polytope->n_edges = n_edges;

edges = (int **) realloc (edges, n_edges * sizeof (int *));
polytope->edges = edges;

for (j = i; j < n_edges; j++)
  {
  edge = static_edge_create (j, polytope);
  edges[j] = edge;
  }

edge = edges[i];
edge[1] = GAMGI_CHEM_POLYTOPE_KEEP;

return i;
}

static void static_edges_mark (int v_new, gamgi_polytope *polytope)
{
int **faces, **edges, **vertices;
int *face_start, *face_end, *face_visible;
int *vertex_new, *edge_visible;
int f_start, f_end, f_visible, e_visible;
int i, j;

faces = polytope->faces;
edges = polytope->edges;
vertices = polytope->vertices;

/*****************************************************************
 * use the vertex conflict data to scan all the visible faces,   *
 * to get all the edges that are susceptible of changes, without *
 * scanning the whole set of edges                               *
 *                                                               *
 * the edges in visible faces are either in the boundary of the  *
 * visible region (only one face is visible), forming new faces  *
 * with the new vertex, or are inside the visible region (both   *
 * faces are visible), and must be removed.                      *
 *****************************************************************/

vertex_new = vertices[v_new];

for (i = 4; i < vertex_new[0] - 1; i += 2)
  {
  if (vertex_new[i] < 0) continue;

  f_visible = vertex_new[i];
  face_visible = faces[f_visible];

  for (j = 4 + 2 * face_visible[2]; j < face_visible[0] - 2; j += 3)
    {
    if (face_visible[j] < 0) continue;

    e_visible = face_visible[j];
    edge_visible = edges[e_visible];
    edge_visible[1] = GAMGI_CHEM_POLYTOPE_MODIFY;

    f_start = edge_visible[2];
    face_start = faces[f_start];
    f_end = edge_visible[3];
    face_end = faces[f_end];

    if (face_start[1] == GAMGI_CHEM_POLYTOPE_REMOVE &&
    face_end[1] == GAMGI_CHEM_POLYTOPE_REMOVE)
      edge_visible[1] = GAMGI_CHEM_POLYTOPE_REMOVE;
    }
  }

}

static void static_face_start (int f, gamgi_polytope *polytope)
{
int *face;
int i;

face = polytope->faces[f];

/*****************************************************
 * this face is no longer in use and can be recycled *
 *****************************************************/

face[1] = GAMGI_CHEM_POLYTOPE_USE;
face[2] = 0;
for (i = 3; i < face[0]; i++) face[i] = -1;
}

static void static_face_remove (int f, gamgi_polytope *polytope)
{
int **vertices;
int *vertex, *face;
int v, i, offset;

vertices = polytope->vertices;

face = polytope->faces[f];

/******************************************
 * remove all conflicts for this face     *
 * from vertices before removing the face *
 ******************************************/

for (i = 3; i < 3 + 2 * face[2] - 1; i += 2)
  {
  v = face[i + 0];
  offset = face[i + 1];
  vertex = vertices[v];
  vertex[offset] = -1;
  }

/*****************************************************
 * this face is no longer in use and can be recycled *
 *****************************************************/

static_face_start (f, polytope);
}

static void static_faces_remove (int v_new, gamgi_polytope *polytope)
{
int *vertex_new;
int f_visible, i;

/************************************************************
 * use the vertex conflict data to remove all visible faces *
 ************************************************************/

vertex_new = polytope->vertices[v_new];

for (i = 4; i < vertex_new[0] - 1; i += 2)
  {
  if (vertex_new[i] < 0) continue;

  f_visible = vertex_new[i];
  static_face_remove (f_visible, polytope);
  }
}

static void static_face_clean (int f, gamgi_polytope *polytope)
{
int **edges;
int *face, *edge, *edge_next;
int before, next, next_next;
int e, e_next, i;

edges = polytope->edges;

/*****************************
 * use old edge if it exists *
 *****************************/

face = polytope->faces[f];

for (i = 3 + 2 * face[2] - 2; i < face[0] - 2; i += 3)
  {
  e = face[i];
  if (e < 0) continue;

  edge = edges[e];
  if (edge[1] == GAMGI_CHEM_POLYTOPE_KEEP && edge[2] == edge[3])
    {
    /*******************************************************
     * each edge has three positions, the first is the     *
     * edge number, the second is the offset of the edge   *
     * before and the third is the offset of the edge next *
     *                                                     *
     * the edges are orientated counter-clockwise in a     *
     * circular way: positions i+1 and i+2 indicate        *
     * the offset for the edge before and next             *
     *                                                     *
     * mark e and e_next to be removed later. This in      *
     * turn will remove the vertex between e and e_next    *
     *******************************************************/

    before = face[i + 1];
    next = face[i + 2];

    next_next = face[next + 2];
    face[before + 2] = next_next;
    face[next_next + 1] = before;

    edge[1] = GAMGI_CHEM_POLYTOPE_REMOVE;
    e_next = face[next];
    edge_next = edges[e_next];
    edge_next[1] = GAMGI_CHEM_POLYTOPE_REMOVE;
    }
  }

face[1] = GAMGI_CHEM_POLYTOPE_KEEP;
}

static int static_face_edge_recreate (int f, gamgi_polytope *polytope)
{
int **faces, **edges;
int *face, *edge;
int e, i, j;

faces = polytope->faces;
edges = polytope->edges;

face = faces[f];

/*****************************
 * use old edge if it exists *
 *****************************/

for (i = 4 + 2 * face[2]; i < face[0] - 2; i += 3)
  {
  e = face[i];
  if (e < 0) return i;

  edge = edges[e];
  if (edge[1] == GAMGI_CHEM_POLYTOPE_USE) return i;
  }

/**********************************************************
 * if that is not available increase and initialize array *
 **********************************************************/

face[0] += GAMGI_CHEM_POLYTOPE_FACE_STEP;
face =  (int *) realloc (face, face[0] * sizeof (int));
faces[f] = face;

for (j = i; j < face[0]; j++) face[j] = -1;

return i;
}

static int static_face_conflict_recreate (int f, gamgi_polytope *polytope)
{
int **faces;
int *face;
int i, j;

faces = polytope->faces;
face = faces[f];

/******************************************
 * look for a old face conflict position  *
 *                                        *
 * the first position for conflicts is 3  *
 *                                        *
 * the last 10 positions are reserved for *
 * n_edges and 3 edges, for this new face *
 ******************************************/

for (i = 3; i < face[0] - 11; i+= 2)
  if (face[i] < 0) return i;

/*******************************************
 * if that is not available increase array *
 *******************************************/

face[0] += GAMGI_CHEM_POLYTOPE_FACE_STEP;
face = (int *) realloc (face, face[0] * sizeof (int));
faces[f] = face;

for (j = i; j < face[0]; j++) face[j] = -1;

return i;
}

static void static_conflict_add (int v1, int v2, int v3, 
int v, int f, int *n_conflicts, gamgi_polytope *polytope)
{
int *face, *vertex;
int offset_f, offset_v;

/***************************************************
 * v1, v2, v3 must be orientated counter-clockwise *
 * when seen from the outside of the polytope!     *
 ***************************************************/

if (static_visible (v1, v2, v3, v, polytope) == TRUE)
  {
  offset_v = static_vertex_conflict_recreate (v, polytope);
  offset_f = static_face_conflict_recreate (f, polytope);

  face = polytope->faces[f];
  face[offset_f + 0] = v;
  face[offset_f + 1] = offset_v;

  vertex = polytope->vertices[v];
  vertex[offset_v + 0] = f;
  vertex[offset_v + 1] = offset_f;

  *n_conflicts += 1;
  }
}

static int static_conflict_face (int f_old, int f, 
int v_start, int v_end, int v_new, gamgi_polytope *polytope)
{
int **faces, **vertices;
int *face, *face_old;
int n_conflicts;
int i, v;

faces = polytope->faces;
vertices = polytope->vertices;

face_old = faces[f_old];
face = faces[f];

n_conflicts = 0;
for (i = 3; i < 3 + 2 * face_old[2] - 1; i += 2)
  {
  v = face_old[i];

  /***************************************************************
   * vertices supplied to static_visible must be orientated      *
   * counter-clockwise when seen from the outside. As v_start,   *
   * v_end are orientated for the first face (face hidden), they *
   * must be reversed for the second face (face new). The third  *
   * vertex (the new vertex) of this triangular face comes next. *
   ***************************************************************/

  static_conflict_add (v_end, v_start, v_new, v, f, &n_conflicts, polytope);
  }

return n_conflicts;
}

static int static_conflict_rebuild (int f, int e, int v, gamgi_polytope *polytope)
{
int *edge;
int f_start, f_end;
int v_start, v_end;
int n_conflicts;

edge = polytope->edges[e];

/************************************************
 * vertex conflicts to new face must come from  *
 * vertex conflicts with faces adjacent to edge *
 ************************************************/

n_conflicts = 0;
f_start = edge[2];
f_end = edge[3];

v_start = edge[4];
v_end = edge[5];

n_conflicts += static_conflict_face (f, f_start, v_start, v_end, v, polytope);
n_conflicts += static_conflict_face (f, f_end, v_start, v_end, v, polytope);

return n_conflicts;
}

static void static_conflict_build (gamgi_polytope *polytope)
{
int **faces, **vertices;
int *face;
int n_conflicts, n_faces, n_vertices;
int v1, v2, v3;
int f, v;

faces = polytope->faces;
n_faces = 4;

vertices = polytope->vertices;
n_vertices = polytope->n_vertices;

for (f = 0; f < n_faces; f++)
  {
  face = faces[f];

  /* add vertex function here */

  n_conflicts = 0;
  for (v = 4; v < n_vertices; v++)
    {
    /***************************************************************
     * vertices supplied to static_visible must be orientated      *
     * counter-clockwise when seen from the outside. As v_start,   *
     * v_end are orientated for the first face (face hidden), they *
     * must be reversed for the second face (face new). The third  *
     * vertex (the new vertex) of this triangular face comes next. *
     ***************************************************************/

    static_conflict_add (v1, v2, v3, v, f, &n_conflicts, polytope);
    }

  face[2] = n_conflicts;
  }

}

static int *static_face_create (int f, gamgi_polytope *polytope)
{
int *face;
int total;

/***************************************************
 * create face array, with 3 edges and 0 conflicts *
 ***************************************************/

total = GAMGI_CHEM_POLYTOPE_FACE_START;
face = (int *) malloc (total * sizeof (int));
face[0] = total;

polytope->faces[f] = face;
static_face_start (f, polytope);

return face;
}

static int static_faces_create (gamgi_polytope *polytope)
{
int **faces;
int *face;
int n_faces;
int i, j;

faces = polytope->faces;
n_faces = polytope->n_faces;

for (i = 0; i < n_faces; i++)
  {
  /*************************
   * reuse old faces array *
   *************************/

  face = faces[i];
  if (face[1] == GAMGI_CHEM_POLYTOPE_USE)
    {
    face[1] = GAMGI_CHEM_POLYTOPE_KEEP;
    return i;
    }
  }

/************************
 * increase faces array *
 ************************/

n_faces += GAMGI_CHEM_POLYTOPE_FACES_STEP;
polytope->n_faces = n_faces;

faces = (int **) realloc (faces, n_faces * sizeof (int *));
polytope->faces = faces;

for (j = i; j < n_faces; j++)
  {
  face = static_face_create (j, polytope);
  faces[j] = face;
  }

face = faces[i];
face[1] = GAMGI_CHEM_POLYTOPE_KEEP;

return i;
}

static void static_faces_mark (int v_new, gamgi_polytope *polytope)
{
int **faces, **vertices;
int *vertex_new, *face_visible;
int f_visible;
int i;

faces = polytope->faces;
vertices = polytope->vertices;

/*******************************************
 * signal all faces that should be removed *
 *******************************************/

vertex_new = vertices[v_new];

for (i = 4; i < vertex_new[0] - 1; i += 2)
  {
  f_visible = vertex_new[i];
  if (f_visible < 0) continue;

  face_visible = faces[f_visible];
  face_visible[2] = GAMGI_CHEM_POLYTOPE_REMOVE;
  }

}

static void static_face_build (int v_new, int v_start, int v_end,
int e_visible, gamgi_polytope *polytope)
{
int **faces, **edges, **vertices;
int *face_new, *vertex_new;
int *vertex_start, *vertex_end;
int *edge_start, *edge_end;
int f_new, e_start, e_end;
int n_conflicts, offset;

faces = polytope->faces;
edges = polytope->edges;
vertices = polytope->vertices;

vertex_new = vertices[v_new];
vertex_start = vertices[v_start];
vertex_end = vertices[v_end];

/******************************************************
 * Create new face, with 3 vertices and 3 edges,      *
 * counter-clockwise: for each edge, the vertices are *
 * orientated counter-clockwise for the first face,   *
 * for the second face the vertices must be reversed  *
 ******************************************************/
 
f_new = static_faces_create (polytope);
n_conflicts = static_conflict_rebuild (f_new, e_visible, v_new, polytope);

face_new = faces[f_new];
face_new[2] = n_conflicts;
offset = 3 + 2 * n_conflicts;

face_new[offset + 0] = 3;
face_new[offset + 1] = e_visible;
face_new[offset + 2] = offset + 7;
face_new[offset + 3] = offset + 4;

if (vertex_start[1] < 0)
  {
  /******************************************************
   * create new edge, pointing to edges before and next *
   ******************************************************/

  e_start = static_edges_create (polytope);
  edge_start = edges[e_start];
  face_new[offset + 4] = e_start;
  face_new[offset + 5] = offset + 1;
  face_new[offset + 6] = offset + 7;
  edge_start[2] = f_new;
  edge_start[3] = -1;
  edge_start[4] = v_start;
  edge_start[5] = v_new;

  /**************************************************
   * increase vertex counting (its number of edges) *
   **************************************************/

  vertex_start[2]++;
  vertex_new[2]++;

  /****************************
   * point vertex to new edge *
   ****************************/

  vertex_start[1] = e_start;
  }
else
  {
  /***************************
   * the edge already exists *
   ***************************/

  e_start = vertex_start[1];
  edge_start = edges[e_start];
  face_new[offset + 4] = e_start;
  face_new[offset + 5] = offset + 1;
  face_new[offset + 6] = offset + 7;
  edge_start[3] = f_new;

  /******************************
   * reset vertex edge pointing *
   ******************************/

  vertex_start[1] = -1;
  }

if (vertex_end[1] < 0)
  {
  /******************************************************
   * create new edge, pointing to edges before and next *
   ******************************************************/

  e_end = static_edges_create (polytope);
  edge_end = edges[e_end];
  face_new[offset + 7] = e_end;
  face_new[offset + 8] = offset + 4;
  face_new[offset + 9] = offset + 1;
  edge_end[2] = f_new;
  edge_end[3] = -1;
  edge_end[4] = v_new;
  edge_end[5] = v_end;

  /**************************************************
   * increase vertex counting (its number of edges) *
   **************************************************/

  vertex_new[2]++;
  vertex_end[2]++;

  /****************************
   * point vertex to new edge *
   ****************************/

  vertex_end[1] = e_end;
  }
else
  {
  /***************************
   * the edge already exists *
   ***************************/

  e_end = vertex_end[1];
  edge_end = edges[e_end];
  face_new[offset + 7] = e_end;
  face_new[offset + 8] = offset + 4;
  face_new[offset + 9] = offset + 1;
  edge_end[4] = f_new;

  /******************************
   * reset vertex edge pointing *
   ******************************/

  vertex_end[1] = -1;
  }

}

static void static_face_merge (int v_new, int v_start, int v_end,
int e_visible, int before, int next, gamgi_polytope *polytope)
{
int **faces, **edges, **vertices;
int *face_hidden, *vertex_new;
int *vertex_start, *vertex_end;
int *edge_start, *edge_end, *edge_visible;
int f_hidden;
int e_start, e_end;
int start, end;

faces = polytope->faces;
edges = polytope->edges;
vertices = polytope->vertices;

vertex_new = vertices[v_new];
vertex_start = vertices[v_start];
vertex_end = vertices[v_end];

edge_visible = edges[e_visible];
f_hidden = edge_visible[2];
face_hidden = faces[f_hidden];

/***********************************************************
 * This is the same face as face_hidden: 1) add vertex_new *
 * to face_hidden, in counter-clockwise direction.         *
 *                                                         *
 * 2) mark face_hidden to be checked again in the end,     *
 * looking for common edges. This happens when v_new is    *
 * coplanar with more than one face, in which case         *
 * redundant edges and vertices must be removed.           *
 ***********************************************************/

face_hidden[2] = GAMGI_CHEM_POLYTOPE_MODIFY;

if (vertex_start[1] < 0)
  {
  /*******************
   * create new edge *
   *******************/

  e_start = static_edges_create (polytope);
  edge_start = edges[e_start];

  start = static_face_edge_recreate (f_hidden, polytope);
  face_hidden[start] = e_start;

  edge_start[2] = f_hidden;
  edge_start[3] = -1;
  edge_start[4] = v_start;
  edge_start[5] = v_new;

  /**************************************************
   * increase vertex counting (its number of edges) *
   **************************************************/

  vertex_start[2]++;
  vertex_new[2]++;

  /****************************
   * point vertex to new edge *
   ****************************/

  vertex_start[1] = e_start;
  }
else
  {
  /***************************
   * the edge already exists *
   ***************************/

  e_start = vertex_start[1];
  edge_start = edges[e_start];
  edge_start[3] = f_hidden;

  start = static_face_edge_recreate (f_hidden, polytope);
  face_hidden[start] = e_start;

  /******************************
   * reset vertex edge pointing *
   ******************************/

  vertex_start[1] = -1;
  }

if (vertex_end[1] < 0)
  {
  /*******************
   * create new edge *
   *******************/

  e_end = static_edges_create (polytope);
  edge_end = edges[e_end];

  end = static_face_edge_recreate (f_hidden, polytope);
  face_hidden[end] = e_end;

  edge_end[2] =  f_hidden;
  edge_end[3] =  -1;
  edge_end[4] = v_new;
  edge_end[5] = v_end;

  /**************************************************
   * increase vertex counting (its number of edges) *
   **************************************************/

  vertex_new[2]++;
  vertex_end[2]++;

  /****************************
   * point vertex to new edge *
   ****************************/

  vertex_end[1] = e_end;
  }
else
  {
  /***************************
   * the edge already exists *
   ***************************/

  e_end = vertex_end[1];
  edge_end = edges[e_end];
  edge_end[3] = f_hidden;

  end = static_face_edge_recreate (f_hidden, polytope);
  face_hidden[end] = e_end;

  /******************************
   * reset vertex edge pointing *
   ******************************/

  vertex_end[1] = -1;
  }

/*************************************************
 * update before and next pointers to all edges, *
 * ordered as before -> start -> end -> next,    *
 * to keep all edges orientated counter-clocwise *
 *************************************************/

face_hidden[start + 1] = before;
face_hidden[start + 2] = e_end;
face_hidden[end + 1] = e_start;
face_hidden[end + 2] = next;

face_hidden[before + 2] = start;
face_hidden[next + 1] = end;
 
edge_visible[1] = GAMGI_CHEM_POLYTOPE_REMOVE;
}

static void static_build_vertex (int v_new, gamgi_polytope *polytope)
{
int **faces, **edges, **vertices;
int *vertex_new, *edge_visible, *face_visible, *face_edge;
int v_before, v_start, v_end, e_visible, f_visible, f_edge;
int offset, before, next;
int i, j, aux;

faces = polytope->faces;
edges = polytope->edges;
vertices = polytope->vertices;

static_faces_mark (v_new, polytope);
static_edges_mark (v_new, polytope);

/*****************************************************************
 * use again the vertex conflict data to scan all the            *
 * visible faces, to get all the edges that are susceptible      *
 * of changes, without scanning the whole set of edges           *
 *                                                               *
 * build a new face for each edge in the boundary of the visible *
 * region (only one face is visible), linking its two vertices   *
 * with the new vertex,                                          *
 *****************************************************************/

for (i = 4; i < vertex_new[0] - 1; i += 2)
  {
  if (vertex_new[i] < 0) continue;

  f_visible = vertex_new[i];
  face_visible = faces[f_visible];

  for (j = 4 + 2 * face_visible[2]; j < face_visible[0] - 2; j += 3)
    {
    if (face_visible[j] < 0) continue;

    e_visible = face_visible[j];
    edge_visible = edges[e_visible];

    before = face_visible[j + 1];
    next = face_visible[j + 2];

    /*************************************************
     * possible cases: 1) both faces are visible,    *
     * and edge will be removed only in the end;     *
     * 2) new coplanar edges have been added to face *
     * and should be ignored; 3) only one face is    *
     * visible and the second face must be created   *
     *************************************************/

    if (edge_visible[1] == GAMGI_CHEM_POLYTOPE_REMOVE) continue;

    if (edge_visible[1] == GAMGI_CHEM_POLYTOPE_KEEP) continue;

    /****************************************************
     * edge format: face1 face2 vertex1 vertex2 where   *
     * vertex1 vertex2 are orientated counter-clockwise *
     * for face1 (and clockwise for face2)              *
     ****************************************************/

    f_edge = edge_visible[2];
    face_edge = faces[f_edge];
    if (face_edge[1] == GAMGI_CHEM_POLYTOPE_REMOVE)
      {
      /************************************************************ 
       * The face that survives, face_hidden, should be first     *
       * and vertices should be ordered according to face_hidden. *
       * When face_hidden is not first, swap faces and vertices.  *
       ************************************************************/

      aux = edge_visible[2];
      edge_visible[2] = edge_visible[3];
      edge_visible[3] = aux;

      aux = edge_visible[4];
      edge_visible[4] = edge_visible[5];
      edge_visible[5] = aux;
      }

    /************************************************
     * get vertices in visible edge and edge before *
     ************************************************/

    v_start = edge_visible[4];
    v_end = edge_visible[5];
    v_before = static_before (offset, edge_visible[2], polytope);

    if (static_coplanar (v_before, v_start, v_end, v_new, polytope) == FALSE)
      static_face_build (v_new, v_start, v_end, e_visible, polytope);
    else
      static_face_merge (v_new, v_start, v_end, e_visible, before, next, polytope);
    }
  }

/***********************************************************
 * use the vertex conflict data to scan all visible faces, *
 * looking for edges, vertices to remove in the same face  *
 ***********************************************************/

for (i = 4; i < vertex_new[0] - 1; i += 2)
  {
  f_visible = vertex_new[i];
  face_visible = faces[f_visible];
  if (face_visible[2] == GAMGI_CHEM_POLYTOPE_MODIFY)
    static_face_clean (f_visible, polytope);
  }

static_edges_remove (v_new, polytope);
static_faces_remove (v_new, polytope);
}

static void static_build_reset (gamgi_polytope *polytope)
{
int n_faces, n_edges, n_vertices;
int i;

/****************
 * reset arrays *
 ****************/

n_edges = polytope->n_edges;
for (i = 0; i < n_edges; i++)
  static_edge_start (i, polytope);

n_vertices = polytope->n_vertices;
for (i = 0; i < n_vertices; i++)
  static_vertex_start (i, polytope);

n_faces = polytope->n_faces;
for (i = 0; i < n_faces; i++)
  static_face_start (i, polytope);
}

static void static_build_vertices (int cell_x, int cell_y, int cell_z,
gamgi_slist *start, gamgi_atom *atom, gamgi_polytope *polytope)
{
/* go around cells and create vertices array, with closer
vertices, in reverse order, so the vertices more far away
are handled first */
}

static void static_build_tetrahedron (gamgi_polytope *polytope)
{
/* build first polygon with conflict lists */

}

static void static_build_central (int cell_x, int cell_y, int cell_z, 
gamgi_slist *start, gamgi_atom *atom, gamgi_polytope *polytope)
{
int i;

static_build_reset (polytope);
static_build_vertices (cell_x, cell_y, cell_z, start, atom, polytope);
static_build_tetrahedron (polytope);
static_conflict_build (polytope);

for (i = 4; i < polytope->n_vertices; i++)
  static_build_vertex (i, polytope);

static_export (atom, polytope);
}

static gamgi_bool static_build (gamgi_polytope *polytope)
{
gamgi_atom *atom;
gamgi_slist **cells;
gamgi_slist *slist, *start;
int cell_x, cell_y, cell_z;
int n_cells, offset;

cells = polytope->cells;
n_cells = polytope->n_cells;

for (cell_z = 0; cell_z < n_cells; cell_z++)
  {
  for (cell_y = 0; cell_y < n_cells; cell_y++)
    {
    for (cell_x = 0; cell_x < n_cells; cell_x++)
      {
      offset = cell_z * n_cells * n_cells + cell_y * n_cells + cell_x;
      start = cells[offset];
      for (slist = start; slist != NULL; slist = slist->next)
        {
        /*****************************************************
         * build polyhedron from distance-ordered neighbours *
         *****************************************************/

        atom = GAMGI_CAST_ATOM slist->data;
        if (atom->element == polytope->central)
          static_build_central (cell_x, cell_y, cell_z, start, atom, polytope);
        }

      }
    }
  }

return TRUE;
}

static void static_start_cells (gamgi_polytope *polytope)
{
gamgi_atom *atom;
gamgi_slist **cells;
gamgi_dlist *dlist;
double min_x, max_x;
double min_y, max_y;
double min_z, max_z; 
double width_x, width_y, width_z;
double r_cells, min;
int cell_x, cell_y, cell_z;
int central, neighbour;
int offset;
int n;

central = polytope->central;
neighbour = polytope->neighbour;

/**************************************************
 * determine: 1) system dimensions; 2) number of  *
 * atoms; 3) maximum radius (radical tesselation) *
 **************************************************/

min_x = DBL_MAX; max_x = -DBL_MAX;
min_y = DBL_MAX; max_y = -DBL_MAX;
min_z = DBL_MAX; max_z = -DBL_MAX;

n = 0;
for (dlist = polytope->atoms; dlist != NULL; dlist = dlist->next)
  {
  /*****************************************************
   * give a different, sequential, number to each atom *
   *                                                   *
   *     find minimum and maximum cell boundaries      *
   *                                                   *
   *     find maximum radius (radical tesselation)     *
   *****************************************************/

  atom = GAMGI_CAST_ATOM dlist->data;
  if (atom->element != central && atom->element != neighbour) continue;
  n++;

  if (atom->position[0] < min_x) min_x = atom->position[0];
  if (atom->position[0] > max_x) max_x = atom->position[0];

  if (atom->position[1] < min_y) min_y = atom->position[1];
  if (atom->position[1] > max_y) max_y = atom->position[1];

  if (atom->position[2] < min_z) min_z = atom->position[2];
  if (atom->position[2] > max_z) max_z = atom->position[2];
  }
polytope->n_atoms = n;

/************************************************************
 * determine number of cells on each direction:             *
 * on average, each cell has approximately one atom.        *
 *                                                          *
 * allocate array of cells:                                 *
 * each cell points to a linked list of atoms in that cell. *
 ************************************************************/

if (modf (pow (n, 1/3.0), &r_cells) > 0.5) r_cells++;
polytope->n_cells = n = r_cells;

cells = (gamgi_slist **) calloc (pow (n, 3), sizeof (gamgi_slist *));
polytope->cells = cells;

/***********************************************************
 * determine cell dimensions: when the width is too small  *
 * (for example when the atoms have the same z coordinate) *
 * we need to enlarge the cell dimensions symmetricaly, to *
 * guarantee that the distance from the seeds at old min,  *
 * max coordinates to the new border planes is the same    *
 ***********************************************************/

min = GAMGI_CHEM_POLYTOPE_WIDTH;

if (max_x - min_x < min)
  {
  min_x = (min_x + max_x - min) / 2.0;
  max_x = (min_x + max_x + min) / 2.0;
  }
width_x = (max_x - min_x) / n;

if (max_y - min_y < min)
  {
  min_y = (min_y + max_y - min) / 2.0;
  max_y = (min_y + max_y + min) / 2.0;
  }
width_y = (max_y - min_y) / n;

if (max_z - min_z < min)
  {
  min_z = (min_z + max_z - min) / 2.0;
  max_z = (min_z + max_z + min) / 2.0;
  }
width_z = (max_z - min_z) / n;

/********************************************
 * 1) find cell where each atom belongs     *
 * 2) add atom to linked list for that cell *
 ********************************************/

for (dlist = polytope->atoms; dlist != NULL; dlist = dlist->next)
  {
  atom = GAMGI_CAST_ATOM dlist->data; 
  if (atom->element != central && atom->element != neighbour) continue;

  cell_x = floor ((atom->position[0] - min_x) / width_x);
  if (cell_x == n) cell_x--;
  cell_y = floor ((atom->position[1] - min_y) / width_y);
  if (cell_y == n) cell_y--;
  cell_z = floor ((atom->position[2] - min_z) / width_z);
  if (cell_z == n) cell_z--;

  offset = cell_z * n * n + cell_y * n + cell_x;
  cells[offset] = gamgi_engine_slist_add_start (cells[offset]);
  cells[offset]->data = atom;
  }

}

static void static_end_cells (gamgi_polytope *polytope)
{
gamgi_slist **cells;
gamgi_slist *slist;
int n_cells;
int cell_x, cell_y, cell_z;
int offset;

/********************************************
 * remove lists of nodes, one for each cell *
 ********************************************/

cells = polytope->cells;
n_cells = polytope->n_cells;
for (cell_z = 0; cell_z < n_cells; cell_z++)
  {
  for (cell_y = 0; cell_y < n_cells; cell_y++)
    {
    for (cell_x = 0; cell_x < n_cells; cell_x++)
      {
      offset = cell_z * n_cells * n_cells + cell_y * n_cells + cell_x;
      slist = cells[offset];

      while (slist != NULL)
        slist = gamgi_engine_slist_remove_start (slist);
      }
    }
  }

}

static gamgi_polytope *static_start (gamgi_object *parent, int central, 
int neighbour, int coordination, double radius, gamgi_bool periodic, 
double offset, gamgi_enum style, float *color)
{
gamgi_polytope *polytope;
gamgi_dlist *atoms;

/*****************************************
 * allocate atom list and main structure *
 *****************************************/

atoms = gamgi_engine_dlist_atom_object (parent, NULL);
if (atoms == NULL) return NULL;

polytope = (gamgi_polytope *) malloc (sizeof (gamgi_polytope));
polytope->atoms = atoms;

/***************************
 * allocate atoms in cells *
 ***************************/

static_start_cells (polytope);

/*******************
 * save input data *
 *******************/

polytope->central = central;
polytope->neighbour = neighbour;
polytope->coordination = coordination;
polytope->radius = radius;
polytope->periodic = periodic;
polytope->offset = offset;
polytope->style = style;
polytope->color = color;

return polytope;
}

static void static_end (gamgi_polytope *polytope)
{
gamgi_dlist *dlist;
int i;

/***************
 * free arrays *
 ***************/

for (i = 0; i < polytope->n_faces; i++) free (polytope->faces[i]);
free (polytope->faces);

for (i = 0; i < polytope->n_edges; i++) free (polytope->edges[i]);
free (polytope->edges);

for (i = 0; i < polytope->n_vertices; i++) free (polytope->vertices[i]);
free (polytope->vertices);

free (polytope->xyz);

/**************
 * free cells *
 **************/

static_end_cells (polytope);

/**************************
 * reset atoms, free list *
 **************************/

dlist = polytope->atoms;
while (dlist != NULL)
  {
  (GAMGI_CAST_ATOM dlist->data)->mark = 0;
  dlist = gamgi_engine_dlist_remove_start (dlist);
  }

free (polytope);
}

gamgi_bool gamgi_chem_polytope (gamgi_object *parent, int central, 
int neighbour, int coordination, double radius, gamgi_bool periodic, 
double offset, gamgi_enum style, float *color)
{
gamgi_polytope *polytope;
gamgi_bool valid;

polytope = static_start (parent, central, neighbour, 
coordination, radius, periodic, offset, style, color);
if (polytope == NULL) return FALSE;

valid = static_build (polytope);

static_end (polytope);

return valid;
}

Generated by  Doxygen 1.6.0   Back to index