!*************************************************************************
Autori: 
Paolo Sordi Ordine degli ingegneri di Roma N° 22979 - Sezione A
Leonardo Santini Ordine degli ingegneri di Roma N° 22757 - Sezione A
 
!*************************************************************************
refinement.f90
 
SUBROUTINE RIFINISCI_LA_MESH_BIDIMENSIONALE_CON_CRITERIO_FRATTALE
USE CONFIG_MESH_2D
USE WINTERACTER
USE IDENTIFIERS
IMPLICIT NONE
TYPE (win_message) :: MESSAGE
INTEGER :: ITYPE 
LOGICAL :: continua
INTEGER :: ind,numlati
INTEGER :: It,i
INTEGER :: Iv1,iv2,iv3,im1,im2,im3,im4,im5,im6
INTEGER, DIMENSION (NTRI_2D*3,3) :: elenco_lati
CALL WDIALOGPUTINTEGER (idf_integer2 , NTRI_2D*8)
 
CALL CALCOLO_DEL_FATTORE_DI_SCALA_PER_LE_PROIEZIONI_ORTOGONALI ()
numlati = 0
DO It=1, NTRI_2D
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO_X_BIDIMENSIONALE (NT_2D(it,1),NT_2D(it,2), elenco_lati, numlati)
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO_X_BIDIMENSIONALE (NT_2D(it,1),NT_2D(it,3), elenco_lati, numlati)
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO_X_BIDIMENSIONALE (NT_2D(it,2),NT_2D(it,3), elenco_lati, numlati)
END DO
 
DO It=1, NTRI_2D
    iv1 =  NT_2D(it,1)
    iv2 =  NT_2D(it,2)
    iv3 =  NT_2D(it,3)
    DO i=1,numlati
        IF (elenco_lati (i,1).eq.iv1.and.elenco_lati (i,2).eq.iv2) im1 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv1.and.elenco_lati (i,1).eq.iV2) im1 = elenco_lati (i,3)
 
        IF (elenco_lati (i,1).eq.iv1.and.elenco_lati (i,2).eq.iv3) im2 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv1.and.elenco_lati (i,1).eq.iV3) im2 = elenco_lati (i,3)
 
        IF (elenco_lati (i,1).eq.iv2.and.elenco_lati (i,2).eq.iv3) im3 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv2.and.elenco_lati (i,1).eq.iV3) im3 = elenco_lati (i,3)
        CALL WMESSAGEPEEK (ITYPE, MESSAGE)
        IF (ITYPE.eq.CLOSEREQUEST) STOP
    END DO
    CALL WMESSAGEPEEK (ITYPE, MESSAGE)
    IF (ITYPE.eq.CLOSEREQUEST) STOP
    NT_2D(it,1) =  iv1
    NT_2D(it,2) =  im1
    NT_2D(it,3) =  im2
    NTRI_2D = NTRI_2D + 1
    NT_2D(NTRI_2D,1) =  iv2
    NT_2D(NTRI_2D,2) =  im3
    NT_2D(NTRI_2D,3) =  im1
    NTRI_2D = NTRI_2D + 1
    NT_2D(NTRI_2D,1) =  iv3
    NT_2D(NTRI_2D,2) =  im2
    NT_2D(NTRI_2D,3) =  im3
    NTRI_2D = NTRI_2D + 1
    NT_2D(NTRI_2D,1) =  im1
    NT_2D(NTRI_2D,2) =  im3
    NT_2D(NTRI_2D,3) =  im2
!    CALL WDIALOGPUTINTEGER (idf_integer3 , NTETRA)
!    CALL WDIALOGPUTINTEGER (idf_integer2 , NLAT)
!    CALL WDIALOGPUTINTEGER (idf_integer1 , NNOD)
    CALL WMESSAGEPEEK      (ITYPE, MESSAGE)
    IF (ITYPE.eq.CLOSEREQUEST) STOP
END DO
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE INSERISCI_UN_NODO_IN_MEZZO_AL_LATO_X_BIDIMENSIONALE (in1, in2, elenco_lati, numlati)
USE CONFIG_MESH_2D
IMPLICIT NONE
INTEGER, DIMENSION (NTRI_2D*3,3) :: elenco_lati
INTEGER :: in1,in2,numlati
INTEGER :: i
DO i=1,numlati
    if (elenco_lati (i,1).eq.in1.and.elenco_lati (i,2).eq.in2) RETURN
    if (elenco_lati (i,2).eq.in1.and.elenco_lati (i,1).eq.in2) RETURN
END DO
NNOD_2D = NNOD_2D + 1
X_2D (NNOD_2D) = (X_2D(in1) + X_2D(in2)) / 2.0
Y_2D (NNOD_2D) = (Y_2D(in1) + Y_2D(in2)) / 2.0
Z_2D (NNOD_2D) = (z_2D(in1) + z_2D(in2)) / 2.0
numlati = numlati + 1
elenco_lati (numlati,1) = in1
elenco_lati (numlati,2) = in2
elenco_lati (numlati,3) = NNOD_2D
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE RIFINISCI_LA_MESH_CON_CRITERIO_FRATTALE (i_esci)
USE CONFIG_MESH
USE CONFIG_FEM
USE WINTERACTER
USE IDENTIFIERS
IMPLICIT NONE
TYPE (win_message) :: MESSAGE
INTEGER :: ITYPE,i_esci
LOGICAL :: continua
INTEGER :: ind,numlati
INTEGER :: It,i
INTEGER :: Iv1,iv2,iv3,iv4,im1,im2,im3,im4,im5,im6
 
 
INTEGER, DIMENSION (MAXTET*3,3) :: elenco_lati
IF (i_esci.ne.1) CALL WDIALOGSELECT     (IDD_DIALOG001)
IF (i_esci.ne.1) CALL WDIALOGPUTINTEGER (idf_integer5 , NTETRA*8)
 
CALL CALCOLO_DEL_FATTORE_DI_SCALA_PER_LE_PROIEZIONI_ORTOGONALI ()
numlati = 0
 
