module metricas

use variaveis

contains

!--------------------------------------------------------------------

subroutine metricas_aproximadas

	raux1 = 0.5d0/d_ksi
	raux2 = 0.5d0/d_eta

	! Nós internos
	do i = 2, Nv2-1
		
		aux1 = (i-1)*Nv1

		do j = 2, Nv1-1

			aux2 = aux1 + j
			aux3 = aux2 - Nv1
			aux4 = aux2 + Nv1
			aux5 = aux2 - 1
			aux6 = aux2 + 1

			x_eta(aux2) = raux2*(xv(aux4) - xv(aux3))
			y_eta(aux2) = raux2*(yv(aux4) - yv(aux3))

			x_ksi(aux2) = raux1*(xv(aux6) - xv(aux5))
			y_ksi(aux2) = raux1*(yv(aux6) - yv(aux5))

			alpha(aux2) = x_eta(aux2)**2 + y_eta(aux2)**2
			beta(aux2)  = x_eta(aux2)*x_ksi(aux2) + y_eta(aux2)*y_ksi(aux2)
			gama(aux2)  = x_ksi(aux2)**2 + y_ksi(aux2)**2

			Jac(aux2)   = 1.0d0 / ( x_ksi(aux2)*y_eta(aux2) - x_eta(aux2)*y_ksi(aux2) )

		end do

	end do

end subroutine metricas_aproximadas

!--------------------------------------------------------------------

subroutine centroides

	! Volumes internos

	raux1 = 1.0d0/3.0d0

	! Área de cada elemento de malha
	Sp = 0

	do i = 1, N2

		aux1 = (i-1)*N1
		aux3 = i*Nf1

		do j = 1, N1

			aux2 = aux1 + j
			aux4 = aux3 + j + 1

			! Centroides - primeiro triângulo
			raux2 = raux1*( xv(vert(aux2,1)) + xv(vert(aux2,2)) + xv(vert(aux2,3)) )
			raux3 = raux1*( yv(vert(aux2,1)) + yv(vert(aux2,2)) + yv(vert(aux2,3)) )

			! Área - primeiro triângulo
			raux4 = ( xv(vert(aux2,1))*yv(vert(aux2,2)) + xv(vert(aux2,2))*yv(vert(aux2,3)) + &
			          xv(vert(aux2,3))*yv(vert(aux2,1)) - xv(vert(aux2,3))*yv(vert(aux2,2)) - &
					  xv(vert(aux2,2))*yv(vert(aux2,1)) - xv(vert(aux2,1))*yv(vert(aux2,3)) ) * 0.5d0

			! Centroides - segundo triangulo
			raux5 = raux1*( xv(vert(aux2,1)) + xv(vert(aux2,3)) + xv(vert(aux2,4)) )
			raux6 = raux1*( yv(vert(aux2,1)) + yv(vert(aux2,3)) + yv(vert(aux2,4)) )

			! Área - segundo triângulo
			raux7 = ( xv(vert(aux2,1))*yv(vert(aux2,3)) + xv(vert(aux2,3))*yv(vert(aux2,4)) + &
			          xv(vert(aux2,4))*yv(vert(aux2,1)) - xv(vert(aux2,4))*yv(vert(aux2,3)) - &
					  xv(vert(aux2,3))*yv(vert(aux2,1)) - xv(vert(aux2,1))*yv(vert(aux2,4)) ) * 0.5d0

			!Área do elemento de área
			Sp(aux2) = raux4 + raux7

			! Posição dos centroides
			raux8 = 1.0d0 / (raux4 + raux7)
			xc(aux2) = (raux2*raux4 + raux5*raux7)*raux8 
			yc(aux2) = (raux3*raux4 + raux6*raux7)*raux8

			! Posição dos centroides (com fictícios)
			xcf(aux4) = xc(aux2)
			ycf(aux4) = yc(aux2)

		end do

	end do

	! Contorno sul
	aux1  = Nf1/2
	raux1 = 1.0d0 / (aux1 - 1.0d0)
	do j = 2, aux1
	
		raux2 = (j-1.5d0)*raux1

		xcf(j) = 0.0d0
		ycf(j) = 1.0d0 - raux2

	end do

	do j = aux1 + 1, Nf1-1

		raux2 = ((j-aux1)-0.5d0)*raux1

		xcf(j) = raux2
		ycf(j) = 0.0d0

	end do

	! Contorno norte
	aux1  = Nf1/2
	aux2  = (Nf2-1)*Nf1
	raux1 = 0.5d0 / (aux1 - 1.0d0)
	do j = 2, aux1

		aux3  = aux2 + j
	
		raux2 = (j-1.5d0)*raux1

		xcf(aux3) = 0.5d0
		ycf(aux3) = 1.0d0 - raux2

	end do

	do j = aux1 + 1, Nf1-1

		aux3 = aux2 + j

		raux2 = ((j-aux1)-0.5d0)*raux1

		xcf(aux3) = 0.5d0 + raux2
		ycf(aux3) = 0.5d0

	end do

	! Contorno oeste
	raux1 = 0.5d0 / N2
	do i = 2, Nf2-1
		
		aux1  = (i-1)*Nf1 + 1
		raux2 =	(i-1.5d0)*raux1

		xcf(aux1) = raux2
		ycf(aux1) = 1.0d0

	end do

	! Contorno leste
	raux1 = 0.5d0 / N2
	do i = 2, Nf2-1

		aux1  = i*Nf1
		raux2 = (i-1.5d0)*raux1

		xcf(aux1) = 1.0d0
		ycf(aux1) = raux2

	end do

	! Canto sudoeste
	xcf(1) = 0.0d0
	ycf(1) = 1.0d0

	! Canto sudeste
	xcf(Nf1) = 1.0d0
	ycf(Nf1) = 0.0d0

	! Canto noroeste
	xcf( (Nf2-1)*Nf1 + 1) = 0.5d0
	ycf( (Nf2-1)*Nf1 + 1) = 1.0d0

	! Canto nordeste
	xcf( Nf2*Nf1 ) = 1.0d0
	ycf( Nf2*Nf1 ) = 0.5d0


