Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reaction Diffusion on Deforming Domain #189

Open
wants to merge 10 commits into
base: develop
Choose a base branch
from
Open

Reaction Diffusion on Deforming Domain #189

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
42ba19b
Merge pull request #1 from OpenCMISS/develop
jaredmcollette Jul 31, 2019
0b4b2c7
Added fitting to Nearly-Incompressible Mooney Rivlin
Aug 2, 2019
0ce479b
Merge pull request #2 from OpenCMISS/develop
jaredmcollette Mar 1, 2020
1b04bdb
Implementation of Linear Elasticity with Stokes Damping
jaredmcollette Mar 2, 2020
0710994
Implementation of Reaction Diffusion on Moving Domain
jaredmcollette Mar 2, 2020
2dd5048
Merge remote-tracking branch 'fork_jared/develop' into github_jared
jaredmcollette Mar 2, 2020
55ed651
Merge pull request #3 from jaredmcollette/github_jared
jaredmcollette Mar 2, 2020
9534342
Code Clean-up and moving Beta to material field
jaredmcollette Mar 2, 2020
d2bedf0
Merging Jareds Changes
jaredmcollette Mar 2, 2020
f382963
Merge pull request #4 from jaredmcollette/github_jared
jaredmcollette Mar 2, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 1 addition & 4 deletions src/coordinate_routines.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3885,10 +3885,7 @@ SUBROUTINE Coordinates_MaterialSystemCalculate(geometricInterpPointMetrics,fibre
ELSE
!No fibre field
numberOfNuDimensions=0
DO xiIdx=1,numberOfXiDimensions
dNudXiTemp(1:numberOfXDimensions,xiIdx)=geometricInterpPointMetrics%DX_DXI(1:numberOfXDimensions,xiIdx)
CALL Normalise(dNudXiTemp(1:numberOfXDimensions,xiIdx),dXdNu(1:numberOfXDimensions,xiIdx),err,error,*999)
ENDDO !xiIdx
CALL IdentityMatrix(dXdNU,err,error,*999)
ENDIF
!Calculate dNu/dX the inverse of dX/dNu (same as transpose due to orthogonality)
CALL MatrixTranspose(dXdNu(1:numberOfXDimensions,1:numberOfXDimensions),dNudX(1:numberOfXDimensions,1: &
Expand Down
3 changes: 2 additions & 1 deletion src/equations_routines.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -660,7 +660,8 @@ SUBROUTINE Equations_InterpolationInitialise(equations,err,error,*)
CALL Field_InterpolatedPointsInitialise(equations%interpolation%independentInterpParameters, &
& equations%interpolation%independentInterpPoint,err,error,*999)
IF(equations%interpolation%independentField%type==FIELD_GEOMETRIC_TYPE.OR. &
& equations%interpolation%independentField%type==FIELD_FIBRE_TYPE) THEN
& equations%interpolation%independentField%type==FIELD_FIBRE_TYPE.OR. &
& equations%interpolation%independentField%type==FIELD_GEOMETRIC_GENERAL_TYPE) THEN
CALL Field_InterpolatedPointsMetricsInitialise(equations%interpolation%independentInterpPoint, &
& equations%interpolation%independentInterpPointMetrics,err,error,*999)
END IF
Expand Down
2 changes: 2 additions & 0 deletions src/equations_set_constants.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,8 @@ MODULE EquationsSetConstants
INTEGER(INTG), PARAMETER :: EQUATIONS_SET_THREE_DIMENSIONAL_SUBTYPE=4
INTEGER(INTG), PARAMETER :: EQUATIONS_SET_PLATE_SUBTYPE=5
INTEGER(INTG), PARAMETER :: EQUATIONS_SET_SHELL_SUBTYPE=6
INTEGER(INTG), PARAMETER :: EQUATIONS_SET_ONE_DIM_STOKES_DAMPING_SUBTYPE=7
INTEGER(INTG), PARAMETER :: EQUATIONS_SET_TWO_DIM_PLANE_STRESS_STOKES_DAMPING_SUBTYPE=8
! Finite elasticity
INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE=1
INTEGER(INTG), PARAMETER :: EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE=2
Expand Down
44 changes: 40 additions & 4 deletions src/finite_elasticity_routines.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3757,7 +3757,7 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field
CALL FieldVariable_FieldGet(fieldVariable,field,err,error,*999)
fieldVarType=fieldVariable%VARIABLE_TYPE

CALL Field_NumberOfComponentsCheck(field,fieldVarType,6,err,error,*999)
CALL Field_NumberOfComponentsCheck(field,fieldVarType,numberOfComponents,err,error,*999)
CALL Field_ComponentInterpolationGet(field,fieldVarType,1,fieldInterpolation,err,error,*999)
!Check the interpolation type
SELECT CASE(fieldInterpolation)
Expand Down Expand Up @@ -4260,7 +4260,8 @@ SUBROUTINE FiniteElasticity_StressStrainPoint(equationsSet,evaluateType,numberOf

SELECT CASE(equationsSet%specification(3))
CASE(EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
& EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE)
& EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, &
& EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE)
!Calculate the Cauchy stress tensor (in Voigt form) at the gauss point.
Jznu=dependentInterpolatedPointMetrics%JACOBIAN/geometricInterpolatedPointMetrics%JACOBIAN
! Note that some problems, e.g. active contraction, require additonal fields to be evaluated at Gauss points. This is
Expand Down Expand Up @@ -6995,17 +6996,20 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
REAL(DP), INTENT(IN) :: Jznu !Determinant of deformation gradient tensor (AZL)
INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER,GAUSS_POINT_NUMBER !<Element/Gauss point number
INTEGER(INTG), INTENT(OUT) :: err !<The error code
INTEGER(INTG) :: i
TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string
!Local Variables
INTEGER(INTG) :: PRESSURE_COMPONENT,component_idx,dof_idx
REAL(DP) :: P
REAL(DP) :: I1 !Invariants, if needed
REAL(DP) :: I1,I2,I3 !Invariants, if needed
REAL(DP) :: TEMPTERM1,TEMPTERM2,VALUE !Temporary variables
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)
REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3),AZU(3,3),PIOLA_TENSOR(3,3),IDENTITY(3,3)
REAL(DP) :: AZL_SQUARED(3,3)
REAL(DP) :: B(6),E(6),DQ_DE(6)
REAL(DP) :: WV_PRIME
REAL(DP), POINTER :: C(:) !Parameters for constitutive laws