CALL WDIALOGRANGEPROGRESSBAR (IDF_PROGRESS1,1,NTETRA)
 
 
DO It = 1, NTETRA
    CALL WDIALOGPUTPROGRESSBAR (IDF_PROGRESS1,it,Absolute)
 
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO (NT(it,1),NT(it,2), elenco_lati, numlati)
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO (NT(it,1),NT(it,3), elenco_lati, numlati)
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO (NT(it,1),NT(it,4), elenco_lati, numlati)
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO (NT(it,2),NT(it,3), elenco_lati, numlati)
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO (NT(it,2),NT(it,4), elenco_lati, numlati)
    CALL INSERISCI_UN_NODO_IN_MEZZO_AL_LATO (NT(it,3),NT(it,4), elenco_lati, numlati)
    CALL WDIALOGPUTINTEGER (idf_integer1, NNOD)
    CALL WMESSAGEPEEK      (ITYPE, MESSAGE)
    IF (ITYPE.eq.CLOSEREQUEST) STOP
END DO
 
IF (i_esci.eq.1) RETURN
 
DO It = 1, NTETRA
    iv1 =  NT(it,1)
    iv2 =  NT(it,2)
    iv3 =  NT(it,3)
    iv4 =  NT(it,4)
    DO i=1,numlati
        IF (elenco_lati (i,1).eq.iv1.and.elenco_lati (i,2).eq.iv2) im1 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv1.and.elenco_lati (i,1).eq.iV2) im1 = elenco_lati (i,3)
 
        IF (elenco_lati (i,1).eq.iv1.and.elenco_lati (i,2).eq.iv3) im2 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv1.and.elenco_lati (i,1).eq.iV3) im2 = elenco_lati (i,3)
 
        IF (elenco_lati (i,1).eq.iv1.and.elenco_lati (i,2).eq.iv4) im3 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv1.and.elenco_lati (i,1).eq.iV4) im3 = elenco_lati (i,3)
 
        IF (elenco_lati (i,1).eq.iv2.and.elenco_lati (i,2).eq.iv3) im4 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv2.and.elenco_lati (i,1).eq.iV3) im4 = elenco_lati (i,3)
 
        IF (elenco_lati (i,1).eq.iv2.and.elenco_lati (i,2).eq.iv4) im5 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv2.and.elenco_lati (i,1).eq.iV4) im5 = elenco_lati (i,3)
 
        IF (elenco_lati (i,1).eq.iv3.and.elenco_lati (i,2).eq.iv4) im6 = elenco_lati (i,3)
        IF (elenco_lati (i,2).eq.iv3.and.elenco_lati (i,1).eq.iV4) im6 = elenco_lati (i,3)
        CALL WMESSAGEPEEK (ITYPE, MESSAGE)
        IF (ITYPE.eq.CLOSEREQUEST) STOP
    END DO
    CALL WMESSAGEPEEK (ITYPE, MESSAGE)
    IF (ITYPE.eq.CLOSEREQUEST) STOP
    NT(it,1) =  iv1
    NT(it,2) =  im1
    NT(it,3) =  im2
    NT(it,4) =  im3
 
    CALL ORDINANO (it)
    NTETRA = NTETRA + 1
    NT(NTETRA,1) =  iv2
    NT(NTETRA,2) =  im1
    NT(NTETRA,3) =  im4
    NT(NTETRA,4) =  im5
    NT_PRIMARIO (NTETRA) =  NT_PRIMARIO (it)
 
 
    CALL ORDINANO (NTETRA)
    MATERIALE_TETRAEDRO (NTETRA) = MATERIALE_TETRAEDRO (it)
    NTETRA = NTETRA + 1
    NT(NTETRA,1) =  iv3
    NT(NTETRA,2) =  im2
    NT(NTETRA,3) =  im4
    NT(NTETRA,4) =  im6
    NT_PRIMARIO (NTETRA) =  NT_PRIMARIO (it)
 
    CALL ORDINANO (NTETRA)
    MATERIALE_TETRAEDRO (NTETRA) = MATERIALE_TETRAEDRO (it)
    NTETRA = NTETRA + 1
    NT(NTETRA,1) =  iv4
    NT(NTETRA,2) =  im3
    NT(NTETRA,3) =  im5
    NT(NTETRA,4) =  im6
    NT_PRIMARIO (NTETRA) = NT_PRIMARIO (it)
    CALL ORDINANO (NTETRA)
    MATERIALE_TETRAEDRO (NTETRA) = MATERIALE_TETRAEDRO (it)
!-------------------------------------------------------------------------
    NTETRA = NTETRA + 1
    NT(NTETRA,1) =  im1
    NT(NTETRA,2) =  im6
    NT(NTETRA,3) =  im2
    NT(NTETRA,4) =  im3
    NT_PRIMARIO (NTETRA) = NT_PRIMARIO (it)
    CALL ORDINANO (NTETRA)
    MATERIALE_TETRAEDRO (NTETRA) = MATERIALE_TETRAEDRO (it)
    NTETRA = NTETRA + 1
    NT(NTETRA,1) =  im1
    NT(NTETRA,2) =  im6
    NT(NTETRA,3) =  im2
    NT(NTETRA,4) =  im4
    NT_PRIMARIO (NTETRA) = NT_PRIMARIO (it)
    CALL ORDINANO (NTETRA)
    MATERIALE_TETRAEDRO (NTETRA) = MATERIALE_TETRAEDRO (it)
    NTETRA = NTETRA + 1
    NT(NTETRA,1) =  im1
    NT(NTETRA,2) =  im6
    NT(NTETRA,3) =  im5
    NT(NTETRA,4) =  im4
    NT_PRIMARIO (NTETRA) = NT_PRIMARIO (it)
    CALL ORDINANO (NTETRA)
    MATERIALE_TETRAEDRO (NTETRA) = MATERIALE_TETRAEDRO (it)
    NTETRA = NTETRA + 1
    NT(NTETRA,1) =  im1
    NT(NTETRA,2) =  im6
    NT(NTETRA,3) =  im3
    NT(NTETRA,4) =  im5
    NT_PRIMARIO (NTETRA) = NT_PRIMARIO (it)
    CALL ORDINANO (NTETRA)
    MATERIALE_TETRAEDRO (NTETRA) = MATERIALE_TETRAEDRO (it)
    CALL WDIALOGPUTINTEGER (idf_integer3 , NTETRA)
    CALL WDIALOGPUTINTEGER (idf_integer2 , NLAT)
    CALL WDIALOGPUTINTEGER (idf_integer1 , NNOD)
    CALL WMESSAGEPEEK      (ITYPE, MESSAGE)
    IF (ITYPE.eq.CLOSEREQUEST) STOP
END DO
CALL DISEGNA_LE_PROIEZIONI_ORTOGONALI_X_CARCIOFO (0)
CALL DISEGNA_LE_PROIEZIONI_ORTOGONALI_DEI_TETRAEDRI ()
CALL WDIALOGHIDE
 
