Add support for concentrated loads

Author: Javier Velo

PRIVATE

@@ -41,7 +41,8 @@ program DAMASK_mesh
     Delta_t_prev = 0.0_pREAL                                                                        !< previous time interval
   logical :: &
     guess, &                                                                                        !< guess along former trajectory
-    stagIterate
+    stagIterate, &
+    printed_bc_type
   logical, allocatable, dimension(:) :: &
     read_BC_entries
   integer :: &
@@ -64,15 +65,16 @@ program DAMASK_mesh
     step_discretization
   type(tList), pointer :: &
     load_steps, &
-    mech_u, &
+    mech_u => NULL(), &
+    mech_f => NULL(), &
     step_mech
   character(len=pSTRLEN) :: &
     incInfo
   integer :: &
-    bc_tag, &
     stagItMax, &                                                                                    !< max number of field level staggered iterations
     maxCutBack, &                                                                                   !< max number of cutbacks
-    skip_T1, skip_T2                                                                                !< number of characters to skip (T descriptor)
+    bc_tag, &                                                                                       !< tag used for BC in YAML load file
+    skip_T1                                                                                         !< number of characters to skip (T descriptor)
 
   type(tLoadCase), allocatable, dimension(:) :: loadCases                                           !< array of all load cases
   type(tSolutionState), allocatable, dimension(:) :: solres
@@ -81,7 +83,7 @@ program DAMASK_mesh
   character(len=:), allocatable :: &
     fileContent, fname
   character(len=pSTRLEN) :: &
-    bc_label_name
+    bc_label                                                                                        !< mesh label or generic label for BC
 
 
 !--------------------------------------------------------------------------------------------------
@@ -126,21 +128,24 @@ program DAMASK_mesh
     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
     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.)
+      allocate(loadCases(l)%mechBC(boundary)%dot_u(dimPlex), source = 0.0_pREAL)
+      allocate(loadCases(l)%mechBC(boundary)%dot_u_active(dimPlex), source = .false.)
+      allocate(loadCases(l)%mechBC(boundary)%dot_f(dimPlex), source = 0.0_pREAL)
+      allocate(loadCases(l)%mechBC(boundary)%dot_f_active(dimPlex), source = .false.)
     end do
 
+!--------------------------------------------------------------------------------------------------
+! check valid tags/labels
     do m = 1, size(step_mech)
       mech_BC => step_mech%get_dict(m)
       if (mech_BC%contains('label')) then
-        bc_label_name = mech_BC%get_asStr('label')
-        boundary = findloc(mesh_BCLabels, bc_label_name, dim = 1)
+        bc_label = mech_BC%get_asStr('label')
+        boundary = findloc(mesh_bcLabels, bc_label, 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])
+          call IO_error(812_pI16, 'label', trim(bc_label), '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', &
+          call IO_error(812_pI16, 'duplicated entries: label', trim(bc_label), 'and tag', &
                         mesh_boundariesIS(boundary), emph = [2,4])
         read_BC_entries(boundary) = .true.
         loadCases(l)%mechBC(boundary)%use_label = .true.
@@ -151,25 +156,50 @@ program DAMASK_mesh
           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])
+                        trim(mesh_bcLabels(boundary)), emph = [2, 4])
         read_BC_entries(boundary) = .true.
         loadCases(l)%mechBC(boundary)%use_label = .false.
       else
         call IO_error(812_pI16, 'neither "label" nor "tag" given for boundary condition', m, emph=[2])
       end if
-      mech_u => mech_BC%get_list('dot_u')
-      do component = 1, dimPlex
-        if (mech_u%get_asStr(component) /= 'x') then
-          loadCases(l)%mechBC(boundary)%active(component) = .true.
-          loadCases(l)%mechBC(boundary)%dot_u(component)  = mech_u%get_asReal(component)
-        end if
-      end do
+! check valid BC definition of dot_u and dot_f, store values
+      if (mech_BC%contains('dot_u')) mech_u => mech_BC%get_list('dot_u')
+      if (mech_BC%contains('u_dot')) mech_u => mech_BC%get_list('u_dot')
+      if (mech_BC%contains('dot_f')) mech_f => mech_BC%get_list('dot_f')
+      if (mech_BC%contains('f_dot')) mech_f => mech_BC%get_list('f_dot')
+
+      if (.not. (associated(mech_u) .or. associated(mech_f))) &
+        call IO_error(812_pI16, 'dot_u/dot_f missing for boundary condition', m, emph = [2])
+
+      associate (BC_mech => loadCases(l)%mechBC(boundary))
+        do component = 1, dimplex
+          if (associated(mech_u)) then
+            if (mech_u%get_asStr(component) /= 'x') then
+              BC_mech%dot_u(component) = mech_u%get_asReal(component)
+              BC_mech%dot_u_active(component) = .true.
+            end if
+          end if
+          if (associated(mech_f)) then
+            if (mech_f%get_asStr(component) /= 'x') then
+              BC_mech%dot_f(component) = mech_f%get_asReal(component)
+              BC_mech%dot_f_active(component) = .true.
+            end if
+          end if
+
+          if (BC_mech%dot_f_active(component) .and. BC_mech%dot_u_active(component)) &
+            call IO_error(812_pI16, 'displacement and force prescribed in the same DOF', &
+                          IO_EOL, 'in boundary condition', m, emph = [4])
+        end do
+      end associate
+      nullify(mech_u, mech_f)
     end do
     read_BC_entries = .false.
 
