/*============================================================================
 * Deal with the exchange of data included in opaque structure fvm_nodal_t
 *============================================================================*/

/*
  This file is part of the "Finite Volume Mesh" library, intended to provide
  finite volume mesh and associated fields I/O and manipulation services.

  Copyright (C) 2006  EDF

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/*----------------------------------------------------------------------------
 * BFT library headers
 *----------------------------------------------------------------------------*/

#include <bft_mem.h>
#include <bft_printf.h>

/*----------------------------------------------------------------------------
 *  Local headers
 *----------------------------------------------------------------------------*/

#include <fvm_config_defs.h>
#include <fvm_defs.h>
#include <fvm_io_num.h>
#include <fvm_parall.h>

/*----------------------------------------------------------------------------
 *  Header for the current file
 *----------------------------------------------------------------------------*/

#include <fvm_nodal.h>
#include <fvm_nodal_priv.h>
#include <fvm_nodal_extract.h>

/*----------------------------------------------------------------------------*/

#ifdef __cplusplus
extern "C" {
#if 0
} /* Fake brace to force back Emacs auto-indentation back to column 0 */
#endif
#endif /* __cplusplus */

/*============================================================================
 * Static global variables
 *============================================================================*/

/*============================================================================
 * Private function definitions
 *============================================================================*/

/*============================================================================
 * Semi-private function definitions (prototypes in fvm_nodal_priv.h)
 *============================================================================*/

/*============================================================================
 * Public function definitions
 *============================================================================*/

/*----------------------------------------------------------------------------
 * Copy global vertex numbers to an array.
 *
 * parameters:
 *   this_nodal <-- pointer to nodal mesh structure
 *   g_vtx_num  --> global vertex numbers (pre-allocated)
 *----------------------------------------------------------------------------*/

void
fvm_nodal_get_global_vertex_num(const fvm_nodal_t  *this_nodal,
                                fvm_gnum_t         *g_vtx_num)
{
  size_t size;
  fvm_lnum_t  vertex_id;
  fvm_io_num_t *global_io_num = this_nodal->global_vertex_num;

  assert(g_vtx_num != NULL || this_nodal->n_vertices == 0);

  if (g_vtx_num == NULL)
    return;

  if (global_io_num == NULL) {

    for (vertex_id = 0; vertex_id < this_nodal->n_vertices; vertex_id++)
      g_vtx_num[vertex_id] = vertex_id + 1;

  }
  else {

    size = sizeof(fvm_gnum_t) * fvm_io_num_get_local_count(global_io_num);
    memcpy(g_vtx_num, fvm_io_num_get_global_num(global_io_num), size);

  }
}

/*----------------------------------------------------------------------------
 * Copy global element numbers of a given element type to an array.
 *
 * Note that if the mesh contains multiple sections of the same element type,
 * the global element numbers are continued from one section to the next,
 * so to the user, all is as if the sections were concatenated.
 *
 * parameters:
 *   this_nodal   <-- pointer to nodal mesh structure
 *   element_type <-- type of elements to deal with
 *   g_elt_num    <-> pointer to global_element_num array (pre-allocated)
 *
 * returns:
 *----------------------------------------------------------------------------*/

