cs_cdo_quantities.h

Go to the documentation of this file.

#ifndef __CS_CDO_QUANTITIES_H__
#define __CS_CDO_QUANTITIES_H__

/*============================================================================
 * Manage geometrical quantities needed in CDO schemes
 *============================================================================*/

/*
  This file is part of Code_Saturne, a general-purpose CFD tool.

  Copyright (C) 1998-2021 EDF S.A.

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

  This program 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 General Public License for more
  details.

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

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

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

#include "cs_base.h"
#include "cs_cdo_connect.h"
#include "cs_flag.h"
#include "cs_math.h"
#include "cs_mesh.h"
#include "cs_mesh_quantities.h"

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

BEGIN_C_DECLS

/*============================================================================
 * Macro definitions
 *============================================================================*/

/*============================================================================
 * Type definitions
 *============================================================================*/

typedef enum {

  /* Set of geometrical quantities related to CDO schemes */

  CS_CDO_QUANTITIES_EB_SCHEME                  = 1<<0,  /* =   1 */
  CS_CDO_QUANTITIES_FB_SCHEME                  = 1<<1,  /* =   2 */
  CS_CDO_QUANTITIES_HHO_SCHEME                 = 1<<2,  /* =   4 */
  CS_CDO_QUANTITIES_VB_SCHEME                  = 1<<3,  /* =   8 */
  CS_CDO_QUANTITIES_VCB_SCHEME                 = 1<<4,  /* =  16 */

  /* How to compute the cell center
     ------------------------------ */

  CS_CDO_QUANTITIES_MEANV_CENTER               = 1<<7,  /* = 128 */
  CS_CDO_QUANTITIES_BARYC_CENTER               = 1<<8,  /* = 256 */
  CS_CDO_QUANTITIES_SATURNE_CENTER             = 1<<9,  /* = 512 */

} cs_cdo_quantities_bit_t;


/* Structure storing information about variation of entities across the
   mesh for a given type of entity (cell, face and edge) */
typedef struct {

  /* Measure is either a volume for cells, a surface for faces or a length
     for edges */
  double   meas_min;  /* Min. value of the entity measure */
  double   meas_max;  /* Max. value of the entity measure */
  double   h_min;     /* Estimation of the min. value of the diameter */
  double   h_max;     /* Estimation of the max. value of the diameter */

} cs_quant_info_t;

/* For primal vector quantities (edge or face) */
typedef struct {

  double  meas;       /* length or area */
  double  unitv[3];   /* unitary vector: tangent or normal to the element */
  double  center[3];

} cs_quant_t;

typedef struct { /* Specific mesh quantities */

  /* Flag storing which quantities to compute and how to define the cell
   * center */
  cs_flag_t        flag;

  /* Global mesh quantities */
  double           vol_tot;

  /* Cell-based quantities */
  /* ===================== */

  cs_lnum_t         n_cells;        /* Local number of cells */
  cs_gnum_t         n_g_cells;      /* Global number of cells */
  cs_real_t        *cell_centers;   /* May be shared according to options */
  const cs_real_t  *cell_vol;       /* Shared with cs_mesh_quantities_t */

  cs_quant_info_t   cell_info;

  /* Face-based quantities */
  /* ===================== */

  cs_lnum_t         n_faces;        /* n_i_faces + n_b_faces */
  cs_gnum_t         n_g_faces;      /* Global number of faces */

  /* cs_quant_t structure attached to a face (interior or border) can be built
     on-the-fly: cs_quant_get(flag, f_id, quant). In order to reduce the memory
     consumption one shares face quantities with the ones defined in the legacy
     part and stored in the cs_mesh_quantities_t structure */

  cs_lnum_t         n_i_faces;      /* Local number of interior faces */
  const cs_real_t  *i_face_normal;  /* Shared with cs_mesh_quantities_t */
  const cs_real_t  *i_face_center;  /* Shared with cs_mesh_quantities_t */
  const cs_real_t  *i_face_surf;    /* Shared with cs_mesh_quantities_t */

  cs_lnum_t         n_b_faces;      /* Local number of border faces */
  const cs_real_t  *b_face_normal;  /* Shared with cs_mesh_quantities_t */
  const cs_real_t  *b_face_center;  /* Shared with cs_mesh_quantities_t */
  const cs_real_t  *b_face_surf;    /* Shared with cs_mesh_quantities_t */

  /* cs_nvec3_t structure attached to a dual edge can be built on-the-fly to
     access to its length and its unit tangential vector. One
     recalls that a dual edge is associated to a primal face and is shared with
     two cells if the face lies inside the domain and shared with one cell if
     the face lies on the domain boundary.
     Scan this quantity with the c2f connectivity
  */

  cs_real_t        *dedge_vector;   /* Allocation to 3*c2f->idx[n_faces] */

  cs_real_t        *pvol_fc;        /* Portion of volume surrounding a face
                                     * in each cell. This is a pyramid of
                                     * base the face and apex the cell center
                                     * Scanned with the c2f adjacency.
                                     * Not always allocated.
                                     */
  cs_quant_info_t   face_info;

  /* Edge-based quantities */
  /* ===================== */

  cs_lnum_t         n_edges;        /* Local number of edges */
  cs_gnum_t         n_g_edges;      /* Global number of edges */

  cs_real_t        *edge_vector;    /* Allocation to 3*n_edges
                                       Norm of the vector is equal to the
                                       distance between two vertices.
                                       Unit vector is the tangential direction
                                       attached to the edge */

  /* For each edge e belonging to a cell c, the dual face is built from the
     contributions of two triangles s(x_c, x_f, x_e) and s(x_c, x_f', x_e) with
     the faces f and f' belonging to F_e \cap F_c
     Scan this quantity with the c2e connectivity */

  cs_real_t        *dface_normal;   /* Vector-valued normal for each dual face
                                     * inside a cell associated to an edge */
  cs_real_t        *pvol_ec;        /* Portion of volume surrounding an edge
                                     * in each cell. Scanned with the c2e
                                     * adjacency.
                                     * Not always allocated. */

  cs_quant_info_t   edge_info;

  /* Vertex-based quantities */
  /* ======================= */

  cs_lnum_t         n_vertices;      /* Local number of vertices */
  cs_gnum_t         n_g_vertices;    /* Global number of vertices */

  cs_real_t        *dcell_vol;       /* Dual volume related to each vertex.
                                      * Scanned with the c2v adjacency.
                                      * Not always allocated.
                                      */
  const cs_real_t  *vtx_coord;       /* Shared with the cs_mesh_t structure */

} cs_cdo_quantities_t;

