RANDOM NUMBER probleme

**BAmal** · 31/10/2013, 15h31

Salut,
je voudrais bien générer des positions aléatoires à l aide d la subroutine Random number ,mais je voudrais aussi imposer une distance minimale qui sépare n importe quels deux points de mon maillage ,quelq'un peut m aider sur ce problème ,voilà ce que j ai bien réussi à faire ,merci

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
30      do   m=1,3
			    call random_number(x)
		        call random_number(y)
          10     do   i=1,n*n
          20     do   j=1,n*n
			    dx=x(j)-x(i)
	            dy=y(j)-y(i)
			    rad=sqrt(dx**2+dy**2)
				write(300,*) i,j,rad
		        enddo
				if(rad<dd) then  ! dd est la distance minimale que j ai imposée 
				!write(*,*) rad
				go to 10
				else
				go to 30
				!write(*,*) x(i),y(j) 
 
				endif
				!write(*,*) rad
			    enddo
		        enddo

**xflr6** · 02/11/2013, 16h23

Bonjour,

Même si j'en ai une idée, j'ai un peu de mal à comprendre ce que tu essaies de faire.
1) Essaies-tu de construire un réseau de points répartis aléatoirement (avec des distances inter-points "correctes"), c'est à dire un pb de type sphères dures en 2D ?
2) Ou alors tu as déjà un réseau, et tu désires tirer certains points de ce réseau aléatoirement ?

J'ai l'impression que c'est l'option 1, right ? Pour celle-ci, il faut
- tirer le premier point
- tirer les suivants en vérifiant pour le n-ième qu'il n'overlap pas les n-1 précédents. Si overlap, recommencer jusqu'à ce qu'il n'est ai plus => attention si le rayon est trop grand et le nombre de points trop grands, ça va boucler indéfiniement, il faut donc choisir le rayon en fonction du nombre de points ou inversement pour avoir une densité relativement faible

Pour avoir une idée, jette un coup d'oeil au code ci-dessous, ça fait le job pour des sphères dures sans conditions de bords. Ce n'est pas forcément la ou une bonne solution, mais l'algo est détaillé :

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
program hs
  implicit none
  integer, parameter          :: n     = 10000
  double precision, parameter :: rad   = 0.50d0/sqrt(dble(n))
  integer, parameter          :: nstop = 10000
  integer            :: i, j, nrejet ! iterateurs
  logical            :: accept       ! condition d'acceptance du tirage
  double precision   :: x(n), y(n)   ! les mailles du reseau
  double precision   :: d2, a, b     ! variables intermediaires pour le tirage
  double precision   :: d2min        ! distance min entre 2 points au carre
  !
  d2min = (2.d0*rad)**2 ! rad est le rayon d'une sphère, d2min distance minimale au carré
  ! Le premier point... car il en faut bien un ;-) !
  call random_number(a)
  call random_number(b)
  x(1)=a
  y(1)=b
  ! Les autres points
  do i=2,n
    nrejet = 0               ! On initialise un compteur quelconque ...
                             !    ... pour ne pas boucler indéfiniement
lm: do                       ! On commence un boucle pour tirer un point valide
      call random_number(a)  ! On tire x ...
      call random_number(b)  !           et y ...
      accept = .true.        !                et on suppose que c'est ok
lt:   do j=1,i-1             ! ET ON TESTE SUR TOUS LES POINTS DEJA VALIDES
        d2 = (x(j)-a)**2 + (y(j)-b)**2
        if ( d2 < d2min ) then
          accept = .false.   ! si on a un et un seul overlap, c'est mauvais !
          exit lt            !    => pas la peine de tester les autres points
        endif
      enddo lt
      if ( accept ) exit lm  ! si on a accept qui est encore true ici : YES !
                             !          ... et on peut sortir de la boucle lm
      nrejet = nrejet+1      ! sinon, on a un rejet supplémentaire et on boucle...
                             ! pas trop quand meme...
      if ( nrejet > nstop ) stop 'trop de rejets successifs...'
    enddo lm
    ! si on en est arrive la sans STOP, on a notre point 'i' 
    x(i) = a
    y(i) = b
  enddo
end program hs

Il faut aussi faire attention aux effets de bords s'il y en a.

**BAmal** · 04/11/2013, 17h15

SALUT,

Un grand merci pour votre intérêt,j essaie de comprendre votre code ,je vois à peu prés ce que vous voulez faire,je vais regarder encore pour en déduire peut être ce que je veux vraiment obtenir!
en pièce joint les nœuds de maillage que je voulais créer ,des nœuds aléatoires ouiiiiiiii mais un maillage homogène c est ce que j arrive pas a faire!!!!

