! ! (c) 2019 Guide Star Engineering, LLC ! This Software was developed for the US Nuclear Regulatory Commission (US NRC) ! under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under ! Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007 ! ! ! NEMO - Numerical Engine (for) Multiphysics Operators ! Copyright (c) 2007, Stefano Toninel ! Gian Marco Bianchi University of Bologna ! David P. Schmidt University of Massachusetts - Amherst ! Salvatore Filippone University of Rome Tor Vergata ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without modification, ! are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! 3. Neither the name of the NEMO project nor the names of its contributors ! may be used to endorse or promote products derived from this software ! without specific prior written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED ! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE ! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! !--------------------------------------------------------------------------------- SUBMODULE (tools_mesh) cmp_mesh_implementation USE class_face IMPLICIT NONE CONTAINS MODULE PROCEDURE cmp_moving_surf USE class_psblas!, ONLY : mypnum_ USE class_connectivity USE class_surface!, ONLY : alloc_surface USE class_vertex IMPLICIT NONE !! $Id: cmp_moving_surf.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: Allocates and calculates shapes of moving surfaces !! INTEGER :: ib, info IF(mypnum_() == 0) THEN WRITE(*,*) 'Autodetecting virtual surfaces:' ! Not necessary since INTENT(OUT) !!$ if (associated (surf) ) then !!$ call abort_psblas !!$ endif ! Allocates pointer of SURFACE objects ALLOCATE(surf(nbc),stat=info) IF (info /= 0) THEN WRITE(*,100) CALL abort_psblas ENDIF ! Identifies surfaces DO ib = 1, nbc CALL alloc_surface(v2b,ib,verts,surf(ib)) END DO WRITE(*,*) ENDIF 100 FORMAT(' ERROR! Failure to allocate memory in CMP_MOVING_SURF') END PROCEDURE cmp_moving_surf MODULE PROCEDURE cmp_mesh_c2c USE class_psblas USE class_connectivity IMPLICIT NONE !! $Id: cmp_mesh_c2c.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: !! INTEGER :: ic, IF, im, is, info INTEGER :: n, ncells, nfl_faces INTEGER, ALLOCATABLE :: itab(:,:), kconn(:) IF(mypnum_() == 0) THEN WRITE(*,*) 'Computing cells adjacency graph C2C' n = f2c%max_conn() ! Maximum connectivity degree ncells = f2c%nel_() ! Number of cells nfl_faces = COUNT(faces%flag_() == 0) ! Number of fluid faces ALLOCATE(itab(n,ncells),kconn(ncells),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF itab = 0 kconn = 0 ! WARNING! CMP_MESH_C2C is called before the global-to-local reallocation. ! Thus faces are still numbered in such a way that fluid ones come first DO IF = 1, nfl_faces im = faces(IF)%master_() is = faces(IF)%slave_() kconn(im) = kconn(im) + 1 kconn(is) = kconn(is) + 1 itab(kconn(im),im) = is itab(kconn(is),is) = im END DO ! Counts total amount of connectivities n = 0 DO ic = 1, ncells n = n + kconn(ic) END DO ! Creates C2C connectivity object CALL alloc_conn(c2c,nel=ncells,nconn=n) ! Sets C2C connectivity DO ic = 1, ncells n = kconn(ic) CALL c2c%set_ith_conn(ic,itab(1:n,ic)) END DO DEALLOCATE(itab,kconn,stat=info) IF(info /= 0) THEN WRITE(*,200) CALL abort_psblas END IF ENDIF 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_C2C') 200 FORMAT(' ERROR! Memory deallocation failure in CMP_MESH_C2C') END PROCEDURE cmp_mesh_c2c MODULE PROCEDURE cmp_mesh_f2b USE class_psblas USE class_connectivity IMPLICIT NONE !! $Id: cmp_mesh_f2b.f90 3175 2008-06-13 12:59:07Z sfilippo $ !! !! Description: !! INTEGER :: ib, IF, info, maxconn, n, nfaces,j INTEGER :: flag(SIZE(faces)), iface(SIZE(faces)), kf(-1:nbc) INTEGER, ALLOCATABLE :: iconn(:) IF(mypnum_() == 0) & & WRITE(*,*) 'Computing face to boundary connectivity F2B' nfaces = SIZE(faces) iface = (/(IF, IF = 1, nfaces)/) flag = faces%flag_() ! Computes maximum connectivity degree IF (.FALSE.) THEN maxconn = 0 DO ib = -1, nbc kf(ib) = COUNT(flag == ib) maxconn = MAX(maxconn,kf(ib)) END DO ELSE CALL psb_msort(flag(1:nfaces),ix=iface(1:nfaces),flag=psb_sort_keep_idx_) IF (nfaces > 0) THEN IF ((flag(1) < -1).OR.(flag(nfaces)>nbc)) THEN WRITE(0,*) 'Error in cmp_mesh_f2b: flags out of bounds ',flag(1),flag(nfaces) CALL abort_psblas END IF END IF maxconn = 0 kf(:) = 0 DO IF = 1, nfaces ib = flag(IF) kf(ib) = kf(ib) + 1 maxconn = MAX(maxconn,kf(ib)) END DO END IF CALL alloc_conn(f2b,nel=nbc+2,nconn=nfaces,lb=-1) ALLOCATE(iconn(maxconn),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF j = 0 DO ib = -1, nbc n = kf(ib) iconn(1:n) = iface(j+1:j+n) !!$ iconn(1:n) = pack(iface,flag == ib) CALL f2b%set_ith_conn(ib,iconn(1:n)) j = j + n END DO DEALLOCATE(iconn) 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_F2B') END PROCEDURE cmp_mesh_f2b MODULE PROCEDURE cmp_mesh_f2f USE class_psblas USE class_connectivity IMPLICIT NONE !! $Id: cmp_mesh_f2f.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: !! Provides the Face To Face adjacency graph. !! The criterium for associating the graph nodes is arbitrarly defined by the !! user. In this case, using a co-located cell-centered strategy for variable !! displacement, the F2F connectivity is obtained by linking each face to all !! the ones belonging to its master and slave cells. !! INTEGER :: i, ic, IF, if1, if2, info, j, n, ncells, ncf INTEGER, POINTER :: if2c(:) => NULL() INTEGER, ALLOCATABLE :: kconn(:), itab(:,:) IF(mypnum_() == 0) THEN WRITE(*,*) 'Computing faces adjacency graph F2F' n = (f2c%max_conn() - 1) * 2 ! F2F maximum connectivity degree ncells = f2c%nel_() ! Number of cells ALLOCATE(itab(n,nfaces),kconn(nfaces),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF kconn = 0 DO ic = 1, ncells CALL f2c%get_ith_conn(if2c,ic) ncf = SIZE(if2c) DO i = 1, ncf - 1 if1 = if2c(i) DO j = i + 1, ncf if2 = if2c(j) ! Add direct arc kconn(if1) = kconn(if1) + 1 itab(kconn(if1),if1) = if2 ! Add indirect arc kconn(if2) = kconn(if2) + 1 itab(kconn(if2),if2) = if1 END DO END DO END DO ! Computes total number of connectivities n = 0 DO IF = 1, nfaces n = n + kconn(IF) END DO CALL alloc_conn(f2f,nel=nfaces,nconn=n) ! Sets F2F connectivity DO IF = 1, nfaces n = kconn(IF) CALL f2f%set_ith_conn(IF,itab(1:n,IF)) END DO ! Frees memory NULLIFY(if2c) DEALLOCATE(itab,kconn) ENDIF 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_F2F') END PROCEDURE cmp_mesh_f2f MODULE PROCEDURE cmp_mesh_part USE class_psblas USE class_connectivity USE tools_part USE part_graph, ONLY : bld_part_graph, get_part_graph USE part_block, ONLY : bld_part_block USE part_random, ONLY : bld_part_random IMPLICIT NONE !! $Id: cmp_mesh_part.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: !! INTEGER :: info, nprocs, ncells INTEGER, ALLOCATABLE :: xadj_glob(:) INTEGER, ALLOCATABLE :: adjncy_glob(:) IF (mypnum_() == 0) THEN CALL sw_par%tic() nprocs = nprocs_() IF(nprocs > 1 ) WRITE(*,*) ncells = c2c%nel_() ! Checks status of PART_CELLS vector IF(ALLOCATED(part_cells)) DEALLOCATE(part_cells) ALLOCATE(part_cells(ncells),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF ! Retrieves connectivity data in CSR-like format IF(ipart > 0) THEN ! ParMetis requires global connectivity in CSR format, therefore ! the call to GET_CONN_CSR is done here, outside BLD_PART_GRAPH. CALL c2c%get_conn_csr(xadj_glob,adjncy_glob) END IF SELECT CASE(ipart) CASE(-1) ! Random partitioning CALL bld_part_random(ncells,nprocs,part_cells) CASE(0) ! HPF Block partitioning CALL bld_part_block(ncells,nprocs,part_cells) CASE(1) ! METIS Partitioning CALL bld_part_graph(xadj_glob,adjncy_glob,ipart=1) CALL get_part_graph(part_cells) CASE DEFAULT WRITE(*,200) CALL abort_psblas END SELECT WRITE(*,*) IF(ipart > 0) THEN DEALLOCATE(xadj_glob,adjncy_glob,stat=info) IF(info /= 0) THEN WRITE(*,300) CALL abort_psblas END IF END IF CALL sw_par%toc() ENDIF 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_PART') 200 FORMAT(' ERROR! Unknown partitioning method in CMP_MESH_PART') 300 FORMAT(' ERROR! Memory deallocation failure in CMP_MESH_PART') END PROCEDURE cmp_mesh_part MODULE PROCEDURE cmp_mesh_renum USE class_psblas USE class_cell USE class_connectivity USE renum IMPLICIT NONE !! $Id: cmp_mesh_renum.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: !! INTEGER :: icontxt, mypnum icontxt = icontxt_() mypnum = mypnum_() CALL sw_ord%tic() SELECT CASE(irenum) CASE(0) IF(mypnum == 0) THEN WRITE(*,*) 'Matrix renumbering: none' WRITE(*,*) END IF CASE(1) ! Builds permutation array on P0, according to IRENUM value. IF(mypnum == 0) THEN CALL start_renum(irenum,c2c) !! call print_renum(0) ! Applies permutation to cell-related structures CALL apply_renum(c2c,to_a_and_b_) CALL apply_renum(v2c,to_b_) CALL apply_renum(f2c,to_b_) CALL apply_renum(c2g,to_a_) CALL apply_renum(cells) CALL apply_renum(faces) ! Deallocates private storage of renum module CALL stop_renum ENDIF CASE DEFAULT WRITE(*,100) CALL abort_psblas END SELECT CALL sw_ord%toc() 100 FORMAT(' ERROR! Unknown renumbering method in CMP_MESH_RENUM') END PROCEDURE cmp_mesh_renum MODULE PROCEDURE cmp_mesh_v2b USE class_psblas USE class_connectivity IMPLICIT NONE !! $Id: cmp_mesh_v2b.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: Sets up vertex to boundary connectivity, so that given a boundary, the !! GET_ITH_CONN method can be used to obtain the list of vertices that belong to the boundary. !! Boundary ID's correspond to those used to flag faces: !! !!---------------------------| !! Vertex Flag | Vertex Type | !!---------------------------| !! 0 | interior | !! < 0 | proc. bndry | !! > 0 | phys. bndry | !!---------------------------| !! LOGICAL :: inspected(0:nbc) INTEGER :: flag, IF, info, iv, k, n, nverts INTEGER :: kconn(0:nbc) INTEGER, ALLOCATABLE :: itab(:,:) INTEGER, POINTER :: if2v(:) => NULL() TYPE(connectivity) :: f2v IF(mypnum_() == 0) THEN WRITE(*,*) 'Computing vertex to boundary connectivity V2B' CALL v2f%get_dual_conn(f2v) nverts = f2v%nel_() ALLOCATE(itab(nverts,0:nbc),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF itab = 0 kconn = 0 DO iv = 1, nverts inspected(:) = .FALSE. CALL f2v%get_ith_conn(if2v,iv) n = SIZE(if2v) DO k = 1, n IF = if2v(k) flag = faces(IF)%flag_() IF(flag > 0 .AND. .NOT.inspected(flag)) THEN inspected(flag) = .TRUE. kconn(flag) = kconn(flag) + 1 itab(kconn(flag),flag) = iv END IF END DO IF(ALL(inspected .EQV. .FALSE.)) THEN kconn(0) = kconn(0) + 1 itab(kconn(0),0) = iv END IF END DO ! Computes number of V2B connectivity n = 0 DO k = 0, nbc n = n + kconn(k) END DO CALL alloc_conn(v2b,nel=nbc+1,nconn=n,lb=0) DO k = 0, nbc n = kconn(k) CALL v2b%set_ith_conn(k,itab(1:n,k)) END DO DEALLOCATE(itab) NULLIFY(if2v) CALL free_conn(f2v) ENDIF 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2B') END PROCEDURE cmp_mesh_v2b MODULE PROCEDURE cmp_mesh_v2e USE class_psblas USE class_connectivity IMPLICIT NONE !! $Id: cmp_mesh_v2e.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: !! Provides the Vertex To Edge connectivity. !! LOGICAL, ALLOCATABLE :: inspected(:), twice(:) INTEGER :: ie, ie1, ie2 INTEGER :: iv, iv1, iv2, iv3, iv4 INTEGER :: i, IF, info, j, n, nedges, nfaces, nverts INTEGER, POINTER :: iv2f(:) => NULL() INTEGER, POINTER :: iv2e(:) => NULL(), ie2v(:) => NULL() TYPE(connectivity) :: work, e2v TYPE(stopwatch) :: sw_v2e INTEGER :: d1, d2, s1, s2 ! V2E makes sense only for 3D meshes! IF(ncd == 2) RETURN IF(mypnum_() == 0) & WRITE(*,*) 'Computing vertex to edge connectivity V2E' ! Constructs stopwatch sw_v2e = stopwatch_(icontxt_()) ! Start timing CALL sw_v2e%tic() IF(mypnum_() == 0) THEN nfaces = v2f%nel_() nedges = 0 DO IF = 1, nfaces CALL v2f%get_ith_conn(iv2f,IF) nedges = nedges + SIZE(iv2f) END DO CALL alloc_conn(work,nel=nedges,nconn=2 * nedges) ie = 1 DO IF = 1, nfaces CALL v2f%get_ith_conn(iv2f,IF) n = SIZE(iv2f) ! # of face vertices is equal to # of face edges N DO i = 1, n - 1 CALL work%set_ith_conn(ie,(/ iv2f(i), iv2f(i+1) /)) ie = ie + 1 END DO CALL work%set_ith_conn(ie,(/ iv2f(n), iv2f(1) /)) ie = ie + 1 END DO CALL work%get_dual_conn(e2v) ALLOCATE(inspected(nedges),twice(nedges),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF nverts = e2v%nel_() inspected = .FALSE. twice = .FALSE. DO iv = 1, nverts CALL e2v%get_ith_conn(ie2v,iv) n = SIZE(ie2v) edge1: DO i = 1, n - 1 ie1 = ie2v(i) IF(twice(ie1) .OR. inspected(ie1)) CYCLE edge1 inspected(ie1) = .TRUE. CALL work%get_ith_conn(iv2e,ie1) iv1 = iv2e(1) iv2 = iv2e(2) s1 = iv1 + iv2 d1 = ABS(iv1 - iv2) edge2: DO j = i + 1, n ie2 = ie2v(j) IF(twice(ie2) .OR. inspected(ie2)) CYCLE edge2 CALL work%get_ith_conn(iv2e,ie2) iv3 = iv2e(1) iv4 = iv2e(2) s2 = iv3 + iv4 d2 = ABS(iv3 - iv4) IF(s1 == s2 .AND. d1 == d2) twice(ie2) = .TRUE. !!$ if((iv4 == iv1 .and. iv3 == iv2) .or.& !!$ iv4 == iv2 .and. iv3 == iv1 ) then !!$ twice(ie2) = .true. !!$ end if END DO edge2 END DO edge1 END DO nedges = COUNT(.NOT.twice) CALL alloc_conn(v2e,nel=nedges,nconn=nedges*2) ie = 0 n = work%nel_() DO i = 1, n IF(twice(i)) CYCLE ie = ie + 1 CALL work%get_ith_conn(iv2e,i) CALL v2e%set_ith_conn(ie,iv2e) END DO ! Frees storage NULLIFY(iv2e,ie2v,iv2f) DEALLOCATE(inspected,twice) CALL free_conn(e2v) CALL free_conn(work) ENDIF ! Stop timing CALL sw_v2e%toc() ! Synchronization CALL sw_v2e%synchro() IF(mypnum_() == 0) WRITE(*,'(a,es13.6)') & &' - Elapsed time for V2E building:', sw_v2e%partial_() 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2E') END PROCEDURE cmp_mesh_v2e MODULE PROCEDURE cmp_mesh_v2v USE class_psblas USE class_connectivity IMPLICIT NONE !! $Id: cmp_mesh_v2v.f90 3093 2008-04-22 14:51:09Z sfilippo $ !! !! Description: !! Provides the Vertex To Vertex adjacency graph. !! The criterium for associating the graph nodes is arbitrarly defined by the !! user. In this case, using a co-located cell-centered strategy for variable !! displacement, the V2V connectivity is obtained by linking each vertex to all !! the ones belonging to the cells which share it. !! INTEGER, PARAMETER :: ncv_max = 50 ! Max # of cells sharing a vertex INTEGER :: i, ic, inb, info, iv, iv1, iv2, j, n, ncv, ncells INTEGER, POINTER :: iv2c(:) => NULL() INTEGER, ALLOCATABLE :: kconn(:), itab(:,:) LOGICAL :: inspected(nverts) TYPE(stopwatch) :: sw_v2v IF(mypnum_() == 0) WRITE(*,*) 'Computing vertices adjacency graph V2V' ! Constructs stopwatch sw_v2v = stopwatch_(icontxt_()) ! Start timing CALL sw_v2v%tic() n = (v2c%max_conn() - 1) * ncv_max ! Maximum conn. degree (overestimated) ncells = v2c%nel_() ALLOCATE(itab(n,nverts),kconn(nverts),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF kconn = 0 DO ic = 1, ncells CALL v2c%get_ith_conn(iv2c,ic) ncv = SIZE(iv2c) DO i = 1, ncv - 1 iv1 = iv2c(i) DO j = i + 1, ncv iv2 = iv2c(j) ! Add direct arc kconn(iv1) = kconn(iv1) + 1 itab(kconn(iv1),iv1) = iv2 ! Add indirect arc kconn(iv2) = kconn(iv2) + 1 itab(kconn(iv2),iv2) = iv1 END DO END DO END DO ! Eliminates duplicates inspected = .FALSE. DO iv = 1, nverts j = 0 n = kconn(iv) DO i = 1, n inb = itab(i,iv) IF(.NOT.inspected(inb)) THEN j = j + 1 itab(j,iv) = inb inspected(inb) = .TRUE. END IF END DO kconn(iv) = j ! Reset inspected DO i = 1, j inb = itab(i,iv) inspected(inb) = .FALSE. END DO END DO ! Computes total number of connectivities n = 0 DO iv = 1, nverts n = n + kconn(iv) END DO CALL alloc_conn(v2v,nel=nverts,nconn=n) ! Sets V2V connectivity DO iv = 1, nverts n = kconn(iv) CALL v2v%set_ith_conn(iv,itab(1:n,iv)) END DO ! Frees memory NULLIFY(iv2c) DEALLOCATE(kconn,itab) ! Stop timing CALL sw_v2v%toc() ! Synchronization CALL sw_v2v%synchro() IF(mypnum_() == 0) WRITE(*,'(a,es13.6)') & & ' - Elapsed time for V2V building:', sw_v2v%partial_() 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2V') END PROCEDURE cmp_mesh_v2v MODULE PROCEDURE cmp_mesh_v2ve USE class_psblas USE class_connectivity USE tools_mesh, ONLY : cmp_mesh_v2e IMPLICIT NONE !! $Id: cmp_mesh_v2ve.f90 8157 2014-10-09 13:02:44Z sfilippo $ !! !! Description: !! Builds a vertex to vertex connectivity based on edge connections. !! INTEGER :: i, ie, info, iv, iv_nb, nconn, nmax, nverts INTEGER, POINTER :: ie2v(:) => NULL(), iv2e(:) => NULL() INTEGER, ALLOCATABLE :: iv2v(:,:), kconn(:) TYPE(connectivity) :: v2e, e2v TYPE(stopwatch) :: sw_v2v IF(mypnum_() == 0) WRITE(*,*) 'Computing vertices adjacency graph V2VE' ! Constructs stopwatch sw_v2v = stopwatch_(icontxt_()) ! Start timing CALL sw_v2v%tic() SELECT CASE(ncd) CASE(2) IF(mypnum_() == 0) CALL v2f%get_dual_conn(e2v) CASE(3) CALL cmp_mesh_v2e(ncd,v2f,v2e) IF(mypnum_() == 0) CALL v2e%get_dual_conn(e2v) END SELECT IF(mypnum_() == 0) THEN nverts = e2v%nel_() nmax = e2v%max_conn() ALLOCATE(iv2v(nverts,nmax),kconn(nverts),stat=info) IF(info /= 0) THEN WRITE(*,100) CALL abort_psblas END IF SELECT CASE(ncd) CASE(2) ! 2D DO iv = 1, nverts CALL e2v%get_ith_conn(ie2v,iv) nconn = SIZE(ie2v) DO i = 1, nconn ie = ie2v(i) ! In 2D edges = faces => use V2F CALL v2f%get_ith_conn(iv2e,ie) iv_nb = iv2e(1) + iv2e(2) - iv iv2v(iv,i) = iv_nb END DO kconn(iv) = nconn END DO CASE(3) ! 3D DO iv = 1, nverts CALL e2v%get_ith_conn(ie2v,iv) nconn = SIZE(ie2v) DO i = 1, nconn ie = ie2v(i) ! In 3D edges /= faces => use V2E CALL v2e%get_ith_conn(iv2e,ie) iv_nb = iv2e(1) + iv2e(2) - iv iv2v(iv,i) = iv_nb END DO kconn(iv) = nconn END DO END SELECT ! Computes total number of connectivities nconn = 0 DO iv = 1, nverts nconn = nconn + kconn(iv) END DO ! Allocate and sets V2V CALL alloc_conn(v2v,nel=nverts,nconn=nconn) DO iv = 1, nverts nconn = kconn(iv) CALL v2v%set_ith_conn(iv,iv2v(iv,1:nconn)) END DO ! Free memory storage NULLIFY(iv2e,ie2v) DEALLOCATE(iv2v,kconn) CALL free_conn(e2v) IF(ncd == 3) CALL free_conn(v2e) ENDIF ! Stop timing CALL sw_v2v%toc() ! Synchronization CALL sw_v2v%synchro() IF(mypnum_() == 0) WRITE(*,'(a,es13.6)') & &' - Elapsed time for V2VE building:', sw_v2v%partial_() 100 FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2VE') END PROCEDURE cmp_mesh_v2ve END SUBMODULE cmp_mesh_implementation