Examples for the implementation of element types

 

In this section we give some examples for the ELEMENT3D_DESCRIPTION. We point out one possible way how to define parts of the ELEMENT3D_DESCRIPTION for a cube, a prism, and a tetrahedron. For all elements we define local_coord, polygon_length, polygon_vertex, and polygon_neighbour. In the case of a cuboidal grid we define the main part of the function neighbour(). The second function boundary() uses the same ideas as neighbour() but is much easier to implement.

As a first example we start with one of the simplest mesh types and consider a set of equidistributed points given by a Finite Difference method. This should underline that beside general meshes containing various types of elements also the most basic grids can easily be handled by this concept. The whole grid is totally described by the constant mesh sizes h_x, h_y, h_z and the numbers of mesh points n_x, n_y, n_z in the x, y, and z-direction. Such information can be handled by the following

typedef
{
  int    n_x, n_y, n_z;     /* number of points in x, y, z-direction      */
  float  h_x, h_y, h_z;     /* distance of points in x, y, z-direction    */
  float  **f_data;          /* pointer to data matrix                     */
} USER_DATA;

This set of points describes a cuboidal grid. All cubes of the grid are of equal shape with diameters h_x, h_y, h_z. A single cube is totally described by the index triple of one vertex, for example the vertex with smallest x, y and z coordinates.

typedef
{
  int i_x, i_y, i_z;        /* user's reference for the actual cube       */
} CUBE;

For a cube with index (i_x, i_y, i_z) the world coordinates of the reference vertex are simple (i_x*h_x, i_y*h_y, i_z*h_z). The world coordinates of the other vertices are calculated in a similar way. This is done by a procedure fill_points(ELEMENT3D *el).

We use the following components of the ELEMENT3D_DESCRIPTION for one single cube :

/* vertex indices of the polygons                                         */
static int c_v0[4] = {0,3,2,1}, c_v1[4] = {4,5,6,7}, c_v2[4] = {0,1,5,4};
static int c_v3[4] = {1,2,6,5}, c_v4[4] = {2,3,7,6}, c_v5[4] = {0,4,7,3};
/* polygon adjacencies                                                    */
static int c_p0[4] = {5,4,3,2}, c_p1[4] = {2,3,4,5}, c_p2[4] = {0,3,1,5};
static int c_p3[4] = {0,4,1,2}, c_p4[4] = {0,5,1,3}, c_p5[4] = {2,1,4,0};
/* local coordinates of the vertices                                      */
static  double c_c0[3] = {0.,0.,0.}, c_c1[3] = {1.,0.,0.};
static  double c_c2[3] = {1.,1.,0.}, c_c3[3] = {0.,1.,0.};
static  double c_c4[3] = {0.,0.,1.}, c_c5[3] = {1.,0.,1.};
static  double c_c6[3] = {1.,1.,1.}, c_c7[3] = {0.,1.,1.};

static double *cube_coord[8] = {c_c0,c_c1,c_c2,c_c3,c_c4,c_c5,c_c6,c_c7};
static int   cube_polygon_length[6] = {4,4,4,4,4,4};
static int   *cube_vertex[6] =        {c_v0,c_v1,c_v2,c_v3,c_v4,c_v5};
static int   *cube_next_polygon[6] =  {c_p0,c_p1,c_p2,c_p3,c_p4,c_p5};

The procedure which gives adjacency information on the grid can be implemented in the following way:

static ELEMENT3D  *cube_neighbour(ELEMENT3D *el, int pn, int flag,
                                   double *coord, float *xyz)
{
  CUBE      *cd;
  USER_DATA *ud;

  cd = (CUBE*)el->user_data;     /* triple index of the actual element */
  ud = (USER_DATA*)el->mesh->user_data;   /* user data         */

  switch(pn) {
    case 0 :
      if (cd->i_z == 0) return(NIL);
      cd->i_z--;
      fill_points(el);
      coord[2] = 1.0;
      return(el);
    case 1:
      if (cd->i_z == ud->n_z-2) return(NIL);
      cd->i_z++;
      fill_points(el);
      coord[2] = 0.0;
      return(el);
    case 2:
      ...
  }
}

The description of a cube then simply is:

static ELEMENT3D_DESCRIPTION cube_descr =
{
  8 , 6 , cube_polygon_length, cube_vertex, cube_next_polygon,
  3, cube_coord, 1, 
  cube_world_to_coord, cube_coord_to_world, cube_check_inside,
  cube_neighbour, cube_boundary
};

The basic functions of MESH3D are first_element() and next_element(). The first one has to allocate memory for one single element and for the triple index of the ``user's element''.

static ELEMENT3D *cube_first_element(MESH3D *mesh)
{
  ELEMENT3D  *el;
  CUBE       *cd;

  if (mesh == NIL)  return(NIL);
  if ((el = (ELEMENT3D*)malloc(sizeof(ELEMENT3D))) == NIL)  return(NIL);

  el->mesh = mesh;
  el->vertex = (float**)malloc(8*sizeof(float*));
  el->vindex = (int*)malloc(8*sizeof(int));
  el->user_data = (void*)(cd = (CUBE*)malloc(sizeof(CUBE)));
  el->size_of_user_data = sizeof(CUBE);
  el->descr = &cube_descr;

  cd->i_x = cd->i_y = cd->i_z = 0;
  fill_points(el);
  return(el);
}

