preprocessing_2Delements

public subroutine preprocessing_2Delements(me, mshFile)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Arguments

Type	Intent	Optional	Attributes		Name
type(mesh2D_prop_t),	intent(out)			::	me
character(len=*),	intent(in)			::	mshFile
Called by

Help
Source Code

subroutine preprocessing_2Delements(me,mshFile)  
    type(mesh2D_prop_t), intent(out) :: me
    character(*), intent(in) :: mshFile

    type(elements_1D_t), allocatable :: elem1D(:) ! Type for 1D element (only nodes are required)
    
    real(dp) :: den, x_face, y_face, check
    real(dp) :: nn(2),cf(2)
 
    integer, allocatable :: cnt(:,:), nod1(:,:)
    integer :: nod_shar(2)
    integer :: cnt_nod_shar,cnt_nod_shar_1D,nod_alone,ii
    integer :: Nb_1D_elements
    integer :: line,nb_n,nb_e,nb_p,nb_1De
    integer :: j_e,nb_v,j_v,j_1De
    integer :: pt1,pt2,pt3,pt4,pt5

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    open(unit=1,file=mshFile) ! msh file generated by GMSH     -->    export - File.msh - Version 2 ASCII  (without any additional option)
    do line=1,3 ! Mesh Format
      read(1,*)
    enddo
    read(1,*)  ! Physical names
    read(1,*) me%Nb_label
    do line=1,me%Nb_label
      read(1,*)
    enddo
    read(1,*)
    read(1,*)  ! nodes
    read(1,*) me%nb_nodes
    do line=1,me%nb_nodes
      read(1,*)
    enddo
    read(1,*)
    read(1,*)
    read(1,*) me%nb_elements  ! Elements (multi-D at this stage)
    do line=1,me%nb_elements
      read(1,*) pt1, pt2
      if(pt2==2) then
        Nb_1D_elements=pt1-1
        exit
      endif
    enddo
    me%nb_elements=me%nb_elements-Nb_1D_elements ! to keep only the number of 2D elements
    close(1)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    allocate(me%LabelName(me%Nb_label))
    allocate(me%node(me%nb_nodes),me%elem(me%nb_elements),elem1D(Nb_1D_elements))
    allocate(cnt(me%nb_elements,me%nb_elements),nod1(me%nb_elements,me%nb_elements))
  
    do nb_n=1,me%nb_nodes
      me%node(nb_n)%x_coord=0.0_dp;me%node(nb_n)%y_coord=0.0_dp
    enddo
    do nb_e=1,me%nb_elements
      me%elem(nb_e)%nod(:)=0
    enddo
    do nb_e=1,Nb_1D_elements
      elem1D(nb_e)%nod(:)=0
    enddo

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    open(unit=1,file=mshFile)
    do line=1,3 ! Mesh Format
      read(1,*)
    enddo
    read(1,*)  ! Physical names
    read(1,*) 
    do line=1,me%Nb_label
      read(1,*) pt1, pt2, me%LabelName(line)  ! In practice, pt2=line
    enddo
    read(1,*)
    read(1,*)  ! nodes
    read(1,*)
    do line=1,me%nb_nodes
      read(1,*) pt1, me%node(line)%x_coord, me%node(line)%y_coord, check
    enddo
    read(1,*)
    read(1,*)  ! Elements
    read(1,*)
    do line=1,Nb_1D_elements
      read(1,*) pt1, pt2, pt3, elem1D(line)%PhysIdx, pt4, elem1D(line)%nod(:)
    enddo
    do line=1,me%nb_elements
      read(1,*) pt1, pt2, pt3, me%elem(line)%PhysIdxE, pt5, me%elem(line)%nod(:)
    enddo
    close(1) ! End of reading
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    cnt(:,:)=0
    me%nnz=0
    do nb_e=1,me%nb_elements

      me%elem(nb_e)%PhysE=trim(me%LabelName(me%elem(nb_e)%PhysIdxE))
  
      me%elem(nb_e)%surface=0.5_dp*&
                  abs(me%node(me%elem(nb_e)%nod(1))%x_coord*(me%node(me%elem(nb_e)%nod(2))%y_coord-&
                      me%node(me%elem(nb_e)%nod(3))%y_coord)+&
                      me%node(me%elem(nb_e)%nod(2))%x_coord*(me%node(me%elem(nb_e)%nod(3))%y_coord-&
                      me%node(me%elem(nb_e)%nod(1))%y_coord)+&
                      me%node(me%elem(nb_e)%nod(3))%x_coord*(me%node(me%elem(nb_e)%nod(1))%y_coord-&
                      me%node(me%elem(nb_e)%nod(2))%y_coord))
  
      me%elem(nb_e)%center%x_coord=(me%node(me%elem(nb_e)%nod(1))%x_coord+me%node(me%&
                                           elem(nb_e)%nod(2))%x_coord+&
                                           me%node(me%elem(nb_e)%nod(3))%x_coord)/3.0_dp
      me%elem(nb_e)%center%y_coord=(me%node(me%elem(nb_e)%nod(1))%y_coord+me%node(me%&
                                           elem(nb_e)%nod(2))%y_coord+&
                                           me%node(me%elem(nb_e)%nod(3))%y_coord)/3.0_dp
  
      me%elem(nb_e)%NbNeigh=0
      do nb_v=1,me%nb_elements
        do j_e=1,3
          do j_v=1,3    
            if((me%elem(nb_e)%nod(j_e)==me%elem(nb_v)%nod(j_v)) .and. (nb_v/=nb_e)) then
              cnt(nb_e,nb_v)=cnt(nb_e,nb_v)+1
              if(cnt(nb_e,nb_v)/=2) nod1(nb_e,nb_v)=me%elem(nb_e)%nod(j_e)
            endif
          enddo
  
          if(cnt(nb_e,nb_v)==2) then
            me%elem(nb_e)%NbNeigh=me%elem(nb_e)%NbNeigh+1            
            me%elem(nb_e)%Neigh(me%elem(nb_e)%NbNeigh)=nb_v
  
            me%elem(nb_e)%NeighCenter(me%elem(nb_e)%NbNeigh)%x_coord=(me%node(me%elem(nb_v)%&
                                                                      nod(1))%x_coord+&
                                                                      me%node(me%elem(nb_v)%nod(2))%x_coord+&
                                                                      me%node(me%elem(nb_v)%nod(3))%x_coord)/3.0_dp
            me%elem(nb_e)%NeighCenter(me%elem(nb_e)%NbNeigh)%y_coord=(me%node(me%elem(nb_v)%&
                                                                      nod(1))%y_coord+&
                                                                      me%node(me%elem(nb_v)%nod(2))%y_coord+&
                                                                      me%node(me%elem(nb_v)%nod(3))%y_coord)/3.0_dp
  
            me%elem(nb_e)%face(me%elem(nb_e)%NbNeigh)=sqrt((me%node(me%elem(nb_e)%nod(j_e))%&
                                               x_coord-me%node(nod1(nb_e,nb_v))%x_coord)**2+&
                                               (me%node(me%elem(nb_e)%nod(j_e))%y_coord-&
                                               me%node(nod1(nb_e,nb_v))%y_coord)**2)
  

            x_face=0.5_dp*(me%node(me%elem(nb_e)%nod(j_e))%x_coord+me%node(nod1(nb_e,nb_v))%x_coord)
            y_face=0.5_dp*(me%node(me%elem(nb_e)%nod(j_e))%y_coord+me%node(nod1(nb_e,nb_v))%y_coord)

            me%elem(nb_e)%centTOface(me%elem(nb_e)%NbNeigh)%x_coord=x_face-me%elem(nb_e)%center%x_coord
            me%elem(nb_e)%centTOface(me%elem(nb_e)%NbNeigh)%y_coord=y_face-me%elem(nb_e)%center%y_coord

            me%elem(nb_e)%centTOfaceN(me%elem(nb_e)%NbNeigh)%x_coord=x_face-&
                                                         me%elem(nb_e)%NeighCenter(me%elem(nb_e)%NbNeigh)%x_coord
            me%elem(nb_e)%centTOfaceN(me%elem(nb_e)%NbNeigh)%y_coord=y_face-&
                                                         me%elem(nb_e)%NeighCenter(me%elem(nb_e)%NbNeigh)%y_coord
  

            den=me%elem(nb_e)%face(me%elem(nb_e)%NbNeigh)
            me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%x_coord=-(me%node(me%&
                                                                elem(nb_e)%nod(j_e))%y_coord-&
                                                                me%node(nod1(nb_e,nb_v))%y_coord)/den
            me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%y_coord= (me%node(me%&
                                                                elem(nb_e)%nod(j_e))%x_coord-&
                                                                me%node(nod1(nb_e,nb_v))%x_coord)/den
  
            ! check for unit normales to be outgoing
            nn(1)=me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%x_coord
            nn(2)=me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%y_coord
            cf(1)=me%elem(nb_e)%centTOface(me%elem(nb_e)%NbNeigh)%x_coord
            cf(2)=me%elem(nb_e)%centTOface(me%elem(nb_e)%NbNeigh)%y_coord
            if(dot_product(nn,cf)<=0.0_dp) then
              me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%x_coord=&
                                               -me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%x_coord
              me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%y_coord=&
                                               -me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%y_coord
            endif
            nn(1)=me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%x_coord
            nn(2)=me%elem(nb_e)%norm(me%elem(nb_e)%NbNeigh)%y_coord
            me%elem(nb_e)%Delta(me%elem(nb_e)%NbNeigh)=abs(dot_product(nn,cf))

            cf(1)=me%elem(nb_e)%centTOfaceN(me%elem(nb_e)%NbNeigh)%x_coord
            cf(2)=me%elem(nb_e)%centTOfaceN(me%elem(nb_e)%NbNeigh)%y_coord
            me%elem(nb_e)%DeltaNeigh(me%elem(nb_e)%NbNeigh)=abs(dot_product(-nn,cf))
  
            cnt(nb_e,nb_v)=-10 ! to avoid duplication
  
          endif
        enddo
      enddo
  
      if(me%elem(nb_e)%NbNeigh==1) then

        me%elem(nb_e)%Neigh(2)=0
        me%elem(nb_e)%Neigh(3)=0

        nod_shar(:)=0
        do j_e=1,3
          cnt_nod_shar=0
          ! Neigh 1 (only one neighbour)
          do j_v=1,3
            if(me%elem(nb_e)%nod(j_e)==me%elem(me%elem(nb_e)%Neigh(1))%nod(j_v)) then
              cnt_nod_shar=cnt_nod_shar+1
            endif
          enddo
  
          if(cnt_nod_shar==1) then
            if(nod_shar(1)==0) then
              nod_shar(1)=me%elem(nb_e)%nod(j_e)
            else
              nod_shar(2)=me%elem(nb_e)%nod(j_e)
            endif
          else if(cnt_nod_shar==0) then
              nod_alone=me%elem(nb_e)%nod(j_e)
          endif
        enddo

        do ii=2,3
          me%elem(nb_e)%face(ii)=sqrt((me%node(nod_shar(ii-1))%x_coord-me%node(nod_alone)%x_coord)**2+&
                                    (me%node(nod_shar(ii-1))%y_coord-&
                                    me%node(nod_alone)%y_coord)**2)
    
          x_face=0.5_dp*(me%node(nod_shar(ii-1))%x_coord+me%node(nod_alone)%x_coord)
          y_face=0.5_dp*(me%node(nod_shar(ii-1))%y_coord+me%node(nod_alone)%y_coord)
          me%elem(nb_e)%centTOface(ii)%x_coord=x_face-me%elem(nb_e)%center%x_coord
          me%elem(nb_e)%centTOface(ii)%y_coord=y_face-me%elem(nb_e)%center%y_coord
    
          den=me%elem(nb_e)%face(ii)
          me%elem(nb_e)%norm(ii)%x_coord=-(me%node(nod_shar(ii-1))%y_coord-me%node(nod_alone)%y_coord)/den
          me%elem(nb_e)%norm(ii)%y_coord= (me%node(nod_shar(ii-1))%x_coord-me%node(nod_alone)%x_coord)/den
    
          ! check for unit normales to be outgoing
          nn(1)=me%elem(nb_e)%norm(ii)%x_coord
          nn(2)=me%elem(nb_e)%norm(ii)%y_coord
          cf(1)=me%elem(nb_e)%centTOface(ii)%x_coord
          cf(2)=me%elem(nb_e)%centTOface(ii)%y_coord
          if(dot_product(nn,cf)<=0.0_dp) then
            me%elem(nb_e)%norm(ii)%x_coord=-me%elem(nb_e)%norm(ii)%x_coord
            me%elem(nb_e)%norm(ii)%y_coord=-me%elem(nb_e)%norm(ii)%y_coord
          endif
          nn(1)=me%elem(nb_e)%norm(ii)%x_coord
          nn(2)=me%elem(nb_e)%norm(ii)%y_coord
          me%elem(nb_e)%Delta(ii)=abs(dot_product(nn,cf))
        enddo

        me%elem(nb_e)%PhysN(1)=trim(me%LabelName(me%elem(me%elem(nb_e)%Neigh(1))%PhysIdxE))

        do ii=2,3
          do nb_1De=1,Nb_1D_elements
            cnt_nod_shar_1D=0
            do j_1De=1,2
              if(nod_shar(ii-1)==elem1D(nb_1De)%nod(j_1De) .or. nod_alone==elem1D(nb_1De)%nod(j_1De)) then
                cnt_nod_shar_1D=cnt_nod_shar_1D+1
              endif
            enddo        
            if(cnt_nod_shar_1D==2) then
              me%elem(nb_e)%PhysN(ii)=trim(me%LabelName(elem1D(nb_1De)%PhysIdx))
            endif
          enddo
        enddo
  
      else if(me%elem(nb_e)%NbNeigh==3) then
        do j_v=1,3
          me%elem(nb_e)%PhysN(j_v)=trim(me%LabelName(me%elem(me%elem(nb_e)%Neigh(j_v))%PhysIdxE))
        enddo

      else if(me%elem(nb_e)%NbNeigh==2) then
        me%elem(nb_e)%Neigh(3)=0
  
        nod_shar(:)=0
        do j_e=1,3
          cnt_nod_shar=0
          ! Neigh 1
          do j_v=1,3
            if(me%elem(nb_e)%nod(j_e)==me%elem(me%elem(nb_e)%Neigh(1))%nod(j_v)) then
              cnt_nod_shar=cnt_nod_shar+1
            endif
          enddo
  
          ! Neigh 2
          do j_v=1,3
            if(me%elem(nb_e)%nod(j_e)==me%elem(me%elem(nb_e)%Neigh(2))%nod(j_v)) then
              cnt_nod_shar=cnt_nod_shar+1
            endif
          enddo
  
          if(cnt_nod_shar==1) then
            if(nod_shar(1)==0) then
              nod_shar(1)=me%elem(nb_e)%nod(j_e)
            else
              nod_shar(2)=me%elem(nb_e)%nod(j_e)
            endif
          endif
        enddo
  
        me%elem(nb_e)%face(3)=sqrt((me%node(nod_shar(2))%x_coord-me%node(nod_shar(1))%x_coord)**2+&
                                  (me%node(nod_shar(2))%y_coord-&
                                  me%node(nod_shar(1))%y_coord)**2)
  
        x_face=0.5_dp*(me%node(nod_shar(2))%x_coord+me%node(nod_shar(1))%x_coord)
        y_face=0.5_dp*(me%node(nod_shar(2))%y_coord+me%node(nod_shar(1))%y_coord)
        me%elem(nb_e)%centTOface(3)%x_coord=x_face-me%elem(nb_e)%center%x_coord
        me%elem(nb_e)%centTOface(3)%y_coord=y_face-me%elem(nb_e)%center%y_coord
  
        den=me%elem(nb_e)%face(3)
        me%elem(nb_e)%norm(3)%x_coord=-(me%node(nod_shar(2))%y_coord-me%node(nod_shar(1))%y_coord)/den
        me%elem(nb_e)%norm(3)%y_coord= (me%node(nod_shar(2))%x_coord-me%node(nod_shar(1))%x_coord)/den
  
        ! check for unit normales to be outgoing
        nn(1)=me%elem(nb_e)%norm(3)%x_coord
        nn(2)=me%elem(nb_e)%norm(3)%y_coord
        cf(1)=me%elem(nb_e)%centTOface(3)%x_coord
        cf(2)=me%elem(nb_e)%centTOface(3)%y_coord
        if(dot_product(nn,cf)<=0.0_dp) then
          me%elem(nb_e)%norm(3)%x_coord=-me%elem(nb_e)%norm(3)%x_coord
          me%elem(nb_e)%norm(3)%y_coord=-me%elem(nb_e)%norm(3)%y_coord
        endif
        nn(1)=me%elem(nb_e)%norm(3)%x_coord
        nn(2)=me%elem(nb_e)%norm(3)%y_coord
        me%elem(nb_e)%Delta(3)=abs(dot_product(nn,cf))
  
        do j_v=1,2
          me%elem(nb_e)%PhysN(j_v)=trim(me%LabelName(me%elem(me%elem(nb_e)%Neigh(j_v))%PhysIdxE))
        enddo

        do nb_1De=1,Nb_1D_elements
          cnt_nod_shar_1D=0
          do j_1De=1,2
            if(nod_shar(1)==elem1D(nb_1De)%nod(j_1De) .or. nod_shar(2)==elem1D(nb_1De)%nod(j_1De)) then
              cnt_nod_shar_1D=cnt_nod_shar_1D+1
            endif
          enddo        
          if(cnt_nod_shar_1D==2) then
            me%elem(nb_e)%PhysN(3)=trim(me%LabelName(elem1D(nb_1De)%PhysIdx))
          endif
        enddo
  
      endif

      me%nnz=me%nnz+1+me%elem(nb_e)%NbNeigh

      if(me%elem(nb_e)%NbNeigh==1) then
        print*,"ACCUMULATION TO BE CONSIDERED FOR 2D CELLS THAT ARE CONNECTED TO 2 DIFFERENT LABELS ?"
        print*," check if it is the same for the 2D cells connected to labels"
        read(*,*)
      endif
 
    enddo

    print*,mshFile,' --->  Preprocessing performed'

end subroutine preprocessing_2Delements
preprocessing_2Delements Subroutine

Contents

Source Code

public subroutine preprocessing_2Delements(me, mshFile)

Arguments

Called by

Source Code
