From b300cc7faa6f5e5a977bfed02c6da0513dfa003b Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Fri, 19 Aug 2011 11:18:56 +0000 Subject: [PATCH] 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 | 3 +-- code/constitutive_nonlocal.f90 | 6 +++--- code/crystallite.f90 | 5 ++++- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index ef19b2d1b..cb97ad9fc 100644 --- a/code/CPFEM.f90 +++ b/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) diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index f676f391c..d258c8a26 100644 --- a/code/constitutive_nonlocal.f90 +++ b/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 diff --git a/code/crystallite.f90 b/code/crystallite.f90 index eaee146e2..2abee3d54 100644 --- a/code/crystallite.f90 +++ b/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