/*============================================================================ * * 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 #include #include #include #include /*---------------------------------------------------------------------------- * Fichiers `include' librairie standard C ou BFT *----------------------------------------------------------------------------*/ #include #include /*---------------------------------------------------------------------------- * 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 */