/*============================================================================
*
* 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
*
*============================================================================*/
/*============================================================================
* Définitions, variables globales, et fonctions de base
*============================================================================*/
/* includes système */
#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
/* Includes librairie BFT et FVM */
#include <bft_error.h>
#include <bft_mem.h>
#include <bft_printf.h>
/* Includes librairie */
#include "cs_base.h"
#include "cs_maillage.h"
#include "cs_parallel.h"
#include "cs_voiset.h"
#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */
/*============================================================================
* Variables statiques locales
*============================================================================*/
/*============================================================================
* Prototypes de fonctions privées
*============================================================================*/
/*============================================================================
* Prototypes de fonctions FORTRAN appelées
*============================================================================*/
/*----------------------------------------------------------------------------
* Calcul des gradients par moindres carrés (standard ou avec voisinage étendu)
*----------------------------------------------------------------------------*/
void CS_PROCF (gradmc, GRADMC)
(
const cs_int_t *const idbia0, /* --> numéro 1ère case libre dans IA */
const cs_int_t *const idbra0, /* --> numéro 1ère case libre dans RA */
const cs_int_t *const ncelet, /* --> nombre de cellules étendu */
const cs_int_t *const ncel, /* --> nombre de cellules */
const cs_int_t *const nfac, /* --> nombre de faces internes */
const cs_int_t *const nfabor, /* --> nombre de faces de bord */
const cs_int_t *const ncelbr, /* --> nombre de cellules de bord */
const cs_int_t *const nituse, /* --> longueur du tableau ituser[] */
const cs_int_t *const nrtuse, /* --> longueur du tableau rtuser[] */
const cs_int_t *const inc, /* --> 0 ou 1 : incrément ou non */
const cs_int_t *const iccocg, /* --> 1 ou 0 : recalcul de COCG ou non */
const cs_int_t *const nswrgp, /* --> >1 : avec reconstruction */
const cs_int_t *const idimte, /* --> 0, 1, 2 : variable
scalaire, vecteur, tenseur */
const cs_int_t *const itenso, /* --> pour la périodicité de rotation */
const cs_int_t *const iphydp, /* --> prise en compte de P hydrostatiqu*/
const cs_int_t *const imrgra, /* --> mode de calcul du gradient */
const cs_int_t *const iwarnp, /* --> niveau d'impression */
const cs_int_t *const nfecra, /* --> unité sortie standard */
const cs_real_t *const epsrgp, /* --> précision pour calcul du gradient
itératif */
const cs_real_t *const extrap, /* --> extrapolation des grad. au bord */
const cs_int_t ifacel[], /* --> liste des faces internes */
const cs_int_t ifabor[], /* --> liste des faces de bord */
const cs_int_t icelbr[], /* --> numéro des cellules ayant au
moins une face de bord */
const cs_int_t ipcvse[], /* --> pour chaque cellule, position
dans ielvse de la première cellule
du voisinage étendu */
const cs_int_t ielvse[], /* --> numéro des cellules du voisinage
étendu */
cs_int_t ituser[], /* <-> tab. complémentaire utilisateur */
cs_int_t ia[], /* <-> macro-tableau entier */
const cs_real_t volume[], /* --> volumes des cellules */
const cs_real_t surfac[], /* --> surfaces des faces internes */
const cs_real_t surfbo[], /* --> surfaces des faces de bord */
const cs_real_t surfbn[], /* --> norme de surfbo */
const cs_real_t pond[], /* --> pondération géométrique faces int*/
const cs_real_t dist[], /* --> distance entre I' et J' faces int*/
const cs_real_t distbr[], /* --> distance du cdgfbo à I' fac. brd */
const cs_real_t dijpf[], /* --> vecteur I'J' aux faces internes */
const cs_real_t diipb[], /* --> vecteur II' aux faces de bord */
const cs_real_t fextx[], /* --> composantes de la force */
const cs_real_t fexty[], /* extérieure générant la pression */
const cs_real_t fextz[], /* hydrostatique */
const cs_real_t xyzcen[], /* --> c.d.g. des cellules */
const cs_real_t cdgfac[], /* --> c.d.g. des faces internes */
const cs_real_t cdgfbo[], /* --> c.d.g. des faces de bord */
const cs_real_t coefap[], /* --> cond. limites aux faces de bord */
const cs_real_t coefbp[], /* --> cond. limites aux faces de bord */
const cs_real_t pvar[], /* --> variable dont on calcule le grad */
cs_real_t cocgb[], /* <-> contribution à COCG des faces
internes des cellules de bord */
cs_real_t cocg[], /* <-> contribution à COCG des faces
internes des cellules de bord */
cs_real_t rtuser[], /* <-> tab. complémentaire utilisateur */
cs_real_t dpdx[], /* <-- composante x du gradient */
cs_real_t dpdy[], /* <-- composante y du gradient */
cs_real_t dpdz[], /* <-- composante z du gradient */
cs_real_t bx[], /* - tab. de travail local */
cs_real_t by[], /* - tab. de travail local */
cs_real_t bz[], /* - tab. de travail local */
cs_real_t ra[] /* <-> macro-tableau réel */
);
/*----------------------------------------------------------------------------
* Limitation des gradients (standard ou avec voisinage étendu)
*----------------------------------------------------------------------------*/
void CS_PROCF (limgrd, LIMGRD)
(
const cs_int_t *const ndim, /* --> dimension de l'espace */
const cs_int_t *const ncelet, /* --> nombre de cellules étendu */
const cs_int_t *const ncel, /* --> nombre de cellules */
const cs_int_t *const nfac, /* --> nombre de faces internes */
const cs_int_t *const imrgra, /* --> mode de calcul du gradient */
const cs_int_t *const imligp, /* --> mode de limitation du gradient */
const cs_int_t *const idimte, /* --> 0, 1, 2 : variable
scalaire, vecteur, tenseur */
const cs_int_t *const itenso, /* --> pour la périodicité de rotation */
const cs_int_t *const iwarnp, /* --> niveau d'impression */
const cs_int_t *const nfecra, /* --> unité sortie standard */
const cs_real_t *const climgp, /* --> coefficient de limitation du grad*/
const cs_int_t ifacel[], /* --> liste des faces internes */
const cs_int_t ipcvse[], /* --> pour chaque cellule, position
dans ielvse de la première cellule
du voisinage étendu */
const cs_int_t ielvse[], /* --> numéro des cellules du voisinage
étendu */
const cs_real_t xyzcen[], /* --> c.d.g. des cellules */
const cs_real_t pvar[], /* --> variable dont on calcule le grad */
cs_real_t dpdx[], /* <-> composante x du gradient */
cs_real_t dpdy[], /* <-> composante y du gradient */
cs_real_t dpdz[], /* <-> composante z du gradient */
cs_real_t faclim[], /* - tab. de travail local */
cs_real_t denom[], /* - tab. de travail local */
cs_real_t denum[] /* - tab. de travail local */
);
void CS_PROCF (filtre, FILTRE)
(
const cs_int_t *const ncelet, /* --> nombre de cellules étendu */
const cs_int_t *const ncel, /* --> nombre de cellules */
const cs_int_t *const nfac, /* --> nombre de faces internes */
const cs_int_t *const nfabor, /* --> nombre de faces de bord */
const cs_int_t ifacel[], /* --> liste des faces internes */
const cs_int_t ifabor[], /* --> liste des faces de bord */
const cs_int_t ipcvse[], /* --> pour chaque cellule, position
dans ielvse de la première cellule
du voisinage étendu */
const cs_int_t ielvse[], /* --> numéro des cellules du voisinage
étendu */
const cs_real_t volume[], /* --> volumes des cellules */
const cs_real_t pvar[], /* --> variable a filtrer */
cs_real_t varfil[], /* --> resultat de la variable filtree */
cs_real_t w1[], /* - tab. de travail local */
cs_real_t w2[] /* - tab. de travail local */
);
/*============================================================================
* Fonctions publiques
*============================================================================*/
/*----------------------------------------------------------------------------
* Definition du voisinage étendu pour le calcul des gradients par moindres
* carrés : au premier passage, on choisit les cellules et on modifie la
* structure maillage qui stocke le voisinage étendu.
* ensuite, on se contente de retourner les pointeurs
*
* Attention : le voisinage étendu permet de toucher toutes les cellules
* qui partagent un noeud avec une cellule donnée, mais il ne
* comporte pas les cellules visibles d'une cellule au travers
* des faces (elles sont accessibles par IFACEL)
*
* Interface Fortran :
*
* SUBROUTINE DEFVSE
* *****************
* & ( NDIM , NCEL , NFAC , NNOD ,
* & IMRGRA , ANOMAX ,
* & IFACEL , IPNFAC , NODFAC , XYZCEN , SURFAC )
*
*----------------------------------------------------------------------------*/
void CS_PROCF (defvse, DEFVSE)
(
const cs_int_t *const ndim, /* --> Dimension = 3 */
const cs_int_t *const ncel, /* --> Nombre de cellules */
const cs_int_t *const nfac, /* --> Nombre de faces internes */
const cs_int_t *const nnod, /* --> Nombre de sommets */
const cs_int_t *const imrgra, /* --> Mode de reconstruction */
const cs_real_t *const anomax, /* --> Angle de non orthogonalité en radian
au-delà duquel on retient dans le
voisinage étendu de chaque cellule
voisine les cellules dont un sommet
est sur la face */
const cs_int_t *const ifacel, /* --> Connectivité face interne -> cellules*/
const cs_int_t *const ipnfac, /* --> Position du premier noeud d'une face
dans nodfac */
const cs_int_t *const nodfac, /* --> Numéro des noeuds des faces internes */
const cs_real_t *const xyzcen, /* --> Coord. du centre des cellules */
const cs_real_t *const surfac /* --> Vecteur surface des faces internes */
)
{
cs_maillage_t *maillage = cs_glob_maillage;
cs_int_t ind_som;
cs_int_t ind_fac;
cs_int_t ind_cel;
cs_int_t ind_cel_i;
cs_int_t ind_cel_j;
cs_int_t ind_bcl_cel;
cs_int_t ind_bcl_cel_vse;
cs_int_t ind_bcl_fac;
cs_int_t num_cel;
cs_int_t num_cel_vse;
cs_int_t nbr_val_connect_som_cel;
cs_int_t nbr_val_connect_som_cel_estime;
cs_int_t nbr_cel_elimine;
cs_int_t ind_val_last;
cs_int_t ind_pos_last;
cs_int_t icoo;
cs_int_t *ielvse;
cs_int_t *ipcvse;
cs_int_t *pos_connect_som_fac;
cs_int_t *val_connect_som_fac;
cs_int_t *pos_connect_som_cel;
cs_int_t *val_connect_som_cel;
cs_real_t vij[CS_DIM_3];
cs_real_t vnf[CS_DIM_3];
cs_real_t normeij;
cs_real_t normenf;
cs_real_t cos_ij_nf;
cs_real_t cos_ij_nf_min;
cs_bool_t bool_cel_vu;
const cs_real_t *xyzcen_avec_voiset;
static cs_int_t cs_loc_premier_appel = 0;
enum {X, Y, Z};
#define CS_LOC_MODULE(vect) \
sqrt(vect[X] * vect[X] + vect[Y] * vect[Y] + vect[Z] * vect[Z])
/* Au premier passage, on sélectionne les cellules */
if (cs_loc_premier_appel == 0) {
/* Pour la suite, on note que l'on est passé une fois */
cs_loc_premier_appel = 1;
/* Si on n'a pas de voisinage étendu, on s'arrête */
ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
ielvse = maillage->val_connect_voiset_cel_cel_dom;
if ( ipcvse == NULL || ielvse == NULL
|| maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SANS) {
bft_error (__FILE__, __LINE__, 0,
_("Lorsque la méthode de calcul des gradients par moindres\n"
"carrés sur support étendu est activée (IMRGRA = <%d>), \n"
"il est indispensable de générer un voisinage étendu \n"
"en utilisant l'Enveloppe en pré-traitement avec \n"
"l'option -ve (COMMANDE_VOISET dans le lanceur)."),
(int)(*imrgra));
}
/* Selon imrgra, on sélectionne le voisinage de diverses façons */
/*
Par défaut : on conserve tout
=============================
*/
/*
imrgra = 3 : sélection sur critère de non orthogonalité
=======================================================
pour chaque face interne, au-dessus d'un angle de non
orthogonalité donné, on retient dans le voisinage étendu
des deux cellules voisines de la face toutes les cellules
qui s'appuient sur un des sommets de la face
*/
if ((*imrgra) == 3) {
/*
On construit d'abord la connectivité sommet -> cellules
------------------------------------------------------
On a déjà face -> sommets, il faut l'inverser
on utilisera ensuite face -> cellules
*/
/*
On construit la connectivité sommet -> faces
en inversant face -> sommets
*/
/* Allocation de pos_connect_som_fac */
BFT_MALLOC(pos_connect_som_fac, (*nnod)+1, cs_int_t);
/* On compte les faces de chaque sommet, on utilise pos_connect_som_fac*/
for (ind_som = 0 ; ind_som < (*nnod)+1 ; ind_som++)
pos_connect_som_fac[ind_som] = 0;
for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++)
for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ;
ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
ind_som = nodfac[ind_bcl_fac] - 1;
pos_connect_som_fac[ind_som+1]++;
}
/* Construction de pos_connect_som_fac */
pos_connect_som_fac[0] = 1;
for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++)
pos_connect_som_fac[ind_som+1] += pos_connect_som_fac[ind_som];
/* Allocations de val_connect_som_fac */
BFT_MALLOC(val_connect_som_fac, pos_connect_som_fac[(*nnod)]-1, cs_int_t);
/* Construction de val_connect_som_fac : pour cela, on utilise
pos_connect_som_fac pour repérer la case à remplir dans
val_connect_som_fac : on l'incrémente donc à chaque nouvelle valeur
ajoutée. A la fin, on a donc décalé les valeurs de pos_connect_som_fac
d'une case et il faut les remettre en place. */
for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++)
for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ;
ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
ind_som = nodfac[ind_bcl_fac] - 1;
val_connect_som_fac[pos_connect_som_fac[ind_som]-1] = ind_fac + 1;
pos_connect_som_fac[ind_som]++;
}
for (ind_som = (*nnod) ; ind_som > 0 ; ind_som--)
pos_connect_som_fac[ind_som] = pos_connect_som_fac[ind_som-1];
pos_connect_som_fac[0] = 1;
#if 0
for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) {
printf(" sommet %d :: ", ind_som+1) ;
for (ind_bcl_cel = pos_connect_som_fac[ind_som ] - 1 ;
ind_bcl_cel < pos_connect_som_fac[ind_som+1] - 1 ;
ind_bcl_cel++)
printf(" %d ", val_connect_som_fac[ind_bcl_cel]) ;
printf("\n") ;
}
#endif
/*
On construit la connectivité sommet -> cellules
avec sommet -> faces et face -> cellules
*/
/* Allocations */
nbr_val_connect_som_cel_estime = 3*(*nnod);
BFT_MALLOC(val_connect_som_cel, nbr_val_connect_som_cel_estime, cs_int_t);
BFT_MALLOC(pos_connect_som_cel, (*nnod)+1, cs_int_t);
/* Initialisations */
nbr_val_connect_som_cel = 0;
pos_connect_som_cel[0] = 1;
/* Pour tous les sommets */
for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) {
/* On balaye les faces */
for (ind_bcl_fac = pos_connect_som_fac[ind_som ] - 1 ;
ind_bcl_fac < pos_connect_som_fac[ind_som+1] - 1 ; ind_bcl_fac++) {
ind_fac = val_connect_som_fac[ind_bcl_fac] - 1;
/* Et on note les cellules (il faut absolument
qu'elles ne soient pas en plusieurs exemplaires */
/* Voisin 1 */
num_cel = ifacel[2*ind_fac ];
/* Test si deja vu */
bool_cel_vu = CS_FALSE;
ind_bcl_cel = pos_connect_som_cel[ind_som] - 1;
while ( (bool_cel_vu == CS_FALSE)
&& (ind_bcl_cel < nbr_val_connect_som_cel)) {
if (num_cel == val_connect_som_cel[ind_bcl_cel])
bool_cel_vu = CS_TRUE;
ind_bcl_cel++;
}
/* Si pas encore vu, on le note */
if (bool_cel_vu == CS_FALSE) {
if (nbr_val_connect_som_cel + 1 > nbr_val_connect_som_cel_estime) {
nbr_val_connect_som_cel_estime *= 2;
BFT_REALLOC(val_connect_som_cel,
nbr_val_connect_som_cel_estime, cs_int_t);
}
val_connect_som_cel[nbr_val_connect_som_cel] = num_cel;
nbr_val_connect_som_cel++;
}
/* Voisin 2 */
num_cel = ifacel[2*ind_fac+1];
/* Test si deja vu */
bool_cel_vu = CS_FALSE;
ind_bcl_cel = pos_connect_som_cel[ind_som] - 1;
while ( (bool_cel_vu == CS_FALSE)
&& (ind_bcl_cel < nbr_val_connect_som_cel)) {
if (num_cel == val_connect_som_cel[ind_bcl_cel])
bool_cel_vu = CS_TRUE;
ind_bcl_cel++;
}
/* Si pas encore vu, on le note */
if (bool_cel_vu == CS_FALSE) {
if (nbr_val_connect_som_cel + 1 > nbr_val_connect_som_cel_estime) {
nbr_val_connect_som_cel_estime *=2;
BFT_REALLOC(val_connect_som_cel,
nbr_val_connect_som_cel_estime, cs_int_t);
}
val_connect_som_cel[nbr_val_connect_som_cel] = num_cel;
nbr_val_connect_som_cel++;
}
}
pos_connect_som_cel[ind_som+1] = nbr_val_connect_som_cel + 1;
}
/* Réallocation */
BFT_REALLOC(val_connect_som_cel, nbr_val_connect_som_cel, cs_int_t);
/* Libération de sommet -> faces */
BFT_FREE(val_connect_som_fac);
BFT_FREE(pos_connect_som_fac);
#if 0
for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) {
printf(" sommet %d :: ", ind_som+1) ;
for (ind_bcl_cel = pos_connect_som_cel[ind_som ] - 1 ;
ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ;
ind_bcl_cel++)
printf(" %d ", val_connect_som_cel[ind_bcl_cel]) ;
printf("\n") ;
}
#endif
/*
On sélectionne enfin les cellules à conserver dans le voisinage
---------------------------------------------------------------
On les marque d'abord, on élimine celles qui sont inutiles ensuite
*/
/*
Marquage des cellules à éliminer : en négatif
*/
assert((*ndim) == CS_DIM_3);
/* Coordonnées du centre */
if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
&& maillage->nbr_cel_fac_par_voiset > 0) {
xyzcen_avec_voiset = maillage->coord_cel_avec_voiset;
}
else
xyzcen_avec_voiset = xyzcen;
/* Sélection par les faces internes */
for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++) {
/* On calcule le cosinus de l'angle de non orthogonalité
des faces internes (angle entre la normale et la droite
passant par le centre des cellules voisines */
/* Vecteur IJ */
ind_cel_i = ifacel[2*ind_fac ] - 1;
ind_cel_j = ifacel[2*ind_fac+1] - 1;
icoo = 0;
vij[icoo]
= xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo]
- xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo];
icoo = 1;
vij[icoo]
= xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo]
- xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo];
icoo = 2;
vij[icoo]
= xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo]
- xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo];
normeij = CS_LOC_MODULE(vij);
assert (normeij > 0.);
/* Vecteur normal à la face */
icoo = 0;
vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo];
icoo = 1;
vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo];
icoo = 2;
vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo];
normenf = CS_LOC_MODULE(vnf);
assert (normenf > 0.);
/* Produit scalaire vij.vnf */
cos_ij_nf = (vij[0]*vnf[0]+vij[1]*vnf[1]+vij[2]*vnf[2])
/ (normeij*normenf);
/* Comparaison à une limite prédéfinie,
En dessous de cette limite, c'est non orthogonal et
on garde la cellule dans le voisinage étendu des deux voisines
de la face (on la marque en négatif, puis on changera tous les
signes, pour que les cellules à éliminer soient en négatif */
cos_ij_nf_min = cos((*anomax));
/* Test sur l'angle de non orthogonalité */
if (cos_ij_nf < cos_ij_nf_min) {
/* Pour chaque cellule voisine de la face :
intersection du groupe des cellules du voisinage étendu
que l'on n'a pas encore décidé de garder avec
le groupe de cellules possédant un sommet sur la face */
/* Cellule voisine 1 (uniquement si elle est réelle) */
if (ind_cel_i < (*ncel)) {
ind_cel = ind_cel_i;
/* Voisinage étendu non encore retenu */
for (ind_bcl_cel_vse = ipcvse[ind_cel ] - 1 ;
ind_bcl_cel_vse < ipcvse[ind_cel+1] - 1 ;
ind_bcl_cel_vse++) {
num_cel_vse = ielvse[ind_bcl_cel_vse];
if (num_cel_vse > 0)
/* Cellules partageant un sommet avec la face */
for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ;
ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
ind_som = nodfac[ind_bcl_fac] - 1;
for (ind_bcl_cel = pos_connect_som_cel[ind_som ] - 1 ;
ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ;
ind_bcl_cel++) {
num_cel = val_connect_som_cel[ind_bcl_cel];
/* Comparaison et sélection */
if (num_cel == num_cel_vse && ielvse[ind_bcl_cel_vse] > 0)
ielvse[ind_bcl_cel_vse] = - ielvse[ind_bcl_cel_vse];
}
}
}
}
/* Cellule voisine 2 (uniquement si elle est réelle) */
if (ind_cel_j < (*ncel)) {
ind_cel = ind_cel_j;
/* Voisinage étendu non encore retenu */
for (ind_bcl_cel_vse = ipcvse[ind_cel ] - 1 ;
ind_bcl_cel_vse < ipcvse[ind_cel+1] - 1 ;
ind_bcl_cel_vse++) {
num_cel_vse = ielvse[ind_bcl_cel_vse];
if (num_cel_vse > 0)
/* Cellules partageant un sommet avec la face */
for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ;
ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
ind_som = nodfac[ind_bcl_fac] - 1;
for (ind_bcl_cel = pos_connect_som_cel[ind_som ] - 1 ;
ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ;
ind_bcl_cel++) {
num_cel = val_connect_som_cel[ind_bcl_cel];
/* Comparaison et sélection */
if (num_cel == num_cel_vse && ielvse[ind_bcl_cel_vse] > 0)
ielvse[ind_bcl_cel_vse] = - ielvse[ind_bcl_cel_vse] ;
}
}
}
}
}
}
/* Libération de sommet -> cellules */
BFT_FREE(val_connect_som_cel);
BFT_FREE(pos_connect_som_cel);
/* On change tous les signes de ielvse, pour que les cellules à éliminer
soient en négatif */
for (ind_bcl_cel = 0 ; ind_bcl_cel < ipcvse[(*ncel)]-1 ; ind_bcl_cel++)
ielvse[ind_bcl_cel] = - ielvse[ind_bcl_cel];
/* On élimine les cellules marquées négativement */
nbr_cel_elimine = 0;
ind_val_last = -1;
ind_pos_last = 0;
for (ind_cel = 0 ; ind_cel < (*ncel) ; ind_cel++) {
for (ind_bcl_cel = ind_pos_last ;
ind_bcl_cel < ipcvse[ind_cel+1] - 1 ; ind_bcl_cel++) {
if (ielvse[ind_bcl_cel] > 0) {
ind_val_last++;
ielvse[ind_val_last] = ielvse[ind_bcl_cel];
}
else {
nbr_cel_elimine++;
}
}
ind_pos_last = ipcvse[ind_cel+1] - 1;
ipcvse[ind_cel+1] = ipcvse[ind_cel+1] - nbr_cel_elimine;
}
maillage->nbr_val_connect_voiset_cel_cel_dom
= maillage->nbr_val_connect_voiset_cel_cel_dom - nbr_cel_elimine;
/* Réallocation avec la dimension exacte */
BFT_REALLOC(maillage->val_connect_voiset_cel_cel_dom,
ipcvse[(*ncel)]-1, cs_int_t);
bft_printf(_("\n %d cellules conservées dans le voisinage étendu\n"),
maillage->nbr_val_connect_voiset_cel_cel_dom);
#if 0
for (ind_cel = 0 ; ind_cel < (*ncel) ; ind_cel++) {
printf(" cellule %d :: ", ind_cel+1);
for (ind_bcl_cel = ipcvse[ind_cel ] - 1 ;
ind_bcl_cel < ipcvse[ind_cel+1] - 1 ;
ind_bcl_cel++)
printf(" %d ", ielvse[ind_bcl_cel]);
printf("\n");
}
#endif
}
}
}
/*----------------------------------------------------------------------------
* Appel de la routine FORTRAN GRADMC pour le calcul des gradients par moindres
* carrés. On ajoute aux arguments le voisinage étendu de la structure
* maillage.
*
* Interface Fortran :
*
* SUBROUTINE CGRDMC
* *****************
*
* & ( IDBIA0 , IDBRA0 ,
* & NCELET , NCEL , NFAC , NFABOR , NCELBR , NITUSE , NRTUSE ,
* & INC , ICCOCG , NSWRGP , IDIMTE , ITENSO , IPHYDP , IMRGRA ,
* & IWARNP , NFECRA , EPSRGP , EXTRAP ,
* & IFACEL , IFABOR , ICELBR , ITUSER , IA ,
* & VOLUME , SURFAC , SURFBO , SURFBN , POND ,
* & DIST , DISTBR , DIJPF , DIIPB ,
* & FEXTX , FEXTY , FEXTZ ,
* & XYZCEN , CDGFAC , CDGFBO , COEFAP , COEFBP , PVAR ,
* & COCGB , COCG , RTUSER ,
* & DPDX , DPDY , DPDZ ,
* & BX , BY , BZ , RA )
*
*----------------------------------------------------------------------------*/
void CS_PROCF (cgrdmc, CGRDMC)
(
const cs_int_t *const idbia0, /* --> numéro 1ère case libre dans IA */
const cs_int_t *const idbra0, /* --> numéro 1ère case libre dans RA */
const cs_int_t *const ncelet, /* --> nombre de cellules étendu */
const cs_int_t *const ncel, /* --> nombre de cellules */
const cs_int_t *const nfac, /* --> nombre de faces internes */
const cs_int_t *const nfabor, /* --> nombre de faces de bord */
const cs_int_t *const ncelbr, /* --> nombre de cellules de bord */
const cs_int_t *const nituse, /* --> longueur du tableau ituser[] */
const cs_int_t *const nrtuse, /* --> longueur du tableau rtuser[] */
const cs_int_t *const inc, /* --> 0 ou 1 : incrément ou non */
const cs_int_t *const iccocg, /* --> 1 ou 0 : recalcul de COCG ou non */
const cs_int_t *const nswrgp, /* --> >1 : avec reconstruction */
const cs_int_t *const idimte, /* --> 0, 1, 2 : variable
scalaire, vecteur, tenseur */
const cs_int_t *const itenso, /* --> pour la périodicité de rotation */
const cs_int_t *const iphydp, /* --> prise en compte de P hydrostatiqu*/
const cs_int_t *const imrgra, /* --> mode de calcul du gradient */
const cs_int_t *const iwarnp, /* --> niveau d'impression */
const cs_int_t *const nfecra, /* --> unité sortie standard */
const cs_real_t *const epsrgp, /* --> précision pour calcul du gradient
itératif */
const cs_real_t *const extrap, /* --> extrapolation des grad. au bord */
const cs_int_t ifacel[], /* --> liste des faces internes */
const cs_int_t ifabor[], /* --> liste des faces de bord */
const cs_int_t icelbr[], /* --> numéro des cellules ayant au
moins une face de bord */
cs_int_t ituser[], /* <-> tab. complémentaire utilisateur */
cs_int_t ia[], /* <-> macro-tableau entier */
const cs_real_t volume[], /* --> volumes des cellules */
const cs_real_t surfac[], /* --> surfaces des faces internes */
const cs_real_t surfbo[], /* --> surfaces des faces de bord */
const cs_real_t surfbn[], /* --> norme de surfbo */
const cs_real_t pond[], /* --> pondération géométrique faces int*/
const cs_real_t dist[], /* --> distance entre I' et J' faces int*/
const cs_real_t distbr[], /* --> distance du cdgfbo à I' fac. brd */
const cs_real_t dijpf[], /* --> vecteur I'J' aux faces internes */
const cs_real_t diipb[], /* --> vecteur II' aux faces de bord */
const cs_real_t fextx[], /* --> composantes de la force */
const cs_real_t fexty[], /* extérieure générant la pression */
const cs_real_t fextz[], /* hydrostatique */
const cs_real_t xyzcen[], /* --> c.d.g. des cellules */
const cs_real_t cdgfac[], /* --> c.d.g. des faces internes */
const cs_real_t cdgfbo[], /* --> c.d.g. des faces de bord */
const cs_real_t coefap[], /* --> cond. limites aux faces de bord */
const cs_real_t coefbp[], /* --> cond. limites aux faces de bord */
const cs_real_t pvar[], /* --> variable dont on calcule le grad */
cs_real_t cocgb[], /* <-> contribution à COCG des faces
internes des cellules de bord */
cs_real_t cocg[], /* <-> contribution à COCG des faces
internes des cellules de bord */
cs_real_t rtuser[], /* <-> tab. complémentaire utilisateur */
cs_real_t dpdx[], /* <-- composante x du gradient */
cs_real_t dpdy[], /* <-- composante y du gradient */
cs_real_t dpdz[], /* <-- composante z du gradient */
cs_real_t bx[], /* - tab. de travail local */
cs_real_t by[], /* - tab. de travail local */
cs_real_t bz[], /* - tab. de travail local */
cs_real_t ra[] /* <-> macro-tableau réel */
)
{
#if defined(_CS_HAVE_MPI)
cs_int_t iel;
cs_int_t nceltot;
#endif /*_CS_HAVE_MPI */
cs_int_t *ipcvse;
cs_int_t *ielvse;
const cs_real_t *xyzcen_avec_voiset;
const cs_real_t *pvar_avec_voiset;
const cs_real_t *fextx_avec_voiset;
const cs_real_t *fexty_avec_voiset;
const cs_real_t *fextz_avec_voiset;
#if defined(_CS_HAVE_MPI)
cs_real_t *pvar_voiset;
cs_real_t *fextx_voiset;
cs_real_t *fexty_voiset;
cs_real_t *fextz_voiset;
cs_real_t *fextx_avec_voiset_tmp;
cs_real_t *fexty_avec_voiset_tmp;
cs_real_t *fextz_avec_voiset_tmp;
#endif /*_CS_HAVE_MPI */
cs_maillage_t *maillage = cs_glob_maillage;
#if defined(_CS_HAVE_MPI)
cs_maillage_tmp_t *maillage_tmp = cs_glob_maillage_tmp;
#endif /*_CS_HAVE_MPI */
/*
Remarque sur les tests faisant intervenir _CS_HAVE_MPI
Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie
definie que si on _CS_HAVE_MPI (si le parallélisme est possible)
Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est
correct car il ne peut y avoir un voisinage separé que en parallele
etant donne le choix fait dans l'enveloppe (un voisinage separe
en non parallele n'apporterait rien car il ne permettrait pas de
reduire le volume des echanges).
On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI
n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t
par des tests sur _CS_HAVE_MPI pour permettre la compilation
lorsque _CS_HAVE_MPI n'est pas defini
Par coherence, on fait de meme avec les variables pvar_voiset
fext?_voiset et fext?_avec_voiset_tmp, nceltot, iel
*/
/* La complexite du traitement dans cette fonction est liée à la prise en
compte d'un voisinage étendu séparé du halo classique */
/* On va transmettre au gradient par moindre carrés une variable
pvar_avec_voiset qui sera :
maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe
et qu'il est traité séparément
pvar sinon (cas classiques)
La variable maillage_tmp->var_cel_avec_voiset devra contenir
pvar dans les ncelet premieres cases et les valeurs actualisees
du voisinage étendu séparé dans le reste, noté pvar_voiset
Il faut de même transmettre la variable allongée correcte pour les
centres de gravité des cellules et,
lorsque la pression hydrostatique est prise en compte, transmettre
également des variables allongées pour la fext.
L'utilisation des pointeurs fext?_avec_voiset_tmp au lieu de
fext?_avec_voiset directement permet d'éviter un warning et de conserver
le caractere const des fext? en argument.
*/
/*
Echanges pour mise à jour du voisinage étendu séparé
----------------------------------------------------
*/
#if defined(_CS_HAVE_MPI)
if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) {
/* Identification des pointeurs sur le voisinage étendu séparé */
if (maillage->nbr_cel_fac_par_voiset > 0) {
pvar_voiset = maillage_tmp->var_cel_avec_voiset + (*ncelet);
}
else {
pvar_voiset = NULL;
}
if (*iphydp != 0) {
if (maillage->nbr_cel_fac_par_voiset > 0) {
nceltot = (*ncelet)+maillage->nbr_cel_fac_par_voiset;
BFT_MALLOC(fextx_avec_voiset_tmp,nceltot,cs_real_t);
BFT_MALLOC(fexty_avec_voiset_tmp,nceltot,cs_real_t);
BFT_MALLOC(fextz_avec_voiset_tmp,nceltot,cs_real_t);
fextx_voiset = fextx_avec_voiset_tmp + (*ncelet);
fexty_voiset = fexty_avec_voiset_tmp + (*ncelet);
fextz_voiset = fextz_avec_voiset_tmp + (*ncelet);
}
else {
fextx_avec_voiset_tmp = NULL;
fexty_avec_voiset_tmp = NULL;
fextz_avec_voiset_tmp = NULL;
fextx_voiset = NULL;
fexty_voiset = NULL;
fextz_voiset = NULL;
}
}
/* Echanges parallèles et périodiques
Attention, un processeur peut n'avoir rien à recevoir
(maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant
des éléments à envoyer : on doit donc entrer dans parcve
pour tous les processeurs lorsque l'on est en parallèle */
if (maillage->nbr_dom > 1)
cs_parallel_parcve(pvar, pvar_voiset);
if (maillage->nbr_per > 0)
cs_perio_percve (pvar, pvar_voiset);
if (*iphydp != 0) {
if (maillage->nbr_dom > 1) {
cs_parallel_parcve(fextx, fextx_voiset);
cs_parallel_parcve(fexty, fexty_voiset);
cs_parallel_parcve(fextz, fextz_voiset);
}
if (maillage->nbr_per > 0) {
cs_perio_percve (fextx, fextx_voiset);
cs_perio_percve (fexty, fexty_voiset);
cs_perio_percve (fextz, fextz_voiset);
}
}
}
#endif /*_CS_HAVE_MPI */
/*
Appel du gradient par moindres carres
-------------------------------------
*/
/* Identification du pointeur sur la variable allongée
(de longueur ncelet + voisinage étendu séparé éventuel)
et copie des valeurs 1-ncelet si nécessaire */
#if defined(_CS_HAVE_MPI)
if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
&& maillage->nbr_cel_fac_par_voiset > 0) {
for (iel = 0 ; iel < (*ncelet) ; iel++)
maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel];
pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset;
if (*iphydp != 0) {
for (iel = 0 ; iel < (*ncelet) ; iel++)
fextx_avec_voiset_tmp[iel] = fextx[iel];
fextx_avec_voiset = fextx_avec_voiset_tmp;
for (iel = 0 ; iel < (*ncelet) ; iel++)
fexty_avec_voiset_tmp[iel] = fexty[iel];
fexty_avec_voiset = fexty_avec_voiset_tmp;
for (iel = 0 ; iel < (*ncelet) ; iel++)
fextz_avec_voiset_tmp[iel] = fextz[iel];
fextz_avec_voiset = fextz_avec_voiset_tmp;
}
else {
fextx_avec_voiset = fextx;
fexty_avec_voiset = fexty;
fextz_avec_voiset = fextz;
}
}
else {
#endif /*_CS_HAVE_MPI */
pvar_avec_voiset = pvar;
fextx_avec_voiset = fextx;
fexty_avec_voiset = fexty;
fextz_avec_voiset = fextz;
#if defined(_CS_HAVE_MPI)
}
#endif /*_CS_HAVE_MPI */
/* Identification du pointeur sur les coordonnées allongées
(de longueur ncelet + voisinage étendu séparé éventuel)*/
if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
&& maillage->nbr_cel_fac_par_voiset > 0) {
xyzcen_avec_voiset = maillage->coord_cel_avec_voiset;
}
else
xyzcen_avec_voiset = xyzcen;
/* Pointeurs sur le voisinage étendu */
ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
ielvse = maillage->val_connect_voiset_cel_cel_dom;
/* Appel du calcul de gradient */
CS_PROCF(gradmc, GRADMC)
(idbia0, idbra0,
ncelet, ncel , nfac , nfabor, ncelbr, nituse, nrtuse,
inc , iccocg, nswrgp, idimte, itenso, iphydp, imrgra,
iwarnp, nfecra, epsrgp, extrap,
ifacel, ifabor, icelbr, ipcvse, ielvse, ituser, ia ,
volume, surfac, surfbo, surfbn, pond ,
dist , distbr, dijpf , diipb ,
fextx_avec_voiset, fexty_avec_voiset, fextz_avec_voiset,
xyzcen_avec_voiset,
cdgfac, cdgfbo, coefap, coefbp, pvar_avec_voiset,
cocgb , cocg , rtuser,
dpdx , dpdy , dpdz ,
bx , by , bz , ra );
/* Libération si on a réservé quelque chose */
#if defined(_CS_HAVE_MPI)
if ( (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE)
&& (*iphydp != 0)
&& (maillage->nbr_cel_fac_par_voiset > 0)) {
BFT_FREE(fextx_avec_voiset_tmp);
BFT_FREE(fexty_avec_voiset_tmp);
BFT_FREE(fextz_avec_voiset_tmp);
}
#endif /*_CS_HAVE_MPI */
}
/*----------------------------------------------------------------------------
* Appel de la routine FORTRAN LIMGRD pour la limitation éventuelle des
* gradients. On ajoute aux arguments le voisinage étendu de la structure
* maillage.
*
* Interface Fortran :
*
* SUBROUTINE CLMGRD
* *****************
*
* & ( NDIM , NCELET , NCEL , NFAC ,
* & IMRGRA , IMLIGP , IDIMTE , ITENSO , IWARNP , NFECRA , CLIMGP ,
* & IFACEL ,
* & XYZCEN , PVAR ,
* & DPDX , DPDY , DPDZ ,
* & FACLIM , DENOM , DENUM )
*
*----------------------------------------------------------------------------*/
void CS_PROCF (clmgrd, CLMGRD)
(
const cs_int_t *const ndim, /* --> dimension de l'espace */
const cs_int_t *const ncelet, /* --> nombre de cellules étendu */
const cs_int_t *const ncel, /* --> nombre de cellules */
const cs_int_t *const nfac, /* --> nombre de faces internes */
const cs_int_t *const imrgra, /* --> mode de calcul du gradient */
const cs_int_t *const imligp, /* --> mode de limitation du gradient */
const cs_int_t *const idimte, /* --> 0, 1, 2 : variable
scalaire, vecteur, tenseur */
const cs_int_t *const itenso, /* --> pour la périodicité de rotation */
const cs_int_t *const iwarnp, /* --> niveau d'impression */
const cs_int_t *const nfecra, /* --> unité sortie standard */
const cs_real_t *const climgp, /* --> coefficient de limitation du grad*/
const cs_int_t ifacel[], /* --> liste des faces internes */
const cs_real_t xyzcen[], /* --> c.d.g. des cellules */
const cs_real_t pvar[], /* --> variable dont on calcule le grad */
cs_real_t dpdx[], /* <-> composante x du gradient */
cs_real_t dpdy[], /* <-> composante y du gradient */
cs_real_t dpdz[], /* <-> composante z du gradient */
cs_real_t faclim[], /* - tab. de travail local */
cs_real_t denom[], /* - tab. de travail local */
cs_real_t denum[] /* - tab. de travail local */
)
{
#if defined(_CS_HAVE_MPI)
cs_int_t iel;
cs_int_t nceltot;
#endif /*_CS_HAVE_MPI */
cs_int_t *ipcvse;
cs_int_t *ielvse;
const cs_real_t *xyzcen_avec_voiset;
const cs_real_t *pvar_avec_voiset;
cs_real_t *dpdx_avec_voiset;
cs_real_t *dpdy_avec_voiset;
cs_real_t *dpdz_avec_voiset;
cs_real_t *denum_avec_voiset;
cs_real_t *denom_avec_voiset;
#if defined(_CS_HAVE_MPI)
cs_real_t *pvar_voiset;
cs_real_t *dpdx_voiset;
cs_real_t *dpdy_voiset;
cs_real_t *dpdz_voiset;
cs_real_t *dpdx_avec_voiset_tmp;
cs_real_t *dpdy_avec_voiset_tmp;
cs_real_t *dpdz_avec_voiset_tmp;
cs_real_t *denum_avec_voiset_tmp;
cs_real_t *denom_avec_voiset_tmp;
#endif /*_CS_HAVE_MPI */
cs_maillage_t *maillage = cs_glob_maillage;
#if defined(_CS_HAVE_MPI)
cs_maillage_tmp_t *maillage_tmp = cs_glob_maillage_tmp;
#endif /*_CS_HAVE_MPI */
/*
Remarque sur les tests faisant intervenir _CS_HAVE_MPI
Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie
definie que si on _CS_HAVE_MPI (si le parallélisme est possible)
Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est
correct car il ne peut y avoir un voisinage separé que en parallele
etant donne le choix fait dans l'enveloppe (un voisinage separe
en non parallele n'apporterait rien car il ne permettrait pas de
reduire le volume des echanges).
On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI
n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t
par des tests sur _CS_HAVE_MPI pour permettre la compilation
lorsque _CS_HAVE_MPI n'est pas defini
Par coherence, on fait de meme avec les variables pvar_voiset
fext?_voiset et fext?_avec_voiset_tmp, nceltot, iel
*/
/* La complexite du traitement dans cette fonction est liée à la prise en
compte d'un voisinage étendu séparé du halo classique */
/* On va transmettre au gradient par moindre carrés une variable
pvar_avec_voiset qui sera :
maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe
et qu'il est traité séparément
pvar sinon (cas classiques)
La variable maillage_tmp->var_cel_avec_voiset devra contenir
pvar dans les ncelet premieres cases et les valeurs actualisees
du voisinage étendu séparé dans le reste, noté pvar_voiset
Il faut de même transmettre la variable allongée correcte pour les
centres de gravité des cellules et,
pour le mode de limitation imligp = 1, transmettre également
une variable allongée pour le gradient et les tableaux de travail
denom et denum.
*/
/*
Echanges pour mise à jour du voisinage étendu séparé
----------------------------------------------------
*/
#if defined(_CS_HAVE_MPI)
if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) {
/* Identification des pointeurs sur le voisinage étendu séparé */
if (maillage->nbr_cel_fac_par_voiset > 0) {
pvar_voiset = maillage_tmp->var_cel_avec_voiset
+ maillage->nbr_cel_etendu;
}
else {
pvar_voiset = NULL;
}
/* Pour le gradient et deux tableaux de travail denom et denum.
L'utilisation d'un pointeur _tmp pour le gradient (au lieu d'un
pointeur dpdx_avec_voiset directement) provient du fait qu'on a recopié
la fonction CGRDMC. Ici ce n'était pas indispensable, dpdx en entrée
n'étant pas const.
Pour denom et denum, pas nécessaire d'initialiser les tableaux
on n'a donc pas besoin de denom_voiset et denum_voiset */
if (*imligp == 1) {
if (maillage->nbr_cel_fac_par_voiset > 0) {
nceltot = maillage->nbr_cel_etendu + maillage->nbr_cel_fac_par_voiset;
BFT_MALLOC(dpdx_avec_voiset_tmp, nceltot, cs_real_t);
BFT_MALLOC(dpdy_avec_voiset_tmp, nceltot, cs_real_t);
BFT_MALLOC(dpdz_avec_voiset_tmp, nceltot, cs_real_t);
BFT_MALLOC(denum_avec_voiset_tmp, nceltot, cs_real_t);
BFT_MALLOC(denom_avec_voiset_tmp, nceltot, cs_real_t);
dpdx_voiset = dpdx_avec_voiset_tmp + maillage->nbr_cel_etendu;
dpdy_voiset = dpdy_avec_voiset_tmp + maillage->nbr_cel_etendu;
dpdz_voiset = dpdz_avec_voiset_tmp + maillage->nbr_cel_etendu;
}
else {
dpdx_avec_voiset_tmp = NULL;
dpdy_avec_voiset_tmp = NULL;
dpdz_avec_voiset_tmp = NULL;
denum_avec_voiset_tmp = NULL;
denom_avec_voiset_tmp = NULL;
dpdx_voiset = NULL;
dpdy_voiset = NULL;
dpdz_voiset = NULL;
}
}
/* Echanges parallèles et périodiques
Attention, un processeur peut n'avoir rien à recevoir
(maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant
des éléments à envoyer : on doit donc entrer dans parcve
pour tous les processeurs lorsque l'on est en parallèle */
if (maillage->nbr_dom > 1)
cs_parallel_parcve(pvar, pvar_voiset);
if (maillage->nbr_per > 0)
cs_perio_percve(pvar, pvar_voiset);
/* Echange pour les gradients
(Pour les tableaux de travail, échange inutile) */
if (*imligp == 1) {
if (maillage->nbr_dom > 1) {
cs_parallel_parcve(dpdx, dpdx_voiset);
cs_parallel_parcve(dpdy, dpdy_voiset);
cs_parallel_parcve(dpdz, dpdz_voiset);
}
if (maillage->nbr_per > 0) {
cs_perio_percve(dpdx, dpdx_voiset);
cs_perio_percve(dpdy, dpdy_voiset);
cs_perio_percve(dpdz, dpdz_voiset);
}
}
}
#endif /*_CS_HAVE_MPI */
/*
Appel de la limitation du gradient
----------------------------------
*/
/* Identification du pointeur sur la variable allongée
(de longueur ncelet + voisinage étendu séparé éventuel)
et copie des valeurs 1-ncelet si nécessaire */
#if defined(_CS_HAVE_MPI)
if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
&& maillage->nbr_cel_fac_par_voiset > 0) {
for (iel = 0 ; iel < (*ncelet) ; iel++)
maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel];
pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset;
if (*imligp == 1) {
for (iel = 0 ; iel < (*ncelet) ; iel++)
dpdx_avec_voiset_tmp[iel] = dpdx[iel];
dpdx_avec_voiset = dpdx_avec_voiset_tmp;
for (iel = 0 ; iel < (*ncelet) ; iel++)
dpdy_avec_voiset_tmp[iel] = dpdy[iel];
dpdy_avec_voiset = dpdy_avec_voiset_tmp;
for (iel = 0 ; iel < (*ncelet) ; iel++)
dpdz_avec_voiset_tmp[iel] = dpdz[iel];
dpdz_avec_voiset = dpdz_avec_voiset_tmp;
denum_avec_voiset = denum_avec_voiset_tmp;
denom_avec_voiset = denom_avec_voiset_tmp;
}
else {
dpdx_avec_voiset = dpdx;
dpdy_avec_voiset = dpdy;
dpdz_avec_voiset = dpdz;
denum_avec_voiset = denum;
denom_avec_voiset = denom;
}
}
else {
#endif /*_CS_HAVE_MPI */
pvar_avec_voiset = pvar;
dpdx_avec_voiset = dpdx;
dpdy_avec_voiset = dpdy;
dpdz_avec_voiset = dpdz;
denum_avec_voiset = denum;
denom_avec_voiset = denom;
#if defined(_CS_HAVE_MPI)
}
#endif /*_CS_HAVE_MPI */
/* Identification du pointeur sur les coordonnées allongées
(de longueur ncelet + voisinage étendu séparé éventuel)*/
if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
&& maillage->nbr_cel_fac_par_voiset > 0) {
xyzcen_avec_voiset = maillage->coord_cel_avec_voiset;
}
else
xyzcen_avec_voiset = xyzcen;
/* Pointeurs sur le voisinage étendu */
ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
ielvse = maillage->val_connect_voiset_cel_cel_dom;
/* Appel de la limitation du gradient */
CS_PROCF(limgrd, LIMGRD)
(ndim, ncelet, ncel, nfac,
imrgra, imligp, idimte, itenso, iwarnp, nfecra, climgp,
ifacel, ipcvse, ielvse,
xyzcen_avec_voiset, pvar_avec_voiset,
dpdx_avec_voiset, dpdy_avec_voiset, dpdz_avec_voiset,
faclim, denom_avec_voiset, denum_avec_voiset);
/* Transfert du gradient limité si on a utilisé un tableau différent
puis libération des tableaux locaux
(voisinage étendu séparé, non vide et avec imligp == 1) */
#if defined(_CS_HAVE_MPI)
if ( (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE)
&& (maillage->nbr_cel_fac_par_voiset > 0)
&& (*imligp == 1)) {
for (iel = 0 ; iel < (*ncelet) ; iel++)
dpdx[iel] = dpdx_avec_voiset[iel];
for (iel = 0 ; iel < (*ncelet) ; iel++)
dpdy[iel] = dpdy_avec_voiset[iel];
for (iel = 0 ; iel < (*ncelet) ; iel++)
dpdz[iel] = dpdz_avec_voiset[iel];
BFT_FREE(dpdx_avec_voiset_tmp);
BFT_FREE(dpdy_avec_voiset_tmp);
BFT_FREE(dpdz_avec_voiset_tmp);
BFT_FREE(denum_avec_voiset_tmp);
BFT_FREE(denom_avec_voiset_tmp);
}
#endif /*_CS_HAVE_MPI */
}
/*----------------------------------------------------------------------------
* Appel de la routine FORTRAN CFILTR pour le calcul des filtres avec voisinages
* etendus pour le modele dynamique.
*
* Interface Fortran :
*
* SUBROUTINE CFILTR
* *****************
*
* & ( NCELET , NCEL , NFAC , NFABOR ,
* & IFACEL , IFABOR ,
* & VOLUME ,
* & PVAR , VARFIL ,
* & W1 , W2 )
*
*----------------------------------------------------------------------------*/
void CS_PROCF (cfiltr, CFILTR)
(
const cs_int_t *const ncelet, /* --> nombre de cellules étendu */
const cs_int_t *const ncel, /* --> nombre de cellules */
const cs_int_t *const nfac, /* --> nombre de faces internes */
const cs_int_t *const nfabor, /* --> nombre de faces de bord */
const cs_int_t ifacel[], /* --> liste des faces internes */
const cs_int_t ifabor[], /* --> liste des faces de bord */
const cs_real_t volume[], /* --> volumes des cellules */
const cs_real_t pvar[], /* --> variable a filter */
cs_real_t varfil[], /* --> resultat de la variable filtree */
cs_real_t w1[], /* - tab. de travail local */
cs_real_t w2[] /* - tab. de travail local */
)
{
#if defined(_CS_HAVE_MPI)
cs_int_t iel;
cs_int_t nceltot;
#endif /*_CS_HAVE_MPI */
cs_int_t *ipcvse;
cs_int_t *ielvse;
const cs_real_t *pvar_avec_voiset;
const cs_real_t *volume_avec_voiset;
#if defined(_CS_HAVE_MPI)
cs_real_t *pvar_voiset;
cs_real_t *volume_voiset;
cs_real_t *volume_avec_voiset_tmp;
#endif /*_CS_HAVE_MPI */
cs_maillage_t *maillage = cs_glob_maillage;
#if defined(_CS_HAVE_MPI)
cs_maillage_tmp_t *maillage_tmp = cs_glob_maillage_tmp;
#endif /*_CS_HAVE_MPI */
/*
Remarque sur les tests faisant intervenir _CS_HAVE_MPI
Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie
definie que si on _CS_HAVE_MPI (si le parallélisme est possible)
Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est
correct car il ne peut y avoir un voisinage separé que en parallele
etant donne le choix fait dans l'enveloppe (un voisinage separe
en non parallele n'apporterait rien car il ne permettrait pas de
reduire le volume des echanges).
On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI
n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t
par des tests sur _CS_HAVE_MPI pour permettre la compilation
lorsque _CS_HAVE_MPI n'est pas defini
Par coherence, on fait de meme avec les variables pvar_voiset
nceltot, iel
*/
/* La complexite du traitement dans cette fonction est liée à la prise en
compte d'un voisinage étendu séparé du halo classique */
/* On va transmettre au gradient par moindre carrés une variable
pvar_avec_voiset qui sera :
maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe
et qu'il est traité séparément
pvar sinon (cas classiques)
La variable maillage_tmp->var_cel_avec_voiset devra contenir
pvar dans les ncelet premieres cases et les valeurs actualisees
du voisinage étendu séparé dans le reste, noté pvar_voiset
Il faut de même transmettre la variable allongée correcte pour les
centres de gravité des cellules et,
lorsque la pression hydrostatique est prise en compte, transmettre
également des variables allongées pour la fext.
L'utilisation des pointeurs fext?_avec_voiset_tmp au lieu de
fext?_avec_voiset directement permet d'éviter un warning et de conserver
le caractere const des fext? en argument.
*/
/*
Echanges pour mise à jour du voisinage étendu séparé
----------------------------------------------------
*/
#if defined(_CS_HAVE_MPI)
if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) {
/* Identification des pointeurs sur le voisinage étendu séparé */
if (maillage->nbr_cel_fac_par_voiset > 0) {
pvar_voiset = maillage_tmp->var_cel_avec_voiset + (*ncelet);
}
else {
pvar_voiset = NULL;
}
if (maillage->nbr_cel_fac_par_voiset > 0) {
nceltot = (*ncelet)+maillage->nbr_cel_fac_par_voiset;
BFT_MALLOC(volume_avec_voiset_tmp,nceltot,cs_real_t);
volume_voiset = volume_avec_voiset_tmp + (*ncelet);
}
else {
volume_avec_voiset_tmp = NULL;
volume_voiset = NULL;
}
/* Echanges parallèles et périodiques
Attention, un processeur peut n'avoir rien à recevoir
(maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant
des éléments à envoyer : on doit donc entrer dans parcve
pour tous les processeurs lorsque l'on est en parallèle */
if (maillage->nbr_dom > 1)
cs_parallel_parcve(pvar, pvar_voiset);
cs_parallel_parcve(volume, volume_voiset);
if (maillage->nbr_per > 0)
cs_perio_percve (pvar, pvar_voiset);
cs_perio_percve (volume, volume_voiset);
}
#endif /*_CS_HAVE_MPI */
#if defined(_CS_HAVE_MPI)
if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
&& maillage->nbr_cel_fac_par_voiset > 0) {
for (iel = 0 ; iel < (*ncelet) ; iel++)
maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel];
pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset;
for (iel = 0 ; iel < (*ncelet) ; iel++)
volume_avec_voiset_tmp[iel] = volume[iel];
volume_avec_voiset = volume_avec_voiset_tmp;
}
else{
#endif /*_CS_HAVE_MPI */
pvar_avec_voiset = pvar;
volume_avec_voiset = volume;
#if defined(_CS_HAVE_MPI)
}
#endif /*_CS_HAVE_MPI */
/* Pointeurs sur le voisinage étendu */
ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
ielvse = maillage->val_connect_voiset_cel_cel_dom;
/* Identification du pointeur sur la variable allongée
(de longueur ncelet + voisinage étendu séparé éventuel)
et copie des valeurs 1-ncelet si nécessaire */
CS_PROCF(filtre, FILTRE)
(ncelet, ncel , nfac , nfabor,
ifacel, ifabor, ipcvse, ielvse,
volume_avec_voiset,
pvar_avec_voiset, varfil,
w1 , w2 );
/* Libération si on a réservé quelque chose */
#if defined(_CS_HAVE_MPI)
if ( (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE)
&& (maillage->nbr_cel_fac_par_voiset > 0)) {
BFT_FREE(volume_avec_voiset_tmp);
}
#endif /*_CS_HAVE_MPI */
}
/*----------------------------------------------------------------------------*/
#ifdef __cplusplus
}
#endif /* __cplusplus */
syntax highlighted by Code2HTML, v. 0.9.1