| 1 |
//------------------------------------------------------------ |
| 2 |
//------------------------------------------------------------ |
| 3 |
// MAGiC |
| 4 |
// Jean Christophe Cuilli�e et Vincent FRANCOIS |
| 5 |
// D�artement de G�ie M�anique - UQTR |
| 6 |
//------------------------------------------------------------ |
| 7 |
// Le projet MAGIC est un projet de recherche du d�artement |
| 8 |
// de g�ie m�anique de l'Universit�du Qu�ec � |
| 9 |
// Trois Rivi�es |
| 10 |
// Les librairies ne peuvent �re utilis�s sans l'accord |
| 11 |
// des auteurs (contact : francois@uqtr.ca) |
| 12 |
//------------------------------------------------------------ |
| 13 |
//------------------------------------------------------------ |
| 14 |
// |
| 15 |
// fem_tetra.cpp |
| 16 |
// |
| 17 |
//------------------------------------------------------------ |
| 18 |
//------------------------------------------------------------ |
| 19 |
// COPYRIGHT 2000 |
| 20 |
// Version du 02/03/2006 �11H22 |
| 21 |
//------------------------------------------------------------ |
| 22 |
//------------------------------------------------------------ |
| 23 |
|
| 24 |
|
| 25 |
#include "gestionversion.h" |
| 26 |
#include <math.h> |
| 27 |
#include "fem_element3.h" |
| 28 |
#include "xfem_element3.h" |
| 29 |
#include "fem_noeud.h" |
| 30 |
#include "fem_maillage.h" |
| 31 |
#include "mg_element_maillage.h" |
| 32 |
#include "ot_boite_3d.h" |
| 33 |
#include "ot_mathematique.h" |
| 34 |
|
| 35 |
FEM_ELEMENT3::FEM_ELEMENT3(unsigned long num,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(num,mai) |
| 36 |
{ |
| 37 |
} |
| 38 |
FEM_ELEMENT3::FEM_ELEMENT3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo):FEM_ELEMENT_MAILLAGE(num,topo) |
| 39 |
{ |
| 40 |
} |
| 41 |
FEM_ELEMENT3::FEM_ELEMENT3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(num,topo,mai) |
| 42 |
{ |
| 43 |
} |
| 44 |
|
| 45 |
FEM_ELEMENT3::FEM_ELEMENT3(class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(mai) |
| 46 |
{ |
| 47 |
} |
| 48 |
|
| 49 |
FEM_ELEMENT3::FEM_ELEMENT3(class MG_ELEMENT_TOPOLOGIQUE* topo):FEM_ELEMENT_MAILLAGE(topo) |
| 50 |
{ |
| 51 |
} |
| 52 |
FEM_ELEMENT3::FEM_ELEMENT3(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai):FEM_ELEMENT_MAILLAGE(topo,mai) |
| 53 |
{ |
| 54 |
} |
| 55 |
|
| 56 |
|
| 57 |
FEM_ELEMENT3::FEM_ELEMENT3(FEM_ELEMENT3& mdd):FEM_ELEMENT_MAILLAGE(mdd) |
| 58 |
{ |
| 59 |
} |
| 60 |
|
| 61 |
|
| 62 |
|
| 63 |
FEM_ELEMENT3::~FEM_ELEMENT3() |
| 64 |
{ |
| 65 |
} |
| 66 |
|
| 67 |
bool FEM_ELEMENT3::get_param_element_fini(double *xyz,double *uvw) |
| 68 |
{ |
| 69 |
OT_MATRICE_3D mat; |
| 70 |
OT_VECTEUR_3D vec; |
| 71 |
uvw[0]=0.; |
| 72 |
uvw[1]=0.; |
| 73 |
uvw[2]=0.; |
| 74 |
|
| 75 |
int compteur=0; |
| 76 |
int ok=0; |
| 77 |
while (ok!=1) |
| 78 |
{ |
| 79 |
ok=0; |
| 80 |
int cas=compteur/100-1; |
| 81 |
int base=4; |
| 82 |
if (compteur!=0) |
| 83 |
{ |
| 84 |
int i=cas%base; |
| 85 |
int j=(cas/base)%base; |
| 86 |
int k=(cas/base/base)%base; |
| 87 |
uvw[0]=-1.+2.*i/(base-1); |
| 88 |
uvw[1]=-1.+2.*j/(base-1); |
| 89 |
uvw[2]=-1.+2.*k/(base-1); |
| 90 |
} |
| 91 |
if (cas>base*base*base) return false; |
| 92 |
while (ok==0) |
| 93 |
{ |
| 94 |
mat(0,0)=0.;mat(0,1)=0.;mat(0,2)=0.; |
| 95 |
mat(1,0)=0.;mat(1,1)=0.;mat(1,2)=0.; |
| 96 |
mat(2,0)=0.;mat(2,1)=0.;mat(2,2)=0.; |
| 97 |
vec(0)=xyz[0];vec(1)=xyz[1];vec(2)=xyz[2]; |
| 98 |
|
| 99 |
for (int i=0;i<get_nb_fem_noeud();i++) |
| 100 |
{ |
| 101 |
mat(0,0)=mat(0,0)+get_fonction_derive_interpolation(i+1,1,uvw)*get_fem_noeud(i)->get_x(); |
| 102 |
mat(0,1)=mat(0,1)+get_fonction_derive_interpolation(i+1,2,uvw)*get_fem_noeud(i)->get_x(); |
| 103 |
mat(0,2)=mat(0,2)+get_fonction_derive_interpolation(i+1,3,uvw)*get_fem_noeud(i)->get_x(); |
| 104 |
mat(1,0)=mat(1,0)+get_fonction_derive_interpolation(i+1,1,uvw)*get_fem_noeud(i)->get_y(); |
| 105 |
mat(1,1)=mat(1,1)+get_fonction_derive_interpolation(i+1,2,uvw)*get_fem_noeud(i)->get_y(); |
| 106 |
mat(1,2)=mat(1,2)+get_fonction_derive_interpolation(i+1,3,uvw)*get_fem_noeud(i)->get_y(); |
| 107 |
mat(2,0)=mat(2,0)+get_fonction_derive_interpolation(i+1,1,uvw)*get_fem_noeud(i)->get_z(); |
| 108 |
mat(2,1)=mat(2,1)+get_fonction_derive_interpolation(i+1,2,uvw)*get_fem_noeud(i)->get_z(); |
| 109 |
mat(2,2)=mat(2,2)+get_fonction_derive_interpolation(i+1,3,uvw)*get_fem_noeud(i)->get_z(); |
| 110 |
vec(0)=vec(0)-get_fonction_interpolation(i+1,uvw)*get_fem_noeud(i)->get_x(); |
| 111 |
vec(1)=vec(1)-get_fonction_interpolation(i+1,uvw)*get_fem_noeud(i)->get_y(); |
| 112 |
vec(2)=vec(2)-get_fonction_interpolation(i+1,uvw)*get_fem_noeud(i)->get_z(); |
| 113 |
} |
| 114 |
double det=mat.get_determinant(); |
| 115 |
if (fabs(det)<1e-12) ok=3; |
| 116 |
OT_MATRICE_3D mat1(vec,mat.get_vecteur2(),mat.get_vecteur3()); |
| 117 |
OT_MATRICE_3D mat2(mat.get_vecteur1(),vec,mat.get_vecteur3()); |
| 118 |
OT_MATRICE_3D mat3(mat.get_vecteur1(),mat.get_vecteur2(),vec); |
| 119 |
double d1=mat1.get_determinant()/det; |
| 120 |
double d2=mat2.get_determinant()/det; |
| 121 |
double d3=mat3.get_determinant()/det; |
| 122 |
compteur++; |
| 123 |
if (compteur%100==0) ok=2; |
| 124 |
if ((fabs(d1)<1e-8)&&(fabs(d2)<1e-8)&&(fabs(d3)<1e-8)) ok=1; |
| 125 |
uvw[0]=uvw[0]+d1; |
| 126 |
uvw[1]=uvw[1]+d2; |
| 127 |
uvw[2]=uvw[2]+d3; |
| 128 |
} |
| 129 |
} |
| 130 |
return true; |
| 131 |
} |
| 132 |
|
| 133 |
|
| 134 |
void FEM_ELEMENT3::get_interpolation_xyz(double* uvw, double* xyz) |
| 135 |
{ |
| 136 |
xyz[0]=0.;xyz[1]=0.;xyz[2]=0; |
| 137 |
for (int i=0;i<get_nb_fem_noeud();i++) |
| 138 |
{ |
| 139 |
xyz[0]=xyz[0]+get_fonction_interpolation(i+1,uvw)*get_fem_noeud(i)->get_x(); |
| 140 |
xyz[1]=xyz[1]+get_fonction_interpolation(i+1,uvw)*get_fem_noeud(i)->get_y(); |
| 141 |
xyz[2]=xyz[2]+get_fonction_interpolation(i+1,uvw)*get_fem_noeud(i)->get_z(); |
| 142 |
} |
| 143 |
} |
| 144 |
|
| 145 |
double FEM_ELEMENT3::get_jacobien(double* jac,double *uv,double unite) |
| 146 |
{ |
| 147 |
int nb=get_nb_fem_noeud(); |
| 148 |
|
| 149 |
OT_MATRICE_3D jacobien; |
| 150 |
for (int i=0;i<3;i++) |
| 151 |
{ |
| 152 |
jacobien(i,0)=0.;jacobien(i,1)=0.;jacobien(i,2)=0.; |
| 153 |
for (int k=0;k<nb;k++) |
| 154 |
{ |
| 155 |
double valderiv=get_fonction_derive_interpolation(k+1,i+1,uv); |
| 156 |
double *xyz=get_fem_noeud(k)->get_coord(); |
| 157 |
jacobien(i,0)=jacobien(i,0)+valderiv*xyz[0]*unite; |
| 158 |
jacobien(i,1)=jacobien(i,1)+valderiv*xyz[1]*unite; |
| 159 |
jacobien(i,2)=jacobien(i,2)+valderiv*xyz[2]*unite; |
| 160 |
} |
| 161 |
} |
| 162 |
jac[0]=jacobien(0,0); |
| 163 |
jac[1]=jacobien(0,1); |
| 164 |
jac[2]=jacobien(0,2); |
| 165 |
|
| 166 |
jac[3]=jacobien(1,0); |
| 167 |
jac[4]=jacobien(1,1); |
| 168 |
jac[5]=jacobien(1,2); |
| 169 |
|
| 170 |
jac[6]=jacobien(2,0); |
| 171 |
jac[7]=jacobien(2,1); |
| 172 |
jac[8]=jacobien(2,2); |
| 173 |
|
| 174 |
double det=jacobien.get_determinant(); |
| 175 |
return det; |
| 176 |
} |
| 177 |
|
| 178 |
void FEM_ELEMENT3::get_inverse_jacob(double* j,double *uv,double unite) |
| 179 |
{ |
| 180 |
double jac[9]; |
| 181 |
get_jacobien(jac,uv,unite); |
| 182 |
OT_MATRICE_3D J; |
| 183 |
J(0,0)=jac[0]; |
| 184 |
J(0,1)=jac[1]; |
| 185 |
J(0,2)=jac[2]; |
| 186 |
J(1,0)=jac[3]; |
| 187 |
J(1,1)=jac[4]; |
| 188 |
J(1,2)=jac[5]; |
| 189 |
J(2,0)=jac[6]; |
| 190 |
J(2,1)=jac[7]; |
| 191 |
J(2,2)=jac[8]; |
| 192 |
OT_MATRICE_3D j_i=J.inverse(); |
| 193 |
|
| 194 |
|
| 195 |
j[0]=j_i(0,0); |
| 196 |
j[1]=j_i(0,1); |
| 197 |
j[2]=j_i(0,2); |
| 198 |
|
| 199 |
j[3]=j_i(1,0); |
| 200 |
j[4]=j_i(1,1); |
| 201 |
j[5]=j_i(1,2); |
| 202 |
|
| 203 |
j[6]=j_i(2,0); |
| 204 |
j[7]=j_i(2,1); |
| 205 |
j[8]=j_i(2,2); |
| 206 |
} |