!call SCRITTURA_DELLE_COORDINATE_SU_FILE()
!pause
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE INSERISCI_UN_NODO_IN_MEZZO_AL_LATO (in1, in2, elenco_lati, numlati)
USE CONFIG_FEM
USE CONFIG_MESH
IMPLICIT NONE
INTEGER, DIMENSION (MAXTET*3,3) :: elenco_lati
INTEGER :: in1,in2,numlati
INTEGER :: i
DO i=1,numlati
    if (elenco_lati (i,1).eq.in1.and.elenco_lati (i,2).eq.in2) RETURN
    if (elenco_lati (i,2).eq.in1.and.elenco_lati (i,1).eq.in2) RETURN
END DO
NNOD = NNOD + 1
CONDIZIONE_INIZIALE (NNOD) = -1.0
X (NNOD) = (X(in1) + X(in2)) / 2.0
Y (NNOD) = (Y(in1) + Y(in2)) / 2.0
Z (NNOD) = (z(in1) + z(in2)) / 2.0
numlati = numlati + 1
elenco_lati (numlati,1) = in1
elenco_lati (numlati,2) = in2
elenco_lati (numlati,3) = nnod
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE SPLITTA_IL_TETRAEDRO_SUL_LATO_PIU_LUNGO ()
USE CONFIG_MESH
IMPLICIT NONE
INTEGER :: k_div
 
INTEGER :: k_migliore
INTEGER :: i_best,j_best,it_best
 
INTEGER :: i_tetra
DOUBLE PRECISION :: lunghezza_best
DOUBLE PRECISION :: x1,y1,z1,x2,y2,z2
 
DO i_tetra = 1,NTETRA
   lunghezza_best = 0.0
   i_best = 0
   j_best = 0
   CALL TROVA_LATO_PIU_LUNGO (i_tetra,1,2,lunghezza_best,i_best,j_best,it_best)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra,1,3,lunghezza_best,i_best,j_best,it_best)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra,1,4,lunghezza_best,i_best,j_best,it_best)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra,2,3,lunghezza_best,i_best,j_best,it_best)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra,2,4,lunghezza_best,i_best,j_best,it_best)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra,3,4,lunghezza_best,i_best,j_best,it_best)
   x1 = X(NT(i_tetra,i_best))
   y1 = Y(NT(i_tetra,i_best))
   z1 = Z(NT(i_tetra,i_best))
   x2 = X(NT(i_tetra,j_best))
   y2 = Y(NT(i_tetra,j_best))
   z2 = Z(NT(i_tetra,j_best))
   NNOD = NNOD + 1
   X (NNOD) = x1 + (x2-x1) / 2.0
   Y (NNOD) = y1 + (y2-y1) / 2.0
   Z (NNOD) = z1 + (z2-z1) / 2.0
   CALL TETRAEDIZZAZIONE_PER_NODI (NNOD, 0, 0, NNOD)
END DO
CALL WDIALOGHIDE
END SUBROUTINE
 
 
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE TROVA_I_TETRAEDRI_CHE_CONTENGONO_LO_SPIGOLO (it_array,i_vert_1,i_vert_2,i_vert_opp1,i_vert_opp2,trovato)
USE CONFIG_MESH
USE CONFIG_PARAMETRI
USE IDENTIFIERS
USE WINTERACTER
IMPLICIT NONE
INTEGER,DIMENSION (0:100) :: it_array
LOGICAL :: trovato
INTEGER :: i
INTEGER :: i_vert_1,i_vert_2
INTEGER :: i_vert_opp1,i_vert_opp2
INTEGER :: num_tetraedri
 
TETRA_IN_STRUCTURE = .true.
 
trovato =.false.
it_array = 0
DO i=1,NTETRA
   IF (i_vert_1.EQ.NT(i,1).OR.i_vert_1.EQ.NT(i,2).OR.i_vert_1.EQ.NT(i,3).OR.i_vert_1.EQ.NT(i,4)) THEN
       IF (i_vert_2.EQ.NT(i,1).OR.i_vert_2.EQ.NT(i,2).OR.i_vert_2.EQ.NT(i,3).OR.i_vert_2.EQ.NT(i,4)) THEN
            TETRA_IN_STRUCTURE (i) = .false.
            trovato = .true.
            it_array(1) = i
            IF (NT(i,1).NE.i_vert_1.AND.NT(i,1).NE.i_vert_2) i_vert_opp1 = NT (i,1)
            IF (NT(i,2).NE.i_vert_1.AND.NT(i,2).NE.i_vert_2) i_vert_opp1 = NT (i,2)
            IF (NT(i,3).NE.i_vert_1.AND.NT(i,3).NE.i_vert_2) i_vert_opp1 = NT (i,3)
            IF (NT(i,4).NE.i_vert_1.AND.NT(i,4).NE.i_vert_2) i_vert_opp1 = NT (i,4)
            IF (NT(i,1).NE.i_vert_1.AND.NT(i,1).NE.i_vert_2.AND.NT(i,1).NE.i_vert_opp1) i_vert_opp2 = NT (i,1)
            IF (NT(i,2).NE.i_vert_1.AND.NT(i,2).NE.i_vert_2.AND.NT(i,2).NE.i_vert_opp1) i_vert_opp2 = NT (i,2)
            IF (NT(i,3).NE.i_vert_1.AND.NT(i,3).NE.i_vert_2.AND.NT(i,3).NE.i_vert_opp1) i_vert_opp2 = NT (i,3)
            IF (NT(i,4).NE.i_vert_1.AND.NT(i,4).NE.i_vert_2.AND.NT(i,4).NE.i_vert_opp1) i_vert_opp2 = NT (i,4)
            GOTO 100
       END IF
   END IF
END DO
100 CONTINUE
num_tetraedri = 1
DO i=1,NTETRA
   IF (i_vert_1.EQ.NT(i,1).OR.i_vert_1.EQ.NT(i,2).OR.i_vert_1.EQ.NT(i,3).OR.i_vert_1.EQ.NT(i,4)) THEN
       IF (i_vert_2.EQ.NT(i,1).OR.i_vert_2.EQ.NT(i,2).OR.i_vert_2.EQ.NT(i,3).OR.i_vert_2.EQ.NT(i,4)) THEN
           IF (i.ne.it_array(1)) THEN
                TETRA_IN_STRUCTURE (i) = .false.
                num_tetraedri = num_tetraedri + 1
                it_array(0) = num_tetraedri
                it_array(num_tetraedri) = i
           END IF
       END IF
   END IF
