/*============================================================================
*
* Code_Saturne version 1.3
* ------------------------
*
*
* This file is part of the Code_Saturne Kernel, element of the
* Code_Saturne CFD tool.
*
* Copyright (C) 1998-2007 EDF S.A., France
*
* contact: saturne-support@edf.fr
*
* The Code_Saturne Kernel 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.
*
* The Code_Saturne Kernel 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 the Code_Saturne Kernel; if not, write to the
* Free Software Foundation, Inc.,
* 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*
*============================================================================*/
/*============================================================================
* Grandeurs associées à un maillage
*============================================================================*/
/* includes système */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
/*----------------------------------------------------------------------------
* Fichiers `include' librairie standard C ou BFT
*----------------------------------------------------------------------------*/
#include <bft_mem.h>
#include <bft_error.h>
/*----------------------------------------------------------------------------
* Fichiers `include' locaux
*----------------------------------------------------------------------------*/
#include "cs_base.h"
#include "cs_maillage.h"
#include "cs_maillage_connect.h"
#include "cs_parallel.h"
#include "cs_perio.h"
/*----------------------------------------------------------------------------
* Fichiers `include' associés au fichier courant
*----------------------------------------------------------------------------*/
#include "cs_maillage_grd.h"
#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */
/*=============================================================================
* Définitions de macros locales
*============================================================================*/
enum {X, Y, Z} ;
#define _CS_PRODUIT_VECTORIEL(prod_vect, vect1, vect2) \
(prod_vect[X] = vect1[Y] * vect2[Z] - vect2[Y] * vect1[Z], \
prod_vect[Y] = vect2[X] * vect1[Z] - vect1[X] * vect2[Z], \
prod_vect[Z] = vect1[X] * vect2[Y] - vect2[X] * vect1[Y])
#define _CS_PRODUIT_SCALAIRE(vect1, vect2) \
(vect1[X] * vect2[X] + vect1[Y] * vect2[Y] + vect1[Z] * vect2[Z])
#define _CS_MODULE(vect) \
sqrt(vect[X] * vect[X] + vect[Y] * vect[Y] + vect[Z] * vect[Z])
/*============================================================================
* Variables globales statiques
*============================================================================*/
/* Pointeur sur les grandeurs du maillage principal */
cs_maillage_grd_t *cs_glob_maillage_grd = NULL;
/* Choix de l'algorithme de calcul des centres de gravité des cellules */
static int cs_glob_maillage_grd_cdg_cel = 0;
/*============================================================================
* Prototypes de fonctions privées
*============================================================================*/
/*----------------------------------------------------------------------------
* Calcul des grandeurs associées à des faces (intérieures ou de bord)
*
* Pi+1
* *---------* B : barycentre du polygone P
* / . . \
* / . . \ Pi : sommet du polygone P
* / . . \
* / . . Ti \ Ti : triangle
* *.........B.........* Pi
* Pn-1 \ . . /
* \ . . /
* \ . . /
* \ . T0 . /
* *---------*
* P0
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_calc_faces
(
const cs_int_t dim, /* --> Dimension de l'espace */
const cs_int_t nbr_fac, /* --> Nombre de faces */
const cs_real_t coo_som[], /* --> Coordonnées des sommets */
const cs_int_t pos_fac_som[], /* --> Positions de départ dans
connect. faces -> sommets */
const cs_int_t val_fac_som[], /* --> Connect. faces -> sommets */
cs_real_t cdg_fac[], /* <-- CDG faces */
cs_real_t surf_fac[] /* <-- Normales faces */
);
/*----------------------------------------------------------------------------
* Calcul des centres des cellules C *
* *
* à partir de leurs n sommets S(i) où i=0,n-1 *
* *
* *
* --> 1 n-1 --> *
* OB(C) = --- Somme OSi *
* n i=0 *
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_cdg_cel_som
(
const cs_maillage_t *const maillage, /* --> Maillage associé */
cs_real_t cen_cel[] /* <-- Barycentres des cellules */
);
/*----------------------------------------------------------------------------
* Calcul des centres G des cellules C *
* *
* à partir de leurs n faces F(i) où i=0,n-1 *
* *
* *
* n-1 *
* Somme Surf(Fi) G(Fi) *
* i=0 *
* G(C) = ----------------------- *
* n-1 *
* Somme Surf(Fi) *
* i=0 *
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_cdg_cel_fac
(
const cs_maillage_t *const maillage, /* --> Maillage associé */
const cs_real_t surf_fac[], /* --> Surfaces des faces internes */
const cs_real_t cdg_fac[], /* --> Centres de gravité des
faces internes */
const cs_real_t surf_fbr[], /* --> Surfaces des faces de bord */
const cs_real_t cdg_fbr[], /* --> Centres de gravité des
faces de bord */
cs_real_t cen_cel[] /* <-- Centres de gravité des
cellules */
);
/*----------------------------------------------------------------------------
* Calcul du volume des cellules C à partir de leurs n faces F(i) *
* et de leur centre de gravité G(Fi) où i=0,n-1 *
* *
* 1 n-1 *
* G(C) = - . Somme Surf(Fi) G(Fi) *
* 3 i = 0 *
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_vol_cel
(
const cs_maillage_t *const maillage, /* --> Maillage associé */
const cs_real_t surf_fac[], /* --> Surfaces des faces internes */
const cs_real_t cdg_fac[], /* --> Centres de gravité des
faces internes */
const cs_real_t surf_fbr[], /* --> Surfaces des faces de bord */
const cs_real_t cdg_fbr[], /* --> Centres de gravité des
faces de bord */
cs_real_t vol_cel[] /* <-- Volume des cellules */
);
/*============================================================================
* Fonctions publiques pour API Fortran
*============================================================================*/
/*----------------------------------------------------------------------------
* Interrogation ou modification du choix de l'algorithme de calcul
* des centres de gravité des cellules ;
* Au retour, IOPT contient la valeur (1 à 2) correspondant à l'algorithme
* choisi
*
* Interface Fortran :
*
* SUBROUTINE ALGCEN (IOPT)
* *****************
*
* INTEGER IOPT : <-> : Choix de l'algorithme
* < 0 : interrogation
* 0 : calcul basé sur les centres et
* surfaces des faces (par défaut)
* 1 : calcul basé sur les sommets
*----------------------------------------------------------------------------*/
void CS_PROCF (algcen, ALGCEN)
(
cs_int_t *const iopt
)
{
int iopt_ret = cs_maillage_grd_choix_cdg_cel((int)(*iopt));
*iopt = iopt_ret;
}
/*============================================================================
* Fonctions publiques
*============================================================================*/
/*----------------------------------------------------------------------------
* Interrogation ou modification du choix de l'algorithme de calcul
* des centres de gravité des cellules ;
* < 0 : interrogation
* 0 : calcul basé sur les centres et surfaces des faces (choix par défaut)
* 1 : calcul basé sur les sommets
* Renvoie la valeur (1 à 2) correspondant à l'algorithme choisi
*----------------------------------------------------------------------------*/
int cs_maillage_grd_choix_cdg_cel
(
const int choix_algo
)
{
if (choix_algo > 1)
bft_error(__FILE__, __LINE__,0,
_("L'indicateur de choix d'algorithme de calcul du centre "
"de gravité\n"
"des cellules peut prendre les valeurs suivantes :\n"
" 0 : calcul basé sur les centres et surfaces des faces\n"
" 1 : calcul basé sur les sommets\n"
"et non %d."), cs_glob_maillage_grd_cdg_cel);
else if (choix_algo >= 0)
cs_glob_maillage_grd_cdg_cel = choix_algo;
return cs_glob_maillage_grd_cdg_cel;
}
/*----------------------------------------------------------------------------
* Création d'une structure grandeurs maillage
*----------------------------------------------------------------------------*/
cs_maillage_grd_t * cs_maillage_grd_cree
(
void
)
{
cs_maillage_grd_t * maillage_grd;
BFT_MALLOC (maillage_grd, 1, cs_maillage_grd_t);
maillage_grd->cen_cel = NULL;
maillage_grd->vol_cel = NULL;
maillage_grd->surf_fac = NULL;
maillage_grd->surf_fbr = NULL;
maillage_grd->cdg_fac = NULL;
maillage_grd->cdg_fbr = NULL;
return (maillage_grd);
}
/*----------------------------------------------------------------------------
* Destruction d'une structure grandeurs maillage
*----------------------------------------------------------------------------*/
cs_maillage_grd_t * cs_maillage_grd_detruit
(
cs_maillage_grd_t * maillage_grd
)
{
BFT_FREE(maillage_grd->cen_cel);
BFT_FREE(maillage_grd->vol_cel);
BFT_FREE(maillage_grd->surf_fac);
BFT_FREE(maillage_grd->surf_fbr);
BFT_FREE(maillage_grd->cdg_fac);
BFT_FREE(maillage_grd->cdg_fbr);
BFT_FREE(maillage_grd);
return (maillage_grd);
}
/*----------------------------------------------------------------------------
* Calcul des grandeurs associées à un maillage
*----------------------------------------------------------------------------*/
void cs_maillage_grd_calc
(
const cs_maillage_t *const maillage, /* --> maillage utilisé */
cs_maillage_grd_t *const maillage_grd /* <-> grandeurs associées */
)
{
cs_int_t dim = maillage->dim;
cs_int_t nbr_fac = maillage->nbr_fac;
cs_int_t nbr_fbr = maillage->nbr_fbr;
cs_int_t nbr_cel = maillage->nbr_cel;
cs_int_t nbr_cel_etendu = maillage->nbr_cel_etendu;
/* Allocation des structures s'il ne s'agit pas d'une mise à jour */
if (maillage_grd->surf_fac == NULL)
BFT_MALLOC(maillage_grd->surf_fac, nbr_fac*dim, cs_real_t);
if (maillage_grd->cdg_fac == NULL)
BFT_MALLOC(maillage_grd->cdg_fac, nbr_fac*dim, cs_real_t);
if (maillage_grd->surf_fbr == NULL)
BFT_MALLOC(maillage_grd->surf_fbr, nbr_fbr*dim, cs_real_t);
if (maillage_grd->cdg_fbr == NULL)
BFT_MALLOC(maillage_grd->cdg_fbr, nbr_fbr*dim, cs_real_t);
if (maillage_grd->cen_cel == NULL)
BFT_MALLOC(maillage_grd->cen_cel, nbr_cel_etendu*dim, cs_real_t);
if (maillage_grd->vol_cel == NULL)
BFT_MALLOC(maillage_grd->vol_cel, nbr_cel_etendu, cs_real_t);
/* Calcul des centres de gravité et des surfaces des faces internes */
_cs_maillage_grd_calc_faces(dim,
nbr_fac,
maillage->coo_som,
maillage->pos_fac_som,
maillage->val_fac_som,
maillage_grd->cdg_fac,
maillage_grd->surf_fac);
/* Calcul des centres de gravité et des surfaces des faces de bord */
_cs_maillage_grd_calc_faces(dim,
nbr_fbr,
maillage->coo_som,
maillage->pos_fbr_som,
maillage->val_fbr_som,
maillage_grd->cdg_fbr,
maillage_grd->surf_fbr);
/* Calcul des centres des cellules basé sur le barycentre des faces
ou sur les sommets */
switch (cs_glob_maillage_grd_cdg_cel) {
case 0:
_cs_maillage_grd_cdg_cel_fac(maillage,
maillage_grd->surf_fac,
maillage_grd->cdg_fac,
maillage_grd->surf_fbr,
maillage_grd->cdg_fbr,
maillage_grd->cen_cel);
break;
case 1:
_cs_maillage_grd_cdg_cel_som(maillage,
maillage_grd->cen_cel);
break;
default:
assert(0);
}
/* Calcul du volume des cellules */
_cs_maillage_grd_vol_cel(maillage,
maillage_grd->surf_fac,
maillage_grd->cdg_fac,
maillage_grd->surf_fbr,
maillage_grd->cdg_fbr,
maillage_grd->vol_cel);
if (cs_glob_base_nbr > 1) {
cs_int_t idim, icel;
cs_real_t *coord_tmp = NULL;
/* Synchronisation des coordonnées des centres des cellules
entre processeurs dans les zones halo */
BFT_MALLOC(coord_tmp, nbr_cel_etendu, cs_real_t);
for (idim = 0 ; idim < dim ; idim++) {
for (icel = 0 ; icel < nbr_cel ; icel++)
coord_tmp[icel] = maillage_grd->cen_cel[icel*dim + idim];
for (icel = nbr_cel ; icel < nbr_cel_etendu ; icel++)
coord_tmp[icel] = 0.;
CS_PROCF (parcom, PARCOM) (coord_tmp);
for (icel = nbr_cel ; icel < nbr_cel_etendu ; icel++)
maillage_grd->cen_cel[icel*dim + idim] = coord_tmp[icel];
}
BFT_FREE(coord_tmp);
CS_PROCF (pargve, PARGVE) (&dim, maillage_grd->cen_cel);
/* Synchronisation du volume des cellules entre processeurs
dans les zones halo */
for (icel = nbr_cel ; icel < nbr_cel_etendu ; icel++)
maillage_grd->vol_cel[icel] = 0.;
CS_PROCF (parcom, PARCOM) (maillage_grd->vol_cel);
}
/* En cas de halo monoprocesseur (périodicité), mise à zéro du halo */
else if (nbr_cel_etendu > nbr_cel) {
cs_int_t idim, icel;
for (icel = nbr_cel ; icel < nbr_cel_etendu ; icel++) {
/* Traitement du volume des cellules */
maillage_grd->vol_cel[icel] = 0.;
/* Traitement des centres cellule */
for (idim = 0 ; idim < dim ; idim++)
maillage_grd->cen_cel[icel*dim + idim] = 0.;
}
}
/* En cas de périodicité, mise à jour du halo */
if (maillage->nbr_per > 0)
cs_perio_sync_geo() ;
}
/*============================================================================
* Fonctions privées
*============================================================================*/
/*----------------------------------------------------------------------------
* Calcul des grandeurs associées à des faces (intérieures ou de bord)
*
* Pi+1
* *---------* B : barycentre du polygone P
* / . . \
* / . . \ Pi : sommet du polygone P
* / . . \
* / . . Ti \ Ti : triangle
* *.........B.........* Pi
* Pn-1 \ . . /
* \ . . /
* \ . . /
* \ . T0 . /
* *---------*
* P0
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_calc_faces
(
const cs_int_t dim, /* --> Dimension de l'espace */
const cs_int_t nbr_fac, /* --> Nombre de faces */
const cs_real_t coo_som[], /* --> Coordonnées des sommets */
const cs_int_t pos_fac_som[], /* --> Positions de départ dans
connect. faces -> sommets */
const cs_int_t val_fac_som[], /* --> Connect. faces -> sommets */
cs_real_t cdg_fac[], /* <-- CDG faces */
cs_real_t surf_fac[] /* <-- Normales faces */
)
{
cs_int_t ifac;
cs_int_t nbr_som_fac, nbr_max_som_fac;
cs_int_t icoo, icoo_inf;
cs_int_t isom, isom_inf, isom_sup;
cs_int_t itri;
cs_point_t *coo_som_fac;
cs_point_t barycentre_fac;
cs_point_t normale_fac;
cs_point_t centre_fac;
cs_point_t centre_tri;
cs_point_t *normale_tri;
cs_point_t vect1, vect2;
cs_real_t surface_fac;
cs_real_t surface_tri;
cs_real_t contrib_vol_fac;
cs_real_t contrib_vol_tri;
cs_real_t correc_cdg;
/* Retour si l'on n'a pas les informations suffisantes
(cas SolCom hors radiatif ou Enveloppe pré 1.2.d sans option "-n") */
if (pos_fac_som == NULL && val_fac_som == NULL)
return;
/* Vérifications */
/*---------------*/
if (dim != 3)
bft_error(__FILE__, __LINE__,0,
_("Le calcul des grandeurs géométriques associées aux faces\n"
"ne sont implantés qu'en 3D."));
assert(cdg_fac != NULL || nbr_fac == 0);
assert(surf_fac != NULL || nbr_fac == 0);
/* Comptages et allocations */
/*--------------------------*/
nbr_max_som_fac = 0;
for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
nbr_som_fac = pos_fac_som[ifac + 1] - pos_fac_som[ifac];
if (nbr_max_som_fac <= nbr_som_fac) nbr_max_som_fac = nbr_som_fac;
}
BFT_MALLOC(coo_som_fac, nbr_max_som_fac + 1, cs_point_t);
BFT_MALLOC(normale_tri, nbr_max_som_fac , cs_point_t);
/*=========================================================================*/
/* Boucle sur les faces */
/*=========================================================================*/
for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
contrib_vol_tri = 0. ;
surface_fac = 0.0 ;
/*------------------------------------------------------------------------*/
/* Définition du polygone P en fonction des sommets Pi de la face */
/*------------------------------------------------------------------------*/
isom_inf = pos_fac_som[ifac] - 1;
isom_sup = pos_fac_som[ifac + 1] - 1;
nbr_som_fac = 0;
for (isom = isom_inf ; isom < isom_sup ; isom++) {
icoo_inf = 3 * (val_fac_som[isom] - 1);
for (icoo = X ; icoo < 3 ; icoo++)
coo_som_fac[nbr_som_fac][icoo] = coo_som[icoo_inf + icoo];
nbr_som_fac++;
}
for (icoo = X ; icoo < 3 ; icoo++)
coo_som_fac[nbr_som_fac][icoo] = coo_som_fac[0][icoo];
/*------------------------------------------------------------------------*/
/* Calcul des coordonnees du barycentre B du polygone P */
/* */
/* --> 1 n-1 --> */
/* OB = - Somme OPi */
/* n i=0 */
/*------------------------------------------------------------------------*/
for (icoo = X ; icoo < 3 ; icoo++) {
barycentre_fac[icoo] = 0.0;
for (isom = 0 ; isom < nbr_som_fac ; isom++)
barycentre_fac[icoo] += coo_som_fac[isom][icoo];
barycentre_fac[icoo] /= nbr_som_fac ;
}
for (icoo = X ; icoo < 3 ; icoo++) {
normale_fac[icoo] = 0.0;
centre_fac [icoo] = 0.0;
}
/* Première boucle sur les triangles de la face (calcul des normales) */
/*========================================================================*/
for (itri = 0 ; itri < nbr_som_fac ; itri++) {
/*----------------------------------------------------------------------*/
/* Calcul de la normale au triangle Ti : */
/* */
/* -> --> --> */
/* N(Ti) = 1/2 ( BPi X BPi+1 ) */
/*----------------------------------------------------------------------*/
for (icoo = X ; icoo < 3 ; icoo++) {
vect1[icoo] = coo_som_fac[itri ][icoo] - barycentre_fac[icoo];
vect2[icoo] = coo_som_fac[itri + 1][icoo] - barycentre_fac[icoo];
}
_CS_PRODUIT_VECTORIEL(normale_tri[itri], vect1, vect2);
for (icoo = X ; icoo < 3 ; icoo++)
normale_tri[itri][icoo] *= 0.5;
/*----------------------------------------------------------------------*/
/* Calcul de la normale au polygone */
/* comme somme vectorielle des normales aux triangles : */
/* */
/* -> n-1 -> */
/* N(P) = Somme( N(Ti) ) */
/* i=0 */
/*----------------------------------------------------------------------*/
for (icoo = X ; icoo < 3 ; icoo++)
normale_fac[icoo] += normale_tri[itri][icoo];
} /* Fin de la premiere boucle sur les triangles de la face */
/* Seconde boucle sur les triangles de la face (pour le barycentre) */
/*==================================================================*/
for (itri = 0 ; itri < nbr_som_fac ; itri++) {
/*----------------------------------------------------------------------*/
/* Calcul du centre de gravite G(Ti) du triangle Ti : */
/* */
/* --> --> --> --> */
/* OG(Ti) = 1/3 ( OB + OPi + OPi+1 ) */
/* */
/* Et de la contribution au volume du triangle Ti : */
/* */
/* --> -> */
/* OG(Ti).N(Ti) */
/*----------------------------------------------------------------------*/
for (icoo = X ; icoo < 3 ; icoo++) {
centre_tri[icoo] = barycentre_fac[icoo]
+ coo_som_fac[itri ][icoo]
+ coo_som_fac[itri + 1][icoo];
centre_tri[icoo] /= 3.0;
contrib_vol_tri += (centre_tri[icoo] * normale_tri[itri][icoo]);
}
/*----------------------------------------------------------------------*/
/* Calcul de la surface du triangle Ti comme module de la normale */
/* */
/* -> */
/* Surf(Ti) = | N(Ti) | */
/*----------------------------------------------------------------------*/
surface_tri = _CS_MODULE(normale_tri[itri]);
if (_CS_PRODUIT_SCALAIRE(normale_tri[itri], normale_fac) < 0.0)
surface_tri *= -1.0;
surface_fac += surface_tri;
/*----------------------------------------------------------------------*/
/* n-1 */
/* Somme Surf(Ti) G(Ti) */
/* i=0 */
/*----------------------------------------------------------------------*/
for (icoo = X ; icoo < 3 ; icoo++)
centre_fac[icoo] += surface_tri * centre_tri[icoo];
} /* Fin de la seconde boucle sur les triangles de la face */
/*------------------------------------------------------------------------*/
/* Calcul du centre de gravite G(P) du polygone P : */
/* */
/* n-1 */
/* Somme Surf(Ti) G(Ti) */
/* i=0 */
/* G(P) = ----------------------- */
/* n-1 */
/* Somme Surf(Ti) */
/* i=0 */
/* */
/* Calcul de la contribution au volume du polygone (avant correction) : */
/* */
/* --> -> */
/* OG(P).N(P) */
/*------------------------------------------------------------------------*/
contrib_vol_fac = 0.0;
for (icoo = X ; icoo < 3 ; icoo++) {
centre_fac[icoo] = centre_fac[icoo] / surface_fac;
contrib_vol_fac += (centre_fac[icoo] * normale_fac[icoo]);
}
correc_cdg = (contrib_vol_tri - contrib_vol_fac)
/ (surface_fac * surface_fac);
for (icoo = X ; icoo < 3 ; icoo++)
centre_fac[icoo] += correc_cdg * normale_fac[icoo];
/* Stockage du résultat dans les structures appropriées */
/*------------------------------------------------------*/
for (icoo = X ; icoo < 3 ; icoo++) {
cdg_fac [ifac*3 + icoo] = centre_fac[icoo];
surf_fac[ifac*3 + icoo] = normale_fac[icoo];
}
} /* Fin : boucle sur les faces */
BFT_FREE(normale_tri);
BFT_FREE(coo_som_fac);
}
/*----------------------------------------------------------------------------
* Calcul des centres G des cellules C *
* *
* à partir de leurs n faces F(i) où i=0,n-1 *
* *
* *
* n-1 *
* Somme Surf(Fi) G(Fi) *
* i=0 *
* G(C) = ----------------------- *
* n-1 *
* Somme Surf(Fi) *
* i=0 *
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_cdg_cel_fac
(
const cs_maillage_t *const maillage, /* --> Maillage associé */
const cs_real_t surf_fac[], /* --> Surfaces des faces internes */
const cs_real_t cdg_fac[], /* --> Centres de gravité des
faces internes */
const cs_real_t surf_fbr[], /* --> Surfaces des faces de bord */
const cs_real_t cdg_fbr[], /* --> Centres de gravité des
faces de bord */
cs_real_t cen_cel[] /* <-- Centres de gravité des
cellules */
)
{
cs_int_t ifac, icoo, ival;
cs_int_t icel, icel1, icel2;
cs_real_t *surf_tot;
cs_real_t surf;
cs_real_t vect_loc_surf[3];
/* Connectivité du maillage */
const cs_int_t dim = maillage->dim;
const cs_int_t nbr_fac = maillage->nbr_fac;
const cs_int_t nbr_fbr = maillage->nbr_fbr;
const cs_int_t nbr_cel = maillage->nbr_cel;
const cs_int_t nbr_cel_etendu = maillage->nbr_cel_etendu;
const cs_int_t *fac_cel = maillage->fac_cel;
const cs_int_t *fbr_cel = maillage->fbr_cel;
/* Retour si l'on n'a pas les informations suffisantes
(cas SolCom hors radiatif ou Enveloppe pré 1.2.d sans option "-n") */
if (maillage->val_fac_som == NULL && maillage->val_fbr_som == NULL)
return;
/* Vérifications */
/* ------------- */
if (dim != 3)
bft_error(__FILE__, __LINE__,0,
_("Le calcul des centres des cellules\n"
"ne sont implantés qu'en 3D."));
assert(cen_cel != NULL);
/* Initialisation */
/* -------------- */
BFT_MALLOC(surf_tot, nbr_cel_etendu, cs_real_t);
for (ival = 0 ; ival < nbr_cel_etendu ; ival++) {
surf_tot[ival] = 0.;
for (icoo = 0; icoo < dim ; icoo++) {
cen_cel[dim*ival + icoo] = 0.;
}
}
/* ----------------------------- */
/* Boucle sur les faces internes */
/* ----------------------------- */
for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
/* ----------------------------------------------------------- */
/* Pour chacune des deux cellules voisines à la face interne, */
/* on met à jour le numérateur de cen_cel et surf_tot */
/* ----------------------------------------------------------- */
icel1 = fac_cel[2*ifac ] - 1;
icel2 = fac_cel[2*ifac + 1] - 1;
/* Calcul de la surface de la face */
for (icoo = 0 ; icoo < dim ; icoo++)
vect_loc_surf[icoo] = surf_fac[dim*ifac + icoo];
surf = _CS_MODULE(vect_loc_surf);
surf_tot[icel1] += surf;
surf_tot[icel2] += surf;
/* Calcul du numérateur */
for (icoo = 0 ; icoo < dim ; icoo++) {
cen_cel[dim*icel1 + icoo] += cdg_fac[dim*ifac + icoo]*surf;
cen_cel[dim*icel2 + icoo] += cdg_fac[dim*ifac + icoo]*surf;
}
} /* Fin de la boucle sur les faces internes */
/* ----------------------------- */
/* Boucle sur les faces de bord */
/* ----------------------------- */
for (ifac = 0 ; ifac < nbr_fbr ; ifac++) {
/* ----------------------------------------------------------- */
/* Pour chaque cellule voisine de la face bord, */
/* on met à jour le numérateur de cen_cel et surf_tot */
/* ----------------------------------------------------------- */
icel1 = fbr_cel[ifac] - 1;
/* Calcul de la surface de la face */
for (icoo = 0 ; icoo < dim ; icoo++)
vect_loc_surf[icoo] = surf_fbr[dim*ifac + icoo];
surf = _CS_MODULE(vect_loc_surf);
surf_tot[icel1] += surf;
/* Calcul du numérateur */
for (icoo = 0 ; icoo < dim ; icoo++)
cen_cel[dim*icel1 + icoo] += cdg_fbr[dim*ifac + icoo]*surf;
} /* Fin de la boucle sur les faces de bord */
/* ------------------------------------------------------------------*/
/* Boucle sur les cellules pour finaliser la détermination du centre */
/* ------------------------------------------------------------------*/
for (icel = 0 ; icel < nbr_cel ; icel++) {
for (icoo = 0 ; icoo < dim ; icoo++ )
cen_cel[icel*dim + icoo] /= surf_tot[icel];
} /* Fin de la boucle sur les cellules */
/* Libération de la mémoire des pointeurs */
BFT_FREE(surf_tot);
}
/*----------------------------------------------------------------------------
* Calcul des centres des cellules C *
* *
* à partir de leurs n sommets S(i) où i=0,n-1 *
* *
* *
* --> 1 n-1 --> *
* OB(C) = --- Somme OSi *
* n i=0 *
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_cdg_cel_som
(
const cs_maillage_t *const maillage, /* --> Maillage associé */
cs_real_t cen_cel[] /* <-- Barycentres des cellules */
)
{
cs_int_t icel, ifac, isom, icoo;
cs_int_t ipos_fac, ipos_som;
cs_int_t num_fac;
cs_int_t cpt_som_cel;
cs_int_t *ind_cel_som;
cs_int_t *pos_fac_som_loc;
cs_int_t *val_fac_som_loc;
cs_int_t *pos_cel_fac = NULL;
cs_int_t *val_cel_fac = NULL;
/* Retour si l'on n'a pas les informations suffisantes
(cas SolCom hors radiatif ou Enveloppe pré 1.2.d sans option "-n") */
if (maillage->val_fac_som == NULL && maillage->val_fbr_som == NULL)
return;
/* Vérifications */
/* ------------- */
if (maillage->dim != 3)
bft_error(__FILE__, __LINE__,0,
_("Le calcul des centres des cellules\n"
"ne sont implantés qu'en 3D."));
assert(cen_cel != NULL);
/* Allocations et Initialisation */
/* ----------------------------- */
BFT_MALLOC(ind_cel_som, maillage->nbr_som, cs_int_t);
for (isom = 0 ; isom < maillage->nbr_som ; isom++)
ind_cel_som[isom] = -1;
/* Extraction de la connectivité "cellules -> faces" */
cs_maillage_ret_cel_fac(maillage,
maillage->nbr_cel,
NULL,
&pos_cel_fac,
&val_cel_fac);
/* Boucle sur les cellules */
/* ----------------------- */
for (icel = 0 ; icel < maillage->nbr_cel ; icel++) {
cpt_som_cel = 0;
for (icoo = 0 ; icoo < 3 ; icoo++)
cen_cel[icel*3 + icoo] = 0.;
/* Boucle sur les faces de la cellule */
for (ipos_fac = pos_cel_fac[icel] ;
ipos_fac < pos_cel_fac[icel + 1] ;
ipos_fac++) {
num_fac = val_cel_fac[ipos_fac - 1];
/* Face interne ou face de bord */
if (num_fac > maillage->nbr_fbr) {
ifac = num_fac - maillage->nbr_fbr - 1;
pos_fac_som_loc = maillage->pos_fac_som;
val_fac_som_loc = maillage->val_fac_som;
}
else {
ifac = num_fac - 1;
pos_fac_som_loc = maillage->pos_fbr_som;
val_fac_som_loc = maillage->val_fbr_som;
}
/* Boucle sur les sommets de la face */
for (ipos_som = pos_fac_som_loc[ifac] ;
ipos_som < pos_fac_som_loc[ifac + 1] ;
ipos_som++) {
isom = val_fac_som_loc[ipos_som] - 1;
if (ind_cel_som[isom] < icel) {
for (icoo = 0 ; icoo < 3 ; icoo++)
cen_cel[icel*3 + icoo] += maillage->coo_som[isom*3 + icoo];
cpt_som_cel += 1;
ind_cel_som[isom] = icel;
}
}
} /* Fin de la boucle sur les faces de la cellule */
for (icoo = 0 ; icoo < 3 ; icoo++)
cen_cel[icel*3 + icoo] /= (double)cpt_som_cel;
} /* Fin de la boucle sur les cellules */
/* Libération de la mémoire allouée */
BFT_FREE(ind_cel_som);
BFT_FREE(pos_cel_fac);
BFT_FREE(val_cel_fac);
}
/*----------------------------------------------------------------------------
* Calcul du volume des cellules C à partir de leurs n faces F(i) *
* et de leur centre de gravité G(Fi) où i=0,n-1 *
* *
* 1 n-1 *
* G(C) = - . Somme Surf(Fi) G(Fi) *
* 3 i = 0 *
*----------------------------------------------------------------------------*/
static void _cs_maillage_grd_vol_cel
(
const cs_maillage_t *const maillage, /* --> Maillage associé */
const cs_real_t surf_fac[], /* --> Surfaces des faces internes */
const cs_real_t cdg_fac[], /* --> Centres de gravité des
faces internes */
const cs_real_t surf_fbr[], /* --> Surfaces des faces de bord */
const cs_real_t cdg_fbr[], /* --> Centres de gravité des
faces de bord */
cs_real_t vol_cel[] /* <-- Volume des cellules */
)
{
cs_int_t i_elt1, i_elt2;
cs_int_t i_fac, i_fbr, i_cel;
cs_int_t i_dim;
cs_real_t flux = 0;
const cs_real_t un_tiers = 1.0/3.0;
const cs_int_t dim = maillage->dim;
/* Initialisation */
/*----------------*/
for (i_cel = 0; i_cel < maillage->nbr_cel_etendu; i_cel++)
vol_cel[i_cel] = 0;
/* Boucle sur les faces internes */
/* ----------------------------- */
for (i_fac = 0; i_fac < maillage->nbr_fac; i_fac++) {
i_elt1 = maillage->fac_cel[2*i_fac] - 1;
i_elt2 = maillage->fac_cel[2*i_fac + 1] - 1;
flux = 0;
for (i_dim = 0; i_dim < dim; i_dim++)
flux += surf_fac[dim*i_fac + i_dim] * cdg_fac[dim*i_fac + i_dim];
vol_cel[i_elt1] += flux;
vol_cel[i_elt2] -= flux;
}
/* Boucle sur les faces de bord */
/* ---------------------------- */
for (i_fbr = 0; i_fbr < maillage->nbr_fbr; i_fbr++) {
i_elt1 = maillage->fbr_cel[i_fbr] - 1;
flux = 0;
for (i_dim = 0; i_dim < dim; i_dim++)
flux += surf_fbr[dim*i_fbr + i_dim] * cdg_fbr[dim*i_fbr + i_dim];
vol_cel[i_elt1] += flux;
}
/* Calcul du volume */
for (i_cel = 0; i_cel < maillage->nbr_cel; i_cel++)
vol_cel[i_cel] *= un_tiers;
}
/*----------------------------------------------------------------------------*/
#if 0 /* Morceau de test bonne orientation faces */
cs_int_t idim, ifac, icel;
cs_real_t cdgfac[3];
cs_real_t cdgcel[3];
cs_real_t normale[3];
cs_real_t pscal;
for (ifac = 0 ; ifac < maillage->nbr_fbr ; ifac++) {
icel = maillage->fbr_cel[ifac] - 1;
pscal = 0;
for (idim = 0 ; idim < 3 ; idim++) {
cdgcel[idim] = cs_glob_maillage_grd->cen_cel[icel*3 + idim];
cdgfac[idim] = cs_glob_maillage_grd->cdg_fbr[ifac*3 + idim];
normale[idim] = cs_glob_maillage_grd->surf_fbr[ifac*3 + idim];
pscal += normale[idim] * (cdgfac[idim] - cdgcel[idim]);
}
if (pscal < 0.0)
printf("num_fac_brd = %d, num_cel = %d, pscal = %f\n",
ifac + 1, icel + 1, pscal);
}
#endif
/*----------------------------------------------------------------------------*/
/* Suppression des macros locales */
#undef _CS_PRODUIT_VECTORIEL
#undef _CS_PRODUIT_SCALAIRE
#undef _CS_MODULE
#ifdef __cplusplus
}
#endif /* __cplusplus */
syntax highlighted by Code2HTML, v. 0.9.1