merci

**xflr6** · 04/11/2013, 19h09

Bonsoir,

Il me faut plus d'explications pour comprendre.

Sur la pièce jointe :
Est-ce ce que vous voulez pour les points noirs ?
Il y a me semble-t-il un traitement spécial pour les bords ?
C'est quoi le truc rouge ?

Pour le reste :
N'y a-t-il pas contradiction conceptuelle entre aléatoire et homogène ?
L’homogénéité arriverait avec passage à la limite (rayon -> 0 et N -> \infty).
Ou alors, le rayon, c'est pour "gagner" en "homogénéité" (ou plutot pour éviter la clusterisation) ?

**BAmal** · 05/11/2013, 09h39

Je vais faire de mon mieux pour vous expliquer ;

Alors:
1- Sur la pièce jointe que je vous ai envoyé la dernière fois ---- le truc en rouge n'est pas important dans cette étape ,j'en aurai besoin après ------- ,ce sont plutôt les points noirs ,que je veux créer --------plus précisément ces points noirs sont bien les nœuds de mon maillage ,ils ont été crées de manière aléatoire mais le maillage dans son ensemble est homogène ,sur ce que je viens de vous envoyer maintenant le même maillage c'est à dire des points noirs aléatoires pour le même nombre de nœuds, avant de passer à l'étape homogène c est à dire mes nœuds avant que j'impose une certaine distance minimale) ---------
2-Oui , vous avez raison ,les nœuds aux bords ,je les ai crées séparément ,mais il leurs faut aussi une condition pour qu'ils ne soient pas aussi non homogène c est à dire (en clustérisation si je peux dire ça)
3-Oui, y a une contradiction conceptuelle entre aléatoire et homogène dans le sens propre de mot ,mais ici, ce qu on veux bien c est de ne pas avoir une clusterisation qui génère bien des matrices singulières dans notre calcul parce que par exemple ,en utilisant une maniéré aléatoire on peut facilement générer deux points qui sont trop proches l'une de l'autre ce qui complique nos mesures après.
4-Pourquoi donc ne pas utiliser un maillage régulier? ,oui,le maillage régulier je l'ai fait et j ai validé mes calcul avec,ce que je veux bien c est d avoir plus de précision dans mes calculs et cela est fait par le maillage irrégulier (aléatoire /homogène) ,ce qui a été déjà prouvé par des chercheurs dans le même domaine.

j’espère que je me fais bien comprendre.
Merci encore une fois pour votre intérêt.

**xflr6** · 05/11/2013, 10h56

Ok, mais alors, modulo éventuellement les points aux bords, qu'est-ce qui ne va pas et que vous voulez changer dans la figure que vous avez mise en PJ hier ?

**LittleOwl** · 05/11/2013, 10h57

Bonjour,

Je ne suis pas sur d'avoir absolument bien saisi ce que vous vouliez faire, mais j'ai l'impression que vous souhaitez générer un maillage cartésien irrégulier de manière aléatoire tel que la distance entre 2 noeuds quelconques soit supérieure à une distance minimale.

Si c'est le cas, l'hypothèse d'un maillage cartésien simplifie considérablement le problème car, pour que la dernière condition soit vérifiée, il faut et il suffit que le pas du réseaux selon les abscisses (resp. les ordonnées) soit supérieur à la distance minimale que l'on souhaite imposer entre les noeuds.

Le code ci-dessous implante cette approche.

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
 
program random_irregular_cartesian_grid
implicit none
real, parameter		::  dmin = 2., &	! Distance minimale entre 2 noeuds (unite arbitraire)
				     dmax = 3.	! Distance maximale entre 2 noeuds (unite arbitraire)
 
integer, parameter	:: n  = 100	! Nbr de points de la grille
 
real, dimension(n)   :: x, &	! grille (abscisses)
			       y	! grille (ordonnees)
real, dimension(2*n-2) :: rnd	! Nombres pseudo-aléatoires 0 <= rnd <= 1
 
integer :: i, j, k, l
real	:: dx, dy, d
real	:: time_begin,time_end
 
!**** Generation de la grille
 
call cpu_time(time_begin)
call random_number(rnd)	
x(1) = 0.; y(1) = 0.	! Le premier point est choisi comme origine
j = 1				! Indice permettant de parcourir le tableau rnd
do i = 2, n
  dx = rnd(j)*dmin   + (1.-rnd(j))*dmax;	x(i) = x(i-1) + dx  ! dmin <= dx <= dmax
  dy = rnd(j+1)*dmin + (1.-rnd(j+1))*dmax;	y(i) = y(i-1) + dy  ! dmin <= dy <= dmax
  j = j + 2