END DO
!IF (num_tetraedri.gt.3) THEN
!    write (1000,*) 'Numero di tetredri per lo spigolo .ne. 3',num_tetraedri
!    CALL WDIALOGPUTINTEGER (idf_integer37,num_tetraedri)
!END IF
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE TROVA_I_DUE_TETRAEDRI_CHE_CONTENGANO_LA_FACCIA (i_tetra1, i_tetra2, i_vert_1, i_vert_2, i_vert_3, i_vert_opp, trovato)
USE CONFIG_MESH
USE CONFIG_PARAMETRI
IMPLICIT NONE
LOGICAL :: trovato
INTEGER :: i
INTEGER :: i_tetra1,i_tetra2
INTEGER :: i_vert_1,i_vert_2,i_vert_3
INTEGER :: i_vert_opp
trovato =.false.
i_tetra1= 0
i_tetra2= 0
TETRA_IN_STRUCTURE = .true.
DO i=1,NTETRA
   IF (i_vert_1.EQ.NT(i,1).OR.i_vert_1.EQ.NT(i,2).OR.i_vert_1.EQ.NT(i,3).OR.i_vert_1.EQ.NT(i,4)) THEN
       IF (i_vert_2.EQ.NT(i,1).OR.i_vert_2.EQ.NT(i,2).OR.i_vert_2.EQ.NT(i,3).OR.i_vert_2.EQ.NT(i,4)) THEN
           IF (i_vert_3.EQ.NT(i,1).OR.i_vert_3.EQ.NT(i,2).OR.i_vert_3.EQ.NT(i,3).OR.i_vert_3.EQ.NT(i,4)) THEN
               trovato = .true.
               i_tetra1 = i
               TETRA_IN_STRUCTURE (i) = .false.
               IF (NT(i_tetra1,1).NE.i_vert_1.AND.NT(i_tetra1,1).NE.i_vert_2.AND.NT(i_tetra1,1).NE.i_vert_3)  i_vert_opp = NT (i_tetra1,1)
               IF (NT(i_tetra1,2).NE.i_vert_1.AND.NT(i_tetra1,2).NE.i_vert_2.AND.NT(i_tetra1,2).NE.i_vert_3)  i_vert_opp = NT (i_tetra1,2)
               IF (NT(i_tetra1,3).NE.i_vert_1.AND.NT(i_tetra1,3).NE.i_vert_2.AND.NT(i_tetra1,3).NE.i_vert_3)  i_vert_opp = NT (i_tetra1,3)
               IF (NT(i_tetra1,4).NE.i_vert_1.AND.NT(i_tetra1,4).NE.i_vert_2.AND.NT(i_tetra1,4).NE.i_vert_3)  i_vert_opp = NT (i_tetra1,4)
               GOTO 100
           END IF
       END IF
   END IF
END DO
100 CONTINUE
DO i=1,NTETRA
   IF (i_vert_1.EQ.NT(i,1).OR.i_vert_1.EQ.NT(i,2).OR.i_vert_1.EQ.NT(i,3).OR.i_vert_1.EQ.NT(i,4)) THEN
       IF (i_vert_2.EQ.NT(i,1).OR.i_vert_2.EQ.NT(i,2).OR.i_vert_2.EQ.NT(i,3).OR.i_vert_2.EQ.NT(i,4)) THEN
           IF (i_vert_3.EQ.NT(i,1).OR.i_vert_3.EQ.NT(i,2).OR.i_vert_3.EQ.NT(i,3).OR.i_vert_3.EQ.NT(i,4)) THEN
               TETRA_IN_STRUCTURE (i) = .false.
               IF (i.ne.i_tetra1) i_tetra2 = i
           END IF
       END IF
   END IF
END DO
END SUBROUTINE 
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE  CONTROLLA_NODO_PRIMARIO_PIU_VICINO (xs,ys,zs,distanza_minima,num_nodo)
USE CONFIG_MESH
IMPLICIT NONE
DOUBLE PRECISION :: xs,ys,zs
DOUBLE PRECISION :: distanza_minima
DOUBLE PRECISION :: lunghezza_lato
INTEGER :: i
INTEGER :: num_nodo
num_nodo = 1
distanza_minima = DSQRT(((X(1)-Xs)**2.0 + (Y(1)-Ys)**2.0 + (Z(1)-Zs)**2.0))
DO i=1 , NNOD_PRIMARI
    lunghezza_lato = DSQRT(((X(i)-Xs)**2.0 + (Y(i)-Ys)**2.0 + (Z(i)-Zs)**2.0))
    IF (lunghezza_lato.le.distanza_minima) THEN
        distanza_minima = lunghezza_lato
        num_nodo = i
    END IF
END DO
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE SPOSTA_NODI_SECONDARI_SUI_NODI_PRIMARI
USE CONFIG_MESH
USE CONFIG_PARAMETRI
IMPLICIT NONE
DOUBLE PRECISION :: lato_minimo_presente
DOUBLE PRECISION :: distanza
INTEGER          :: num_nodo
INTEGER          :: in
lato_minimo_presente = DMax1 (XMAX-XMIN,YMAX-YMIN,ZMAX-ZMIN) / FATTORE_REFINIMENT
WRITE (1000,*) 'Lato_Minimo_Presente',lato_minimo_presente
WRITE (1000,*) 'NNOD_PRIMARI ',NNOD_PRIMARI
 
DO In = NNOD_PRIMARI + 1  ,NNOD
    CALL CONTROLLA_NODO_PRIMARIO_PIU_VICINO (x(IN),y(IN),z(IN),distanza,num_nodo)
    WRITE (1000,*) 'distanza',distanza
    IF (distanza.LT.(lato_minimo_presente/15.0)) THEN
         WRITE (1000,*) 'Scambio Effettuato'
         X (in) = x (num_nodo)
         Y (in) = y (num_nodo)
         Z (in) = z (num_nodo)
    END IF
