| 1 |
#include "gestionversion.h" |
| 2 |
#include "toibrep.h" |
| 3 |
#include "mg_file.h" |
| 4 |
#include "ot_mathematique.h" |
| 5 |
#include "tpl_map_entite.h" |
| 6 |
#include "tpl_relation_entite.h" |
| 7 |
#include "tpl_octree.h" |
| 8 |
#include "toibrep_point.h" |
| 9 |
#include "vct.h" |
| 10 |
#include <math.h> |
| 11 |
#include <fstream> |
| 12 |
#include <string> |
| 13 |
#include "IBrep.h" |
| 14 |
#include "ot_cpu.h" |
| 15 |
|
| 16 |
|
| 17 |
|
| 18 |
TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm,int nbpas,OT_CPU* compt):geo(ge),gest(g),mai(femm),NPAS(nbpas),affichageactif(0),compteur(compt) |
| 19 |
{ |
| 20 |
} |
| 21 |
|
| 22 |
TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm,int nbpas):geo(ge),gest(g),mai(femm),NPAS(nbpas),affichageactif(0),compteur(NULL) |
| 23 |
{ |
| 24 |
} |
| 25 |
TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm):geo(ge),gest(g),mai(femm),NPAS(50),affichageactif(0),compteur(NULL) |
| 26 |
{ |
| 27 |
} |
| 28 |
TOIBREP::TOIBREP(class MG_GESTIONNAIRE *g,class MG_GEOMETRIE *ge,class FEM_MAILLAGE* femm,class OT_CPU* compt):geo(ge),gest(g),mai(femm),NPAS(50),affichageactif(0),compteur(compt) |
| 29 |
{ |
| 30 |
} |
| 31 |
|
| 32 |
TOIBREP::~TOIBREP() |
| 33 |
{ |
| 34 |
} |
| 35 |
|
| 36 |
|
| 37 |
|
| 38 |
void TOIBREP::active_affichage(void (*fonc)(char*)) |
| 39 |
{ |
| 40 |
affiche=fonc; |
| 41 |
affichageactif=1; |
| 42 |
} |
| 43 |
|
| 44 |
int TOIBREP::compare_etat(int etat,int valeur) |
| 45 |
{ |
| 46 |
if ((etat==1) && (valeur==INTERIEUR)) return 1; |
| 47 |
if ((etat==3) && (valeur==INTERIEUR)) return 1; |
| 48 |
if ((etat==3) && (valeur==STRICTINTERIEUR)) return 1; |
| 49 |
if ((etat==1) && (valeur==SUR_FACE)) return 1; |
| 50 |
return 0; |
| 51 |
} |
| 52 |
|
| 53 |
|
| 54 |
|
| 55 |
int TOIBREP::estdansletetra(FEM_TETRA *tet,double x,double y, double z) |
| 56 |
{ |
| 57 |
double *xyz1,*xyz2,*xyz3,*xyz4; |
| 58 |
if (tet->get_nb_fem_noeud()==4) |
| 59 |
{ |
| 60 |
xyz1=tet->get_fem_noeud(0)->get_coord(); |
| 61 |
xyz2=tet->get_fem_noeud(1)->get_coord(); |
| 62 |
xyz3=tet->get_fem_noeud(2)->get_coord(); |
| 63 |
xyz4=tet->get_fem_noeud(3)->get_coord(); |
| 64 |
} |
| 65 |
if (tet->get_nb_fem_noeud()==10) |
| 66 |
{ |
| 67 |
xyz1=tet->get_fem_noeud(0)->get_coord(); |
| 68 |
xyz2=tet->get_fem_noeud(2)->get_coord(); |
| 69 |
xyz3=tet->get_fem_noeud(4)->get_coord(); |
| 70 |
xyz4=tet->get_fem_noeud(9)->get_coord(); |
| 71 |
} |
| 72 |
OT_VECTEUR_3D v1(xyz2[0]-xyz1[0],xyz2[1]-xyz1[1],xyz2[2]-xyz1[2]); |
| 73 |
OT_VECTEUR_3D v2(xyz3[0]-xyz1[0],xyz3[1]-xyz1[1],xyz3[2]-xyz1[2]); |
| 74 |
OT_VECTEUR_3D v3(xyz4[0]-xyz1[0],xyz4[1]-xyz1[1],xyz4[2]-xyz1[2]); |
| 75 |
OT_VECTEUR_3D v4(x-xyz1[0],y-xyz1[1],z-xyz1[2]); |
| 76 |
OT_MATRICE_3D mat(v1,v2,v3); |
| 77 |
OT_MATRICE_3D mat1(v4,v2,v3); |
| 78 |
OT_MATRICE_3D mat2(v1,v4,v3); |
| 79 |
OT_MATRICE_3D mat3(v1,v2,v4); |
| 80 |
double det=mat.get_determinant(); |
| 81 |
double xsi=mat1.get_determinant()/det; |
| 82 |
double eta=mat2.get_determinant()/det; |
| 83 |
double dseta=mat3.get_determinant()/det; |
| 84 |
int reponse1=1; |
| 85 |
int reponse2=1; |
| 86 |
double eps=0.0000001; |
| 87 |
if (xsi<-eps) reponse1=0; |
| 88 |
if (eta<-eps) reponse1=0; |
| 89 |
if (dseta<-eps) reponse1=0; |
| 90 |
if (xsi+eta+dseta>1.+eps) reponse1=0; |
| 91 |
if (xsi<eps) reponse2=0; |
| 92 |
if (eta<eps) reponse2=0; |
| 93 |
if (dseta<eps) reponse2=0; |
| 94 |
if (xsi+eta+dseta>1.-eps) reponse2=0; |
| 95 |
return reponse1+2*reponse2; |
| 96 |
} |
| 97 |
|
| 98 |
|
| 99 |
|
| 100 |
double TOIBREP::calculdist(double *n,double x,double y,double z,FEM_NOEUD* noeud) |
| 101 |
{ |
| 102 |
double* xyz=noeud->get_coord(); |
| 103 |
double dist=sqrt((xyz[0]-x)*(xyz[0]-x)+(xyz[1]-y)*(xyz[1]-y)+(xyz[2]-z)*(xyz[2]-z)); |
| 104 |
OT_VECTEUR_3D vec1(n[0],n[1],n[2]); |
| 105 |
OT_VECTEUR_3D vec2(xyz[0]-x,xyz[1]-y,xyz[2]-z); |
| 106 |
double ps=vec1*vec2; |
| 107 |
if (ps<0.) dist=-dist; |
| 108 |
return dist; |
| 109 |
} |
| 110 |
|
| 111 |
void TOIBREP::recherche_arete_tangeante(TPL_MAP_ENTITE<MG_ARETE*> &lst,TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> &lsttopo) |
| 112 |
{ |
| 113 |
|
| 114 |
|
| 115 |
LISTE_MG_ARETE::iterator it; |
| 116 |
for (MG_ARETE* are=geo->get_premier_arete(it);are!=NULL;are=geo->get_suivant_arete(it)) |
| 117 |
{ |
| 118 |
if (lsttopo.existe(are)) |
| 119 |
if (ot.arete_virtuelle(are)) lst.ajouter(are); |
| 120 |
} |
| 121 |
|
| 122 |
} |
| 123 |
|
| 124 |
|
| 125 |
//void TOIBREP::importer(std::string nomfichier,class MagXchange* data,std::string nomfichierout) |
| 126 |
IBrep TOIBREP::importer(std::string nomfichier,std::string nomfichieribrep,MG_GROUPE_TOPOLOGIQUE* mggt) |
| 127 |
{ |
| 128 |
TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> lst; |
| 129 |
if (mggt!=NULL) |
| 130 |
{ |
| 131 |
int nb=mggt->get_nb(); |
| 132 |
for (int i=0;i<nb;i++) |
| 133 |
{ |
| 134 |
lst.ajouter(mggt->get(i)); |
| 135 |
mggt->get(i)->get_topologie_sousjacente(&lst); |
| 136 |
} |
| 137 |
} |
| 138 |
else |
| 139 |
{ |
| 140 |
int nbvol=geo->get_nb_mg_volume(); |
| 141 |
for (int i=0;i<nbvol;i++) |
| 142 |
{ |
| 143 |
lst.ajouter(geo->get_mg_volume(i)); |
| 144 |
geo->get_mg_volume(i)->get_topologie_sousjacente(&lst); |
| 145 |
} |
| 146 |
} |
| 147 |
affiche((char*)"Recherche aretes tangeantes"); |
| 148 |
TPL_MAP_ENTITE<MG_ARETE*> lstaretetan; |
| 149 |
recherche_arete_tangeante(lstaretetan,lst); |
| 150 |
char mess[300]; |
| 151 |
sprintf(mess," Nombre d'aretes tangeantes %d",lstaretetan.get_nb()); |
| 152 |
affiche(mess); |
| 153 |
if (compteur!=NULL) compteur->ajouter_etape("Aretes tangeantes"); |
| 154 |
affiche((char*)"Traitement de base"); |
| 155 |
int nbface=geo->get_nb_mg_face(); |
| 156 |
int nbvraiface=0; |
| 157 |
TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*>::ITERATEUR it; |
| 158 |
for (MG_ELEMENT_TOPOLOGIQUE* ele=lst.get_premier(it);ele!=NULL;ele=lst.get_suivant(it)) |
| 159 |
if (ele->get_dimension()==2) nbvraiface++; |
| 160 |
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base selection"); |
| 161 |
int nb_noeud=mai->get_nb_fem_noeud(); |
| 162 |
int nb_tetra=mai->get_nb_fem_tetra(); |
| 163 |
std::string nomfichier2=nomfichier+"_dis.sol"; |
| 164 |
FEM_SOLUTION* solution=new FEM_SOLUTION(mai,nbvraiface+lstaretetan.get_nb(),(char*)nomfichier2.c_str(),nb_noeud,"LS",ENTITE_NOEUD); |
| 165 |
int i=0; |
| 166 |
gest->ajouter_fem_solution(solution); |
| 167 |
nomfichier2=nomfichier+"_ele.sol"; |
| 168 |
FEM_SOLUTION* solution_ele=new FEM_SOLUTION(mai,nbvraiface+lstaretetan.get_nb(),(char*)nomfichier2.c_str(),nb_tetra,"A",ENTITE_TETRA); |
| 169 |
gest->ajouter_fem_solution(solution_ele); |
| 170 |
for (MG_ELEMENT_TOPOLOGIQUE* ele=lst.get_premier(it);ele!=NULL;ele=lst.get_suivant(it)) |
| 171 |
{ |
| 172 |
if (ele->get_dimension()!=2) continue; |
| 173 |
char mess[255]; |
| 174 |
sprintf(mess,"F %lu",ele->get_id()); |
| 175 |
solution->change_legende(i,mess); |
| 176 |
solution_ele->change_legende(i,mess); |
| 177 |
++i; |
| 178 |
} |
| 179 |
TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR itare; |
| 180 |
for (MG_ARETE* are=lstaretetan.get_premier(itare);are!=NULL;are=lstaretetan.get_suivant(itare)) |
| 181 |
{ |
| 182 |
char mess[255]; |
| 183 |
sprintf(mess,"A %lu",are->get_id()); |
| 184 |
solution->change_legende(i,mess); |
| 185 |
solution_ele->change_legende(i,mess); |
| 186 |
++i; |
| 187 |
} |
| 188 |
|
| 189 |
|
| 190 |
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base creation solution"); |
| 191 |
double xmin=1e300,ymin=1e300,zmin=1e308; |
| 192 |
double xmax=-1e300,ymax=-1e300,zmax=-1e308; |
| 193 |
TPL_MAP_ENTITE<FEM_NOEUD*> lst_noeud; |
| 194 |
LISTE_FEM_NOEUD::iterator it2; |
| 195 |
i=0; |
| 196 |
for (FEM_NOEUD* noeud=mai->get_premier_noeud(it2);noeud!=NULL;noeud=mai->get_suivant_noeud(it2)) |
| 197 |
{ |
| 198 |
double* xyz=noeud->get_coord(); |
| 199 |
if (xyz[0]<xmin) xmin=xyz[0]; |
| 200 |
if (xyz[1]<ymin) ymin=xyz[1]; |
| 201 |
if (xyz[2]<zmin) zmin=xyz[2]; |
| 202 |
if (xyz[0]>xmax) xmax=xyz[0]; |
| 203 |
if (xyz[1]>ymax) ymax=xyz[1]; |
| 204 |
if (xyz[2]>zmax) zmax=xyz[2]; |
| 205 |
for (int j=0;j<nbvraiface+lstaretetan.get_nb();j++) |
| 206 |
solution->ecrire(i,j,1e300); |
| 207 |
lst_noeud.ajouter(noeud); |
| 208 |
i++; |
| 209 |
} |
| 210 |
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base recherche boite"); |
| 211 |
octree.initialiser(&lst_noeud,1,xmin,ymin,zmin,xmax,ymax,zmax); |
| 212 |
//if (compteur!=NULL) compteur->ajouter_etape("Traitement de base ini octree"); |
| 213 |
LISTE_FEM_TETRA::iterator it3; |
| 214 |
i=0; |
| 215 |
for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(it3)) |
| 216 |
{ |
| 217 |
octree.inserer(tet); |
| 218 |
for (int j=0;j<nbvraiface+lstaretetan.get_nb();j++) |
| 219 |
solution_ele->ecrire(i,j,0); |
| 220 |
i++; |
| 221 |
} |
| 222 |
if (compteur!=NULL) compteur->ajouter_etape("Traitement de base"); |
| 223 |
for (int i=0;i<nbface;i++) |
| 224 |
{ |
| 225 |
MG_FACE* face=geo->get_mg_face(i); |
| 226 |
if (affichageactif) |
| 227 |
{ |
| 228 |
char mess[300]; |
| 229 |
sprintf(mess," Face %d numero %lu",i,face->get_id()); |
| 230 |
affiche(mess); |
| 231 |
} |
| 232 |
if (!lst.existe(face)) continue; |
| 233 |
vector<MG_FACE*> lstface; |
| 234 |
lstface.push_back(face); |
| 235 |
/// ici topologie virtuelle |
| 236 |
|
| 237 |
|
| 238 |
/*TPL_LISTE_ENTITE<double> lst; |
| 239 |
int type=face->get_surface()->get_type_geometrique(lst); |
| 240 |
int idem=0; |
| 241 |
for (int j=0;j<i;j++) |
| 242 |
{ |
| 243 |
TPL_LISTE_ENTITE<double> lst2; |
| 244 |
MG_FACE* face2=geo->get_mg_face(j); |
| 245 |
int type2=face2->get_surface()->get_type_geometrique(lst2); |
| 246 |
if (type==type2) |
| 247 |
if (lst.get_nb()==lst2.get_nb()) |
| 248 |
{ |
| 249 |
int diff=0; |
| 250 |
for (int i=0;i<lst.get_nb();i++) |
| 251 |
if (fabs(lst.get(i)-lst2.get(i))>0.000001) diff=1; |
| 252 |
if (diff==0) idem=1; |
| 253 |
} |
| 254 |
} |
| 255 |
if (!idem) lstface.push_back(geo->get_mg_face(i)); |
| 256 |
for (int j=i+1;j<nbface;j++) |
| 257 |
{ |
| 258 |
TPL_LISTE_ENTITE<double> lst2; |
| 259 |
MG_FACE* face2=geo->get_mg_face(j); |
| 260 |
int type2=face2->get_surface()->get_type_geometrique(lst2); |
| 261 |
if (type==type2) |
| 262 |
if (lst.get_nb()==lst2.get_nb()) |
| 263 |
{ |
| 264 |
int diff=0; |
| 265 |
for (int i=0;i<lst.get_nb();i++) |
| 266 |
if (fabs(lst.get(i)-lst2.get(i))>0.000001) diff=1; |
| 267 |
if (diff==0) lstface.push_back(face2); |
| 268 |
} |
| 269 |
}*/ |
| 270 |
levelsetn(&lst,&lstface,solution,solution_ele,i); |
| 271 |
} |
| 272 |
|
| 273 |
TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR itare2; |
| 274 |
i=0; |
| 275 |
for (MG_ARETE* are=lstaretetan.get_premier(itare2);are!=NULL;are=lstaretetan.get_suivant(itare2)) |
| 276 |
{ |
| 277 |
if (affichageactif) |
| 278 |
{ |
| 279 |
char mess[300]; |
| 280 |
sprintf(mess," Arete tangeante %d numero %lu",i,are->get_id()); |
| 281 |
affiche(mess); |
| 282 |
} |
| 283 |
traite_arete_tangeante(are,nbface+i,solution,solution_ele); |
| 284 |
i++; |
| 285 |
} |
| 286 |
|
| 287 |
|
| 288 |
affiche((char*)"Exportation IBREP"); |
| 289 |
IBrep brep=exporter_IBrep(nomfichieribrep,solution,solution_ele,mggt); |
| 290 |
if (compteur!=NULL) compteur->ajouter_etape("Exportation IBREP"); |
| 291 |
return brep; |
| 292 |
} |
| 293 |
|
| 294 |
void TOIBREP::traite_arete_tangeante(MG_ARETE* are,int numsol,FEM_SOLUTION* solution,FEM_SOLUTION* solution_ele) |
| 295 |
{ |
| 296 |
//NPAS=150; |
| 297 |
solution->active_solution(numsol); |
| 298 |
solution_ele->active_solution(numsol); |
| 299 |
vector<TOIBREP_POINT*> lst; |
| 300 |
LISTE_FEM_TETRA::iterator ittetra; |
| 301 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra)) |
| 302 |
tet->change_etat(0); |
| 303 |
LISTE_FEM_NOEUD::iterator itnoeud; |
| 304 |
for (FEM_NOEUD* no=mai->get_premier_noeud(itnoeud);no!=NULL;no=mai->get_suivant_noeud(itnoeud)) |
| 305 |
no->change_etat(0); |
| 306 |
double tmin=are->get_tmin(); |
| 307 |
double tmax=are->get_tmax(); |
| 308 |
double eps=1e-3*(tmax-tmin); |
| 309 |
for (int i=0;i<NPAS+1;i++) |
| 310 |
{ |
| 311 |
double t=tmin+i*1.0/NPAS*(tmax-tmin); |
| 312 |
double xyz[3],xyzplus[3],xyzmoins[3]; |
| 313 |
are->evaluer(t,xyz); |
| 314 |
TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,NULL); |
| 315 |
lst.push_back(pt); |
| 316 |
if (t+eps<tmax) |
| 317 |
{ |
| 318 |
are->evaluer(t+eps,xyzplus); |
| 319 |
pt->ajoute_point_eps(xyzplus); |
| 320 |
} |
| 321 |
if (t-eps>tmin) |
| 322 |
{ |
| 323 |
are->evaluer(t-eps,xyzmoins); |
| 324 |
pt->ajoute_point_eps(xyzmoins); |
| 325 |
} |
| 326 |
|
| 327 |
} |
| 328 |
for (int i=0;i<lst.size();i++) |
| 329 |
{ |
| 330 |
TOIBREP_POINT* pt=lst[i]; |
| 331 |
double xyz[3]; |
| 332 |
pt->get_coord3(xyz); |
| 333 |
TPL_MAP_ENTITE<FEM_TETRA*> liste; |
| 334 |
octree.rechercher(xyz[0],xyz[1],xyz[2],0.,liste); |
| 335 |
TPL_MAP_ENTITE<FEM_TETRA*>::ITERATEUR it; |
| 336 |
for (FEM_TETRA* tet=liste.get_premier(it);tet!=NULL;tet=liste.get_suivant(it)) |
| 337 |
{ |
| 338 |
if (tet->get_etat()==1) continue; |
| 339 |
int etat=estdansletetra(tet,xyz[0],xyz[1],xyz[2]); |
| 340 |
if ((compare_etat(etat,SUR_FACE)) && (pt->get_interieur())) |
| 341 |
{ |
| 342 |
int nbeps=pt->get_nb_point_eps(); |
| 343 |
for (int k=0;k<nbeps;k++) |
| 344 |
{ |
| 345 |
double xyztmp[3]; |
| 346 |
pt->get_point_eps(k,xyztmp); |
| 347 |
int etattmp=estdansletetra(tet,xyztmp[0],xyztmp[1],xyztmp[2]); |
| 348 |
if (compare_etat(etattmp,SUR_FACE)) {etat=STRICTINTERIEUR;break;} |
| 349 |
if (compare_etat(etattmp,INTERIEUR)) {etat=STRICTINTERIEUR;break;} |
| 350 |
if (compare_etat(etattmp,STRICTINTERIEUR)) {etat=STRICTINTERIEUR;break;} |
| 351 |
} |
| 352 |
} |
| 353 |
if (compare_etat(etat,STRICTINTERIEUR)) |
| 354 |
{ |
| 355 |
for (int k=0;k<tet->get_nb_fem_noeud();k++) |
| 356 |
{ |
| 357 |
double dist=calcul_distance_level_ortho(tet->get_fem_noeud(k),are,pt); |
| 358 |
tet->get_fem_noeud(k)->change_solution(dist); |
| 359 |
tet->get_fem_noeud(k)->change_numero(pt->get_id()); |
| 360 |
tet->change_etat(1); //temp |
| 361 |
tet->get_fem_noeud(k)->change_etat(1); //temp |
| 362 |
} |
| 363 |
} |
| 364 |
} |
| 365 |
} |
| 366 |
remplir_trou_tangeant(&lst,are,solution_ele); |
| 367 |
|
| 368 |
|
| 369 |
|
| 370 |
int i=0; |
| 371 |
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud)) |
| 372 |
{ |
| 373 |
if (noeud->get_etat()==1) solution->ecrire(i,numsol,noeud->get_solution()); |
| 374 |
else solution->ecrire(i,numsol,1e300); |
| 375 |
++i; |
| 376 |
} |
| 377 |
i=0; |
| 378 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra)) |
| 379 |
{ |
| 380 |
if (tet->get_etat()==1) solution_ele->ecrire(i,numsol,1.); |
| 381 |
++i; |
| 382 |
} |
| 383 |
int nbpt=lst.size(); |
| 384 |
for (int i=0;i<nbpt;i++) delete lst[i]; |
| 385 |
TOIBREP_POINT::remisecompteurid(); |
| 386 |
} |
| 387 |
|
| 388 |
|
| 389 |
void TOIBREP::echantillonne_sommets(std::vector<TOIBREP_POINT*> &lst,MG_FACE* face) |
| 390 |
{ |
| 391 |
TPL_MAP_ENTITE<MG_SOMMET*> lstsom; |
| 392 |
int nbbou=face->get_nb_mg_boucle(); |
| 393 |
for (int i=0;i<nbbou;i++) |
| 394 |
{ |
| 395 |
MG_BOUCLE* boucle=face->get_mg_boucle(i); |
| 396 |
int nb_arete=boucle->get_nb_mg_coarete(); |
| 397 |
for (int j=0;j<nb_arete;j++) |
| 398 |
{ |
| 399 |
MG_ARETE* are=boucle->get_mg_coarete(j)->get_arete(); |
| 400 |
lstsom.ajouter(are->get_cosommet1()->get_sommet()); |
| 401 |
lstsom.ajouter(are->get_cosommet2()->get_sommet()); |
| 402 |
} |
| 403 |
} |
| 404 |
TPL_MAP_ENTITE<MG_SOMMET*>::ITERATEUR it; |
| 405 |
for (MG_SOMMET* som=lstsom.get_premier(it);som!=NULL;som=lstsom.get_suivant(it)) |
| 406 |
{ |
| 407 |
double xyz[3]; |
| 408 |
som->get_point()->evaluer(xyz); |
| 409 |
double uv[2]; |
| 410 |
face->inverser(uv,xyz); |
| 411 |
int N=NPAS; |
| 412 |
double rayon=1e-2; |
| 413 |
for (int k=0;k<N+1;k++) |
| 414 |
{ |
| 415 |
double uv2[2]; |
| 416 |
uv2[0]=uv[0]+rayon*cos(k*1.0/N*2.*M_PI); |
| 417 |
uv2[1]=uv[1]+rayon*sin(k*1.0/N*2.*M_PI); |
| 418 |
if ((face->valide_parametre_u(uv2[0])) && (face->valide_parametre_v(uv2[1]))) |
| 419 |
{ |
| 420 |
double distance=ot.calcule_distance_contour_face_uv(uv2,face); |
| 421 |
if (distance>0.) |
| 422 |
{ |
| 423 |
double xyz[3]; |
| 424 |
face->evaluer(uv2,xyz); |
| 425 |
TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],1,face); |
| 426 |
lst.push_back(pt); |
| 427 |
} |
| 428 |
} |
| 429 |
} |
| 430 |
|
| 431 |
} |
| 432 |
} |
| 433 |
void TOIBREP::echantillonne_aretes(std::vector<TOIBREP_POINT*> &lst,MG_FACE* face) |
| 434 |
{ |
| 435 |
echantillonne_sommets(lst,face); |
| 436 |
int nbbou=face->get_nb_mg_boucle(); |
| 437 |
for (int i=0;i<nbbou;i++) |
| 438 |
{ |
| 439 |
MG_BOUCLE* boucle=face->get_mg_boucle(i); |
| 440 |
int nb_arete=boucle->get_nb_mg_coarete(); |
| 441 |
for (int j=0;j<nb_arete;j++) |
| 442 |
{ |
| 443 |
MG_ARETE* are=boucle->get_mg_coarete(j)->get_arete(); |
| 444 |
double tmin=are->get_tmin(); |
| 445 |
double tmax=are->get_tmax(); |
| 446 |
double N=NPAS; |
| 447 |
for (int k=0;k<N+1;k++) |
| 448 |
{ |
| 449 |
double t=tmin+k*1.0/N*(tmax-tmin); |
| 450 |
double xyz[3]; |
| 451 |
are->evaluer(t,xyz); |
| 452 |
TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face); |
| 453 |
lst.push_back(pt); |
| 454 |
// octree.inserer(pt); |
| 455 |
} |
| 456 |
/* double t=tmin+0.2/NPAS*(tmax-tmin); |
| 457 |
double xyz[3]; |
| 458 |
are->evaluer(t,xyz); |
| 459 |
TOIBREP_POINT *pt1=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face); |
| 460 |
lst.push_back(pt1); |
| 461 |
t=tmax-0.2/NPAS*(tmax-tmin); |
| 462 |
are->evaluer(t,xyz); |
| 463 |
TOIBREP_POINT *pt2=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,1,face); |
| 464 |
lst.push_back(pt2);*/ |
| 465 |
|
| 466 |
|
| 467 |
|
| 468 |
} |
| 469 |
} |
| 470 |
} |
| 471 |
|
| 472 |
void TOIBREP::levelsetn(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lsttopo,vector<MG_FACE*> *lstface,class FEM_SOLUTION* solution,FEM_SOLUTION* solution_ele,int numsol) |
| 473 |
{ |
| 474 |
if (affichageactif) |
| 475 |
{ |
| 476 |
char mess[300]; |
| 477 |
sprintf(mess," Echantillonnage"); |
| 478 |
affiche(mess); |
| 479 |
} |
| 480 |
|
| 481 |
solution->active_solution(numsol); |
| 482 |
solution_ele->active_solution(numsol); |
| 483 |
int nbface=lstface->size(); |
| 484 |
if (nbface==0) return; |
| 485 |
vector<TOIBREP_POINT*> lst; |
| 486 |
BOITE_2D boite=ot.get_boite_2D((*lstface)[0]); |
| 487 |
echantillonne_aretes(lst,(*lstface)[0]); |
| 488 |
for (int iface=1;iface<nbface;iface++) |
| 489 |
{ |
| 490 |
MG_FACE* face=(*lstface)[iface]; |
| 491 |
BOITE_2D boiteface=ot.get_boite_2D(face); |
| 492 |
boite=boite+boiteface; |
| 493 |
echantillonne_aretes(lst,(*lstface)[iface]); |
| 494 |
} |
| 495 |
MG_FACE* face=(*lstface)[0]; |
| 496 |
int trouve=0; |
| 497 |
int orientation; |
| 498 |
int j=0; |
| 499 |
do |
| 500 |
{ |
| 501 |
MG_COFACE* coface=face->get_mg_coface(j); |
| 502 |
MG_VOLUME* vol=coface->get_coquille()->get_mg_volume(); |
| 503 |
if (lsttopo->existe(vol)) |
| 504 |
{ |
| 505 |
orientation=coface->get_orientation(); |
| 506 |
trouve=1; |
| 507 |
} |
| 508 |
j++; |
| 509 |
} |
| 510 |
while (trouve==0); |
| 511 |
LISTE_FEM_TETRA::iterator ittetra; |
| 512 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra)) |
| 513 |
tet->change_etat(0); |
| 514 |
LISTE_FEM_NOEUD::iterator itnoeud; |
| 515 |
for (FEM_NOEUD* no=mai->get_premier_noeud(itnoeud);no!=NULL;no=mai->get_suivant_noeud(itnoeud)) |
| 516 |
no->change_etat(0); |
| 517 |
double umin=boite.get_xmin(); |
| 518 |
double umax=boite.get_xmax(); |
| 519 |
double vmin=boite.get_ymin(); |
| 520 |
double vmax=boite.get_ymax(); |
| 521 |
if (face->get_surface()->est_periodique_u()) |
| 522 |
{ |
| 523 |
umin=0; |
| 524 |
umax=face->get_surface()->get_periode_u(); |
| 525 |
} |
| 526 |
if (face->get_surface()->est_periodique_v()) |
| 527 |
{ |
| 528 |
vmin=0; |
| 529 |
vmax=face->get_surface()->get_periode_v(); |
| 530 |
} |
| 531 |
double epsu=1e-3*(umax-umin); |
| 532 |
double epsv=1e-3*(vmax-vmin); |
| 533 |
for (int i=0;i<NPAS;i++) |
| 534 |
for (int j=0;j<NPAS;j++) |
| 535 |
{ |
| 536 |
double uv[2]; |
| 537 |
uv[0]=umin+(i*1.0+0.5)/NPAS*(umax-umin); |
| 538 |
uv[1]=vmin+(j*1.0+0.5)/NPAS*(vmax-vmin); |
| 539 |
double xyz[3]; |
| 540 |
if ((face->valide_parametre_u(uv[0])) && (face->valide_parametre_v(uv[1]))) |
| 541 |
{ |
| 542 |
double distance=ot.calcule_distance_contour_face_uv(uv,face); |
| 543 |
int inte=0; |
| 544 |
if (distance>0.) inte=1.; else continue; |
| 545 |
face->evaluer(uv,xyz); |
| 546 |
TOIBREP_POINT *pt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],inte,face); |
| 547 |
lst.push_back(pt); |
| 548 |
if (uv[0]+epsu<umax) |
| 549 |
{ |
| 550 |
double uvtmp[2]={uv[0]+epsu,uv[1]}; |
| 551 |
double distance=ot.calcule_distance_contour_face_uv(uvtmp,face); |
| 552 |
if (distance>1e-8) |
| 553 |
{ |
| 554 |
face->evaluer(uvtmp,xyz); |
| 555 |
pt->ajoute_point_eps(xyz); |
| 556 |
} |
| 557 |
} |
| 558 |
if (uv[0]-epsu>umin) |
| 559 |
{ |
| 560 |
double uvtmp[2]={uv[0]-epsu,uv[1]}; |
| 561 |
double distance=ot.calcule_distance_contour_face_uv(uvtmp,face); |
| 562 |
if (distance>1e-8) |
| 563 |
{ |
| 564 |
face->evaluer(uvtmp,xyz); |
| 565 |
pt->ajoute_point_eps(xyz); |
| 566 |
} |
| 567 |
} |
| 568 |
if (uv[1]+epsv<vmax) |
| 569 |
{ |
| 570 |
double uvtmp[2]={uv[0],uv[1]+epsv}; |
| 571 |
double distance=ot.calcule_distance_contour_face_uv(uvtmp,face); |
| 572 |
if (distance>1e-8) |
| 573 |
{ |
| 574 |
face->evaluer(uvtmp,xyz); |
| 575 |
pt->ajoute_point_eps(xyz); |
| 576 |
} |
| 577 |
} |
| 578 |
if (uv[1]-epsv>vmin) |
| 579 |
{ |
| 580 |
double uvtmp[2]={uv[0],uv[1]-epsv}; |
| 581 |
double distance=ot.calcule_distance_contour_face_uv(uvtmp,face); |
| 582 |
if (distance>1e-8) |
| 583 |
{ |
| 584 |
face->evaluer(uvtmp,xyz); |
| 585 |
pt->ajoute_point_eps(xyz); |
| 586 |
} |
| 587 |
} |
| 588 |
// octree.inserer(pt); |
| 589 |
} |
| 590 |
} |
| 591 |
// debut temp |
| 592 |
/* |
| 593 |
MG_GESTIONNAIRE gestemp; |
| 594 |
MG_MAILLAGE *mai=new MG_MAILLAGE(NULL); |
| 595 |
gestemp.ajouter_mg_maillage(mai); |
| 596 |
for (int i=0;i<lst.size();i++) |
| 597 |
{ |
| 598 |
double xyz[3]; |
| 599 |
lst[i]->get_coord3(xyz); |
| 600 |
MG_NOEUD* no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],TRIANGULATION); |
| 601 |
mai->ajouter_mg_noeud(no); |
| 602 |
if (i!=0) {MG_SEGMENT* seg=new MG_SEGMENT(NULL,no,mai->get_mg_noeud(mai->get_nb_mg_noeud()-2),TRIANGULATION);mai->ajouter_mg_segment(seg);} |
| 603 |
} |
| 604 |
|
| 605 |
gestemp.enregistrer("debug.magic"); |
| 606 |
*/ |
| 607 |
// fin temp |
| 608 |
if (compteur!=NULL) compteur->ajouter_etape("Echantillonnage"); |
| 609 |
if (affichageactif) |
| 610 |
{ |
| 611 |
char mess[300]; |
| 612 |
sprintf(mess," Preparation newton"); |
| 613 |
affiche(mess); |
| 614 |
} |
| 615 |
calcullevelsetpremierepasse(face,orientation,&lst,0,lst.size()); |
| 616 |
if (compteur!=NULL) compteur->ajouter_etape("Preparation newton"); |
| 617 |
|
| 618 |
if (affichageactif) |
| 619 |
{ |
| 620 |
char mess[300]; |
| 621 |
sprintf(mess," Newton"); |
| 622 |
affiche(mess); |
| 623 |
} |
| 624 |
calcullevelsetdeuxiemepasse(&lst); |
| 625 |
if (compteur!=NULL) compteur->ajouter_etape("Newton"); |
| 626 |
if (affichageactif) |
| 627 |
{ |
| 628 |
char mess[300]; |
| 629 |
sprintf(mess," Remplissage des trous"); |
| 630 |
affiche(mess); |
| 631 |
} |
| 632 |
remplir_trou(&lst,face,orientation); |
| 633 |
if (compteur!=NULL) compteur->ajouter_etape("Remplissage des trous"); |
| 634 |
//etendrelevelset(solution,numsol); |
| 635 |
int i=0; |
| 636 |
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud)) |
| 637 |
{ |
| 638 |
if (noeud->get_etat()==1) solution->ecrire(i,numsol,noeud->get_solution()); |
| 639 |
else solution->ecrire(i,numsol,1e300); |
| 640 |
++i; |
| 641 |
} |
| 642 |
i=0; |
| 643 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra)) |
| 644 |
{ |
| 645 |
if (tet->get_etat()==1) solution_ele->ecrire(i,numsol,1.); |
| 646 |
++i; |
| 647 |
} |
| 648 |
int nbpt=lst.size(); |
| 649 |
for (int i=0;i<nbpt;i++) delete lst[i]; |
| 650 |
TOIBREP_POINT::remisecompteurid(); |
| 651 |
|
| 652 |
|
| 653 |
|
| 654 |
} |
| 655 |
void TOIBREP::remplir_trou_tangeant(std::vector<TOIBREP_POINT*> *lst,MG_ARETE* are,FEM_SOLUTION *solution_ele) |
| 656 |
{ |
| 657 |
MG_FACE* face1=are->get_mg_coarete(0)->get_boucle()->get_mg_face(); |
| 658 |
MG_FACE* face2=are->get_mg_coarete(1)->get_boucle()->get_mg_face(); |
| 659 |
TPL_MAP_ENTITE<FEM_TETRA*> lst_tetra_concerne; |
| 660 |
int num1,num2; |
| 661 |
int nbchamp=solution_ele->get_nb_champ(); |
| 662 |
for (int i=0;i<nbchamp;i++) |
| 663 |
{ |
| 664 |
std::string nom=solution_ele->get_legende(i); |
| 665 |
char nom2[2]; |
| 666 |
unsigned long id; |
| 667 |
sscanf(nom.c_str(),"%s %lu",nom2,&id); |
| 668 |
if (id==face1->get_id()) num1=i; |
| 669 |
if (id==face2->get_id()) num2=i; |
| 670 |
} |
| 671 |
LISTE_FEM_TETRA::iterator ittet; |
| 672 |
int i=0; |
| 673 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet)) |
| 674 |
{ |
| 675 |
if (solution_ele->lire(i,num1)>0.000001) |
| 676 |
lst_tetra_concerne.ajouter(tet); |
| 677 |
if (solution_ele->lire(i,num2)>0.000001) |
| 678 |
lst_tetra_concerne.ajouter(tet); |
| 679 |
i++; |
| 680 |
} |
| 681 |
|
| 682 |
|
| 683 |
int couche_active; |
| 684 |
int nb_couche=0; |
| 685 |
do |
| 686 |
{ |
| 687 |
nb_couche++; |
| 688 |
TPL_MAP_ENTITE<FEM_TETRA*>::ITERATEUR ittetra; |
| 689 |
std::vector<FEM_TETRA*> couche; |
| 690 |
couche_active=0; |
| 691 |
for (FEM_TETRA* tet=lst_tetra_concerne.get_premier(ittetra);tet!=NULL;tet=lst_tetra_concerne.get_suivant(ittetra)) |
| 692 |
{ |
| 693 |
if (tet->get_etat()!=0) continue; |
| 694 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 695 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 696 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 697 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 698 |
int nbactif=0; |
| 699 |
int nbpositif=0; |
| 700 |
if ((no1->get_etat()==1) || (no1->get_etat()==-1)) |
| 701 |
{ |
| 702 |
nbactif++; |
| 703 |
if (no1->get_solution()>0.) nbpositif++; |
| 704 |
} |
| 705 |
if ((no2->get_etat()==1) || (no2->get_etat()==-1)) |
| 706 |
{ |
| 707 |
nbactif++; |
| 708 |
if (no2->get_solution()>0.) nbpositif++; |
| 709 |
} |
| 710 |
if ((no3->get_etat()==1) || (no3->get_etat()==-1)) |
| 711 |
{ |
| 712 |
nbactif++; |
| 713 |
if (no3->get_solution()>0.) nbpositif++; |
| 714 |
} |
| 715 |
if ((no4->get_etat()==1) || (no4->get_etat()==-1)) |
| 716 |
{ |
| 717 |
nbactif++; |
| 718 |
if (no4->get_solution()>0.) nbpositif++; |
| 719 |
} |
| 720 |
if ((nbactif==3) && (nbpositif<3) && (nbpositif>0)) |
| 721 |
{ |
| 722 |
couche.push_back(tet); |
| 723 |
} |
| 724 |
} |
| 725 |
for (int i=0;i<couche.size();i++) |
| 726 |
{ |
| 727 |
FEM_TETRA* tet=couche[i]; |
| 728 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 729 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 730 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 731 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 732 |
TOIBREP_POINT* pt; |
| 733 |
FEM_NOEUD* noeudref; |
| 734 |
if (no1->get_etat()==0) {noeudref=no1;pt=(*lst)[no2->get_numero()];} |
| 735 |
if (no2->get_etat()==0) {noeudref=no2;pt=(*lst)[no1->get_numero()];} |
| 736 |
if (no3->get_etat()==0) {noeudref=no3;pt=(*lst)[no1->get_numero()];} |
| 737 |
if (no4->get_etat()==0) {noeudref=no4;pt=(*lst)[no1->get_numero()];} |
| 738 |
MG_FACE* face=pt->get_mg_face(); |
| 739 |
double t,xyz[3]; |
| 740 |
double dist=calcul_distance_level_ortho(noeudref,are,pt,t,xyz); |
| 741 |
noeudref->change_solution(dist); |
| 742 |
noeudref->change_etat(2); |
| 743 |
TOIBREP_POINT *nvpt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],t,0.,1,NULL); |
| 744 |
lst->push_back(nvpt); |
| 745 |
noeudref->change_numero(nvpt->get_id()); |
| 746 |
} |
| 747 |
for (FEM_TETRA* tet=lst_tetra_concerne.get_premier(ittetra);tet!=NULL;tet=lst_tetra_concerne.get_suivant(ittetra)) |
| 748 |
{ |
| 749 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 750 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 751 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 752 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 753 |
if (no1->get_etat()!=0) |
| 754 |
if (no2->get_etat()!=0) |
| 755 |
if (no3->get_etat()!=0) |
| 756 |
if (no4->get_etat()!=0) |
| 757 |
if (tet->get_etat()==0) |
| 758 |
{ |
| 759 |
int reactive=0; |
| 760 |
for (int i=0;i<tet->get_nb_fem_noeud();i++) |
| 761 |
{ |
| 762 |
if (reactive) continue; |
| 763 |
FEM_NOEUD* no=tet->get_fem_noeud(i); |
| 764 |
TOIBREP_POINT *pt=(*lst)[no->get_numero()]; |
| 765 |
double t; |
| 766 |
pt->get_coord1(t); |
| 767 |
double eps=1e-4*(are->get_tmax()-are->get_tmin()); |
| 768 |
if ((t>=are->get_tmin()-eps) && (t<are->get_tmax()+eps)) |
| 769 |
{ |
| 770 |
tet->change_etat(1); |
| 771 |
no1->change_etat(1); |
| 772 |
no2->change_etat(1); |
| 773 |
no3->change_etat(1); |
| 774 |
no4->change_etat(1); |
| 775 |
reactive=1; |
| 776 |
couche_active=1; |
| 777 |
|
| 778 |
} |
| 779 |
|
| 780 |
} |
| 781 |
if (!reactive) |
| 782 |
{ |
| 783 |
if (no1->get_etat()==2) no1->change_etat(-1); |
| 784 |
if (no2->get_etat()==2) no2->change_etat(-1); |
| 785 |
if (no3->get_etat()==2) no3->change_etat(-1); |
| 786 |
if (no4->get_etat()==2) no4->change_etat(-1); |
| 787 |
tet->change_etat(-1); |
| 788 |
} |
| 789 |
} |
| 790 |
} |
| 791 |
} |
| 792 |
while (couche_active==1); |
| 793 |
TPL_MAP_ENTITE<FEM_TETRA*>::ITERATEUR ittetra; |
| 794 |
for (FEM_TETRA* tet=lst_tetra_concerne.get_premier(ittetra);tet!=NULL;tet=lst_tetra_concerne.get_suivant(ittetra)) |
| 795 |
{ |
| 796 |
if (tet->get_etat()!=1) continue; |
| 797 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 798 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 799 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 800 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 801 |
int nbpositif=0; |
| 802 |
int nbnegatif=0; |
| 803 |
double epsdist=1e-6*tet->get_boite_3D().get_rayon(); |
| 804 |
if (no1->get_solution()>epsdist) nbpositif++; |
| 805 |
if (no1->get_solution()<-epsdist) nbnegatif++; |
| 806 |
if (no2->get_solution()>epsdist) nbpositif++; |
| 807 |
if (no2->get_solution()<-epsdist) nbnegatif++; |
| 808 |
if (no3->get_solution()>epsdist) nbpositif++; |
| 809 |
if (no3->get_solution()<-epsdist) nbnegatif++; |
| 810 |
if (no4->get_solution()>epsdist) nbpositif++; |
| 811 |
if (no4->get_solution()<-epsdist) nbnegatif++; |
| 812 |
if (!(((nbnegatif>0) && (nbpositif>0)) ||(nbnegatif+nbpositif==1))) |
| 813 |
tet->change_etat(0); |
| 814 |
} |
| 815 |
} |
| 816 |
|
| 817 |
|
| 818 |
void TOIBREP::remplir_trou(std::vector<TOIBREP_POINT*> *lst,MG_FACE* face,int sens) |
| 819 |
{ |
| 820 |
int castot=0; |
| 821 |
int cas=0; |
| 822 |
int couche_active; |
| 823 |
int nb_couche=0; |
| 824 |
LISTE_FEM_TETRA::iterator ittetra; |
| 825 |
do |
| 826 |
{ |
| 827 |
nb_couche++; |
| 828 |
|
| 829 |
std::vector<FEM_TETRA*> couche; |
| 830 |
couche_active=0; |
| 831 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra)) |
| 832 |
{ |
| 833 |
if (tet->get_etat()!=0) continue; |
| 834 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 835 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 836 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 837 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 838 |
int nbactif=0; |
| 839 |
int nbpositif=0; |
| 840 |
if ((no1->get_etat()==1) || (no1->get_etat()==-1)) |
| 841 |
{ |
| 842 |
nbactif++; |
| 843 |
if (no1->get_solution()>0.) nbpositif++; |
| 844 |
} |
| 845 |
if ((no2->get_etat()==1) || (no2->get_etat()==-1)) |
| 846 |
{ |
| 847 |
nbactif++; |
| 848 |
if (no2->get_solution()>0.) nbpositif++; |
| 849 |
} |
| 850 |
if ((no3->get_etat()==1) || (no3->get_etat()==-1)) |
| 851 |
{ |
| 852 |
nbactif++; |
| 853 |
if (no3->get_solution()>0.) nbpositif++; |
| 854 |
} |
| 855 |
if ((no4->get_etat()==1) || (no4->get_etat()==-1)) |
| 856 |
{ |
| 857 |
nbactif++; |
| 858 |
if (no4->get_solution()>0.) nbpositif++; |
| 859 |
} |
| 860 |
if ((nbactif==3) && (nbpositif<3) && (nbpositif>0)) |
| 861 |
{ |
| 862 |
couche.push_back(tet); |
| 863 |
} |
| 864 |
} |
| 865 |
for (int i=0;i<couche.size();i++) |
| 866 |
{ |
| 867 |
FEM_TETRA* tet=couche[i]; |
| 868 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 869 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 870 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 871 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 872 |
TOIBREP_POINT* pt; |
| 873 |
FEM_NOEUD* noeudref; |
| 874 |
if (no1->get_etat()==0) {noeudref=no1;pt=(*lst)[no2->get_numero()];} |
| 875 |
if (no2->get_etat()==0) {noeudref=no2;pt=(*lst)[no1->get_numero()];} |
| 876 |
if (no3->get_etat()==0) {noeudref=no3;pt=(*lst)[no1->get_numero()];} |
| 877 |
if (no4->get_etat()==0) {noeudref=no4;pt=(*lst)[no1->get_numero()];} |
| 878 |
MG_FACE* face=pt->get_mg_face(); |
| 879 |
double uv[2]; |
| 880 |
double xyz[3]; |
| 881 |
double dist=calcul_distance(noeudref,face,pt,xyz,uv); |
| 882 |
double normal[3]; |
| 883 |
face->calcul_normale_unitaire(uv,normal); |
| 884 |
normal[0]=normal[0]*sens*(-1.); |
| 885 |
normal[1]=normal[1]*sens*(-1.); |
| 886 |
normal[2]=normal[2]*sens*(-1.); |
| 887 |
double dist2=calculdist(normal,xyz[0],xyz[1],xyz[2],noeudref); |
| 888 |
if (dist2<0.) dist=dist*(-1); |
| 889 |
noeudref->change_solution(dist); |
| 890 |
noeudref->change_etat(2); |
| 891 |
TOIBREP_POINT *nvpt=new TOIBREP_POINT(xyz[0],xyz[1],xyz[2],uv[0],uv[1],1,face); |
| 892 |
lst->push_back(nvpt); |
| 893 |
noeudref->change_numero(nvpt->get_id()); |
| 894 |
|
| 895 |
} |
| 896 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra)) |
| 897 |
{ |
| 898 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 899 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 900 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 901 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 902 |
if (no1->get_etat()!=0) |
| 903 |
if (no2->get_etat()!=0) |
| 904 |
if (no3->get_etat()!=0) |
| 905 |
if (no4->get_etat()!=0) |
| 906 |
if (tet->get_etat()==0) |
| 907 |
{ |
| 908 |
castot++; |
| 909 |
double lstptdecoup[5*3]; |
| 910 |
int nb=0; |
| 911 |
decoupe_tetra(tet,lstptdecoup,&nb); |
| 912 |
int reactive=0; |
| 913 |
for (int i=0;i<nb;i++) |
| 914 |
{ |
| 915 |
if (reactive) continue; |
| 916 |
double uv[2]; |
| 917 |
face->inverser(uv,lstptdecoup+3*i); |
| 918 |
if ((face->valide_parametre_u(uv[0])) && (face->valide_parametre_v(uv[1]))) |
| 919 |
{ |
| 920 |
double distance=ot.calcule_distance_contour_face_uv(uv,face); |
| 921 |
if (distance>1e-6) |
| 922 |
{ |
| 923 |
cas++; |
| 924 |
tet->change_etat(1); |
| 925 |
no1->change_etat(1); |
| 926 |
no2->change_etat(1); |
| 927 |
no3->change_etat(1); |
| 928 |
no4->change_etat(1); |
| 929 |
reactive=1; |
| 930 |
couche_active=1; |
| 931 |
} |
| 932 |
} |
| 933 |
} |
| 934 |
if (!reactive) |
| 935 |
{ |
| 936 |
if (no1->get_etat()==2) no1->change_etat(-1); |
| 937 |
if (no2->get_etat()==2) no2->change_etat(-1); |
| 938 |
if (no3->get_etat()==2) no3->change_etat(-1); |
| 939 |
if (no4->get_etat()==2) no4->change_etat(-1); |
| 940 |
tet->change_etat(-1); |
| 941 |
} |
| 942 |
} |
| 943 |
} |
| 944 |
} |
| 945 |
while (couche_active==1); |
| 946 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittetra);tet!=NULL;tet=mai->get_suivant_tetra(ittetra)) |
| 947 |
{ |
| 948 |
if (tet->get_etat()!=1) continue; |
| 949 |
FEM_NOEUD *no1=tet->get_fem_noeud(0); |
| 950 |
FEM_NOEUD *no2=tet->get_fem_noeud(1); |
| 951 |
FEM_NOEUD *no3=tet->get_fem_noeud(2); |
| 952 |
FEM_NOEUD *no4=tet->get_fem_noeud(3); |
| 953 |
int nbpositif=0; |
| 954 |
int nbnegatif=0; |
| 955 |
double epsdist=1e-6*tet->get_boite_3D().get_rayon(); |
| 956 |
if (no1->get_solution()>epsdist) nbpositif++; |
| 957 |
if (no1->get_solution()<-epsdist) nbnegatif++; |
| 958 |
if (no2->get_solution()>epsdist) nbpositif++; |
| 959 |
if (no2->get_solution()<-epsdist) nbnegatif++; |
| 960 |
if (no3->get_solution()>epsdist) nbpositif++; |
| 961 |
if (no3->get_solution()<-epsdist) nbnegatif++; |
| 962 |
if (no4->get_solution()>epsdist) nbpositif++; |
| 963 |
if (no4->get_solution()<-epsdist) nbnegatif++; |
| 964 |
if (!(((nbnegatif>0) && (nbpositif>0)) ||(nbnegatif+nbpositif==1))) |
| 965 |
tet->change_etat(0); |
| 966 |
} |
| 967 |
double res=cas*100.0/castot; |
| 968 |
char mess[255]; |
| 969 |
sprintf(mess," %lf%% en %d couches",res,nb_couche); |
| 970 |
affiche(mess); |
| 971 |
} |
| 972 |
|
| 973 |
void TOIBREP::decoupe_tetra(FEM_TETRA* tet,double* lst,int *nb) |
| 974 |
{ |
| 975 |
FEM_NOEUD* no1=tet->get_fem_noeud(0); |
| 976 |
FEM_NOEUD* no2=tet->get_fem_noeud(1); |
| 977 |
FEM_NOEUD* no3=tet->get_fem_noeud(2); |
| 978 |
FEM_NOEUD* no4=tet->get_fem_noeud(3); |
| 979 |
if (no1->get_solution()*no2->get_solution()<-0.00000001) |
| 980 |
{ |
| 981 |
double t=no1->get_solution()/(no1->get_solution()-no2->get_solution()); |
| 982 |
double x=no1->get_x()+t*(no2->get_x()-no1->get_x()); |
| 983 |
double y=no1->get_y()+t*(no2->get_y()-no1->get_y()); |
| 984 |
double z=no1->get_z()+t*(no2->get_z()-no1->get_z()); |
| 985 |
lst[3*(*nb)]=x; |
| 986 |
lst[3*(*nb)+1]=y; |
| 987 |
lst[3*(*nb)+2]=z; |
| 988 |
(*nb)++; |
| 989 |
} |
| 990 |
if (no1->get_solution()*no3->get_solution()<-0.00000001) |
| 991 |
{ |
| 992 |
double t=no1->get_solution()/(no1->get_solution()-no3->get_solution()); |
| 993 |
double x=no1->get_x()+t*(no3->get_x()-no1->get_x()); |
| 994 |
double y=no1->get_y()+t*(no3->get_y()-no1->get_y()); |
| 995 |
double z=no1->get_z()+t*(no3->get_z()-no1->get_z()); |
| 996 |
lst[3*(*nb)]=x; |
| 997 |
lst[3*(*nb)+1]=y; |
| 998 |
lst[3*(*nb)+2]=z; |
| 999 |
(*nb)++; |
| 1000 |
} |
| 1001 |
if (no1->get_solution()*no4->get_solution()<-0.00000001) |
| 1002 |
{ |
| 1003 |
double t=no1->get_solution()/(no1->get_solution()-no4->get_solution()); |
| 1004 |
double x=no1->get_x()+t*(no4->get_x()-no1->get_x()); |
| 1005 |
double y=no1->get_y()+t*(no4->get_y()-no1->get_y()); |
| 1006 |
double z=no1->get_z()+t*(no4->get_z()-no1->get_z()); |
| 1007 |
lst[3*(*nb)]=x; |
| 1008 |
lst[3*(*nb)+1]=y; |
| 1009 |
lst[3*(*nb)+2]=z; |
| 1010 |
(*nb)++; |
| 1011 |
} |
| 1012 |
if (no2->get_solution()*no3->get_solution()<-0.00000001) |
| 1013 |
{ |
| 1014 |
double t=no2->get_solution()/(no2->get_solution()-no3->get_solution()); |
| 1015 |
double x=no2->get_x()+t*(no3->get_x()-no2->get_x()); |
| 1016 |
double y=no2->get_y()+t*(no3->get_y()-no2->get_y()); |
| 1017 |
double z=no2->get_z()+t*(no3->get_z()-no2->get_z()); |
| 1018 |
lst[3*(*nb)]=x; |
| 1019 |
lst[3*(*nb)+1]=y; |
| 1020 |
lst[3*(*nb)+2]=z; |
| 1021 |
(*nb)++; |
| 1022 |
} |
| 1023 |
if (no2->get_solution()*no4->get_solution()<-0.00000001) |
| 1024 |
{ |
| 1025 |
double t=no2->get_solution()/(no2->get_solution()-no4->get_solution()); |
| 1026 |
double x=no2->get_x()+t*(no4->get_x()-no2->get_x()); |
| 1027 |
double y=no2->get_y()+t*(no4->get_y()-no2->get_y()); |
| 1028 |
double z=no2->get_z()+t*(no4->get_z()-no2->get_z()); |
| 1029 |
lst[3*(*nb)]=x; |
| 1030 |
lst[3*(*nb)+1]=y; |
| 1031 |
lst[3*(*nb)+2]=z; |
| 1032 |
(*nb)++; |
| 1033 |
} |
| 1034 |
if (no3->get_solution()*no4->get_solution()<-0.00000001) |
| 1035 |
{ |
| 1036 |
double t=no3->get_solution()/(no3->get_solution()-no4->get_solution()); |
| 1037 |
double x=no3->get_x()+t*(no4->get_x()-no3->get_x()); |
| 1038 |
double y=no3->get_y()+t*(no4->get_y()-no3->get_y()); |
| 1039 |
double z=no3->get_z()+t*(no4->get_z()-no3->get_z()); |
| 1040 |
lst[3*(*nb)]=x; |
| 1041 |
lst[3*(*nb)+1]=y; |
| 1042 |
lst[3*(*nb)+2]=z; |
| 1043 |
(*nb)++; |
| 1044 |
} |
| 1045 |
} |
| 1046 |
|
| 1047 |
void TOIBREP::calcullevelsetpremierepasse(MG_FACE* face,int sens,vector<TOIBREP_POINT*> *lst,int n1,int n2) |
| 1048 |
{ |
| 1049 |
/*LISTE_FEM_TETRA::iterator it3; |
| 1050 |
for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(it3)) |
| 1051 |
{ |
| 1052 |
double *xyz1,*xyz2,*xyz3,*xyz4; |
| 1053 |
if (tet->get_nb_fem_noeud()==4) |
| 1054 |
{ |
| 1055 |
xyz1=tet->get_fem_noeud(0)->get_coord(); |
| 1056 |
xyz2=tet->get_fem_noeud(1)->get_coord(); |
| 1057 |
xyz3=tet->get_fem_noeud(2)->get_coord(); |
| 1058 |
xyz4=tet->get_fem_noeud(3)->get_coord(); |
| 1059 |
} |
| 1060 |
if (tet->get_nb_fem_noeud()==10) |
| 1061 |
{ |
| 1062 |
xyz1=tet->get_fem_noeud(0)->get_coord(); |
| 1063 |
xyz2=tet->get_fem_noeud(2)->get_coord(); |
| 1064 |
xyz3=tet->get_fem_noeud(4)->get_coord(); |
| 1065 |
xyz4=tet->get_fem_noeud(9)->get_coord(); |
| 1066 |
} |
| 1067 |
double xcentre=0.25*(xyz1[0]+xyz2[0]+xyz3[0]+xyz4[0]); |
| 1068 |
double ycentre=0.25*(xyz1[1]+xyz2[1]+xyz3[1]+xyz4[1]); |
| 1069 |
double zcentre=0.25*(xyz1[2]+xyz2[2]+xyz3[2]+xyz4[2]); |
| 1070 |
double rayon1=(xyz1[0]-xcentre)*(xyz1[0]-xcentre)+(xyz1[1]-ycentre)*(xyz1[1]-ycentre)+(xyz1[2]-zcentre)*(xyz1[2]-zcentre); |
| 1071 |
double rayon2=(xyz2[0]-xcentre)*(xyz2[0]-xcentre)+(xyz2[1]-ycentre)*(xyz2[1]-ycentre)+(xyz2[2]-zcentre)*(xyz2[2]-zcentre); |
| 1072 |
double rayon3=(xyz3[0]-xcentre)*(xyz3[0]-xcentre)+(xyz3[1]-ycentre)*(xyz3[1]-ycentre)+(xyz3[2]-zcentre)*(xyz3[2]-zcentre); |
| 1073 |
double rayon4=(xyz4[0]-xcentre)*(xyz4[0]-xcentre)+(xyz4[1]-ycentre)*(xyz4[1]-ycentre)+(xyz4[2]-zcentre)*(xyz4[2]-zcentre); |
| 1074 |
double rayon=std::max(rayonETRA_ACTIF1,std::max(rayon2,std::max(rayon3,rayon4))); |
| 1075 |
rayon=sqrt(rayon); |
| 1076 |
TPL_MAP_ENTITE<TOIBREP_POINT*> liste; |
| 1077 |
octree.rechercher(xcentre,ycentre,zcentre,rayon,liste); |
| 1078 |
TPL_MAP_ENTITE<TOIBREP_POINT*>::ITERATEUR it; |
| 1079 |
for (TOIBREP_POINT* pt=liste.get_premier(it);pt!=NULL;pt=liste.get_suivant(it)) |
| 1080 |
{ |
| 1081 |
double xyzpt[3]; |
| 1082 |
pt->get_coord3(xyzpt); |
| 1083 |
int etat=estdansletetra(tet,xyzpt[0],xyzpt[1],xyzpt[2]); |
| 1084 |
if (compare_etat(etat,INTERIEUR)) |
| 1085 |
{ |
| 1086 |
for (int k=0;k<tet->get_nb_fem_noeud();k++) |
| 1087 |
{ |
| 1088 |
double normal[3]; |
| 1089 |
double uv[2]; |
| 1090 |
face->inverser(uv,xyzpt); |
| 1091 |
face->calcul_normale_unitaire(uv,normal); |
| 1092 |
normal[0]=normal[0]*sens*(-1.); |
| 1093 |
normal[1]=normal[1]*sens*(-1.); |
| 1094 |
normal[2]=normal[2]*sens*(-1.); |
| 1095 |
double dist=calculdist(normal,xyzpt[0],xyzpt[1],xyzpt[2],tet->get_fem_noeud(k)); |
| 1096 |
if (fabs(dist)<fabs(tet->get_fem_noeud(k)->get_solution())) |
| 1097 |
{ |
| 1098 |
tet->get_fem_noeud(k)->change_solution(dist); |
| 1099 |
tet->get_fem_noeud(k)->change_numero(pt->get_id()); |
| 1100 |
pt->change_coord2(uv); |
| 1101 |
} |
| 1102 |
} |
| 1103 |
} |
| 1104 |
if ((compare_etat(etat,STRICTINTERIEUR)) && (pt->get_interieur())) |
| 1105 |
{ |
| 1106 |
tet->change_etat(1); |
| 1107 |
for (int k=0;k<tet->get_nb_fem_noeud();k++) |
| 1108 |
tet->get_fem_noeud(k)->change_etat(1); |
| 1109 |
} |
| 1110 |
if (tet->get_etat()==1) break; |
| 1111 |
} |
| 1112 |
|
| 1113 |
} |
| 1114 |
octree.vide(); |
| 1115 |
it}*/ |
| 1116 |
/*LISTE_FEM_TETRA::iterator it3; |
| 1117 |
for (FEM_TETRA* tet=mai->get_premier_tetra(it3);tet!=NULL;tet=mai->get_suivant_tetra(it3)) |
| 1118 |
octree.inserer(tet);*/ |
| 1119 |
for (int i=n1;i<n2;i++) |
| 1120 |
{ |
| 1121 |
TOIBREP_POINT* pt=(*lst)[i]; |
| 1122 |
double uv[2]; |
| 1123 |
double xyz[3]; |
| 1124 |
pt->get_coord3(xyz); |
| 1125 |
double normal[3]; |
| 1126 |
face->inverser(uv,xyz); |
| 1127 |
face->calcul_normale_unitaire(uv,normal); |
| 1128 |
normal[0]=normal[0]*sens*(-1.); |
| 1129 |
normal[1]=normal[1]*sens*(-1.); |
| 1130 |
normal[2]=normal[2]*sens*(-1.); |
| 1131 |
TPL_MAP_ENTITE<FEM_TETRA*> liste; |
| 1132 |
octree.rechercher(xyz[0],xyz[1],xyz[2],0.,liste); |
| 1133 |
TPL_MAP_ENTITE<FEM_TETRA*>::ITERATEUR it; |
| 1134 |
for (FEM_TETRA* tet=liste.get_premier(it);tet!=NULL;tet=liste.get_suivant(it)) |
| 1135 |
{ |
| 1136 |
if (tet->get_etat()==1) continue; |
| 1137 |
int etat=estdansletetra(tet,xyz[0],xyz[1],xyz[2]); |
| 1138 |
if ((compare_etat(etat,SUR_FACE)) && (pt->get_interieur())) |
| 1139 |
{ |
| 1140 |
int nbeps=pt->get_nb_point_eps(); |
| 1141 |
for (int k=0;k<nbeps;k++) |
| 1142 |
{ |
| 1143 |
double xyztmp[3]; |
| 1144 |
pt->get_point_eps(k,xyztmp); |
| 1145 |
int etattmp=estdansletetra(tet,xyztmp[0],xyztmp[1],xyztmp[2]); |
| 1146 |
if (compare_etat(etattmp,SUR_FACE)) {etat=STRICTINTERIEUR;break;} |
| 1147 |
if (compare_etat(etattmp,INTERIEUR)) {etat=STRICTINTERIEUR;break;} |
| 1148 |
if (compare_etat(etattmp,STRICTINTERIEUR)) {etat=STRICTINTERIEUR;break;} |
| 1149 |
} |
| 1150 |
} |
| 1151 |
if (compare_etat(etat,INTERIEUR)) |
| 1152 |
{ |
| 1153 |
for (int k=0;k<tet->get_nb_fem_noeud();k++) |
| 1154 |
{ |
| 1155 |
double dist=calculdist(normal,xyz[0],xyz[1],xyz[2],tet->get_fem_noeud(k)); |
| 1156 |
if (fabs(dist)<fabs(tet->get_fem_noeud(k)->get_solution())) |
| 1157 |
{ |
| 1158 |
tet->get_fem_noeud(k)->change_solution(dist); |
| 1159 |
tet->get_fem_noeud(k)->change_numero(pt->get_id()); |
| 1160 |
pt->change_coord2(uv); |
| 1161 |
} |
| 1162 |
} |
| 1163 |
} |
| 1164 |
if ((compare_etat(etat,STRICTINTERIEUR)) && (pt->get_interieur())) |
| 1165 |
{ |
| 1166 |
tet->change_etat(1); |
| 1167 |
for (int k=0;k<tet->get_nb_fem_noeud();k++) |
| 1168 |
{ |
| 1169 |
tet->get_fem_noeud(k)->change_etat(1); |
| 1170 |
//octree.supprimer(tet); |
| 1171 |
} |
| 1172 |
} |
| 1173 |
} |
| 1174 |
} |
| 1175 |
//octree.vide(); |
| 1176 |
|
| 1177 |
|
| 1178 |
|
| 1179 |
|
| 1180 |
} |
| 1181 |
|
| 1182 |
void TOIBREP::calcullevelsetdeuxiemepasse(vector<TOIBREP_POINT*> *lst) |
| 1183 |
{ |
| 1184 |
LISTE_FEM_NOEUD::iterator itno; |
| 1185 |
for (FEM_NOEUD* no=mai->get_premier_noeud(itno);no!=NULL;no=mai->get_suivant_noeud(itno)) |
| 1186 |
if (no->get_solution()<1e250) |
| 1187 |
{ |
| 1188 |
TOIBREP_POINT* pt=(*lst)[no->get_numero()]; |
| 1189 |
MG_FACE* face=pt->get_mg_face(); |
| 1190 |
double xyz[3],uv[2]; |
| 1191 |
double dist=calcul_distance(no,face,pt,xyz,uv); |
| 1192 |
int signe=1; |
| 1193 |
if (no->get_solution()<0) signe=-1; |
| 1194 |
no->change_solution(signe*dist); |
| 1195 |
} |
| 1196 |
|
| 1197 |
|
| 1198 |
|
| 1199 |
|
| 1200 |
} |
| 1201 |
|
| 1202 |
|
| 1203 |
|
| 1204 |
double TOIBREP::calcul_distance(FEM_NOEUD* noeud,MG_ARETE* are,TOIBREP_POINT* pt,double &tii,double *xyz,double precision) |
| 1205 |
{ |
| 1206 |
double eps; |
| 1207 |
are->inverser(tii,noeud->get_coord()); |
| 1208 |
int compteur=0; |
| 1209 |
OT_VECTEUR_3D Pt(noeud->get_x(),noeud->get_y(),noeud->get_z()); |
| 1210 |
do |
| 1211 |
{ |
| 1212 |
compteur++; |
| 1213 |
double ti=tii; |
| 1214 |
double xyz[3],dxyz[3],ddxyz[3]; |
| 1215 |
are->deriver_seconde(ti,ddxyz,dxyz,xyz); |
| 1216 |
OT_VECTEUR_3D Ct(xyz[0],xyz[1],xyz[2]); |
| 1217 |
OT_VECTEUR_3D Ct_deriver(dxyz[0],dxyz[1],dxyz[2]); |
| 1218 |
OT_VECTEUR_3D Ct_deriver_seconde(ddxyz[0],ddxyz[1],ddxyz[2]); |
| 1219 |
OT_VECTEUR_3D Distance = Ct-Pt; |
| 1220 |
tii=ti-Ct_deriver*Distance/(Ct_deriver_seconde*Distance+Ct_deriver.get_longueur2()); |
| 1221 |
eps=fabs(tii-ti); |
| 1222 |
if (compteur>500) return 1e300; |
| 1223 |
if (tii<are->get_tmin()) |
| 1224 |
{ |
| 1225 |
tii=are->get_tmin(); |
| 1226 |
eps=0.; |
| 1227 |
} |
| 1228 |
if (tii>are->get_tmax()) |
| 1229 |
{ |
| 1230 |
tii=are->get_tmax(); |
| 1231 |
eps=0.; |
| 1232 |
} |
| 1233 |
} |
| 1234 |
while (eps>precision); |
| 1235 |
are->evaluer(tii,xyz); |
| 1236 |
OT_VECTEUR_3D Ct(xyz[0],xyz[1],xyz[2]); |
| 1237 |
double distance=(Ct-Pt).get_longueur(); |
| 1238 |
return distance; |
| 1239 |
} |
| 1240 |
|
| 1241 |
|
| 1242 |
double TOIBREP::calcul_distance_level_ortho(FEM_NOEUD* no,MG_ARETE* are,TOIBREP_POINT* pt,double precision) |
| 1243 |
{ |
| 1244 |
double t; |
| 1245 |
double xyz[3]; |
| 1246 |
return calcul_distance_level_ortho(no,are,pt,t,xyz,precision); |
| 1247 |
} |
| 1248 |
|
| 1249 |
double TOIBREP::calcul_distance_level_ortho(FEM_NOEUD* no,MG_ARETE* are,TOIBREP_POINT* pt,double &t,double *xyz,double precision) |
| 1250 |
{ |
| 1251 |
calcul_distance(no,are,pt,t,xyz,precision); |
| 1252 |
double dir[3]; |
| 1253 |
are->deriver(t,dir); |
| 1254 |
double uv[2]; |
| 1255 |
MG_FACE* face=are->get_mg_coarete(0)->get_boucle()->get_mg_face(); |
| 1256 |
face->inverser(uv,xyz); |
| 1257 |
double nor[3]; |
| 1258 |
face->calcul_normale_unitaire(uv,nor); |
| 1259 |
OT_VECTEUR_3D dirv(dir); |
| 1260 |
OT_VECTEUR_3D norv(nor); |
| 1261 |
OT_VECTEUR_3D sensv=norv&dirv; |
| 1262 |
double d=-(sensv.get_x()*xyz[0]+sensv.get_y()*xyz[1],sensv.get_z()*xyz[2]); |
| 1263 |
double dist=fabs(sensv.get_x()*no->get_x()+sensv.get_y()*no->get_y()+sensv.get_z()*no->get_z()+d); |
| 1264 |
OT_VECTEUR_3D distv(no->get_x()-xyz[0],no->get_y()-xyz[1],no->get_z()-xyz[2]); |
| 1265 |
if (sensv*distv<0.) dist=-dist; |
| 1266 |
if (fabs(fabs(dist)-fabs(no->get_y()))>1e-10) |
| 1267 |
std::cout << "pb calcul" << std::endl; |
| 1268 |
return dist; |
| 1269 |
} |
| 1270 |
|
| 1271 |
|
| 1272 |
double TOIBREP::calcul_distance(FEM_NOEUD* noeud,MG_FACE* face,TOIBREP_POINT* pt,double *xyz,double *uv,double precision) |
| 1273 |
{ |
| 1274 |
double uvii[2],eps; |
| 1275 |
pt->get_coord2(uvii); |
| 1276 |
int compteur=0; |
| 1277 |
OT_VECTEUR_3D Pt(noeud->get_x(),noeud->get_y(),noeud->get_z()); |
| 1278 |
double delta_u,delta_v; |
| 1279 |
do |
| 1280 |
{ |
| 1281 |
compteur++; |
| 1282 |
double uvi[2]; |
| 1283 |
uvi[0]=uvii[0]; |
| 1284 |
uvi[1]=uvii[1]; |
| 1285 |
double xyzduu[3],xyzdvv[3],xyzduv[3],xyzdu[3],xyzdv[3],xyz[3]; |
| 1286 |
face->deriver_seconde(uvi,xyzduu,xyzduv,xyzdvv,xyz,xyzdu,xyzdv); |
| 1287 |
OT_VECTEUR_3D S(xyz[0],xyz[1],xyz[2]); |
| 1288 |
OT_VECTEUR_3D Su(xyzdu[0],xyzdu[1],xyzdu[2]); |
| 1289 |
OT_VECTEUR_3D Sv(xyzdv[0],xyzdv[1],xyzdv[2]); |
| 1290 |
OT_VECTEUR_3D Suu(xyzduu[0],xyzduu[1],xyzduu[2]); |
| 1291 |
OT_VECTEUR_3D Suv(xyzduv[0],xyzduv[1],xyzduv[2]); |
| 1292 |
OT_VECTEUR_3D Svv(xyzdvv[0],xyzdvv[1],xyzdvv[2]); |
| 1293 |
OT_VECTEUR_3D Distance = S-Pt; |
| 1294 |
double a[4],b[2]; |
| 1295 |
a[0]=Su.get_longueur2()+Distance*Suu; |
| 1296 |
a[1]=Su*Sv+Distance*Suv; |
| 1297 |
a[2]=Su*Sv+Distance*Suv; |
| 1298 |
a[3]=Sv.get_longueur2()+Distance*Svv; |
| 1299 |
b[0]=Distance*Su;b[0]=-b[0]; |
| 1300 |
b[1]=Distance*Sv;b[1]=-b[1]; |
| 1301 |
double deltau,deltav; |
| 1302 |
double denominateur_delta=(a[0]*a[3]-a[2]*a[1]); |
| 1303 |
if (a[0]<1E-12) |
| 1304 |
deltau=0; |
| 1305 |
else delta_u=(b[0]*a[3]-b[1]*a[1])/denominateur_delta; |
| 1306 |
if (a[3]<1E-12) |
| 1307 |
deltav=0; |
| 1308 |
else delta_v=(a[0]*b[1]-a[2]*b[0])/denominateur_delta; |
| 1309 |
/*if (fabs(denominateur_delta) < ( (fabs(a[0])+fabs(a[1])+fabs(a[2])+fabs(a[3]))*1e-12 ) ) |
| 1310 |
return 1e300;*/ |
| 1311 |
uvii[0]=uvi[0]+delta_u; |
| 1312 |
uvii[1]=uvi[1]+delta_v; |
| 1313 |
if (face->get_surface()->est_periodique_u()==1) |
| 1314 |
{ |
| 1315 |
if(uvii[0]<0.) uvii[0]=face->get_surface()->get_periode_u()-uvii[0]; |
| 1316 |
if(uvii[0]>face->get_surface()->get_periode_u()) uvii[0]=uvii[0]-face->get_surface()->get_periode_u(); |
| 1317 |
} |
| 1318 |
if (face->get_surface()->est_periodique_v()==1) |
| 1319 |
{ |
| 1320 |
if(uvii[1]<0.) uvii[0]=face->get_surface()->get_periode_v()-uvii[1]; |
| 1321 |
if(uvii[1]>face->get_surface()->get_periode_v()) uvii[1]=uvii[1]-face->get_surface()->get_periode_v(); |
| 1322 |
} |
| 1323 |
delta_u=uvii[0]-uvi[0]; |
| 1324 |
delta_v=uvii[1]-uvi[1]; |
| 1325 |
if (compteur>500) return 1e300; |
| 1326 |
} |
| 1327 |
|
| 1328 |
while ((fabs(delta_u)>precision)||(fabs(delta_v)>precision)); |
| 1329 |
face->evaluer(uvii,xyz); |
| 1330 |
OT_VECTEUR_3D S(xyz[0],xyz[1],xyz[2]); |
| 1331 |
double distance=(S-Pt).get_longueur(); |
| 1332 |
uv[0]=uvii[0]; |
| 1333 |
uv[1]=uvii[1]; |
| 1334 |
return distance; |
| 1335 |
} |
| 1336 |
|
| 1337 |
|
| 1338 |
|
| 1339 |
void TOIBREP::etendrelevelset(FEM_SOLUTION* sol,int numsol) |
| 1340 |
{ |
| 1341 |
sol->active_solution(numsol); |
| 1342 |
|
| 1343 |
|
| 1344 |
LISTE_FM know; |
| 1345 |
LISTE_FM trial; |
| 1346 |
LISTE_FM far; |
| 1347 |
LISTE_FM exterieur; |
| 1348 |
LISTE_FM_TRI trialtri; |
| 1349 |
LISTE_FM_TRI_ID trialtriid; |
| 1350 |
|
| 1351 |
|
| 1352 |
LISTE_FEM_TETRA::iterator ittet; |
| 1353 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet)) |
| 1354 |
{ |
| 1355 |
tet->change_solution(0.); |
| 1356 |
int numsol=0; |
| 1357 |
if (tet->get_fem_noeud(0)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(0)->change_numero(1);} else tet->get_fem_noeud(0)->change_numero(0); |
| 1358 |
if (tet->get_fem_noeud(1)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(1)->change_numero(1);} else tet->get_fem_noeud(1)->change_numero(0); |
| 1359 |
if (tet->get_fem_noeud(2)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(2)->change_numero(1);} else tet->get_fem_noeud(2)->change_numero(0); |
| 1360 |
if (tet->get_fem_noeud(3)->get_solution()<1e200) {numsol++;tet->get_fem_noeud(3)->change_numero(1);} else tet->get_fem_noeud(3)->change_numero(0); |
| 1361 |
if (numsol==4) |
| 1362 |
ajouter_liste(know,tet); |
| 1363 |
|
| 1364 |
else if (numsol==3) |
| 1365 |
{ |
| 1366 |
if (tet->get_fem_noeud(0)->get_numero()==0) tet->change_numero(0); |
| 1367 |
else if (tet->get_fem_noeud(1)->get_numero()==0) tet->change_numero(1); |
| 1368 |
else if (tet->get_fem_noeud(2)->get_numero()==0) tet->change_numero(2); |
| 1369 |
else tet->change_numero(3); |
| 1370 |
ajouter_liste(trial,tet); |
| 1371 |
} |
| 1372 |
else |
| 1373 |
ajouter_liste(far,tet); |
| 1374 |
|
| 1375 |
} |
| 1376 |
for (LISTE_FM::iterator i=trial.begin();i!=trial.end();i++) |
| 1377 |
{ |
| 1378 |
FEM_TETRA* tet=*i; |
| 1379 |
int signe; |
| 1380 |
double sol=resoudgradT(tet,&signe); |
| 1381 |
if (fabs(sol)>0.00000001) |
| 1382 |
{ |
| 1383 |
if (fabs(sol)<fabs(tet->get_fem_noeud(tet->get_numero())->get_solution())) |
| 1384 |
{ |
| 1385 |
tet->get_fem_noeud(tet->get_numero())->change_solution(sol); |
| 1386 |
ajouter_liste(trialtri,trialtriid,tet,fabs(tet->get_fem_noeud(tet->get_numero())->get_solution())); |
| 1387 |
} |
| 1388 |
} |
| 1389 |
else |
| 1390 |
{ |
| 1391 |
ajouter_liste(exterieur,tet); |
| 1392 |
tet->change_solution(1e300); |
| 1393 |
} |
| 1394 |
} |
| 1395 |
int fin=0; |
| 1396 |
LISTE_FM_TRI::iterator itfin=trialtri.end(); |
| 1397 |
itfin--; |
| 1398 |
double longref=(*itfin).first; |
| 1399 |
do |
| 1400 |
{ |
| 1401 |
LISTE_FM_TRI::iterator it=trialtri.begin(); |
| 1402 |
FEM_TETRA* tet=(*it).second; |
| 1403 |
double longcourant=(*it).first; |
| 1404 |
supprimer_liste(trialtri,trialtriid,tet); |
| 1405 |
ajouter_liste(know,tet); |
| 1406 |
FEM_NOEUD* noeud=tet->get_fem_noeud(tet->get_numero()); |
| 1407 |
noeud->change_numero(1); |
| 1408 |
if (noeud->get_solution()>20000) |
| 1409 |
cout << "BUGGGGGGG" <<endl; |
| 1410 |
int nbtetra=noeud->get_lien_tetra()->get_nb(); |
| 1411 |
for (int i=0;i<nbtetra;i++) |
| 1412 |
{ |
| 1413 |
FEM_TETRA* tet2=noeud->get_lien_tetra()->get(i); |
| 1414 |
if (tet2==tet) continue; |
| 1415 |
LISTE_FM_TRI_ID::iterator it=trialtriid.find(tet2->get_id()); |
| 1416 |
if (it!=trialtriid.end()) |
| 1417 |
{ |
| 1418 |
int signe; |
| 1419 |
double sol=resoudgradT(tet2,&signe); |
| 1420 |
double solution=tet2->get_fem_noeud(tet2->get_numero())->get_solution(); |
| 1421 |
if (fabs(sol)>0.00000001) |
| 1422 |
if (!((solution<1e200)&&(sol*solution<0))) |
| 1423 |
if (fabs(sol)<fabs(solution)) |
| 1424 |
{ |
| 1425 |
supprimer_liste(trialtri,trialtriid,tet2); |
| 1426 |
ajouter_liste(trialtri,trialtriid,tet2,fabs(sol)); |
| 1427 |
tet2->get_fem_noeud(tet2->get_numero())->change_solution(sol); |
| 1428 |
} |
| 1429 |
} |
| 1430 |
LISTE_FM::iterator it2=find(far.begin(),far.end(),tet2); |
| 1431 |
if (it2!=far.end()) |
| 1432 |
{ |
| 1433 |
int numsol=0; |
| 1434 |
if (tet2->get_fem_noeud(0)->get_numero()==1) numsol++; |
| 1435 |
if (tet2->get_fem_noeud(1)->get_numero()==1) numsol++; |
| 1436 |
if (tet2->get_fem_noeud(2)->get_numero()==1) numsol++; |
| 1437 |
if (tet2->get_fem_noeud(3)->get_numero()==1) numsol++; |
| 1438 |
//if (numsol==4) |
| 1439 |
//cout << " BUG " <<endl; |
| 1440 |
if (numsol==3) |
| 1441 |
{ |
| 1442 |
if (tet2->get_fem_noeud(0)->get_numero()==0) tet2->change_numero(0); |
| 1443 |
else if (tet2->get_fem_noeud(1)->get_numero()==0) tet2->change_numero(1); |
| 1444 |
else if (tet2->get_fem_noeud(2)->get_numero()==0) tet2->change_numero(2); |
| 1445 |
else tet2->change_numero(3); |
| 1446 |
int signe; |
| 1447 |
double sol=resoudgradT(tet2,&signe); |
| 1448 |
double ancsol=tet2->get_fem_noeud(tet2->get_numero())->get_solution(); |
| 1449 |
if (fabs(sol)>0.00000001) |
| 1450 |
{ |
| 1451 |
if (!((ancsol<1e200) && (ancsol*sol<0.) )) |
| 1452 |
{ |
| 1453 |
tet2->get_fem_noeud(tet2->get_numero())->change_solution(sol); |
| 1454 |
supprimer_liste(far,tet2); |
| 1455 |
ajouter_liste(trialtri,trialtriid,tet2,fabs(sol)); |
| 1456 |
} |
| 1457 |
} |
| 1458 |
else |
| 1459 |
{ |
| 1460 |
tet2->change_solution(1e300); |
| 1461 |
ajouter_liste(exterieur,tet2); |
| 1462 |
} |
| 1463 |
} |
| 1464 |
} |
| 1465 |
} |
| 1466 |
if (trialtri.size()==0) fin=1; |
| 1467 |
if (exterieur.size()>0) |
| 1468 |
if (fin==0) |
| 1469 |
if (longcourant>longref) |
| 1470 |
{ |
| 1471 |
int nombre=exterieur.size(); |
| 1472 |
for (int i=0;i<nombre;i++) |
| 1473 |
{ |
| 1474 |
FEM_TETRA* tet2=(*(exterieur.begin())); |
| 1475 |
supprimer_liste(exterieur,tet2); |
| 1476 |
ajouter_liste(know,tet2); |
| 1477 |
for (int nd=0;nd<4;nd++) |
| 1478 |
{ |
| 1479 |
FEM_NOEUD* noeud=tet2->get_fem_noeud(nd); |
| 1480 |
int nbtetra=noeud->get_lien_tetra()->get_nb(); |
| 1481 |
for (int i=0;i<nbtetra;i++) |
| 1482 |
{ |
| 1483 |
FEM_TETRA* tet3=noeud->get_lien_tetra()->get(i); |
| 1484 |
if (tet2==tet3) continue; |
| 1485 |
LISTE_FM::iterator it2=find(far.begin(),far.end(),tet3); |
| 1486 |
if (it2!=far.end()) |
| 1487 |
{ |
| 1488 |
supprimer_liste(far,tet3); |
| 1489 |
ajouter_liste(exterieur,tet3); |
| 1490 |
tet3->change_solution(1e300); |
| 1491 |
} |
| 1492 |
} |
| 1493 |
} |
| 1494 |
} |
| 1495 |
LISTE_FM_TRI::iterator itfin=trialtri.end(); |
| 1496 |
itfin--; |
| 1497 |
longref=(*itfin).first; |
| 1498 |
} |
| 1499 |
} |
| 1500 |
while (fin==0); |
| 1501 |
LISTE_FEM_NOEUD::iterator itnoeud; |
| 1502 |
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud)) |
| 1503 |
noeud->change_numero(0); |
| 1504 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet)) |
| 1505 |
{ |
| 1506 |
if (tet->get_solution()>1e200) continue; |
| 1507 |
tet->get_fem_noeud(0)->change_numero(1); |
| 1508 |
tet->get_fem_noeud(1)->change_numero(1); |
| 1509 |
tet->get_fem_noeud(2)->change_numero(1); |
| 1510 |
tet->get_fem_noeud(3)->change_numero(1); |
| 1511 |
} |
| 1512 |
int i=0; |
| 1513 |
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud)) |
| 1514 |
{ |
| 1515 |
if (noeud->get_numero()==1) sol->ecrire(i,numsol,noeud->get_solution()); else sol->ecrire(i,numsol,1e300); |
| 1516 |
++i; |
| 1517 |
} |
| 1518 |
} |
| 1519 |
IBrep TOIBREP::exporter_IBrep(string chemin,FEM_SOLUTION* solution,FEM_SOLUTION *solution_ele,MG_GROUPE_TOPOLOGIQUE* mggt) |
| 1520 |
{ |
| 1521 |
TPL_MAP_ENTITE<MG_VOLUME*> lst; |
| 1522 |
if (mggt==NULL) |
| 1523 |
{ |
| 1524 |
int nbvol=geo->get_nb_mg_volume(); |
| 1525 |
for (int i=0;i<nbvol;i++) |
| 1526 |
lst.ajouter(geo->get_mg_volume(i)); |
| 1527 |
} |
| 1528 |
else |
| 1529 |
{ |
| 1530 |
int nbmggt=mggt->get_nb(); |
| 1531 |
for (int i=0;i<nbmggt;i++) |
| 1532 |
if (mggt->get(i)->get_dimension()==3) |
| 1533 |
lst.ajouter((MG_VOLUME*)mggt->get(i)); |
| 1534 |
} |
| 1535 |
IBrep brep(100); |
| 1536 |
int nbvolexp=lst.get_nb(); |
| 1537 |
for (int i=0;i<nbvolexp;i++) |
| 1538 |
{ |
| 1539 |
MG_VOLUME* vol=lst.get(i); |
| 1540 |
int nbcoq=vol->get_nb_mg_coquille(); |
| 1541 |
IVolumeN ivoltmp(vol->get_id()); |
| 1542 |
pair<IVolume*,bool> pair_ivol=brep.AddVolume(ivoltmp); |
| 1543 |
if (vol->get_idoriginal()!="") brep.SetName(pair_ivol.first,vol->get_idoriginal()); |
| 1544 |
IVolume* ivol=pair_ivol.first; |
| 1545 |
for (int j=0;j<nbcoq;j++) |
| 1546 |
{ |
| 1547 |
MG_COQUILLE* coq=vol->get_mg_coquille(j); |
| 1548 |
int nbface=coq->get_nb_mg_coface(); |
| 1549 |
IShell ishell(nbface); |
| 1550 |
for (int k=0;k<nbface;k++) |
| 1551 |
{ |
| 1552 |
MG_COFACE* coface=coq->get_mg_coface(k); |
| 1553 |
ishell[k]=coface->get_id(); |
| 1554 |
MG_FACE* face=coface->get_face(); |
| 1555 |
int sens=coface->get_orientation(); |
| 1556 |
ICoFaceN icoface(coface->get_id(),face->get_id(),sens); |
| 1557 |
brep.AddCoFace(icoface); |
| 1558 |
IFace* iface=brep.GetFace(face->get_id()); |
| 1559 |
if (iface==NULL) |
| 1560 |
{ |
| 1561 |
IFaceN ifacenew(face->get_id()); |
| 1562 |
pair<IFace*,bool> pair_iface=brep.AddFace(ifacenew); |
| 1563 |
if (face->get_idoriginal()!="") brep.SetName(pair_iface.first,face->get_idoriginal()); |
| 1564 |
iface=pair_iface.first; |
| 1565 |
int nbbou=face->get_nb_mg_boucle(); |
| 1566 |
for (int l=0;l<nbbou;l++) |
| 1567 |
{ |
| 1568 |
MG_BOUCLE* bou=face->get_mg_boucle(l); |
| 1569 |
int nbare=bou->get_nb_mg_coarete(); |
| 1570 |
ILoop iloop(nbare); |
| 1571 |
for (int m=0;m<nbare;m++) |
| 1572 |
{ |
| 1573 |
MG_COARETE* coare=bou->get_mg_coarete(m); |
| 1574 |
iloop[m]=coare->get_id(); |
| 1575 |
MG_ARETE* are=coare->get_arete(); |
| 1576 |
int sens=coare->get_orientation(); |
| 1577 |
ICoEdgeN icoedge(coare->get_id(),are->get_id(),sens,face->get_id()); |
| 1578 |
brep.AddCoEdge(icoedge); |
| 1579 |
IEdge* iedge=brep.GetEdge(are->get_id()); |
| 1580 |
IVertex *ivertex1,*ivertex2; |
| 1581 |
ICoVertex *icover1,*icover2; |
| 1582 |
if (iedge==NULL) |
| 1583 |
{ |
| 1584 |
MG_COSOMMET* cover1=are->get_cosommet1(); |
| 1585 |
MG_COSOMMET* cover2=are->get_cosommet2(); |
| 1586 |
MG_SOMMET* ver1=cover1->get_sommet(); |
| 1587 |
MG_SOMMET* ver2=cover2->get_sommet(); |
| 1588 |
ICoVertexN icovertmp1(cover1->get_id(),are->get_id(),ver1->get_id(),cover1->get_t()); |
| 1589 |
ICoVertexN icovertmp2(cover2->get_id(),are->get_id(),ver2->get_id(),cover2->get_t()); |
| 1590 |
pair<ICoVertex*,bool> pair_icover1=brep.AddCoVertex(icovertmp1); |
| 1591 |
pair<ICoVertex*,bool> pair_icover2=brep.AddCoVertex(icovertmp2); |
| 1592 |
icover1=pair_icover1.first; |
| 1593 |
icover2=pair_icover2.first; |
| 1594 |
IEdgeN iedgetmp(are->get_id(),cover1->get_id(),cover2->get_id()); |
| 1595 |
ivertex1=brep.GetVertex(ver1->get_id()); |
| 1596 |
ivertex2=brep.GetVertex(ver2->get_id()); |
| 1597 |
if (ivertex1==NULL) |
| 1598 |
{ |
| 1599 |
MG_POINT* pt=ver1->get_point(); |
| 1600 |
double xyz[3]; |
| 1601 |
pt->evaluer(xyz); |
| 1602 |
IVertexN ivertex(ver1->get_id(),xyz); |
| 1603 |
pair<IVertex*,bool> pair_ivertex1=brep.AddVertex(ivertex); |
| 1604 |
if (ver1->get_idoriginal()!="") brep.SetName(pair_ivertex1.first,ver1->get_idoriginal()); |
| 1605 |
ivertex1=pair_ivertex1.first; |
| 1606 |
} |
| 1607 |
if (ivertex2==NULL) |
| 1608 |
{ |
| 1609 |
MG_POINT* pt=ver2->get_point(); |
| 1610 |
double xyz[3]; |
| 1611 |
pt->evaluer(xyz); |
| 1612 |
IVertexN ivertex(ver2->get_id(),xyz); |
| 1613 |
pair<IVertex*,bool> pair_ivertex2=brep.AddVertex(ivertex); |
| 1614 |
if (ver2->get_idoriginal()!="") brep.SetName(pair_ivertex2.first,ver2->get_idoriginal()); |
| 1615 |
ivertex2=pair_ivertex2.first; |
| 1616 |
} |
| 1617 |
ivertex1->AddCoVertex(cover1->get_id()); |
| 1618 |
ivertex2->AddCoVertex(cover2->get_id()); |
| 1619 |
pair<IEdge*,bool> pair_iedge=brep.AddEdge(iedgetmp); |
| 1620 |
if (are->get_idoriginal()!="") brep.SetName(pair_iedge.first,are->get_idoriginal()); |
| 1621 |
iedge=pair_iedge.first; |
| 1622 |
} |
| 1623 |
iedge->AddCoEdge(coare->get_id(),brep); |
| 1624 |
if (sens==1) |
| 1625 |
{ |
| 1626 |
icover1=brep.GetCoVertex(iedge->FromCoVertex()); |
| 1627 |
ivertex1=brep.GetVertex(icover1->Vertex()); |
| 1628 |
ivertex1->AddFace(face->get_id()); |
| 1629 |
} |
| 1630 |
else |
| 1631 |
{ |
| 1632 |
icover2=brep.GetCoVertex(iedge->ToCoVertex()); |
| 1633 |
ivertex2=brep.GetVertex(icover2->Vertex()); |
| 1634 |
ivertex2->AddFace(face->get_id()); |
| 1635 |
} |
| 1636 |
} |
| 1637 |
iface->AddLoop(iloop); |
| 1638 |
} |
| 1639 |
} |
| 1640 |
iface->AddCoFace(coface->get_id()); |
| 1641 |
} |
| 1642 |
ivol->AddShell(ishell); |
| 1643 |
} |
| 1644 |
} |
| 1645 |
int nbsol=solution->get_nb_champ(); |
| 1646 |
//unsigned long *correspondface=new unsigned long[nbsol]; |
| 1647 |
for (int i=0;i<nbsol;i++) |
| 1648 |
{ |
| 1649 |
std::string nom=solution->get_legende(i); |
| 1650 |
char nom2[2]; |
| 1651 |
unsigned long id; |
| 1652 |
sscanf(nom.c_str(),"%s %lu",nom2,&id); |
| 1653 |
std::pair<IField *,bool> fld=brep.AddField(id); |
| 1654 |
if (nom2[0]!='F') |
| 1655 |
{ |
| 1656 |
IField *f=fld.first; |
| 1657 |
f->info.tag=1; |
| 1658 |
} |
| 1659 |
} |
| 1660 |
LISTE_FEM_NOEUD::iterator itnoeud; |
| 1661 |
int i=0; |
| 1662 |
for (FEM_NOEUD* noeud=mai->get_premier_noeud(itnoeud);noeud!=NULL;noeud=mai->get_suivant_noeud(itnoeud)) |
| 1663 |
{ |
| 1664 |
int nbsolnoeud=0; |
| 1665 |
for (int j=0;j<nbsol;j++) |
| 1666 |
{ |
| 1667 |
if (solution->lire(i,j)<1e200) nbsolnoeud++; |
| 1668 |
} |
| 1669 |
double *xyz=noeud->get_coord(); |
| 1670 |
INodeN inoeudtmp(noeud->get_id(),xyz[0],xyz[1],xyz[2],nbsolnoeud); |
| 1671 |
pair<INode*,bool> pair_inoeud=brep.AddNode(inoeudtmp); |
| 1672 |
INode* inoeud=pair_inoeud.first; |
| 1673 |
int num=0; |
| 1674 |
IBrep::iterator_fields itf=brep.begin_fields(); |
| 1675 |
for (int j=0;j<nbsol;j++) |
| 1676 |
{ |
| 1677 |
if (solution->lire(i,j)<1e200) |
| 1678 |
{ |
| 1679 |
inoeud->Fields()[num]=(*itf).num; |
| 1680 |
inoeud->Values()[num]=solution->lire(i,j); |
| 1681 |
num++; |
| 1682 |
} |
| 1683 |
itf++; |
| 1684 |
} |
| 1685 |
i++; |
| 1686 |
} |
| 1687 |
LISTE_FEM_TETRA::iterator ittet; |
| 1688 |
int num=0; |
| 1689 |
i=0; |
| 1690 |
for (FEM_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet)) |
| 1691 |
{ |
| 1692 |
int nbsolele=0; |
| 1693 |
for (int j=0;j<nbsol;j++) |
| 1694 |
{ |
| 1695 |
if (fabs(solution_ele->lire(i,j)-1.)<0.0000000001) nbsolele++; |
| 1696 |
} |
| 1697 |
IElementN xfemele(IElement::TETRAHEDRON,nbsolele); |
| 1698 |
num++; |
| 1699 |
xfemele.num=num; |
| 1700 |
xfemele.GetNode(0)=tet->get_fem_noeud(0)->get_id(); |
| 1701 |
xfemele.GetNode(1)=tet->get_fem_noeud(1)->get_id(); |
| 1702 |
xfemele.GetNode(2)=tet->get_fem_noeud(2)->get_id(); |
| 1703 |
xfemele.GetNode(3)=tet->get_fem_noeud(3)->get_id(); |
| 1704 |
int num=0; |
| 1705 |
IBrep::iterator_fields itf=brep.begin_fields(); |
| 1706 |
for (int j=0;j<nbsol;j++) |
| 1707 |
{ |
| 1708 |
if (fabs(solution_ele->lire(i,j)-1.)<0.0000000001) |
| 1709 |
{ |
| 1710 |
xfemele.Fields()[num]=(*itf).num; |
| 1711 |
num++; |
| 1712 |
} |
| 1713 |
itf++; |
| 1714 |
} |
| 1715 |
brep.AddElement(xfemele); |
| 1716 |
i++; |
| 1717 |
} |
| 1718 |
|
| 1719 |
|
| 1720 |
|
| 1721 |
|
| 1722 |
std::ofstream output(chemin.c_str()); |
| 1723 |
output << brep << std::endl; |
| 1724 |
output.close(); |
| 1725 |
|
| 1726 |
|
| 1727 |
|
| 1728 |
|
| 1729 |
return brep; |
| 1730 |
} |
| 1731 |
|
| 1732 |
|
| 1733 |
|
| 1734 |
void TOIBREP::ajouter_liste(LISTE_FM_TRI &lst,LISTE_FM_TRI_ID &lstid,FEM_TETRA* tet,double val) |
| 1735 |
{ |
| 1736 |
pair<double,FEM_TETRA*> p(val,tet); |
| 1737 |
LISTE_FM_TRI::iterator it=lst.insert(p); |
| 1738 |
lstid[tet->get_id()]=it; |
| 1739 |
} |
| 1740 |
|
| 1741 |
void TOIBREP::supprimer_liste(LISTE_FM_TRI &lst,LISTE_FM_TRI_ID &lstid,FEM_TETRA* tet) |
| 1742 |
{ |
| 1743 |
LISTE_FM_TRI::iterator it2=lstid[tet->get_id()]; |
| 1744 |
LISTE_FM_TRI_ID::iterator it3=lstid.find(tet->get_id()); |
| 1745 |
lstid.erase(it3); |
| 1746 |
lst.erase(it2); |
| 1747 |
} |
| 1748 |
|
| 1749 |
void TOIBREP::ajouter_liste(LISTE_FM& lst,FEM_TETRA* tet) |
| 1750 |
{ |
| 1751 |
lst.push_back(tet); |
| 1752 |
} |
| 1753 |
|
| 1754 |
void TOIBREP::supprimer_liste(LISTE_FM& lst,FEM_TETRA* tet) |
| 1755 |
{ |
| 1756 |
LISTE_FM::iterator it=find(lst.begin(),lst.end(),tet); |
| 1757 |
if (it!=lst.end()) lst.erase(it); |
| 1758 |
} |
| 1759 |
|
| 1760 |
|
| 1761 |
|
| 1762 |
double TOIBREP::resoudgradT(FEM_TETRA* tet,int *signe) |
| 1763 |
{ |
| 1764 |
double j[9]; |
| 1765 |
double N[12]; |
| 1766 |
double jN[12]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; |
| 1767 |
double T[4]; |
| 1768 |
double uv[2]={0.25,0.25}; |
| 1769 |
|
| 1770 |
|
| 1771 |
tet->get_inverse_jacob(j,uv); |
| 1772 |
for (int i=0;i<3;i++) |
| 1773 |
for (int k=0;k<4;k++) |
| 1774 |
N[i*4+k]=tet->get_fonction_derive_interpolation(k+1,i+1,uv); |
| 1775 |
int premier=0; |
| 1776 |
double tmin=1e300; |
| 1777 |
double tmax=-1e300; |
| 1778 |
for (int i=0;i<4;i++) |
| 1779 |
{ |
| 1780 |
if (i==tet->get_numero()) continue; |
| 1781 |
T[i]=tet->get_fem_noeud(i)->get_solution(); |
| 1782 |
if (fabs(T[i])>0.000001) |
| 1783 |
{ |
| 1784 |
if (premier==0) |
| 1785 |
if (T[i]>0) (*signe)=1; else (*signe)=-1; |
| 1786 |
else if (T[i]*(*signe)<0) (*signe)=0; |
| 1787 |
} |
| 1788 |
T[i]=fabs(T[i]); |
| 1789 |
if (tet->get_numero()!=i) |
| 1790 |
{ |
| 1791 |
if (T[i]<tmin) tmin=T[i]; |
| 1792 |
if (T[i]>tmax) tmax=T[i]; |
| 1793 |
} |
| 1794 |
premier=1; |
| 1795 |
} |
| 1796 |
for (int i=0;i<3;i++) |
| 1797 |
for (int k=0;k<4;k++) |
| 1798 |
for (int l=0;l<3;l++) |
| 1799 |
jN[i*4+k]=jN[i*4+k]+j[i*3+l]*N[l*4+k]; |
| 1800 |
double a=0.,b=0.,c=-1.; |
| 1801 |
for (int i=0;i<3;i++) |
| 1802 |
{ |
| 1803 |
double coef=0.; |
| 1804 |
double coefinc=0.; |
| 1805 |
for (int l=0;l<4;l++) |
| 1806 |
{ |
| 1807 |
if (tet->get_numero()!=l) coef=coef+jN[i*4+l]*T[l]; |
| 1808 |
else coefinc=coefinc+jN[i*4+l]; |
| 1809 |
} |
| 1810 |
c=c+coef*coef; |
| 1811 |
a=a+coefinc*coefinc; |
| 1812 |
b=b+2*coef*coefinc; |
| 1813 |
} |
| 1814 |
/*if (*signe==0) |
| 1815 |
cout << "attention " <<endl;*/ |
| 1816 |
double det=b*b-4.*a*c; |
| 1817 |
if (det<0.) det=0.; else det=sqrt(det); |
| 1818 |
double sol1=(-b-det)/2./a; |
| 1819 |
double sol2=(-b+det)/2./a; |
| 1820 |
double sol=sol1; |
| 1821 |
if (sol2>sol1) sol=sol2; |
| 1822 |
if (sol<tmin*0.99) |
| 1823 |
sol=0.; |
| 1824 |
sol=sol*(*signe); |
| 1825 |
return sol; |
| 1826 |
} |
