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.

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

Implement the rotation of the stiffness tensor from the local frame to the global frame (+ other fixes) #3

Open
wants to merge 17 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
72 changes: 72 additions & 0 deletions MEF90/m_MEF90_LinAlg.F90
Original file line number Diff line number Diff line change
Expand Up @@ -2843,6 +2843,78 @@ Function Tens4OS3DTransform(T,M)
Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
End Function Tens4OS3DTransform

Function Tens4OS2DXTens4OS2D(T1,T2)
!!! Compute the double contraction of 4th order tensors
!!! i.e. C_{ijkl} = T1_{ijpq}:T2_{pqkl}
!!! (c) 2022 Jean-Michel Scherer [email protected]
Type(Tens4OS2D),Intent(IN) :: T1,T2

Type(Tens4OS2D) :: Tens4OS2DXTens4OS2D

PetscReal,Dimension(2,2,2,2) :: TT1,TT2,C
Integer :: i,j,k,l
Integer :: p,q
PetscLogDouble :: flops
PetscErrorCode :: ierr

TT1 = T1
TT2 = T2
C = 0.0_Kr
Do i = 1,2
Do j = 1,2
Do k = 1,2
Do l = 1,2
Do p = 1,2
Do q = 1,2
C(i,j,k,l) = C(i,j,k,l) + TT1(i,j,p,q) * TT2(p,q,k,l)
End Do
End Do
End Do
End Do
End Do
End Do

Tens4OS2DXTens4OS2D = C
flops = 192 ! 2**6 * 3
Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
End Function Tens4OS2DXTens4OS2D

Function Tens4OS3DXTens4OS3D(T1,T2)
!!! Compute the double contraction of 4th order tensors
!!! i.e. C_{ijkl} = T1_{ijpq}:T2_{pqkl}
!!! (c) 2022 Jean-Michel Scherer [email protected]
Type(Tens4OS3D),Intent(IN) :: T1,T2

Type(Tens4OS3D) :: Tens4OS3DXTens4OS3D

PetscReal,Dimension(3,3,3,3) :: TT1,TT2,C
Integer :: i,j,k,l
Integer :: p,q
PetscLogDouble :: flops
PetscErrorCode :: ierr

TT1 = T1
TT2 = T2
C = 0.0_Kr
Do i = 1,3
Do j = 1,3
Do k = 1,3
Do l = 1,3
Do p = 1,3
Do q = 1,3
C(i,j,k,l) = C(i,j,k,l) + TT1(i,j,p,q) * TT2(p,q,k,l)
End Do
End Do
End Do
End Do
End Do
End Do

Tens4OS3DXTens4OS3D = C
flops = 2187 ! 3**6 * 3
Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
End Function Tens4OS3DXTens4OS3D

Function Tens4OS2DSquareRoot(T)
Type(Tens4OS2D),Intent(IN) :: T
Type(Tens4OS2D) :: Tens4OS2DSquareRoot
Expand Down
4 changes: 2 additions & 2 deletions Makefile.include
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@ endif
MEF90_INCLUDE=${MEF90_DIR}/objs/${PETSC_ARCH}

ifneq (${SNLP_DIR},)
FC_FLAGS+= -DMEF90_HAVE_SNLP -I${SNLP_DIR}/include
C_FLAGS+= -I${SNLP_DIR}/include
FC_FLAGS+= -DMEF90_HAVE_SNLP -I${SNLP_DIR}/include -I${SNLP_DIR}/obj
C_FLAGS+= -I${SNLP_DIR}/include -I${SNLP_DIR}/obj
SNLP_LIB:= -L${SNLP_DIR}/lib -lsnlp -lsnlpF90 -Wl,-rpath,${SNLP_DIR}/lib
endif

Expand Down
1 change: 1 addition & 0 deletions m_DefMech/m_MEF90_DefMechPlasticity.F90
Original file line number Diff line number Diff line change
Expand Up @@ -304,6 +304,7 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast
elemDisplacement,elemDisplacementType,elemScal,elemScalType,ierr)
Allocate(damageloc(elemScalType%numDof))
Do cell = 1,size(cellID)
!print *,"cell = ", cell
!! actualiser le ctx ( HookesLaw ,InelasticStrainSec, plasticStrainStrainSec, plasticStrainOldSec )
Call SectionRealRestrict(plasticStrainSec,cellID(cell),plasticStrainLoc,ierr);CHKERRQ(ierr)
Call SectionRealRestrict(plasticStrainOldSec,cellID(cell),plasticStrainOldLoc,ierr);CHKERRQ(ierr)
Expand Down