END DO
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE MIGLIORA_QUALITA_MUOVENDO_LATI_SECONDARI
USE CONFIG_MESH
USE CONFIG_PARAMETRI
USE WINTERACTER
IMPLICIT NONE
TYPE (win_message) :: MESSAGE
INTEGER :: itype
INTEGER :: I,J,k
DOUBLE PRECISION :: x_aux,y_aux,z_aux
DOUBLE PRECISION :: qualita
DOUBLE PRECISION :: volume_minimo_per_cento,volume_parallalelepipedo
DOUBLE PRECISION :: volume_minimo_old
INTEGER          :: i_minimo
CALL TROVA_TETRAEDRO_VOLUME_MINIMO(i_minimo,volume_minimo_per_cento,volume_parallalelepipedo)
volume_minimo_old = volume_minimo_per_cento
WRITE (1000,*) 'Nodi primari', NNOD_PRIMARI
WRITE (1000,*) 'Nodi        ', NNOD
DO i = NNOD_PRIMARI + 1, NNOD
    CALL STAMPA_INFORMAZIONI_SU_SCHERMO(0,0)
    CALL CALCOLO_COEFFICIENTE_QUALITA_GLOBALE ()
    QUALITA = QUALITA_MEDIA
    IF (X(i).eq.XMIN.OR.X(i).eq.XMAX) GOTO 2000
    IF (Y(i).eq.YMIN.OR.Y(i).eq.YMAX) GOTO 2000
    IF (Z(i).eq.ZMIN.OR.Z(i).eq.ZMAX) GOTO 2000
    x_aux = x(i)
    y_aux = y(i)
    z_aux = z(i)
    DO j=1,NTETRA
         CALL WMESSAGEPEEK (ITYPE, MESSAGE)
         IF (ITYPE.eq.CLOSEREQUEST) STOP
         IF (i.EQ.NT(J,1).OR.i.EQ.NT(J,2).OR.i.EQ.NT(J,3).OR.i.EQ.NT(J,4)) THEN
              DO k = 1,4
                   IF (i.NE.NT(j,k)) THEN
                         WRITE (1000,*) 'Numero nodo  ',i
                         X(i) = (X_aux * 8 + X(NT(j,k))) / 9.0
                         Y(i) = (Y_aux * 8 + Y(NT(j,k))) / 9.0
                         Z(i) = (Z_aux * 8 + Z(NT(j,k))) / 9.0
                         WRITE (1000,*) 'x_aux, y_aux, z_aux  ',x_aux,y_aux,z_aux
                         WRITE (1000,*) 'nt                   ',X(NT(j,k)),Y(NT(j,k)),Z(NT(j,k))
                         WRITE (1000,*) 'x(i) y(i) z(i)       ',X(i) , Y(i) , Z(i)
                         CALL TROVA_TETRAEDRO_VOLUME_MINIMO(i_minimo,volume_minimo_per_cento,volume_parallalelepipedo)
                         WRITE (1000,*) 'Volume minimo per cento', volume_minimo_per_cento
                         WRITE (1000,*) 'volume minimo old      ', volume_minimo_old
                         IF (volume_minimo_per_cento.lt.volume_minimo_old) THEN
                               x(i) = x_aux
                               y(i) = y_aux
                               z(i) = z_aux
                               WRITE (1000,*) '----- Volume minimo peggiorato -----'
                               GOTO 100
                         ELSE
                               volume_minimo_old = volume_minimo_per_cento
                         END IF
                         CALL CALCOLO_COEFFICIENTE_QUALITA_GLOBALE ()
                         WRITE (1000,*) 'Qualita_media',QUALITA_MEDIA
                         IF (QUALITA_MEDIA.lt.QUALITA) THEN
                               x(i) = x_aux
                               y(i) = y_aux
                               z(i) = z_aux
                         ELSE
                               QUALITA = QUALITA_MEDIA
                         ENDIF
                         100 CONTINUE
                   END IF
              END DO
         END IF
    END DO
2000 CONTINUE
END DO
CALL WDIALOGHIDE
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE  TROVA_LATO_PIU_LUNGO_CON_REGOLE (it,i,j,lunghezza_best,lunghezza_worse,i_best,j_best,it_best)
USE CONFIG_MESH 
IMPLICIT NONE
DOUBLE PRECISION :: lunghezza_best,lunghezza_lato
DOUBLE PRECISION :: lunghezza_worse
INTEGER :: it,i,j,i_best,j_best,it_best
lunghezza_lato = DSQRT(((X(NT(It,i))-X(NT(It,j)))**2.0 + (Y(NT(It,i))-Y(NT(It,j)))**2.0 + (Z(NT(It,i))-Z(NT(It,j)))**2.0))
IF (lunghezza_lato.gt.lunghezza_best.and.lunghezza_lato.lt.lunghezza_worse-1.E-2)  then
   lunghezza_best=lunghezza_lato
   i_best = i
   j_best = j
   it_best = it
END IF
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE  TROVA_LATO_PIU_LUNGO (it,i,j,lunghezza_best,i_best,j_best,it_best)
USE CONFIG_MESH 
IMPLICIT NONE
DOUBLE PRECISION ::   lunghezza_best,lunghezza_lato
INTEGER :: it,i,j,i_best,j_best,it_best
lunghezza_lato = DSQRT(((X(NT(It,i))-X(NT(It,j)))**2.0 + (Y(NT(It,i))-Y(NT(It,j)))**2.0 + (Z(NT(It,i))-Z(NT(It,j)))**2.0))
IF (lunghezza_lato.gt.lunghezza_best)  then
   lunghezza_best=lunghezza_lato
   i_best = i
   j_best = j
   it_best = it
END IF
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE RIFINISCI_LA_MESH_AL_MASSIMO
USE CONFIG_MESH
USE CONFIG_PARAMETRI
USE CONFIG_FEM
USE CONFIG_ZERI
USE WINTERACTER
USE IDENTIFIERS
IMPLICIT NONE
TYPE (win_message) :: MESSAGE
INTEGER :: ITYPE 
LOGICAL :: continua
INTEGER :: ind
INTEGER :: k_div
INTEGER :: It
INTEGER :: k_migliore
INTEGER :: i_best,j_best,it_best
DOUBLE PRECISION :: lunghezza_best
DOUBLE PRECISION :: lato_minimo_parall
DOUBLE PRECISION :: x1,y1,z1,x2,y2,z2
CALL WDIALOGSELECT     (IDD_DIALOG001)
CALL CALCOLO_DEL_FATTORE_DI_SCALA_PER_LE_PROIEZIONI_ORTOGONALI ()
CALL WMESSAGEBOX (0,0,1,'Questa subroutine deve essere stoppata','')
lato_minimo_parall = DMax1 (XMAX-XMIN,YMAX-YMIN,ZMAX-ZMIN)
WRITE (1000,*) 'lato_minimo_parall',lato_minimo_parall
lunghezza_best = lato_minimo_parall
FATTORE_REFINIMENT = FATTORE_REFINIMENT + 0.5
WRITE (1000,*) 'FATTORE_REFINIMENT',FATTORE_REFINIMENT
continua = .TRUE.
DO WHILE (CONTINUA)
continua =.FALSE.
   DO It=1, NTETRA
