ip volume is now based on the determinant of F. "mesh_ipVolume" represents the initial volume and is multiplied with det(F) wherever the current volume is needed. Since this works for all solver types, the "volume" output in crystallite is now also correct for spectral method and abaqus.

code/CPFEM.f90

@@ -533,10 +533,9 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
               write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
             !$OMP END CRITICAL (write2out)
           endif
-          if (FEsolver == 'Marc') then                                    ! marc updates nodal coordinates, whereas Abaqus and spectral solver directly update ip coordinates. In the latter case it is not possible to get the current ip volume, since the current nodal positions are unknown
+          if (FEsolver == 'Marc') then                                    ! marc returns nodal coordinates, whereas Abaqus and spectral solver return ip coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates.
             call mesh_build_subNodeCoords()                               ! update subnodal coordinates
             call mesh_build_ipCoordinates()                               ! update ip coordinates
-            call mesh_build_ipVolumes()                                   ! update ip volumes
           endif
           if (debug_verbosity > 0) then
             !$OMP CRITICAL (write2out)

code/constitutive_nonlocal.f90

@@ -1008,7 +1008,7 @@ if (.not. phase_localConstitution(phase)) then
       neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance)
       neighboring_ns = constitutive_nonlocal_totalNslip(neighboring_instance)
       neighboring_invFe = math_inv3x3(Fe(1:3,1:3,1,neighboring_ip,neighboring_el))
-      neighboring_ipVolumeSideLength = mesh_ipVolume(neighboring_ip,neighboring_el) ** (1.0_pReal/3.0_pReal)
+      neighboring_ipVolumeSideLength = mesh_ipVolume(neighboring_ip,neighboring_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here
       forall (s = 1:neighboring_ns, c = 1:2) &
         neighboring_rhoExcess(c,1,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*neighboring_ns+s) &  ! positive mobiles
                                      - state(g,neighboring_ip,neighboring_el)%p((2*c-1)*neighboring_ns+s)    ! negative mobiles
@@ -1178,7 +1178,7 @@ if (.not. phase_localConstitution(phase)) then
                 sigma = sigma * constitutive_nonlocal_Gmod(neighboring_instance) &
                               * constitutive_nonlocal_burgersPerSlipSystem(s,neighboring_instance) &
                               / (4.0_pReal * pi * (1.0_pReal - nu)) &
-                              * mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength
+                              * mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength      ! reference volume is used here (according to the segment length calculation)
                 Tdislo_neighboringLattice = Tdislo_neighboringLattice &
                       + math_mul33x33(math_transpose3x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)), &
                         math_mul33x33(sigma, constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)))
@@ -1678,7 +1678,7 @@ forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el)
 !*** check CFL condition for flux
 
 if (any(1.2_pReal * constitutive_nonlocal_v(1:ns,1:4,g,ip,el) * timestep &                    ! security factor 1.2
-        > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then
+        > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then                          ! reference volume and area
   dotState%p = NaN(3)
   return
 endif

code/crystallite.f90

@@ -3063,6 +3063,7 @@ function crystallite_postResults(&
                                       math_QuaternionToAxisAngle, &
                                       math_mul33x33, &
                                       math_transpose3x3, &
+                                      math_det3x3, &
                                       math_I3, &
                                       inDeg, &
                                       math_Mandel6to33
@@ -3090,6 +3091,7 @@ function crystallite_postResults(&
  
  !*** local variables ***!
  real(pReal), dimension(3,3) ::       Ee
+ real(pReal)                          detF
  integer(pInt)                        o,c,crystID,mySize
 
  crystID = microstructure_crystallite(mesh_element(4,e))
@@ -3110,7 +3112,8 @@ function crystallite_postResults(&
        crystallite_postResults(c+1) = material_texture(g,i,e)                   ! textureID of grain
      case ('volume')
        mySize = 1_pInt
-       crystallite_postResults(c+1) = mesh_ipVolume(i,e) / homogenization_Ngrains(mesh_element(3,e)) ! grain volume (not fraction but absolute)
+       detF = math_det3x3(crystallite_partionedF(1:3,1:3,g,i,e))                ! V_current = det(F) * V_reference
+       crystallite_postResults(c+1) = detF * mesh_ipVolume(i,e) / homogenization_Ngrains(mesh_element(3,e)) ! grain volume (not fraction but absolute)
      case ('orientation')
        mySize = 4_pInt
        crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,g,i,e)    ! grain orientation as quaternion