text: /*****************************************************************************
   * Copyright (c) 2019 FrontISTR Commons
   * This software is released under the MIT License, see LICENSE.txt
   *****************************************************************************/
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  #include <assert.h>
  #include <errno.h>
  
  #include "hecmw_msgno.h"
  #include "hecmw_struct.h"
  #include "hecmw_error.h"
  #include "hecmw_comm.h"
  
  #include "hecmw_couple_define.h"
  #include "hecmw_couple_control.h"
  #include "hecmw_couple_boundary_info.h"
  #include "hecmw_couple_bounding_box.h"
  #include "hecmw_couple_background_cell.h"
  
  #define INFINITE (1.0E+37)
  
  #define EPS (1.0E-06)
  
  /*================================================================================================*/
  
  extern void HECMW_couple_free_background_cell(
      struct hecmw_couple_background_cell *bgcell) {
    if (bgcell == NULL) return;
  
    HECMW_free(bgcell);
    bgcell = NULL;
  }
  
  static struct hecmw_couple_background_cell *alloc_struct_bgcell(void) {
    struct hecmw_couple_background_cell *bgcell = NULL;
    int size;
  
    size   = sizeof(struct hecmw_couple_background_cell);
    bgcell = (struct hecmw_couple_background_cell *)HECMW_malloc(size);
    if (bgcell == NULL) {
      HECMW_set_error(errno, "");
      return NULL;
    }
  
    bgcell->n    = 0;
    bgcell->coef = 0.0;
    bgcell->nx   = 0;
    bgcell->ny   = 0;
    bgcell->nz   = 0;
    bgcell->dx   = 0.0;
    bgcell->dy   = 0.0;
    bgcell->dz   = 0.0;
  
    return bgcell;
  }
  
  static int elem_size_by_elem(const struct hecmwST_local_mesh *mesh,
                               const struct hecmw_couple_boundary *boundary,
                               double *min_dx, double *min_dy, double *min_dz,
                               double *max_dx, double *max_dy, double *max_dz) {
    double min_x, min_y, min_z, max_x, max_y, max_z, dx, dy, dz, coord_x, coord_y,
        coord_z;
    int elem, node, i, j;
  
    *min_dx = *min_dy = *min_dz = INFINITE;
    *max_dx = *max_dy = *max_dz = 0.0;
  
    for (i = 0; i < boundary->elem->n; i++) {
      elem = boundary->elem->item[i];
  
      min_x = min_y = min_z = +INFINITE;
      max_x = max_y = max_z = -INFINITE;
  
      for (j = mesh->elem_node_index[elem - 1]; j < mesh->elem_node_index[elem];
           j++) {
        node    = mesh->elem_node_item[j];
        coord_x = mesh->node[3 * (node - 1)];
        coord_y = mesh->node[3 * (node - 1) + 1];
        coord_z = mesh->node[3 * (node - 1) + 2];
  
        if (coord_x < min_x) min_x = coord_x;
        if (coord_y < min_y) min_y = coord_y;
        if (coord_z < min_z) min_z = coord_z;
        if (coord_x > max_x) max_x = coord_x;
        if (coord_y > max_y) max_y = coord_y;
        if (coord_z > max_z) max_z = coord_z;
      }
  
      dx = max_x - min_x;
      dy = max_y - min_y;
      dz = max_z - min_z;
  
      if (dx < *min_dx) *min_dx = dx;
      if (dy < *min_dy) *min_dy = dy;
      if (dz < *min_dz) *min_dz = dz;
      if (dx > *max_dx) *max_dx = dx;
      if (dy > *max_dy) *max_dy = dy;
      if (dz > *max_dz) *max_dz = dz;
    }
  
    return 0;
  }
  
  static int elem_size_by_surf(const struct hecmwST_local_mesh *mesh,
                               const struct hecmw_couple_boundary *boundary,
                               double *min_dx, double *min_dy, double *min_dz,
                               double *max_dx, double *max_dy, double *max_dz) {
    double min_x, min_y, min_z, max_x, max_y, max_z, dx, dy, dz, coord_x, coord_y,
        coord_z;
    int elem, node, i, j;
  
    *min_dx = *min_dy = *min_dz = INFINITE;
    *max_dx = *max_dy = *max_dz = 0.0;
  
    for (i = 0; i < boundary->surf->n; i++) {
      elem = boundary->surf->item[2 * i];
  
      min_x = min_y = min_z = +INFINITE;
      max_x = max_y = max_z = -INFINITE;
  
      for (j = mesh->elem_node_index[elem - 1]; j < mesh->elem_node_index[elem];
           j++) {
        node    = mesh->elem_node_item[j];
        coord_x = mesh->node[3 * (node - 1)];
        coord_y = mesh->node[3 * (node - 1) + 1];
        coord_z = mesh->node[3 * (node - 1) + 2];
  
        if (coord_x < min_x) min_x = coord_x;
        if (coord_y < min_y) min_y = coord_y;
        if (coord_z < min_z) min_z = coord_z;
        if (coord_x > max_x) max_x = coord_x;
        if (coord_y > max_y) max_y = coord_y;
        if (coord_z > max_z) max_z = coord_z;
      }
  
      dx = max_x - min_x;
      dy = max_y - min_y;
      dz = max_z - min_z;
  
      if (dx < *min_dx) *min_dx = dx;
      if (dy < *min_dy) *min_dy = dy;
      if (dz < *min_dz) *min_dz = dz;
      if (dx > *max_dx) *max_dx = dx;
      if (dy > *max_dy) *max_dy = dy;
      if (dz > *max_dz) *max_dz = dz;
    }
  
    return 0;
  }
  
  /*================================================================================================*/
  
  extern struct hecmw_couple_background_cell *HECMW_couple_set_background_cell(
      const char *boundary_id, const struct hecmwST_local_mesh *mesh,
      const struct hecmw_couple_bounding_box *bbox,
      const struct hecmw_couple_boundary *boundary) {
    struct hecmw_couple_background_cell *bgcell;
    double min_dx, min_dy, min_dz, max_dx, max_dy, max_dz, dx, dy, dz;
    double bbox_size_x, bbox_size_y, bbox_size_z;
    int nx, ny, nz;
  
    if (boundary_id == NULL) {
      HECMW_set_error(
          HECMWCPL_E_INVALID_ARG,
          "HECMW_couple_set_background_cell(): 'boundary_id' is NULL");
      return NULL;
    }
    if (mesh == NULL) {
      HECMW_set_error(HECMWCPL_E_INVALID_ARG,
                      "HECMW_couple_set_background_cell(): 'mesh' is NULL");
      return NULL;
    }
    if (bbox == NULL) {
      HECMW_set_error(HECMWCPL_E_INVALID_ARG,
                      "HECMW_couple_set_background_cell(): 'bbox' is NULL");
      return NULL;
    }
    if (boundary == NULL) {
      HECMW_set_error(HECMWCPL_E_INVALID_ARG,
                      "HECMW_couple_set_background_cell(): 'boundary' is NULL");
      return NULL;
    }
  
    if ((bgcell = alloc_struct_bgcell()) == NULL) return NULL;
  
    if (HECMW_couple_ctrl_get_bgcoef(boundary_id, &bgcell->coef) != HECMW_SUCCESS)
      goto error;
  
    if (boundary->geom_type ==
        HECMW_COUPLE_NODE_GROUP) { /* node group		*/
      HECMW_set_error(HECMWCPL_E_NONSUPPORT_GEOMTYPE,
                      "In current version, node group is not supported.");
      goto error;
  
    } else if (boundary->geom_type ==
               HECMW_COUPLE_ELEMENT_GROUP) { /* element group	*/
      HECMW_set_error(HECMWCPL_E_NONSUPPORT_GEOMTYPE,
                      "In current version, element group is not supported.");
      goto error;
      /*		elem_size_by_elem(mesh, boundary, &min_dx, &min_dy, &min_dz, &max_dx,
       * &max_dy, &max_dz); */
  
    } else if (boundary->geom_type ==
               HECMW_COUPLE_SURFACE_GROUP) { /* surface group	*/
      elem_size_by_surf(mesh, boundary, &min_dx, &min_dy, &min_dz, &max_dx,
                        &max_dy, &max_dz);
  
    } else { /* error			*/
      HECMW_set_error(HECMWCPL_E_INVALID_GEOMTYPE, "");
      goto error;
    }
  
    bbox_size_x = bbox->enlarged->max_x - bbox->enlarged->min_x;
    bbox_size_y = bbox->enlarged->max_y - bbox->enlarged->min_y;
    bbox_size_z = bbox->enlarged->max_z - bbox->enlarged->min_z;
  
    if (max_dx * bgcell->coef > bbox->tolerance + EPS) {
      dx = max_dx * bgcell->coef;
    } else {
      dx = bbox->tolerance + EPS;
    }
    if (max_dy * bgcell->coef > bbox->tolerance + EPS) {
      dy = max_dy * bgcell->coef;
    } else {
      dy = bbox->tolerance + EPS;
    }
    if (max_dz * bgcell->coef > bbox->tolerance + EPS) {
      dz = max_dz * bgcell->coef;
    } else {
      dz = bbox->tolerance + EPS;
    }
  
    nx = bbox_size_x / dx;
    ny = bbox_size_y / dy;
    nz = bbox_size_z / dz;
  
    bgcell->nx = (nx > 0) ? nx : 1;
    bgcell->ny = (ny > 0) ? ny : 1;
    bgcell->nz = (nz > 0) ? nz : 1;
  
    bgcell->dx = bbox_size_x / bgcell->nx;
    bgcell->dy = bbox_size_y / bgcell->ny;
    bgcell->dz = bbox_size_z / bgcell->nz;
  
    bgcell->n = bgcell->nx * bgcell->ny * bgcell->nz;
  
    return bgcell;
  
  error:
    HECMW_couple_free_background_cell(bgcell);
    return NULL;
  }
  