From 2ffb0eade5a842da6d423db7431a7a287fc8e4f3 Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Tue, 27 Mar 2018 17:24:32 +0200 Subject: [PATCH 01/11] comment 1 --- src/opencmiss_iron.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index 944490b7..1a9bc9b8 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -27676,6 +27676,8 @@ SUBROUTINE cmfe_Field_CreateFinishObj(field,err) CALL TAU_STATIC_PHASE_STOP('field Create') #endif +! my comment + EXITS("cmfe_Field_CreateFinishObj") RETURN 999 ERRORSEXITS("cmfe_Field_CreateFinishObj",err,error) From c83a531c7a9805b4c4b8c7b51884e8744544039b Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Tue, 27 Mar 2018 17:31:39 +0200 Subject: [PATCH 02/11] comment gone --- src/opencmiss_iron.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index 1a9bc9b8..944490b7 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -27676,8 +27676,6 @@ SUBROUTINE cmfe_Field_CreateFinishObj(field,err) CALL TAU_STATIC_PHASE_STOP('field Create') #endif -! my comment - EXITS("cmfe_Field_CreateFinishObj") RETURN 999 ERRORSEXITS("cmfe_Field_CreateFinishObj",err,error) From f156df320c77a2b2c61781d4785c0014a4f8828d Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Tue, 27 Mar 2018 17:32:46 +0200 Subject: [PATCH 03/11] Another comment... --- src/opencmiss_iron.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index 944490b7..c167d89d 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -27676,6 +27676,8 @@ SUBROUTINE cmfe_Field_CreateFinishObj(field,err) CALL TAU_STATIC_PHASE_STOP('field Create') #endif +! another comment + EXITS("cmfe_Field_CreateFinishObj") RETURN 999 ERRORSEXITS("cmfe_Field_CreateFinishObj",err,error) From 184db9644cab0f2d506173902fd6422bf792bc56 Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Mon, 6 Aug 2018 13:26:24 +0200 Subject: [PATCH 04/11] Simplex fix --- src/basis_routines.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/basis_routines.F90 b/src/basis_routines.F90 index 98f14871..46a95b76 100644 --- a/src/basis_routines.F90 +++ b/src/basis_routines.F90 @@ -6936,7 +6936,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*) CALL FlagError(localError,err,error,*999) ENDIF !Gauss point 1 - x(1,1)=1.0_DP/2.0 + x(1,1)=1.0_DP/2.0_DP x(2,1)=1.0_DP/2.0_DP w(1)=1.0_DP CASE(2) @@ -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 @@ -7017,7 +7017,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 @@ -7099,8 +7099,8 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*) CALL FlagError(localError,err,error,*999) ENDIF lC=1.0_DP/3.0_DP - wC=-9.0_DP/16.0 - alpha1=25.0_DP + wC=-3.0_DP/4.0 + alpha1=2.0_DP/5.0_DP wAlpha1=25.0_DP/48.0_DP l1Alpha1=(1.0_DP+2.0_DP*alpha1)/3.0_DP l2Alpha1=(1.0_DP-alpha1)/3.0_DP From 3bfcd8aefdb644f4d43261bf2869ffc4c3d93bda Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Mon, 6 Aug 2018 13:41:36 +0200 Subject: [PATCH 05/11] Simplex fix --- src/basis_routines.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/basis_routines.F90 b/src/basis_routines.F90 index 46a95b76..b26d709e 100644 --- a/src/basis_routines.F90 +++ b/src/basis_routines.F90 @@ -7099,7 +7099,7 @@ SUBROUTINE Gauss_Simplex(order,numberOfVertices,n,x,w,err,error,*) CALL FlagError(localError,err,error,*999) ENDIF lC=1.0_DP/3.0_DP - wC=-3.0_DP/4.0 + wC=-3.0_DP/4.0_DP alpha1=2.0_DP/5.0_DP wAlpha1=25.0_DP/48.0_DP l1Alpha1=(1.0_DP+2.0_DP*alpha1)/3.0_DP From 9afb4c24029901030453154a1fb46d8fe74bbef7 Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Thu, 18 Oct 2018 16:49:26 +0200 Subject: [PATCH 06/11] Comments in EQUATIONS_SET_BACKSUBSTITUTE cf. deludeln wrong in Laplace --- src/equations_set_routines.F90 | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/equations_set_routines.F90 b/src/equations_set_routines.F90 index c37c09c7..cef93296 100644 --- a/src/equations_set_routines.F90 +++ b/src/equations_set_routines.F90 @@ -1793,6 +1793,7 @@ SUBROUTINE EQUATIONS_SET_BACKSUBSTITUTE(equationsSet,BOUNDARY_CONDITIONS,err,err & ROW_INDICES,COLUMN_INDICES,err,error,*999) !Loop over the non-ghosted rows in the equations set DO equations_row_number=1,vectorMapping%numberOfRows + ! Same as other case (BLOCK) RHS_VALUE=0.0_DP rhs_variable_dof=rhsMapping%equationsRowToRHSDOFMap(equations_row_number) rhs_global_dof=RHS_DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(rhs_variable_dof) @@ -1801,6 +1802,7 @@ SUBROUTINE EQUATIONS_SET_BACKSUBSTITUTE(equationsSet,BOUNDARY_CONDITIONS,err,err CASE(BOUNDARY_CONDITION_DOF_FREE) !Back substitute !Loop over the local columns of the equations matrix + ! Different vs. BLOCK! DO equations_column_idx=ROW_INDICES(equations_row_number), & ROW_INDICES(equations_row_number+1)-1 equations_column_number=COLUMN_INDICES(equations_column_idx) @@ -1809,6 +1811,20 @@ SUBROUTINE EQUATIONS_SET_BACKSUBSTITUTE(equationsSet,BOUNDARY_CONDITIONS,err,err DEPENDENT_VALUE=DEPENDENT_PARAMETERS(variable_dof) RHS_VALUE=RHS_VALUE+MATRIX_VALUE*DEPENDENT_VALUE ENDDO !equations_column_idx + + ! CASE block storage above + !Back substitute + !Loop over the local columns of the equations matrix + !DO equations_column_idx=1,COLUMN_DOMAIN_MAPPING%TOTAL_NUMBER_OF_LOCAL + ! equations_column_number=COLUMN_DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP( & + ! & equations_column_idx) + ! variable_dof=equations_column_idx + ! MATRIX_VALUE=equationsMatrixData(equations_row_number+ & + ! & (equations_column_number-1)*vectorMatrices%totalNumberOfRows) + ! DEPENDENT_VALUE=DEPENDENT_PARAMETERS(variable_dof) + ! RHS_VALUE=RHS_VALUE+MATRIX_VALUE*DEPENDENT_VALUE + !ENDDO !equations_column_idx + CASE(BOUNDARY_CONDITION_DOF_FIXED) !Do nothing CASE(BOUNDARY_CONDITION_DOF_MIXED) From 6b12bdd61987d0949112c810a3c722abf5ef4cda Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Fri, 16 Nov 2018 13:31:08 +0100 Subject: [PATCH 07/11] Non-affecting changes towards solving Neumann_point incorrect behavior. --- src/boundary_condition_routines.F90 | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/boundary_condition_routines.F90 b/src/boundary_condition_routines.F90 index 2a9e9368..87412b08 100644 --- a/src/boundary_condition_routines.F90 +++ b/src/boundary_condition_routines.F90 @@ -2230,7 +2230,7 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab INTEGER(INTG) :: numberOfPointDofs, numberNonZeros, numberRowEntries, neumannConditionNumber, localNeumannConditionIdx INTEGER(INTG) :: neumannIdx, globalDof, localDof, localDofNyy, domainIdx, numberOfDomains, domainNumber, componentNumber INTEGER(INTG) :: nodeIdx, derivIdx, nodeNumber, versionNumber, derivativeNumber, columnNodeNumber, lineIdx, faceIdx, columnDof - INTEGER(INTG), ALLOCATABLE :: rowIndices(:), columnIndices(:), localDofNumbers(:) + INTEGER(INTG), ALLOCATABLE :: rowIndices(:), columnIndices(:), localDofNumbers(:), tempArray(:) REAL(DP) :: pointValue INTEGER(INTG) :: dummyErr TYPE(VARYING_STRING) :: dummyError @@ -2485,7 +2485,9 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab END DO !local DOFs CALL LIST_DESTROY(rowColumnIndicesList,err,error,*999) - CALL LIST_DETACH_AND_DESTROY(columnIndicesList,numberNonZeros,columnIndices,err,error,*999) + CALL LIST_DETACH_AND_DESTROY(columnIndicesList,numberNonZeros,tempArray,err,error,*999) + columnIndices=tempArray(1:numberNonZeros) + IF(ALLOCATED(tempArray)) DEALLOCATE(tempArray) IF(DIAGNOSTICS1) THEN CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Neumann integration matrix sparsity",err,error,*999) CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Number non-zeros = ", numberNonZeros,err,error,*999) @@ -2654,7 +2656,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,* INTEGER(INTG) :: faceNumber,lineNumber INTEGER(INTG) :: ms,os,nodeNumber,derivativeNumber,versionNumber LOGICAL :: dependentGeometry - REAL(DP) :: integratedValue,phim,phio + REAL(DP) :: integratedValue,phim,phio, integratedValueSum TYPE(BoundaryConditionsNeumannType), POINTER :: neumannConditions TYPE(BASIS_TYPE), POINTER :: basis TYPE(FIELD_TYPE), POINTER :: geometricField @@ -2679,6 +2681,8 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,* NULLIFY(interpolatedPointMetrics) NULLIFY(integratedValues) + integratedValueSum = 0.0_DP + neumannConditions=>rhsBoundaryConditions%neumannBoundaryConditions !Check that Neumann conditions are associated, otherwise do nothing IF(ASSOCIATED(neumannConditions)) THEN @@ -2913,9 +2917,16 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,* ! Add integral term to N matrix CALL DistributedMatrix_ValuesAdd(neumannConditions%integrationMatrix,localDof,neumannDofIdx, & & integratedValue,err,error,*999) + integratedValueSum = integratedValueSum+integratedValue END DO END DO END DO facesLoop + WRITE(*,*) "neumannGlobalDof" + WRITE(*,*) neumannGlobalDof + WRITE(*,*) "Value" + WRITE(*,*) integratedValueSum + integratedValueSum = 0.0_DP + CASE DEFAULT CALL FlagError("The dimension is invalid for point Neumann conditions",err,error,*999) END SELECT From 189395e1bc7bbb002acb23cbd9c5a5b423fccef4 Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Fri, 21 Dec 2018 14:00:00 +0100 Subject: [PATCH 08/11] Fixed 2D at FINITE_ELASTICITY_GAUSS_STRESS_TENSOR case Mooney-Rivlin. --- src/finite_elasticity_routines.F90 | 98 ++++++++++++++++++++++++------ 1 file changed, 78 insertions(+), 20 deletions(-) diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index ab1ebcf6..2f9a6030 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -1683,13 +1683,18 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN & GEOMETRIC_INTERPOLATED_POINT_METRICS,FIBRE_INTERPOLATED_POINT,DZDNU,ERR,ERROR,*999) Jznu=DEPENDENT_INTERPOLATED_POINT_METRICS%JACOBIAN/GEOMETRIC_INTERPOLATED_POINT_METRICS%JACOBIAN + WRITE (*,*) "Jznu" + WRITE (*,*) Jznu JGW=DEPENDENT_INTERPOLATED_POINT_METRICS%JACOBIAN*DEPENDENT_QUADRATURE_SCHEME%GAUSS_WEIGHTS(gauss_idx) !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) + WRITE (*,*) "Number of dimensions" + WRITE (*,*) numberOfDimensions + ! Convert from Voigt form to tensor form and multiply with Jacobian and Gauss weight. DO nh=1,numberOfDimensions DO mh=1,numberOfDimensions @@ -2038,7 +2043,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) @@ -3714,6 +3719,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) @@ -4244,8 +4250,11 @@ 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) + + WRITE (*,*) "Number of dimensions" + WRITE (*,*) numberOfDimensions + !Convert from Voigt form to tensor form. \TODO needs to be generalised for 2D DO nh=1,numberOfDimensions DO mh=1,numberOfDimensions @@ -4495,8 +4504,11 @@ 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) + WRITE (*,*) "Number of dimensions" + WRITE (*,*) numberOfDimensions + !Convert from Voigt form to tensor form. \TODO needs to be generalised for 2D DO nh=1,numberOfDimensions DO mh=1,numberOfDimensions @@ -4707,11 +4719,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 ! numberOfDimensions?? + DO mh=1,3 ! numberOfDimensions?? cauchyStressTensor(mh,nh)=cauchyStressVoigt(TENSOR_TO_VOIGT3(mh,nh)) ENDDO ENDDO @@ -6944,7 +6956,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) @@ -6990,15 +7006,19 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO !W=c1*(I1-3)+c2*(I2-3)+p/2*(I3-1) !Calculate isochoric fictitious 2nd Piola tensor (in Voigt form) - I1=AZL(1,1)+AZL(2,2)+AZL(3,3) + ! Needs a 2D test-case for Mooney-Rivlin!!! (i.e. c2=/0) + 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 @@ -7017,8 +7037,41 @@ 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 + 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 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)) + END DO + END DO + 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 + + ! 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 (*,*) "Cauchy Stress Tensor" + WRITE (*,*) STRESS_TENSOR CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) @@ -7053,6 +7106,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))// & From 9a1834dd90e58c4d89fed25eac9e84a71798dffe Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Fri, 21 Dec 2018 14:57:22 +0100 Subject: [PATCH 09/11] Some more refinements at finite_elasticity. --- src/finite_elasticity_routines.F90 | 47 +++++++++++++++++------------- 1 file changed, 26 insertions(+), 21 deletions(-) diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index 4ef4722b..35bf58ce 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -1709,8 +1709,6 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN & GEOMETRIC_INTERPOLATED_POINT_METRICS,FIBRE_INTERPOLATED_POINT,DZDNU,ERR,ERROR,*999) Jznu=DEPENDENT_INTERPOLATED_POINT_METRICS%JACOBIAN/GEOMETRIC_INTERPOLATED_POINT_METRICS%JACOBIAN - WRITE (*,*) "Jznu" - WRITE (*,*) Jznu JGW=DEPENDENT_INTERPOLATED_POINT_METRICS%JACOBIAN*DEPENDENT_QUADRATURE_SCHEME%GAUSS_WEIGHTS(gauss_idx) !Calculate the Cauchy stress tensor (in Voigt form) at the gauss point. @@ -1718,9 +1716,6 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,STRESS_TENSOR,DZDNU,Jznu, & & elementNumber,gauss_idx,numberOfDimensions,err,error,*999) - WRITE (*,*) "Number of dimensions" - WRITE (*,*) numberOfDimensions - ! Convert from Voigt form to tensor form and multiply with Jacobian and Gauss weight. DO nh=1,numberOfDimensions DO mh=1,numberOfDimensions @@ -4278,9 +4273,6 @@ SUBROUTINE FiniteElasticity_StressStrainPoint(equationsSet,evaluateType,numberOf CALL FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(equationsSet,dependentInterpolatedPoint, & & materialsInterpolatedPoint,geometricInterpolatedPoint,cauchyStressVoigt,dZdNu,Jznu, & & elementNumber,0,numberOfDimensions,ERR,ERROR,*999) - - WRITE (*,*) "Number of dimensions" - WRITE (*,*) numberOfDimensions !Convert from Voigt form to tensor form. \TODO needs to be generalised for 2D DO nh=1,numberOfDimensions @@ -4532,9 +4524,6 @@ SUBROUTINE FiniteElasticity_TensorInterpolateGaussPoint(equationsSet,tensorEvalu CALL FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(equationsSet,dependentInterpolatedPoint, & & materialsInterpolatedPoint,geometricInterpolatedPoint,cauchyStressVoigt,dZdNu,Jznu, & & localElementNumber,0,numberOfDimensions,ERR,ERROR,*999) - - WRITE (*,*) "Number of dimensions" - WRITE (*,*) numberOfDimensions !Convert from Voigt form to tensor form. \TODO needs to be generalised for 2D DO nh=1,numberOfDimensions @@ -7018,7 +7007,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO REAL(DP) :: ONETHIRD_TRACE TYPE(VARYING_STRING) :: LOCAL_ERROR TYPE(FIELD_VARIABLE_TYPE), POINTER :: FIELD_VARIABLE - REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3), DZDNUInv(3,3), PK2Tensor(3,3) + REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3), DZDNUInv(3,3), PK2Tensor(3,3), sigmaTensor(3,3) REAL(DP) :: B(6),E(6),DQ_DE(6), PK2(6), Jznu_dummy REAL(DP), POINTER :: C(:) !Parameters for constitutive laws @@ -7045,8 +7034,13 @@ 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) - ! Needs a 2D test-case for Mooney-Rivlin!!! (i.e. c2=/0) I1 = 0.0_DP DO component_idx=1,numberOfDimensions I1=I1+AZL(component_idx,component_idx) @@ -7080,20 +7074,19 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:numberOfDimensions))/numberOfDimensions_DP STRESS_TENSOR(1:numberOfDimensions)=STRESS_TENSOR(1:numberOfDimensions)-ONETHIRD_TRACE+P - ! Compute PK2 + ! 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 from Voigt to tensor + + ! 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 - 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 ! Write out PK2 IF(DIAGNOSTICS1) THEN @@ -7109,9 +7102,21 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO & '(16X,2(X,E13.6))',err,error,*999) END IF END IF - - WRITE (*,*) "Cauchy Stress Tensor" - WRITE (*,*) STRESS_TENSOR + + ! 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) From be11cd2efd59928144c9ee234d168b475c4df7eb Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Fri, 21 Dec 2018 15:13:25 +0100 Subject: [PATCH 10/11] Pull request just finite elasticity (and one small at integration) --- src/boundary_condition_routines.F90 | 18 ++++-------------- src/equations_set_routines.F90 | 17 +---------------- src/opencmiss_iron.F90 | 2 -- 3 files changed, 5 insertions(+), 32 deletions(-) diff --git a/src/boundary_condition_routines.F90 b/src/boundary_condition_routines.F90 index 36f1eb01..8e23b509 100644 --- a/src/boundary_condition_routines.F90 +++ b/src/boundary_condition_routines.F90 @@ -2242,7 +2242,7 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab INTEGER(INTG) :: numberOfPointDofs, numberNonZeros, numberRowEntries, neumannConditionNumber, localNeumannConditionIdx INTEGER(INTG) :: neumannIdx, globalDof, localDof, localDofNyy, domainIdx, numberOfDomains, domainNumber, componentNumber INTEGER(INTG) :: nodeIdx, derivIdx, nodeNumber, versionNumber, derivativeNumber, columnNodeNumber, lineIdx, faceIdx, columnDof - INTEGER(INTG), ALLOCATABLE :: rowIndices(:), columnIndices(:), localDofNumbers(:), tempArray(:) + INTEGER(INTG), ALLOCATABLE :: rowIndices(:), columnIndices(:), localDofNumbers(:) REAL(DP) :: pointValue INTEGER(INTG) :: dummyErr TYPE(VARYING_STRING) :: dummyError @@ -2497,9 +2497,7 @@ SUBROUTINE BoundaryConditions_NeumannMatricesInitialise(boundaryConditionsVariab END DO !local DOFs CALL LIST_DESTROY(rowColumnIndicesList,err,error,*999) - CALL LIST_DETACH_AND_DESTROY(columnIndicesList,numberNonZeros,tempArray,err,error,*999) - columnIndices=tempArray(1:numberNonZeros) - IF(ALLOCATED(tempArray)) DEALLOCATE(tempArray) + CALL LIST_DETACH_AND_DESTROY(columnIndicesList,numberNonZeros,columnIndices,err,error,*999) IF(DIAGNOSTICS1) THEN CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Neumann integration matrix sparsity",err,error,*999) CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Number non-zeros = ", numberNonZeros,err,error,*999) @@ -2668,7 +2666,7 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,* INTEGER(INTG) :: faceNumber,lineNumber INTEGER(INTG) :: ms,os,nodeNumber,derivativeNumber,versionNumber LOGICAL :: dependentGeometry - REAL(DP) :: integratedValue,phim,phio, integratedValueSum + REAL(DP) :: integratedValue,phim,phio TYPE(BoundaryConditionsNeumannType), POINTER :: neumannConditions TYPE(BASIS_TYPE), POINTER :: basis TYPE(FIELD_TYPE), POINTER :: geometricField @@ -2693,8 +2691,6 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,* NULLIFY(interpolatedPointMetrics) NULLIFY(integratedValues) - integratedValueSum = 0.0_DP - neumannConditions=>rhsBoundaryConditions%neumannBoundaryConditions !Check that Neumann conditions are associated, otherwise do nothing IF(ASSOCIATED(neumannConditions)) THEN @@ -2929,16 +2925,9 @@ SUBROUTINE BoundaryConditions_NeumannIntegrate(rhsBoundaryConditions,err,error,* ! Add integral term to N matrix CALL DistributedMatrix_ValuesAdd(neumannConditions%integrationMatrix,localDof,neumannDofIdx, & & integratedValue,err,error,*999) - integratedValueSum = integratedValueSum+integratedValue END DO END DO END DO facesLoop - WRITE(*,*) "neumannGlobalDof" - WRITE(*,*) neumannGlobalDof - WRITE(*,*) "Value" - WRITE(*,*) integratedValueSum - integratedValueSum = 0.0_DP - CASE DEFAULT CALL FlagError("The dimension is invalid for point Neumann conditions",err,error,*999) END SELECT @@ -3988,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 5735f56d..37d385f0 100644 --- a/src/equations_set_routines.F90 +++ b/src/equations_set_routines.F90 @@ -1792,7 +1792,6 @@ SUBROUTINE EQUATIONS_SET_BACKSUBSTITUTE(equationsSet,BOUNDARY_CONDITIONS,err,err & ROW_INDICES,COLUMN_INDICES,err,error,*999) !Loop over the non-ghosted rows in the equations set DO equations_row_number=1,vectorMapping%numberOfRows - ! Same as other case (BLOCK) RHS_VALUE=0.0_DP rhs_variable_dof=rhsMapping%equationsRowToRHSDOFMap(equations_row_number) rhs_global_dof=RHS_DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(rhs_variable_dof) @@ -1801,7 +1800,6 @@ SUBROUTINE EQUATIONS_SET_BACKSUBSTITUTE(equationsSet,BOUNDARY_CONDITIONS,err,err CASE(BOUNDARY_CONDITION_DOF_FREE) !Back substitute !Loop over the local columns of the equations matrix - ! Different vs. BLOCK! DO equations_column_idx=ROW_INDICES(equations_row_number), & ROW_INDICES(equations_row_number+1)-1 equations_column_number=COLUMN_INDICES(equations_column_idx) @@ -1810,20 +1808,6 @@ SUBROUTINE EQUATIONS_SET_BACKSUBSTITUTE(equationsSet,BOUNDARY_CONDITIONS,err,err DEPENDENT_VALUE=DEPENDENT_PARAMETERS(variable_dof) RHS_VALUE=RHS_VALUE+MATRIX_VALUE*DEPENDENT_VALUE ENDDO !equations_column_idx - - ! CASE block storage above - !Back substitute - !Loop over the local columns of the equations matrix - !DO equations_column_idx=1,COLUMN_DOMAIN_MAPPING%TOTAL_NUMBER_OF_LOCAL - ! equations_column_number=COLUMN_DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP( & - ! & equations_column_idx) - ! variable_dof=equations_column_idx - ! MATRIX_VALUE=equationsMatrixData(equations_row_number+ & - ! & (equations_column_number-1)*vectorMatrices%totalNumberOfRows) - ! DEPENDENT_VALUE=DEPENDENT_PARAMETERS(variable_dof) - ! RHS_VALUE=RHS_VALUE+MATRIX_VALUE*DEPENDENT_VALUE - !ENDDO !equations_column_idx - CASE(BOUNDARY_CONDITION_DOF_FIXED) !Do nothing CASE(BOUNDARY_CONDITION_DOF_MIXED) @@ -7231,3 +7215,4 @@ END SUBROUTINE EquationsSet_ResidualEvaluateStaticNodal ! END MODULE EQUATIONS_SET_ROUTINES + diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90 index 2cf69121..300cf38f 100644 --- a/src/opencmiss_iron.F90 +++ b/src/opencmiss_iron.F90 @@ -29397,8 +29397,6 @@ SUBROUTINE cmfe_Field_CreateFinishObj(field,err) CALL TAU_STATIC_PHASE_STOP('field Create') #endif -! another comment - EXITS("cmfe_Field_CreateFinishObj") RETURN 999 ERRORSEXITS("cmfe_Field_CreateFinishObj",err,error) From 040db5afa90ac0a07e6e421fea9a8600bf1213dc Mon Sep 17 00:00:00 2001 From: lorenzo-mechbau Date: Fri, 21 Dec 2018 15:25:11 +0100 Subject: [PATCH 11/11] Small adjustment finite elasticity --- src/finite_elasticity_routines.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index 35bf58ce..d89747d3 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -4738,8 +4738,8 @@ SUBROUTINE FiniteElasticity_TensorInterpolateXi(equationsSet,tensorEvaluateType, & localElementNumber,0,numberOfDimensions,ERR,ERROR,*999) !Convert from Voigt form to tensor form. - DO nh=1,3 ! numberOfDimensions?? - DO mh=1,3 ! numberOfDimensions?? + DO nh=1,3 + DO mh=1,3 cauchyStressTensor(mh,nh)=cauchyStressVoigt(TENSOR_TO_VOIGT3(mh,nh)) ENDDO ENDDO