end subroutine centroides

!--------------------------------------------------------------------

subroutine centroides_fict_nest

	! Canto sudoeste
	xcf(1) = 0.0d0
	ycf(1) = 1.0d0

	! Canto sudeste
	xcf(Nf1) = 1.0d0
	ycf(Nf1) = 0.0d0

	! Canto noroeste
	xcf( (Nf2-1)*Nf1 + 1) = 0.5d0
	ycf( (Nf2-1)*Nf1 + 1) = 1.0d0

	! Canto nordeste
	xcf( Nf2*Nf1 ) = 1.0d0
	ycf( Nf2*Nf1 ) = 0.5d0

	! Contorno sul
	aux1  = Nf1/2
	raux1 = 1.0d0 / (aux1 - 1.0d0)
	do j = 2, aux1
	
		raux2 = (j-1.5d0)*raux1
		raux3 = 1.0d0 - raux2
		raux4 = raux3 - ycf(j+Nf1)

		xcf(j) = -xcf(j+Nf1)
		ycf(j) = raux3 + raux4

	end do

	do j = aux1 + 1, Nf1-1

		raux2 = ((j-aux1)-0.5d0)*raux1
		raux3 = raux2 - xcf(j+Nf1)

		xcf(j) = raux2 + raux3
		ycf(j) = -ycf(j+Nf1)

	end do

	! Contorno norte
	aux1  = Nf1/2
	aux2  = (Nf2-1)*Nf1
	raux1 = 0.5d0 / (aux1 - 1.0d0)
	do j = 2, aux1

		aux3  = aux2 + j
	
		raux2 = (j-1.5d0)*raux1
		raux3 = 1.0d0 - raux2
		raux4 = raux3 - ycf(aux3-Nf1)

		xcf(aux3) = 1.0d0 - xcf(aux3-Nf1)
		ycf(aux3) = raux3 + raux4

	end do

	do j = aux1 + 1, Nf1-1

		aux3 = aux2 + j

		raux2 = ((j-aux1)-0.5d0)*raux1
		raux3 = 0.5d0 + raux2
		raux4 = raux3 - xcf(aux3-Nf1)

		xcf(aux3) = raux3 + raux4
		ycf(aux3) = 1.0d0 - ycf(aux3-Nf1)

	end do

	! Contorno oeste
	raux1 = 0.5d0 / N2
	do i = 2, Nf2-1
		
		aux1  = (i-1)*Nf1 + 1
		raux2 =	(i-1.5d0)*raux1
		raux3 = raux2 - xcf(aux1+1)

		xcf(aux1) = raux2 + raux3
		ycf(aux1) = 2.0d0 - ycf(aux1+1)

	end do

	! Contorno leste
	raux1 = 0.5d0 / N2
	do i = 2, Nf2-1

		aux1  = i*Nf1
		raux2 = (i-1.5d0)*raux1
		raux3 = raux2 - xcf(aux3-1)

		xcf(aux1) = 2.0d0 - xcf(aux1-1)
		ycf(aux1) = raux2 + raux3

	end do

end subroutine centroides_fict_nest

!--------------------------------------------------------------------

subroutine conectividades

conec = 0

! Volumes internos
do i = 2, N2-1
	
	aux1 = (i-1)*N1

	do j = 2, N1-1

		aux2 = aux1 + j

		conec(1,aux2) = aux2 - N1
		conec(2,aux2) = aux2 + 1
		conec(3,aux2) = aux2 + N1
		conec(4,aux2) = aux2 - 1

	end do

end do

! Volumes adjacentes ao contorno sul (exceto cantos)
do j = 2, N1-1
	
	conec(2,j) = j + 1
	conec(3,j) = j + N1
	conec(4,j) = j - 1

	if ( (conec(2,j) > Nt).or.(conec(2,j) < 1) ) conec(2,j) = 0
	if ( (conec(3,j) > Nt).or.(conec(3,j) < 1) ) conec(3,j) = 0
	if ( (conec(4,j) > Nt).or.(conec(4,j) < 1) ) conec(4,j) = 0

