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

# Content
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 }