MAGiC  V5.0
Mailleurs Automatiques de Géometries intégrés à la Cao
mg_geodesic.cpp
Aller à la documentation de ce fichier.
1 //####//------------------------------------------------------------
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 //####// mg_geodesic.cpp
15 //####//
16 //####//------------------------------------------------------------
17 //####//------------------------------------------------------------
18 //####// COPYRIGHT 2000-2024
19 //####// jeu 13 jun 2024 11:58:56 EDT
20 //####//------------------------------------------------------------
21 //####//------------------------------------------------------------
22 
23 
24 #include "gestionversion.h"
25 #include "mg_geodesic.h"
26 #include "mg_fast_marching2D.h"
27 #include "mg_gestionnaire.h"
28 #include "fem_maillage.h"
29 #include "fem_triangle3.h"
30 #include "fem_solution.h"
31 #include "tpl_map_entite.h"
33 #include "ot_mathematique.h"
34 #include <math.h>
35 #include "pars_argument.h"
36 #include "parse.h"
37 #include "fem_maillage_outils.h"
38 #include "ot_boite_3d.h"
39 
40 
41 namespace MAGIC
42 {
43 #define INFINI 1e300
44 
45 
46 MG_GEODESIC::MG_GEODESIC(char* chaine,MG_GESTIONNAIRE *g,FEM_MAILLAGE *f,int nbpt):MAGIC_AFFICHE(),fem(f),gest(g)
47 {
48  char chaine2[500];
49  sprintf(chaine2,"%s0_fm.sol",chaine);
50  sol=gest->get_fem_solution(fem,nbpt+2,chaine2,(char*)"Geodesic");
51  if (sol==NULL)
52  {
53  sol=new FEM_SOLUTION(fem,nbpt+2,chaine2,1,"Geodesic");
54  for (int i=0;i<nbpt;i++)
55  {
56  char numlevel[500];
57  sprintf(numlevel,"L%d",i+1);
58  sol->change_legende(i,numlevel);
59  }
60  sol->change_legende(nbpt,(char*)"Li+Lj");
61  sol->change_legende(nbpt+1,(char*)"Li-Lj");
63  }
64 /* sprintf(chaine2,"%s1_fm.sol",chaine);
65  solgrad=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
66  if (solgrad==NULL)
67  {
68  solgrad=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::VECTEUR);
69  solgrad->change_legende(0,(char*)"LA'");
70  solgrad->change_legende(1,(char*)"LB'");
71  gest->ajouter_fem_solution(solgrad);
72  }
73  sprintf(chaine2,"%s2_fm.sol",chaine);
74  solgradnoeud=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
75  if (solgradnoeud==NULL)
76  {
77  solgradnoeud=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::VECTEUR);
78  solgradnoeud->change_legende(0,(char*)"LA'_ext");
79  solgradnoeud->change_legende(1,(char*)"LB'_ext");
80  gest->ajouter_fem_solution(solgradnoeud);
81  }
82  sprintf(chaine2,"%s3_fm.sol",chaine);
83  solgradgrad=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
84  if (solgradgrad==NULL)
85  {
86  solgradgrad=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::TENSEUR);
87  solgradgrad->change_legende(0,(char*)"LA''");
88  solgradgrad->change_legende(1,(char*)"LB''");
89  gest->ajouter_fem_solution(solgradgrad);
90  }
91  sprintf(chaine2,"%s4_fm.sol",chaine);
92  courburemax=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
93  if (courburemax==NULL)
94  {
95  courburemax=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::SCALAIRE);
96  courburemax->change_legende(0,(char*)"CAmax");
97  courburemax->change_legende(1,(char*)"CBmax");
98  gest->ajouter_fem_solution(courburemax);
99  }*/
100 
101 }
102 
103 
104 
105 MG_GEODESIC::MG_GEODESIC(MG_GEODESIC &mdd):sol(mdd.sol),fem(mdd.fem)
106 {
107 }
108 
110 {
111 }
112 
113 
114 
115 
116 bool MG_GEODESIC::init_noeud(char* entite)
117 {
118 PARS_ARGUMENT param[100];
119 PARSE parse;
120 parse.decode(entite,"&",param);
121 int nb=param[0].argument.size();
123 for (int i=0;i<nb;i++)
124  {
125  int num=atoi(param[0].argument[i].c_str());
126  MG_ELEMENT_TOPOLOGIQUE *ele=NULL;
127  ele=geo->get_mg_sommetid(num);
128  if (ele!=NULL)
129  {
130  int nbele=ele->get_lien_fem_maillage()->get_nb();
131  for (int j=0;j<nbele;j++)
132  {
133  FEM_NOEUD* no=(FEM_NOEUD*)ele->get_lien_fem_maillage()->get(j);
134  if (fem->get_fem_noeudid(no->get_id())==no) lst.push_back(no);
135  }
136  }
137 
138  }
139 initdouble=false;
140 int nbnoeudinit=lst.size();
141 if (nbnoeudinit!=2) return false;
142 return true;
143 }
144 bool MG_GEODESIC::init_coord(char* entite)
145 {
146 PARS_ARGUMENT param[100];
147 PARSE parse;
148 parse.decode(entite,"@,@,@,@,@,@",param);
149 lstdouble.push_back(atof(param[0].argument[0].c_str()));
150 lstdouble.push_back(atof(param[1].argument[0].c_str()));
151 lstdouble.push_back(atof(param[2].argument[0].c_str()));
152 lstdouble.push_back(atof(param[3].argument[0].c_str()));
153 lstdouble.push_back(atof(param[4].argument[0].c_str()));
154 lstdouble.push_back(atof(param[5].argument[0].c_str()));
155 initdouble=true;
156 return true;
157 }
158 
159 
160 
161 
162 double MG_GEODESIC::calcul_fm(int num,double xdep,double ydep,double zdep,std::vector<double> &lstpoints)
163 {
164 MG_FAST_MARCHING2D fm(sol,num);
166 fm.init();
167 if (initdouble==false)
168  {
169  std::vector<FEM_NOEUD*> lst2;
170  lst2.push_back(lst[num]);
171  fm.init_noeud(lst2);
172  }
173 else
174  fm.init_coord(xdep,ydep,zdep);
175 fm.propage(1.0);
176 fm.finalise();
177 sol->active_solution(num);
178 double valA;
179 if (initdouble==false) valA=lst[(num+1)%2]->get_solution();
180 else
181  {
182  int nbval=lstpoints.size()/3;
183  for (int i=0;i<nbval;i++)
184  {
185  valA=fm.get_valeur(lstpoints[3*i+0],lstpoints[3*i+1],lstpoints[3*i+2]);
186  lstpoints.push_back(valA);
187  }
188  }
189 return valA;
190 }
191 
192 void MG_GEODESIC::calcul_reseau(std::vector<double> lstpoint,double *matdist)
193 {
194 initdouble=true;
195 int nb=lstpoint.size()/3;
196 for (int i=0;i<nb;i++)
197  for (int j=0;j<nb;j++)
198  matdist[i*nb+j]=0;
199 for (int i=0;i<nb;i++)
200  {
201  std::vector<double> lst=lstpoint;
202  lst.erase(lst.begin()+3*i,lst.begin()+3*i+3);
203  calcul_fm(i,lstpoint[3*i],lstpoint[3*i+1],lstpoint[3*i+2],lst);
204  int ii=0;
205  for (int j=0;j<nb;j++)
206  if (j!=i)
207  {
208  double val=lst[3*nb-3+ii];
209  matdist[i*nb+j]=val;
210  //matdist[j*nb+i]=matdist[j*nb+i]+lst[3*nb-3+ii];
211  ii++;
212  }
213  }
214 for (int i=0;i<nb;i++)
215  for (int j=i+1;j<nb;j++)
216  if (i!=j)
217  {
218  double vala,valb;
219  calcul_milieux(i,j,nb,nb+1,vala,valb);
220  matdist[i*nb+j]=matdist[i*nb+j]+matdist[j*nb+i]-vala-valb;
221  matdist[j*nb+i]=matdist[i*nb+j];
222  }
223 }
224 
225 double MG_GEODESIC::calcul(double *tab)
226 {
227 std::vector<double> lst;
228 double xdep,ydep,zdep;
229 if (initdouble==true)
230  {
231  lst.push_back(lstdouble[3]);
232  lst.push_back(lstdouble[4]),
233  lst.push_back(lstdouble[5]);
234  xdep=lstdouble[0];
235  ydep=lstdouble[1];
236  zdep=lstdouble[2];
237  }
238 double valA=calcul_fm(0,xdep,ydep,zdep,lst);
239 if (initdouble==true)
240  {
241  lst.clear();
242  lst.push_back(lstdouble[0]);
243  lst.push_back(lstdouble[1]),
244  lst.push_back(lstdouble[2]);
245  xdep=lstdouble[3];
246  ydep=lstdouble[4];
247  zdep=lstdouble[5];
248  }
249 double valB=calcul_fm(1,xdep,ydep,zdep,lst);
250 double vala,valb;
251 calcul_milieux(0,1,2,3,vala,valb);
252 double longueur=valA+valB-vala-valb;
253 if (tab!=NULL)
254  {
255  tab[0]=valA;
256  tab[1]=valB;
257  tab[2]=longueur;
258  }
259 return longueur;
260 }
261 
262 void MG_GEODESIC::calcul_milieux(int num1,int num2,int num3, int num4,double &vala,double &valb)
263 {
267 sol->active_solution(num4);
268 
269 int ligne=0;
270 LISTE_FEM_NOEUD::iterator itn;
271 for (FEM_NOEUD* no=fem->get_premier_noeud(itn);no!=NULL;no=fem->get_suivant_noeud(itn))
272  {
273  double val=sol->lire(ligne,num3);
274  no->change_solution(val,1);
275  val=sol->lire(ligne,num1);
276  no->change_solution(val,2);
277  val=sol->lire(ligne,num2);
278  no->change_solution(val,3);
279  ligne++;
280  }
281 double limit=0.;
282 double minApB=INFINI;
283 LISTE_FEM_ELEMENT2::iterator it;
284 double xm,ym,zm;
285 for (FEM_ELEMENT2 *ele=fem->get_premier_element2(it);ele!=NULL;ele=fem->get_suivant_element2(it))
286  {
287  double val1=ele->get_fem_noeud(0)->get_solution();
288  double val2=ele->get_fem_noeud(1)->get_solution();
289  double val3=ele->get_fem_noeud(2)->get_solution();
290  if ((val1-limit)*(val2-limit)<1e-12)
291  {
292  double t=(limit-val1)/(val2-val1);
293  double val1apb=ele->get_fem_noeud(0)->get_solution(1);
294  double val2apb=ele->get_fem_noeud(1)->get_solution(1);
295  double valApB=val1apb+t*(val2apb-val1apb);
296  if (valApB<minApB)
297  {
298  minApB=valApB;
299  double val1a=ele->get_fem_noeud(0)->get_solution(2);
300  double val2a=ele->get_fem_noeud(1)->get_solution(2);
301  vala=val1a+t*(val2a-val1a);
302  double val1b=ele->get_fem_noeud(0)->get_solution(3);
303  double val2b=ele->get_fem_noeud(1)->get_solution(3);
304  valb=val1b+t*(val2b-val1b);
305  xm=ele->get_fem_noeud(0)->get_x()+t*(ele->get_fem_noeud(1)->get_x()-ele->get_fem_noeud(0)->get_x());
306  ym=ele->get_fem_noeud(0)->get_y()+t*(ele->get_fem_noeud(1)->get_y()-ele->get_fem_noeud(0)->get_y());
307  zm=ele->get_fem_noeud(0)->get_z()+t*(ele->get_fem_noeud(1)->get_z()-ele->get_fem_noeud(0)->get_z());
308  }
309  }
310  if ((val2-limit)*(val3-limit)<1e-12)
311  {
312  double t=(limit-val2)/(val3-val2);
313  double val2apb=ele->get_fem_noeud(1)->get_solution(1);
314  double val3apb=ele->get_fem_noeud(2)->get_solution(1);
315  double valApB=val2apb+t*(val3apb-val2apb);
316  if (valApB<minApB)
317  {
318  minApB=valApB;
319  double val2a=ele->get_fem_noeud(1)->get_solution(2);
320  double val3a=ele->get_fem_noeud(2)->get_solution(2);
321  vala=val2a+t*(val3a-val2a);
322  double val2b=ele->get_fem_noeud(1)->get_solution(3);
323  double val3b=ele->get_fem_noeud(2)->get_solution(3);
324  valb=val2b+t*(val3b-val2b);
325  xm=ele->get_fem_noeud(1)->get_x()+t*(ele->get_fem_noeud(2)->get_x()-ele->get_fem_noeud(1)->get_x());
326  ym=ele->get_fem_noeud(1)->get_y()+t*(ele->get_fem_noeud(2)->get_y()-ele->get_fem_noeud(1)->get_y());
327  zm=ele->get_fem_noeud(1)->get_z()+t*(ele->get_fem_noeud(2)->get_z()-ele->get_fem_noeud(1)->get_z());
328  }
329  }
330  if ((val1-limit)*(val3-limit)<1e-12)
331  {
332  double t=(limit-val1)/(val3-val1);
333  double val1apb=ele->get_fem_noeud(0)->get_solution(1);
334  double val3apb=ele->get_fem_noeud(2)->get_solution(1);
335  double valApB=val1apb+t*(val3apb-val1apb);
336  if (valApB<minApB)
337  {
338  minApB=valApB;
339  double val1a=ele->get_fem_noeud(0)->get_solution(2);
340  double val3a=ele->get_fem_noeud(2)->get_solution(2);
341  vala=val1a+t*(val3a-val1a);
342  double val1b=ele->get_fem_noeud(0)->get_solution(3);
343  double val3b=ele->get_fem_noeud(2)->get_solution(3);
344  valb=val1b+t*(val3b-val1b);
345  xm=ele->get_fem_noeud(0)->get_x()+t*(ele->get_fem_noeud(2)->get_x()-ele->get_fem_noeud(0)->get_x());
346  ym=ele->get_fem_noeud(0)->get_y()+t*(ele->get_fem_noeud(2)->get_y()-ele->get_fem_noeud(0)->get_y());
347  zm=ele->get_fem_noeud(0)->get_z()+t*(ele->get_fem_noeud(2)->get_z()-ele->get_fem_noeud(0)->get_z());
348  }
349  }
350 
351  }
352 char mess[500];
353 sprintf(mess," Le point milieu x=%lf, y=%lf, z=%lf la demi longueur de la LS1=%f la demi longueur de la LS2=%f ",xm,ym,zm,vala,valb) ;
354 affiche(mess);
355 }
356 
357 };
MAGIC::MG_GEODESIC::lst
std::vector< FEM_NOEUD * > lst
Definition: mg_geodesic.h:61
FEM_SOLUTION
Definition: fem_solution.h:40
MAGIC::MG_GEODESIC
Definition: mg_geodesic.h:42
gestionversion.h
MAGIC::MG_FAST_MARCHING2D::propage
virtual void propage(double vitesse=1.0)
Definition: mg_fast_marching2D.cpp:148
MAGIC::MG_GEODESIC::fem
FEM_MAILLAGE * fem
Definition: mg_geodesic.h:59
MAGIC_AFFICHE::active_affichage
virtual void active_affichage(fonction_affiche *fonc)
Definition: magic_affiche.cpp:49
mg_fast_marching2D.h
MAGIC::MG_GEODESIC::calcul
virtual double calcul(double *tab=NULL)
Definition: mg_geodesic.cpp:225
fem_solution.h
mg_gestionnaire.h
MAGIC::MG_GEODESIC::init_noeud
virtual bool init_noeud(char *entite)
Definition: mg_geodesic.cpp:116
MAGIC::MG_GEODESIC::MG_GEODESIC
MG_GEODESIC(char *fichiermagic, MG_GESTIONNAIRE *g, FEM_MAILLAGE *f, int nbpt=2)
Definition: mg_geodesic.cpp:46
MAGIC::MG_GEODESIC::gest
MG_GESTIONNAIRE * gest
Definition: mg_geodesic.h:64
MG_IDENTIFICATEUR::get_id
unsigned long get_id()
Definition: mg_identificateur.cpp:53
fem_maillage.h
MAGIC::MG_FAST_MARCHING::init_noeud
virtual void init_noeud(std::vector< FEM_NOEUD * > &lst)
Definition: mg_fast_marching.cpp:75
FEM_MAILLAGE::get_premier_noeud
FEM_NOEUD * get_premier_noeud(LISTE_FEM_NOEUD::iterator &it)
Definition: fem_maillage.cpp:174
FEM_MAILLAGE::get_mg_geometrie
MG_GEOMETRIE * get_mg_geometrie(void)
Definition: fem_maillage.cpp:88
FEM_MAILLAGE_OUTILS::operation_champs_solution
void operation_champs_solution(class FEM_SOLUTION *sol1, int numchamp1, FEM_SOLUTION *sol2, int numchamp2, FEM_SOLUTION *sol3, int numchamp3, int operation)
Definition: fem_maillage_outils.cpp:622
MAGIC_AFFICHE::affiche
virtual void affiche(char *mess)
Definition: magic_affiche.cpp:43
MG_GESTIONNAIRE
Definition: mg_gestionnaire.h:57
MG_GESTIONNAIRE::ajouter_fem_solution
int ajouter_fem_solution(FEM_SOLUTION *mgsol)
Definition: mg_gestionnaire.cpp:902
MAGIC_AFFICHE::affiche2
fonction_affiche * affiche2
Definition: magic_affiche.h:41
FEM_MAILLAGE::get_fem_noeudid
FEM_NOEUD * get_fem_noeudid(unsigned long num)
Definition: fem_maillage.cpp:150
FEM_ELEMENT2
Definition: fem_element2.h:34
MAGIC::MG_FAST_MARCHING2D
Definition: mg_fast_marching2D.h:42
parse.h
MG_GESTIONNAIRE::get_fem_solution
FEM_SOLUTION * get_fem_solution(unsigned int num)
Definition: mg_gestionnaire.cpp:930
FEM_MAILLAGE::get_premier_element2
FEM_ELEMENT2 * get_premier_element2(LISTE_FEM_ELEMENT2::iterator &it)
Definition: fem_maillage.cpp:561
MG_ELEMENT_TOPOLOGIQUE
Definition: mg_element_topologique.h:51
MAGIC::MG_GEODESIC::init_coord
virtual bool init_coord(char *entite)
Definition: mg_geodesic.cpp:144
MAGIC
Definition: mg_fast_marching.cpp:40
f
double f(double x, long nb, double *xfonc, double *fonc, double eng, double eni, double lambda, double nor, double *fonc2)
Definition: fct_generateur_calibrage.cpp:96
FEM_SOLUTION::active_solution
void active_solution(int num)
Definition: fem_solution.cpp:490
FEM_MAILLAGE
Definition: fem_maillage.h:66
mg_geodesic.h
MAGIC::MG_FAST_MARCHING2D::init_coord
virtual void init_coord(double x, double y, double z)
Definition: mg_fast_marching2D.cpp:65
tpl_liste_unique_entite.h
INFINI
#define INFINI
Definition: mg_geodesic.cpp:43
MAGIC_AFFICHE
Definition: magic_affiche.h:30
PARSE::decode
void decode(char *code, std::string masque, class PARS_ARGUMENT *arg)
Definition: parse.cpp:71
tpl_map_entite.h
MAGIC::OPERATION_FEM_SOLUTION::SOUSTRACTION
@ SOUSTRACTION
Definition: mg_definition.h:128
PARSE
Definition: parse.h:32
TPL_LISTE_ENTITE::get_nb
virtual int get_nb(void)
Definition: tpl_liste_entite.h:67
FEM_MAILLAGE::get_suivant_element2
FEM_ELEMENT2 * get_suivant_element2(LISTE_FEM_ELEMENT2::iterator &it)
Definition: fem_maillage.cpp:569
FEM_SOLUTION::lire
double lire(int i, int j, int coord=0, int num_no=0)
Definition: fem_solution.cpp:398
TPL_LISTE_ENTITE::get
virtual X get(int num)
Definition: tpl_liste_entite.h:72
ot_mathematique.h
FEM_NOEUD
Definition: fem_noeud.h:35
MAGIC::MG_FAST_MARCHING2D::get_valeur
virtual double get_valeur(double x, double y, double z)
Definition: mg_fast_marching2D.cpp:116
MAGIC::MG_GEODESIC::calcul_reseau
void calcul_reseau(std::vector< double > lstpoint, double *matdist)
Definition: mg_geodesic.cpp:192
PARS_ARGUMENT::argument
std::vector< std::string > argument
Definition: pars_argument.h:44
MAGIC::MG_GEODESIC::initdouble
bool initdouble
Definition: mg_geodesic.h:63
FEM_MAILLAGE_OUTILS
Definition: fem_maillage_outils.h:56
MAGIC::MG_FAST_MARCHING::finalise
virtual void finalise(void)
Definition: mg_fast_marching.cpp:240
MAGIC::MG_GEODESIC::~MG_GEODESIC
virtual ~MG_GEODESIC()
Definition: mg_geodesic.cpp:109
ot_boite_3d.h
PARS_ARGUMENT
Definition: pars_argument.h:37
MAGIC_AFFICHE::affichageactif
int affichageactif
Definition: magic_affiche.h:42
MAGIC::MG_GEODESIC::calcul_milieux
void calcul_milieux(int num1, int num2, int num3, int num4, double &vala, double &valb)
Definition: mg_geodesic.cpp:262
MAGIC::MG_FAST_MARCHING::init
virtual void init(void)
Definition: mg_fast_marching.cpp:61
fem_triangle3.h
MG_ELEMENT_TOPOLOGIQUE::get_lien_fem_maillage
virtual TPL_LISTE_ENTITE< FEM_ELEMENT_MAILLAGE * > * get_lien_fem_maillage(void)
Definition: mg_element_topologique.cpp:67
MG_GEOMETRIE
Definition: mg_geometrie.h:84
MAGIC::MG_GEODESIC::sol
FEM_SOLUTION * sol
Definition: mg_geodesic.h:60
FEM_MAILLAGE::get_suivant_noeud
FEM_NOEUD * get_suivant_noeud(LISTE_FEM_NOEUD::iterator &it)
Definition: fem_maillage.cpp:182
MAGIC::MG_GEODESIC::lstdouble
std::vector< double > lstdouble
Definition: mg_geodesic.h:62
MAGIC::OPERATION_FEM_SOLUTION::ADDITION
@ ADDITION
Definition: mg_definition.h:128
pars_argument.h
MAGIC::MG_GEODESIC::calcul_fm
double calcul_fm(int num, double xdep, double ydep, double zdep, std::vector< double > &lstpoints)
Definition: mg_geodesic.cpp:162
MG_GEOMETRIE::get_mg_sommetid
MG_SOMMET * get_mg_sommetid(unsigned long num)
Definition: mg_geometrie.cpp:513
FEM_SOLUTION::change_legende
void change_legende(int num, std::string val)
Definition: fem_solution.cpp:457
fem_maillage_outils.h