!      IF (MATERIALE_TETRAEDRO (it).EQ.1) goto 1000
 
       lunghezza_best = 0.0
       CALL TROVA_LATO_PIU_LUNGO (it,1,2,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,1,3,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,1,4,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,2,3,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,2,4,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,3,4,lunghezza_best,i_best,j_best,it_best)
!       WRITE (1000,*) 'lunghezza_best',lunghezza_best
       IF (lunghezza_best.GE.(lato_minimo_parall)/FATTORE_REFINIMENT) THEN
           continua = .TRUE.
           x1 = X(NT(it,i_best))
           y1 = Y(NT(it,i_best))
           z1 = Z(NT(it,i_best))
           x2 = X(NT(it,j_best))
           y2 = Y(NT(it,j_best))
           z2 = Z(NT(it,j_best))
           k_div = 80
           k_migliore = k_div / 2
           NNOD = NNOD + 1
           X (NNOD) = x1 + (x2-x1) * k_migliore / k_div
           Y (NNOD) = y1 + (y2-y1) * k_migliore / k_div
           Z (NNOD) = z1 + (z2-z1) * k_migliore / k_div
           CALL POSIZIONE_PUNTO_IE_RISPETTO_TETRAEDRO (it, X(NNOD), y(NNOD), Z(NNOD), ind, ZERO_NODO_IN_TETRAEDRO)
           CALL IMMISSIONE_NODO_NEL_TETRAEDRO (it,NNOD,ind)
           CONDIZIONE_INIZIALE (NNOD) = -1.0
           CALL DISEGNA_LE_PROIEZIONI_ORTOGONALI_X_CARCIOFO (0)
           CALL DISEGNA_LE_PROIEZIONI_ORTOGONALI_DEI_TETRAEDRI ()
           CALL WDIALOGPUTINTEGER (idf_integer3 , NTETRA)
           CALL WDIALOGPUTINTEGER (idf_integer2 , NLAT)
           CALL WDIALOGPUTINTEGER (idf_integer1 , NNOD)
           CALL WDIALOGPUTREAL    (IDF_REAL8, REAL (FATTORE_REFINIMENT))
           CALL WMESSAGEPEEK      (ITYPE, MESSAGE)
           IF (ITYPE.eq.CLOSEREQUEST) STOP
       END IF
       1000 CONTINUE
   END DO
END DO
PAUSE
CALL WDIALOGHIDE
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE CALCOLA_QUALITA_PEGGIORE_DEL_TETRAEDRO (qualita_peggiore)
USE CONFIG_MESH
IMPLICIT NONE
INTEGER :: I
INTEGER :: It
INTEGER :: i_best,j_best,k_best
INTEGER :: i_tetra_peggiore
DOUBLE PRECISION :: x1,y1,z1,x2,y2,z2,x3,y3,z3
DOUBLE PRECISION :: area_best
DOUBLE PRECISION :: qualita_peggiore
qualita_peggiore = 1.0
DO It = 1, NTETRA
    CALL TROVA_TETRA_CON_QUALITA_PEGGIORE (it,qualita_peggiore,i_tetra_peggiore)
END DO
WRITE (1000,*) 'qualita_peggiore',i_tetra_peggiore,qualita_peggiore
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE  TROVA_TETRA_CON_QUALITA_PEGGIORE (i,qualita_peggiore,itetra_peggiore)
USE CONFIG_MESH
USE CONFIG_PARAMETRI
IMPLICIT NONE
INTEGER :: i
DOUBLE PRECISION :: qt
DOUBLE PRECISION :: qualita_tetra
DOUBLE PRECISION :: qualita_peggiore
INTEGER :: itetra_peggiore
qualita_tetra = QT (I)
IF (qualita_tetra.LE.qualita_peggiore) THEN
     qualita_peggiore = qualita_tetra
     itetra_peggiore = I
END IF
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
SUBROUTINE RIFINISCI_LA_MESH_AL_MASSIMO_IN_BASE_ALLA_QUALITA
USE CONFIG_MESH
USE CONFIG_PARAMETRI
IMPLICIT NONE
INTEGER :: I,k,k_div
INTEGER :: It
INTEGER :: k_migliore
INTEGER :: i_best,j_best
INTEGER :: i_tetra_peggiore
DOUBLE PRECISION :: lunghezza_best
DOUBLE PRECISION :: x1,y1,z1,x2,y2,z2
DOUBLE PRECISION :: xs,ys,zs
DOUBLE PRECISION :: distanza_migliore,distanza_minima
DOUBLE PRECISION :: qualita_peggiore
lunghezza_best = yMAX-yMIN
FATTORE_REFINIMENT = FATTORE_REFINIMENT + 0.2
DO I= 1, 10
   qualita_peggiore = 1.0
   DO It = 1, NTETRA
        CALL TROVA_TETRA_CON_QUALITA_PEGGIORE (it,qualita_peggiore,i_tetra_peggiore)
   END DO
   WRITE (1000,*) 'qualita_peggiore',i_tetra_peggiore,qualita_peggiore
   lunghezza_best = 0.0
   CALL TROVA_LATO_PIU_LUNGO (i_tetra_peggiore,1,2,lunghezza_best,i_best,j_best,i_tetra_peggiore)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra_peggiore,1,3,lunghezza_best,i_best,j_best,i_tetra_peggiore)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra_peggiore,1,4,lunghezza_best,i_best,j_best,i_tetra_peggiore)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra_peggiore,2,3,lunghezza_best,i_best,j_best,i_tetra_peggiore)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra_peggiore,2,4,lunghezza_best,i_best,j_best,i_tetra_peggiore)
   CALL TROVA_LATO_PIU_LUNGO (i_tetra_peggiore,3,4,lunghezza_best,i_best,j_best,i_tetra_peggiore)
   x1 = X(NT(i_tetra_peggiore,i_best))
   y1 = Y(NT(i_tetra_peggiore,i_best))
   z1 = Z(NT(i_tetra_peggiore,i_best))
   x2 = X(NT(i_tetra_peggiore,j_best))
   y2 = Y(NT(i_tetra_peggiore,j_best))
   z2 = Z(NT(i_tetra_peggiore,j_best))
   distanza_migliore = 0.0
   k_div = 40
   DO k= 18, k_div-18
       xs = x1 + (x2-x1) * k / k_div
       ys = y1 + (y2-y1) * k / k_div
       zs = z1 + (z2-z1) * k / k_div
       CALL CONTROLLO_DISTANZA_MINIMA_CON_ALTRI_NODI (xs,ys,zs,distanza_minima)
       IF (distanza_minima.gt.distanza_migliore) THEN
             distanza_migliore = distanza_minima
             k_migliore = k
       END IF 
   END DO
   NNOD = NNOD + 1
   k_migliore=20
   X (NNOD) = x1 + (x2-x1) * k_migliore / k_div
   Y (NNOD) = y1 + (y2-y1) * k_migliore / k_div
   Z (NNOD) = z1 + (z2-z1) * k_migliore / k_div
   CALL TETRAEDIZZAZIONE_PER_NODI (NNOD, 0, 0, NNOD)