/*============================================================================
 * Global variables
 *============================================================================*/

/*============================================================================
 * Public function prototypes
 *============================================================================*/

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

static inline double
cs_compute_area_from_quant(const cs_quant_t    qa,
                           const cs_real_t    *xb)
{
  const double  xab[3] = {xb[0] - qa.center[0],
                          xb[1] - qa.center[1],
                          xb[2] - qa.center[2]};
  const double  cp[3] = {qa.unitv[1]*xab[2] - qa.unitv[2]*xab[1],
                         qa.unitv[2]*xab[0] - qa.unitv[0]*xab[2],
                         qa.unitv[0]*xab[1] - qa.unitv[1]*xab[0]};

  return 0.5 * qa.meas * cs_math_3_norm(cp);
}

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

void
cs_cdo_quantities_set(cs_flag_t   option_flag);

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

void
cs_cdo_quantities_set_algo_ccenter(cs_cdo_quantities_bit_t   algo);

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

cs_cdo_quantities_t *
cs_cdo_quantities_build(const cs_mesh_t             *m,
                        const cs_mesh_quantities_t  *mq,
                        const cs_cdo_connect_t      *topo);

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

cs_cdo_quantities_t *
cs_cdo_quantities_free(cs_cdo_quantities_t   *q);

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

void
cs_cdo_quantities_summary(const cs_cdo_quantities_t  *quant);

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

void
cs_cdo_quantities_dump(const cs_cdo_quantities_t  *cdoq);

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

void
cs_cdo_quantities_compute_pvol_fc(const cs_cdo_quantities_t    *cdoq,
                                  const cs_adjacency_t         *c2f,
                                  cs_real_t                   **p_pvol_fc);

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

void
cs_cdo_quantities_compute_pvol_ec(const cs_cdo_quantities_t   *cdoq,
                                  const cs_adjacency_t        *c2e,
                                  cs_real_t                  **p_pvol_ec);

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

void
cs_cdo_quantities_compute_dual_volumes(const cs_cdo_quantities_t   *cdoq,
                                       const cs_adjacency_t        *c2v,
                                       cs_real_t                   *dual_vol);

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

void
cs_cdo_quantities_compute_i_tef(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     f_id,
                                cs_real_t                     tef[]);

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

void
cs_cdo_quantities_compute_b_tef(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     bf_id,
                                cs_real_t                     tef[]);

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

void
cs_cdo_quantities_compute_i_wvf(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     f_id,
                                cs_real_t                     wvf[]);

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

void
cs_cdo_quantities_compute_b_wvf(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     bf_id,
                                cs_real_t                     wvf[]);

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

inline static const cs_real_t *
cs_quant_get_face_vector_area(cs_lnum_t                    f_id,
                              const cs_cdo_quantities_t   *cdoq)
{
  if (f_id < cdoq->n_i_faces)   /* Interior face */
    return cdoq->i_face_normal + 3*f_id;
  else                          /* Border face */
    return cdoq->b_face_normal + 3*(f_id - cdoq->n_i_faces);
}

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

inline static const cs_real_t *
cs_quant_get_face_center(cs_lnum_t                    f_id,
                         const cs_cdo_quantities_t   *cdoq)
{
  if (f_id < cdoq->n_i_faces)   /* Interior face */
    return cdoq->i_face_center + 3*f_id;
  else                          /* Border face */
    return cdoq->b_face_center + 3*(f_id - cdoq->n_i_faces);
}

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

cs_quant_t
cs_quant_set_face(cs_lnum_t                    f_id,
                  const cs_cdo_quantities_t   *cdoq);

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

cs_nvec3_t
cs_quant_set_face_nvec(cs_lnum_t                    f_id,
                       const cs_cdo_quantities_t   *cdoq);

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

cs_nvec3_t
cs_quant_set_edge_nvec(cs_lnum_t                    e_id,
                       const cs_cdo_quantities_t   *cdoq);

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

cs_nvec3_t
cs_quant_set_dedge_nvec(cs_lnum_t                     f_shift,
                        const cs_cdo_quantities_t    *cdoq);

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

void
cs_quant_dump(FILE             *f,
              cs_lnum_t         num,
              const cs_quant_t  q);

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

END_C_DECLS

#endif /* __CS_CDO_QUANTITIES_H__ */