void
fvm_nodal_get_global_element_num(const fvm_nodal_t  *this_nodal,
                                 fvm_element_t       element_type,
                                 fvm_gnum_t         *g_elt_num)
{
  int element_id, section_id;

  fvm_lnum_t n_elements = 0, element_count = 0;
  fvm_gnum_t n_g_elements = 0, global_count = 0;
  const fvm_gnum_t *g_num = NULL;

  /* Define global element numbers */

  for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {

    const fvm_nodal_section_t *section = this_nodal->sections[section_id];

    if (section->type == element_type) {

      const fvm_io_num_t *io_num = section->global_element_num;

      if (io_num == NULL) {

        for (element_id = 0; element_id < section->n_elements; element_id++)
          g_elt_num[element_count + element_id] = element_id + global_count + 1;

        element_count += section->n_elements;
        global_count += section->n_elements;

      }

      else { /* io_num != NULL */

        n_elements = fvm_io_num_get_local_count(io_num);
        n_g_elements = fvm_io_num_get_global_count(io_num);
        g_num = fvm_io_num_get_global_num(io_num);

        if (element_count == 0) {
          memcpy((char *)g_elt_num,
                 g_num,
                 sizeof(fvm_gnum_t) * n_elements);

        }
        else {

          for (element_id = 0; element_id < n_elements; element_id++)
            g_elt_num[element_count + element_id]
              = g_num[element_id] + global_count;

          element_count += n_elements;
          global_count += n_g_elements;

        }

      }

    } /* section->type == element_type */

  } /* End of loop on sections */

  return;
}

/*----------------------------------------------------------------------------
 * Copy vertex coordinates to an array.
 *
 * parameters:
 *   this_nodal     <-- pointer to nodal mesh structure
 *   interlace      <-- indicates if destination array is interlaced
 *   vertex_coords  --> vertices coordinates (pre-allocated)
 *----------------------------------------------------------------------------*/

void
fvm_nodal_get_vertex_coords(const fvm_nodal_t  *this_nodal,
                            fvm_interlace_t     interlace,
                            fvm_coord_t        *vertex_coords)
{
  int i;
  fvm_lnum_t vertex_id;

  const int dim = this_nodal->dim;
  const fvm_lnum_t n_vertices = this_nodal->n_vertices;
  const fvm_coord_t *coords = this_nodal->vertex_coords;
  const fvm_lnum_t  *parent_num = this_nodal->parent_vertex_num;

  if (this_nodal->parent_vertex_num == NULL) {

    if (interlace == FVM_INTERLACE)
      memcpy(vertex_coords, coords, sizeof(fvm_coord_t) * n_vertices * dim);

    else {
      for (i = 0; i < dim; i++) {
        for (vertex_id = 0; vertex_id < n_vertices; vertex_id++)
          vertex_coords[n_vertices*i + vertex_id]
            = coords[vertex_id*dim + i];
      }

    }

  }
  else { /* parent_vertex_num != NULL */

    if (interlace == FVM_INTERLACE) {

      for (i = 0; i < dim; i++) {
        for (vertex_id = 0; vertex_id < n_vertices; vertex_id++)
          vertex_coords[vertex_id * dim + i]
            = coords[(parent_num[vertex_id]-1) * dim + i];
      }

    }
    else {

      for (i = 0; i < dim; i++) {
        for (vertex_id = 0; vertex_id < n_vertices; vertex_id++)
          vertex_coords[n_vertices*i + vertex_id]
            = coords[(parent_num[vertex_id]-1) * dim + i];
      }

    }

  }

}

/*----------------------------------------------------------------------------
 * Copy element centers to an array.
 *
 * Note that if the mesh contains multiple cell element sections of, the
 * cell_centers array spans all sections, so to the user, all is as if the
 * sections were concatenated.
 *
 * parameters:
 *   this_nodal     <-- pointer to nodal mesh structure
 *   interlace      <-- indicates if destination array is interlaced
 *   entity_dim     <-- dimension of entities we want to count (0 to 3)
 *   cell_centers   --> cell centers coordinates (pre-allocated)
 *----------------------------------------------------------------------------*/