end do
call cpu_time(time_end)
 
print *,"*** Grille (Pt n°|x|y)",n*n," points."
do i = 1, n
	do j = 1, n
	 print *,j+n*(i-1),x(i),y(j)
	end do
end do
print '("Generation",t20,": ",es7.1," s.")',(time_end-time_begin)
 
!*** Verification de la condition ||M_kl.M_ij|| >= 
 
call cpu_time(time_begin)
do i = 1, n
 do j = 1, n
   do k = 1, n
     do l = 1, n
	   dx = x(k) - x(i)
	   dy = y(l) - y(j)
	   d  = sqrt(dx*dx + dy*dy)
	   if ((d < dmin).AND.((k /= i).OR.(l /= j))) then
		print '("Error : Pt N°",i7," - ",i7,"d = ",g9.3)',j+n*(i-1),l+n*(l-1),d
	   end if
	 end do
   end do
 end do
end do 
call cpu_time(time_end)
print '("Verification",t20,": ",es7.1," s.")',(time_end-time_begin)
 
end program random_irregular_cartesian_grid

1. l19 - l26 : Génération du maillage. Le premier point (x(1),y(1)) étant choisi comme origine, il est nécessaire de générer aléatoirement n-1 abscisses et n-1 ordonnées. Pour ce faire, le générateur de nombres peuso-aléatoires est appelé l19 pour générer 2n-2 nombres aléatoires, tous compris entre 0 et 1. Les nombres d'indices pairs (resp. impairs) sont utilisés pour générer un pas aléatoire compris entre dmin et dmax selon l'axe des abscisses (l23) et des ordonnées (l24).

2. l29 - l35 : Affichage de la liste des noeuds. Les noeuds sont numérotés et leurs coordonnées affichées en commençant par l'origine et en suivant l'ordre des lignes.

3. l40 - l53 : Vérification de la condition sur la distance entre les noeuds du réseau. La distance entre 2 noeuds quelconques du réseau est calculée de manière systématique et on vérifie qu'elle est bien supérieure à la distance minimale souhaitée. Le temps de calcul de cette étape peut être réduit lorsque le nombre de noeuds est élevé en se limitant aux premiers voisins.

Ce problème est intermédiaire entre la génération d'un maillage régulier et celui d'un maillage non structuré aléatoire, type sphère dure, tel que le décrit xflr6.

Bonne continuation.

**BAmal** · 05/11/2013, 14h35

à xflr6

la PJ que je vous ai envoyé hier c est ce que j aimerais bien obtenir ,ce sont pas mes résultats

, mes résultats sont sur la pj que je vous ai envoyé aujourd’hui (mon maillage est toujours aléatoires).

à LittleOwl

un grand merci à vous je viens de visualiser le maillage généré par votre code, c'est pas exactement un maillage cartésien irrégulier que je voulais ,mais je vois très bien que le maillage généré peut m aider beaucoup dans mes calculs, un grand merci à vous,cependant, comme je travaille sur une cavité de dimension (1m*1m) ,comment je fais pour imposer cette dimension "l",les dimensions de ma géométrie car sur votre code les coordonnées x(n) et y(n) dépendent de dx et de dy qui sont générés aléatoirement et moi je veux bien des nœuds qui ont des abscisses (des ordonnées ) égales à "l" (des nœuds sur les bords et qui depassent pas la longueur de la géométrie (la largeur) quel que soit la dmin imposé ) les dimensions de (ma géométrie) ,j espère que je me fais comprendre
un grand merci à vous encore une fois

**xflr6** · 05/11/2013, 16h05

Envoyé par BAmal

à xflr6
la PJ que je vous ai envoyé hier c est ce que j aimerais bien obtenir ,ce sont pas mes résultats

, mes résultats sont sur la pj que je vous ai envoyé aujourd’hui (mon maillage est toujours aléatoires).

Très bien, excusez moi, je n'avais pas bien compris.

J'ai 'impression que vous désirez un réseau non structuré vu la PJ d'hier. Je vous mets en PJ des exemples de ce qui est obtenu avec le code que je vous avais donné [modulo les lignes de paramètres

Code :

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

1
2
3
  integer, parameter          :: n     = 10000
  double precision, parameter :: rad   = 0.25d0/sqrt(dble(n))
  integer, parameter          :: nstop = 100000

