ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/acismesh/m3d_angle_front2.cpp
Revision: 1
Committed: Mon Jun 11 22:53:07 2007 UTC (17 years, 11 months ago)
File size: 3288 byte(s)
Log Message:

File Contents

# User Rev Content
1 1 /*****************************************************************
2    
3     m3d_angle_front2.cpp Type:Func
4    
5     Calcul l angle entre un element du front 3d et un noeud lie
6    
7     Date de creation : 10-3-1999 11 :32 :3
8     Derniere version : 10-3-1999 11 :32 :3
9    
10     Vincent FRANCOIS
11    
12     *****************************************************************/
13    
14    
15    
16    
17    
18     /**************************/
19     /* include */
20     #include <stdio.h>
21     #include <string.h>
22     #include <stdlib.h>
23     #include <math.h>
24     #include "const.h"
25     #include "memoire.h"
26     #include "struct.h"
27     #include "prototype.h"
28    
29     /**************************/
30     /* variables globales */
31     extern struct environnement env;
32     extern struct s_mesh *mesh;
33    
34    
35    
36     /**************************/
37     /* programme principal */
38    
39     float m3d_angle_front2(struct s_front3d *ft,struct s_segment *seg)
40     {
41     int n1,n2,n3,n4;
42     struct s_noeud *no1,*no2,*no3,*no4;
43     float xproj,yproj,zproj,det,xsi,eta,psi;
44     float n1n2[4],n1n3[4],n1n4[4],n1n[4],n1nproj[4];
45     float angle,ps;
46    
47     if (seg->n1==ft->n1) {n1=ft->n1;n2=ft->n2;n3=ft->n3;n4=seg->n2;}
48     else if (seg->n1==ft->n2) {n1=ft->n2;n2=ft->n3;n3=ft->n1;n4=seg->n2;}
49     else if (seg->n1==ft->n3) {n1=ft->n3;n2=ft->n1;n3=ft->n2;n4=seg->n2;}
50     else if (seg->n2==ft->n1) {n1=ft->n1;n2=ft->n2;n3=ft->n3;n4=seg->n1;}
51     else if (seg->n2==ft->n2) {n1=ft->n2;n2=ft->n3;n3=ft->n1;n4=seg->n1;}
52     else if (seg->n2==ft->n3) {n1=ft->n3;n2=ft->n1;n3=ft->n2;n4=seg->n1;}
53    
54    
55    
56     no1=ADRESSE(n1,noeud,mesh->);
57     no2=ADRESSE(n2,noeud,mesh->);
58     no3=ADRESSE(n3,noeud,mesh->);
59     no4=ADRESSE(n4,noeud,mesh->);
60    
61    
62     VEC(n1n2,no1,no2);
63     VEC(n1n3,no1,no3);
64     VEC(n1n4,no1,no4);
65     PVEC(n1n,n1n2,n1n3);
66    
67     det=n1n2[0]*n1n3[1]*n1n[2]+n1n2[1]*n1n3[2]*n1n[0]+n1n2[2]*n1n3[0]*n1n[1]
68     -n1n2[2]*n1n3[1]*n1n[0]-n1n2[0]*n1n3[2]*n1n[1]-n1n2[1]*n1n3[0]*n1n[2];
69     xsi=n1n4[0]*n1n3[1]*n1n[2]+n1n4[1]*n1n3[2]*n1n[0]+n1n4[2]*n1n3[0]*n1n[1]
70     -n1n4[2]*n1n3[1]*n1n[0]-n1n4[0]*n1n3[2]*n1n[1]-n1n4[1]*n1n3[0]*n1n[2];
71     eta=n1n2[0]*n1n4[1]*n1n[2]+n1n2[1]*n1n4[2]*n1n[0]+n1n2[2]*n1n4[0]*n1n[1]
72     -n1n2[2]*n1n4[1]*n1n[0]-n1n2[0]*n1n4[2]*n1n[1]-n1n2[1]*n1n4[0]*n1n[2];
73     psi=n1n2[0]*n1n3[1]*n1n4[2]+n1n2[1]*n1n3[2]*n1n4[0]+n1n2[2]*n1n3[0]*n1n4[1]
74     -n1n2[2]*n1n3[1]*n1n4[0]-n1n2[0]*n1n3[2]*n1n4[1]-n1n2[1]*n1n3[0]*n1n4[2];
75     xsi=xsi/det;
76     eta=eta/det;
77     psi=psi/det;
78    
79     xproj=(1-xsi-eta)*no1->x+xsi*no2->x+eta*no3->x;
80     yproj=(1-xsi-eta)*no1->y+xsi*no2->y+eta*no3->y;
81     zproj=(1-xsi-eta)*no1->z+xsi*no2->z+eta*no3->z;
82    
83     n1nproj[0]=xproj-no1->x;
84     n1nproj[1]=yproj-no1->y;
85     n1nproj[2]=zproj-no1->z;
86     NORME(n1n4);
87     NORME(n1nproj);
88     ps=PSCA(n1n4,n1nproj);
89     if (ps>1.) ps=1.;
90     if (ps<-1.) ps=-1.;
91     angle=(float)acos((double)ps);
92     if ((xsi<0.)||(eta<0.)) angle=angle+PI/2.;
93     if (psi<0.) angle=2*PI-angle;
94     return(angle);
95     }