Support named mesh element groups for BC

Allow the use of named groups (labels) from the mesh file to
apply BC and/or material ID, instead of numbers (tag).

Using labels is still allowed for backwards compatibility.
If a group has a label, either of the tag or the label can be used to
reference it (but not both).

Add error checks for:

- Inexistent labels (i.e. trying to apply a BC onto a label not defined
  in the mesh file).
- No groups defined for material ID assignment.
  - Groups must be created; they may or may not be named.

Use the label for printing loadcase information, when appropriate.

Update/improve some loadcase/boundary condition error messages.

Author: Javier Velo

PRIVATE

@@ -508,6 +508,8 @@ subroutine IO_error_new(error_ID, &
 ! errors related to the mesh solver
     case (800)
       msg = 'invalid mesh'
+    case (812)
+      msg = 'invalid boundary conditions'
 
 !-------------------------------------------------------------------------------------------------
 ! errors related to the grid solver

src/IO.f90

@@ -11,6 +11,7 @@ program DAMASK_mesh
   use PetscDM
   use prec
   use CLI
+  use IO
   use parallelization
   use materialpoint
   use discretization_mesh
@@ -40,10 +41,11 @@ program DAMASK_mesh
   logical :: &
     guess, &                                                                                        !< guess along former trajectory
     stagIterate
+  logical, allocatable, dimension(:) :: &
+    read_BC_entries
   integer :: &
     l, &
     m, &
-    errorID, &
     cutBackLevel = 0, &                                                                             !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
     stepFraction = 0, &                                                                             !< fraction of current time interval
     inc, &                                                                                          !< current increment in current load case
@@ -64,32 +66,29 @@ program DAMASK_mesh
     mech_u, &
     step_mech
   character(len=pSTRLEN) :: &
-    incInfo, &
-    loadcase_string
+    incInfo
   integer :: &
+    bc_tag, &
     stagItMax, &                                                                                    !< max number of field level staggered iterations
     maxCutBack, &                                                                                   !< max number of cutbacks
-    skipBCTagDigits, &                                                                              !< number of characters to skip to print BC tag (T descriptor)
-    BCTag                                                                                           !< tag value read from YAML load file
-  integer, allocatable, dimension(:) :: &
-    knownBCTags                                                                                     !< array of read BC tags (check for duplicates)
-
+    skip_T1, skip_T2                                                                                !< number of characters to skip (T descriptor)
 
   type(tLoadCase), allocatable, dimension(:) :: loadCases                                           !< array of all load cases
   type(tSolutionState), allocatable, dimension(:) :: solres
   PetscInt :: boundary, dimPlex
   PetscErrorCode :: err_PETSc
   character(len=:), allocatable :: &
-    fileContent, fname, tagPrintFormat
-  character(len=6) :: BC_elem
+    fileContent, fname
+  character(len=pSTRLEN) :: &
+    bc_label_name
 
 
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
   call materialpoint_initAll()
   print'(/,1x,a)', '<<<+-  DAMASK_mesh init  -+>>>'; flush(IO_STDOUT)
 
-!---------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 ! reading field information from numerics file and do sanity checks
   num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
   num_mesh   => num_solver%get_dict('mesh',defaultVal=emptyDict)
@@ -99,8 +98,10 @@ program DAMASK_mesh
   if (stagItMax < 0)  call IO_error(301,ext_msg='N_staggered_iter_max')
   if (maxCutBack < 0) call IO_error(301,ext_msg='N_cutback_max')
 
-! reading basic information from load case file and allocate data structure containing load cases
-  call DMGetDimension(geomMesh,dimPlex,err_PETSc)                                                   !< dimension of mesh (2D or 3D)
+!--------------------------------------------------------------------------------------------------
+! reading basic information from load case file, allocate data structure containing load cases
+! and some checks for invalid tags/labels use
+  call DMGetDimension(geomMesh,dimPlex,err_PETSc)                                                   ! dimension of mesh (2D or 3D)
   CHKERRA(err_PETSc)
   allocate(solres(1))
 
@@ -117,31 +118,45 @@ program DAMASK_mesh
   load => YAML_str_asDict(fileContent)
   load_steps => load%get_list('loadstep')
 
+  allocate(read_BC_entries(mesh_nBoundaries), source = .false.)
   allocate(loadCases(size(load_steps)))
   do l = 1, size(load_steps)
     load_step => load_steps%get_dict(l)
     step_bc   => load_step%get_dict('boundary_conditions')
     step_mech => step_bc%get_list('mechanical')
     allocate(loadCases(l)%mechBC(mesh_Nboundaries))
-    loadCases(l)%mechBC(:)%nComponents = dimPlex                                                   !< X, Y (, Z) displacements
+    loadCases(l)%mechBC(:)%nComponents = dimPlex                                                    ! X, Y (, Z) displacements
     do boundary = 1_pPETSCINT, mesh_Nboundaries
       allocate(loadCases(l)%mechBC(boundary)%dot_u(dimPlex),  source = 0.0_pREAL)
       allocate(loadCases(l)%mechBC(boundary)%active(dimPlex), source = .false.)
     end do
 
-    allocate(knownBCTags(size(step_mech)), source = -1)
     do m = 1, size(step_mech)
       mech_BC => step_mech%get_dict(m)
-      BCTag = mech_BC%get_asInt('tag')
-      boundary = findloc(mesh_boundariesIS, BCtag, dim = 1)
-      if (boundary == 0) then                                                                      !< tag not defined in mesh file
-        call IO_error(error_ID = 837, ext_msg = 'BC tag '//mech_BC%get_asStr('tag')// &
-                                                ' refers to inexistent node/edge/face mesh group')
-      else if (findloc(knownBCTags, BCTag, dim = 1) /= 0) then                                     !< duplicated tag
-        call IO_error(error_ID = 837, ext_msg = 'BC redefinition not allowed (tag '// &
-                                                mech_BC%get_asStr('tag')//')')
+      if (mech_BC%contains('label')) then
+        bc_label_name = mech_BC%get_asStr('label')
+        boundary = findloc(mesh_BCLabels, bc_label_name, dim = 1)
+        if (boundary == 0) &                                                                        ! label not defined in mesh file
+          call IO_error(812_pI16, 'label', trim(bc_label_name), 'not defined', emph = [2])
+        if (read_BC_entries(boundary)) &                                                            ! duplicated label/tag
+          call IO_error(812_pI16, 'duplicated entries: label', trim(bc_label_name), "and tag", &
+                        mesh_boundariesIS(boundary), emph = [2,4])
+        read_BC_entries(boundary) = .true.
+        loadCases(l)%mechBC(boundary)%use_label = .true.
+      else if (mech_BC%contains('tag')) then
+        bc_tag = mech_BC%get_asInt('tag')
+        boundary = findloc(mesh_boundariesIS, bc_tag, dim = 1)
+        if (boundary == 0) &                                                                        ! tag not defined in mesh file
+          call IO_error(812_pI16, 'tag', bc_tag, 'not defined', emph = [2])
+        if (read_BC_entries(boundary)) &                                                            ! duplicated tag/label
+          call IO_error(812_pI16, 'duplicated entries: tag', bc_tag, 'and label', &
+                        trim(mesh_BCLabels(boundary)), emph = [2, 4])
+        read_BC_entries(boundary) = .true.
+        loadCases(l)%mechBC(boundary)%use_label = .false.
+      else
+        call IO_error(837_pI16, 'entry'//IO_intAsStr(m)//': at least one of   &
+                      & "tag/label"   required')
       end if
-      knownBCTags(m) = BCTag
       mech_u => mech_BC%get_list('dot_u')
       do component = 1, dimPlex
         if (mech_u%get_asStr(component) /= 'x') then
@@ -150,7 +165,7 @@ program DAMASK_mesh
         end if
       end do
     end do
-    deallocate(knownBCTags)
+    read_BC_entries = .false.
 
     step_discretization => load_step%get_dict('discretization')
     loadCases(l)%t = step_discretization%get_asReal('t')
@@ -163,38 +178,45 @@ program DAMASK_mesh
     end if
     loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
   end do
+  deallocate(read_BC_entries)
 
 !--------------------------------------------------------------------------------------------------
-! consistency checks and output of load case
-  errorID = 0
-  skipBCTagDigits = 4+6+1+floor(log10(real(maxval(mesh_boundariesIS))))+1+2                        !< Indentation(4)+BC_elem(6)+blank(1)+descriptor T(1)+NumDigits(floor(..)+1)+2(start printing after 1 blank)
-  allocate(character(len=skipBCTagDigits) :: tagPrintFormat)
-  write(tagPrintFormat,'(i0)') skipBCTagDigits
+! loadcase parameters checks and output of load case information
+  skip_T1 = 4+max(len(PETSC_GENERIC_LABELS), maxval(len_trim(mesh_BCLabels)))+2                     ! indentation(4)+length_longest_label+blank
+  skip_T2 = skip_T1+1+floor(log10(real(maxval(mesh_boundariesIS))))+1+1+1                           ! T1+"("+NumDigits(floor()+1)+")"+blank
   checkLoadcases: do l = 1, size(load_steps)
-    write (loadcase_string, '(i0)' ) l
+    if (loadCases(l)%N < 1) &
+      call IO_error(813_pI16, 'loadcase', IO_intAsStr(l)//':', 'N', &
+                    'must be positive', emph = [3])
+    if (loadCases(l)%f_out < 1) &
+      call IO_error(813_pI16, 'loadcase', IO_intAsStr(l)//':', 'f_out', &
+                    'must be positive', emph = [3])
+
     print'(/,1x,a,1x,i0)', 'load case:', l
-    if (.not. loadCases(l)%estimate_rate) &
-      print'(2x,a)', 'drop guessing along trajectory'
+    if (.not. loadCases(l)%estimate_rate) print'(2x,a)', 'drop guessing along trajectory'
     print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
 
     do boundary = 1_pPETSCINT, mesh_Nboundaries
-      BC_elem = merge('Vertex', merge('Edge  ', 'Face  ', dimplex == 2_pPETSCINT), &
-                      mesh_boundariesIdx(boundary) == mesh_BCTypeVertex)
+      if (loadCases(l)%mechBC(boundary)%use_label) then
+        bc_label_name = trim(mesh_BCLabels(boundary))
+      else
+        m = mesh_boundariesIdx(boundary)
+        if (dimPlex == 2_pPETSCINT .and. m < 4) m = m + 1
+        bc_label_name = PETSC_GENERIC_LABELS(m)
+      end if
       do component = 1_pPETSCINT, dimPlex
         if (loadCases(l)%mechBC(boundary)%active(component)) &
-          print'(4x,a,1x,i0,T'//tagPrintFormat//',a,1x,i1,1x,a,1x,f12.7)', &
-          BC_elem, mesh_boundariesIS(boundary), &
-          'Component', component, &
-          'Value', loadCases(l)%mechBC(boundary)%dot_u(component)
+          print'(4x,a,T'//IO_intAsStr(skip_T1)//',a,i0,a,'// &
+                'T'//IO_intAsStr(skip_T2)//',a,1x,i1,1x,a,1x,f12.7)', &
+            trim(bc_label_name), "(", mesh_boundariesIS(boundary), ")", &
+            'Component', component, &
+            'Value', loadCases(l)%mechBC(boundary)%dot_u(component)
       end do
     end do
 
-    print'(2x,a,T21,g0.6)', 'time:', loadCases(l)%t
-    if (loadCases(l)%N < 1) errorID = 835                                                          ! non-positive incs count
-    print'(2x,a,T21,i0)',   'increments:', loadCases(l)%N
-    if (loadCases(l)%f_out < 1) errorID = 836                                                      ! non-positive result frequency
+    print'(2x,a,T21,g0.6)', 'time:',             loadCases(l)%t
+    print'(2x,a,T21,i0)',   'increments:',       loadCases(l)%N
     print'(2x,a,T21,i0)',   'output frequency:', loadCases(l)%f_out
-    if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string)                  ! exit with error message
   end do checkLoadcases
 
 !--------------------------------------------------------------------------------------------------
@@ -213,8 +235,8 @@ program DAMASK_mesh
   call materialpoint_result(0,0.0_pREAL)
 
   loadCaseLooping: do l = 1, size(load_steps)
-    t_0 = t                                                                                        ! load case start time
-    guess = loadCases(l)%estimate_rate                                                             ! change of load case? homogeneous guess for the first inc
+    t_0 = t                                                                                         ! load case start time
+    guess = loadCases(l)%estimate_rate                                                              ! change of load case? homogeneous guess for the first inc
 
     incLooping: do inc = 1, loadCases(l)%N
       totalIncsCounter = totalIncsCounter + 1
@@ -292,7 +314,7 @@ program DAMASK_mesh
         print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged'
       end if; flush(IO_STDOUT)
 
-      if (mod(inc,loadCases(l)%f_out) == 0) then                            ! at output frequency
+      if (mod(inc,loadCases(l)%f_out) == 0) then                                                    ! at output frequency
         print'(/,1x,a)', '... saving results ........................................................'
         call FEM_mechanical_updateCoords()
         call materialpoint_result(totalIncsCounter,t)

src/mesh/DAMASK_mesh.f90

@@ -50,6 +50,7 @@ module FEM_utilities
     integer :: nComponents = 0
     real(pREAL), allocatable, dimension(:) :: dot_u
     logical,     allocatable, dimension(:) :: active
+    logical                                :: use_label
   end type tMechBC
 
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
@@ -78,11 +79,11 @@ subroutine FEM_utilities_init(num_mesh)
   type(tDict), pointer :: &
     num_mech
   character(len=:), allocatable :: &
-    petsc_options
+    PETSc_options
   integer :: &
-    p_s, &                                                                                          !< order of shape functions
-    p_i                                                                                             !< integration order (quadrature rule)
-  PetscErrorCode  :: err_PETSc
+    p_s, &                                                                                          ! order of shape functions
+    p_i                                                                                             ! integration order (quadrature rule)
+  PetscErrorCode :: err_PETSc
 
 
   print'(/,1x,a)',   '<<<+-  FEM_utilities init  -+>>>'
@@ -151,61 +152,62 @@ end subroutine utilities_constitutiveResponse
 !--------------------------------------------------------------------------------------------------
 !> @brief Project BC values to local vector
 !--------------------------------------------------------------------------------------------------
-subroutine utilities_projectBCValues(dm_local,solution_local,section,mechBC,Delta_t,dimPlex)
+subroutine utilities_projectBCValues(dm_local,solution_local_vec,section,mechBC,Delta_t,dimPlex)
   DM  :: dm_local
-  Vec :: solution_local
+  Vec :: solution_local_vec
   PetscSection   :: section
   type(tMechBC),  dimension(:), intent(in) :: &
     mechBC
   real(pREAL),    intent(in) :: Delta_t
   PetscInt, intent(in) :: dimPlex
 
-  PetscInt       :: component, boundary, bcSize, &
-                    nBcPoints, point, dof, numDof, numComp, offset
-  IS             :: bcPointsIS
+  PetscInt       :: component, boundary, bc_n_points, &
+                    points_IS_size, point, dof, n_field_dof, n_field_cmp, offset
+  IS             :: bc_points_IS
   PetscErrorCode :: err_PETSc
-  PetscInt,    pointer :: bcPoints(:)
-  real(pREAL), pointer :: localArray(:)
-
-  character(len=11) :: setLabel
-
-
-  do boundary = 1_pPETSCINT, mesh_Nboundaries; do component = 1_pPETSCINT, dimPlex
-   if (mechBC(boundary)%active(component)) then
-     setLabel = mesh_BCTypeLabel(mesh_boundariesIdx(boundary))
-     call DMGetStratumSize(dm_local,setLabel,mesh_boundariesIS(boundary),bcSize,err_PETSc)
-     if (bcSize > 0_pPETSCINT) then
-       call DMGetStratumIS(dm_local,setLabel,mesh_boundariesIS(boundary),bcPointsIS,err_PETSc)
-       CHKERRQ(err_PETSc)
-       call PetscSectionGetFieldComponents(section,0_pPETSCINT,numComp,err_PETSc)
-       CHKERRQ(err_PETSc)
-       call ISGetSize(bcPointsIS,nBcPoints,err_PETSc)
-       CHKERRQ(err_PETSc)
-       if (nBcPoints > 0_pPETSCINT) call ISGetIndices(bcPointsIS,bcPoints,err_PETSc)
-       call VecGetArray(solution_local,localArray,err_PETSc)
-       CHKERRQ(err_PETSc)
-       do point = 1_pPETSCINT, nBcPoints
-         call PetscSectionGetFieldDof(section,bcPoints(point),0_pPETSCINT,numDof,err_PETSc)
-         CHKERRQ(err_PETSc)
-         call PetscSectionGetFieldOffset(section,bcPoints(point),0_pPETSCINT,offset,err_PETSc)
-         CHKERRQ(err_PETSc)
-         do dof = offset+component, offset+numDof, numComp
-           localArray(dof) = localArray(dof) + mechBC(boundary)%dot_u(component)*Delta_t
-         end do
-       end do
-       call VecRestoreArray(solution_local,localArray,err_PETSc)
-       CHKERRQ(err_PETSc)
-       call VecAssemblyBegin(solution_local, err_PETSc)
-       CHKERRQ(err_PETSc)
-       call VecAssemblyEnd(solution_local, err_PETSc)
-       CHKERRQ(err_PETSc)
-       if (nBcPoints > 0_pPETSCINT) call ISRestoreIndices(bcPointsIS,bcPoints,err_PETSc)
-       CHKERRQ(err_PETSc)
-       call ISDestroy(bcPointsIS,err_PETSc)
-       CHKERRQ(err_PETSc)
-     end if
-   end if
-  end do; end do
+  PetscInt,    pointer   :: bc_points(:)
+  real(pREAL), pointer   :: solution_local(:)
+  character(len=pSTRLEN) :: bc_label
+
+
+  do boundary = 1, mesh_nBoundaries;
+    bc_label = PETSC_GENERIC_LABELS(mesh_boundariesIdx(boundary))
+    do component = 1_pPETSCINT, dimPlex
+      if (mechBC(boundary)%active(component)) then
+        call DMGetStratumSize(dm_local,bc_label,mesh_boundariesIS(boundary),bc_n_points,err_PETSc)
+        if (bc_n_points > 0_pPETSCINT) then
+          call DMGetStratumIS(dm_local,bc_label,mesh_boundariesIS(boundary),bc_points_IS,err_PETSc)
+          CHKERRQ(err_PETSc)
+          call PetscSectionGetFieldComponents(section,0_pPETSCINT,n_field_cmp,err_PETSc)
+          CHKERRQ(err_PETSc)
+          call ISGetSize(bc_points_IS,points_IS_size,err_PETSc)
+          CHKERRQ(err_PETSc)
+          if (points_IS_size > 0_pPETSCINT) call ISGetIndices(bc_points_IS,bc_points,err_PETSc)
+          call VecGetArray(solution_local_vec,solution_local,err_PETSc)
+          CHKERRQ(err_PETSc)
+          do point = 1_pPETSCINT, points_IS_size
+            call PetscSectionGetFieldDof(section,bc_points(point),0_pPETSCINT,n_field_dof,err_PETSc)
+            CHKERRQ(err_PETSc)
+            call PetscSectionGetFieldOffset(section,bc_points(point),0_pPETSCINT,offset,err_PETSc)
+            CHKERRQ(err_PETSc)
+            do dof = offset+component, offset+n_field_dof, n_field_cmp
+              solution_local(dof) = solution_local(dof) + mechBC(boundary)%dot_u(component)*Delta_t
+            end do
+          end do
+          call VecRestoreArray(solution_local_vec,solution_local,err_PETSc)
+          CHKERRQ(err_PETSc)
+          call VecAssemblyBegin(solution_local_vec, err_PETSc)
+          CHKERRQ(err_PETSc)
+          call VecAssemblyEnd(solution_local_vec, err_PETSc)
+          CHKERRQ(err_PETSc)
+          if (points_IS_size > 0_pPETSCINT) call ISRestoreIndices(bc_points_IS,bc_points,err_PETSc)
+          CHKERRQ(err_PETSc)
+          call ISDestroy(bc_points_IS,err_PETSc)
+          CHKERRQ(err_PETSc)
+        end if
+      end if
+    end do
+  end do
 
 end subroutine utilities_projectBCValues