+!--------------------------------------------------------------------------------------------------
+! store discretization, time and frequency
     step_discretization => load_step%get_dict('discretization')
     loadCases(l)%t = step_discretization%get_asReal('t')
-    loadCases(l)%N = step_discretization%get_asInt  ('N')
+    loadCases(l)%N = step_discretization%get_asInt('N')
 
     if (load_step%get_asStr('f_out',defaultVal='n/a') == 'none') then
       loadCases(l)%f_out = huge(0)
@@ -182,39 +212,59 @@ program DAMASK_mesh
 
 !--------------------------------------------------------------------------------------------------
 ! 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
+  skip_T1 = 4+max(len(PETSC_GENERIC_LABELS), maxval(len_trim(mesh_bcLabels)))+2                     ! indentation(4)+length_longest_label+blank
   checkLoadcases: do l = 1, size(load_steps)
     if (loadCases(l)%N < 1) &
-      call IO_error(813_pI16, 'loadcase', l, 'has non-positive number of steps ("N")', emph = [2])
+      call IO_error(301_pI16, 'loadcase', l, 'has non-positive number of steps ("N")', emph = [2])
     if (loadCases(l)%f_out < 1) &
-      call IO_error(813_pI16, 'loadcase', l, 'has non-positive output frequency ("f_out")', emph = [2])
+      call IO_error(301_pI16, 'loadcase', l, 'has non-positive output frequency ("f_out")', emph = [2])
 
     print'(/,1x,a,1x,i0)', 'load case:', l
     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
       if (loadCases(l)%mechBC(boundary)%use_label) then
-        bc_label_name = trim(mesh_BCLabels(boundary))
+        bc_label = 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)
+        if (dimPlex == 2_pPETSCINT .and. m < size(PETSC_GENERIC_LABELS)) m = m + 1                  ! adjust for 2D (cells -> faces)
+        bc_label = PETSC_GENERIC_LABELS(m)
       end if
-      do component = 1_pPETSCINT, dimPlex
-        if (loadCases(l)%mechBC(boundary)%active(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
-    print'(2x,a,T21,i0)',   'increments:',       loadCases(l)%N
-    print'(2x,a,T21,i0)',   'output frequency:', loadCases(l)%f_out
+      printed_bc_type = .false.
+      if (count(loadCases(l)%mechBC(boundary)%dot_u_active) > 0) then
+        if (.not. printed_bc_type) then
+          print'(3x,a,T'//IO_intAsStr(skip_T1)//',a,i0,a)', &
+            trim(bc_label), '(', mesh_boundariesIS(boundary),')'
+          printed_bc_type = .true.
+        end if
+        print'(4x,a)', 'Rates of Displacement'
+        do component = 1_pPETSCINT, dimPlex
+          if (loadCases(l)%mechBC(boundary)%dot_u_active(component)) &
+            print'(5x,a,1x,i1,a,1x,en12.3e2,2x,a)', &
+            'Component', component, ':', &
+            loadCases(l)%mechBC(boundary)%dot_u(component), 'm/s'
+        end do
+      end if
+      if (count(loadCases(l)%mechBC(boundary)%dot_f_active) > 0) then
+        if (.not. printed_bc_type) then
+          print'(3x,a,T'//IO_intAsStr(skip_T1)//',a,i0,a)', &
+            trim(bc_label), '(', mesh_boundariesIS(boundary), ')'
+          printed_bc_type = .true.
+        end if
+        print'(4x,a)', 'Rates of Force'
+        do component = 1_pPETSCINT, dimPlex
+          if (loadCases(l)%mechBC(boundary)%dot_f_active(component)) &
+            print'(5x,a,1x,i1,a,1x,en12.3e2,2x,a)', &
+            'Component', component, ':', &
+            loadCases(l)%mechBC(boundary)%dot_f(component), 'N/s'
+        end do
+      end if
+    end do
+    print'(2x,a,T19,en12.3e2,2x,a)', 'time:',             loadCases(l)%t, 's'
+    print'(2x,a,T22,i0)',            'increments:',       loadCases(l)%N
+    print'(2x,a,T22,i0)',            'output frequency:', loadCases(l)%f_out
   end do checkLoadcases
 
 !--------------------------------------------------------------------------------------------------
@@ -317,10 +367,7 @@ program DAMASK_mesh
         call FEM_mechanical_updateCoords()
         call materialpoint_result(totalIncsCounter,t)
       end if
-
-
     end do incLooping
-
   end do loadCaseLooping