static ELEMENT3D *cube_next_element(ELEMENT3D *el)
{
  CUBE      *cd;
  USER_DATA *ud;
  if (el == NIL) return(NIL);
  cd = (CUBE*)el->user_data;
  ud = (USER_DATA*)el->mesh->user_data;

  if(++cd->i_x == ud->n_x-1) {
    cd->i_x = 0;
    if(++cd->i_y == ud->n_y-1) {
      cd->i_y = 0;
      if(++cd->i_z == ud->n_z-1) {
        cube_free_element(el);   /* el was last element of the grid       */
        return(NIL);
      }
    }
  }
  fill_points(el);
  return(el);
}

For the visualization of given data from a numerical method we have to implement a procedure f() with corresponding f_el_info() in an F_DATA3D structure. Usually, Finite Difference methods only give values at the nodes of the grid. Such data can be interpolated by a piecewise trilinear function. In the concept of local coordinates it is very easy to evaluate such a function. Let u[0]u[7] be those values at the vertices of one cube of the grid (for the sake of simplicity we assume that these are scalar values). The main part of the function f(ELEMENT3D *el, int ind, double coord[], double val[]) would be

if (coord == NIL)
  *val = u[ind];
else {
  *val = (1-coord[2]) * ((1-coord[1])*((1-coord[0])*u[0]+coord[0]*u[1])+
                             coord[1]*((1-coord[0])*u[3]+coord[0]*u[2]));
  *val +=   coord[2]  * ((1-coord[1])*((1-coord[0])*u[4]+coord[0]*u[5])+
                             coord[1]*((1-coord[0])*u[7]+coord[0]*u[6]));
}

The values of u[0]u[7] can easily be obtained from the user arrays (possibly stored at el->mesh->user_data) and some index arithmetic.

In the second example we give definitions for local_coord, polygon_length, polygon_vertex, and polygon_neighbour for tetrahedral and prismatic elements. Since usually these elements occur in methods using unstructured grids the adjacency connectivity of the elements usually can not be computed by indexing operations as described above. But for the Finite Element / Finite Volume methods itself the underlying user data structure usually provides such information. Thus the functions *neighbour() and *boundary() can be implemented very easily using the user's data structure.

/* vertex indices of the polygons for a prism and a tetrahedron           */
static int p_v0[4] = {1,2,5,4}, p_v1[4] = {0,3,5,2}, p_v2[4] = {0,1,4,3};
static int p_v3[3] = {0,2,1},   p_v4[3] = {3,4,5};
static int t_v0[3] = {1,2,3},   t_v1[3] = {2,0,3};
static int t_v2[3] = {0,1,3},   t_v3[3] = {2,1,0};
/* polygon adjacencies  for a prism and a tetrahedron                     */
static int p_p0[4] = {3,1,4,2}, p_p1[4] = {2,4,0,3}, p_p2[4] = {3,0,4,1};
static int p_p3[3] = {1,0,2},   p_p4[3] = {2,0,1};
static int t_p0[3] = {3,1,2},   t_p1[3] = {3,2,0};
static int t_p2[3] = {3,0,1},   t_p3[3] = {0,2,1};
/* local coordinates of the vertices for a prism and a tetrahedron        */
static  double p_c0[4] = {1.,0.,0.,0.}, p_c1[4] = {0.,1.,0.,0.};
static  double p_c2[4] = {0.,0.,1.,0.}, p_c3[4] = {1.,0.,0.,1.};
static  double p_c4[4] = {0.,1.,0.,1.}, p_c5[4] = {0.,0.,1.,1.};
static  double t_c0[4] = {1.,0.,0.,0.}, t_c1[4] = {0.,1.,0.,0.};
static  double t_c2[4] = {0.,0.,1.,0.}, t_c3[4] = {0.,0.,0.,1.};

static int   prism_polygon_length[5] = {4,4,4,3,3};
static int   *prism_vertex[5] =       {p_v0,p_v1,p_v2,p_v3,p_v4};
static int   *prism_next_polygon[5] = {p_p0,p_p1,p_p2,p_p3,p_p4,};
static double *prism_coord[6] =        {p_c0,p_c1,p_c2,p_c3,p_c4,p_c5};

static int   tetra_polygon_length[4] = {3, 3, 3, 3};
static int   *tetra_vertex[4] =       {t_v0,t_v1,t_v2,t_v3};
static int   *tetra_next_polygon[4] = {t_p0,t_p1,t_p2,t_p3};
static double *tetra_coord[4] =        {t_c0,t_c1,t_c2,t_c3};

The ELEMENT3D_DESCRIPTION for a prism and a tetrahedron then would be

static ELEMENT3D_DESCRIPTION prism_descr =
{
  6 , 5 , prism_polygon_length, prism_vertex, prism_next_polygon,
  4, prism_coord, 1, 
  prism_world_to_coord, prism_coord_to_world, prism_check_inside,
  prism_neighbour,  prism_boundary
};
static ELEMENT3D_DESCRIPTION tetra_descr =
{
  4 , 4 , tetra_polygon_length, tetra_vertex, tetra_next_polygon,
  4, tetra_coord, 1,
  tetra_world_to_coord, tetra_coord_to_world, tetra_check_inside,
  tetra_neighbour,  tetra_boundary
};