que j'ai changé (nstop à 1000000, n=1000, et rad varie de 0.2 à 0.8/sqrt(n). Vous avez le script test.ksh qui génère les résultats, le code hardspheres.f90 et le fichier test.pdf avec les résultats (certaines fois il y a trop de rejets -ce qui peut dire que c'est impossible d'en mettre plus, ou alors que c'est trop long de tirer un truc dans le peu de place qui reste-), il n'y a donc pas 1000 points.
- à 0.2/sqrt(n) on a une distrib clusterisée
- à 0.4/sqrt(n) on a une distrib avec 1000 points pas mal
- à 0.6/sqrt(n) et plus on a une distrib bien homogène mais on n'est pas arrivé à mettre 1000 points...
- idem pour 0.8

Pour chaque test, il y a le # de rejet en fonction de 'indice du tirage (le numéro du point), puis la répartition des sphères dures (avec leur rayon cohérent), et la répartition des points pour se rendre compte de la distribution des points.

Cordialement.

PS1: test.ksh est renommé en test.txt pour l'upload
PS2: hardspheres.f90 renommé en hardspheres.txt pour l'upload aussi.

**LittleOwl** · 05/11/2013, 23h45

comme je travaille sur une cavité de dimension (1m*1m) ,comment je fais pour imposer cette dimension "l",les dimensions de ma géométrie car sur votre code les coordonnées x(n) et y(n) dépendent de dx et de dy qui sont générés aléatoirement et moi je veux bien des nœuds qui ont des abscisses (des ordonnées ) égales à "l" (des nœuds sur les bords et qui depassent pas la longueur de la géométrie (la largeur) quel que soit la dmin imposé ) les dimensions de (ma géométrie) ,j espère que je me fais comprendre
un grand merci à vous encore une fois

J'ai modifié le code de manière à ce que :
1 - Vous puissiez définir les dimensions (lx) et (ly) du domaine rectangulaire
2 - Les bords du domaine comportent des noeuds.

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
 
program random_irregular_cartesian_grid_v2
 
implicit none
 
real, parameter		::	lx = 1., &	! Dimension du domaine selon l'axe Ox
					ly = 1.		! Dimension du domaine selon l'axe Oy
 
real, parameter		::  dmin = 0.05	! Distance minimale entre 2 noeuds
 
integer, parameter	::  nx = 10, &	! Nbr de colonnes
				     ny = 15		! Nbr de lignes
 
real, dimension(nx)   :: x	 ! grille (abscisses)
real, dimension(ny)   :: y	 ! grille (ordonnees)
real, dimension(nx+ny-4) :: rnd	! Nombres pseudo-aléatoires 0 <= rnd <= 1
 
integer :: i, j, k, l
integer :: nxmax, nymax
real	:: dx, dy, d, hx, hy
real	:: time_begin,time_end
 
!**** Generation de la grille
 
nxmax = lx/dmin + 1; nymax = ly/dmin + 1
if (nx > nxmax) then
	print *,"Contrainte impossible à satisfaire. Trop de noeuds selon Ox (nx > nxmax)"
	stop
end if
 
if (ny > nxmax) then
	print *,"Contrainte impossible à satisfaire. Trop de noeuds selon Oy (ny > nymax)"
	stop
end if
 
call cpu_time(time_begin)
hx = lx/(nx-1); hy = ly/(ny-1)	! Pas du maillage régulier du domaine
x(1) = 0.; x(nx) = lx	! Les bords inf. et sup. du domaine sont mailles
y(1) = 0.; y(ny) = ly	! Les bords gauche et droit du domaine sont mailles
call random_number(rnd)	
do i = 2, nx-1
  x(i) = (i-1)*hx + (rnd(i-1)-0.5)*(hx-dmin)
end do
 
do i = 2, ny-1
  y(i) = (i-1)*hy + (rnd(nx+i-3)-0.5)*(hy-dmin)
end do
call cpu_time(time_end)
 
print *,"*** Grille (Pt n°|x|y)",nx*ny," points."
do i = 1, nx
	do j = 1, ny
	 print *,j+ny*(i-1),x(i),y(j)
	end do
end do
print '("Generation",t20,": ",es7.1," s.")',(time_end-time_begin)
 
!*** Verification de la condition ||M_kl.M_ij|| >= 
 
call cpu_time(time_begin)
do i = 1, nx
 do j = 1, ny
   do k = 1, nx
     do l = 1, ny
	   dx = x(k) - x(i)
	   dy = y(l) - y(j)
	   d  = sqrt(dx*dx + dy*dy)
	   if ((d < dmin).AND.((k /= i).OR.(l /= j))) then
		print '("Error : Pt N°",i7," - ",i7," d = ",g9.3)',j+ny*(i-1),l+ny*(l-1),d
	   end if
	 end do
   end do
 end do
