// 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-2022 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 "Util.h"
#include "MathUtil.h"
#include "MathUtil2.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
// =========================================================================================
// routine (4) generale de calcul de la courbure
// cas 4: calcul de la courbure à partir d'un polynôme passant par les 6 points
Tenseur_ns2BB Met_Sfe1::Courbure4
( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH & aiH
,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ // calcul de la normale centrale
// on utilise des expressions strictement identique ici et dans dcourbure
// de manière à éviter les erreurs d'arrondis
//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;
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
Vecteur sm(3); // le second membre pour le calcul de la courbure
// 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
Coordonnee3 U1,CB,V,OpA;
Coordonnee3 Mm,OpMm;
Coordonnee3 OpAp; // le vecteur horizontal correspondant à OpA
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// 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);
// examen de condition limite éventuelle
// ncot % 3 +1: = le coté en face du noeud ncot
if (tabTypeCL(ncot % 3 +1) != TANGENTE_CL)
{ // cas sans conditions limites particulières
if (D != A) // cas ou le noeud externe existe
{OpA = A - Op;
double theta3 = OpA * N; // cote du point externe
OpAp = OpA - theta3 * N; // projection de OpA sur la facette
double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2); // les coordonnées locales de A
// construction de la ligne de la matrice
mat(ncot,1) = theta1 * (theta1 - 1.);
mat(ncot,2) = theta2 * (theta2 - 1.);
mat(ncot,3) = theta1 * theta2;
// et du second membre correspondant
sm(ncot) = theta3;
}
else // cas ou il n'y a pas de noeud externe
{// On va calculer la normale non normée au coté externe
// ce vecteur est dans le plan de la facette
CB = B - C;
U1 = CB ;
V = Util::ProdVec_coor(N,U1);
// maintenant on calcul les coordonnées locales de V
double theta1 = V * aiH.Coordo(1); double theta2 = V * aiH.Coordo(2);
// construction de la ligne de la matrice
mat(ncot,1) = theta1 * theta1;
mat(ncot,2) = theta2 * theta2;
mat(ncot,3) = theta1 * theta2;
// et du second membre correspondant
sm(ncot) = 0;
};
}
else
{ // cas avec une condition de tangente imposée
// on commence par définir la direction de la tangente
CB = B - C;
V = Util::ProdVec_coor(CB,vplan(ncot % 3 +1));
// calcul des coordonnées locales du vecteur tangent
double d1 = V*aiH.Coordo(1); double d2=V*aiH.Coordo(2); double d3 = V*N;
// on calcul les coordonnées du point milieu de BC que l'on nome Mm
Mm = 0.5*(B+C);
OpMm = Mm - Op;
// les coordonnées locales de Mm
double theta1 = OpMm * aiH.Coordo(1); double theta2 = OpMm * aiH.Coordo(2);
// construction de la ligne de la matrice
mat(ncot,1) = (2.*theta1 - 1.) * d1;
mat(ncot,2) = (2.*theta2 - 1.) * d2;
mat(ncot,3) = (theta1 * d2 + theta2 * d1);
// et du second membre correspondant
sm(ncot) = d3;
};
};
// calcul du tenseur de courbure dans la base ai
Mat_pleine mati = mat.Inverse();
Vecteur curbi = mati * sm;
// retour de la courbure
Tenseur_ns2BB curb;
curb.Coor(1,1)=2. * curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(3); curb.Coor(2,2)= 2. * curbi(2);
return curb;
};
// routine (4) 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 4: calcul de la courbure à partir d'un polynôme passant par les 6 points
void Met_Sfe1::Dcourbure4
( 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 & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ int nbddl = (nbNoeudCentral+3)*3; //18; // nombre de ddl total
Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
Tableau <Mat_pleine> dmat(nbddl,mat); // sa variation
Vecteur sm(3); // le second membre pour le calcul de la courbure
Tableau <Vecteur> dsm(nbddl,sm); // sa variation
// 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 : 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);
};
// 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 autres dérivées sont nulles
// 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 CB,U1,V,dV;
Coordonnee3 OpA,OpAp;
Coordonnee3 dA,dOpA,dOpAp,dB,dC;
Coordonnee3 Mm,dMm,OpMm,dOpMm;
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// 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);
// examen de condition limite éventuelle
int nbcote = ncot % 3 +1; // le numéro du coté en face du noeud ncot
if (tabTypeCL(nbcote) != TANGENTE_CL)
{ // cas sans conditions limites particulières
if (D != A) // cas ou le noeud externe existe
{OpA = A - Op;
double theta3 = OpA * N; // cote du point externe
OpAp = OpA - theta3 * N; // projection de OpA sur la facette
// les coordonnées locales de A
double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2);
// construction de la ligne de la matrice
mat(ncot,1) = theta1 * (theta1 - 1.);
mat(ncot,2) = theta2 * (theta2 - 1.);
mat(ncot,3) = theta1 * theta2;
// et du second membre correspondant
sm(ncot) = theta3;
// calcul des variations
for (int ia=1;ia<= 3;ia++)
{{// indice de A: jnA, noeud extérieur
int jnA = 9+(ncot-1)*3+ia;
dA = Ia(ia);
dOpA = dA; // dOp(jnA)=0 // OpA = A - Op;
// la coordonnée en 3
double dtheta3 = dOpA * N; // dN(jnA)=0 car noeud externe // theta3 = OpA * N
dOpAp = dOpA - dtheta3 * N ; //dN(jnA)=0; // OpAp = OpA - theta3 * N
// les coordonnées locales
double dtheta1 = dOpAp * aiH.Coordo(1); // DaiH(jnA).Coordo(1)=0; // theta1 = OpAp * aiH.Coordo(1)
double dtheta2 = dOpAp * aiH.Coordo(2);// DaiH(jnA).Coordo(2)=0; // theta2 = OpAp * aiH.Coordo(2)
// la matrice
dmat(jnA)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
dmat(jnA)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
dmat(jnA)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
// le second membre
dsm(jnA)(ncot) = dtheta3; // sm(ncot) = theta3
}
{// indice de B: jnB, noeud de la facette centrale
int jnB = (indi(ncot+1)-1)*3+ia;
// dA = 0 ;
dOpA = - dOp(jnB); // OpA = A - Op;
// la coordonnée en 3
double dtheta3 = dOpA * N + OpA * dN(jnB); // theta3 = OpA * N
dOpAp = dOpA - dtheta3 * N - theta3 * dN(jnB); // OpAp = OpA - theta3 * N
// les coordonnées locales
// 1) à noter ici que OpAp est un vecteur de la facette, donc normale à N, et que DaiH ne contient que les
// variations plane de aiH (il manque les variations suivant la normale) mais cette partie suffit car la partie
// qui manque est suivant N et donc normale à OpAp (donc elle disparait au niveau du produit scalaire)
double dtheta1 = dOpAp * aiH.Coordo(1) + OpAp * DaiH(jnB).Coordo(1); // theta1 = OpAp * aiH.Coordo(1)
double dtheta2 = dOpAp * aiH.Coordo(2) + OpAp * DaiH(jnB).Coordo(2); // theta2 = OpAp * aiH.Coordo(2)
// la matrice
dmat(jnB)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
dmat(jnB)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
dmat(jnB)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
// le second membre
dsm(jnB)(ncot) = dtheta3; // sm(ncot) = theta3
}
{// indice de C: jnC
int jnC = (indi(ncot+2)-1)*3+ia;
// dA = 0 ;
dOpA = - dOp(jnC); // OpA = A - Op;
// la coordonnée en 3
double dtheta3 = dOpA * N + OpA * dN(jnC); // theta3 = OpA * N
dOpAp = dOpA - dtheta3 * N - theta3 * dN(jnC); // OpAp = OpA - theta3 * N
// les coordonnées locales
// mêmes remarques 1) que plus haut
double dtheta1 = dOpAp * aiH.Coordo(1) + OpAp * DaiH(jnC).Coordo(1); // theta1 = OpAp * aiH.Coordo(1)
double dtheta2 = dOpAp * aiH.Coordo(2) + OpAp * DaiH(jnC).Coordo(2); // theta2 = OpAp * aiH.Coordo(2)
// la matrice
dmat(jnC)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
dmat(jnC)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
dmat(jnC)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
// le second membre
dsm(jnC)(ncot) = dtheta3; // sm(ncot) = theta3
}
{// indice de D: jnD
int jnD = (indi(ncot)-1)*3+ia;
// dA = 0 ;
dOpA = - dOp(jnD); // OpA = A - Op;
// la coordonnée en 3
double dtheta3 = dOpA * N + OpA * dN(jnD); // theta3 = OpA * N
dOpAp = dOpA - dtheta3 * N - theta3 * dN(jnD); // OpAp = OpA - theta3 * N
// les coordonnées locales
// mêmes remarques 1) que plus haut
double dtheta1 = dOpAp * aiH.Coordo(1) + OpAp * DaiH(jnD).Coordo(1); // theta1 = OpAp * aiH.Coordo(1)
double dtheta2 = dOpAp * aiH.Coordo(2) + OpAp * DaiH(jnD).Coordo(2); // theta2 = OpAp * aiH.Coordo(2)
// la matrice
dmat(jnD)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
dmat(jnD)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
dmat(jnD)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
// le second membre
dsm(jnD)(ncot) = dtheta3; // sm(ncot) = theta3
}
}; // fin du for sur ia
}
else // cas ou il n'y a pas de noeud externe
{ // On va calculer la normale non normée au coté externe
// ce vecteur est dans le plan de la facette
CB = B - C;
U1 = CB ;
V = Util::ProdVec_coor(N,U1);
// maintenant on calcul les coordonnées locales de V
double theta1 = V * aiH.Coordo(1); double theta2 = V * aiH.Coordo(2);
// construction de la ligne de la matrice
mat(ncot,1) = theta1 * theta1;
mat(ncot,2) = theta2 * theta2;
mat(ncot,3) = theta1 * theta2;
// et du second membre correspondant
sm(ncot) = 0;
// calcul des variations
for (int ia=1;ia<= 3;ia++)
{
// {// indice de A: jnA, noeud extérieur
// comme il n'y a pas de noeud extérieur et que N ne dépend pas des ddl extérieur,
// normalement jnA n'intervient pas !!
// int jnA = 9+(ncot-1)*3+ia;
// // dCB = 0.; // CB = B - C
// // dU1 = 0.; // U1 = CB
// dV = Util::ProdVec_coor(dN(jnA),U1); // V = Util::ProdVec_coor(N,U1)
// double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnA).Coordo(1); // theta1 = V * aiH.Coordo(1);
// double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnA).Coordo(2); // theta2 = V * aiH.Coordo(2);
// dmat(jnA)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
// dmat(jnA)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
// dmat(jnA)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
// dsm(jnA)(ncot) = 0.; // sm(ncot) = 0
// }
{// indice de B: jnB, noeud de la facette centrale
int jnB = (indi(ncot+1)-1)*3+ia;
dB = Ia(ia);
// d_CB = dB // CB = B - C
// dU1 = dB; // U1 = CB
dV = Util::ProdVec_coor(dN(jnB),U1) + Util::ProdVec_coor(N,dB); // V = Util::ProdVec_coor(N,U1)
// 2) à noter ici que V est un vecteur de la facette, donc normale à N, et que DaiH ne contient que les
// variations plane de aiH (il manque les variations suivant la normale) mais cette partie suffit car la partie
// qui manque est suivant N et donc normale à V (donc elle disparait au niveau du produit scalaire)
double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnB).Coordo(1); // theta1 = V * aiH.Coordo(1);
double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnB).Coordo(2); // theta2 = V * aiH.Coordo(2);
dmat(jnB)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
dmat(jnB)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
dmat(jnB)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
dsm(jnB)(ncot) = 0.; // sm(ncot) = 0
}
{// indice de C: jnC
int jnC = (indi(ncot+2)-1)*3+ia;
dC = Ia(ia);
//dCB = -dC; // CB = B - C
// dU1 = -dC; // U1 = CB
dV = Util::ProdVec_coor(dN(jnC),U1) + Util::ProdVec_coor(N,-dC); // V = Util::ProdVec_coor(N,U1)
// mêmes remarques 2) que plus haut
double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnC).Coordo(1); // theta1 = V * aiH.Coordo(1);
double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnC).Coordo(2); // theta2 = V * aiH.Coordo(2);
dmat(jnC)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
dmat(jnC)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
dmat(jnC)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
dsm(jnC)(ncot) = 0.; // sm(ncot) = 0
}
{// indice de D: jnD
int jnD = (indi(ncot)-1)*3+ia;
// dCB = 0.; // CB = B - C
// dU1 = 0.; // U1 = CB
dV = Util::ProdVec_coor(dN(jnD),U1); // V = Util::ProdVec_coor(N,U1)
// mêmes remarques 2) que plus haut
double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnD).Coordo(1); // theta1 = V * aiH.Coordo(1);
double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnD).Coordo(2); // theta2 = V * aiH.Coordo(2);
dmat(jnD)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
dmat(jnD)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
dmat(jnD)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
dsm(jnD)(ncot) = 0.; // sm(ncot) = 0
}
}; // fin du for sur ia
}; // fin du cas où il n'y a pas de noeud externe
}
else
{ // cas avec une condition de tangente imposée
// on commence par définir la direction de la tangente
CB = B - C;
V = Util::ProdVec_coor(CB,vplan(nbcote));
// on calcul les coordonnées du point milieu de BC que l'on nome Mm
Mm = 0.5*(B+C);
OpMm = Mm - Op;
// calcul des coordonnées locales du vecteur tangent
double d1 = V*aiH.Coordo(1); double d2=V*aiH.Coordo(2); double d3 = V*N;
// les coordonnées locales de Mm
double theta1 = OpMm * aiH.Coordo(1); double theta2 = OpMm * aiH.Coordo(2);
// les coordonnée sont ici des coordonnées matérielles qui dans le cas linéaires sont constantes (0.5,0) ou
// (0.5,0.5) ou (0,0.5) leur variations est donc nulles
// ***** donc on peut simplifier les expressions qui suivent en annulant les variations de theta_alpha
// construction de la ligne de la matrice
mat(ncot,1) = (2.*theta1 - 1.) * d1;
mat(ncot,2) = (2.*theta2 - 1.) * d2;
mat(ncot,3) = (theta1 * d2 + theta2 * d1);
// et du second membre correspondant
sm(ncot) = d3;
// calcul des variations
for (int ia=1;ia<= 3;ia++)
{ {// indice de B: jnB, noeud de la facette centrale
int jnB = (indi(ncot+1)-1)*3+ia;
dB = Ia(ia);
// d_CB = dB // CB = B - C
dV = Util::ProdVec_coor(dB,vplan(nbcote)); // V = Util::ProdVec_coor(CB,vplan(nbcote))
// mêmes remarques 2) que plus haut
double dd1=dV*aiH.Coordo(1)+V*DaiH(jnB).Coordo(1); // d1 = V*aiH.Coordo(1)
double dd2=dV*aiH.Coordo(2)+V*DaiH(jnB).Coordo(2); // d2=V*aiH.Coordo(2)
double dd3=dV*N+V*dN(jnB); // d3 = V*N
// dMm = 0.5 * dB ; // Mm = 0.5*(B+C);
dOpMm = 0.5 * dB - dOp(jnB); // OpMm = Mm - Op;
// les coordonnées locales
// ici également, OpMm est dans le plan de la facette d'où mêmes remarques 2) que plus haut
double dtheta1 = dOpMm * aiH.Coordo(1) + OpMm * DaiH(jnB).Coordo(1); // theta1 = OpMm * aiH.Coordo(1)
double dtheta2 = dOpMm * aiH.Coordo(2) + OpMm * DaiH(jnB).Coordo(2); // theta2 = OpMm * aiH.Coordo(2)
dmat(jnB)(ncot,1) = (2.*dtheta1-1.)*d1+(2.*theta1-1.)*dd1; // mat(ncot,1) = (2.*theta1 - 1.) * d1
dmat(jnB)(ncot,2) = (2.*dtheta2-1.)*d2+(2.*theta2-1.)*dd2; // mat(ncot,1) = (2.*theta1 - 1.) * d1
dmat(jnB)(ncot,3) = (dtheta1*d2+dtheta2*d1)+(theta1*dd2+theta2*dd1); // mat(ncot,3) = (theta1 * d2 + theta2 * d1);
dsm(jnB)(ncot) = dd3; // sm(ncot) = d3
}
{// indice de C: jnC
int jnC = (indi(ncot+2)-1)*3+ia;
dC = Ia(ia);
//dCB = -dC; // CB = B - C
dV = Util::ProdVec_coor(-dC,vplan(nbcote)); // V = Util::ProdVec_coor(CB,vplan(nbcote))
double dd1=dV*aiH.Coordo(1)+V*DaiH(jnC).Coordo(1); // d1 = V*aiH.Coordo(1)
double dd2=dV*aiH.Coordo(2)+V*DaiH(jnC).Coordo(2); // d2=V*aiH.Coordo(2)
double dd3=dV*N+V*dN(jnC); // d3 = V*N
// dMm = 0.5 * dC ; // Mm = 0.5*(B+C);
dOpMm = 0.5 * dC - dOp(jnC); // OpMm = Mm - Op;
// les coordonnées locales
// ici également, OpMm est dans le plan de la facette d'où mêmes remarques 2) que plus haut
double dtheta1 = dOpMm * aiH.Coordo(1) + OpMm * DaiH(jnC).Coordo(1); // theta1 = OpMm * aiH.Coordo(1)
double dtheta2 = dOpMm * aiH.Coordo(2) + OpMm * DaiH(jnC).Coordo(2); // theta2 = OpMm * aiH.Coordo(2)
dmat(jnC)(ncot,1) = (2.*dtheta1-1.)*d1+(2.*theta1-1.)*dd1; // mat(ncot,1) = (2.*theta1 - 1.) * d1
dmat(jnC)(ncot,2) = (2.*dtheta2-1.)*d2+(2.*theta2-1.)*dd2; // mat(ncot,1) = (2.*theta1 - 1.) * d1
dmat(jnC)(ncot,3) = (dtheta1*d2+dtheta2*d1)+(theta1*dd2+theta2*dd1); // mat(ncot,3) = (theta1 * d2 + theta2 * d1);
dsm(jnC)(ncot) = dd3; // sm(ncot) = d3
}
{// indice de D: jnD
int jnD = (indi(ncot)-1)*3+ia;
// dCB = 0.; // CB = B - C
dV = 0. ; // V = Util::ProdVec_coor(CB,vplan(nbcote))
double dd1=V*DaiH(jnD).Coordo(1); // d1 = V*aiH.Coordo(1)
double dd2=V*DaiH(jnD).Coordo(2); // d2=V*aiH.Coordo(2)
double dd3=V*dN(jnD); // d3 = V*N
// dMm = 0.5 * dB ; // Mm = 0.5*(B+C);
dOpMm = - dOp(jnD); // OpMm = Mm - Op;
// les coordonnées locales
// ici également, OpMm est dans le plan de la facette d'où mêmes remarques 2) que plus haut
double dtheta1 = dOpMm * aiH.Coordo(1) + OpMm * DaiH(jnD).Coordo(1); // theta1 = OpMm * aiH.Coordo(1)
double dtheta2 = dOpMm * aiH.Coordo(2) + OpMm * DaiH(jnD).Coordo(2); // theta2 = OpMm * aiH.Coordo(2)
dmat(jnD)(ncot,1) = (2.*dtheta1-1.)*d1+(2.*theta1-1.)*dd1; // mat(ncot,1) = (2.*theta1 - 1.) * d1
dmat(jnD)(ncot,2) = (2.*dtheta2-1.)*d2+(2.*theta2-1.)*dd2; // mat(ncot,1) = (2.*theta1 - 1.) * d1
dmat(jnD)(ncot,3) = (dtheta1*d2+dtheta2*d1)+(theta1*dd2+theta2*dd1); // mat(ncot,3) = (theta1 * d2 + theta2 * d1);
dsm(jnD)(ncot) = dd3; // sm(ncot) = d3
}
}; // fin du for sur ia
}; // fin du test sur les CL: condition de tangence
}; // fin de la boucle sur les 3 noeuds externes
// calcul du tenseur de courbure dans la base ai
Mat_pleine mati = mat.Inverse();
//---- pour le debug des cas particuliers ----
#ifdef MISE_AU_POINT
if (ParaGlob::NiveauImpression() > 3)
{ int i,j;
if (mati.MaxiValAbs(i,j) > ConstMath::grand)
{ cout << "\n *** Met_Sfe1::Dcourbure4, mati.MaxiValAbs(i,j) est trop grand :\n mat= ";
mat.Affiche();
cout << "\n inverse de mat: mati= ";
mati.Affiche();
// affichage des éléments du triangle
cout << "\n pour chaque cote : on utilise les notations : triangle principal DBC et secondaire BAC";
cout << "\n Op= "<< Op << "\n N= "<< N ;
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// simplification de l'ecriture
cout << "\n cote : "<< ncot ;
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);
cout << "\n D= "<< D << "\n B= " << B << "\n C= "<< C << "\n A= "<< A;
cout << "\n aiH(1)= "<<aiH.Coordo(1) << "\n aiH(2)= "<<aiH.Coordo(2);
int nbcote = ncot % 3 +1; // le numéro du coté en face du noeud ncot
if (tabTypeCL(nbcote) != TANGENTE_CL) // seul cas envisagé pour l'instant
{ // cas sans conditions limites particulières
if (D != A) // cas ou le noeud externe existe
{OpA = A - Op;
double theta3 = OpA * N; // cote du point externe
OpAp = OpA - theta3 * N; // projection de OpA sur la facette
cout << "\n OpA = A - Op= "<< OpA
<< "\n OpAp = OpA - theta3 * N= " << OpAp ;
// les coordonnées locales de A
double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2);
cout << "\n theta1= OpAp * aiH(1)= "<< theta1 << " theta2= OpAp * aiH(2)= " << theta2 << " ";
cout << "\n mat(ncot,1) = theta1 * (theta1 - 1.) ";
cout << "\n mat(ncot,2) = theta2 * (theta2 - 1.)";
cout << "\n mat(ncot,3) = theta1 * theta2";
cout << "\n sm(ncot) = theta3";
};
};
};
cout << endl;
};
};
#endif
//---- fin pour le debug des cas particuliers ----
Vecteur curbi = mati * sm;
// retour de la courbure
curb.Coor(1,1)=2. * curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(3); curb.Coor(2,2)= 2. * curbi(2);
// idem pour les variations
Coordonnee3 dcurbi;
// Mat_pleine dmata(mat); // une matrice intermédiaire
for (int ddl=1;ddl<=nbddl;ddl++)
{ dcurbi = mati * (dsm(ddl) - dmat(ddl) * curbi );
Tenseur_ns2BB& dcur = dcurb(ddl);
// ----- pour le débug
// cout << "\n dcurbi(1),(" << ddl << ")= " << dcurbi(1)
// << " dcurbi(2),(" << ddl << ")= " << dcurbi(2)
// << " dcurbi(3),(" << ddl << ")= " << dcurbi(3) << endl;
// ----- fin pour le débug
dcur.Coor(1,1)=2. * dcurbi(1); dcur.Coor(1,2)=dcur.Coor(2,1)=dcurbi(3); dcur.Coor(2,2)= 2. * dcurbi(2);
};
// ----- pour le débug
// Sortie(1);
// ----- fin pour le débug
};
// 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_anormale4(double inf_normale, const Tableau<Coordonnee3>& tab_coor
, const BaseB & aiB , const BaseH & aiH)
{ // calcul de la normale centrale
// on utilise des expressions strictement identique ici et dans dcourbure
// de manière à éviter les erreurs d'arrondis
//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;
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
Vecteur sm(3); // le second membre pour le calcul de la courbure
// 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
Coordonnee3 U1,CB,V,OpA;
Coordonnee3 Mm,OpMm;
Coordonnee3 OpAp; // le vecteur horizontal correspondant à OpA
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// 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);
// pas d'examen de condition limite éventuelle
// on ne considère que le cas sans conditions limites particulières
// ncot % 3 +1: = le coté en face du noeud ncot
{ // cas sans conditions limites particulières
if (D != A) // cas ou le noeud externe existe
{OpA = A - Op;
double theta3 = OpA * N; // cote du point externe
OpAp = OpA - theta3 * N; // projection de OpA sur la facette
double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2); // les coordonnées locales de A
// construction de la ligne de la matrice
mat(ncot,1) = theta1 * (theta1 - 1.);
mat(ncot,2) = theta2 * (theta2 - 1.);
mat(ncot,3) = theta1 * theta2;
// et du second membre correspondant
sm(ncot) = theta3;
}
else // cas ou il n'y a pas de noeud externe
{// On va calculer la normale non normée au coté externe
// ce vecteur est dans le plan de la facette
CB = B - C;
U1 = CB ;
V = Util::ProdVec_coor(N,U1);
// maintenant on calcul les coordonnées locales de V
double theta1 = V * aiH.Coordo(1); double theta2 = V * aiH.Coordo(2);
// construction de la ligne de la matrice
mat(ncot,1) = theta1 * theta1;
mat(ncot,2) = theta2 * theta2;
mat(ncot,3) = theta1 * theta2;
// et du second membre correspondant
sm(ncot) = 0;
};
}
};
// le principe du test est de regarder les valeurs de mat qui doit-être une matrice
// non singulière
double det = mat.Determinant();
//--- debug
#ifdef MISE_AU_POINT
//cout << "\n debug Met_Sfe1::Test_courbure_anormale4 "
// << " determinant= "<< det;
//if (Dabs(det) < ConstMath::petit)
if ((Dabs(det) < inf_normale)&& (ParaGlob::NiveauImpression() > 5))
{ cout << "\n *** Met_Sfe1::Test_courbure_anormale4, determinant trop petit = "<< det
<< " \n mat= ";
mat.Affiche();
// affichage des éléments du triangle
cout << "\n pour chaque cote : on utilise les notations : triangle principal DBC et secondaire BAC";
cout << "\n Op= "<< Op << "\n N= "<< N ;
for (int ncot=1;ncot<=3;ncot++)
{ // on utilise les notations : triangle principal DBC et secondaire BAC
// simplification de l'ecriture
cout << "\n cote : "<< ncot ;
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);
cout << "\n D= "<< D << "\n B= " << B << "\n C= "<< C << "\n A= "<< A;
cout << "\n aiH(1)= "<<aiH.Coordo(1) << "\n aiH(2)= "<<aiH.Coordo(2);
int nbcote = ncot % 3 +1; // le numéro du coté en face du noeud ncot
{ // cas sans conditions limites particulières
if (D != A) // cas ou le noeud externe existe
{OpA = A - Op;
double theta3 = OpA * N; // cote du point externe
OpAp = OpA - theta3 * N; // projection de OpA sur la facette
cout << "\n OpA = A - Op= "<< OpA
<< "\n OpAp = OpA - theta3 * N= " << OpAp ;
// les coordonnées locales de A
double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2);
cout << "\n theta1= OpAp * aiH(1)= "<< theta1 << " theta2= OpAp * aiH(2)= " << theta2 << " ";
cout << "\n mat(ncot,1) = theta1 * (theta1 - 1.) ";
cout << "\n mat(ncot,2) = theta2 * (theta2 - 1.)";
cout << "\n mat(ncot,3) = theta1 * theta2";
cout << "\n sm(ncot) = theta3";
};
};
};
cout << endl;
};
#endif
//--- fin debug
int retour = 1; // init par défaut
// if (det < ConstMath::petit)
if (Dabs(det) < inf_normale)
retour = 0;
// retour de l'indicateur
return retour;
};
// programme de vérification par différences finies des dérivées
void Met_Sfe1::Verif_Dcourbure
( const Tableau<Noeud *>& tab_noeud,bool gradV_instantane
,Mat_pleine const & dphiS,int nombre_noeud,Vecteur const & phiS
,Mat_pleine const & tabdphiH,Vecteur const& phiH,const Epai* epas
,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
{ // l'idée est de faire une vérification des dérivées à l'aide d'une méthode de différence finie
int dim = ParaGlob::Dimension();
// dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_DcourbureSFE == 0)
{ cout << "\n ****verification du calcul de la variation de courbure ****";
cout << "\n Met_Sfe1::Verif_Dcourbure(.. \n";
}
indic_Verif_DcourbureSFE++;
// on définit 2 seconde métriques, car il faudra modifier les vecteurs ai et les vecteurs gi, donc
// il vaut mieux travailler sur une copie
Met_Sfe1 metrique_bis_1(*this);Met_Sfe1 metrique_bis_2(*this);
// ici on considère que l'on a même nombre de ddl par noeud = dim
// on va modifier chaque ddl de chaque noeud systématiquement
int nbnoeud = tab_noeud.Taille();
// le deltat pour les différences finis
double delta = ConstMath::unpeupetit;
double mini_val = ConstMath::pasmalpetit;
int numddl = 1; // le compteur de ddl
int nberr = 10;
Tableau<int> erreur(nberr);
double diff_admi = 0.1;
bool premier_calcul=true;
for (int inoeud=1;inoeud<=nbnoeud;inoeud++)
// on ne fait la vérification que si le noeud externe existe
if (( inoeud <= 3) || ( (tab_noeud(inoeud-3)->Num_noeud()) != (tab_noeud(inoeud)->Num_noeud())))
{// on récupère les coordonnées du noeud
Coordonnee coordtdt = tab_noeud(inoeud)->Coord2();
for (int ix= 1;ix<=dim;ix++,numddl++)
{ Coordonnee X(dim); X(ix) += delta;
tab_noeud(inoeud)->Ajout_coord2(X);
// appel de la métrique
metrique_bis_1.CalSfe1_explicit_tdt( tab_noeud,gradV_instantane,dphiS,nombre_noeud,phiS,premier_calcul
,tabdphiH,phiH,epas,tabTypeCL,vplan);
// maintenant on remet les coordonnées du noeud à l'état initial
tab_noeud(inoeud)->Change_coord2(coordtdt);
// et on décale dégativement
X(ix) -= 2*delta;
tab_noeud(inoeud)->Ajout_coord2(X);
// appel de la métrique
metrique_bis_2.CalSfe1_explicit_tdt( tab_noeud,gradV_instantane,dphiS,nombre_noeud,phiS,premier_calcul
,tabdphiH,phiH,epas,tabTypeCL,vplan);
// maintenant on remet les coordonnées du noeud à l'état initial
tab_noeud(inoeud)->Change_coord2(coordtdt);
premier_calcul=false;
// calcul des dérivées numériques et vérification
int nb_vecteur = 2;
for (int j=1;j<=nb_vecteur;j++)
{ // variation des vecteurs aiB_tdt
CoordonneeB daiB = ((*metrique_bis_1.aiB_tdt)(j) -(*metrique_bis_2.aiB_tdt)(j))/(2*delta);
CoordonneeB toto = ((*metrique_bis_1.aiB_tdt)(j) -(*metrique_bis_2.aiB_tdt)(j));
for (int i=1;i<=dim;i++)
if (diffpourcent(daiB(i),(*d_aiB_tdt)(numddl)(j)(i),MaX(Dabs(daiB(i)),Dabs((*d_aiB_tdt)(numddl)(j)(i))),diff_admi))
{ if (DabsMiN(daiB(i),(*d_aiB_tdt)(numddl)(j)(i)) <= mini_val)
{if ( DabsMaX(daiB(i),(*d_aiB_tdt)(numddl)(j)(i)) > 50.*delta)
{erreur(1) += 1;
cout << "\ndaiB("<<i<<")="<<daiB(i)<<",(*d_aiB_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiB_tdt)(numddl)(j)(i);
cout << "\n aiB(xi+delta)="<<(*metrique_bis_1.aiB_tdt)(j)(i) << " aiB(xi-delta)="<<(*metrique_bis_2.aiB_tdt)(j)(i);
cout << "\n toto= " << toto << " " << (*aiB_tdt)(j) << " " << (*metrique_bis_1.aiB_tdt)(j);
}
}
else
{erreur(2) += 1;
cout << "\ndaiB("<<i<<")="<<daiB(i)<<",(*d_aiB_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiB_tdt)(numddl)(j)(i);
};
};
};
for (int j=1;j<=nb_vecteur;j++)
{ // variation des vecteurs aiH_tdt
CoordonneeH d_aiH = ((*metrique_bis_1.aiH_tdt)(j) -(*metrique_bis_2.aiH_tdt)(j))/(2*delta);
Coordonnee daiH = d_aiH.Coor_const(); // un coordonné à la même place
// en fait la variation de aiH calculée est uniquement celle dans le plan de la facette, il faut donc
// retirer la partie hors facette pour la comparaison
Met_Sfe1::Calcul_N_tdt(); // calcul de la normale
daiH -= ((d_aiH * N_tdt) * N_tdt).Coor_const(); // on retire la partie portée par la normale
for (int i=1;i<=dim;i++)
if (diffpourcent(daiH(i),(*d_aiH_tdt)(numddl)(j)(i),MaX(Dabs(daiH(i)),Dabs((*d_aiH_tdt)(numddl)(j)(i))),diff_admi))
{if (DabsMiN(daiH(i),(*d_aiH_tdt)(numddl)(j)(i)) <= mini_val)
{if ( DabsMaX(daiH(i),(*d_aiH_tdt)(numddl)(j)(i)) > 50.*delta)
{erreur(3) += 1;
cout << "\ndaiH("<<i<<")="<<daiH(i)<<",(*d_aiH_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiH_tdt)(numddl)(j)(i);
cout << "\n aiH_(xi+delta)="<<(*metrique_bis_1.aiH_tdt)(j)(i) << " aiH_(xi-delta)="<<(*metrique_bis_2.aiH_tdt)(j)(i);
}
}
else
{erreur(4) += 1;
cout << "\ndaiH("<<i<<")="<<daiH(i)<<",(*d_aiH_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiH_tdt)(numddl)(j)(i);
};
};
};
// on ne peut pas faire la vérif des métriques car elles ne sont pas calculées !!
/*
// variation du tenseur aijBB_tdt
TenseurBB * aijBB = NevezTenseurBB(*metrique_bis_1.aijBB_tdt);
*aijBB = (*metrique_bis_1.aijBB_tdt - *metrique_bis_2.aijBB_tdt) / (2*delta);
for (int i1=1;i1<=nb_vecteur;i1++)
for (int j1=1;j1<=nb_vecteur;j1++)
if (diffpourcent((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1),
DabsMaX((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1)),diff_admi))
if (DabsMiN((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1)) <= mini_val)
{if( DabsMaX((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1)) > 50.*delta)
{erreur(5) += 1;
cout << "\ndaijBB("<<i1<<","<<j1<<")="<<(*aijBB)(i1,j1)<<",(*d_aijBB_tdt)("<<numddl<<")("<<i1<<","<<j1<<")="<< (*(d_aijBB_tdt(numddl)))(i1,j1);
cout << "\n(*d_aiB_tdt)("<<numddl<<")(1)="<< (*d_aiB_tdt)(numddl)(1) << " aiB(1)=" << (*aiB_tdt)(1);
}}
else erreur(6) += 1;
// variation du tenseur aijHH_tdt
TenseurHH * aijHH = NevezTenseurHH(*metrique_bis_1.aijHH_tdt);
*aijHH = (*metrique_bis_1.aijHH_tdt - *metrique_bis_2.aijHH_tdt) / (2*delta);
for (int i1=1;i1<=nb_vecteur;i1++)
for (int j1=1;j1<=nb_vecteur;j1++)
if (diffpourcent((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1),
DabsMaX((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1)),diff_admi))
if (DabsMiN((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1)) <= mini_val)
{if( DabsMaX((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1)) > 50.*delta) erreur(7) += 1;}
else erreur(8) += 1;
*/
// variation du tenseur de courbure
Tenseur_ns2BB d_curb(curb_tdt);
// calcul de la variation de la métrique dans le repère absolu
d_curb = (metrique_bis_1.curb_tdt - metrique_bis_2.curb_tdt) / (2*delta);
// on calcul the same pour la variation calculée
// Tenseur_ns2BB d_curb_glob;
// dcurb_tdt(numddl).BaseAbsolue(d_curb_glob,*aiH_tdt);
for (int i1=1;i1<=nb_vecteur;i1++)
for (int j1=1;j1<=nb_vecteur;j1++)
if (diffpourcent(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1),
DabsMaX(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1)),diff_admi))
{if (DabsMiN(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1)) <= mini_val)
{if( DabsMaX(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1)) > 50.*delta)
{erreur(9) += 1;
cout << "\nd_curb("<<i1<<","<<j1<<")="<<d_curb(i1,j1)<<",dcurb_tdt("<<numddl<<")("<<i1<<","<<j1<<")="<< dcurb_tdt(numddl)(i1,j1);
cout << "\n curbe= " << curb_tdt << ",dcurb_tdt("<<numddl<<")("<< dcurb_tdt(numddl) ;
}
}
else
{erreur(10) += 1;
cout << "\nd_curb("<<i1<<","<<j1<<")="<<d_curb(i1,j1)<<",dcurb_tdt("<<numddl<<")("<<i1<<","<<j1<<")="<< dcurb_tdt(numddl)(i1,j1);
}
};
// effacement des tenseurs intermédiaires
// delete aijHH; delete aijBB;
};// fin de boucle sur la dimension de coordonnée
} // fin du if: s'il y a un noeud externe
else
{numddl += 3;}; // on passe les ddl
// fin de boucle sur les noeuds
// message d'erreur si besoin
bool une_erreur=false;
for (int l=1;l<=nberr;l++) if (erreur(l) != 0) une_erreur=true;
if (une_erreur)
{ cout << "\n erreur dans le calcul analytique des derivees de la courbure";
cout << "\n Met_Sfe1::Verif_Dcourbure(.."
<< " , numero de passage = " << indic_Verif_DcourbureSFE << " erreur = " << erreur << endl;
// Sortie(1);
}
};