diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 66ed96df7..d30f2d3a3 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3478,7 +3478,7 @@ function crystallite_postResults(dt, gr, ip, el) crystallite_postResults = 0.0_pReal c = 0_pInt - crystallite_postResults(c+1) = real(crystallite_sizePostResults(crystID),pReal) ! size of results from cryst + crystallite_postResults(c+1) = real(crystallite_sizePostResults(crystID),pReal) ! size of results from cryst c = c + 1_pInt do o = 1_pInt,crystallite_Noutput(crystID) @@ -3486,52 +3486,55 @@ function crystallite_postResults(dt, gr, ip, el) select case(crystallite_output(o,crystID)) case ('phase') mySize = 1_pInt - crystallite_postResults(c+1) = real(material_phase(gr,ip,el),pReal) ! phaseID of grain + crystallite_postResults(c+1) = real(material_phase(gr,ip,el),pReal) ! phaseID of grain case ('texture') mySize = 1_pInt - crystallite_postResults(c+1) = real(material_texture(gr,ip,el),pReal) ! textureID of grain + crystallite_postResults(c+1) = real(material_texture(gr,ip,el),pReal) ! textureID of grain case ('volume') mySize = 1_pInt - detF = math_det33(crystallite_partionedF(1:3,1:3,gr,ip,el)) ! V_current = det(F) * V_reference - crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) / homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute) + detF = math_det33(crystallite_partionedF(1:3,1:3,gr,ip,el)) ! V_current = det(F) * V_reference + crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) / & + homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute) case ('orientation') mySize = 4_pInt - crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,gr,ip,el) ! grain orientation as quaternion + crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,gr,ip,el) ! grain orientation as quaternion case ('eulerangles') mySize = 3_pInt - crystallite_postResults(c+1:c+mySize) = inDeg * math_qToEuler(crystallite_orientation(1:4,gr,ip,el)) ! grain orientation as Euler angles in degree + crystallite_postResults(c+1:c+mySize) = inDeg * & + math_qToEuler(crystallite_orientation(1:4,gr,ip,el)) ! grain orientation as Euler angles in degree case ('grainrotation') mySize = 4_pInt - crystallite_postResults(c+1:c+mySize) = math_qToAxisAngle(crystallite_rotation(1:4,gr,ip,el)) ! grain rotation away from initial orientation as axis-angle in crystal reference coordinates - crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree + crystallite_postResults(c+1:c+mySize) = math_qToAxisAngle(crystallite_rotation(1:4,gr,ip,el))! grain rotation away from initial orientation as axis-angle in crystal reference coordinates + crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree case ('grainrotationx') mySize = 1_pInt rotation = math_qToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,gr,ip,el), & - crystallite_rotation(1:4,gr,ip,el)), & - math_qConj(crystallite_orientation(1:4,gr,ip,el)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree + crystallite_rotation(1:4,gr,ip,el)), & + math_qConj(crystallite_orientation(1:4,gr,ip,el)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree case ('grainrotationy') mySize = 1_pInt rotation = math_qToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,gr,ip,el), & - crystallite_rotation(1:4,gr,ip,el)), & - math_qConj(crystallite_orientation(1:4,gr,ip,el)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree + crystallite_rotation(1:4,gr,ip,el)), & + math_qConj(crystallite_orientation(1:4,gr,ip,el)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree case ('grainrotationz') mySize = 1_pInt rotation = math_qToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,gr,ip,el), & - crystallite_rotation(1:4,gr,ip,el)), & - math_qConj(crystallite_orientation(1:4,gr,ip,el)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree + crystallite_rotation(1:4,gr,ip,el)), & + math_qConj(crystallite_orientation(1:4,gr,ip,el)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates + crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree case ('ipcoords') mySize = 3_pInt - crystallite_postResults(c+1:c+mySize) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates + crystallite_postResults(c+1:c+mySize) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates ! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33 ! thus row index i is slow, while column index j is fast. reminder: "row is slow" case ('defgrad','f') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_partionedF(1:3,1:3,gr,ip,el)),[mySize]) + crystallite_postResults(c+1:c+mySize) = & + reshape(math_transpose33(crystallite_partionedF(1:3,1:3,gr,ip,el)),[mySize]) case ('e') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( & @@ -3539,23 +3542,29 @@ function crystallite_postResults(dt, gr, ip, el) crystallite_partionedF(1:3,1:3,gr,ip,el)) - math_I3),[mySize]) case ('fe') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fe(1:3,1:3,gr,ip,el)),[mySize]) + crystallite_postResults(c+1:c+mySize) = & + reshape(math_transpose33(crystallite_Fe(1:3,1:3,gr,ip,el)),[mySize]) case ('ee') - Ee = 0.5_pReal * (math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,gr,ip,el)), crystallite_Fe(1:3,1:3,gr,ip,el)) - math_I3) + Ee = 0.5_pReal *(math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,gr,ip,el)), & + crystallite_Fe(1:3,1:3,gr,ip,el)) - math_I3) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize]) case ('fp') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fp(1:3,1:3,gr,ip,el)),[mySize]) + crystallite_postResults(c+1:c+mySize) = & + reshape(math_transpose33(crystallite_Fp(1:3,1:3,gr,ip,el)),[mySize]) case ('lp') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Lp(1:3,1:3,gr,ip,el)),[mySize]) + crystallite_postResults(c+1:c+mySize) = & + reshape(math_transpose33(crystallite_Lp(1:3,1:3,gr,ip,el)),[mySize]) case ('p','firstpiola','1stpiola') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_P(1:3,1:3,gr,ip,el)),[mySize]) + crystallite_postResults(c+1:c+mySize) = & + reshape(math_transpose33(crystallite_P(1:3,1:3,gr,ip,el)),[mySize]) case ('s','tstar','secondpiola','2ndpiola') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,gr,ip,el)),[mySize]) + crystallite_postResults(c+1:c+mySize) = & + reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,gr,ip,el)),[mySize]) case ('elasmatrix') mySize = 36_pInt crystallite_postResults(c+1:c+mySize) = reshape(constitutive_homogenizedC(gr,ip,el),[mySize]) @@ -3576,10 +3585,9 @@ function crystallite_postResults(dt, gr, ip, el) crystallite_postResults(c+1) = real(constitutive_sizePostResults(gr,ip,el),pReal) ! size of constitutive results c = c + 1_pInt if (constitutive_sizePostResults(gr,ip,el) > 0_pInt) & - crystallite_postResults(c+1:c+constitutive_sizePostResults(gr,ip,el)) = constitutive_postResults(crystallite_Tstar_v(1:6,gr,ip,el), & - crystallite_Fe, & - crystallite_Temperature(gr,ip,el), & - dt, gr, ip, el) + crystallite_postResults(c+1:c+constitutive_sizePostResults(gr,ip,el)) = & + constitutive_postResults(crystallite_Tstar_v(1:6,gr,ip,el), crystallite_Fe, & + crystallite_Temperature(gr,ip,el), dt, gr, ip, el) c = c + constitutive_sizePostResults(gr,ip,el) end function crystallite_postResults diff --git a/code/mesh.f90 b/code/mesh.f90 index b3f46a228..5c01db930 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -5121,9 +5121,9 @@ subroutine mesh_write_cellGeom filename = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'_ipbased.vtk', & mesh_topology = 'UNSTRUCTURED_GRID') err = VTK_geo(NN = mesh_Ncellnodes, & - X = mesh_cellnode(1,:), & - Y = mesh_cellnode(2,:), & - Z = mesh_cellnode(3,:)) + X = mesh_cellnode(1,1:mesh_Ncellnodes), & + Y = mesh_cellnode(2,1:mesh_Ncellnodes), & + Z = mesh_cellnode(3,1:mesh_Ncellnodes)) err = VTK_con(NC = mesh_Ncells, & connect = cellconnection(1:j), & cell_type = celltype) @@ -5161,18 +5161,18 @@ subroutine mesh_write_elemGeom i = i + 1_pInt + FE_Nnodes(t) enddo - err = VTK_ini(output_format = 'ASCII', & - title=trim(getSolverJobName())//' element mesh', & - filename = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName()//'_nodebased.vtk', & - mesh_topology = 'UNSTRUCTURED_GRID') - err = VTK_geo(NN = mesh_Nnodes, & - X = mesh_node0(1,1:mesh_Nnodes), & - Y = mesh_node0(2,1:mesh_Nnodes), & - Z = mesh_node0(3,1:mesh_Nnodes)) - err = VTK_con(NC = mesh_Nelems, & - connect = elementconnection(1:i), & - cell_type = elemtype) - err = VTK_end() + err =VTK_ini(output_format = 'ASCII', & + title=trim(getSolverJobName())//' element mesh', & + filename = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'_nodebased.vtk', & + mesh_topology = 'UNSTRUCTURED_GRID') + err =VTK_geo(NN = mesh_Nnodes, & + X = mesh_node0(1,1:mesh_Nnodes), & + Y = mesh_node0(2,1:mesh_Nnodes), & + Z = mesh_node0(3,1:mesh_Nnodes)) + err =VTK_con(NC = mesh_Nelems, & + connect = elementconnection(1:i), & + cell_type = elemtype) + err =VTK_end() end subroutine mesh_write_elemGeom diff --git a/lib/Lib_VTK_IO.f90 b/lib/Lib_VTK_IO.f90 index d38313223..85cd70970 100644 --- a/lib/Lib_VTK_IO.f90 +++ b/lib/Lib_VTK_IO.f90 @@ -3435,6 +3435,7 @@ contains !--------------------------------------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------------------------------------- + rf = 1 if (present(cf)) rf = cf close(unit=vtk(rf)%u,iostat=E_IO) call vtk_update(act='remove',cf=rf,Nvtk=Nvtk,vtk=vtk)