END DO 
CALL WDIALOGHIDE
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE RIFINISCI_LA_MESH_AL_MASSIMO_IN_BASE_ALLA_QUALITA_CON_BARICENTRO
USE CONFIG_MESH
USE CONFIG_PARAMETRI
IMPLICIT NONE
INTEGER :: I
INTEGER :: It
INTEGER :: i_tetra_peggiore
DOUBLE PRECISION :: x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4
DOUBLE PRECISION :: qualita_peggiore
FATTORE_REFINIMENT = FATTORE_REFINIMENT + 0.2
DO I= 1, 10
   qualita_peggiore = 1.0
   DO It = 1, NTETRA
        CALL TROVA_TETRA_CON_QUALITA_PEGGIORE (it,qualita_peggiore,i_tetra_peggiore)
   END DO
   WRITE (1000,*) 'qualita_peggiore',i_tetra_peggiore,qualita_peggiore
   x1 = X(NT(i_tetra_peggiore,1))
   y1 = Y(NT(i_tetra_peggiore,1))
   z1 = Z(NT(i_tetra_peggiore,1))
   x2 = X(NT(i_tetra_peggiore,2))
   y2 = Y(NT(i_tetra_peggiore,2))
   z2 = Z(NT(i_tetra_peggiore,2))
   x3 = X(NT(i_tetra_peggiore,3))
   y3 = Y(NT(i_tetra_peggiore,3))
   z3 = Z(NT(i_tetra_peggiore,3))
   x4 = X(NT(i_tetra_peggiore,4))
   y4 = Y(NT(i_tetra_peggiore,4))
   z4 = Z(NT(i_tetra_peggiore,4))
   NNOD = NNOD + 1
   X (NNOD) = (x1 + x2 + x3 + x4) / 4.0
   Y (NNOD) = (y1 + y2 + y3 + y4) / 4.0
   Z (NNOD) = (z1 + z2 + z3 + z4) / 4.0
   CALL TETRAEDIZZAZIONE_PER_NODI (NNOD, 0, 0, NNOD)
END DO 
CALL WDIALOGHIDE
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
SUBROUTINE RIFINISCI_LA_MESH_AL_MASSIMO_IN_BASE_ALLA_QUALITA_CON_AREA
USE CONFIG_MESH
IMPLICIT NONE
INTEGER :: I
INTEGER :: It
INTEGER :: i_best,j_best,k_best
INTEGER :: i_tetra_peggiore
DOUBLE PRECISION :: x1,y1,z1,x2,y2,z2,x3,y3,z3
DOUBLE PRECISION :: area_best
DOUBLE PRECISION :: qualita_peggiore
DO I= 1, 10
   qualita_peggiore = 1.0
   DO It = 1, NTETRA
        CALL TROVA_TETRA_CON_QUALITA_PEGGIORE (it,qualita_peggiore,i_tetra_peggiore)
   END DO
   WRITE (1000,*) 'qualita_peggiore',i_tetra_peggiore,qualita_peggiore
   area_best = 0.0
   CALL TROVA_AREA_PIU_GRANDE (i_tetra_peggiore,1,2,3,area_best,i_best,j_best,k_best)
   CALL TROVA_AREA_PIU_GRANDE (i_tetra_peggiore,1,2,4,area_best,i_best,j_best,k_best)
   CALL TROVA_AREA_PIU_GRANDE (i_tetra_peggiore,2,3,4,area_best,i_best,j_best,k_best)
   CALL TROVA_AREA_PIU_GRANDE (i_tetra_peggiore,1,3,4,area_best,i_best,j_best,k_best)
   x1 = X(NT(i_tetra_peggiore,i_best))
   y1 = Y(NT(i_tetra_peggiore,i_best))
   z1 = Z(NT(i_tetra_peggiore,i_best))
   x2 = X(NT(i_tetra_peggiore,j_best))
   y2 = Y(NT(i_tetra_peggiore,j_best))
   z2 = Z(NT(i_tetra_peggiore,j_best))
   x3 = X(NT(i_tetra_peggiore,k_best))
   y3 = Y(NT(i_tetra_peggiore,k_best))
   z3 = Z(NT(i_tetra_peggiore,k_best))
   NNOD = NNOD + 1
   X (NNOD) = (x1 + x2 + x3) / 3.0
   Y (NNOD) = (y1 + y2 + y3) / 3.0
   Z (NNOD) = (z1 + z2 + z3) / 3.0
   CALL TETRAEDIZZAZIONE_PER_NODI (NNOD, 0, 0, NNOD)
END DO 
CALL WDIALOGHIDE
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE TROVA_AREA_PIU_GRANDE (i_tetra,i,j,k,area_best,i_best,j_best,k_best)
USE CONFIG_MESH
IMPLICIT NONE
INTEGER :: i,j,k
INTEGER :: i_best,j_best,k_best
INTEGER :: i_tetra
DOUBLE PRECISION :: area_best,area
CALL CALCOLO_AREA_DEL_TRIANGOLO (NT(i_tetra,i), NT(i_tetra,j), NT(i_tetra,k), area)
IF (area.gt.area_best) THEN
    i_best = i
    j_best = j
    k_best = k
    area_best = area
