cs_cdo_local.h

Go to the documentation of this file.

#ifndef __CS_CDO_LOCAL_H__
#define __CS_CDO_LOCAL_H__

/*============================================================================
 * Routines to handle low-level routines related to CDO local quantities:
 * - local matrices (stored in dense format),
 * - local quantities related to a cell.
 *============================================================================*/

/*
  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_cdo_connect.h"
#include "cs_cdo_quantities.h"
#include "cs_flag.h"
#include "cs_param_cdo.h"
#include "cs_sdm.h"

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

BEGIN_C_DECLS

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

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

typedef struct {

  cs_real_t     t_pty_eval;
  cs_real_t     t_bc_eval;
  cs_real_t     t_st_eval;
  cs_flag_t     cell_flag;

  /* Store the cellwise value for the grad-div, the time and reaction properties
   * since the associated Hodge operator is linked to the unity as related
   * property */

  double        gpty_val;  
  double        tpty_val;  
  double        rpty_vals[CS_CDO_N_MAX_REACTIONS];
  double        rpty_val;  
  /* Advection-related values */
  double       *adv_fluxes;

  /* Temporary buffers (erase and updated several times during the system
     build */
  int          *ids;     
  double       *values;  
  cs_real_3_t  *vectors; 
  /* Structures used to build specific terms composing the algebraic system */
  cs_sdm_t     *loc;   
  cs_sdm_t     *aux;   
} cs_cell_builder_t;

typedef struct {

  cs_lnum_t   c_id;       
  int         n_dofs;   
  cs_lnum_t  *dof_ids;  
  cs_flag_t  *dof_flag; 
  cs_sdm_t   *mat;      
  double     *rhs;      
  double     *source;   
  double     *val_n;    
  double     *val_nm1;  
  /* Boundary conditions for the local system */
  short int   n_bc_faces;  
  short int  *_f_ids;      
  cs_lnum_t  *bf_ids;      
  cs_flag_t  *bf_flag;     
  /* Dirichlet BCs */
  bool        has_dirichlet; 
  double     *dir_values;    
  /* Neumann BCs */
  bool        has_nhmg_neumann; 
  double     *neu_values;       
  /* Robin BCs */
  bool        has_robin;     
  double     *rob_values;    
  /* Sliding BCs */
  bool        has_sliding;   
  /* Internal enforcement of DoFs */
  bool        has_internal_enforcement;  
  cs_lnum_t  *intern_forced_ids;         
} cs_cell_sys_t;

typedef struct {

  cs_eflag_t     flag;    
  fvm_element_t  type;    
  /* Sizes used to allocate buffers */
  short int      n_max_vbyc;
  short int      n_max_ebyc;
  short int      n_max_fbyc;

  /* Cell information */
  cs_lnum_t      c_id;    
  cs_real_3_t    xc;      
  double         vol_c;   
  double         diam_c;  
  /* Vertex information */
  short int    n_vc;  
  cs_lnum_t   *v_ids; 
  double      *xv;    
  double      *wvc;   
  /* Edge information */
  short int    n_ec;   
  cs_lnum_t   *e_ids;  
  cs_quant_t  *edge;   
  cs_nvec3_t  *dface;  
  cs_real_t   *pvol_e; 
  /* Face information */
  short int    n_fc;        
  cs_lnum_t    bface_shift; 
  cs_lnum_t   *f_ids;       
  short int   *f_sgn;       
  double      *f_diam;      
  double      *hfc;         
  cs_quant_t  *face;        
  cs_nvec3_t  *dedge;      
  cs_real_t   *pvol_f;      
  /* Local e2v connectivity: size 2*n_ec (allocated to 2*n_max_ebyc) */
  short int   *e2v_ids;  
  short int   *e2v_sgn;  
  /* Local f2v connectivity: size = 2*n_max_ebyc */
  short int   *f2v_idx;  
  short int   *f2v_ids;  
  /* Local f2e connectivity: size = 2*n_max_ebyc */
  short int   *f2e_idx;  
  short int   *f2e_ids;  
  short int   *f2e_sgn;  
  double      *tef;      
  /* Local e2f connectivity: size 2*n_ec (allocated to 2*n_max_ebyc) */
  short int   *e2f_ids;  
  cs_nvec3_t  *sefc;     
} cs_cell_mesh_t;