end do

! Volumes adjacentes ao contorno norte (exceto cantos)
aux1 = (N2-1)*N1
do j = 2, N1-1

	aux2 = aux1 + j

	conec(1,aux2) = aux2 - N1
	conec(2,aux2) = aux2 + 1
	conec(4,aux2) = aux2 - 1

	if ( (conec(1,aux2) > Nt).or.(conec(1,aux2) < 1) ) conec(1,aux2) = 0
	if ( (conec(2,aux2) > Nt).or.(conec(2,aux2) < 1) ) conec(2,aux2) = 0
	if ( (conec(4,aux2) > Nt).or.(conec(4,aux2) < 1) ) conec(4,aux2) = 0

end do

! Volumes adjacentes ao contorno oeste (exceto cantos)
do i = 2, N2-1

	aux2 = (i-1)*N1 + 1

	conec(1,aux2) = aux2 - N1
	conec(2,aux2) = aux2 + 1
	conec(3,aux2) = aux2 + N1

	if ( (conec(1,aux2) > Nt).or.(conec(1,aux2) < 1) ) conec(1,aux2) = 0
	if ( (conec(2,aux2) > Nt).or.(conec(2,aux2) < 1) ) conec(2,aux2) = 0
	if ( (conec(3,aux2) > Nt).or.(conec(3,aux2) < 1) ) conec(3,aux2) = 0

end do

! Volumes adjacentes ao contorno leste (exceto cantos)
do i = 2, N2-1

	aux2 = i*N1

	conec(1,aux2) = aux2 - N1
	conec(3,aux2) = aux2 + N1
	conec(4,aux2) = aux2 - 1

	if ( (conec(1,aux2) > Nt).or.(conec(1,aux2) < 1) ) conec(1,aux2) = 0
	if ( (conec(3,aux2) > Nt).or.(conec(3,aux2) < 1) ) conec(3,aux2) = 0
	if ( (conec(4,aux2) > Nt).or.(conec(4,aux2) < 1) ) conec(4,aux2) = 0

end do

! Canto SW
conec(2,1) = 2
conec(3,1) = 1 + N1

if ( (conec(2,1) > Nt).or.(conec(2,1) < 1) ) conec(2,1) = 0
if ( (conec(3,1) > Nt).or.(conec(3,1) < 1) ) conec(3,1) = 0

! Canto SE
conec(3,N1) = 2*N1
conec(4,N1) = N1 - 1

if ( (conec(3,N1) > Nt).or.(conec(3,N1) < 1) ) conec(3,N1) = 0
if ( (conec(4,N1) > Nt).or.(conec(4,N1) < 1) ) conec(4,N1) = 0

! Canto NW
aux2 = (N2-1)*N1 + 1
conec(1,aux2) = aux2 - N1
conec(2,aux2) = aux2 + 1

if ( (conec(1,aux2) > Nt).or.(conec(1,aux2) < 1) ) conec(1,aux2) = 0
if ( (conec(2,aux2) > Nt).or.(conec(2,aux2) < 1) ) conec(2,aux2) = 0

! Canto NE
aux2 = N1*N2
conec(1,aux2) = aux2 - N1
conec(4,aux2) = aux2 - 1

if ( (conec(1,aux2) > Nt).or.(conec(1,aux2) < 1) ) conec(1,aux2) = 0
if ( (conec(4,aux2) > Nt).or.(conec(4,aux2) < 1) ) conec(4,aux2) = 0

end subroutine conectividades

!--------------------------------------------------------------------

subroutine conectividades_fict

	conec = 0

	! Volumes internos
	do i = 2, Nf2-1

		aux1 = (i-1)*Nf1

		do j = 2, Nf1-1

			aux2 = aux1 + j

			conec(1,aux2) = 4
			conec(2,aux2) = aux2 - Nf1
			conec(3,aux2) = aux2 + 1
			conec(4,aux2) = aux2 + Nf1
			conec(5,aux2) = aux2 - 1

		end do

	end do

	! Contorno sul
	do j = 2, Nf1-1
		
		conec(1,j) = 1
		conec(2,j) = j + Nf1

	end do

	! Contorno norte
	aux1  = (Nf2-1)*Nf1
	do j = 2, Nf1-1

		aux2 = aux1 + j

		conec(1,aux2) = 1
		conec(2,aux2) = aux2 - Nf1

	end do

	! Contorno oeste
	do i = 2, Nf2-1

		aux1 = (i-1)*Nf1 + 1

		conec(1,aux1) = 1
		conec(2,aux1) = aux1 + 1

	end do

	! Contorno leste
	do i = 2, Nf2-1

		aux1 = i*Nf1

		conec(1,aux1) = 1
		conec(2,aux1) = aux1 - 1

	end do

end subroutine conectividades_fict

!--------------------------------------------------------------------

end module metricas