ENTERS("FINITE_ELASTICITY_GAUSS_STRESS_TENSOR",err,error,*999)
Expand All @@ -7021,6 +7025,38 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
C=>MATERIALS_INTERPOLATED_POINT%VALUES(:,NO_PART_DERIV)

SELECT CASE(EQUATIONS_SET%specification(3))
CASE(EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE)

CALL INVERT(AZL,AZU,I3,ERR,ERROR,*999)

IDENTITY=0.0_DP
DO i=1,3
IDENTITY(i,i)=1.0_DP
ENDDO

!Form of constitutive model is:
! W_hat=c1*(I1_hat-3)+c2*(I2_hat-3)+p*J*C^(-1) + W^v(J)
! take W^v(J) = 1/2 * kappa * (J-1)^2
WV_PRIME = C(3)*(Jznu - 1.0_DP)
!compute the invariants, I3 a few lines up
I1 = AZL(1,1) + AZL(2,2) + AZL(3,3)
CALL MatrixProduct(AZL,AZL,AZL_SQUARED,err,error,*999)
I2 = 0.5_DP * (I1**2 - AZL_SQUARED(1,1) - AZL_SQUARED(2,2) - AZL_SQUARED(3,3))

PIOLA_TENSOR=2.0_DP*Jznu**(-2.0_DP/3.0_DP)*((C(1)+C(2)*I1)*IDENTITY-C(2)*AZL &
& -(C(1)*I1+2.0_DP*C(2)*I2-1.5_DP*WV_PRIME*Jznu**(5.0_DP/3.0_DP))/3.0_DP*AZU)

!Calculate isochoric fictitious 2nd Piola tensor (in Voigt form)
STRESS_TENSOR(1)=PIOLA_TENSOR(1,1)
STRESS_TENSOR(2)=PIOLA_TENSOR(2,2)
STRESS_TENSOR(3)=PIOLA_TENSOR(3,3)
STRESS_TENSOR(4)=PIOLA_TENSOR(2,1)
STRESS_TENSOR(5)=PIOLA_TENSOR(3,1)
STRESS_TENSOR(6)=PIOLA_TENSOR(3,2)

!Do push-forward of 2nd Piola tensor.
CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,err,error,*999)

CASE(EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
& EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
& EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE)
Expand Down
Loading