diamesh/src/m3d_move.cpp

/* BOUGE DE POINTS SAUCE P.L.GEORGE Article INRIA STRUCOME 1994 */
#include <stdio.h>
#include <math.h>
#include "m3d_struct.h"
#include "m3d_const.h"
#include "m3d_hotes.h"
#include "m3d_erreur.h"
#include "prototype.h"
extern GEST_MEM *gest ;
extern int debug ;
extern float crit_min ;
int m3d_move(int nump)
{
    NOEUD *n1 ;
    CONNEC *connec ;
    TETRAEDRE *tab_tetra[1000] ;
    int nb_ele, i, j, k, *numele  ;
    float  tab_crit[1000], vec_opt[3], cmin ;
    float tab_opt[3000] ;
    int tabele[3000], sens ;
    float gamma, vn[3], xi, yi, zi, vab[3], vac[3], vbc[3], pvec[3], hauteur, vap[3] ;
    float dab, dac, dbc, perimetre, norme, coeff, alpha  ;
    float delta[3] ;
    float vec_save[3], bmin, bmax, *coord, eps, crit_save ;
    float calpha, calpha_moins_eps, vresu[3] ;

    /* coordonnees et connectivites */
    coord = gest->coord ;
    numele = gest->numele ;

    if (nump<gest->nb_init) return(VRAI) ; /* on ne bouge pas un noeud sur la peau */
    ADRESSE(noeud,nump,n1)
    if (n1 == NULL)
    {
        m3d_erreur(ERR_SYST) ;
        return(FAUX) ;
    }
    if (n1->mark == KILLED) return(VRAI) ;

    /* sauvegarde des coordonnees du noeud */
    vec_save[0] = coord[x(nump)] ;
    vec_save[1] = coord[y(nump)] ;
    vec_save[2] = coord[z(nump)] ;

    /* tetraedres connectes au noeud */
    nb_ele = 0 ;
    connec = n1->lis_con ;
    while (connec != NULL)
    {
        if ((connec->tetra)->mark != KILLED)
        {
            tab_tetra[nb_ele] = connec->tetra ;
            /* determination des noeuds de la face de la boule */
            k = 0 ;
            for (j=0;j<4;j++)
                if (gest->numele[4*(connec->tetra->num)+ j] != n1->num)
                {
                    tabele[3*nb_ele + k] = gest->numele[4*(connec->tetra->num)+ j] ;
                    k++ ;
                }
            if (k!=3) return(FAUX) ;
            /* determination du point optimal de la face nb_ele */
            xi = (coord[x(tabele[3*nb_ele])] + coord[x(tabele[3*nb_ele+1])] + coord[x(tabele[3*nb_ele+2])])/3. ;
            yi = (coord[y(tabele[3*nb_ele])] + coord[y(tabele[3*nb_ele+1])] + coord[y(tabele[3*nb_ele+2])])/3. ;
            zi = (coord[z(tabele[3*nb_ele])] + coord[z(tabele[3*nb_ele+1])] + coord[z(tabele[3*nb_ele+2])])/3. ;

            /* calcul de la normale a la face */
            vab[0] = coord[x(tabele[3*nb_ele+1])] - coord[x(tabele[3*nb_ele])] ;
            vab[1] = coord[y(tabele[3*nb_ele+1])] - coord[y(tabele[3*nb_ele])] ;
            vab[2] = coord[z(tabele[3*nb_ele+1])] - coord[z(tabele[3*nb_ele])] ;
            dab = NORME(vab) ;
            vac[0] = coord[x(tabele[3*nb_ele+2])] - coord[x(tabele[3*nb_ele])] ;
            vac[1] = coord[y(tabele[3*nb_ele+2])] - coord[y(tabele[3*nb_ele])] ;
            vac[2] = coord[z(tabele[3*nb_ele+2])] - coord[z(tabele[3*nb_ele])] ;
            dac = NORME(vac) ;
            vbc[0] = coord[x(tabele[3*nb_ele+2])] - coord[x(tabele[3*nb_ele+1])] ;
            vbc[1] = coord[y(tabele[3*nb_ele+2])] - coord[y(tabele[3*nb_ele+1])] ;
            vbc[2] = coord[z(tabele[3*nb_ele+2])] - coord[z(tabele[3*nb_ele+1])] ;
            dbc = NORME(vbc) ;

            /* calcul du critere 2 D */
            pvec[0] = PVECX(vab,vac) ;
            pvec[1] = PVECY(vab,vac) ;
            pvec[2] = PVECZ(vab,vac) ;
            perimetre = dab + dbc + dac ;
            hauteur = (perimetre/3.) * 0.8 ;

            vn[0] = PVECX(vab,vac) ;
            vn[1] = PVECY(vab,vac) ;
            vn[2] = PVECZ(vab,vac) ;
            norme = NORME(vn) ;
            if (norme<EPSILON)
            {
                //      if (debug) aff_text("normale nulle M3D_MOVE\n ") ;
                m3d_erreur(ERR_SYST) ;
                return(FAUX) ;
            }
            for (i=0;i<3;i++) vn[i] = vn[i]/norme ;
            /* calcul du vecteur AP et orientation de la face */
            vap[0] = coord[x(nump)] - coord[x(tabele[3*nb_ele])] ;
            vap[1] = coord[y(nump)] - coord[y(tabele[3*nb_ele])] ;
            vap[2] = coord[z(nump)] - coord[z(tabele[3*nb_ele])] ;
            if (PROSCA(vap,vn) < EPSILON)
            {
                for (i=0;i<3;i++) vn[i] = -1. * vn[i] ;
                i = tabele[3*nb_ele+1] ;
                tabele[3*nb_ele+1] = tabele[3*nb_ele+2] ;
                tabele[3*nb_ele+2] = i ;
            }
            /* point optimal */
            tab_opt[x(nb_ele)] = xi + hauteur * vn[0] ;
            tab_opt[y(nb_ele)] = yi + hauteur * vn[1] ;
            tab_opt[z(nb_ele)] = zi + hauteur * vn[2] ;

            nb_ele ++ ;
        }
        connec = connec->suivant ;
    }
    if (nb_ele == 0) return(VRAI) ;
    /* calcul de la qualite de chaque tetraedre */
    cmin = 1. ;
    for (i=0;i<nb_ele;i++)
    {
        /* calcul du critere */
        tab_crit[i] = m3d_e_qual2(gest->coord,tabele[3*i],tabele[3*i+1],tabele[3*i+2],nump) ;
        if (tab_crit[i]<EPSILON)
        {
            //  if (debug) aff_text("critere nul M3D_MOVE\n ") ;
            m3d_erreur(ERR_SYST) ;
            return(FAUX) ;
        }
        if (tab_crit[i] < cmin) cmin = tab_crit[i] ;
    }
    /* calcul des coordonnees du point ideal */
    gamma = 0. ;
    for (i=0;i<nb_ele;i++)
    {
        gamma = gamma + 1./(tab_crit[i]*tab_crit[i]) ;/* somme des carres des inverses */
    }
    /* calcul du point optimal pour la boule */
    gamma = 1./gamma ;
    vec_opt[0] = 0. ;
    vec_opt[1] = 0. ;
    vec_opt[2] = 0. ;

    for (i=0;i<nb_ele;i++)
    {
        coeff = gamma/(tab_crit[i]*tab_crit[i]) ;
        vec_opt[0] = vec_opt[0] + coeff * tab_opt[x(i)] ;
        vec_opt[1] = vec_opt[1] + coeff * tab_opt[y(i)] ;
        vec_opt[2] = vec_opt[2] + coeff * tab_opt[z(i)] ;
    }

    delta[0] = vec_opt[0] - coord[x(nump)] ;
    delta[1] = vec_opt[1] - coord[y(nump)] ;
    delta[2] = vec_opt[2] - coord[z(nump)] ;
    /* sauvegarde du critere mini */
    crit_save = cmin ;

    /* initialisation des bornes du parametre */
    bmin = 0. ;
    bmax = 1. ;
    sens = 1 ;/* augmnentation */

    vresu[0] = vec_save[0] ;
    vresu[1] = vec_save[1] ;
    vresu[2] = vec_save[2] ;

    for (i=0;i<5;i++)
    {
        alpha = 0.5 * (bmin + bmax) ;
        coord[x(nump)] = vec_save[0] + alpha * delta[0] ;
        coord[y(nump)] = vec_save[1] + alpha * delta[1] ;
        coord[z(nump)] = vec_save[2] + alpha * delta[2] ;

        /* evaluation de la qualite pour la valeur alpha */
        /* on evalue la qualite cmin pour alpha */
        calpha = 1. ;
        for (j=0;j<nb_ele;j++)
        {
            /* calcul du critere */
            tab_crit[j] = m3d_e_qual2(gest->coord,tabele[3*j],tabele[3*j+1],tabele[3*j+2],nump) ;
            if (tab_crit[j] < calpha) calpha = tab_crit[j] ;
        }
        /* la valeur du critere pour alpha est superieure, il faut la stocker */
        if (calpha > crit_save)
        {
            vresu[0] = coord[x(nump)] ;
            vresu[1] = coord[y(nump)] ;
            vresu[2] = coord[z(nump)] ;
            crit_save = calpha ;
        }

        /* on evalue la qualite c1 pour alpha - eps */
        eps = (bmax-bmin)/50. ;
        coord[x(nump)] = vec_save[0] + (alpha - eps) * delta[0] ;
        coord[y(nump)] = vec_save[1] + (alpha - eps) * delta[1] ;
        coord[z(nump)] = vec_save[2] + (alpha - eps) * delta[2] ;
        calpha_moins_eps = 1. ;
        for (j=0;j<nb_ele;j++)
        {
            /* calcul du critere */
            tab_crit[j] = m3d_e_qual2(gest->coord,tabele[3*j],tabele[3*j+1],tabele[3*j+2],nump) ;
            if (tab_crit[j] < calpha_moins_eps) calpha_moins_eps = tab_crit[j] ;
        }

        coord[x(nump)] = vec_save[0] + alpha * delta[0] ;
        coord[y(nump)] = vec_save[1] + alpha * delta[1] ;
        coord[z(nump)] = vec_save[2] + alpha * delta[2] ;

        if (calpha_moins_eps < calpha)
        {
            bmin = alpha ;
            /* bmax est inchange */
        }
        else
        {
            /* bmin est inchange */
            bmax = alpha ;
        }
    }
    coord[x(nump)] = vresu[0] ;
    coord[y(nump)] = vresu[1] ;
    coord[z(nump)] = vresu[2] ;

    return(VRAI) ;
}


