Trouver surface de triangles de maillage tétraédrique

**mari2** · 20/08/2016, 21h32

je veut convertir cette soubritine ecrit en matlab en fortran
elle trouve les surface de triangles de maillage tétraédrique (SURFTRI)

Code matlab :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
function tri=surftri(p,t)
%SURFTRI Find surface triangles from tetrahedra mesh
%   TRI=SURFTRI(P,T)
 
%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.
 
% Form all faces, non-duplicates are surface triangles
faces=[t(:,[1,2,3]);
       t(:,[1,2,4]);
       t(:,[1,3,4]);
       t(:,[2,3,4])];
node4=[t(:,4);t(:,3);t(:,2);t(:,1)];
faces=sort(faces,2);
[foo,ix,jx]=unique(faces,'rows');
vec=histc(jx,1:max(jx));
qx=find(vec==1);
tri=faces(ix(qx),:);
node4=node4(ix(qx));
 
% Orientation
v1=p(tri(:,2),:)-p(tri(:,1),:);
v2=p(tri(:,3),:)-p(tri(:,1),:);
v3=p(node4,:)-p(tri(:,1),:);
ix=find(dot(cross(v1,v2,2),v3,2)>0);
tri(ix,[2,3])=tri(ix,[3,2]);

merci d'avance

**dva2tlse** · 24/08/2016, 09h49

Bonjour,
moi je n'y connais presque rien en matlab, mais pas trop mal en fortran, donc si ça t'intéresse toujours et si tu es prêt à y passer du temps et des tâtonnements (donc pas du "TOUT CUIT") je suis partant.
David

**mari2** · 24/08/2016, 13h00

merci pour votre réponse et pour Votre volonté , je suis prête aussi pour des tâtonnements avec vous

**dva2tlse** · 24/08/2016, 13h38

rebonjour,
il faut me donner le plus de détails possible sur ce que font toute les fonctions matlab utilisées pour qu'on puisse les reconstruire en fortran.
Dans quel code éléments-finis cela se passe t'il.
Envoyez moi les carte de définition d'un ou deux tétraèdres et des nœuds.
Je ferai un petit bout de fortran et vous laisserai terminer dès que ce sera réaliste.
David

**mari2** · 24/08/2016, 14h07

j'utilise la Méthode des Volumes de Contrôle à base d’Elément Finis (MVCEF) pour la discrétisation tridimensionnelle

P =matrice des coordonées des noeux
T = matrice des élement qui contient numéro des noeux formant (4 noeux)

pour les fonctions matlab utilisées:
B = sort(A,dim) retourne les éléments triés de A le long de la dimension dim. Par exemple, si A est une matrice, puis sort(A,2) trie les éléments dans chaque ligne

[C,ia,ic] = unique(A,'rows') traite chaque ligne de A comme une entité unique et retourne les lignes uniques de A. Les lignes du tableau C sont dans un ordre trié.
L’option 'rows' ne supporte pas les baies de la cellule , aussi indice de retours vecteurs ia et ic, telle que C = A(ia,: ) et A = C(ic,: ).

bincounts = histc(x,binranges) compte le nombre de valeurs de x qui se trouvent dans chaque gamme de bac spécifié. L’entrée, binranges, détermine les points de terminaison pour chaque bac. La sortie, bincounts, contient le nombre d’éléments de x dans chaque emplacement.

k = find(X)renvoie un vecteur contenant les indices linéaires de chaque élément différent de zéro dans le tableau X.

C = cross(A,B) retourne le produit Croix de A et B. Nom : téléchargement.png
Affichages : 499
Taille : 2,2 Ko

Nom : téléchargement.png
Affichages : 499
Taille : 2,2 Ko

C = dot(A,B) calcule le produit scalaire dot de A et B.

et merci d'avance

**mari2** · 24/08/2016, 14h54

