src/Common/periodic.c

/*--------------------------------------------------------------------------*/
/* ALBERTA:  an Adaptive multi Level finite element toolbox using           */
/*           Bisectioning refinement and Error control by Residual          */
/*           Techniques for scientific Applications                         */
/*                                                                          */
/* www.alberta-fem.de                                                       */
/*                                                                          */
/*--------------------------------------------------------------------------*/
/*                                                                          */
/* file:     periodic.c                                                     */
/*                                                                          */
/* description:  support routines for periodic meshes: compute the orbit    */
/*               of macro element vertices under the action of the group    */
/*               of wall transformations. In 3d, we need also the orbits    */
/*               of the macro element edges.                                */
/*                                                                          */
/*--------------------------------------------------------------------------*/
/*                                                                          */
/* This file's authors: Claus-Justus Heine                                  */
/*                      Abteilung f�r Angewandte Mathematik                 */
/*                      Universit�t Freiburg                                */
/*                      Hermann-Herder-Stra�e 10                            */
/*                      79104 Freiburg, Germany                             */
/*                                                                          */
/* (c) by C.-J. Heine (2006).                                               */
/*                                                                          */
/*--------------------------------------------------------------------------*/

#include "alberta_intern.h"
#include "alberta.h"

/* Compute the orbit of a vertex (of the macro triangulatin) under the
 * action of the group of wall transformations.
 *
 * @param[in] dim The dimension of the mesh.
 *
 * @param[in] wall_trafos The wall transformations as permutation group acting
 *                  on the vertex indices of the macro triangulation.
 *
 * @param[in] nwt Number of Wall Transformations.
 *
 * @param[in] v The vertex to compute the orbit of.
 *
 * @param[out] orbit The orbit of e under the action of wall_trafos.
 *
 * @param[in] nv Number of Vertices.
 *
 * @return The length of the orbit.
 *
 */
int _AI_wall_trafo_vertex_orbit(int dim,
				int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2],
				int nwt,
				int v, int *orbit, int nv)
{
  char in_orbit[nv]; /* flag vectors, true if v is contained in the orbit */
  int i, j, k, l, new_v;
  int orb_len;

  /* initialisation: orbit is empty */
  for (i = 0; i < nv; i++) {
    in_orbit[i]   = false;
  }
  orb_len = 0;

  /* stuff v itself into the orbit */
  in_orbit[v]      = true;
  orbit[orb_len++] = v;

  for (i = 0; i < orb_len; i++) {
    v = orbit[i];

    for (j = 0; j < nwt; j++) {
      for (k = 0; k < N_VERTICES(dim-1); k++) {
	/* compute the image of v w.r.t. the wall trafo and its
	 * inverse. We would not nead the inverse, as the group of
	 * wall-trafos is finite, but this should speed up things.
	 */
	for (l = 0; l < 2; l++) {
	  if (v == wall_trafos[j][k][l]) {
	    new_v = wall_trafos[j][k][1-l];
	    if (!in_orbit[new_v]) {
	      orbit[orb_len++] = new_v;
	      in_orbit[new_v]  = true;
	    }
	    /* go to the next wall transformation, there can only be
	     * one match per transformation
	     */
	    goto next_wt;
	  }
	}
      }
    next_wt:
      continue;
    }
  }
  return orb_len;
}

/**Compute the number of orbits under the action of the group of wall
 * transformations on the vertices of the macro triangulation and map
 * each vertex to its orbit number. The map is returned in
 * "orbit_map". If orbit_map[e] == -1 then e is a fix-point.
 *
 * @param[in] wall_trafos The wall transformations as permutation group acting
 *                  on the global vertex indices of the macro triangulation.
 *
 * @param[n] nwt Number of Wall Transformations.
 *
 * @param[out] orbit_map The mapping of vertex numbers to orbit
 * numbers. Fix-points have the virtual orbit number -1. @a orbit_map
 * may be NULL.
 *
 * @param[in,out] nv On entry "*nv" holds the number of vertices, on
 * exit "*nv" holds the size of the set of orbits (i.e. the number of
 * vertices after identifying periodic walls).
 *
 * @return The number of non-trivial orbits, i.e. without counting
 * fix-points.
 *
 */
int _AI_wall_trafo_vertex_orbits(int dim,
				 int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2],
				 int nwt,
				 int *orbit_map, int *nvp)
{
  int v, i, nv = *nvp;
  int orb_len, orbit[nv];
  int n_orbits, nv_in_orbits;
  int orbit_map_space[nv];

  if (!orbit_map) {
    orbit_map = orbit_map_space;
  }
  for (v = 0; v < nv; v++) {
    orbit_map[v] = -1;
  }
  n_orbits = 0;
  nv_in_orbits = 0;
  *nvp = 0;
  for (v = 0; v < nv && nv_in_orbits < nv; v++) {
    if (orbit_map[v] < 0) {
      orb_len =
	_AI_wall_trafo_vertex_orbit(dim, wall_trafos, nwt, v, orbit, nv);
      nv_in_orbits += orb_len;
      ++(*nvp);
      if (orb_len > 1) {
	for (i = 0; i < orb_len; i++) {
	  orbit_map[orbit[i]] = n_orbits;
	}
	n_orbits++;
      }
    }
  }
  return n_orbits;
}

