1 line
45 KiB
C++
Executable file
1 line
45 KiB
C++
Executable file
// 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 "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);
|
|
}
|
|
};
|
|
|