Revision:	283
Committed:	Tue Sep 13 21:11:20 2011 UTC (13 years, 8 months ago) by francois
File size:	8119 byte(s)
Log Message:	structure de l'écriture
#	User	Rev	Content
1	francois	283	/* BOUGE DE POINTS SAUCE P.L.GEORGE Article INRIA STRUCOME 1994 */
2			#include <stdio.h>
3			#include <math.h>
4			#include "m3d_struct.h"
5			#include "m3d_const.h"
6			#include "m3d_hotes.h"
7			#include "m3d_erreur.h"
8			#include "prototype.h"
9			extern GEST_MEM *gest ;
10			extern int debug ;
11			extern float crit_min ;
12			int m3d_move(int nump)
13			{
14			NOEUD *n1 ;
15			CONNEC *connec ;
16			TETRAEDRE *tab_tetra[1000] ;
17			int nb_ele, i, j, k, *numele ;
18			float tab_crit[1000], vec_opt[3], cmin ;
19			float tab_opt[3000] ;
20			int tabele[3000], sens ;
21			float gamma, vn[3], xi, yi, zi, vab[3], vac[3], vbc[3], pvec[3], hauteur, vap[3] ;
22			float dab, dac, dbc, perimetre, norme, coeff, alpha ;
23			float delta[3] ;
24			float vec_save[3], bmin, bmax, *coord, eps, crit_save ;
25			float calpha, calpha_moins_eps, vresu[3] ;
26
27			/* coordonnees et connectivites */
28			coord = gest->coord ;
29			numele = gest->numele ;
30
31			if (nump<gest->nb_init) return(VRAI) ; /* on ne bouge pas un noeud sur la peau */
32			ADRESSE(noeud,nump,n1)
33			if (n1 == NULL)
34			{
35			m3d_erreur(ERR_SYST) ;
36			return(FAUX) ;
37			}
38			if (n1->mark == KILLED) return(VRAI) ;
39
40			/* sauvegarde des coordonnees du noeud */
41			vec_save[0] = coord[x(nump)] ;
42			vec_save[1] = coord[y(nump)] ;
43			vec_save[2] = coord[z(nump)] ;
44
45			/* tetraedres connectes au noeud */
46			nb_ele = 0 ;
47			connec = n1->lis_con ;
48			while (connec != NULL)
49			{
50			if ((connec->tetra)->mark != KILLED)
51			{
52			tab_tetra[nb_ele] = connec->tetra ;
53			/* determination des noeuds de la face de la boule */
54			k = 0 ;
55			for (j=0;j<4;j++)
56			if (gest->numele[4*(connec->tetra->num)+ j] != n1->num)
57			{
58			tabele[3nb_ele + k] = gest->numele[4(connec->tetra->num)+ j] ;
59			k++ ;
60			}
61			if (k!=3) return(FAUX) ;
62			/* determination du point optimal de la face nb_ele */
63			xi = (coord[x(tabele[3nb_ele])] + coord[x(tabele[3nb_ele+1])] + coord[x(tabele[3*nb_ele+2])])/3. ;
64			yi = (coord[y(tabele[3nb_ele])] + coord[y(tabele[3nb_ele+1])] + coord[y(tabele[3*nb_ele+2])])/3. ;
65			zi = (coord[z(tabele[3nb_ele])] + coord[z(tabele[3nb_ele+1])] + coord[z(tabele[3*nb_ele+2])])/3. ;
66
67			/* calcul de la normale a la face */
68			vab[0] = coord[x(tabele[3nb_ele+1])] - coord[x(tabele[3nb_ele])] ;
69			vab[1] = coord[y(tabele[3nb_ele+1])] - coord[y(tabele[3nb_ele])] ;
70			vab[2] = coord[z(tabele[3nb_ele+1])] - coord[z(tabele[3nb_ele])] ;
71			dab = NORME(vab) ;
72			vac[0] = coord[x(tabele[3nb_ele+2])] - coord[x(tabele[3nb_ele])] ;
73			vac[1] = coord[y(tabele[3nb_ele+2])] - coord[y(tabele[3nb_ele])] ;
74			vac[2] = coord[z(tabele[3nb_ele+2])] - coord[z(tabele[3nb_ele])] ;
75			dac = NORME(vac) ;
76			vbc[0] = coord[x(tabele[3nb_ele+2])] - coord[x(tabele[3nb_ele+1])] ;
77			vbc[1] = coord[y(tabele[3nb_ele+2])] - coord[y(tabele[3nb_ele+1])] ;
78			vbc[2] = coord[z(tabele[3nb_ele+2])] - coord[z(tabele[3nb_ele+1])] ;
79			dbc = NORME(vbc) ;
80
81			/* calcul du critere 2 D */
82			pvec[0] = PVECX(vab,vac) ;
83			pvec[1] = PVECY(vab,vac) ;
84			pvec[2] = PVECZ(vab,vac) ;
85			perimetre = dab + dbc + dac ;
86			hauteur = (perimetre/3.) * 0.8 ;
87
88			vn[0] = PVECX(vab,vac) ;
89			vn[1] = PVECY(vab,vac) ;
90			vn[2] = PVECZ(vab,vac) ;
91			norme = NORME(vn) ;
92			if (norme<EPSILON)
93			{
94			// if (debug) aff_text("normale nulle M3D_MOVE\n ") ;
95			m3d_erreur(ERR_SYST) ;
96			return(FAUX) ;
97			}
98			for (i=0;i<3;i++) vn[i] = vn[i]/norme ;
99			/* calcul du vecteur AP et orientation de la face */
100			vap[0] = coord[x(nump)] - coord[x(tabele[3*nb_ele])] ;
101			vap[1] = coord[y(nump)] - coord[y(tabele[3*nb_ele])] ;
102			vap[2] = coord[z(nump)] - coord[z(tabele[3*nb_ele])] ;
103			if (PROSCA(vap,vn) < EPSILON)
104			{
105			for (i=0;i<3;i++) vn[i] = -1. * vn[i] ;
106			i = tabele[3*nb_ele+1] ;
107			tabele[3nb_ele+1] = tabele[3nb_ele+2] ;
108			tabele[3*nb_ele+2] = i ;
109			}
110			/* point optimal */
111			tab_opt[x(nb_ele)] = xi + hauteur * vn[0] ;
112			tab_opt[y(nb_ele)] = yi + hauteur * vn[1] ;
113			tab_opt[z(nb_ele)] = zi + hauteur * vn[2] ;
114
115			nb_ele ++ ;
116			}
117			connec = connec->suivant ;
118			}
119			if (nb_ele == 0) return(VRAI) ;
120			/* calcul de la qualite de chaque tetraedre */
121			cmin = 1. ;
122			for (i=0;i<nb_ele;i++)
123			{
124			/* calcul du critere */
125			tab_crit[i] = m3d_e_qual2(gest->coord,tabele[3i],tabele[3i+1],tabele[3*i+2],nump) ;
126			if (tab_crit[i]<EPSILON)
127			{
128			// if (debug) aff_text("critere nul M3D_MOVE\n ") ;
129			m3d_erreur(ERR_SYST) ;
130			return(FAUX) ;
131			}
132			if (tab_crit[i] < cmin) cmin = tab_crit[i] ;
133			}
134			/* calcul des coordonnees du point ideal */
135			gamma = 0. ;
136			for (i=0;i<nb_ele;i++)
137			{
138			gamma = gamma + 1./(tab_crit[i]tab_crit[i]) ;/ somme des carres des inverses */
139			}
140			/* calcul du point optimal pour la boule */
141			gamma = 1./gamma ;
142			vec_opt[0] = 0. ;
143			vec_opt[1] = 0. ;
144			vec_opt[2] = 0. ;
145
146			for (i=0;i<nb_ele;i++)
147			{
148			coeff = gamma/(tab_crit[i]*tab_crit[i]) ;
149			vec_opt[0] = vec_opt[0] + coeff * tab_opt[x(i)] ;
150			vec_opt[1] = vec_opt[1] + coeff * tab_opt[y(i)] ;
151			vec_opt[2] = vec_opt[2] + coeff * tab_opt[z(i)] ;
152			}
153
154			delta[0] = vec_opt[0] - coord[x(nump)] ;
155			delta[1] = vec_opt[1] - coord[y(nump)] ;
156			delta[2] = vec_opt[2] - coord[z(nump)] ;
157			/* sauvegarde du critere mini */
158			crit_save = cmin ;
159
160			/* initialisation des bornes du parametre */
161			bmin = 0. ;
162			bmax = 1. ;
163			sens = 1 ;/* augmnentation */
164
165			vresu[0] = vec_save[0] ;
166			vresu[1] = vec_save[1] ;
167			vresu[2] = vec_save[2] ;
168
169			for (i=0;i<5;i++)
170			{
171			alpha = 0.5 * (bmin + bmax) ;
172			coord[x(nump)] = vec_save[0] + alpha * delta[0] ;
173			coord[y(nump)] = vec_save[1] + alpha * delta[1] ;
174			coord[z(nump)] = vec_save[2] + alpha * delta[2] ;
175
176			/* evaluation de la qualite pour la valeur alpha */
177			/* on evalue la qualite cmin pour alpha */
178			calpha = 1. ;
179			for (j=0;j<nb_ele;j++)
180			{
181			/* calcul du critere */
182			tab_crit[j] = m3d_e_qual2(gest->coord,tabele[3j],tabele[3j+1],tabele[3*j+2],nump) ;
183			if (tab_crit[j] < calpha) calpha = tab_crit[j] ;
184			}
185			/* la valeur du critere pour alpha est superieure, il faut la stocker */
186			if (calpha > crit_save)
187			{
188			vresu[0] = coord[x(nump)] ;
189			vresu[1] = coord[y(nump)] ;
190			vresu[2] = coord[z(nump)] ;
191			crit_save = calpha ;
192			}
193
194			/* on evalue la qualite c1 pour alpha - eps */
195			eps = (bmax-bmin)/50. ;
196			coord[x(nump)] = vec_save[0] + (alpha - eps) * delta[0] ;
197			coord[y(nump)] = vec_save[1] + (alpha - eps) * delta[1] ;
198			coord[z(nump)] = vec_save[2] + (alpha - eps) * delta[2] ;
199			calpha_moins_eps = 1. ;
200			for (j=0;j<nb_ele;j++)
201			{
202			/* calcul du critere */
203			tab_crit[j] = m3d_e_qual2(gest->coord,tabele[3j],tabele[3j+1],tabele[3*j+2],nump) ;
204			if (tab_crit[j] < calpha_moins_eps) calpha_moins_eps = tab_crit[j] ;
205			}
206
207			coord[x(nump)] = vec_save[0] + alpha * delta[0] ;
208			coord[y(nump)] = vec_save[1] + alpha * delta[1] ;
209			coord[z(nump)] = vec_save[2] + alpha * delta[2] ;
210
211			if (calpha_moins_eps < calpha)
212			{
213			bmin = alpha ;
214			/* bmax est inchange */
215			}
216			else
217			{
218			/* bmin est inchange */
219			bmax = alpha ;
220			}
221			}
222			coord[x(nump)] = vresu[0] ;
223			coord[y(nump)] = vresu[1] ;
224			coord[z(nump)] = vresu[2] ;
225
226			return(VRAI) ;
227			}
228
229
230
231
232
233
234