typedef struct {

  short int    n_max_vbyf; 
  cs_lnum_t    c_id;    
  cs_real_3_t  xc;      
  /* Face information */
  cs_lnum_t    f_id;   
  short int    f_sgn;  
  cs_quant_t   face;   
  cs_nvec3_t   dedge;  
  double       pvol;   
  double       hfc;    
  /* Vertex information */
  short int    n_vf;    
  cs_lnum_t   *v_ids;   
  double      *xv;      
  double      *wvf;     
  /* Edge information */
  short int    n_ef;    
  cs_lnum_t   *e_ids;   
  cs_quant_t  *edge;    
  double      *tef;     
  /* Local e2v connectivity: size 2*n_ec (allocated to 2*n_max_ebyf) */
  short int   *e2v_ids;  
} cs_face_mesh_t;

/*
   A cs_face_mesh_light_t structure is close to a cs_face_mesh_t structure
   There are less members to be buildt quicker.
   Such structure is always associated to a cs_cell_mesh_t structure
*/

typedef struct {

  short int    n_max_vbyf;  /* Max number of vertices belonging to a face
                               (= n_max_ebyf) */

  cs_lnum_t    c_id;    /* id of related cell in the mesh numbering */
  short int    f;       /* id of the face in the cell mesh numbering */

  /* Vertex information */
  short int    n_vf;    /* local number of vertices on this face */
  short int   *v_ids;   /* vertex ids in the cellwise numbering */
  double      *wvf;     /* weights related to each vertex */

  /* Edge information */
  short int    n_ef;    /* local number of edges on this face (= n_vf) */
  short int   *e_ids;   /* edge ids in the cellwise numbering */
  double      *tef;     /* area of the triangle of base e and apex xf */

} cs_face_mesh_light_t;

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

/* Auxiliary buffers for extra-operations related to local problems. These
 * buffers are also used for computing quantities related to a cs_cell_mesh_t
 * (there are as many buffers as threads since a call to these buffers can be
 * inside an OpenMP directive */
extern int                     cs_cdo_local_d_buffer_size;
extern double                **cs_cdo_local_d_buffer;

/* Pointer of pointers to global structures */
extern cs_cell_mesh_t        **cs_cdo_local_cell_meshes;
extern cs_face_mesh_t        **cs_cdo_local_face_meshes;
extern cs_face_mesh_light_t  **cs_cdo_local_face_meshes_light;

/*============================================================================
 * Static inline function prototypes
 *============================================================================*/

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

static inline int
cs_cdo_local_get_d_buffer_size(void)
{
  return cs_cdo_local_d_buffer_size;
}

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

static inline double *
cs_cdo_local_get_d_buffer(int   thr_id)
{
  assert(thr_id > -1 && thr_id < cs_glob_n_threads);
  return cs_cdo_local_d_buffer[thr_id];
}

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

static inline bool
cs_cell_has_boundary_elements(const cs_cell_builder_t  *const cb)
{
  if ((cb->cell_flag & (CS_FLAG_BOUNDARY_CELL_BY_FACE   |
                        CS_FLAG_BOUNDARY_CELL_BY_VERTEX |
                        CS_FLAG_BOUNDARY_CELL_BY_EDGE)) > 0)
    return true;
  else
    return false;
}

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

static inline short int
cs_cell_mesh_get_v(const cs_lnum_t              v_id,
                   const cs_cell_mesh_t  *const cm)
{
  if (cm == NULL)
    return -1;
  for (short int v = 0; v < cm->n_vc; v++)
    if (cm->v_ids[v] == v_id)
      return v;
  return -1;
}

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

static inline short int
cs_cell_mesh_get_e(const cs_lnum_t              e_id,
                   const cs_cell_mesh_t  *const cm)
{
  if (cm == NULL)
    return -1;
  for (short int e = 0; e < cm->n_ec; e++)
    if (cm->e_ids[e] == e_id)
      return e;
  return -1;
}

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

static inline short int
cs_cell_mesh_get_f(const cs_lnum_t              f_id,
                   const cs_cell_mesh_t  *const cm)
{
  if (cm == NULL)
    return -1;
  for (short int f = 0; f < cm->n_fc; f++)
    if (cm->f_ids[f] == f_id)
      return f;
  return -1;
}

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