j'ai essayé les 2 premiers lignes

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
open(17,File="c:tetareda")
       open(16,File="c:faces")
       do j=1,ntet	 
	 read(17,*)element(j),element_type(j),Nnub1(j),Nnub2(j),Nnub3(j),
     & Nnub4(j)
	 end do        
	 do j=1,ntet      
	 write(16,*)Nnub1(j),Nnub2(j),Nnub3(j)
       end do 
       do j=1,ntet     
	 write(16,*)Nnub1(j),Nnub2(j),Nnub4(j)
       end do 
       do j=1,ntet       
	 write(16,*)Nnub1(j),Nnub3(j),Nnub4(j)
       end do 
       do j=1,ntet       
	 write(16,*)Nnub2(j),Nnub3(j),Nnub4(j)
       end do 
       close(17)
        close(16)
************************************************************************************		
 
       open(17,File="c:tetareda")
       open(20,File="c:node4")
       do j=1,ntet	 
	 read(17,*)element(j),element_type(j),Nnub1(j),Nnub2(j),Nnub3(j),
     & Nnub4(j)
	 end do        
	 do j=1,ntet      
	 write(20,*)Nnub4(j)
       end do 
       do j=1,ntet     
	 write(20,*)Nnub3(j)
       end do 
       do j=1,ntet       
	 write(20,*)Nnub2(j)
       end do 
       do j=1,ntet       
	 write(20,*)Nnub1(j)
       end do 
       close(17)
        close(20)

**dva2tlse** · 24/08/2016, 19h48

Ouais, super et bravo; c'est bien d'avoir un peu commencé; maintenant j'aimerais avoir les fichiers c:tetareda et c:faces, ou tout au moins des lignes décrivant deux ou quatre tetras adjacents et leur faces; comme ça je pourrais essayer des trucs avec votre début de source...
à+,
David

**mari2** · 24/08/2016, 21h05

c:tetareda

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
19688           4         652       11346       13331       13451
       19689           4       10010       12017       10711       13732
       19690           4        5133       10853       10399       11715
       19691           4       10159        9394       11112       12316
       19692           4        6552        6921       10429       12764
       19693           4        9747       12102       11074       12685
       19694           4       11043       13329       10386       15430
       19695           4        5913       12018       10460       13083
       19696           4       10622       11056        9862       11371
       19697           4       10719       10740        9920       12115
       19698           4        6826       11584       11465       13992
       19699           4       13665       13674       10634       15945
       19700           4        9529       11569       10449       14501
       19701           4        8137       12196       10150       13846
       19702           4       10190       11214        9657       11614
       19703           4        7778       10537       12120       12488
       19704           4       11434       14098       11206       14902
       19705           4        9518       10507       11161       15327
       19706           4        7195       13400       13214       13991
       19707           4        3816       10174        4646       12504
       19708           4       10500       12564       10349       15910
       19709           4       10309       10773        9704       11788

c:faces

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
652       11346       13331
       10010       12017       10711
        5133       10853       10399
       10159        9394       11112
        6552        6921       10429
        9747       12102       11074
       11043       13329       10386
        5913       12018       10460
       10622       11056        9862
       10719       10740        9920
        6826       11584       11465
       13665       13674       10634
        9529       11569       10449
        8137       12196       10150
       10190       11214        9657
        7778       10537       12120
       11434       14098       11206
        9518       10507       11161
        7195       13400       13214
        3816       10174        4646
       10500       12564       10349
       10309       10773        9704

node4

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Merci beaucoup pour votre aide et pour votre intention

**dva2tlse** · 24/08/2016, 21h29

j'ai un tout petit peu avancé en arrangeant quelques lignes du source pour qu'il lise les fichiers d'entrée dès que vous m'en aurez quelques lignes.
David
PS: demain au boulot j'ai accès au site. (si j'ai du temps)

**dva2tlse** · 25/08/2016, 08h24

Bonjour,
je suis arrivé hyper-tôt au boulot ce matin pour faire tourner le programme ici et avancer.
J'arrive à lire les éléments, mais pour calculer les surfaces des faces, il faudrait les coordonnées des nœuds;
Cf le commentaire à la fin du source :

c Au lieu d'écrire les Numéros des noeuds de chaque face, il suffit de calculer l'aire double de la face,
c en multipliant entre eux deux vecteurs côtés d'une face; mais pour ça il faut avoir les coord's des noeuds.

