From 93e3b3db8e63775fdcfb54e18b6d3ade55f5df1c Mon Sep 17 00:00:00 2001 From: Cooper Trucks Date: Thu, 1 Aug 2024 18:16:43 -0500 Subject: [PATCH 1/2] Added ability for physical groups to be recognized as regions in THOR --- src/globals.f90 | 6 ++++++ src/output_thrm.f90 | 7 ++++++- src/read_gmsh.f90 | 46 ++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 55 insertions(+), 4 deletions(-) diff --git a/src/globals.f90 b/src/globals.f90 index d5a82a7..cc64c62 100644 --- a/src/globals.f90 +++ b/src/globals.f90 @@ -37,6 +37,12 @@ MODULE globals !element region tags INTEGER, ALLOCATABLE :: el_tag(:) + !physical volumes + INTEGER, ALLOCATABLE :: phys_vol_tag(:) + + !element physical volume tag + INTEGER, ALLOCATABLE :: el_phys_group(:) + !adjacency list INTEGER, ALLOCATABLE :: adj_list(:,:) diff --git a/src/output_thrm.f90 b/src/output_thrm.f90 index 611c304..d29b326 100644 --- a/src/output_thrm.f90 +++ b/src/output_thrm.f90 @@ -28,7 +28,11 @@ SUBROUTINE output_thrm_file() ENDDO !print out tet regions and source. Assume each region has its own source. User can change this later DO i=1,num_tets - WRITE(30,'(I0,A,I0,A,I0)')i,' ',el_tag(i),' ',el_tag(i) + IF (el_phys_group(i) .EQ. 0) THEN + WRITE(30,'(I0,A,I0,A,I0,A,I0)')i,' ',el_tag(i),' ',el_tag(i) + ELSE + WRITE(30,'(I0,A,I0,A,I0)')i,' ',el_phys_group(i),' ',el_tag(i) + END IF ENDDO !print out tet composition DO i=1,num_tets @@ -101,6 +105,7 @@ SUBROUTINE calcvols() IF(tets_in_reg(i) .GT. 0)THEN WRITE(*,'(A,I0,A,I0)')'Region ',i,' tets: ',tets_in_reg(i) WRITE(*,'(A,I0,A,ES24.16)')'Region ',i,' volume: ',regvol(i) + WRITE(*,'(A,I0,A,I0)')'Region ',i,' physical group: ',phys_vol_tag(i) ENDIF ENDDO WRITE(*,'(A,I0)')'Total number of tets: ',SUM(tets_in_reg) diff --git a/src/read_gmsh.f90 b/src/read_gmsh.f90 index e706b77..6ade2d8 100644 --- a/src/read_gmsh.f90 +++ b/src/read_gmsh.f90 @@ -26,6 +26,18 @@ SUBROUTINE read_gmsh_file() temp_str=TRIM(ADJUSTL(temp_str)) IF(temp_str(1:2) .NE. '4.')STOP 'gmesh format is wrong, should be 4 version' + !find the entities + REWIND(20) + DO + READ(20,*,iostat=ios)temp_str + IF(temp_str .EQ. '$Entities')THEN + !and read them + CALL read_entities() + EXIT + ENDIF + IF(ios .NE. 0)STOP 'did not find Entities in gmesh file' + ENDDO + !find the nodes REWIND(20) DO @@ -78,12 +90,12 @@ SUBROUTINE read_nodes() SUBROUTINE read_elements() INTEGER :: num_entities,numel,minel,maxel - INTEGER :: temp_int,loc_num_el,i,j,el_indx,el_dim,ent_tag + INTEGER :: temp_int,loc_num_el,i,j,el_indx,el_dim,ent_tag,ent_phys INTEGER, ALLOCATABLE :: temp_array(:,:) !get the elements info READ(20,*)num_entities,numel,minel,maxel !last value of array is tag (region/material) - ALLOCATE(temp_array(numel,5)) + ALLOCATE(temp_array(numel,6)) IF(minel .NE. 1 .OR. maxel .NE. numel)STOP 'elements should be indexed 1 to amount' !read in all node data num_tets=0 @@ -95,11 +107,13 @@ SUBROUTINE read_elements() !actually counting the tets IF(el_dim .EQ. 3)THEN num_tets=num_tets+loc_num_el + ent_phys=phys_vol_tag(ent_tag) !get data since this is a tet set DO j=1,loc_num_el el_indx=el_indx+1 READ(20,*)temp_int,temp_array(el_indx,1:4) temp_array(el_indx,5)=ent_tag + temp_array(el_indx,6)=ent_phys IF(MOD(el_indx,CEILING(numel*1.0/(max_prog-1.0))) .EQ. 0)THEN WRITE(*,'(A)',ADVANCE='NO')'*' prog=prog+1 @@ -113,12 +127,13 @@ SUBROUTINE read_elements() ENDIF ENDDO !allocate elements and element tags - ALLOCATE(element(num_tets,4),el_tag(num_tets)) + ALLOCATE(element(num_tets,4),el_tag(num_tets),el_phys_group(num_tets)) element=0 el_tag=0 DO i=1,num_tets element(i,:)=temp_array(i,1:4) el_tag(i)=temp_array(i,5) + el_phys_group(i)=temp_array(i,6) ENDDO DEALLOCATE(temp_array) DO i=prog,max_prog @@ -126,4 +141,29 @@ SUBROUTINE read_elements() ENDDO WRITE(*,*) ENDSUBROUTINE read_elements + + SUBROUTINE read_entities() + INTEGER :: num_points, num_curves, num_surfaces, num_volumes + INTEGER :: num_skiplines, i, j + REAL :: junk(8) + + READ(20,*) num_points, num_curves, num_surfaces, num_volumes + num_skiplines = num_points + num_curves + num_surfaces + ! skip to volume entities + DO i=1,num_skiplines + READ(20,*) + END DO + ! read first physical group of each volume + ALLOCATE(phys_vol_tag(num_volumes)) + phys_vol_tag = 0 + DO i=1,num_volumes + READ(20,*) junk(1:7), j, phys_vol_tag(i) + IF (j==0) THEN + phys_vol_tag(i) = 0 + WRITE(*,'(A,I0,A)')'Region ',i,' has no physical groups associated.' + ELSE IF (j .NE. 1) THEN + WRITE(*,'(A,I0,A)')'Warning: Region ',i,' has more than one physical group associated!' + END IF + END DO + ENDSUBROUTINE END MODULE read_gmsh From bef96e2b3f9e4b3c342b3f4e09c315a29a679c13 Mon Sep 17 00:00:00 2001 From: Cooper Trucks Date: Tue, 6 Aug 2024 11:40:36 -0500 Subject: [PATCH 2/2] Updated README.md to explain physical group assignment to THOR --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index ca5e024..28c6952 100644 --- a/README.md +++ b/README.md @@ -33,6 +33,8 @@ It should be noted that if reflective boundary conditions are specified, then th If the user tries to assign reflective boundary conditions to a direction with a non-flat boundary, then the converter utility will throw an error and terminate. This check is ignored if all boundary conditions are reflective. +If physical volumes are specified by the Gmsh input file, they will be used as the region indices for elements in the THOR mesh. Elements should not have more than one physical volume assigned to them, as only one physical volume will be read in by the mesh converter per element. Assignment of physical volumes *by material* is recommended in order to simplify region mapping in the THOR input file. + --- ## Examples ---