Herezh_dev/Elements/Mecanique/SFE/Met_Sfe1s4.cc

// 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) != TANGEN