END IF 
END SUBROUTINE 
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE CREA_NODI_SECONDARI_PER_MIGLIORARE
USE CONFIG_MESH
USE CONFIG_PARAMETRI
IMPLICIT NONE
INTEGER :: i,k,k_div
INTEGER :: k_migliore
INTEGER :: It
INTEGER :: i_best,j_best,it_best
DOUBLE PRECISION :: lunghezza_best
DOUBLE PRECISION :: x1, y1, z1, x2, y2, z2
DOUBLE PRECISION :: xs, ys, zs
DOUBLE PRECISION :: distanza_migliore
DOUBLE PRECISION :: distanza_minima
DO I=1,10
   lunghezza_best = 0.0
   DO It=1, NTETRA
       CALL TROVA_LATO_PIU_LUNGO (it,1,2,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,1,3,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,1,4,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,2,3,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,2,4,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,3,4,lunghezza_best,i_best,j_best,it_best)
   END DO
   x1 = X(NT(it_best,i_best))
   y1 = Y(NT(it_best,i_best))
   z1 = Z(NT(it_best,i_best))
   x2 = X(NT(it_best,j_best))
   y2 = Y(NT(it_best,j_best))
   z2 = Z(NT(it_best,j_best))
   distanza_migliore = 0.0
   k_div = 40
   DO k= 14, k_div-14
       xs = x1 + (x2-x1) * k / k_div
       ys = y1 + (y2-y1) * k / k_div
       zs = z1 + (z2-z1) * k / k_div
       CALL CONTROLLO_DISTANZA_MINIMA_CON_ALTRI_NODI (xs,ys,zs,distanza_minima)
       IF (distanza_minima.gt.distanza_migliore) THEN
             distanza_migliore = distanza_minima
             k_migliore = k
       END IF 
   END DO
   NNOD = NNOD + 1
   X (NNOD) = x1 + (x2-x1) * k_migliore / k_div
   Y (NNOD) = y1 + (y2-y1) * k_migliore / k_div
   Z (NNOD) = z1 + (z2-z1) * k_migliore / k_div
   CALL TETRAEDIZZAZIONE_PER_NODI (NNOD, 1, 0, NNOD)
END DO
CALL WDIALOGHIDE
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE  CONTROLLO_DISTANZA_MINIMA_CON_ALTRI_NODI (xs,ys,zs,distanza_minima)
USE CONFIG_MESH
IMPLICIT NONE
DOUBLE PRECISION :: xs,ys,zs
DOUBLE PRECISION :: distanza_minima
DOUBLE PRECISION :: lunghezza_lato
INTEGER :: i
distanza_minima = DSQRT(((X(1)-Xs)**2.0 + (Y(1)-Ys)**2.0 + (Z(1)-Zs)**2.0))
DO i=1,NNOD
    lunghezza_lato = DSQRT(((X(i)-Xs)**2.0 + (Y(i)-Ys)**2.0 + (Z(i)-Zs)**2.0))
    IF (lunghezza_lato.le.distanza_minima) distanza_minima = lunghezza_lato
END DO
DO i = 1, 8
    lunghezza_lato = DSQRT(((X(MaxNODI + i)-Xs)**2.0 + (Y(MaxNODI + i)-Ys)**2.0 + (Z(MaxNODI + i)-Zs)**2.0))
    IF (lunghezza_lato/1.0.le.distanza_minima) distanza_minima = lunghezza_lato
END DO
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE RIFINISCI_LA_MESH_CON_NODI_SECONDARI_LONTANI
USE CONFIG_MESH
USE CONFIG_PARAMETRI
 
IMPLICIT NONE
LOGICAL :: continua
INTEGER :: k_div
INTEGER :: It,K
INTEGER :: k_migliore
INTEGER :: i_best,j_best,it_best
DOUBLE PRECISION :: lunghezza_best
DOUBLE PRECISION :: lato_minimo_parall
DOUBLE PRECISION :: x1,y1,z1,x2,y2,z2
DOUBLE PRECISION :: distanza_migliore,distanza_minima
DOUBLE PRECISION :: xs,ys,zs
 
lato_minimo_parall = DMax1 (XMAX-XMIN,YMAX-YMIN,ZMAX-ZMIN)
WRITE (1000,*) 'lato_minimo_parall',lato_minimo_parall
lunghezza_best = lato_minimo_parall
FATTORE_REFINIMENT = FATTORE_REFINIMENT + 0.3
WRITE (1000,*) 'FATTORE_REFINIMENT',FATTORE_REFINIMENT
continua = .TRUE.
DO WHILE (CONTINUA)
continua =.FALSE.
   DO It=1, NTETRA
       lunghezza_best = 0.0
       CALL TROVA_LATO_PIU_LUNGO (it,1,2,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,1,3,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,1,4,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,2,3,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,2,4,lunghezza_best,i_best,j_best,it_best)
       CALL TROVA_LATO_PIU_LUNGO (it,3,4,lunghezza_best,i_best,j_best,it_best)
       IF (lunghezza_best.GE.(lato_minimo_parall)/FATTORE_REFINIMENT) THEN
           continua = .TRUE.
           x1 = X(NT(it,i_best))
           y1 = Y(NT(it,i_best))
           z1 = Z(NT(it,i_best))
           x2 = X(NT(it,j_best))
           y2 = Y(NT(it,j_best))
           z2 = Z(NT(it,j_best))
           k_div = 80
           k_migliore = k_div / 2
           DO k= 35, k_div-35
               distanza_migliore = 0.0
               xs = x1 + (x2-x1) * k / k_div
               ys = y1 + (y2-y1) * k / k_div
               zs = z1 + (z2-z1) * k / k_div
               CALL CONTROLLO_DISTANZA_MINIMA_CON_ALTRI_NODI (xs,ys,zs,distanza_minima)
               WRITE (1000,*) 'distanza_minima ',distanza_minima
               IF (distanza_minima.gt.distanza_migliore) THEN
                   distanza_migliore = distanza_minima
                   k_migliore = k
               END IF
           END DO
           NNOD = NNOD + 1
           X (NNOD) = x1 + (x2-x1) * k_migliore / k_div
           Y (NNOD) = y1 + (y2-y1) * k_migliore / k_div
           Z (NNOD) = z1 + (z2-z1) * k_migliore / k_div
           CALL TETRAEDIZZAZIONE_PER_NODI (NNOD, 0, 0, NNOD)
       END IF
   END DO
END DO
CALL WDIALOGHIDE
END SUBROUTINE