rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/herezh_pp/Elements/Mecanique/SFE/Met_Sfe1s3.cc

1439 lines
69 KiB
C++
Executable file
Raw Permalink Blame History

// FICHIER : Met_Sfe1s3.cc
// CLASSE : Met_Sfe1
// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2021 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "Met_Sfe1.h"
// =========================================================================================
// vu la taille des executables le fichier est decompose en quatre
// le premier : Met_Sfe1s1.cp concerne les constructeurs, destructeur et
// la gestion des variables
// le second : Met_Sfe1s2.cp concerne le calcul des grandeurs publics
// le troisieme : Met_Sfe1s3.cp concerne le calcul des grandeurs protegees, et des courbures pour SFE1 et SFE2 et QSFE1
// le quatrième : Met_Sfe1s4.cp concerne le calcul des courbures pour les éléments SFE3 et QSFE3
// =========================================================================================
//==================== METHODES PROTEGEES===============================
// calcul des normales a la facette
void Met_Sfe1::Calcul_N_0 ()
{ N_0 = (Util::ProdVec_coorB((*aiB_0)(1),(*aiB_0)(2))).Normer();};
void Met_Sfe1::Calcul_N_t ()
{ N_t = (Util::ProdVec_coorB((*aiB_t)(1),(*aiB_t)(2)));
norme_N_t = N_t.Norme();
N_t /= norme_N_t;
};
void Met_Sfe1::Calcul_N_tdt ()
{ N_tdt = (Util::ProdVec_coorB((*aiB_tdt)(1),(*aiB_tdt)(2)));
norme_N_tdt = N_tdt.Norme();
N_tdt /= norme_N_tdt;
};
//== calcul des variations de la normales
// il faut avant le calcul des ap_alpha_t et de épaisseur à t
void Met_Sfe1::Cal_d_N_t()
{ int nbddlx = tab_ddl.NbDdl_famille(X1);
CoordonneeB d_NN; // un vecteur de travail
for (int iddl=1;iddl<= nbddlx;iddl++)
{ // tout d'abord on calcul la variation de la normale non normée
d_NN = (Util::ProdVec_coorB((*d_aiB_t)(iddl)(1),(*aiB_t)(2)))
+ (Util::ProdVec_coorB((*aiB_t)(1),(*d_aiB_t)(iddl)(2)));
// maintenant on calcul la variation du vecteur normal unitaire
(*d_N_t)(iddl) = Util::VarUnVect_coorB(N_t,d_NN,norme_N_t);
};
};
// il faut avant le calcul des ap_alpha_tdt et de épaisseur à tdt
void Met_Sfe1::Cal_d_N_tdt()
{ int nbddlx = tab_ddl.NbDdl_famille(X1);
CoordonneeB d_NN; // un vecteur de travail
for (int iddl=1;iddl<= nbddlx;iddl++)
{ // tout d'abord on calcul la variation de la normale non normée
d_NN = (Util::ProdVec_coorB((*d_aiB_tdt)(iddl)(1),(*aiB_tdt)(2)))
+ (Util::ProdVec_coorB((*aiB_tdt)(1),(*d_aiB_tdt)(iddl)(2)));
// maintenant on calcul la variation du vecteur normal unitaire
(*d_N_tdt)(iddl) = Util::VarUnVect_coorB(N_tdt,d_NN,norme_N_tdt);
};
};
// ---- calcul des points
// calcul du point a t0
void Met_Sfe1::Calcul_M0( const Tableau<Noeud *>& , const Vecteur& phiH, int )
{ // on considère qu'il y a deux points dans l'épaisseur, un a -h/2 et l'autre à h/2
*M0 = *P0 + (N_0 * (epais.epaisseur0 * 0.5 * (phiH(2)-phiH(1)))).Coor_const();};
void Met_Sfe1::Calcul_Mt( const Tableau<Noeud *>& , const Vecteur& phiH, int )
{ // on considère qu'il y a deux points dans l'épaisseur, un a -h/2 et l'autre à h/2
*Mt = *Pt + (N_t * (epais.epaisseur_t * 0.5 * (phiH(2)-phiH(1)))).Coor_const();};
void Met_Sfe1::Calcul_Mtdt( const Tableau<Noeud *>& , const Vecteur& phiH, int )
{ // on considère qu'il y a deux points dans l'épaisseur, un a -h/2 et l'autre à h/2
*Mtdt = *Ptdt + (N_tdt * (epais.epaisseur_tdt * 0.5 * (phiH(2)-phiH(1)))).Coor_const();};
// cas de l'épaisseur interpolée à 0
void Met_Sfe1::Calcul_epaisseur_0(const Tableau<Noeud *>& tab_noeud, const Vecteur& phiS)
{ double & epaisseur = epais.epaisseur0;
epaisseur = 0.; // init
for (int r=1;r<=nbNoeudCentral;r++)
epaisseur += tab_noeud(r)->Valeur_0(EPAIS) * phiS(r);
};
// cas de l'épaisseur interpolée à t
void Met_Sfe1::Calcul_epaisseur_t(const Tableau<Noeud *>& tab_noeud, const Vecteur& phiS)
{ double & epaisseur = epais.epaisseur_t;
epaisseur = 0.; // init
for (int r=1;r<=nbNoeudCentral;r++)
epaisseur += tab_noeud(r)->Valeur_t(EPAIS) * phiS(r);
};
// cas de l'épaisseur interpolée à tdt
void Met_Sfe1::Calcul_epaisseur_tdt(const Tableau<Noeud *>& tab_noeud, const Vecteur& phiS)
{ double & epaisseur = epais.epaisseur_tdt;
epaisseur = 0.; // init
for (int r=1;r<=nbNoeudCentral;r++)
epaisseur += tab_noeud(r)->Valeur_tdt(EPAIS) * phiS(r);
};
// calcul de la variation de l'épaisseur par rapport aux ddl
void Met_Sfe1::D_epaisseur_t(const Vecteur& phiS)
{ // l'épaisseur ne dépend que des ddl d'épaisseur
for (int r=1;r<=nbNoeudCentral;r++)
(*d_epais)(r).epaisseur_t = phiS(r);
};
// calcul de la variation de l'épaisseur par rapport aux ddl
void Met_Sfe1::D_epaisseur_tdt(const Vecteur& phiS)
{ // l'épaisseur ne dépend que des ddl d'épaisseur
for (int r=1;r<=nbNoeudCentral;r++)
(*d_epais)(r).epaisseur_tdt = phiS(r);
};
//== calcul de la variation de N,alpha
void Met_Sfe1::Cal_d_N_alpha_t()
{ int nbddlx = tab_ddl.NbDdl_famille(X1);
for (int iddl=1;iddl<= nbddlx;iddl++)
{ // variation de la derivee du vecteur normal
(*d_N_alpha_t_1)(iddl).ConstructionAPartirDe_H(
-(dcurb_t(iddl)(1,1) * (*aiH_t)(1) + curb_t(1,1) * (*d_aiH_t)(iddl)(1)
+dcurb_t(iddl)(1,2) * (*aiH_t)(2) + curb_t(1,2) * (*d_aiH_t)(iddl)(2)));
(*d_N_alpha_t_2)(iddl).ConstructionAPartirDe_H(
-(dcurb_t(iddl)(2,1) * (*aiH_t)(1) + curb_t(2,1) * (*d_aiH_t)(iddl)(1)
+ dcurb_t(iddl)(2,2) * (*aiH_t)(2) + curb_t(2,2) * (*d_aiH_t)(iddl)(2)));
};
};
void Met_Sfe1::Cal_d_N_alpha_tdt()
{ int nbddlx = tab_ddl.NbDdl_famille(X1);
for (int iddl=1;iddl<= nbddlx;iddl++)
{ // variation de la derivee du vecteur normal
(*d_N_alpha_tdt_1)(iddl).ConstructionAPartirDe_H(
-(dcurb_tdt(iddl)(1,1) * (*aiH_tdt)(1) + curb_tdt(1,1) * (*d_aiH_tdt)(iddl)(1)
+dcurb_tdt(iddl)(1,2) * (*aiH_tdt)(2) + curb_tdt(1,2) * (*d_aiH_tdt)(iddl)(2)));
(*d_N_alpha_tdt_2)(iddl).ConstructionAPartirDe_H(
-(dcurb_tdt(iddl)(2,1) * (*aiH_tdt)(1) + curb_tdt(2,1) * (*d_aiH_tdt)(iddl)(1)
+ dcurb_tdt(iddl)(2,2) * (*aiH_tdt)(2) + curb_tdt(2,2) * (*d_aiH_tdt)(iddl)(2)));
};
};
// ------ calcul des bases pour la surface médiane -------
// calcul de la base naturel a t0 pour la surface médiane: partie a_alpha
void Met_Sfe1::Calcul_a_alpha_B_0
( const Tableau<Noeud *>& tab_noeud, const Mat_pleine& dphi, int ,const Vecteur&)
{
#ifdef MISE_AU_POINT
if (aiB_0 == NULL)
{ cout << "\nErreur : la base a t=0 n'est pas dimensionne !\n";
cout << "void Met_Sfe1::Calcul_aiB_0 \n";
Sortie(1);
};
#endif
for (int i=1;i<= 2;i++)
{for (int a=1;a<= dim_base;a++)
{ double & aib0_i_a = aiB_0->CoordoB(i)(a);
aib0_i_a = 0.;
for (int r=1;r<=nbNoeudCentral;r++)
{ aib0_i_a += tab_noeud(r)->Coord0()(a) * dphi(i,r);
};
};
};
};
// calcul de la base naturel a t pour la surface médiane: partie a_alpha
void Met_Sfe1::Calcul_a_alpha_B_t
( const Tableau<Noeud *>& tab_noeud, const Mat_pleine& dphi, int ,const Vecteur&)
{
#ifdef MISE_AU_POINT
if (aiB_t == NULL)
{ cout << "\nErreur : la base a t n'est pas dimensionne !\n";
cout << "void Met_Sfe1::Calcul_aiB_t \n";
Sortie(1);
};
#endif
for (int i=1;i<= 2;i++)
{for (int a=1;a<= dim_base;a++)
{ double & aibt_i_a = aiB_t->CoordoB(i)(a);
aibt_i_a = 0.;
for (int r=1;r<=nbNoeudCentral;r++)
{ aibt_i_a += tab_noeud(r)->Coord1()(a) * dphi(i,r);
};
};
};
};
// calcul de la base naturel a tdt pour la surface médiane: partie a_alpha
void Met_Sfe1::Calcul_a_alpha_B_tdt
( const Tableau<Noeud *>& tab_noeud, const Mat_pleine& dphi, int ,const Vecteur&)
{
#ifdef MISE_AU_POINT
if (aiB_tdt == NULL)
{ cout << "\nErreur : la base a t+dt n'est pas dimensionne !\n";
cout << "void Met_Sfe1::Calcul_aiB_tdt \n";
Sortie(1);
};
#endif
for (int i=1;i<= 2;i++)
{for (int a=1;a<= dim_base;a++)
{ double & aibtdt_i_a = aiB_tdt->CoordoB(i)(a);
aibtdt_i_a = 0.;
for (int r=1;r<=nbNoeudCentral;r++)
{ aibtdt_i_a += tab_noeud(r)->Coord2()(a) * dphi(i,r);
};
};
};
};
//== calcul de la variation des bases pour la surface médiane: partie a_alpha
// ici il ne s'agit uniquement que de la variation de a_alpha et pas a_3
void Met_Sfe1::D_a_alpha_B_t( const Mat_pleine& dphi, int ,const Vecteur & )
{ int indice;
// derivees de aBi par rapport a XHbr
for (int b = 1; b<= dim_base; b++)
for (int r = 1; r <= nbNoeudCentral; r++)
{ indice = (r-1)*dim_base + b;
BaseB & d_aiB_t_indice=(*d_aiB_t)(indice);
for (int i =1;i<=2;i++) // alpha = 1 et 2 uniquement
{ d_aiB_t_indice.CoordoB(i).Zero();
d_aiB_t_indice.CoordoB(i)(b)=dphi(i,r);
}
};
};
// ici il ne s'agit uniquement que de la variation de a_alpha et pas a_3
void Met_Sfe1::D_a_alpha_B_tdt( const Mat_pleine& dphi, int ,const Vecteur &)
{ int indice;
for (int b = 1; b<= dim_base; b++)
for (int r = 1; r <= nbNoeudCentral; r++)
{ indice = (r-1)*dim_base + b;
BaseB & d_aiB_tdt_indice=(*d_aiB_tdt)(indice);
for (int i =1;i<=2;i++) // alpha = 1 et 2 uniquement
{ d_aiB_tdt_indice.CoordoB(i).Zero();
d_aiB_tdt_indice.CoordoB(i)(b)=dphi(i,r);
}
}
};
//== calcul de la variation de a_3_B
void Met_Sfe1::D_a_3_B_t()
{ int nbddlx = tab_ddl.NbDdl_famille(X1);
// calcul de la variation de a3B par rapport au ddlxi
for (int iddl=1;iddl<= nbddlx;iddl++)
(*d_aiB_t)(iddl).CoordoB(3) = (epais.epaisseur_t * 0.5) * (*d_N_t)(iddl);
// calcul de la variation en fonction des ddl d'épaisseur
int nbddl_total = nbddlx + nbNoeudCentral;int ih=1;
for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
(*d_aiB_t)(iddl).CoordoB(3) = ((*d_epais)(ih).epaisseur_t * 0.5 ) * N_t;
};
void Met_Sfe1::D_a_3_B_tdt()
{ int nbddlx = tab_ddl.NbDdl_famille(X1);
// calcul de la variation de a3B par rapport au ddlxi
for (int iddl=1;iddl<= nbddlx;iddl++)
(*d_aiB_tdt)(iddl).CoordoB(3) = (epais.epaisseur_tdt * 0.5) * (*d_N_tdt)(iddl);
// calcul de la variation en fonction des ddl d'épaisseur
int nbddl_total = nbddlx + nbNoeudCentral;int ih=1;
for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
(*d_aiB_tdt)(iddl).CoordoB(3) = ((*d_epais)(ih).epaisseur_tdt * 0.5 ) * N_tdt;
};
// --------- calcul de la base naturel a t0 ------
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de la facette
// également l'épaisseur interpolée si on est avec des éléments 3D
void Met_Sfe1::Calcul_giB_0
( const Tableau<Noeud *>& , const Mat_pleine& , int , const Vecteur& phiH)
{
#ifdef MISE_AU_POINT
if (giB_0 == NULL)
{ cout << "\nErreur : la base a t=0 n'est pas dimensionne !\n";
cout << "void Met_Sfe1::Calcul_giB_0 \n";
Sortie(1);
};
#endif
// vecteur de base gi,
// theta^3 = (phiH(2)-phiH(1));
double z0 = epais.epaisseur0 * 0.5 * (phiH(2)-phiH(1));
// ici il n'y a pas le terme de variation de h par rapport au theta^alpha
giB_0->CoordoB(1) = (*aiB_0)(1) + N_alpha_0_1 * z0;
giB_0->CoordoB(2) = (*aiB_0)(2) + N_alpha_0_2 * z0;
if (tCalEpais) // cas 3D
giB_0->CoordoB(3) = (*aiB_0)(3);
};
// calcul de la base naturel a t
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de la facette
void Met_Sfe1::Calcul_giB_t
( const Tableau<Noeud *>& , const Mat_pleine& , int, const Vecteur& phiH)
{
#ifdef MISE_AU_POINT
if (giB_t == NULL)
{ cout << "\nErreur : la base a t n'est pas dimensionne !\n";
cout << "void Met_Sfe1::Calcul_giB_t \n";
Sortie(1);
};
#endif
// vecteur de base gi
double z_t = epais.epaisseur_t * 0.5 * (phiH(2)-phiH(1));
// ici il n'y a pas le terme de variation de h par rapport au theta^alpha
giB_t->CoordoB(1) = (*aiB_t)(1) + N_alpha_t_1 * z_t;
giB_t->CoordoB(2) = (*aiB_t)(2) + N_alpha_t_2 * z_t;
if (tCalEpais) // cas 3D
giB_t->CoordoB(3) = (*aiB_t)(3);
};
// calcul de la base naturel a tdt
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de la facette
void Met_Sfe1::Calcul_giB_tdt
( const Tableau<Noeud *>& , const Mat_pleine& , int, const Vecteur& phiH)
{
#ifdef MISE_AU_POINT
if (giB_tdt == NULL)
{ cout << "\nErreur : la base a t+dt n'est pas dimensionne !\n";
cout << "void Met_Sfe1::Calcul_giB_tdt \n";
Sortie(1);
};
#endif
// vecteur de base gi
double z_tdt = epais.epaisseur_tdt * 0.5 * (phiH(2)-phiH(1));
// ici il n'y a pas le terme de variation de h par rapport au theta^alpha
giB_tdt->CoordoB(1) = (*aiB_tdt)(1) + N_alpha_tdt_1 * z_tdt;
giB_tdt->CoordoB(2) = (*aiB_tdt)(2) + N_alpha_tdt_2 * z_tdt;
if (tCalEpais) // cas 3D
giB_tdt->CoordoB(3) = (*aiB_tdt)(3);
};
//------------// variation des vecteurs de base
// il faut au par avant avoir calculé la variation de l'épaisseur si on est en 3D
void Met_Sfe1::D_giB_t( const Mat_pleine& , int , const Vecteur & phiH)
{ double z_t = epais.epaisseur_t * 0.5 * (phiH(2)-phiH(1));
int nbddlx = tab_ddl.NbDdl_famille(X1);
for (int iddl=1;iddl<= nbddlx;iddl++)
{ // vecteur de base gi
(*d_giB_t)(iddl).CoordoB(1) = (*d_aiB_t)(iddl)(1) + (*d_N_alpha_t_1)(iddl) * z_t;
(*d_giB_t)(iddl).CoordoB(2) = (*d_aiB_t)(iddl)(2) + (*d_N_alpha_t_2)(iddl) * z_t;
};
// cas 3D
if (tCalEpais)
{// variation de giB3 = en fait celle de aiB3
int nbddl_total = nbddlx + nbNoeudCentral;
for (int iddl=1;iddl<= nbddl_total;iddl++)
(*d_giB_t)(iddl).CoordoB(3) = (*d_aiB_t)(iddl)(3);
// ensuite les vecteurs g1B et g2B dépendent aussi des ddl d'épaisseur
// on calcul donc cette variation, et on en profite pour calculer également la variation de g3B
// en fonction des ddl d'épaisseur
int ih=1;
for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
{(*d_giB_t)(iddl).CoordoB(1) = ((*d_epais)(ih).epaisseur_t * 0.5 * (phiH(2)-phiH(1))) * N_alpha_t_1;
(*d_giB_t)(iddl).CoordoB(2) = ((*d_epais)(ih).epaisseur_t * 0.5 * (phiH(2)-phiH(1))) * N_alpha_t_2;
};
};
};
void Met_Sfe1::D_giB_tdt( const Mat_pleine& , int , const Vecteur & phiH)
{ double z_tdt = epais.epaisseur_tdt * 0.5 * (phiH(2)-phiH(1));
int nbddlx = tab_ddl.NbDdl_famille(X1);
for (int iddl=1;iddl<= nbddlx;iddl++)
{// vecteur de base gi
(*d_giB_tdt)(iddl).CoordoB(1) = (*d_aiB_tdt)(iddl)(1) + (*d_N_alpha_tdt_1)(iddl) * z_tdt;
(*d_giB_tdt)(iddl).CoordoB(2) = (*d_aiB_tdt)(iddl)(2) + (*d_N_alpha_tdt_2)(iddl) * z_tdt;
#ifdef MISE_AU_POINT
if (ParaGlob::NiveauImpression() > 5)
{if (((*d_giB_tdt)(iddl)(1).Norme() > 10000.) || ((*d_giB_tdt)(iddl)(1).Norme() < -10000.)
|| ((*d_giB_tdt)(iddl)(2).Norme() > 10000.) || ((*d_giB_tdt)(iddl)(2).Norme() < -10000.) )
{cout << "\n attention erreur , les d_giB_tdt sont trop grand "
<< "\n Met_Sfe1::D_giB_tdt(..."
<< " (*d_giB_tdt)(iddl)(1).Norme()=" << (*d_giB_tdt)(iddl)(1).Norme()
<< " (*d_giB_tdt)(iddl)(2).Norme()=" << (*d_giB_tdt)(iddl)(2).Norme()
<< " (*d_aiB_tdt)(iddl)(1).Norme()=" << (*d_aiB_tdt)(iddl)(1).Norme()
<< " (*d_aiB_tdt)(iddl)(2).Norme()=" << (*d_aiB_tdt)(iddl)(2).Norme()
<< " (*d_N_alpha_tdt_1)(iddl).Norme()=" << (*d_N_alpha_tdt_1)(iddl).Norme()
<< " (*d_N_alpha_tdt_2)(iddl).Norme()=" << (*d_N_alpha_tdt_2)(iddl).Norme()
<< endl ;
};
};
#endif
};
// cas 3D
if (tCalEpais)
{// variation de giB3 = en fait celle de aiB3
int nbddl_total = nbddlx + nbNoeudCentral;
for (int iddl=1;iddl<= nbddl_total;iddl++)
(*d_giB_tdt)(iddl).CoordoB(3) = (*d_aiB_tdt)(iddl)(3);
// ensuite les vecteurs g1B et g2B dépendent aussi des ddl d'épaisseur
// on calcul donc cette variation, et on en profite pour calculer également la variation de g3B
// en fonction des ddl d'épaisseur
int ih=1;
for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
{(*d_giB_tdt)(iddl).CoordoB(1) = ((*d_epais)(ih).epaisseur_tdt * 0.5 * (phiH(2)-phiH(1))) * N_alpha_tdt_1;
(*d_giB_tdt)(iddl).CoordoB(2) = ((*d_epais)(ih).epaisseur_tdt * 0.5 * (phiH(2)-phiH(1))) * N_alpha_tdt_2;
};
};
};
//-----------// calcul du tenseur de courbure dans la base naturelle
// on utilise un tenseur non symétrique, bien qu'en théorie la courbure est symétrique
// suivant les différentes méthodes de calcul, le tenseur en sortie est symétrique ou pas
// on peut toujours le symétriser, si l'on veut un tenseur symétrique
// on a : N,al = - b_{al be} . a^{be}
// plusieurs cas sont etudies suivant l'instant considere
// a l'instant t = 0
Tenseur_ns2BB& Met_Sfe1::courbure_0 ( const Tableau<Noeud *>& tab_noeud
,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ Tableau<Coordonnee3> tab_coor(6);
// actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
// du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
tab_coor(i) = tab_noeud(i)->Coord0();
for (int i=1;i< 4;i++) // les noeuds externes
tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord0();
// for (int i=1;i<=6;i++)
// tab_coor(i) = tab_noeud(i)->Coord0();
curb_0 = (this->*PtCourbure)(tab_coor,(*aiB_0),(*aiH_0),tabTypeCL,vplan);
return curb_0;
};
// a l'instant t
Tenseur_ns2BB& Met_Sfe1::courbure_t ( const Tableau<Noeud *>& tab_noeud
,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ Tableau<Coordonnee3> tab_coor(6);
// actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
// du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
tab_coor(i) = tab_noeud(i)->Coord1();
for (int i=1;i< 4;i++) // les noeuds externes
tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord1();
// for (int i=1;i<=6;i++)
// tab_coor(i) = tab_noeud(i)->Coord1();
curb_t = (this->*PtCourbure)(tab_coor,(*aiB_t),(*aiH_t),tabTypeCL,vplan);
return curb_t;
};
// a l'instant t+dt
Tenseur_ns2BB& Met_Sfe1::courbure_tdt ( const Tableau<Noeud *>& tab_noeud
,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ Tableau<Coordonnee3> tab_coor(6);
// actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
// du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
tab_coor(i) = tab_noeud(i)->Coord2();
for (int i=1;i< 4;i++) // les noeuds externes
tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord2();
// for (int i=1;i<=6;i++)
// tab_coor(i) = tab_noeud(i)->Coord2();
curb_tdt = (this->*PtCourbure)(tab_coor,(*aiB_tdt),(*aiH_tdt),tabTypeCL,vplan);
return curb_tdt;
};
//--------------// calcul du tenseur de courbure et de sa variation
// plusieurs cas sont etudies suivant l'instant considere
// a l'instant t
void Met_Sfe1::Dcourbure_t
( const Tableau<Noeud *>& tab_noeud,Tenseur_ns2BB& curb,TabOper<Tenseur_ns2BB>& dcurb
,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ Tableau<Coordonnee3> tab_coor(6);
// actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
// du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
tab_coor(i) = tab_noeud(i)->Coord1();
for (int i=1;i< 4;i++) // les noeuds externes
tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord1();
// for (int i=1;i<=6;i++)
// tab_coor(i) = tab_noeud(i)->Coord1();
(this->*PtDcourbure) (tab_coor,(*aiB_t),*d_aiB_t,(*aiH_t),*d_aiH_t,curb,dcurb,tabTypeCL,vplan);
};
// a l'instant t+dt
void Met_Sfe1::Dcourbure_tdt
( const Tableau<Noeud *>& tab_noeud,Tenseur_ns2BB& curb,TabOper<Tenseur_ns2BB>& dcurb
,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ Tableau<Coordonnee3> tab_coor(6);
// actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
// du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
tab_coor(i) = tab_noeud(i)->Coord2();
for (int i=1;i< 4;i++) // les noeuds externes
tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord2();
// for (int i=1;i<=6;i++)
// tab_coor(i) = tab_noeud(i)->Coord2();
(this->*PtDcourbure) (tab_coor,(*aiB_tdt),*d_aiB_tdt,(*aiH_tdt),*d_aiH_tdt,curb,dcurb,tabTypeCL,vplan);
};
// ----------- routines internes pour calculer les différentes formes de courbures
// routine generale de calcul de la courbure
// cas 1: calcul des 3 courbures suivant les 3 cotés
Tenseur_ns2BB Met_Sfe1::Courbure1
( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH & aiH
,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
{ Vecteur curbp(3),curbi(3); // : les tenseurs de courbures d'abord dans les axes
// locaux des arretes du triangle puis dans le repere naturel: b11,b12,b22
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
// calcul de la normale centrale
Coordonnee3 N(Util::ProdVec_coorBN(aiB(1),aiB(2)).Normer());
// calcul du centre de gravite
Coordonnee3 G((tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.);
// le meme calcul est effectue pour les trois cotes du triangle principale
//1- on utilise indi un tableau d'indice pour ne pas etre embete par les modulos 3
// 2- on boucle sur le cotes du triangle
Coordonnee3 Ni(3); // normale du triangle exterieur
// def des différentes grandeurs
Coordonnee3 G2;
Coordonnee3 H2G2,CB,U1,CG2,C_G,CA,dNdV;
Coordonnee3 U2p;
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// on cherche a calculer la longueur H2G2
// simplification de l'ecriture (vérifiée)
const Coordonnee3& D = tab_coor(indi(ncot));
const Coordonnee3& B = tab_coor(indi(ncot+1));
const Coordonnee3& C = tab_coor(indi(ncot+2));
const Coordonnee3& A = tab_coor(indi(ncot)+3);
if (D != A) // cas ou le noeud externe existe
{ // calcul de G2, centre de gravite du triangle exterieur
G2 = (B + C + A) / 3.;
// calcul de U1
CB = B - C;
U1 = CB / CB.Norme();
// vecteur CG2
CG2 = G2 - C;
// longueur h2G2 et g_h1 et calcul du vecteur H2G2 norme
// tout d'abord H2G2 ce situe dans le plan de la facette externe
H2G2 = CG2 - (CG2 * U1) * U1;
double h2G2 = H2G2.Norme();
double unSurh2G2= 1./h2G2;
// maintenant on considere la direction H2G2, ramenée dans le plan de la facette centrale: U2p
U2p = Util::ProdVec_coor(N,U1);
C_G = G - C; // on met un _ car ce nom existe pour une autre routine
double g_h1 = sqrt ( C_G * C_G - ( C_G * U1)* ( C_G * U1)) ;
// calcul de la normale du triangle exterieur
CA = A - C;
Ni = (Util::ProdVec_coor(CB,CA)).Normer();
// calcul de la derivee de la normale
dNdV = ( Ni - N) / ( g_h1 + h2G2);
// courbure dans la direction U2p
curbp(ncot) = - dNdV * U2p;
}
else // cas ou il n'y a pas de noeud externe
{ // la courbure est nulle
curbp(ncot) = 0.;
// la direction de la courbure nulle est normale au cote
// c-a-d le produit vectoriel de la normale avec le cote
CB = B - C;
double norCB = CB.Norme();
U1 = CB / norCB;
U2p = Util::ProdVec_coor(N,U1);
}
// calcul de la matrice de passage de Ualpha a aalpha
double Beta21 = U2p * aiH.Coordo(1);
double Beta22 = U2p * aiH.Coordo(2);
// remplissage de la matrice mat
mat(ncot,1) = Beta21*Beta21;
mat(ncot,2) = 2.* Beta21*Beta22;
mat(ncot,3) = Beta22*Beta22;
};
// calcul du tenseur de courbure dans la base ai
Mat_pleine mati = mat.Inverse();
curbi = mati * curbp;
// retour de la courbure
Tenseur_ns2BB curb;
curb.Coor(1,1)=curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(2); curb.Coor(2,2)=curbi(3);
return curb;
};
// routine generale de calcul de la courbure et de sa variation
// en sortie : curb , la courbure b11,b12,b22
// dcurb , la variation de courbure.
// On ne considere dans les variations que les termes relatif aux ddl qui font varier
// les normales
// cas 1: calcul des 3 courbures suivant les 3 cotés
void Met_Sfe1::Dcourbure1
( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB
, const Tableau <BaseB>& DaiB,const BaseH & aiH
, const Tableau <BaseH>& DaiH,Tenseur_ns2BB& curb,TabOper <Tenseur_ns2BB>& dcurb
,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
{ int nbddl = (nbNoeudCentral+3)*3; //18; // nombre de ddl total
Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations
Vecteur curbi(3); // la courbure sous forme de vecteur
Vecteur curbp(3); // : les tenseurs de courbures dans les axes
// locaux des arretes du triangle
double untiers=1./3.;
Vecteur vecNul(3); // un vecteur nul pour l'initialisation
Tableau <Vecteur> dcurbp(nbddl,vecNul); // variation des courbures
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
Tableau <Mat_pleine> dmat(nbddl,mat);
// calcul de la normale centrale
Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
double nor = NN.Norme(); // la norme
Coordonnee3 N = NN / nor;
Coordonnee3 DNN; // vecteur intermediaire
// calcul du centre de gravite
Coordonnee3 G ((tab_coor(1) + tab_coor(2) + tab_coor(3)) * untiers);
// et le tableau de ses variations
Tableau <Coordonnee3> dG(nbddl,Nul);
// calcul de la variation de la normale centrale et du centre de gravité: vérifiée
Tableau <Coordonnee3> dN(nbddl,Nul);
// on s'intéresse ici qu'aux ddl de la face centrale
for (int inc=1;inc<=3;inc++) // indice du numéro de noeud de la face centrale
for (int ib=1;ib<=3;ib++)
{ // cas de la normale
int i1=(inc-1)*3+ib;
DNN = Util::ProdVec_coorBN(DaiB(i1)(1),aiB(2)) + Util::ProdVec_coorBN(aiB(1),DaiB(i1)(2));
dN(i1) = Util::VarUnVect_coor(NN,DNN,nor);
// cas du centre de gravité
dG(i1)=untiers*Ia(ib);
};
// le meme calcul est effectue pour les trois cotes du triangle principale
//1- on utilise un tableau d'indice pour ne pas etre embete par les modulos 3
// 2- on boucle sur le cotes du triangle
Coordonnee3 NNi,Ni; // normale du triangle exterieur
// -------- declaration pour la boucle sur les 3 cotes ------------
Coordonnee3 G2(3); // init a zero et a 3 composantes
Tableau <Coordonnee3> dB(nbddl,Nul),dC(nbddl,Nul),dA(nbddl,Nul),dG2(nbddl,Nul) ; // variation des points
Vecteur d_g_h1(nbddl);
Coordonnee3 CB,U1,CG2,H2G2,U2p;
Tableau <Coordonnee3> dCB(nbddl,Nul),dU1(nbddl,Nul);
Coordonnee3 C_G,CA,dNdV;
Tableau <Coordonnee3> dCA(nbddl,Nul),dNi(nbddl,Nul);
Coordonnee3 ddNdv;
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// on cherche a calculer la longueur H2G2
// simplification de l'ecriture
const Coordonnee3& D = tab_coor(indi(ncot));
const Coordonnee3& B = tab_coor(indi(ncot+1));
const Coordonnee3& C = tab_coor(indi(ncot+2));
const Coordonnee3& A = tab_coor(indi(ncot)+3);
if (D != A) // cas ou le noeud externe existe
{// calcul de G2, centre de gravite du triangle exterieur
G2 = (B + C + A) * untiers;
// calcul de U1
CB = B - C;
double norCB = CB.Norme();
U1 = CB / norCB;
// vecteur CG2
CG2 = G2 - C;
// longueur h2G2 et g_h1 et calcul du vecteur H2G2 norme
// tout d'abord H2G2 ce citue dans le plan de la facette externe
H2G2 = CG2 - (CG2 * U1) * U1;
double h2G2 = H2G2.Norme();
double unSurh2G2= 1./h2G2;
// maintenant on considere la direction H2G2, ramenée dans le plan de la facette centrale: U2p
U2p = Util::ProdVec_coor(N,U1);
C_G = G - C;
double g_h1 = sqrt ( C_G * C_G - ( C_G * U1) * ( C_G * U1) ) ;
// calcul de la normale du triangle exterieur
CA = A - C;
NNi = Util::ProdVec_coor(CB,CA);
double norNi = NNi.Norme();
Ni = NNi/norNi;
// calcul de la derivee de la normale
dNdV = ( Ni - N) / ( g_h1 + h2G2);
// courbure dans la direction U2p
curbp(ncot) = - dNdV * U2p;
// calcul de la matrice de passage de Ualpha a aalpha
double Beta21 = U2p * aiH.Coordo(1);
double Beta22 = U2p * aiH.Coordo(2);
// remplissage de la matrice mat
mat(ncot,1) = Beta21*Beta21;
mat(ncot,2) = 2.* Beta21*Beta22;
mat(ncot,3) = Beta22*Beta22;
// calcul des variations
// 1 - initialisation des vecteurs inter
ddNdv.Zero();
// 2 - init en fct des ddl
for (int idi=1; idi<= nbddl; idi++)
{ dCA(idi).Zero(); dCB(idi).Zero();dNi(idi).Zero();
dB(idi).Zero();dC(idi).Zero();dA(idi).Zero();
dU1(idi).Zero();d_g_h1.Zero();dG2(idi).Zero();
};
// 3 - variation non nulle des Ni pour certaine valeur d'indice
for (int ia=1;ia<= 3;ia++)
{ // on a 3 indices pour chaque ia, pour lesquels la variation est non nulle
// dans cette boucle on n'utilise pas systematiquement les 4 indices
int jnB = (indi(ncot+1)-1)*3+ia; // vérifié
int jnC = (indi(ncot+2)-1)*3+ia;
int jnA = 9+(ncot-1)*3+ia; //vérifié
// tout d'abord les variations des points
dB(jnB) = Ia(ia);
dC(jnC) = Ia(ia);dA(jnA) = Ia(ia);
dG2(jnA) = untiers * Ia(ia);dG2(jnB) = untiers * Ia(ia);dG2(jnC) = untiers * Ia(ia);
// vecteur CB
dCB(jnB) = dB(jnB); dCB(jnC) = - dC(jnC);
// Vecteur CA = A - C
dCA(jnC) = -dC(jnC); dCA(jnA) = dA(jnA);
// Vecteur Ni = (Util::ProdVec_coor(CB,CA)).Normer()
dNi(jnC) = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dCB(jnC),CA)+Util::ProdVec_coor(CB,dCA(jnC))),norNi);
dNi(jnB) = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dCB(jnB),CA)+Util::ProdVec_coor(CB,dCA(jnB))),norNi);
dNi(jnA) = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dCB(jnA),CA)+Util::ProdVec_coor(CB,dCA(jnA))),norNi);
// Vecteur U1 = CB / norCB;
dU1(jnC) = Util::VarUnVect_coor(CB,dCB(jnC),norCB);
dU1(jnB) = Util::VarUnVect_coor(CB,dCB(jnB),norCB);
};
for (int ib=1;ib<= 3;ib++) // indice de coordonnées
{ // on a 4 indices pour chaque ia, pour lesquels la variation est non nulle
// mais ici on utilise systematiquement les 4 indices d'ou une seconde boucle
// egalement moins de stockage intermediaire
// jnb[0] -> D, jnb[1] -> B, jnb[2] -> C, jnb[3] -> A,
int jnb[4] = {(indi(ncot)-1)*3+ib,(indi(ncot+1)-1)*3+ib,
(indi(ncot+2)-1)*3+ib,9+(ncot-1)*3+ib};
for (int jb=0;jb<=3;jb++)
{ // variation de g_h1
Coordonnee3 d_CG =dC(jnb[jb])-dG(jnb[jb]);
double d_CGsur2 = (d_CG * C_G);
double d_CG_U1 = d_CG * U1 + C_G * dU1(jnb[jb]);
double d_g_h1 = (d_CGsur2 - d_CG_U1)/g_h1;
// variation de h2G2
Coordonnee3 d_CG2 =dC(jnb[jb])-dG2(jnb[jb]);
Coordonnee3 d_H2G2 =d_CG2 - (d_CG2 * U1 + CG2 * dU1(jnb[jb])) * U1
- (CG2 * U1) *dU1(jnb[jb]);
double d_h2G2 = (H2G2 * d_H2G2)/h2G2;
// vecteur dNdV = ( Ni - N) / ( g_h1 + h2G2)
ddNdv = ((dNi(jnb[jb]) - dN(jnb[jb]))
- (Ni-N) * (d_g_h1 + d_h2G2) / ( g_h1 + h2G2)) / ( g_h1 + h2G2);
// variation de U2p
Coordonnee3 d_U2p = Util::ProdVec_coor(dN(jnb[jb]),U1)+Util::ProdVec_coor(N,dU1(jnb[jb]));
// courbure curbp(ncot) = dNdV * U2p;
dcurbp(jnb[jb])(ncot) = - (ddNdv * U2p + dNdV * d_U2p);
// -- variation des composantes de mat
double d_Beta21 = d_U2p * aiH.Coordo(1) + U2p * DaiH(jnb[jb]).Coordo(1);
double d_Beta22 = d_U2p * aiH.Coordo(2) + U2p * DaiH(jnb[jb]).Coordo(2);
dmat(jnb[jb])(ncot,1)= 2. * Beta21 * d_Beta21;
dmat(jnb[jb])(ncot,2)= 2. * (d_Beta21 * Beta22 + Beta21 * d_Beta22);
dmat(jnb[jb])(ncot,3)= 2. * Beta22 * d_Beta22;
};
}
}
else // cas ou il n'y a pas de noeud externe
{ // la courbure est nulle
curbp(ncot) = 0.;
// la direction de la courbure nulle est normale au cote
// c-a-d le produit vectoriel de la normale avec le cote
CB = B - C;
double norCB = CB.Norme();
U1 = CB / norCB;
H2G2 = Util::ProdVec_coor(N,U1);
// calcul de la matrice de passage de Ualpha a aalpha
double Beta21 = H2G2 * aiH.Coordo(1);
double Beta22 = H2G2 * aiH.Coordo(2);
// remplissage de la matrice mat
mat(ncot,1) = Beta21*Beta21;
mat(ncot,2) = 2.* Beta21*Beta22;
mat(ncot,3) = Beta22*Beta22;
}
}; //-- fin de la boucle sur les trois cotés
// calcul du tenseur de courbure dans la base ai
Mat_pleine mati = mat.Inverse();
curbi = mati * curbp;
curb.Coor(1,1)=curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(2); curb.Coor(2,2)=curbi(3);
// idem pour les variations
Mat_pleine dmata(mat); // une matrice intermédiaire
for (int ddl=1;ddl<=nbddl;ddl++)
{ // tout d'abord le calcul de la variation de l'inverse de mat
dmata = -(mati * dmat(ddl) * mati);
// puis la variation de la courbure
curbi = mati * dcurbp(ddl) + dmata * curbp;
Tenseur_ns2BB& dcur = dcurb(ddl);
dcur.Coor(1,1)=curbi(1); dcur.Coor(1,2)=dcur.Coor(2,1)=curbi(2); dcur.Coor(2,2)=curbi(3);
};
};
// test si la courbure est anormalement trop grande
// inf_normale : indique en entrée le det mini pour la courbure en locale
// retour 1: si tout est ok,
// 0: une courbure trop grande a été détecté
int Met_Sfe1::Test_courbure_anormale1(double inf_normale,const Tableau<Coordonnee3>& tab_coor
, const BaseB & aiB , const BaseH & aiH)
{// pour l'instant cette méthode n'est pas alimentée on le signale
if (ParaGlob::NiveauImpression()> 8)
cout << "\n *** attention, cette methode n'est pas implante pour la courbure 1 "
<< "\n Met_Sfe1::Test_courbure_anormale1(... "<<endl;
return 1; // retour par défaut
};
// routine generale de calcul de la courbure
// cas 2: moyenne des normales de part et autres des arrêtes
Tenseur_ns2BB Met_Sfe1::Courbure2
( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH &
,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
{ Tenseur_ns2BB curb; // : les tenseurs de courbures d'abord dans les axes
// au lieu d'utiliser des expressions condensées (en commentaire)
// on utilise exactement la même procédure que pour dcurb
// car sinon on a de très petites différences qui ensuite
// conduisent à des différences non négligeables sur les def
// calcul de la normale centrale
// Coordonnee3 N((Util::ProdVec_coorBN(aiB(1),aiB(2))).Normer());
// calcul de la normale centrale
Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
double nor = NN.Norme(); // la norme
Coordonnee3 N = NN / nor;
// le meme calcul est effectue pour les trois cotes du triangle principale
//1- on utilise indi un tableau d'indice pour ne pas etre embete par les modulos 3
// 2- on boucle sur le cotes du triangle
Tableau <Coordonnee3> t_Ni(3); // les 3 normales du triangle exterieur
//cout << "\n N= "<<N;
Coordonnee3 NNi,Ni,NietN; // normale du triangle exterieur
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// on cherche a calculer la longueur H2G2
// simplification de l'ecriture (vérifiée)
const Coordonnee3& D = tab_coor(indi(ncot));
const Coordonnee3& B = tab_coor(indi(ncot+1));
const Coordonnee3& C = tab_coor(indi(ncot+2));
const Coordonnee3& A = tab_coor(indi(ncot)+3);
if (D != A) // cas ou le noeud externe existe
{// calcul de la normale du triangle exterieur
Coordonnee3 CB = B - C;
Coordonnee3 CA = A - C;
//t_Ni(ncot) = ((Util::ProdVec_coor(CB,CA)).Normer() + N).Normer();
NNi = Util::ProdVec_coor(CB,CA);
double norNi = NNi.Norme();
Ni = NNi/norNi;
NietN=(Ni + N);
double norNietN = NietN.Norme();
t_Ni(ncot) = NietN / norNietN;
}
else // cas ou il n'y a pas de noeud externe
{ // on utilise la normale actuelle comme normale sur l'arrête
t_Ni(ncot) = N;
};
};
// calcul de la variation de la normale
Coordonnee3 dN_1 = 2.*(t_Ni(1)-t_Ni(2));
Coordonnee3 dN_2 = 2.*(t_Ni(1)-t_Ni(3));
// courbure
curb.Coor(1,1)=-dN_1 * aiB.Coordo(1);
curb.Coor(1,2)=-dN_1 * aiB.Coordo(2);
curb.Coor(2,1)=-dN_2 * aiB.Coordo(1);
curb.Coor(2,2)=-dN_2 * aiB.Coordo(2);
// retour de la courbure
return curb;
};
// routine generale de calcul de la courbure et de sa variation
// en sortie : curb , la courbure b11,b12,b22
// dcurb , la variation de courbure.
// On ne considere dans les variations que les termes relatif aux ddl qui font varier
// les normales
// cas 2: moyenne des normales de part et autres des arrêtes
void Met_Sfe1::Dcourbure2
( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB
, const Tableau <BaseB>& DaiB,const BaseH &
, const Tableau <BaseH>& ,Tenseur_ns2BB& curb,TabOper <Tenseur_ns2BB>& dcurb
,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
{ int nbddl = (nbNoeudCentral+3)*3; //18; // nombre de ddl total
Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations
double untiers=1./3.;
// calcul de la normale centrale
Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
double nor = NN.Norme(); // la norme
Coordonnee3 N = NN / nor;
Coordonnee3 DNN; // vecteur intermediaire
// calcul de la variation de la normale centrale et du centre de gravité: vérifiée
Tableau <Coordonnee3> dN(nbddl,Nul);
// on s'intéresse ici qu'aux ddl de la face centrale
for (int inc=1;inc<=3;inc++) // indice du numéro de noeud de la face centrale
for (int ib=1;ib<=3;ib++)
{ // cas de la normale
int i1=(inc-1)*3+ib;
DNN = Util::ProdVec_coorBN(DaiB(i1)(1),aiB(2))
+ Util::ProdVec_coorBN(aiB(1),DaiB(i1)(2));
dN(i1) = Util::VarUnVect_coor(NN,DNN,nor);
};
// le meme calcul est effectue pour les trois cotes du triangle principale
//1- on utilise un tableau d'indice pour ne pas etre embete par les modulos 3
// 2- on boucle sur le cotes du triangle
Coordonnee3 NNi,Ni,NietN; // normale du triangle exterieur
// -------- declaration pour la boucle sur les 3 cotes ------------
Coordonnee3 CB,CA;
Tableau <Coordonnee3> t_Ni(3); // les 3 normales du triangle exterieur
Tableau <Tableau <Coordonnee3> > d_t_Ni(nbddl,t_Ni); // le tableau des variations
Coordonnee3 dA,dB,dC,dNi;
//cout << "\n NN= "<<NN <<" nor= "<< nor <<" N= "<<N;
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// on cherche a calculer la longueur H2G2
// simplification de l'ecriture
const Coordonnee3& D = tab_coor(indi(ncot));
const Coordonnee3& B = tab_coor(indi(ncot+1));
const Coordonnee3& C = tab_coor(indi(ncot+2));
const Coordonnee3& A = tab_coor(indi(ncot)+3);
if (D != A) // cas ou le noeud externe existe
{CB = B - C;
double norCB = CB.Norme();
// calcul de la normale du triangle exterieur
CA = A - C;
NNi = Util::ProdVec_coor(CB,CA);
double norNi = NNi.Norme();
Ni = NNi/norNi;
NietN=(Ni + N);
double norNietN = NietN.Norme();
t_Ni(ncot) = NietN / norNietN;
//cout << "\n CB= "<<CB <<" CA= "<< CA
// << "t_Ni("<<ncot<<")= "<<t_Ni(ncot);
//cout << "\n norNi= "<< norNi<< " norNietN= "<< norNietN;
// calcul des variations
for (int ia=1;ia<= 3;ia++)
{{// indice de A: jnA
// int jnA = 9+(ncot-1)*3+ia;
//modif pour généraliser int jnA = 9+(ncot-1)*3+ia; //vérifié
int jnA = nbNoeudCentral*3+(ncot-1)*3+ia;
dA = Ia(ia);
//dCA = dA;
// Vecteur Ni = (Util::ProdVec_coor(CB,CA)).Normer()
dNi = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(CB,dA)),norNi);
d_t_Ni(jnA)(ncot) = Util::VarUnVect_coor(NietN,(dNi+dN(jnA)),norNietN);
}
{// indice de B: jnB
int jnB = (indi(ncot+1)-1)*3+ia;
dB = Ia(ia);
//dCB = dB;
// Vecteur Ni = (Util::ProdVec_coor(CB,CA)).Normer()
dNi = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dB,CA)),norNi);
d_t_Ni(jnB)(ncot) = Util::VarUnVect_coor(NietN,(dNi+dN(jnB)),norNietN);
}
{// indice de C: jnC
int jnC = (indi(ncot+2)-1)*3+ia;
dC = Ia(ia);
//dCB = -dC; dCA= -dC;
// Vecteur Ni = (Util::ProdVec_coor(CB,CA)).Normer()
dNi = Util::VarUnVect_coor(NNi,-(Util::ProdVec_coor(dC,CA)+Util::ProdVec_coor(CB,dC)),norNi);
d_t_Ni(jnC)(ncot) = Util::VarUnVect_coor(NietN,(dNi+dN(jnC)),norNietN);
}
{// indice de D: jnD
int jnD = (indi(ncot)-1)*3+ia;
// Vecteur Ni = (Util::ProdVec_coor(CB,CA)).Normer()
// dNi = vecteur nul;
d_t_Ni(jnD)(ncot) = Util::VarUnVect_coor(NietN,(dN(jnD)),norNietN);
}
}; // fin du for sur ia
}
else // cas ou il n'y a pas de noeud externe
{ // on utilise la normale actuelle comme normale sur l'arrête
t_Ni(ncot) = N;
//cout << "\n t_Ni("<<ncot<<")= "<<t_Ni(ncot);
// et pour les variations, celle de N
for (int ia=1;ia<= 3;ia++)
{ // indice de B: jnB
int jnB = (indi(ncot+1)-1)*3+ia;
d_t_Ni(jnB)(ncot) = dN(jnB);
// indice de C: jnC
int jnC = (indi(ncot+2)-1)*3+ia;
d_t_Ni(jnC)(ncot) = dN(jnC);
// indice de D: jnD
int jnD = (indi(ncot)-1)*3+ia;
d_t_Ni(jnD)(ncot) = dN(jnD);
}; // fin du for sur ia
};
}; //-- fin de la boucle sur les trois cotés
// calcul de la variation de la normale
Coordonnee3 dN_1 = 2.*(t_Ni(1)-t_Ni(2));
Coordonnee3 dN_2 = 2.*(t_Ni(1)-t_Ni(3));
// courbure
curb.Coor(1,1)=-dN_1 * aiB.Coordo(1);
curb.Coor(1,2)=-dN_1 * aiB.Coordo(2);
curb.Coor(2,1)=-dN_2 * aiB.Coordo(1);
curb.Coor(2,2)=-dN_2 * aiB.Coordo(2);
//cout << "\n 2) curb= "; curb.Ecriture(cout);
//cout << endl;
// calcul des variations de la courbure par rapport au ddl
for (int ddl=1;ddl<=nbddl;ddl++)
{Tenseur_ns2BB& dcur = dcurb(ddl);
dcur.Coor(1,1) = (-2.) * ((d_t_Ni(ddl)(1)-d_t_Ni(ddl)(2)) * aiB.Coordo(1)) - dN_1 * DaiB(ddl).Coordo(1);
dcur.Coor(1,2)= -(2.* (d_t_Ni(ddl)(1)-d_t_Ni(ddl)(2)) * aiB.Coordo(2) + (dN_1 * DaiB(ddl).Coordo(2)));
dcur.Coor(2,1)= -(2.*(d_t_Ni(ddl)(1)-d_t_Ni(ddl)(3)) * aiB.Coordo(1) + (dN_2 * DaiB(ddl).Coordo(1)));
dcur.Coor(2,2) = (-2.) * ((d_t_Ni(ddl)(1)-d_t_Ni(ddl)(3)) * aiB.Coordo(2)) - dN_2 * DaiB(ddl).Coordo(2);
};
};
// test si la courbure est anormalement trop grande
// inf_normale : indique en entrée le det mini pour la courbure en locale
// retour 1: si tout est ok,
// 0: une courbure trop grande a été détecté
int Met_Sfe1::Test_courbure_anormale2(double inf_normale,const Tableau<Coordonnee3>& tab_coor
, const BaseB & aiB , const BaseH & aiH)
{// pour l'instant cette méthode n'est pas alimentée on le signale
if (ParaGlob::NiveauImpression()> 8)
cout << "\n *** attention, cette methode n'est pas implante pour la courbure 2 "
<< "\n Met_Sfe1::Test_courbure_anormale2(... "<<endl;
return 1; // retour par défaut
};
// routine (3) generale de calcul de la courbure
// cas 3: moyenne pondérée des normales de part et autres des arrêtes, et localisation du point sur l'arrête
Tenseur_ns2BB Met_Sfe1::Courbure3
( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH & aiH
,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
{ // calcul de la normale centrale
Coordonnee3 N((Util::ProdVec_coorBN(aiB(1),aiB(2))).Normer());
// calcul du centre de gravite
Coordonnee3 G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.;
// les coordonnées locales non nulles des points origines des normales sur les arrêtes
// on les initialise comme si les points étaient au milieu des arrêtes
double theta1_H1=0.5; double theta2_H1=0.5;
double theta2_H2=0.5; double theta1_H3=0.5;
// origine
const Coordonnee3& Op=tab_coor(1);// origine du repère local
// le meme calcul est effectue pour les trois cotes du triangle principale
//1- on utilise indi un tableau d'indice pour ne pas etre embete par les modulos 3
// 2- on boucle sur le cotes du triangle
Tableau <Coordonnee3> t_Ni(3); // normale du triangle exterieur
// def des différentes grandeurs
Coordonnee3 Ne; // normale courante du triangle extérieur
Coordonnee3 CB,CA,CD; Coordonnee3 G2;
Coordonnee3 U1,HaA,OpHi;
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// on cherche a calculer la longueur H2G2
// simplification de l'ecriture (vérifiée)
const Coordonnee3& D = tab_coor(indi(ncot));
const Coordonnee3& B = tab_coor(indi(ncot+1));
const Coordonnee3& C = tab_coor(indi(ncot+2));
const Coordonnee3& A = tab_coor(indi(ncot)+3);
if (D != A) // cas ou le noeud externe existe
{// calcul de G2, centre de gravite du triangle exterieur
G2 = (B + C + A) / 3.;
// calcul de U1
CB = B - C;
CA = A - C;
U1 = CB / CB.Norme();
// longueur HaA et HdD et calcul du vecteur HaA norme
// tout d'abord HaA ce situe dans le plan de la facette externe
HaA = CA - (CA * U1) * U1;
double haA = HaA.Norme();
double unSurhaA= 1./haA;
CD = D - C;
double dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ;
// calcul de la normale du triangle exterieur
Ne = (Util::ProdVec_coor(CB,CA)).Normer();
// calcul de la normale
t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
// calcul des coordonnées de la position des normales sur les arrêtes
double beta = 0.5 * (G + G2 - 2. * B) * U1;
OpHi=B+beta*U1-Op;
switch (ncot)
{case 1: theta1_H1=OpHi*aiH.Coordo(1);theta2_H1=OpHi*aiH.Coordo(2); break;
case 2: theta2_H2=OpHi*aiH.Coordo(2); break;
case 3: theta1_H3=OpHi*aiH.Coordo(1); break;
};
}
else // cas ou il n'y a pas de noeud externe
{// on utilise comme normale extérieure la normale du triangle central
t_Ni(ncot) = N;
// on utilise le point milieu de l'arrête pour le positionnement de la normale
// normalement les points sont déjà initialisés
};
};
// calcul des coefficients des fonctions d'interpolation
double s= theta1_H1 * theta2_H2 - theta1_H3 * theta2_H2 + theta2_H1 * theta1_H3;
double a_r1= theta2_H2 / s; double b_r1 = theta1_H3 / s; double c_r1 = -theta1_H3 * theta2_H2/s;
double a_r2= -theta2_H1 / s; double b_r2 = 1./theta2_H2 - theta1_H3 * theta2_H1 /theta2_H2/ s;
double c_r2 = -theta1_H3 * theta2_H1/s;
double a_r3= 1./theta1_H3 - theta1_H1 * theta2_H2 /theta1_H3/ s; double b_r3 = - theta1_H1 / s;
double c_r3 = theta2_H2 * theta1_H1 /s;
// calcul de la variation de la normale
Coordonnee3 d_N_1= t_Ni(1) * a_r1 + t_Ni(2) * a_r2 + t_Ni(3) * a_r3 ;
Coordonnee3 d_N_2= t_Ni(1) * b_r1 + t_Ni(2) * b_r2 + t_Ni(3) * b_r3 ;
// calcul de la courbure
Tenseur_ns2BB curb;
curb.Coor(1,1) = - d_N_1 * aiB.Coordo(1); curb.Coor(1,2) = - d_N_1 * aiB.Coordo(2);
curb.Coor(2,1) = - d_N_2 * aiB.Coordo(1); curb.Coor(2,2) = - d_N_2 * aiB.Coordo(2);
// retour de la courbure
return curb;
};
// test si la courbure est anormalement trop grande
// inf_normale : indique en entrée le det mini pour la courbure en locale
// retour 1: si tout est ok,
// 0: une courbure trop grande a été détecté
int Met_Sfe1::Test_courbure_anormale3(double inf_normale,const Tableau<Coordonnee3>& tab_coor
, const BaseB & aiB , const BaseH & aiH)
{// pour l'instant cette méthode n'est pas alimentée on le signale
if (ParaGlob::NiveauImpression()> 8)
cout << "\n *** attention, cette methode n'est pas implante pour la courbure 3 "
<< "\n Met_Sfe1::Test_courbure_anormale1(... "<<endl;
return 1; // retour par défaut
};
// routine (3) generale de calcul de la courbure et de sa variation
// en sortie : curb , la courbure b11,b12,b22
// dcurb , la variation de courbure.
// On ne considere dans les variations que les termes relatif aux ddl qui font varier
// les normales
// cas 3: moyenne pondérée des normales de part et autres des arrêtes, et localisation du point sur l'arrête
void Met_Sfe1::Dcourbure3
( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB
, const Tableau <BaseB>& DaiB,const BaseH & aiH
, const Tableau <BaseH>& DaiH,Tenseur_ns2BB& curb,TabOper <Tenseur_ns2BB>& dcurb
,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
{ int nbddl = (nbNoeudCentral+3)*3; //18; // nombre de ddl total
Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations
// les coordonnées locales non nulles des points origines des normales sur les arrêtes
// on les initialise comme si les points étaient au milieu des arrêtes
double theta1_H1=0.5; double theta2_H1=0.5;
double theta2_H2=0.5; double theta1_H3=0.5;
// origine
const Coordonnee3& Op=tab_coor(1);// origine du repère local
Tableau <Coordonnee3> dOp(nbddl,Nul);
dOp(1)=Ia(1);dOp(2)=Ia(2);dOp(3)=Ia(3);
// et les variations
Vecteur d_theta1_H1(nbddl), d_theta2_H1(nbddl), d_theta2_H2(nbddl), d_theta1_H3(nbddl);
double untiers=1./3.;
// calcul du centre de gravite
Coordonnee3 G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) * untiers;
// et le tableau de ses variations
Tableau <Coordonnee3> dG(nbddl,Nul);
// calcul de la normale centrale
Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
double nor = NN.Norme(); // la norme
Coordonnee3 N = NN / nor;
Coordonnee3 DNN; // vecteur intermediaire
// calcul de la variation de la normale centrale et du centre de gravité: vérifiée
Tableau <Coordonnee3> dN(nbddl,Nul);
// on s'intéresse ici qu'aux ddl de la face centrale
for (int inc=1;inc<=3;inc++) // indice du numéro de noeud de la face centrale
for (int ib=1;ib<=3;ib++)
{ // cas de la normale
int i1=(inc-1)*3+ib;
DNN = Util::ProdVec_coorBN(DaiB(i1)(1),aiB(2))
+ Util::ProdVec_coorBN(aiB(1),DaiB(i1)(2));
dN(i1) = Util::VarUnVect_coor(NN,DNN,nor);
// cas du centre de gravité
dG(i1)=untiers*Ia(ib);
};
// le meme calcul est effectue pour les trois cotes du triangle principale
//1- on utilise un tableau d'indice pour ne pas etre embete par les modulos 3
// 2- on boucle sur le cotes du triangle
Tableau <Coordonnee3> t_Ni(3); // normale du triangle exterieur
Tableau < Tableau <Coordonnee3> > d_t_Ni(nbddl,t_Ni); // les dérivées, initialisées à celle de dN
Coordonnee3 NNi; // normale du triangle exterieur
// def des différentes grandeurs
Coordonnee3 Ne,NNe; // normale courante du triangle extérieur et sa valeur non normée
Coordonnee3 CB,CA,CD; Coordonnee3 G2,dG2;
Coordonnee3 U1,HaA,OpHi; Coordonnee3 BgimoinsGB;
// -------- declaration pour la boucle sur les 3 cotes ------------
Coordonnee3 dA,dB,dC,dD,d_U1,dNe,d_HaA,d_NNi,d_Op,d_OpHi;
Coordonnee3 t_NNi; // normale courante non normée sur l'arrête
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// on cherche a calculer la longueur H2G2
// simplification de l'ecriture
const Coordonnee3& D = tab_coor(indi(ncot));
const Coordonnee3& B = tab_coor(indi(ncot+1));
const Coordonnee3& C = tab_coor(indi(ncot+2));
const Coordonnee3& A = tab_coor(indi(ncot)+3);
if (D != A) // cas ou le noeud externe existe
{// calcul de G2, centre de gravite du triangle exterieur
G2 = (B + C + A) * untiers;
// calcul de U1
CB = B - C;
CA = A - C;
double norCB = CB.Norme();
U1 = CB / norCB;
// longueur HaA et HdD et calcul du vecteur HaA norme
// tout d'abord HaA ce situe dans le plan de la facette externe
HaA = CA - (CA * U1) * U1;
double haA = HaA.Norme();
double unSurhaA= 1./haA;
CD = D - C;
double dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ;
// double distH = dHd + haA;
// double inv_distH = 1./distH; double inv_distH2 = inv_distH * inv_distH;
// calcul de la normale du triangle exterieur
NNe = (Util::ProdVec_coor(CB,CA));
double norNe = NNe.Norme();
Ne = NNe/norNe;
// calcul de la normale
t_NNi = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA)));
double nor_t_NNi = t_NNi.Norme();
t_Ni(ncot) = t_NNi / nor_t_NNi;
// calcul des coordonnées de la position des normales sur les arrêtes
BgimoinsGB = (G + G2 - 2. * B);
double beta = 0.5 * BgimoinsGB * U1;
OpHi=B+beta*U1-Op;
switch (ncot)
{case 1: theta1_H1=OpHi*aiH.Coordo(1);theta2_H1=OpHi*aiH.Coordo(2); break;
case 2: theta2_H2=OpHi*aiH.Coordo(2); break;
case 3: theta1_H3=OpHi*aiH.Coordo(1); break;
};
// calcul des variations
for (int ia=1;ia<= 3;ia++)
{{// indice de A: jnA, noeud extérieur
// modif pour généraliser int jnA = 9+(ncot-1)*3+ia;
int jnA = nbNoeudCentral*3+(ncot-1)*3+ia;
dA = Ia(ia);
//dCA = dA;
// Vecteur Ne = (Util::ProdVec_coor(CB,CA)).Normer()
dNe = Util::VarUnVect_coor(NNe,(Util::ProdVec_coor(CB,dA)),norNe);
// ici d_CB=0 donc d_U1=0
// ici d_CD=0 donc d_dHd=0
d_HaA = dA - (dA * U1) * U1;
double d_haA = HaA * d_HaA * unSurhaA ;
// t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
d_NNi = (dNe*(dHd/(dHd+haA)) + dN(jnA)*(haA/(dHd+haA)))
+ (Ne*(-dHd*d_haA/(dHd+haA)/(dHd+haA)) + N*(d_haA/(dHd+haA)))
+ ( N*(-haA*d_haA/(dHd+haA)/(dHd+haA)));
d_t_Ni(jnA)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
// maintenant la variation des theta
// G2 = (B + C + A) * untiers;
dG2 = dA * untiers;
// beta = 0.5 * (G + G2 - 2. * B) * U1;
// ici dG(jnA) = 0
double d_beta = 0.5 * (dG2 * U1);
d_OpHi=d_beta*U1; // OpHi=B+beta*U1-Op;
switch (ncot)
{case 1: {d_theta1_H1(jnA)=d_OpHi*aiH.Coordo(1);d_theta2_H1(jnA)=d_OpHi*aiH.Coordo(2); break;}
case 2: d_theta2_H2(jnA)=d_OpHi*aiH.Coordo(2); break;
case 3: d_theta1_H3(jnA)=d_OpHi*aiH.Coordo(1); break;
};
}
{// indice de B: jnB, noeud de la facette centrale
int jnB = (indi(ncot+1)-1)*3+ia;
dB = Ia(ia); // d_CB = dB, d_CA=0
d_U1= Util::VarUnVect_coor(CB,dB,norCB);
// HaA = CA - (CA * U1) * U1;
d_HaA = - ((CA*d_U1) * U1 +(CA * U1) * d_U1);
double d_haA = HaA * d_HaA * unSurhaA ;
// d_CD = 0; dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ;
double d_dHd = - ( CD * d_U1)* ( CD * U1) / dHd ;
// Vecteur Ne = (Util::ProdVec_coor(CB,CA)).Normer()
dNe = Util::VarUnVect_coor(NNe,(Util::ProdVec_coor(dB,CA)),norNe);
// t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
d_NNi = (dNe*(dHd/(dHd+haA)) + dN(jnB)*(haA/(dHd+haA)))
+ (Ne*(d_dHd/(dHd+haA)) + N*(d_haA/(dHd+haA)))
+ (Ne*(-dHd*d_haA/(dHd+haA)/(dHd+haA))+ N*(-haA*d_haA/(dHd+haA)/(dHd+haA)))
+ (Ne*(-dHd*d_dHd/(dHd+haA)/(dHd+haA)) + N*(-haA*d_dHd/(dHd+haA)/(dHd+haA)));
d_t_Ni(jnB)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
// maintenant la variation des theta
// G2 = (B + C + A) * untiers;
dG2 = dB * untiers;
// beta = 0.5 * (G + G2 - 2. * B) * U1;
double d_beta = 0.5 * ((dG(jnB) + dG2 - 2. * dB) * U1 + BgimoinsGB * d_U1);
// if (ncot == 3) {d_Op=dB;} else { d_Op= Nul;};
// d_OpHi=dB + d_beta*U1+beta*d_U1-d_Op;
d_OpHi=dB + d_beta*U1+beta*d_U1-dOp(jnB);
switch (ncot)
{case 1: {d_theta1_H1(jnB)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnB).Coordo(1);
d_theta2_H1(jnB)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnB).Coordo(2); break; }
case 2: d_theta2_H2(jnB)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnB).Coordo(2); break;
case 3: d_theta1_H3(jnB)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnB).Coordo(1); break;
};
}
{// indice de C: jnC
int jnC = (indi(ncot+2)-1)*3+ia;
dC = Ia(ia);
//dCB = -dC; dCA= -dC; U1 = CB / norCB;
d_U1= Util::VarUnVect_coor(CB,-dC,norCB);
// HaA = CA - (CA * U1) * U1;
d_HaA = -dC - ((-dC * U1)+(CA*d_U1)) * U1 -(CA * U1) * d_U1;
double d_haA = HaA * d_HaA * unSurhaA ;
// d_CD = -dC; dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ;
double d_dHd = - (CD*dC + ( CD * d_U1- dC*U1) * ( CD * U1) ) / dHd ;
// Vecteur Ne = (Util::ProdVec_coor(CB,CA)).Normer()
dNe = Util::VarUnVect_coor(NNe,(-(Util::ProdVec_coor(dC,CA)+Util::ProdVec_coor(CB,dC))),norNe);
// t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
d_NNi = (dNe*(dHd/(dHd+haA)) + dN(jnC)*(haA/(dHd+haA)))
+ (Ne*(d_dHd/(dHd+haA)) + N*(d_haA/(dHd+haA)))
+ (Ne*(-dHd*d_haA/(dHd+haA)/(dHd+haA)) + N*(-haA*d_haA/(dHd+haA)/(dHd+haA)))
+ (Ne*(-dHd*d_dHd/(dHd+haA)/(dHd+haA)) + N*(-haA*d_dHd/(dHd+haA)/(dHd+haA)));
d_t_Ni(jnC)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
// maintenant la variation des theta
// G2 = (B + C + A) * untiers;
dG2 = dC * untiers;
// beta = 0.5 * (G + G2 - 2. * B) * U1;
double d_beta = 0.5 * ((dG(jnC) + dG2 ) * U1 + BgimoinsGB * d_U1);
// if (ncot == 2) {d_Op=dC;} else { d_Op= Nul;};
// d_OpHi=d_beta*U1+beta*d_U1-d_Op;
d_OpHi=d_beta*U1+beta*d_U1-dOp(jnC);
switch (ncot)
{case 1: {d_theta1_H1(jnC)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnC).Coordo(1);
d_theta2_H1(jnC)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnC).Coordo(2); break; }
case 2: d_theta2_H2(jnC)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnC).Coordo(2); break;
case 3: d_theta1_H3(jnC)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnC).Coordo(1); break;
};
}
{// indice de D: jnD
int jnD = (indi(ncot)-1)*3+ia;
// d_CB = 0; d_CA = 0 ==> d_U1=0
// d_HaA = 0, et d_haA = 0
// d_CD = dD;
// dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ;
double d_dHd = (CD*dD - (dD * U1) * ( CD * U1) ) / dHd ;
// Vecteur Ne = (Util::ProdVec_coor(CB,CA)).Normer()
//dNe = 0
// t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
d_NNi = (dN(jnD)*(haA/(dHd+haA)))
+ (Ne*(d_dHd/(dHd+haA)))
+ (Ne*(-dHd*d_dHd/(dHd+haA)/(dHd+haA)) + N*(-haA*d_dHd/(dHd+haA)/(dHd+haA)));
d_t_Ni(jnD)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
// maintenant la variation des theta
// G2 = (B + C + A) * untiers;
//dG2 = 0;
// beta = 0.5 * (G + G2 - 2. * B) * U1;
double d_beta = 0.5 * (dG(jnD) * U1);
// if (ncot == 1) {d_Op=dD;} else { d_Op= Nul;};
// dB=0
// d_OpHi=d_beta*U1-d_Op;
d_OpHi=d_beta*U1-dOp(jnD);
switch (ncot)
{case 1: {d_theta1_H1(jnD)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnD).Coordo(1);
d_theta2_H1(jnD)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnD).Coordo(2); break; }
case 2: d_theta2_H2(jnD)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnD).Coordo(2); break;
case 3: d_theta1_H3(jnD)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnD).Coordo(1); break;
};
}
}; // fin du for sur ia
}
else // cas ou il n'y a pas de noeud externe
{ // on utilise la normale actuelle comme normale sur l'arrête
t_Ni(ncot) = N;
// et pour les variations, celle de N
for (int ia=1;ia<= 3;ia++)
{ // indice de B: jnB
int jnB = (indi(ncot+1)-1)*3+ia;
d_t_Ni(jnB)(ncot) = dN(jnB);
// indice de C: jnC
int jnC = (indi(ncot+2)-1)*3+ia;
d_t_Ni(jnC)(ncot) = dN(jnC);
// indice de D: jnD
int jnD = (indi(ncot)-1)*3+ia;
d_t_Ni(jnD)(ncot) = dN(jnD);
}; // fin du for sur ia
// concernant les variations des coordonnées des positions des normales,
// ces variations sont nulles, on n'y touche pas
};
}; //-- fin de la boucle sur les trois cotés
// calcul des coefficients des fonctions d'interpolation
double s= theta1_H1 * theta2_H2 - theta1_H3 * theta2_H2 + theta2_H1 * theta1_H3;
double a_r1= theta2_H2 / s; double b_r1 = theta1_H3 / s; double c_r1 = -theta1_H3 * theta2_H2/s;
double a_r2= -theta2_H1 / s; double b_r2 = 1./theta2_H2 - theta1_H3 * theta2_H1 /theta2_H2/ s;
double c_r2 = -theta1_H3 * theta2_H1/s;
double a_r3= 1./theta1_H3 - theta1_H1 * theta2_H2 /theta1_H3/ s; double b_r3 = - theta1_H1 / s;
double c_r3 = theta2_H2 * theta1_H1 /s;
// calcul de la variation de la normale
Coordonnee3 d_N_1= t_Ni(1) * a_r1 + t_Ni(2) * a_r2 + t_Ni(3) * a_r3 ;
Coordonnee3 d_N_2= t_Ni(1) * b_r1 + t_Ni(2) * b_r2 + t_Ni(3) * b_r3 ;
// calcul de la courbure (non symétrique)
curb.Coor(1,1) = - d_N_1 * aiB.Coordo(1); curb.Coor(1,2) = - d_N_1 * aiB.Coordo(2);
curb.Coor(2,1) = - d_N_2 * aiB.Coordo(1); curb.Coor(2,2) = - d_N_2 * aiB.Coordo(2);
// retour de la courbure
double unsurs = 1./s; double unsurs2=unsurs * unsurs ;
for (int ddl=1;ddl<=nbddl;ddl++)
{ // tout d'abord le calcul de la variation des coeff des fonctions d'interpolation
double d_s = d_theta1_H1(ddl) * theta2_H2 + theta1_H1 * d_theta2_H2(ddl)
- d_theta1_H3(ddl) * theta2_H2 - theta1_H3 * d_theta2_H2(ddl)
+ d_theta2_H1(ddl) * theta1_H3 + theta2_H1 * d_theta1_H3(ddl);
double d_a_r1 = d_theta2_H2(ddl) * unsurs - theta2_H2 * d_s * unsurs2;
double d_b_r1 = d_theta1_H3(ddl) * unsurs - theta1_H3 * d_s * unsurs2;
double d_a_r2 = -d_theta2_H1(ddl) * unsurs + theta2_H1 * d_s * unsurs2;
double d_b_r2 = - d_theta2_H2(ddl)/theta2_H2/theta2_H2
- d_theta1_H3(ddl) * theta2_H1 /theta2_H2* unsurs
- theta1_H3 * d_theta2_H1(ddl) /theta2_H2* unsurs
+ theta1_H3 * theta2_H1 * d_theta2_H2(ddl)/theta2_H2/theta2_H2 * unsurs
+ theta1_H3 * theta2_H1 * d_s /theta2_H2 * unsurs2;
double d_a_r3 = -d_theta1_H3(ddl)/theta1_H3/theta1_H3
- d_theta1_H1(ddl) * theta2_H2 /theta1_H3* unsurs
- theta1_H1 * d_theta2_H2(ddl) /theta1_H3* unsurs
+ theta1_H1 * theta2_H2 * d_theta1_H3(ddl)/theta1_H3/theta1_H3 * unsurs
+ theta1_H1 * theta2_H2 * d_s/theta1_H3 * unsurs2;
double d_b_r3 = - d_theta1_H1(ddl) * unsurs + theta1_H1 * d_s * unsurs2;
// calcul de la variation de la normale
Coordonnee3 d_d_N_1= d_t_Ni(ddl)(1) * a_r1 + t_Ni(1) * d_a_r1
+ d_t_Ni(ddl)(2) * a_r2 + t_Ni(2) * d_a_r2
+ d_t_Ni(ddl)(3) * a_r3 + t_Ni(3) * d_a_r3;
Coordonnee3 d_d_N_2= d_t_Ni(ddl)(1) * b_r1 + t_Ni(1) * d_b_r1
+ d_t_Ni(ddl)(2) * b_r2 + t_Ni(2) * d_b_r2
+ d_t_Ni(ddl)(3) * b_r3 + t_Ni(3) * d_b_r3;
// variation de la courbure
Tenseur_ns2BB& dcur = dcurb(ddl);
dcur.Coor(1,1) = - (d_d_N_1 * aiB.Coordo(1)+ d_N_1 * DaiB(ddl).Coordo(1));
dcur.Coor(1,2) = - (d_d_N_1 * aiB.Coordo(2)- d_N_1 * DaiB(ddl).Coordo(2));
dcur.Coor(2,1) = - (d_d_N_2 * aiB.Coordo(1)- d_N_2 * DaiB(ddl).Coordo(2));
dcur.Coor(2,2) = - (d_d_N_2 * aiB.Coordo(2)- d_N_2 * DaiB(ddl).Coordo(1));
};
};
Reference in a new issue View git blame Copy permalink