module vector implicit none contains !========================================================== !calculate cross product !--------------------------- function cross_product(a,b) result (c) real(8), intent(in) :: a(:),b(:) real(8), allocatable :: c(:) if(size(a) /= size(b)) then stop "wrong number on size a, b" endif allocate(c(size(a,1))) if(size(c,1)==3) then c(1) = a(2)*b(3) - a(3)*b(2) c(2) = a(3)*b(1) - a(1)*b(3) c(3) = a(1)*b(2) - a(2)*b(1) else stop "wrong number at cross_product" endif end function cross_product !========================================================= !calculate diadic !---------------------- function diadic(a,b) result(c) real(8), intent(in) :: a(:), b(:) real(8), allocatable :: c(:,:) integer n,i,j allocate(c(size(a),size(b) ) ) do i=1,size(a) do j=1,size(b) c(i,j)=a(i)*b(j) enddo enddo end function diadic !========================================================== !calculate gz !-------------- subroutine calcgz(x2,x11,x12,nod_coord,gzi) real(8), intent(in) :: nod_coord(:,:) real(8),intent(out) :: gzi real(8) l real(8),allocatable::avec(:) integer,intent(in):: x2,x11,x12 allocate(avec(2)) l = dot_product( nod_coord(x12,1:2) - nod_coord(x11,1:2), & nod_coord(x12,1:2) - nod_coord(x11,1:2) ) l=l**(1.0d0/2.0d0) avec(1:2) = ( nod_coord(x12,1:2) - nod_coord(x11,1:2) )/l if(l==0.0d0)then print *, "calcgz l=0" gzi=0.0d0 else gzi=1.0d0/l*dot_product( nod_coord(x2,1:2) -nod_coord(x11,1:2),avec(1:2) ) endif deallocate(avec) end subroutine calcgz !========================================================== subroutine eigen_2d(Amat,eigenvector) real(8),intent(in)::Amat(:,:) real(8),intent(inout)::eigenvector(:,:) real(8)::b,c,phy,eigenvalue(2) integer i,j eigenvalue(:)=0.0d0 eigenvector(:,:)=0.0d0 b=-1.0d0*(Amat(1,1)+Amat(2,2)) c=Amat(1,1)*Amat(2,2)-Amat(1,2)*Amat(1,2) if(Amat(1,2)/=Amat(2,1) )then stop "input matrice is not symmetric" endif do i=1,2 eigenvalue(i)=(-1.0d0*b+((-1.0d0)**dble(i))*(b*b-4.0d0*c)**(1.0d0/2.0d0))*(0.50d0) enddo do i=1,2 if(Amat(1,2)==0 )then cycle elseif(Amat(1,2)/=0 )then phy=atan( (eigenvalue(i)-Amat(1,1))/Amat(1,2) ) do j=1,2 eigenvector(i,1:2)=(/cos(phy),sin(phy)/) enddo endif enddo do i=1,2 eigenvector(i,:)=eigenvalue(i)*eigenvector(i,:) enddo end subroutine eigen_2d !========================================================== function signmm(a) result(b) real(8),intent(in)::a real(8) b if(a>0)then b=1.0d0 elseif(a<0)then b=-1.0d0 elseif(a==0)then b=0.0d0 else stop "ERROR: Invalid Real(8) in function_signm" endif end function signmm !========================================================== end module vector