#if DIM_MAX == 3
/* We assign each unordered pair (i,j) with i < j a unique number
 * according to the lexicographical ordering (0,1), (0,2), ...,
 * (1,2),...
 *
 * The formula is just the difference of two arithmetic sums, taking
 * into account that we start the numbering at 0.
 *
 * (i,j) with i < j or j < i, nv is the number of vertices
 */
#define _PAIR_NR(i, j, nv) (((j)-(i)-1)+(i)*(2*(nv)-(i)-1)/2)
#define PAIR_NR(i, j, nv) ((i) < (j) ? _PAIR_NR(i, j, nv) : _PAIR_NR(j, i, nv))

/* Compute the orbit of an edge (of the macro triangulatin) under the
 * action of the group of wall transformations. Edges are represented
 * as (unordered) pairs of vertex numbers.
 *
 * @param[in] wall_trafos The wall transformations as permutation group acting
 *                  on the vertex indices of the macro triangulation.
 *
 * @param[in] nwt Number of Wall Transformations.
 *
 * @param[in] e The edge to compute the orbit of.
 *
 * @param[out] orbit The orbit of e under the action of wall_trafos.
 *
 * @param[in] edges The edges of the macro triangulation as unordered pairs of
 *                  of vertex numbers.
 *
 * @param[in] ne Number of Edges.
 *
 * @return The length of the orbit.
 *
 */
int _AI_wall_trafo_edge_orbit(int (*wall_trafos)[N_VERTICES_2D][2],
			      int nwt,
			      int e, int *orbit,
			      int (*edges)[2], int ne)
{
  const int edge_ind[N_EDGES_2D][2] = { { 1, 2 }, { 2, 0 }, {0, 1} };
  int i, j, nv;
  int edge_wall_trafos[nwt][N_EDGES_2D][2];

  nv = 0;
  for (i = 0; i < ne; i++) {
    nv = MAX(edges[i][0], MAX(edges[i][1], nv));
  }
  ++nv;

  {
    int np = nv*(nv-1)/2;
    int pair2edge[np];

    /* first build up the mappings between memory_3d.c edge numbering
     * and the numbering of edges interpreted as unordered pair
     */
    for (i = 0; i < np; i++) {
      pair2edge[i] = -1;
    }
    for (i = 0; i < ne; i++) {
      pair2edge[PAIR_NR(edges[i][0], edges[i][1], nv)] = i;
    }

    /* Now build a representation of the wall-transformations action
     * as permutations on the edge indices. Each wall transformation
     * is the product of three (N_EDGES_2D) disjoint transpositions.
     */
    for (i = 0; i < nwt; i++) {
      for (j = 0; j < N_EDGES_2D; j++) {
	int frompair, topair;

	frompair =  PAIR_NR(wall_trafos[i][edge_ind[j][0]][0],
			    wall_trafos[i][edge_ind[j][1]][0], nv);
	topair   =  PAIR_NR(wall_trafos[i][edge_ind[j][0]][1],
			    wall_trafos[i][edge_ind[j][1]][1], nv);
	DEBUG_TEST_EXIT(pair2edge[frompair] >= 0 && pair2edge[topair] >= 0,
			"Wall transformations do not seem to operate on "
			"the set of edges of the triangulation.\n");
	edge_wall_trafos[i][j][0] = pair2edge[frompair];
	edge_wall_trafos[i][j][1] = pair2edge[topair];
      }
    }
  }

  /******************************* BIG FAT NOTE ****************************/
  /*                                                                       */
  /* AT THIS POINT WE EXPLOIT THE FACT THAT N_EDGES_2D == N_VERTICES_2D,   */
  /* THIS MEANS WE CAN USE _AI_wall_trafo_vertex_orbit(): edge_wall_trafos */
  /* ACT ON THE SET OF UNSORTED TRIPLETS OF INTEGERS, THIS IS JUST WHAT    */
  /* wall_trafos DOES, FROM THE GROUP THEORITIC POINT OF VIEW THERE IS NO  */
  /* DIFFERENCE.                                                           */
  /*                                                                       */
  /*************************************************************************/

  return _AI_wall_trafo_vertex_orbit(DIM_MAX /* 3 */,
				     edge_wall_trafos, nwt,
				     e, orbit, ne);
}

