| 1 |
francois |
104 |
#include "toxfem.h" |
| 2 |
|
|
#include "MagXchange.h" |
| 3 |
|
|
#include "mg_file.h" |
| 4 |
|
|
#include "ot_mathematique.h" |
| 5 |
|
|
#include "tpl_map_entite.h" |
| 6 |
|
|
#include "tpl_octree.h" |
| 7 |
|
|
#include <math.h> |
| 8 |
|
|
|
| 9 |
|
|
|
| 10 |
|
|
|
| 11 |
|
|
toxfem::toxfem() |
| 12 |
|
|
{ |
| 13 |
|
|
} |
| 14 |
|
|
|
| 15 |
|
|
|
| 16 |
|
|
toxfem::~toxfem() |
| 17 |
|
|
{ |
| 18 |
|
|
} |
| 19 |
|
|
|
| 20 |
|
|
|
| 21 |
|
|
|
| 22 |
|
|
int toxfem::estdansletetra(MG_TETRA *tet,double x,double y, double z) |
| 23 |
|
|
{ |
| 24 |
|
|
double* xyz1=tet->get_noeud1()->get_coord(); |
| 25 |
|
|
double* xyz2=tet->get_noeud2()->get_coord(); |
| 26 |
|
|
double* xyz3=tet->get_noeud3()->get_coord(); |
| 27 |
|
|
double* xyz4=tet->get_noeud4()->get_coord(); |
| 28 |
|
|
OT_VECTEUR_3D v1(xyz2[0]-xyz1[0],xyz2[1]-xyz1[1],xyz2[2]-xyz1[2]); |
| 29 |
|
|
OT_VECTEUR_3D v2(xyz3[0]-xyz1[0],xyz3[1]-xyz1[1],xyz3[2]-xyz1[2]); |
| 30 |
|
|
OT_VECTEUR_3D v3(xyz4[0]-xyz1[0],xyz4[1]-xyz1[1],xyz4[2]-xyz1[2]); |
| 31 |
|
|
OT_VECTEUR_3D v4(x-xyz1[0],y-xyz1[1],z-xyz1[2]); |
| 32 |
|
|
OT_MATRICE_3D mat(v1,v2,v3); |
| 33 |
|
|
OT_MATRICE_3D mat1(v4,v2,v3); |
| 34 |
|
|
OT_MATRICE_3D mat2(v1,v4,v3); |
| 35 |
|
|
OT_MATRICE_3D mat3(v1,v2,v4); |
| 36 |
|
|
double det=mat.get_determinant(); |
| 37 |
|
|
double xsi=mat1.get_determinant()/det; |
| 38 |
|
|
double eta=mat2.get_determinant()/det; |
| 39 |
|
|
double dseta=mat3.get_determinant()/det; |
| 40 |
|
|
int reponse=1; |
| 41 |
|
|
if (xsi<-0.000001) reponse=0; |
| 42 |
|
|
if (eta<-0.000001) reponse=0; |
| 43 |
|
|
if (dseta<-0.000001) reponse=0; |
| 44 |
|
|
if (xsi+eta+dseta>1.000001) reponse=0; |
| 45 |
|
|
return reponse; |
| 46 |
|
|
|
| 47 |
|
|
} |
| 48 |
|
|
|
| 49 |
|
|
double toxfem::calculdist(double *n,double x,double y,double z,MG_NOEUD* noeud,MG_TETRA* tet) |
| 50 |
|
|
{ |
| 51 |
|
|
double* xyz=noeud->get_coord(); |
| 52 |
|
|
/*double t=(-n[0]*xyz[0]-n[1]*xyz[1]-n[2]*xyz[2]-n[3])/(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]); |
| 53 |
|
|
double xyz1[3]; |
| 54 |
|
|
xyz1[0]=xyz[0]+t*n[0]; |
| 55 |
|
|
xyz1[1]=xyz[1]+t*n[1]; |
| 56 |
|
|
xyz1[2]=xyz[2]+t*n[2]; |
| 57 |
|
|
double dist; |
| 58 |
|
|
if (estdansletetra(tet,xyz1[0],xyz1[1],xyz1[2])) |
| 59 |
|
|
dist=sqrt((xyz[0]-xyz1[0])*(xyz[0]-xyz1[0])+(xyz[1]-xyz1[1])*(xyz[1]-xyz1[1])+(xyz[2]-xyz1[2])*(xyz[2]-xyz1[2])); |
| 60 |
|
|
else*/ |
| 61 |
|
|
double dist=sqrt((xyz[0]-x)*(xyz[0]-x)+(xyz[1]-y)*(xyz[1]-y)+(xyz[2]-z)*(xyz[2]-z)); |
| 62 |
|
|
OT_VECTEUR_3D vec1(n[0],n[1],n[2]); |
| 63 |
|
|
OT_VECTEUR_3D vec2(xyz[0]-x,xyz[1]-y,xyz[2]-z); |
| 64 |
|
|
double ps=vec1*vec2; |
| 65 |
|
|
if (ps<0.) dist=-dist; |
| 66 |
|
|
return dist; |
| 67 |
|
|
} |
| 68 |
|
|
|
| 69 |
|
|
void toxfem::importer(std::string nomfichier,class MagXchange* data) |
| 70 |
|
|
{ |
| 71 |
|
|
MG_FILE gest((char*)nomfichier.c_str()); |
| 72 |
|
|
MG_MAILLAGE* mai=gest.get_mg_maillage(1); |
| 73 |
|
|
MG_GEOMETRIE* geo=gest.get_mg_geometrie(0); |
| 74 |
|
|
int nb_noeud=mai->get_nb_mg_noeud(); |
| 75 |
|
|
std::string nomfichier2=nomfichier+".sol"; |
| 76 |
|
|
MG_SOLUTION* solution=new MG_SOLUTION(mai,1,(char*)nomfichier.c_str(),nb_noeud,"level set"); |
| 77 |
|
|
solution->change_legende(0,"level_0"); |
| 78 |
|
|
gest.ajouter_mg_solution(solution); |
| 79 |
|
|
double xmin=1e300,ymin=1e300,zmin=1e308; |
| 80 |
|
|
double xmax=-1e300,ymax=-1e300,zmax=-1e308; |
| 81 |
|
|
TPL_MAP_ENTITE<MG_NOEUD*> lst_noeud; |
| 82 |
|
|
for (int i=0;i<nb_noeud;i++) |
| 83 |
|
|
{ |
| 84 |
|
|
MG_NOEUD* noeud=mai->get_mg_noeud(i); |
| 85 |
|
|
double* xyz=noeud->get_coord(); |
| 86 |
|
|
if (xyz[0]<xmin) xmin=xyz[0]; |
| 87 |
|
|
if (xyz[1]<ymin) ymin=xyz[1]; |
| 88 |
|
|
if (xyz[2]<zmin) zmin=xyz[2]; |
| 89 |
|
|
if (xyz[0]>xmax) xmax=xyz[0]; |
| 90 |
|
|
if (xyz[1]>ymax) ymax=xyz[1]; |
| 91 |
|
|
if (xyz[2]>zmax) zmax=xyz[2]; |
| 92 |
|
|
noeud->change_solution(1e300); |
| 93 |
|
|
lst_noeud.ajouter(noeud); |
| 94 |
|
|
} |
| 95 |
|
|
TPL_OCTREE<MG_TETRA*,MG_NOEUD*> octree; |
| 96 |
|
|
octree.initialiser(&lst_noeud,1,xmin,ymin,zmin,xmax,ymax,zmax); |
| 97 |
|
|
int nb_tetra=mai->get_nb_mg_tetra(); |
| 98 |
|
|
for (int i=0;i<nb_tetra;i++) |
| 99 |
|
|
{ |
| 100 |
|
|
MG_TETRA* tet=mai->get_mg_tetra(i); |
| 101 |
|
|
octree.inserer(tet); |
| 102 |
|
|
} |
| 103 |
|
|
int nbface=geo->get_nb_mg_face(); |
| 104 |
|
|
for (int iface=0;iface<nbface;iface++) |
| 105 |
|
|
{ |
| 106 |
|
|
MG_FACE* face=geo->get_mg_face(iface);TPL_MAP_ENTITE<MG_NOEUD*> lstnd; |
| 107 |
|
|
TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it; |
| 108 |
|
|
for (MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get_premier(it);ele!=NULL;ele=face->get_lien_maillage()->get_suivant(it)) |
| 109 |
|
|
{ |
| 110 |
|
|
MG_TRIANGLE *tri=(MG_TRIANGLE*)ele; |
| 111 |
|
|
lstnd.ajouter(tri->get_noeud1()); |
| 112 |
|
|
lstnd.ajouter(tri->get_noeud2()); |
| 113 |
|
|
lstnd.ajouter(tri->get_noeud3()); |
| 114 |
|
|
} |
| 115 |
|
|
TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR it2; |
| 116 |
|
|
TPL_MAP_ENTITE<MG_TETRA*> listcoupe; |
| 117 |
|
|
for (MG_NOEUD* noeud=lstnd.get_premier(it2);noeud!=NULL;noeud=lstnd.get_suivant(it2)) |
| 118 |
|
|
{ |
| 119 |
|
|
double uv[2]; |
| 120 |
|
|
double *xyz=noeud->get_coord(); |
| 121 |
|
|
double normal[4]; |
| 122 |
|
|
face->inverser(uv,xyz); |
| 123 |
|
|
face->calcul_normale_unitaire(uv,normal); |
| 124 |
|
|
TPL_MAP_ENTITE<MG_TETRA*> liste; |
| 125 |
|
|
octree.rechercher(xyz[0],xyz[1],xyz[2],0.,liste); |
| 126 |
|
|
int nb=liste.get_nb(); |
| 127 |
|
|
for (int k=0;k<nb;k++) |
| 128 |
|
|
{ |
| 129 |
|
|
MG_TETRA* tet=liste.get(k); |
| 130 |
|
|
if ( estdansletetra(tet,xyz[0],xyz[1],xyz[2])) |
| 131 |
|
|
{ |
| 132 |
|
|
listcoupe.ajouter(tet); |
| 133 |
|
|
normal[3]=-normal[0]*xyz[0]-normal[1]*xyz[1]-normal[2]*xyz[2]; |
| 134 |
|
|
double dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_noeud1(),tet); |
| 135 |
|
|
if (fabs(dist)<fabs(tet->get_noeud1()->get_solution())) |
| 136 |
|
|
tet->get_noeud1()->change_solution(dist); |
| 137 |
|
|
dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_noeud2(),tet); |
| 138 |
|
|
if (fabs(dist)<fabs(tet->get_noeud2()->get_solution())) |
| 139 |
|
|
tet->get_noeud2()->change_solution(dist); |
| 140 |
|
|
dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_noeud3(),tet); |
| 141 |
|
|
if (fabs(dist)<fabs(tet->get_noeud3()->get_solution())) |
| 142 |
|
|
tet->get_noeud3()->change_solution(dist); |
| 143 |
|
|
dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_noeud4(),tet); |
| 144 |
|
|
if (fabs(dist)<fabs(tet->get_noeud4()->get_solution())) |
| 145 |
|
|
tet->get_noeud4()->change_solution(dist); |
| 146 |
|
|
|
| 147 |
|
|
} |
| 148 |
|
|
} |
| 149 |
|
|
} |
| 150 |
|
|
} |
| 151 |
|
|
for (int i=0;i<nb_noeud;i++) |
| 152 |
|
|
{ |
| 153 |
|
|
MG_NOEUD* noeud=mai->get_mg_noeud(i); |
| 154 |
|
|
if (i==866) noeud->change_solution(noeud->get_solution()*(-1)); |
| 155 |
|
|
if (i==921) noeud->change_solution(noeud->get_solution()*(-1)); |
| 156 |
|
|
if (i==1213) noeud->change_solution(noeud->get_solution()*(-1)); |
| 157 |
|
|
if (i==1158) noeud->change_solution(noeud->get_solution()*(-1)); |
| 158 |
|
|
solution->ecrire(i,0,noeud->get_solution()); |
| 159 |
|
|
} |
| 160 |
|
|
|
| 161 |
|
|
for (int i=0;i<nb_noeud;i++) |
| 162 |
|
|
{ |
| 163 |
|
|
MG_NOEUD* noeud=mai->get_mg_noeud(i); |
| 164 |
|
|
MagNode *xfemnode; |
| 165 |
|
|
if (noeud->get_solution()<1e200) |
| 166 |
|
|
{ |
| 167 |
|
|
xfemnode=new MagNode(1); |
| 168 |
|
|
xfemnode->values()[0]=noeud->get_solution(); |
| 169 |
|
|
xfemnode->idvalues()[0]=0; |
| 170 |
|
|
} // 1 valeur de levelset |
| 171 |
|
|
else |
| 172 |
|
|
xfemnode=new MagNode(0); // 0 valeur de levelset |
| 173 |
|
|
double * posxfem=xfemnode->position(); |
| 174 |
|
|
for (int i=0;i<3;++i) posxfem[i]=noeud->get_coord()[i]; |
| 175 |
|
|
xfemnode->id(noeud->get_id()); |
| 176 |
|
|
data->addnode(*xfemnode,true); |
| 177 |
|
|
delete xfemnode; |
| 178 |
|
|
} |
| 179 |
|
|
int nbtetra=mai->get_nb_mg_tetra(); |
| 180 |
|
|
for (int i=0;i<nbtetra;i++) |
| 181 |
|
|
{ |
| 182 |
|
|
MG_TETRA* tet=mai->get_mg_tetra(i); |
| 183 |
|
|
MagElement xfemele(MagElement::TETRAHEDRON); |
| 184 |
|
|
xfemele.nodes()[0]=tet->get_noeud1()->get_id(); |
| 185 |
|
|
xfemele.nodes()[1]=tet->get_noeud2()->get_id(); |
| 186 |
|
|
xfemele.nodes()[2]=tet->get_noeud3()->get_id(); |
| 187 |
|
|
xfemele.nodes()[3]=tet->get_noeud4()->get_id(); |
| 188 |
|
|
data->addelement(xfemele); |
| 189 |
|
|
} |
| 190 |
|
|
} |