end do 
call cpu_time(time_end)
print '("Verification",t20,": ",es7.1," s.")',(time_end-time_begin)
 
end program random_irregular_cartesian_grid_v2

L'idée consiste à partir d'un maillage cartésien régulier du domaine, défini par les pas (hx) et (hy) et à appliquer à ces noeuds une petite perturbation, dont l'amplitude est calculée de manière à ce que la distance entre 2 noeuds consécutifs reste toujours supérieure à une distance minimale fixée (dmin).

L'existence de cette distance minimale limite le nombre de noeuds admissibles du maillage. Les tests des l24 - l33 affiche un message d'erreur si le nombre de noeuds selon l'un des axes est incompatible avec la distance minimale définie.

Bonne continuation.

**BAmal** · 06/11/2013, 10h41

Bonjour,

à xflr6
Vraiment

,je trouve pas comment vous remercier, c est vraiment ce que je voulais faire(je mets en pièce jointe les résultats que j ai pu visualisé) ,en revanche

,moi pour mon maillage les nœuds aux bords je les crée séparément
ici sur votre maillage j'arrive pas à trouver une méthode pour pouvoir appliquer la condition à ces nœuds aussi????

sur ce code la méthode que j ai faite pour crée ces nœuds (elle est très simple et sans condition les nœuds sont placés sur les bords aléatoirement )

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
           call random_number(x)
		     call random_number(y)
		     do i=1,n*n 
		     front(i)=100
		     do k=1,n-2	 
		     if (i==k.or.i==n-1) then                            !la frontière y=0
		     front(i)=1
		     y(i)=0		  
		     endif 
		     if (i==k*n.or.i==n.or.i==(n-1)*n) then              !la frontière x=l
		     front(i)=4
		     x(i)=lx		  
		     endif
		     if (i==((k*n)+1).or.i==((n-1)*n+1).or.i==1) then    !la frontière x=0
		     front(i)=2 
		     x(i)=0
		     endif
		     if (i==n*n-k.or.i==n*n) then                        !la frontière y=ly
		     front(i)=3
		     y(i)=ly  
		     endif
		     if (i==n*n) then                                    !la frontière x=lx
		     front(i)=4
		     x(i)=lx
		     y(i)=ly
		     endif

**BAmal** · 06/11/2013, 10h58

à LittleOwl

Un grand merci à vous aussi

,je vais travailler sur les deux maillages celui de xflr6 et le votre,et je pense que je vais avoir de bons résultats,je vous tiens au courant

Bon courage à vous

**xflr6** · 06/11/2013, 18h48

Bonsoir,

Je ne vais pas essayer de chercher à comprendre les lignes de code que tu proposes, pas ce soir. Concernant ta question sur comment faire pour ajouter des points sur les bords, plusieurs questions me viennent à l'esprit :

- est-ce vraiment nécessaire ? car ça brise l'esprit de la méthode... A moins qu'il y ai quelque chose de très spécial sur les bords que vous voudriez imposer.
- combien de points sur les bords, combien de points à l'extérieur ?

**BAmal** · 06/11/2013, 19h24

Bonsoir,

OUI,c'est aux bords (aux nœuds des bords) de notre cavité qu'on impose nos conditions aux limites et pour notre problème, si nous n'avons pas de conditions aux limites nous n'arriverons pas à résoudre l’équation qui décrit le phénomène physique

.

POUR,le nombre de nœuds nécessaire ,par exemple un nombre de int(sqrt(n)) suffit largement pour bien présenter les conditions aux limites ,(31 points par chaque coté de notre cavité carrée----pour un total de n=1000)

Bonne soirée

**xflr6** · 06/11/2013, 19h43

Super.

Ce que j'ai fait (à vérifier car je suis pressé et j'ai pris une solution "Quick&Dirty" avec des copier-coller) :
- je tire les points des bords aléatoirement (y compris l'axe auquel il appartient) en supposant des sphères dures (pour les détails, quick and dirty à la mode 2D alors que c'est 1D, osef)
- sphères dures avec bords durs (c'est à dire pas de points à moins de "rad" des bords) pour les points intérieurs. [bords durs : C'est une contrainte qu'il peut être judicieux de lever]. Les sphères des bords sont prises en comptes bien sûr.

Pour le 31 par côtés, c'est trop pour la distance minimale imposée, j'ai donc 3sqrt(n) points sur les 4 axes. Attention à l'usage trop intuitif des racines...

Je mets le code et le résultat en PJ.
En cas de pb, please take it easy, j'ai fait ça vraiment rapidos.