static inline void
cs_cell_mesh_get_f2v(short int                    f,
                     const cs_cell_mesh_t        *cm,
                     short int                   *n_vf,
                     short int                   *v_ids)
{
  /* Reset */
  *n_vf = 0;
  for (short int v = 0; v < cm->n_vc; v++) v_ids[v] = -1;

  /* Tag vertices belonging to the current face f */
  for (short int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {

    const int  shift_e = 2*cm->f2e_ids[i];
    v_ids[cm->e2v_ids[shift_e]] = 1;
    v_ids[cm->e2v_ids[shift_e+1]] = 1;

  } /* Loop on face edges */

  for (short int v = 0; v < cm->n_vc; v++) {
    if (v_ids[v] > 0)
      v_ids[*n_vf] = v, *n_vf += 1;
  }

}

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

static inline void
cs_cell_mesh_get_next_3_vertices(const short int   *f2e_ids,
                                 const short int   *e2v_ids,
                                 short int         *v0,
                                 short int         *v1,
                                 short int         *v2)
{
  const short int e0  = f2e_ids[0];
  const short int e1  = f2e_ids[1];
  const short int tmp = e2v_ids[2*e1];

  *v0 = e2v_ids[2*e0];
  *v1 = e2v_ids[2*e0+1];
  *v2 = ((tmp != *v0) && (tmp != *v1)) ? tmp : e2v_ids[2*e1+1];
}

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

static inline bool
cs_cell_mesh_is_boundary_face(const cs_cell_mesh_t    *cm,
                              const short int          f)
{
  if (cm->f_ids[f] - cm->bface_shift > -1)
    return true;
  else
    return false;
}

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

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

void
cs_cdo_local_initialize(const cs_cdo_connect_t     *connect);

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

void
cs_cdo_local_finalize(void);

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

cs_cell_sys_t *
cs_cell_sys_create(int    n_max_dofbyc,
                   int    n_max_fbyc,
                   int    n_blocks,
                   int   *block_sizes);

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

void
cs_cell_sys_reset(int              n_fbyc,
                  cs_cell_sys_t   *csys);

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

void
cs_cell_sys_free(cs_cell_sys_t     **p_csys);

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

void
cs_cell_sys_dump(const char              msg[],
                 const cs_cell_sys_t    *csys);

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

cs_cell_builder_t *
cs_cell_builder_create(void);

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

void
cs_cell_builder_free(cs_cell_builder_t     **p_cb);

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

cs_cell_mesh_t *
cs_cell_mesh_create(const cs_cdo_connect_t   *connect);

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

cs_cell_mesh_t *
cs_cdo_local_get_cell_mesh(int    mesh_id);

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

void
cs_cell_mesh_reset(cs_cell_mesh_t   *cm);

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

void
cs_cell_mesh_dump(const cs_cell_mesh_t     *cm);

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

void
cs_cell_mesh_free(cs_cell_mesh_t     **p_cm);

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

void
cs_cell_mesh_build(cs_lnum_t                    c_id,
                   cs_eflag_t                   build_flag,
                   const cs_cdo_connect_t      *connect,
                   const cs_cdo_quantities_t   *quant,
                   cs_cell_mesh_t              *cm);

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

cs_face_mesh_t *
cs_face_mesh_create(short int   n_max_vbyf);

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

cs_face_mesh_t *
cs_cdo_local_get_face_mesh(int    mesh_id);

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

void
cs_face_mesh_free(cs_face_mesh_t     **p_fm);

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

void
cs_face_mesh_build(cs_lnum_t                    c_id,
                   cs_lnum_t                    f_id,
                   const cs_cdo_connect_t      *connect,
                   const cs_cdo_quantities_t   *quant,
                   cs_face_mesh_t              *fm);

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

void
cs_face_mesh_build_from_cell_mesh(const cs_cell_mesh_t    *cm,
                                  short int                f,
                                  cs_face_mesh_t          *fm);

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

cs_face_mesh_light_t *
cs_face_mesh_light_create(short int   n_max_vbyf,
                          short int   n_max_vbyc);

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

cs_face_mesh_light_t *
cs_cdo_local_get_face_mesh_light(int    mesh_id);

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

void
cs_face_mesh_light_free(cs_face_mesh_light_t     **p_fm);

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

void
cs_face_mesh_light_build(const cs_cell_mesh_t    *cm,
                         short int                f,
                         cs_face_mesh_light_t    *fm);

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

END_C_DECLS

#endif /* __CS_CDO_LOCAL_H__ */