/**Compute the number of orbits under the action of the group of wall
 * transformations and map each (macro triangulation) edge to its
 * proper orbit. The map is returned in "orbit_map". Negative entries
 * in "orbit_map" mean that the orbit has only one element, i.e. if
 * orbit_map[e] < 0 then e is a fix-point.
 *
 * @param[in] wall_trafos The wall transformations as permutation group acting
 *                  on the vertex indices of the macro triangulation.
 *
 * @param[n] nwt Number of Wall Transformations.
 *
 * @param[out] orbit_map The mapping of orbit numbers to edge
 *                  numbers. Fix-points have the virtual orbit number -1.
 *
 * @param[in] edges The edges of the macro triangulation as unordered pairs of
 *                  of vertex numbers.
 *
 * @param[in] ne Number of Edges.
 *
 * @return The number of non-trivial orbits, i.e. without counting
 * fix-points.
 *
 */
int _AI_wall_trafo_edge_orbits(int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2],
			       int nwt,
			       int *orbit_map,
			       int (*edges)[2], int ne)
{
  char in_orbit[ne]; /* flag vectors, true if e is contained in any orbit */
  int e, i;
  int orb_len, orbit[ne];
  int n_orbits, ne_in_orbits;

  for (e = 0; e < ne; e++) {
    in_orbit[e] = false;
  }
  n_orbits = 0;
  ne_in_orbits = 0;
  for (e = 0; e < ne && ne_in_orbits < ne; e++) {
    if (in_orbit[e] == false) {
      orb_len = _AI_wall_trafo_edge_orbit(wall_trafos, nwt,
					  e, orbit, edges, ne);
      ne_in_orbits += orb_len;
      if (orb_len > 1) {
	for (i = 0; i < orb_len; i++) {
	  in_orbit[orbit[i]] = true;
	  orbit_map[orbit[i]] = n_orbits;
	}
	n_orbits++;
      } else {
	orbit_map[e] = -1;
      }
    }
  }
  return n_orbits;
}
#endif

/* Given a mesh compute the generators for the group of wall
 * transformations for the macro-level. Vertices are numbered
 * according to the address of the coordinate pointers; we assume that
 * the coordinates are stored contiguously in memory (as is assumed in
 * the memory_Xd.c files).
 *
 * Return value is the number of wall-transformations.
 */
int
_AI_compute_macro_wall_trafos(MESH *mesh,
			      int (**wall_trafos_ptr)[N_VERTICES(DIM_MAX-1)][2])
{
  int dim          = mesh->dim;
  REAL_D   *coords = ((MESH_MEM_INFO *)mesh->mem_info)->coords;
  MACRO_EL *mel, *mel_n;
  int i, j, k;
  int el_wall_trafos[mesh->n_macro_el][N_WALLS(dim)];
  int (*wall_trafos)[N_VERTICES(DIM_MAX-1)][2] = NULL;
  int nwt = 0;

  if (!mesh->is_periodic) {
    *wall_trafos_ptr = NULL;
    return 0;
  }

  memset(el_wall_trafos, 0, sizeof(el_wall_trafos));
  for (i = 0; i < mesh->n_macro_el; i++) {
    mel = mesh->macro_els + i;
    for (j = 0; j < N_WALLS(dim); j++) {
      if (el_wall_trafos[i][j] != 0 || mel->neigh_vertices[j][0] == -1) {
	continue;
      }
      /* Need to create a new wall transformation */
      if (nwt % 100 == 0) {
	wall_trafos =(int (*)[N_VERTICES(DIM_MAX-1)][2])
	  MEM_REALLOC(wall_trafos,
		      nwt * sizeof(*wall_trafos),
		      (nwt+100) * sizeof(*wall_trafos), char);
      }
      mel_n = mel->neigh[j];
      for (k = 0; k < N_VERTICES(dim-1); k++) {
	wall_trafos[nwt][k][0] =
	  (int)((REAL_D *)mel->coord[(j+k+1) % N_VERTICES(dim)] - coords);
	wall_trafos[nwt][k][1] =
	  (int)((REAL_D *)mel_n->coord[mel->neigh_vertices[j][k]] - coords);
      }
      el_wall_trafos[mel->index][j] = nwt+1;
      el_wall_trafos[mel_n->index][mel->opp_vertex[j]] = -(nwt+1);
      nwt++;
    }
  }
  wall_trafos = (int (*)[N_VERTICES(DIM_MAX-1)][2])
    MEM_REALLOC(wall_trafos,
		(nwt + 100 - 1) / 100 * 100 * sizeof(*wall_trafos),
		nwt * sizeof(*wall_trafos),
		char);
  *wall_trafos_ptr = wall_trafos;
  return nwt;
}