à+,
David

Pourquoi veux tu passer de matlab au fortran ?
Quelle sorte d'études et où ?
D.

**mari2** · 25/08/2016, 14h41

Bonjour David , merci beaucoup pour votre aide
Mais malheureusement, je n'arrive pas a lire le fichier que vous m'avez envoyé,
Est-ce que vous pouvez le renvoyer sous forme d'un code.

Ce travail entre dans mon travail de thèse de doctorat en Tunisie, j'ai besoin de cette subroutine ,mais je l'ai trouvé qu'en matlab et mon code est en fortran

merci encore une fois

**dva2tlse** · 25/08/2016, 17h01

Bonjour Mari, (c'est quoi votre vrai prénom ?)
voici le code dans l'état actuel; il faudra ajouter la lecture des coordonnées des nœuds dans un autre fichier, après quoi il sera possible de faire le produit vectoriel de deux côtés de chaque face, et la norme du vecteur obtenu sera le double de la surface voulue.
Voyez si vous pouvez continuer seule, sinon avec les coordonnées je pourrais avancer encore un peu.
David
PS: Quel sujet de thèse ?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
c23456 cp tetra_f.txt tetra.f
c23456 cp tetra.f tetra_f.txt
c23456 gfortran tetra.f -o tetra
c23456 ./tetra
c23456----------------------------------------------------------------|
      program tetra
      implicit none
      integer j, ntet
      integer element(22)
      integer element_type(22)
      integer Nnub1(22)
      integer Nnub2(22)
      integer Nnub3(22)
      integer Nnub4(22)
c23456----------------------------------------------------------------|
      print'(/A)', '--------------------------------- debut '//
     +             'de tetra ------------------------------|'
      ntet=22
c23456----------------------------------------------------------------|
c Ouvertures sous linux à la maison :
      open(17, File="/home/david/F90/tetareda")
c
c23456----------------------------------------------------------------|
c Ouvertures sous linux au boulot :
c     open(17, File="/data/home/dvanaelst/F90/FDEV/tetareda")
c
c23456----------------------------------------------------------------|
c Ouvertures sous DOS/windows :
c     open(17, File="c:\tetareda")
c     open(16, File="c:\faces")
c23456----------------------------------------------------------------|
c Lecture des éléments :
      print'(/A)', 'j, element(j), element_type(j), '//
     +'Nnub1(j), Nnub2(j), Nnub3(j), Nnub4(j)'
      do j=1, ntet	
	 read(17, *)
     +element(j), element_type(j),
     +Nnub1(j), Nnub2(j), Nnub3(j), Nnub4(j)
          print'(8I8)', j, element(j), element_type(j),
     +                  Nnub1(j), Nnub2(j), Nnub3(j), Nnub4(j)
      end do
c23456----------------------------------------------------------------|
c Lecture (ou écriture ???) des noeuds des faces :
      do j=1, ntet
c???     write(16, * )Nnub1(j), Nnub2(j), Nnub3(j) ! face opposée au noueud 4
      end do
c23456
      do j=1, ntet
c???     write(16, *) Nnub1(j), Nnub2(j), Nnub4(j) ! face opposée au noueud 3
      end do
c23456
      do j=1, ntet
c???     write(16, *) Nnub1(j), Nnub3(j), Nnub4(j) ! face opposée au noueud 2
      end do
c23456
      do j=1, ntet
c???     write(16, *) Nnub2(j), Nnub3(j), Nnub4(j) ! face opposée au noueud 1
      end do
c23456----------------------------------------------------------------|
c Au lieu d'écrire les Nos des noeuds de chaque face, il suffit de calculer l'aire double de la face,
c en multipliant entre eux deux vecteurs côtés d'une face; mais pour ça il faut avoir les coord's des noeuds.
c
      close(17)
c     close(16)
      print'(A/)', '----------------------------------- fin '//
     +             'de tetra ------------------------------|'
      stop
      end

**mari2** · 26/08/2016, 18h11

Merci David pour votre aide, et pour mon vrai nom, c'est Mariem

Pour moi, c'était ses lignes qui sont un peu difficiles pour les convertir puisque je suis débutante en fortran et je ne trouve pas les fonctions intrinsèques exactes.

Code MATLAB :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
[foo,ix,jx]=unique(faces,'rows');
vec=histc(jx,1:max(jx));
qx=find(vec==1);
tri=faces(ix(qx),:);
node4=node4(ix(qx));

**dva2tlse** · 26/08/2016, 19h59

Bonjour Mariem,
moi je pense qu'il suffit d'avoir les coordonnées des nœuds, comme celles-ci, fictives que j'ai inventées pour quelques nœuds dans un fichier appelé coords, où j'ai mis les coordonnées des nœuds des trois premiers éléments :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
      652        0.	1.	2.
      5133       3.	4.	5.
      10010      6.	7.	8.
      10711      9.	10.	11.
      10853      12.	13.	14.
      10399      15.	16.	17.
      11346      18.	19.	20.
      12017      21.	22.	23.
      13331      24.	25.	26.

puis on les appelle XN(1) à XN(9) puisqu'il y a neuf nœuds dans mon exemple, et YN(1) à YN(9) et ZN(1) à ZN(9)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
c donc il faut les déclarer :
      integer N(9) ! ça c'est pour les numéros de nœuds
      real XN(9), YN(9), ZN(9)
c et il faut les lire au début du programme :
      open(1, file="/home/david/F90/coords")
      do i=1, 9
          read(1, N(i), XN(N(i)), YN(N(i)), ZN(N(i)))
      enddo
c et ensuite d'après les faces on peut créer les vecteurs des arêtes.

Bon courage si vous arrivez à avancer avec ça.
David

Bon bin j'ai un peu avancé moi même, mais il reste à faire une fonction qui donne la surface en fonction des trois noeuds dont on a les coordonnées. Bon courage encore !

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
 
c23456 cp tetra_f.txt tetra.f
c23456 cp tetra.f tetra_f.txt
c23456 gfortran tetra.f -o tetra
c23456 ./tetra
c23456----------------------------------------------------------------|
      program tetra
      implicit none
      character*27 line
      integer i, j, ntet
      integer element(22)
      integer element_type(22)
      integer Nnub1(22)
      integer Nnub2(22)
      integer Nnub3(22)
      integer Nnub4(22)
      integer N(9) ! ça c'est pour les numéros de nœuds
      real XN(9), YN(9), ZN(9) ! ça c'est pour les coordonnées
c23456----------------------------------------------------------------|
      print'(/A)', '--------------------------------- debut '//
     +             'de tetra ------------------------------|'
      ntet=22
c23456----------------------------------------------------------------|
c Ouvertures sous linux à la maison :
      open(1, file="/home/david/F90/coords")
      open(17, File="/home/david/F90/tetareda")
c
c23456----------------------------------------------------------------|
c Ouvertures sous linux au boulot :
c     open(1, File="/data/home/dvanaelst/F90/FDEV/coords")
c     open(17, File="/data/home/dvanaelst/F90/FDEV/tetareda")
c
c23456----------------------------------------------------------------|
c Ouvertures sous DOS/windows :
c     open(1, file="c:\coords")
c     open(17, File="c:\tetareda")
c23456----------------------------------------------------------------|
c Lecture des éléments :
      print'(/A)', 'j, element(j), element_type(j), '//
     +'Nnub1(j), Nnub2(j), Nnub3(j), Nnub4(j)'
      do j=1, ntet	
	 read(17, *)
     +element(j), element_type(j),
     +Nnub1(j), Nnub2(j), Nnub3(j), Nnub4(j)
          print'(8I8)', j, element(j), element_type(j),
     +                  Nnub1(j), Nnub2(j), Nnub3(j), Nnub4(j)
      end do
c23456----------------------------------------------------------------|
c Lecture des noeuds des faces :
      do i=1, 9
          read(1, '(A27)') line
          read(line, *) N(i), XN(i), YN(i), ZN(i)
          write(*, *) N(i), XN(i), YN(i), ZN(i) ! pour contrôler
      enddo
      close(1)
c23456----------------------------------------------------------------|
c Lecture des noeuds des faces :
      do j=1, ntet
c??? Nnub1(j), Nnub2(j), Nnub3(j) ! face opposée au noeud 4
c arête 1=Nnub1(j), Nnub2(j)
c arête 2=Nnub2(j), Nnub3(j)
      end do
c23456
      do j=1, ntet
c??? Nnub1(j), Nnub2(j), Nnub4(j) ! face opposée au noeud 3
c arête 1=Nnub1(j), Nnub2(j)
c arête 2=Nnub2(j), Nnub4(j)
      end do
c23456
      do j=1, ntet
c??? Nnub1(j), Nnub3(j), Nnub4(j) ! face opposée au noeud 2
c arête 1=Nnub1(j), Nnub3(j)
c arête 2=Nnub3(j), Nnub4(j)
      end do
c23456
      do j=1, ntet
c??? Nnub2(j), Nnub3(j), Nnub4(j) ! face opposée au noeud 1
c arête 1=Nnub2(j), Nnub3(j)
c arête 2=Nnub3(j), Nnub4(j)
      end do
c23456----------------------------------------------------------------|
c Au lieu d'écrire les Nos des noeuds de chaque face, il suffit de calculer l'aire double de la face,
c en multipliant entre eux deux vecteurs côtés d'une face; mais pour ça il faut avoir les coord's des noeuds.
c
      close(17)
c     close(16)
      print'(A/)', '----------------------------------- fin '//
     +             'de tetra ------------------------------|'
      stop
      end
c

**dva2tlse** · 27/08/2016, 16h05

Bonjour Mariem,
je suis en train d'avancer sur le programme, mais c'est en train de devenir une véritable usine à gaz. Mais ça marche pas mal et j'ai presque fini.
Par contre j'ai besoin d'un fichier où il y a les coordonnées des nœuds, comme celui que j'ai créé ci dessus.
à+,
David

**dva2tlse** · 28/08/2016, 16h34

Bonjour Mariem,
voici de nouveau le programme que j'ai fait, et qui utilise le fichier de coordonnées des nœuds que j'ai mentionné plus haut.
Il marche bien, mais pour une seule face du tétraèdre, qui fait face au nœud 4; pour les autres faces, il faut d'abord faire recalculer les longueurs des arêtes.
Bon courage,
David

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
c23456 cp tetra_f.txt tetra.f
c23456 cp tetra.f tetra_f.txt
c23456 gfortran tetra.f -o tetra
c23456 echo; echo; echo; ./tetra
c23456 echo; echo; echo; gdb ./tetra<RUN
c23456----------------------------------------------------------------|
      program tetra
      implicit none
      integer Nnoeud
      integer Nelem
      parameter(Nnoeud=12)
      parameter(Nelem=22)
      character*80 line
      integer i
      integer j
      integer k
      integer element(Nelem) ! il y a Nelem elements
      integer element_type(Nelem)
      integer N1(Nelem)
      integer N2(Nelem)
      integer N3(Nelem)
      integer N4(Nelem)
      integer N(Nnoeud) ! ça c'est pour les numéros de nœuds
      integer tmp
      real a12, a14, a23 ! longueurs des côtés N1-N2
      real p1, p2, p3, p4 ! demi-sommes des longueurs des côtés de la face opposée a N
      real s1, s2, s3, s4 ! surface de la face opposée a N
      real XN(Nnoeud), YN(Nnoeud), ZN(Nnoeud) ! ça c'est pour les coordonnées du noeud N=1->9
      common//N
c23456----------------------------------------------------------------|
      print'(/A)', '--------------------------------- debut '//
     +             'de tetra ------------------------------|'
c23456----------------------------------------------------------------|
c Ouvertures sous linux à la maison :
      open(1, file="/home/david/F90/coords")
      open(2, File="/home/david/F90/tetareda")
c
c23456----------------------------------------------------------------|
c Ouvertures sous linux au boulot :
c     open(1, File="/data/home/dvanaelst/F90/FDEV/coords")
c     open(2, File="/data/home/dvanaelst/F90/FDEV/tetareda")
c
c23456----------------------------------------------------------------|
c Ouvertures sous DOS/windows :
c     open(1, file="c:\coords")
c     open(2, File="c:\tetareda")
c23456----------------------------------------------------------------|
c Lecture des coordonnées des noeuds : (9/10 noeuds seulement, pour qq elts)
      do i=1, Nnoeud
 10       read(1, '(A)', end=15) line
          if(line(1:1)=='#') print'(/A)', line
          if(line(1:1)=='#') goto 10
          write(6, '(3X,A)') trim(line)
          read(line, *) N(i), XN(i), YN(i), ZN(i)
          write(6, '(I2,I6,F7.2,F12.2,F12.2)')
     +i, N(i), XN(i), YN(i), ZN(i) ! pour contrôler
      enddo
 15   close(1)
          print*, '  N(12)=',  N(12), ', XN(12)=', ZN(12)
          print*, ' YN(12)=', YN(12), ', ZN(12)=', ZN(12)
c23456----------------------------------------------------------------|
c Lecture des éléments :
      do j=1, 3 ! Nelem
 20       read(2, '(A)', end=25) line
          if(line(1:1)=='#') print'(/A)', line
          if(line(1:1)=='#') goto 20
          print'(A)', line
          read(line, *) element(j), element_type(j),
     +N1(j), N2(j), N3(j), N4(j)
          print'(2I6, 5I12)', j, element(j), element_type(j),
     +N1(j), N2(j), N3(j), N4(j) ! pour contrôler
      enddo
 25   close(2)
         print*, ' N4(2)=', N4(2), ', N(12)=', N(12)
c23456----------------------------------------------------------------|
c Formule de Héron, permet le calcul de l'aire S, connaissant les longueurs
c des trois côtés a, b et c d'un triangle, et donc aussi leur demi-somme p :
c     S =sqrt(p x ( p − a ) x ( p − b ) x ( p − c ))
c23456----------------------------------------------------------------|
c Traitement des noeuds des faces :
c Il me faut connaître le N° d'ordre i d'un noeud N, que je vais appeler,
c K(N), avec N=N1(j) ou N2(j) ou N3(j) ou N4(j), j étant un N° d'élément.
      do j=1, 3 ! Nelem
          print'(/A,I2)', 'Elt j=', j
          print'(4(A,I6))', 'Noeuds : N1(j)=', N1(j), 
     +', N2(j)=', N2(j), 
     +', N3(j)=', N3(j), 
     +', opposés à N4(j)=', N4(j)
          print'(A,I2, A,I6, A,I6, 3(/6X,A,I6, A,I6 ,X))',' j=', j,
     +', N1(j)=', N1(j), ', K(N1(j))=', K(N1(j)),
     +' N2(j)=', N2(j), ', K(N2(j))=', K(N2(j)),
     +' N3(j)=', N3(j), ', K(N3(j))=', K(N3(j)),
     +' N4(j)=', N4(j), ', K(N4(j))=', K(N4(j))
          print'(3(A, F6.2, X))', 
     +' XN(K(N1(j)))=', XN(K(N1(j))), 
     +', YN(K(N1(j)))=', YN(K(N1(j))),
     +', ZN(K(N1(j)))=', ZN(K(N1(j)))
          print'(3(A, F6.2, X))', 
     +' XN(K(N2(j)))=', XN(K(N2(j))), 
     +', YN(K(N2(j)))=', YN(K(N2(j))),
     +', ZN(K(N2(j)))=', ZN(K(N2(j)))
          print'(3(A, F6.2, X))', 
     +' XN(K(N3(j)))=', XN(K(N3(j))), 
     +', YN(K(N3(j)))=', YN(K(N3(j))),
     +', ZN(K(N3(j)))=', ZN(K(N3(j)))
c arête 12=N1(j), N2(j)
          a12=sqrt( (XN(K(N2(j)))-XN(K(N1(j))) )**2
     +            + (YN(K(N2(j)))-YN(K(N1(j))) )**2
     +            + (ZN(K(N2(j)))-ZN(K(N1(j))) )**2 )
      print'(A,I2, 2(A,I6))', 'Elt j=', j, 
     +', arête a12, de N1(j)=', N1(j), ' à N2(j)=', N2(j)
      print'(2(A,I6), A,F8.3)', 'soit de K(N1(j))=', K(N1(j)), 
     +' à K(N2(j))=', K(N2(j)), '; long.=', a12
c arête 14=N1(j), N4(j) (+)
          a14=sqrt( (XN(K(N4(j)))-XN(K(N1(j))) )**2
     +            + (YN(K(N4(j)))-YN(K(N1(j))) )**2
     +            + (ZN(K(N4(j)))-ZN(K(N1(j))) )**2 )
      print'(A,I2, 2(A,I6))', 'Elt j=', j, 
     +', arête a14, de N1(j)=', N1(j), ' à N4(j)=', N4(j)
      print'(2(A,I6), A,F8.3)', 'soit de K(N1(j))=', K(N1(j)), 
     +' à K(N4(j))=', K(N4(j)), '; long.=', a14
c arête 23=N2(j), N3(j)
          a23=sqrt( (XN(K(N3(j)))-XN(K(N2(j))) )**2
     +            + (YN(K(N3(j)))-YN(K(N2(j))) )**2
     +            + (ZN(K(N3(j)))-ZN(K(N2(j))) )**2 )
      print'(A,I2, 2(A,I6))', 'Elt j=', j, 
     +', arête a23, de N2(j)=', N2(j), ' à N3(j)=', N3(j)
      print'(2(A,I6), A,F8.3)', 'soit de K(N2(j))=', K(N2(j)), 
     +' à K(N3(j))=', K(N3(j)), '; long.=', a23
          p4=(a12+a14+a23)/2
          s4=sqrt(p4 * (p4-a12) * (p4-a14) * (p4-a23) )/2.
      print'(A,I2, 3(A,I6))', 'Elt j=', j, 
     +', face opposée à N4, N1(j)=', N1(j), ', N2(j)=', N2(j), 
     +' et N3(j))=', N3(j)
      print'(3(A,I6), A,F8.3)', ' soit de K(N3(j))=', N3(j), 
     +' à K(N2(j))=', K(N2(j)), ' et K(N3(j))=', K(N3(j)), 
     +'; surf.=', s4
          print*, ' surf.=', s4
c
c??? N1(j), N2(j), N4(j) ! face opposée au noeud 3
c arête 12=N1(j), N2(j) (-)
c arête 14=N1(j), N4(j) (+/-)
c arête 24=N2(j), N4(j)
c
c??? N1(j), N3(j), N4(j) ! face opposée au noeud 2
c arête 13=N1(j), N3(j)
c arête 14=N1(j), N4(j) (+/-)
c arête 34=N3(j), N4(j)
c
c??? N2(j), N3(j), N4(j) ! face opposée au noeud 1
c arête 23=N2(j), N3(j) (-)
c arête 24=N2(j), N4(j) (+/-)
c arête 34=N3(j), N4(j) (-)
      end do
c23456----------------------------------------------------------------|
      print'(A/)', '----------------------------------- fin '//
     +             'de tetra ------------------------------|'
      stop
      end
c
c23456----------------------------------------------------------------|
      function K(NN)
      implicit none
      integer Nnoeud
      parameter(Nnoeud=12)
      integer i
      integer K ! indice du noeud passé en argument
      integer N(Nnoeud)
      integer NN ! N° du noeud dont on cherche l'indice
      common//N
      do i=1, Nnoeud
          if(NN.eq.N(i))then
              K=i
              return
          endif
      enddo
      print'(A,I2,  A,I6)', 'Err. fonc.K(NN), i=', i, ', NN=', NN
      stop
      end
c23456----------------------------------------------------------------|

Trouver surface de triangles de maillage tétraédrique

Fortran

Discussions similaires

Partager

Partager