void
fvm_nodal_get_element_centers(const fvm_nodal_t  *this_nodal,
                              fvm_interlace_t     interlace,
                              int                 entity_dim,
                              fvm_coord_t        *cell_centers)
{
  int i, j;
  int section_id;
  fvm_lnum_t element_id;

  fvm_lnum_t element_count = 0;

  const int dim = this_nodal->dim;
  const fvm_coord_t *coords = this_nodal->vertex_coords;
  const fvm_lnum_t  *parent_num = this_nodal->parent_vertex_num;
  const fvm_lnum_t n_elements = fvm_nodal_get_n_entities(this_nodal,
                                                         entity_dim);

  for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {

    const fvm_nodal_section_t *section = this_nodal->sections[section_id];

    if (section->entity_dim == entity_dim) {

      if (section->stride != 0) {

        const int stride = section->stride;

        for (element_id = 0;
             element_id < section->n_elements;
             element_id++, element_count++) {

          double  cell_center[3] = {0., 0., 0.};
          double  denom = 0.;

          if (this_nodal->parent_vertex_num == NULL) {
            for (j = 0; j < stride; j++) {
              const fvm_lnum_t vertex_id
                = section->vertex_num[(element_id * stride) + j] - 1;
              for (i = 0; i < dim; i++)
                cell_center[i] += coords[vertex_id*dim + i];
              denom += 1.;
            }
          }
          else { /* if (this_nodal->parent_vertex_num != NULL) */
            for (j = 0; j < stride; j++) {
              const fvm_lnum_t vertex_id
                = section->vertex_num[(element_id * stride) + j] - 1;
              for (i = 0; i < dim; i++)
                cell_center[i] += coords[(parent_num[vertex_id]-1)*dim + i];
              denom += 1.;
            }
          }

          if (interlace == FVM_INTERLACE) {
            for (i = 0; i < dim; i++)
              cell_centers[element_count*dim + i] = cell_center[i] / denom;
          }
          else {
            for (i = 0; i < dim; i++)
              cell_centers[n_elements*i + element_count]
                = cell_center[i] / denom;
          }

        }

      }

      assert(section->stride != 0); /* Not implemented yet for non strided elements */

    } /* section->entity_dim == entity_dim */

  } /* section_id */

}

/*----------------------------------------------------------------------------
 * Copy element -> vertex connectivity of a given element type to an array.
 *
 * Note that if the mesh contains multiple sections of the same element type,
 * the connectivity spans all sections, so to the user, all is as if the
 * sections were concatenated.
 *
 * parameters:
 *   this_nodal    <-- pointer to nodal mesh structure
 *   element_type  <-- type of elements of the section to deal with
 *   connectivity  <-> pointer to connectvity (pre-allocated)
 *----------------------------------------------------------------------------*/

void
fvm_nodal_get_strided_connect(const fvm_nodal_t  *this_nodal,
                              fvm_element_t       element_type,
                              fvm_lnum_t         *connectivity)
{
  int i, section_id;
  fvm_lnum_t element_id;
  fvm_lnum_t element_count = 0;
  fvm_gnum_t global_count = 0;

  /* Verify section with "element_type" type exists and is strided */

  if (element_type == FVM_FACE_POLY || element_type == FVM_CELL_POLY)
    bft_error(__FILE__, __LINE__, 0,
              _("Elements of type : \"%s\" are not strided elements.\n"
                "Incorrect use with fvm_nodal_get_strided_connect()\n"
                "Associated nodal mesh : \"%s\"\n"),
              fvm_elements_type_name[element_type], this_nodal->name);

  /* Retrieve connectivity */

  for (section_id = 0; section_id < this_nodal->n_sections; section_id++) {

    const fvm_nodal_section_t *section = this_nodal->sections[section_id];

    if (section->type == element_type) {

      const int stride = section->stride;
      const fvm_lnum_t *num = section->vertex_num;

      for (element_id = 0; element_id < section->n_elements; element_id++) {
        for (i = 0; i < stride; i++) {

          connectivity[element_id * stride + i + element_count]
            = global_count + num[element_id * stride + i];

        }
      }

      global_count += section->n_elements;
      element_count += section->n_elements * stride;

    } /* Section->type == element_type */

  } /* End of loop on sections */

}

/*----------------------------------------------------------------------------*/

#ifdef __cplusplus
}
#endif /* __cplusplus */


syntax highlighted by Code2HTML, v. 0.9.1