diff --git a/src/basis_routines.F90 b/src/basis_routines.F90 index 8700eb64..58f8f46b 100644 --- a/src/basis_routines.F90 +++ b/src/basis_routines.F90 @@ -6993,7 +6993,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*) ENDIF !Gauss point 1 x(1,1)=(1.0_DP-SQRT(0.6_DP))/2.0_DP - x(2,1)=1-x(1,1) + x(2,1)=1.0_DP-x(1,1) w(1)=5.0_DP/18.0_DP !Gauss point 2 x(1,2)=1.0_DP/2.0_DP diff --git a/src/boundary_condition_routines.F90 b/src/boundary_condition_routines.F90 index 8defcb08..8e23b509 100644 --- a/src/boundary_condition_routines.F90 +++ b/src/boundary_condition_routines.F90 @@ -3977,3 +3977,4 @@ END SUBROUTINE BOUNDARY_CONDITIONS_PRESSURE_INCREMENTED_INITIALISE ! END MODULE BOUNDARY_CONDITIONS_ROUTINES + diff --git a/src/equations_set_routines.F90 b/src/equations_set_routines.F90 index 10437fa8..37d385f0 100644 --- a/src/equations_set_routines.F90 +++ b/src/equations_set_routines.F90 @@ -7215,3 +7215,4 @@ END SUBROUTINE EquationsSet_ResidualEvaluateStaticNodal ! END MODULE EQUATIONS_SET_ROUTINES + diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index c5b8c450..d89747d3 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -1714,7 +1714,7 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Calculate the Cauchy stress tensor (in Voigt form) at the gauss point. CALL FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPOLATED_POINT, & & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,STRESS_TENSOR,DZDNU,Jznu, & - & elementNumber,gauss_idx,err,error,*999) + & elementNumber,gauss_idx,numberOfDimensions,err,error,*999) ! Convert from Voigt form to tensor form and multiply with Jacobian and Gauss weight. DO nh=1,numberOfDimensions @@ -2065,7 +2065,7 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Calculate the Cauchy stress tensor (in Voigt form) at the gauss point. CALL FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPOLATED_POINT, & & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,STRESS_TENSOR,Fe,Jznu, & - & elementNumber,gauss_idx,err,error,*999) + & elementNumber,gauss_idx,numberOfDimensions,err,error,*999) ! Convert from Voigt form to tensor form and multiply with Jacobian and Gauss weight. JGW=Jzxi*DEPENDENT_QUADRATURE_SCHEME%GAUSS_WEIGHTS(gauss_idx) @@ -3741,6 +3741,7 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field NULLIFY(coordinateSystem) CALL EquationsSet_CoordinateSystemGet(equationsSet,coordinateSystem,err,error,*999) numberOfDimensions=coordinateSystem%NUMBER_OF_DIMENSIONS + !Check the provided strain field variable has appropriate components and interpolation SELECT CASE(numberOfDimensions) CASE(3) @@ -4271,8 +4272,8 @@ SUBROUTINE FiniteElasticity_StressStrainPoint(equationsSet,evaluateType,numberOf ! of stress at any xi position and the GaussPoint number argument needs to be replace with a set of xi coordinates. CALL FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(equationsSet,dependentInterpolatedPoint, & & materialsInterpolatedPoint,geometricInterpolatedPoint,cauchyStressVoigt,dZdNu,Jznu, & - & elementNumber,0,ERR,ERROR,*999) - + & elementNumber,0,numberOfDimensions,ERR,ERROR,*999) + !Convert from Voigt form to tensor form. \TODO needs to be generalised for 2D DO nh=1,numberOfDimensions DO mh=1,numberOfDimensions @@ -4522,8 +4523,8 @@ SUBROUTINE FiniteElasticity_TensorInterpolateGaussPoint(equationsSet,tensorEvalu ! of stress at any xi position and the GaussPoint number argument needs to be replace with a set of xi coordinates. CALL FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(equationsSet,dependentInterpolatedPoint, & & materialsInterpolatedPoint,geometricInterpolatedPoint,cauchyStressVoigt,dZdNu,Jznu, & - & localElementNumber,0,ERR,ERROR,*999) - + & localElementNumber,0,numberOfDimensions,ERR,ERROR,*999) + !Convert from Voigt form to tensor form. \TODO needs to be generalised for 2D DO nh=1,numberOfDimensions DO mh=1,numberOfDimensions @@ -4734,11 +4735,11 @@ SUBROUTINE FiniteElasticity_TensorInterpolateXi(equationsSet,tensorEvaluateType, ! of stress at any xi position and the GaussPoint number argument needs to be replace with a set of xi coordinates. CALL FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(equationsSet,dependentInterpolatedPoint, & & materialsInterpolatedPoint,geometricInterpolatedPoint,cauchyStressVoigt,dZdNu,Jznu, & - & localElementNumber,0,ERR,ERROR,*999) + & localElementNumber,0,numberOfDimensions,ERR,ERROR,*999) !Convert from Voigt form to tensor form. - DO nh=1,3 - DO mh=1,3 + DO nh=1,3 + DO mh=1,3 cauchyStressTensor(mh,nh)=cauchyStressVoigt(TENSOR_TO_VOIGT3(mh,nh)) ENDDO ENDDO @@ -6984,7 +6985,7 @@ END SUBROUTINE FiniteElasticity_StrainTensor !>Evaluates the Cauchy stress tensor at a given Gauss point SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPOLATED_POINT, & & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,STRESS_TENSOR,DZDNU,Jznu, & - & ELEMENT_NUMBER,GAUSS_POINT_NUMBER,err,error,*) + & ELEMENT_NUMBER,GAUSS_POINT_NUMBER,numberOfDimensions,err,error,*) !Argument variables TYPE(EQUATIONS_SET_TYPE), POINTER, INTENT(IN) :: EQUATIONS_SET !MATERIALS_INTERPOLATED_POINT%VALUES(:,NO_PART_DERIV) @@ -7029,16 +7034,25 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO !Form of constitutive model is: !W=c1*(I1-3)+c2*(I2-3)+p/2*(I3-1) + ! This function needs a 2D test-case for Mooney-Rivlin!!! (i.e. c2=/0) + ! Therefore, add a warning. + IF ((CEILING(C(2))/=0) .AND. numberOfDimensions==2) THEN + CALL FlagWarning("Mooney-Rivlin material (c1,c2\=0) not verified for 2D-case.",err,error,*999) + END IF + !Calculate isochoric fictitious 2nd Piola tensor (in Voigt form) - I1=AZL(1,1)+AZL(2,2)+AZL(3,3) + I1 = 0.0_DP + DO component_idx=1,numberOfDimensions + I1=I1+AZL(component_idx,component_idx) + END DO TEMPTERM1=-2.0_DP*C(2) TEMPTERM2=2.0_DP*(C(1)+I1*C(2)) STRESS_TENSOR(1)=TEMPTERM1*AZL(1,1)+TEMPTERM2 STRESS_TENSOR(2)=TEMPTERM1*AZL(2,2)+TEMPTERM2 - STRESS_TENSOR(3)=TEMPTERM1*AZL(3,3)+TEMPTERM2 - STRESS_TENSOR(4)=TEMPTERM1*AZL(2,1) - STRESS_TENSOR(5)=TEMPTERM1*AZL(3,1) - STRESS_TENSOR(6)=TEMPTERM1*AZL(3,2) + STRESS_TENSOR(3)=TEMPTERM1*AZL(3,3)+TEMPTERM2 ! meaningless if 2D + STRESS_TENSOR(4)=TEMPTERM1*AZL(2,1) + STRESS_TENSOR(5)=TEMPTERM1*AZL(3,1) ! meaningless if 2D + STRESS_TENSOR(6)=TEMPTERM1*AZL(3,2) ! meaningless if 2D IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE) THEN !add active contraction stress values @@ -7057,8 +7071,52 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO !Do push-forward of 2nd Piola tensor. CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,err,error,*999) !Calculate isochoric Cauchy tensor (the deviatoric part) and add the volumetric part (the hydrostatic pressure). - ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP - STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P + ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:numberOfDimensions))/numberOfDimensions_DP + STRESS_TENSOR(1:numberOfDimensions)=STRESS_TENSOR(1:numberOfDimensions)-ONETHIRD_TRACE+P + + ! Compute PK2 (just for comparison) + PK2 = STRESS_TENSOR + Jznu_dummy = Jznu + CALL INVERT(DZDNU,DZDNUInv,Jznu_dummy,err,error,*999) + CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(PK2,DZDNUInv,1.0_DP/Jznu,err,error,*999) + + ! PK2 and sigma from Voigt to tensor + DO component_idx =1,numberOfDimensions + DO dof_idx =1, numberOfDimensions + PK2Tensor(component_idx, dof_idx) = PK2(TENSOR_TO_VOIGT(component_idx,dof_idx, 3)) + sigmaTensor(component_idx, dof_idx) = STRESS_TENSOR(TENSOR_TO_VOIGT(component_idx,dof_idx, 3)) + END DO + END DO + + ! Write out PK2 + IF(DIAGNOSTICS1) THEN + CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"",err,error,*999) + CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Second PK Tensor (S):",err,error,*999) + IF (numberOfDimensions == 3) THEN + CALL WriteStringMatrix(DIAGNOSTIC_OUTPUT_TYPE,1,1,3,1,1,3, & + & 3,3,PK2Tensor,WRITE_STRING_MATRIX_NAME_AND_INDICES,'(" S','(",I1,",:)',' :",3(X,E13.6))', & + & '(16X,3(X,E13.6))',err,error,*999) + ELSE ! must be 2 + CALL WriteStringMatrix(DIAGNOSTIC_OUTPUT_TYPE,1,1,2,1,1,2, & + & 3,3,PK2Tensor(1:2,1:2),WRITE_STRING_MATRIX_NAME_AND_INDICES,'(" S','(",I1,",:)',' :",2(X,E13.6))', & + & '(16X,2(X,E13.6))',err,error,*999) + END IF + END IF + + ! Write out Cauchy Stress Tensor + IF(DIAGNOSTICS1) THEN + CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"",err,error,*999) + CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Cauchy Stress Tensor (\sigma):",err,error,*999) + IF (numberOfDimensions == 3) THEN + CALL WriteStringMatrix(DIAGNOSTIC_OUTPUT_TYPE,1,1,3,1,1,3, & + & 3,3,sigmaTensor,WRITE_STRING_MATRIX_NAME_AND_INDICES, '(" \sigma','(",I1,",:)',' :",3(X,E13.6))', & + & '(16X,3(X,E13.6))',err,error,*999) + ELSE ! must be 2 + CALL WriteStringMatrix(DIAGNOSTIC_OUTPUT_TYPE,1,1,2,1,1,2, & + & 3,3,sigmaTensor(1:2,1:2),WRITE_STRING_MATRIX_NAME_AND_INDICES,'(" \sigma','(",I1,",:)',' :",2(X,E13.6))', & + & '(16X,2(X,E13.6))',err,error,*999) + END IF + END IF CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) @@ -7093,6 +7151,11 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO !Calculate isochoric Cauchy tensor (the deviatoric part) and add the volumetric part (the hydrostatic pressure). ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P + + IF (numberOfDimensions==2) THEN + CALL FlagWarning("This case has not been verified for the 2D-case!!!",err,error,*999) + END IF + CASE DEFAULT LOCAL_ERROR="The third equations set specification of "// & & TRIM(NumberToVString(EQUATIONS_SET%specification(3),"*",err,error))// &