mtu/src/fem_element2.cpp

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       fem_element2.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:54 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------


#include "gestionversion.h"
#include <math.h>
#include "fem_element2.h"
#include "fem_noeud.h"
#include "mg_element_maillage.h"
#include "ot_mathematique.h"

FEM_ELEMENT2::FEM_ELEMENT2(unsigned long num,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(num,mai)
{
}
FEM_ELEMENT2::FEM_ELEMENT2(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo):FEM_ELEMENT_MAILLAGE(num,topo)
{
}
FEM_ELEMENT2::FEM_ELEMENT2(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(num,topo,mai)
{
}

FEM_ELEMENT2::FEM_ELEMENT2(class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(mai)
{
}
FEM_ELEMENT2::FEM_ELEMENT2(class MG_ELEMENT_TOPOLOGIQUE* topo):FEM_ELEMENT_MAILLAGE(topo)
{
}
FEM_ELEMENT2::FEM_ELEMENT2(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(topo,mai)
{
}

FEM_ELEMENT2::FEM_ELEMENT2(FEM_ELEMENT2& mdd):FEM_ELEMENT_MAILLAGE(mdd)
{
}

int FEM_ELEMENT2::get_etat(void)
{
    return etat;  
}
void FEM_ELEMENT2::change_etat(int num)
{
    etat=num;  
}

FEM_ELEMENT2::~FEM_ELEMENT2()
{
}


void FEM_ELEMENT2::extrapoler_solution_noeud(void)
{
    int nb=get_nb_fem_noeud();
    for (int i=0;i<nb;i++)
        for (int j=0;j<MAX_TYPE_SOLUTION;j++)  
          get_fem_noeud(i)->change_solution(solution[j],j);
}

bool FEM_ELEMENT2::get_param_element_fini_2D(double *xyz,double *uv)
{
OT_MATRICE_3D mat;
OT_VECTEUR_3D vec;
uv[0]=0.;
uv[1]=0.;

int compteur=0;
int ok=0;

while (ok==0)
{
mat(0,0)=0.;mat(0,1)=0.;mat(0,2)=0.;
mat(1,0)=0.;mat(1,1)=0.;mat(1,2)=0.;
mat(2,0)=0.;mat(2,1)=0.;mat(2,2)=0.;
vec(0)=xyz[0];vec(1)=xyz[1];vec(2)=xyz[2];

for (int i=0;i<get_nb_fem_noeud();i++)
  {
  mat(0,0)=mat(0,0)+get_fonction_derive_interpolation(i+1,1,uv)*get_fem_noeud(i)->get_x();
  mat(0,1)=mat(0,1)+get_fonction_derive_interpolation(i+1,2,uv)*get_fem_noeud(i)->get_x();
  mat(1,0)=mat(1,0)+get_fonction_derive_interpolation(i+1,1,uv)*get_fem_noeud(i)->get_y();
  mat(1,1)=mat(1,1)+get_fonction_derive_interpolation(i+1,2,uv)*get_fem_noeud(i)->get_y();
  mat(2,0)=mat(2,0)+get_fonction_derive_interpolation(i+1,1,uv)*get_fem_noeud(i)->get_z();
  mat(2,1)=mat(2,1)+get_fonction_derive_interpolation(i+1,2,uv)*get_fem_noeud(i)->get_z();
  vec(0)=vec(0)-get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_x();
  vec(1)=vec(1)-get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_y();
  vec(2)=vec(2)-get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_z();
  }
OT_VECTEUR_3D vtmp1(mat(0,0),mat(1,0),mat(2,0));
OT_VECTEUR_3D vtmp2(mat(0,1),mat(1,1),mat(2,1));
OT_VECTEUR_3D vtmp3=vtmp1&vtmp2;
mat(0,2)=vtmp3.get_x();
mat(1,2)=vtmp3.get_y();
mat(2,2)=vtmp3.get_z();
double det=mat.get_determinant();
if (fabs(det)<1e-12) return false;

OT_MATRICE_3D mat1(vec,vtmp2,vtmp3);
OT_MATRICE_3D mat2(vtmp1,vec,vtmp3);
double d1=mat1.get_determinant()/det;
double d2=mat2.get_determinant()/det;
compteur++;
if (compteur>100) return false;
if ((fabs(d1)<1e-8)&&(fabs(d2)<1e-8))  ok=1;
uv[0]=uv[0]+d1;
uv[1]=uv[1]+d2;
}
return true;
  
}


void FEM_ELEMENT2::get_interpolation_xyz(double* uv, double* xyz)
{
xyz[0]=0.;xyz[1]=0.;xyz[2]=0;
for (int i=0;i<get_nb_fem_noeud();i++)
  {
  xyz[0]=xyz[0]+get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_x();
  xyz[1]=xyz[1]+get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_y();
  xyz[2]=xyz[2]+get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_z();
  }
}

double FEM_ELEMENT2::get_jacobien(double* jac,double *uv,double unite)
{   
    int nb=get_nb_fem_noeud();
      
    OT_MATRICE_3D jacobien;
    for (int i=0;i<2;i++)
    {
     jacobien(i,0)=0.;jacobien(i,1)=0.;jacobien(i,2)=0.;
      for (int k=0;k<nb;k++)
          {
      double valderiv=get_fonction_derive_interpolation(k+1,i+1,uv);
      double *xyz=get_fem_noeud(k)->get_coord();
          jacobien(i,0)=jacobien(i,0)+valderiv*xyz[0]*unite;
          jacobien(i,1)=jacobien(i,1)+valderiv*xyz[1]*unite;
          jacobien(i,2)=jacobien(i,2)+valderiv*xyz[2]*unite;
          }
    }
    OT_VECTEUR_3D j1(jacobien(0,0),jacobien(0,1),jacobien(0,2));
    OT_VECTEUR_3D j2(jacobien(1,0),jacobien(1,1),jacobien(1,2));
    OT_VECTEUR_3D j3=j1&j2;
    j3.norme();
    jacobien(2,0)=j3(0);
    jacobien(2,1)=j3(1);
    jacobien(2,2)=j3(2);
    
    jac[0]=jacobien(0,0);
    jac[1]=jacobien(0,1);
    jac[2]=jacobien(0,2);
    
    jac[3]=jacobien(1,0);
    jac[4]=jacobien(1,1);
    jac[5]=jacobien(1,2);
    
    jac[6]=jacobien(2,0);
    jac[7]=jacobien(2,1);
    jac[8]=jacobien(2,2);
    
    double det=jacobien.get_determinant();
    return det;
}

void FEM_ELEMENT2::get_inverse_jacob(double* j,double *uv,double unite)
{
    double jac[9];
    get_jacobien(jac,uv,unite);
    OT_MATRICE_3D J;
    J(0,0)=jac[0];
    J(0,1)=jac[1];
    J(0,2)=jac[2];
    J(1,0)=jac[3];
    J(1,1)=jac[4];
    J(1,2)=jac[5];
    J(2,0)=jac[6];
    J(2,1)=jac[7];
    J(2,2)=jac[8];
    OT_MATRICE_3D j_i=J.inverse();
  
    
    j[0]=j_i(0,0);
    j[1]=j_i(0,1);
    j[2]=j_i(0,2);
    
    j[3]=j_i(1,0);
    j[4]=j_i(1,1);
    j[5]=j_i(1,2);
    
    j[6]=j_i(2,0);
    j[7]=j_i(2,1);
    j[8]=j_i(2,2);
}
Revision:	1158
Committed:	Thu Jun 13 22:18:49 2024 UTC (11 months, 1 week ago) by francois
File size:	6141 byte(s)
Log Message:	compatibilité Ubuntu 22.04 Suppression des refeences à Windows Ajout d'une banière
#	User	Rev	Content
1	francois	1158	//####//------------------------------------------------------------
2			//####//------------------------------------------------------------
3			//####// MAGiC
4			//####// Jean Christophe Cuilliere et Vincent FRANCOIS
5			//####// Departement de Genie Mecanique - UQTR
6			//####//------------------------------------------------------------
7			//####// MAGIC est un projet de recherche de l equipe ERICCA
8			//####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9			//####// http://www.uqtr.ca/ericca
10			//####// http://www.uqtr.ca/
11			//####//------------------------------------------------------------
12			//####//------------------------------------------------------------
13			//####//
14			//####// fem_element2.cpp
15			//####//
16			//####//------------------------------------------------------------
17			//####//------------------------------------------------------------
18			//####// COPYRIGHT 2000-2024
19			//####// jeu 13 jun 2024 11:58:54 EDT
20			//####//------------------------------------------------------------
21			//####//------------------------------------------------------------
22	francois	283
23
24			#include "gestionversion.h"
25			#include <math.h>
26	francois	309	#include "fem_element2.h"
27	francois	283	#include "fem_noeud.h"
28			#include "mg_element_maillage.h"
29	francois	628	#include "ot_mathematique.h"
30	francois	283
31	francois	309	FEM_ELEMENT2::FEM_ELEMENT2(unsigned long num,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(num,mai)
32	francois	283	{
33			}
34	francois	378	FEM_ELEMENT2::FEM_ELEMENT2(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo):FEM_ELEMENT_MAILLAGE(num,topo)
35			{
36			}
37			FEM_ELEMENT2::FEM_ELEMENT2(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(num,topo,mai)
38			{
39			}
40	francois	283
41	francois	309	FEM_ELEMENT2::FEM_ELEMENT2(class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(mai)
42	francois	283	{
43			}
44	francois	378	FEM_ELEMENT2::FEM_ELEMENT2(class MG_ELEMENT_TOPOLOGIQUE* topo):FEM_ELEMENT_MAILLAGE(topo)
45			{
46			}
47			FEM_ELEMENT2::FEM_ELEMENT2(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(topo,mai)
48			{
49			}
50	francois	283
51	francois	309	FEM_ELEMENT2::FEM_ELEMENT2(FEM_ELEMENT2& mdd):FEM_ELEMENT_MAILLAGE(mdd)
52	francois	283	{
53			}
54
55	gervaislavoie	385	int FEM_ELEMENT2::get_etat(void)
56			{
57			return etat;
58			}
59			void FEM_ELEMENT2::change_etat(int num)
60			{
61			etat=num;
62			}
63	francois	283
64	francois	309	FEM_ELEMENT2::~FEM_ELEMENT2()
65	francois	283	{
66			}
67
68
69	francois	309	void FEM_ELEMENT2::extrapoler_solution_noeud(void)
70	francois	283	{
71			int nb=get_nb_fem_noeud();
72			for (int i=0;i<nb;i++)
73	francois	377	for (int j=0;j<MAX_TYPE_SOLUTION;j++)
74	francois	375	get_fem_noeud(i)->change_solution(solution[j],j);
75	francois	283	}
76
77	francois	628	bool FEM_ELEMENT2::get_param_element_fini_2D(double xyz,double uv)
78			{
79			OT_MATRICE_3D mat;
80			OT_VECTEUR_3D vec;
81			uv[0]=0.;
82			uv[1]=0.;
83	francois	283
84	francois	628	int compteur=0;
85			int ok=0;
86
87			while (ok==0)
88			{
89			mat(0,0)=0.;mat(0,1)=0.;mat(0,2)=0.;
90			mat(1,0)=0.;mat(1,1)=0.;mat(1,2)=0.;
91			mat(2,0)=0.;mat(2,1)=0.;mat(2,2)=0.;
92			vec(0)=xyz[0];vec(1)=xyz[1];vec(2)=xyz[2];
93
94			for (int i=0;i<get_nb_fem_noeud();i++)
95			{
96			mat(0,0)=mat(0,0)+get_fonction_derive_interpolation(i+1,1,uv)*get_fem_noeud(i)->get_x();
97			mat(0,1)=mat(0,1)+get_fonction_derive_interpolation(i+1,2,uv)*get_fem_noeud(i)->get_x();
98			mat(1,0)=mat(1,0)+get_fonction_derive_interpolation(i+1,1,uv)*get_fem_noeud(i)->get_y();
99			mat(1,1)=mat(1,1)+get_fonction_derive_interpolation(i+1,2,uv)*get_fem_noeud(i)->get_y();
100	francois	1118	mat(2,0)=mat(2,0)+get_fonction_derive_interpolation(i+1,1,uv)*get_fem_noeud(i)->get_z();
101			mat(2,1)=mat(2,1)+get_fonction_derive_interpolation(i+1,2,uv)*get_fem_noeud(i)->get_z();
102	francois	628	vec(0)=vec(0)-get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_x();
103			vec(1)=vec(1)-get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_y();
104	francois	1118	vec(2)=vec(2)-get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_z();
105	francois	628	}
106	francois	1118	OT_VECTEUR_3D vtmp1(mat(0,0),mat(1,0),mat(2,0));
107			OT_VECTEUR_3D vtmp2(mat(0,1),mat(1,1),mat(2,1));
108			OT_VECTEUR_3D vtmp3=vtmp1&vtmp2;
109			mat(0,2)=vtmp3.get_x();
110			mat(1,2)=vtmp3.get_y();
111			mat(2,2)=vtmp3.get_z();
112			double det=mat.get_determinant();
113	francois	628	if (fabs(det)<1e-12) return false;
114	francois	1118
115			OT_MATRICE_3D mat1(vec,vtmp2,vtmp3);
116			OT_MATRICE_3D mat2(vtmp1,vec,vtmp3);
117			double d1=mat1.get_determinant()/det;
118			double d2=mat2.get_determinant()/det;
119	francois	628	compteur++;
120			if (compteur>100) return false;
121			if ((fabs(d1)<1e-8)&&(fabs(d2)<1e-8)) ok=1;
122			uv[0]=uv[0]+d1;
123			uv[1]=uv[1]+d2;
124			}
125			return true;
126
127			}
128
129
130			void FEM_ELEMENT2::get_interpolation_xyz(double* uv, double* xyz)
131			{
132			xyz[0]=0.;xyz[1]=0.;xyz[2]=0;
133			for (int i=0;i<get_nb_fem_noeud();i++)
134			{
135			xyz[0]=xyz[0]+get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_x();
136			xyz[1]=xyz[1]+get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_y();
137			xyz[2]=xyz[2]+get_fonction_interpolation(i+1,uv)*get_fem_noeud(i)->get_z();
138			}
139			}
140
141	francois	1104	double FEM_ELEMENT2::get_jacobien(double* jac,double *uv,double unite)
142			{
143			int nb=get_nb_fem_noeud();
144
145			OT_MATRICE_3D jacobien;
146			for (int i=0;i<2;i++)
147			{
148			jacobien(i,0)=0.;jacobien(i,1)=0.;jacobien(i,2)=0.;
149			for (int k=0;k<nb;k++)
150			{
151			double valderiv=get_fonction_derive_interpolation(k+1,i+1,uv);
152			double *xyz=get_fem_noeud(k)->get_coord();
153			jacobien(i,0)=jacobien(i,0)+valderivxyz[0]unite;
154			jacobien(i,1)=jacobien(i,1)+valderivxyz[1]unite;
155			jacobien(i,2)=jacobien(i,2)+valderivxyz[2]unite;
156			}
157			}
158			OT_VECTEUR_3D j1(jacobien(0,0),jacobien(0,1),jacobien(0,2));
159			OT_VECTEUR_3D j2(jacobien(1,0),jacobien(1,1),jacobien(1,2));
160			OT_VECTEUR_3D j3=j1&j2;
161			j3.norme();
162			jacobien(2,0)=j3(0);
163			jacobien(2,1)=j3(1);
164			jacobien(2,2)=j3(2);
165
166			jac[0]=jacobien(0,0);
167			jac[1]=jacobien(0,1);
168			jac[2]=jacobien(0,2);
169
170			jac[3]=jacobien(1,0);
171			jac[4]=jacobien(1,1);
172			jac[5]=jacobien(1,2);
173
174			jac[6]=jacobien(2,0);
175			jac[7]=jacobien(2,1);
176			jac[8]=jacobien(2,2);
177
178			double det=jacobien.get_determinant();
179			return det;
180			}
181
182			void FEM_ELEMENT2::get_inverse_jacob(double* j,double *uv,double unite)
183			{
184			double jac[9];
185			get_jacobien(jac,uv,unite);
186			OT_MATRICE_3D J;
187			J(0,0)=jac[0];
188			J(0,1)=jac[1];
189			J(0,2)=jac[2];
190			J(1,0)=jac[3];
191			J(1,1)=jac[4];
192			J(1,2)=jac[5];
193			J(2,0)=jac[6];
194			J(2,1)=jac[7];
195			J(2,2)=jac[8];
196			OT_MATRICE_3D j_i=J.inverse();
197
198
199			j[0]=j_i(0,0);
200			j[1]=j_i(0,1);
201			j[2]=j_i(0,2);
202
203			j[3]=j_i(1,0);
204			j[4]=j_i(1,1);
205			j[5]=j_i(1,2);
206
207			j[6]=j_i(2,0);
208			j[7]=j_i(2,1);
209			j[8]=j_i(2,2);
210			}