module FiniteDeformationClass use, intrinsic :: iso_fortran_env use MathClass use LinearSolverClass use FEMDomainClass use PostProcessingClass use ConstitutiveModelClass implicit none type:: FiniteDeform_ type(FEMDomain_),pointer ::FEMDomain real(real64),allocatable ::DeformStress(:,:,:) real(real64),allocatable ::DeformStrain(:,:,:) real(real64),allocatable ::DeformStressInit(:,:,:) real(real64),allocatable ::DeformStressinc(:,:,:) real(real64),allocatable ::DeformStressMat(:,:,:) real(real64),allocatable ::DeformStressRHS(:,:) real(real64),allocatable ::DeformVecEBETot(:,:) real(real64),allocatable ::DeformVecEBEInc(:,:) real(real64),allocatable ::DeformVecGloTot(:) real(real64),allocatable ::DeformVecGloInc(:) real(real64),allocatable ::TractionVecGlo(:) real(real64),allocatable ::ResidualVecGlo(:) real(real64),allocatable ::InternalVecGlo(:) real(real64),allocatable ::VolInitCurrEBE(:,:) real(real64),allocatable ::YoungsModulus(:) ! directly give parameter #1 real(real64),allocatable ::PoissonsRatio(:) ! directly give parameter #2 real(real64),allocatable :: PorePressure(:) ! directly give parameter #2 real(real64) ::dt,error,reactionforce real(real64) ::nr_tol=1.0e-8 logical :: ReducedIntegration = .false. logical :: infinitesimal = .false. integer(int32) :: itr integer(int32) :: Step=0 contains procedure :: Solve => SolveFiniteDeformNewton procedure :: UpdateSolution => SolveFiniteDeform procedure :: DivideBC => DevideBCIntoTimestep procedure :: UpdateBC => UpdateBCInTimestep procedure :: UpdateInitConfig => UpdateInitConfig procedure :: Setup => SetupFiniteDeform procedure :: Update => UpdateFiniteDeform procedure :: Display => DisplayDeformStress procedure :: getDBCVector => getDBCVectorDeform procedure :: getDispVector => getDispVectorDeform procedure :: getVolume => getVolumeDeform procedure :: check => checkFiniteDeform procedure :: import => importFiniteDeform procedure :: export => exportFiniteDeform procedure :: remove => removeFiniteDeform procedure :: save => saveFiniteDeform procedure :: open => openFiniteDeform procedure :: link => linkFiniteDeform end type contains subroutine removeFiniteDeform(obj) class(FiniteDeform_),intent(inout) :: obj if(associated(obj%FEMDomain) )then nullify(obj%FEMDomain) endif if( allocated(obj%DeformStress)) deallocate(obj%DeformStress) if( allocated(obj%DeformStrain)) deallocate(obj%DeformStrain) if( allocated(obj%DeformStressInit)) deallocate(obj%DeformStressInit) if( allocated(obj%DeformStressinc)) deallocate(obj%DeformStressinc) if( allocated(obj%DeformStressMat)) deallocate(obj%DeformStressMat) if( allocated(obj%DeformStressRHS)) deallocate(obj%DeformStressRHS) if( allocated(obj%DeformVecEBETot)) deallocate(obj%DeformVecEBETot) if( allocated(obj%DeformVecEBEInc)) deallocate(obj%DeformVecEBEInc) if( allocated(obj%DeformVecGloTot)) deallocate(obj%DeformVecGloTot) if( allocated(obj%DeformVecGloInc)) deallocate(obj%DeformVecGloInc) if( allocated(obj%TractionVecGlo)) deallocate(obj%TractionVecGlo) if( allocated(obj%ResidualVecGlo)) deallocate(obj%ResidualVecGlo) if( allocated(obj%InternalVecGlo)) deallocate(obj%InternalVecGlo) if( allocated(obj%VolInitCurrEBE)) deallocate(obj%VolInitCurrEBE) if( allocated(obj%YoungsModulus)) deallocate(obj%YoungsModulus) if( allocated(obj%PoissonsRatio)) deallocate(obj%PoissonsRatio) if( allocated(obj% PorePressure)) deallocate(obj% PorePressure) obj%dt=1.0d0 obj%error=0.0d0 obj%reactionforce=0.0d0 obj%nr_tol=1.0e-8 obj%ReducedIntegration = .false. obj%infinitesimal = .false. obj%itr=0 obj%Step=0 end subroutine ! ####################################################################### subroutine linkFiniteDeform(obj,FEMDomain) class(FiniteDeform_),intent(inout) :: obj type(FEMDomain_),target,intent(in) :: FEMDomain if(associated(obj%FEMDomain) )then nullify(obj%FEMDomain) obj%FEMDomain => FEMDOmain endif end subroutine ! ####################################################################### ! ####################################################################### subroutine openFiniteDeform(obj,path,name) class(FiniteDeform_),intent(inout) :: obj character(*),intent(in) :: path character(*),optional,intent(in) :: name character(200) :: pathi type(IO_) :: f integer(int32) :: n if(present(name) )then pathi=path !if( index(path, "/", back=.true.) == len(path) )then ! n=index(path, "/", back=.true.) ! pathi(n:n)= " " !endif call execute_command_line("mkdir -p "//trim(pathi)) call execute_command_line("mkdir -p "//trim(pathi)//"/"//trim(adjustl(name)) ) call f%open(trim(pathi)//"/"//trim(adjustl(name)) ,"/"//"FiniteDeform",".prop" ) else pathi=path !if( index(path, "/", back=.true.) == len(path) )then ! n=index(path, "/", back=.true.) ! pathi(n:n)= " " !endif call execute_command_line("mkdir -p "//trim(pathi)) call execute_command_line("mkdir -p "//trim(pathi)//"/FiniteDeform") call f%open(trim(pathi)//"/FiniteDeform","/FiniteDeform",".prop" ) endif ! write smt at here! call openArray(f%fh, obj%DeformStress) call openArray(f%fh, obj%DeformStrain) call openArray(f%fh, obj%DeformStressInit) call openArray(f%fh, obj%DeformStressinc) call openArray(f%fh, obj%DeformStressMat) call openArray(f%fh, obj%DeformStressRHS) call openArray(f%fh, obj%DeformVecEBETot) call openArray(f%fh, obj%DeformVecEBEInc) call openArray(f%fh, obj%DeformVecGloTot) call openArray(f%fh, obj%DeformVecGloInc) call openArray(f%fh, obj%TractionVecGlo) call openArray(f%fh, obj%ResidualVecGlo) call openArray(f%fh, obj%InternalVecGlo) call openArray(f%fh, obj%VolInitCurrEBE) call openArray(f%fh, obj%YoungsModulus) ! directly give parameter #1 call openArray(f%fh, obj%PoissonsRatio) ! directly give parameter #2 call openArray(f%fh, obj%PorePressure) ! directly give parameter #2 read(f%fh,*) obj%dt,obj%error,obj%reactionforce read(f%fh,*) obj%nr_tol read(f%fh,*) obj%ReducedIntegration read(f%fh,*) obj%infinitesimal call f%close() end subroutine ! ####################################################################### subroutine saveFiniteDeform(obj,path,name) class(FiniteDeform_),intent(inout) :: obj character(*),intent(in) :: path character(*),optional,intent(in) :: name character(200) :: pathi type(IO_) :: f integer(int32) :: n if(present(name) )then pathi=path !if( index(path, "/", back=.true.) == len(path) )then ! n=index(path, "/", back=.true.) ! pathi(n:n)= " " !endif call execute_command_line("mkdir -p "//trim(pathi)) call execute_command_line("mkdir -p "//trim(pathi)//"/"//trim(adjustl(name)) ) call f%open(trim(pathi)//"/"//trim(adjustl(name)) ,"/"//"FiniteDeform",".prop" ) else pathi=path !if( index(path, "/", back=.true.) == len(path) )then ! n=index(path, "/", back=.true.) ! pathi(n:n)= " " !endif call execute_command_line("mkdir -p "//trim(pathi)) call execute_command_line("mkdir -p "//trim(pathi)//"/FiniteDeform") call f%open(trim(pathi)//"/FiniteDeform","/FiniteDeform",".prop" ) endif ! write smt at here! call writeArray(f%fh, obj%DeformStress) call writeArray(f%fh, obj%DeformStrain) call writeArray(f%fh, obj%DeformStressInit) call writeArray(f%fh, obj%DeformStressinc) call writeArray(f%fh, obj%DeformStressMat) call writeArray(f%fh, obj%DeformStressRHS) call writeArray(f%fh, obj%DeformVecEBETot) call writeArray(f%fh, obj%DeformVecEBEInc) call writeArray(f%fh, obj%DeformVecGloTot) call writeArray(f%fh, obj%DeformVecGloInc) call writeArray(f%fh, obj%TractionVecGlo) call writeArray(f%fh, obj%ResidualVecGlo) call writeArray(f%fh, obj%InternalVecGlo) call writeArray(f%fh, obj%VolInitCurrEBE) call writeArray(f%fh, obj%YoungsModulus) ! directly give parameter #1 call writeArray(f%fh, obj%PoissonsRatio) ! directly give parameter #2 call writeArray(f%fh, obj%PorePressure) ! directly give parameter #2 write(f%fh,*) obj%dt,obj%error,obj%reactionforce write(f%fh,*) obj%nr_tol write(f%fh,*) obj%ReducedIntegration write(f%fh,*) obj%infinitesimal call f%close() end subroutine ! ####################################################################### ! ############################################################### subroutine importFiniteDeform(obj, YoungsModulus, PoissonsRatio, PorePressure) class(FiniteDeform_),intent(inout) :: obj real(real64),optional,intent(in) :: YoungsModulus(:), PoissonsRatio(:),& PorePressure(:) integer(int32) :: i,j,n if(present(YoungsModulus) )then if(allocated(obj%YoungsModulus) )then deallocate(obj%YoungsModulus) endif n=size(YoungsModulus) allocate(obj%YoungsModulus(n) ) obj%YoungsModulus(:)=YoungsModulus(:) endif if(present(PoissonsRatio) )then if(allocated(obj%PoissonsRatio) )then deallocate(obj%PoissonsRatio) endif n=size(PoissonsRatio) allocate(obj%PoissonsRatio(n) ) obj%PoissonsRatio(:)=PoissonsRatio(:) endif if(present(PorePressure) )then if(allocated(obj%PorePressure) )then deallocate(obj%PorePressure) endif n=size(PorePressure) allocate(obj%PorePressure(n) ) obj%PorePressure(:)=PorePressure(:) endif end subroutine importFiniteDeform ! ############################################################### ! ####################################################### subroutine checkFiniteDeform(obj) class(FiniteDeform_),intent(in) :: obj integer(int32) :: i,j,error_counter ! check conditions and return alartes. error_counter=0 print *, "checkFiniteDeform :: Analyzing imported data and checking consistency..." if(.not.associated(obj%FEMDomain) )then print *, "Alert! checkFiniteDeform >> FEMDomain is not improted." stop endif if(size(obj%FEMDomain%Boundary%DBoundNodID,2)/=size(obj%FEMDomain%Mesh%NodCoord,2))then print *, "size(obj%FEMDomain%Boundary%DBoundNodID,2)/=size(obj%FEMDomain%Mesh%NodCoord,2)" stop endif end subroutine ! ####################################################### !######################## Solve deformation by Netwon's method ######################## subroutine SolveFiniteDeformNewton(obj,OptionItr,Solvertype,nr_tol,infinitesimal,restart) class(FiniteDeform_),intent(inout)::obj integer(int32),optional,intent(in)::OptionItr character(*),optional,intent(in)::Solvertype real(real64),optional,intent(in)::nr_tol character*70 ::solver,defaultsolver logical,optional,intent(in) :: infinitesimal,restart real(real64),allocatable::Amat(:,:),bvec(:),xvec(:) real(real64)::val,er,residual,tolerance integer(int32) ::i,j,n,m,k,l,dim1,dim2,nodeid1,nodeid2,localid,itrmax,SetBC,itr_tol,itr integer(int32) :: dim_num,node_num,elem_num,node_num_elmtl character*70 :: gmsh,GaussJordan logical :: skip=.false. gmsh="Gmsh" GaussJordan="GaussJordan" if(present(infinitesimal) )then obj%infinitesimal = infinitesimal endif if(present(OptionItr) )then itr_tol = OptionItr else itr_tol = obj%FEMDomain%ControlPara%ItrTol endif if(present(nr_tol) )then obj%nr_tol=nr_tol endif if(present(restart) )then if(restart .eqv. .true.)then skip=.true. endif endif do itr=itr + 1 obj%Itr=itr if(itr==1 .and. skip .eqv. .false.)then call SetupFiniteDeform(obj) if(obj%Step==1)then !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=0) endif call SolveFiniteDeform(obj,SolverType=Solvertype) !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=itr) if(obj%infinitesimal .eqv. .true.)then exit endif else call UpdateFiniteDeform(obj) !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=0) call SolveFiniteDeform(obj,OptionItr=itr,SolverType=Solvertype) !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=itr) !print *, "Residual r*r = ",dot_product(obj%ResidualVecGlo,obj%ResidualVecGlo) print *, "Step, Itr, ERROR :: ",obj%step ,obj%itr,obj%error if(obj%error<obj%nr_tol )then print *,"itr=",itr, "Netwton's Method is converged !" exit endif endif if(itr == itr_tol)then print *, "SolveFiniteDeformNewton >> Newton's method did not converge" stop "debug" exit endif enddo call UpdateStressMeasure(obj) call UpdateStrainMeasure(obj) end subroutine !######################## Solve deformation by Netwon's method ######################## !########################## Initialize Boundary Conditions of Finite Deformation ################################### subroutine DevideBCIntoTimestep(obj) class(FiniteDeform_),intent(inout)::obj integer(int32) ::n,m, timestep !debug Display Dirichlet Boundary Condition !call obj%FEMDomain%GmshPlotMesh(onlyDirichletBC=.true.) timestep=obj%FEMDomain%ControlPara%Timestep if(obj%FEMDomain%ControlPara%SimMode==1)then n=size(obj%FEMDomain%Boundary%DBoundVal,1) m=size(obj%FEMDomain%Boundary%DBoundVal,2) if(.not.allocated(obj%FEMDomain%Boundary%DBoundValInc ) )then allocate(obj%FEMDomain%Boundary%DBoundValInc(n,m) ) else if( size(obj%FEMDomain%Boundary%DBoundValInc,1)/=n .or. & size(obj%FEMDomain%Boundary%DBoundValInc,2)/=m)then deallocate(obj%FEMDomain%Boundary%DBoundValInc) endif endif obj%FEMDomain%Boundary%DBoundValInc(:,:)=1.0d0/dble(timestep)*obj%FEMDomain%Boundary%DBoundVal(:,:) !obj%FEMDomain%Boundary%DBoundVal(:,:)=obj%FEMDomain%Boundary%DBoundValInc(:,:) elseif(obj%FEMDomain%ControlPara%SimMode==2)then n=size(obj%FEMDomain%Boundary%NBoundVal,1) m=size(obj%FEMDomain%Boundary%NBoundVal,2) if(.not.allocated(obj%FEMDomain%Boundary%NBoundValInc ) )then allocate(obj%FEMDomain%Boundary%NBoundValInc(n,m) ) else if( size(obj%FEMDomain%Boundary%NBoundValInc,1)/=n .or. & size(obj%FEMDomain%Boundary%NBoundValInc,2)/=m)then deallocate(obj%FEMDomain%Boundary%DBoundValInc) endif endif obj%FEMDomain%Boundary%NBoundValInc(:,:)=1.0d0/dble(timestep)*obj%FEMDomain%Boundary%NBoundVal(:,:) obj%FEMDomain%Boundary%NBoundVal(:,:)=obj%FEMDomain%Boundary%NBoundValInc(:,:) else print *, "ERROR :: Displacement Control or Force Control?" endif end subroutine !########################## Initialize Boundary Conditions of Finite Deformation ################################### !########################## Initialize Boundary Conditions of Finite Deformation ################################### subroutine UpdateBCInTimestep(obj) class(FiniteDeform_),intent(inout)::obj integer(int32) :: timestep if(obj%FEMDomain%ControlPara%SimMode==1)then !obj%FEMDomain%Boundary%DBoundVal(:,:)=obj%FEMDomain%Boundary%DBoundValInc(:,:) elseif(obj%FEMDomain%ControlPara%SimMode==2)then !obj%FEMDomain%Boundary%NBoundVal(:,:)=obj%FEMDomain%Boundary%NBoundValInc(:,:) else print *, "ERROR :: Displacement Control or Force Control?" endif end subroutine !########################## Initialize Boundary Conditions of Finite Deformation ################################### !############################################################# subroutine ImportFEMDomainFiDe(obj,OptionalFileFormat,OptionalProjectName) class(FEMDomain_),intent(inout)::obj character*4,optional,intent(in)::OptionalFileFormat character(*),optional,intent(in)::OptionalProjectName character*4::FileFormat character*70::ProjectName character*74 ::FileName integer(int32),allocatable::IntMat(:,:) real(real64),allocatable::RealMat(:,:) integer(int32) :: fh,i,j,NumOfDomain,n,m,DimNum character*70 Msg call DeallocateFEMDomain(obj) fh =104 if(present(OptionalFileFormat) )then FileFormat=trim(OptionalFileFormat) else FileFormat=".scf" endif if(present(OptionalProjectName) )then ProjectName=trim(OptionalProjectName) else ProjectName="untitled" endif FileName = trim(ProjectName)//trim(FileFormat) obj%FileName = trim(FileName) obj%FilePath = "./" !!print *, "Project : ",ProjectName !!print *, "is Exported as : ",FileFormat," format" !!print *, "File Name is : ",FileName open(fh,file=FileName,status="old") if(trim(FileFormat)==".scf" )then read(fh,*) NumOfDomain allocate(IntMat(NumOfDomain,2)) allocate(obj%Mesh%SubMeshNodFromTo(NumOfDomain,3) ) allocate(obj%Mesh%SubMeshElemFromTo(NumOfDomain,3) ) do i=1,NumOfDomain obj%Mesh%SubMeshNodFromTo(i,1)=i read(fh,*) obj%Mesh%SubMeshNodFromTo(i,2),obj%Mesh%SubMeshNodFromTo(i,3) enddo do i=1,NumOfDomain obj%Mesh%SubMeshElemFromTo(i,1)=i read(fh,*) obj%Mesh%SubMeshElemFromTo(i,3) if(i==1)then obj%Mesh%SubMeshElemFromTo(i,2)=1 else obj%Mesh%SubMeshElemFromTo(i,2)=obj%Mesh%SubMeshElemFromTo(i-1,3) + 1 endif enddo read(fh,*) n,m DimNum=m allocate(obj%Mesh%NodCoord(n,m) ) call ImportArray(obj%Mesh%NodCoord,OptionalFileHandle=fh) call CopyArray(obj%Mesh%NodCoord,obj%Mesh%NodCoordInit ) read(fh,*) n,m read(fh,*)obj%Mesh%ElemType obj%ShapeFunction%ElemType=obj%Mesh%ElemType allocate(obj%Mesh%ElemNod(n,m) ) allocate(obj%Mesh%ElemMat(n ) ) call ImportArray(obj%Mesh%ElemNod,OptionalFileHandle=fh) do i=1,n read(fh,*) obj%Mesh%ElemMat(i) enddo read(fh,*) n,m allocate(obj%MaterialProp%MatPara(n,m) ) call ImportArray(obj%MaterialProp%MatPara,OptionalFileHandle=fh) !######### Dirichlet boundary conditions ################# DimNum=size(obj%Mesh%NodCoord,2) allocate(obj%Boundary%DBoundNum(DimNum )) read(fh,*) obj%Boundary%DBoundNum(:) allocate(obj%Boundary%DBoundNodID( maxval(obj%Boundary%DBoundNum), size(obj%Boundary%DBoundNum) ) ) allocate(obj%Boundary%DBoundVal( maxval(obj%Boundary%DBoundNum), size(obj%Boundary%DBoundNum) ) ) obj%Boundary%DBoundNodID(:,:)=-1 obj%Boundary%DBoundVal(:,:) =0.0d0 do i=1,size(obj%Boundary%DBoundNum,1) do j=1,obj%Boundary%DBoundNum(i) read(fh,*) obj%Boundary%DBoundNodID(j,i) enddo do j=1,obj%Boundary%DBoundNum(i) read(fh,*) obj%Boundary%DBoundVal(j,i) enddo enddo !######### Dirichlet boundary conditions ################# !######### Neumann boundary conditions ################# read(fh,*) n if(n==0)then allocate( obj%Boundary%NBoundNum(0)) allocate(obj%Boundary%NBoundNodID(0,0) ) allocate(obj%Boundary%NBoundVal( 0,0) ) else allocate( obj%Boundary%NBoundNum(DimNum)) allocate(obj%Boundary%NBoundNodID(n, size(obj%Boundary%NBoundNum) ) ) allocate(obj%Boundary%NBoundVal( n, size(obj%Boundary%NBoundNum) ) ) obj%Boundary%NBoundNodID(:,:)=-1 obj%Boundary%NBoundVal(:,:) =0.0d0 obj%Boundary%NBoundNum(:)=n do i=1,n read(fh,*) m obj%Boundary%NBoundNodID(i,:)=m enddo do i=1,n read(fh,*) obj%Boundary%NBoundVal(i,:) enddo endif !######### Neumann boundary conditions ################# !######### Initial conditions ################# read(fh,*) n if(n==0)then allocate( obj%Boundary%TBoundNum(0)) allocate(obj%Boundary%TBoundNodID(0,0) ) allocate(obj%Boundary%TBoundVal( 0,0) ) else allocate(obj%Boundary%TBoundNodID(n,1) ) allocate(obj%Boundary%TBoundVal( n,1) ) allocate(obj%Boundary%TBoundNum( 1) ) obj%Boundary%TBoundNum=n if(n/=0)then if(n<0)then print *, "ERROR :: number of initial conditions are to be zero" else do i=1,n read(fh,*) obj%Boundary%TBoundNodID(i,1) enddo do i=1,n read(fh,*) obj%Boundary%TBoundVal(i,1) enddo endif endif endif !######### Initial conditions ################# read(fh,*) obj%ControlPara%SimMode, obj%ControlPara%ItrTol,obj%ControlPara%Timestep close(fh) else !!print *, "ERROR :: ExportFEMDomain >> only .scf file can be exported." endif end subroutine !############################################################# !############################################################# subroutine SetupFiniteDeform(obj,tol) class(FiniteDeform_),intent(inout)::obj integer(int32),optional,intent(in) :: tol if(obj%dt==0.0d0 .or. obj%dt/=obj%dt)then obj%dt=1.0d0 endif !if(present(tol) )then ! obj%nr_tol = tol !else ! obj%nr_tol = 1.0e-08 !endif call UpdateCurrConfig(obj) call GetDeformStressMatAndVector(obj) end subroutine !############################################################# !############################################################# subroutine UpdateFiniteDeform(obj,restart) class(FiniteDeform_),intent(inout)::obj logical,optional,intent(in) :: restart call UpdateCurrConfig(obj) call GetDeformStressMatAndVector(obj) end subroutine !############################################################# !############################################################# subroutine UpdateCurrConfig(obj,restart) class(FiniteDeform_),intent(inout)::obj logical,optional,intent(in) :: restart integer(int32) :: i,j,n,m integer(int32) :: num_node,num_elem,num_dim num_node=size(obj%FEMDomain%Mesh%NodCoord,1) num_dim =size(obj%FEMDomain%Mesh%NodCoord,2) num_elem=size(obj%FEMDomain%Mesh%ElemNod,1) ! Displacement :: u if(.not.allocated(obj%DeformVecGloTot) ) then allocate(obj%DeformVecGloTot(num_node*num_dim) ) obj%DeformVecGloTot(:)=0.0d0 endif ! ⊿u = v if(.not.allocated(obj%DeformVecGloInc) ) then allocate(obj%DeformVecGloInc(num_node*num_dim) ) obj%DeformVecGloInc(:)=0.0d0 endif if(.not. allocated(obj%FEMDomain%Mesh%NodCoordInit) )then allocate(obj%FEMDomain%Mesh%NodCoordInit(num_node,num_dim) ) obj%FEMDomain%Mesh%NodCoordInit(:,:)=obj%FEMDomain%Mesh%NodCoord(:,:) endif !call showArraySIze(obj%FEMDomain%Mesh%NodCoord) !call showArraySIze(obj%FEMDomain%Mesh%NodCoordInit) !call showArraySize(obj%DeformVecGloInc) !call showArraySize(obj%DeformVecGloTot) !stop do i=1,num_node do j=1,num_dim obj%FEMDomain%Mesh%NodCoord(i,j )=obj%FEMDomain%Mesh%NodCoordInit(i,j) + & obj%DeformVecGloInc(num_dim*(i-1) + j) + obj%DeformVecGloTot(num_dim*(i-1) + j) enddo enddo end subroutine !############################################################# !############################################################# subroutine UpdateInitConfig(obj) class(FiniteDeform_),intent(inout)::obj integer(int32) :: i,j,n,m integer(int32) :: num_node,num_elem,num_dim num_node=size(obj%FEMDomain%Mesh%NodCoord,1) num_dim =size(obj%FEMDomain%Mesh%NodCoord,2) num_elem=size(obj%FEMDomain%Mesh%ElemNod,1) if(.not.allocated(obj%DeformVecGloTot) ) then allocate(obj%DeformVecGloTot(num_node*num_dim) ) obj%DeformVecGloTot(:)=0.0d0 endif if(.not.allocated(obj%DeformVecGloInc) ) then allocate(obj%DeformVecGloInc(num_node*num_dim) ) obj%DeformVecGloInc(:)=0.0d0 endif do i=1,num_node do j=1,num_dim ! Updated Lagrangian obj%FEMDomain%Mesh%NodCoordInit(i,j )=obj%FEMDomain%Mesh%NodCoord(i,j) enddo enddo end subroutine !############################################################# !================================================================================ ! �S�̍����}�g���N�X�̌v�Z subroutine GetDeformStressMatAndVector(obj,OptionalStep,restart) class(FiniteDeform_),intent(inout)::obj type(ConstModel_) ::mdl logical,optional,intent(in) :: restart integer(int32), optional,intent(in) :: OptionalStep integer(int32) :: itr_rm,itr,itr_contact integer(int32) :: nod_num,dim_num,elemnod_num,elem_num real(real64), allocatable :: g(:,:), Bmat(:,:), Ce_neoHK(:,:), BTmat(:,:), Psymat(:,:), & xymat(:,:),xymat_c(:,:), Jmat(:,:), s(:), Kmat_e(:,:),c_nod_coord(:,:),cc_nod_coord(:,:),& dNdxi(:,:),F_iJ_n(:,:),F_iJ(:,:),C_IJ(:,:),Cp_IJ_n(:,:),Cp_IJ(:,:),b_ij(:,:),F_inv(:,:),& F_T_inv(:,:),F_T(:,:),dNdx(:,:),M_IJ(:,:),Cp_IJ_inv(:,:),gvec_e(:) integer(int32),allocatable::ij(:,:) integer(int32) i, j,k, m, p, q,gp_num,NumOfStrainMeasure real(real64) detJ,Lamda,mu,c,phy,psy,E,v,xx,Tol real(real64) ,allocatable::MatPara(:) itr_rm =obj%FEMDomain%ControlPara%ItrTol itr =obj%itr itr_contact=obj%FEMDomain%ControlPara%ItrTol Tol =obj%FEMDomain%ControlPara%Tol !if( abs(tol)<1.0e-15 )then tol = 1.0e-15 !endif nod_num = size(obj%FEMDomain%Mesh%NodCoord,1) elem_num = size(obj%FEMDomain%Mesh%ElemNod,1) elemnod_num = size(obj%FEMDomain%Mesh%ElemNod,2) dim_num = size(obj%FEMDomain%Mesh%NodCoord,2) if( .not. allocated( obj%DeformVecEBEInc) ) then allocate(obj%DeformVecEBEInc(elem_num,elemnod_num*dim_num) ) obj%DeformVecEBEInc(:,:)=0.0d0 endif if( .not. allocated( obj%DeformVecEBETot) ) then allocate(obj%DeformVecEBETot(elem_num,elemnod_num*dim_num) ) obj%DeformVecEBETot=0.0d0 endif if( .not. allocated( obj%DeformStressRHS) )then allocate(obj%DeformStressRHS(elem_num,elemnod_num*dim_num) ) obj%DeformStressRHS(:,:)=0.0d0 endif if( .not. allocated( obj%DeformStressMat) ) then allocate(obj%DeformStressMat(elem_num,elemnod_num*dim_num,elemnod_num*dim_num) ) obj%DeformStressMat(:,:,:)=0.0d0 endif if( .not. allocated( obj%VolInitCurrEBE) ) then allocate(obj%VolInitCurrEBE(elem_num,3) ) !initial, current obj%VolInitCurrEBE(:,:)=1.0d0 endif obj%DeformStressRHS(:,:) =0.0d0 obj%DeformStressMat(:,:,:)=0.0d0 obj%FEMDomain%ShapeFunction%ElemType=obj%FEMDomain%Mesh%GetElemType() call SetShapeFuncType(obj%FEMDomain%ShapeFunction, ReducedIntegration=obj%ReducedIntegration) gp_num=obj%FEMDomain%ShapeFunction%NumOfGp if(dim_num<=2)then NumOfStrainMeasure=15 else NumOfStrainMeasure=49 endif allocate(MatPara(7)) if( .not. allocated( obj%DeformStrain )) then allocate(obj%DeformStrain(elem_num,gp_num,NumOfStrainMeasure) ) obj%DeformStrain(:,:,:) =0.0d0 obj%DeformStrain(:,:,1) =1.0d0 ! Cp_iJ_n(1,1) obj%DeformStrain(:,:,2) =1.0d0 ! Cp_iJ_n(2,2) obj%DeformStrain(:,:,4) =1.0d0 ! Cp_iJ(1,1) obj%DeformStrain(:,:,5) =1.0d0 ! Cp_iJ(2,2) obj%DeformStrain(:,:,7) =1.0d0 ! F_iJ_n(1,1) obj%DeformStrain(:,:,8) =1.0d0 ! F_iJ_n(2,2) obj%DeformStrain(:,:,11)=1.0d0 ! F_iJ(1,1) obj%DeformStrain(:,:,12)=1.0d0 ! F_iJ(2,2) if(NumOfStrainMeasure==49)then obj%DeformStrain(:,:,16) = 1.0d0 ! Cp_iJ_n(3,3) obj%DeformStrain(:,:,19) = 1.0d0 ! Cp_iJ(3,3) obj%DeformStrain(:,:,22) = 1.0d0 ! F_iJ_n(3,3) obj%DeformStrain(:,:,27) = 1.0d0 ! F_iJ(3,3) endif endif !initialize if(.not.allocated(obj%DeformStress))then allocate(obj%DeformStress(elem_num,gp_num,6) ) obj%DeformStress(:,:,:)=0.0d0 endif if(.not.allocated(obj%DeformStressinc))then allocate(obj%DeformStressinc(elem_num,gp_num,6) ) obj%DeformStressinc(:,:,:)=0.0d0 endif m = elemnod_num*dim_num allocate (Kmat_e(elemnod_num*dim_num,elemnod_num*dim_num),gvec_e(elemnod_num*dim_num)) obj%VolInitCurrEBE(:,2)=0.0d0 do i = 1, elem_num !�v�f���ƃ��[� Kmat_e(:,:) = 0.0d0 gvec_e(:) = 0.0d0 E = obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i),1) v = obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i),2) ! allow direct-import of Youngs modulus and Poissons ratio if(allocated(obj%YoungsModulus) )then if(size(obj%YoungsModulus)/=elem_num) then print *, "ERROR :: FiniteDeform :: size(obj%YoungsModulus/=elem_num)" else E = obj%YoungsModulus(i) endif endif ! allow direct-import of Youngs modulus and Poissons ratio if(allocated(obj%PoissonsRatio) )then if(size(obj%PoissonsRatio)/=elem_num) then print *, "ERROR :: FiniteDeform :: size(obj%PoissonsRatio/=elem_num)" else v = obj%PoissonsRatio(i) endif endif mdl%Lamda=v*E/(1.0d0 + v)/(1.0d0-2.0d0*v) mdl%mu=E/2.0d0/(1.0d0 + v) c=obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i),4) phy=obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i),5) psy=obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i),6) MatPara(1)=E MatPara(2)=v MatPara(3)=Lamda MatPara(4)=mu MatPara(5)=c MatPara(6)=phy MatPara(7)=psy do j = 1, obj%FEMDomain%ShapeFunction%NumOfGp !�K�E�X�ϕ����ƃ��[�v ! -----J�}�g���N�X�̌v�Z----------------------------------------- call GetAllShapeFunc(obj%FEMDomain%ShapeFunction,elem_id=i,nod_coord=obj%FEMDomain%Mesh%NodCoord,& elem_nod=obj%FEMDomain%Mesh%ElemNod,OptionalGpID=j,ReducedIntegration=obj%ReducedIntegration) if(obj%infinitesimal .eqv. .false.)then ! ---Compute stress/strain measures of the Finite Elasto-Plast- call F_tensor_ICU(obj,i,j,mdl%F_iJ_n,mdl%F_iJ) call trans_rank_2(mdl%F_iJ,mdl%F_T) call C_tensor(mdl%F_iJ,mdl%C_IJ,mdl%b_ij,itr,dim_num) !call Cp_tensor(i,j,obj%DeformStrain,mdl%Cp_IJ_n,mdl%Cp_IJ,mdl%Cp_IJ_inv,dim_num) call inverse_rank_2(mdl%F_iJ,mdl%F_inv_iJ) call inverse_rank_2(mdl%F_T,mdl%F_T_inv_iJ) !call M_neo_Hookean(mdl%C_IJ,mdl%Cp_IJ,mdl%Cp_IJ_inv,mdl%M_IJ,mdl%Lamda,mdl%mu,i,j) !----Return-Mapping algorithm---------------------------------- !call Return_Mapping_MCDP(dim_num,i,j,mdl%C_IJ,mdl%Cp_IJ,mdl%Cp_IJ_n,mdl%Cp_IJ_inv,mdl%M_IJ,MatPara,& !itr_rm,tol,obj%DeformStress,mdl%F_T,mdl%F_T_inv_ij,itr,itr_contact,obj%DeformStrain,OptionalStep) ! -----current config. Ce matrix----------------------------------------- !call Ce_neoHK_current(dim_num, i,j,mdl%Lamda,mdl%mu,mdl%C_IJ,mdl%Cp_IJ,mdl%b_ij,mdl%M_IJ,Ce_neoHK,mdl%F_T,mdl%F_T_inv,ij) ! -----B,W�}�g���N�X�̌v�Z---------------------------------------- call B_mat(dim_num,obj%FEMDomain%ShapeFunction%dNdgzi, & obj%FEMDomain%ShapeFunction%Jmat, obj%FEMDomain%ShapeFunction%detJ, mdl%Bmat,m) ! -----BT�}�g���N�X�̌v�Z----------------------------------------- !call trans_rank_2(Bmat,BTmat) !-----2-D neo-Hookean current stifness matrix------------------- ! Get Volumetric change obj%VolInitCurrEBE(i,2)=obj%VolInitCurrEBE(i,2)+& det_mat(mdl%F_iJ,size(mdl%F_iJ,1) )/dble(obj%FEMDomain%ShapeFunction%NumOfGp) mdl%ModelType="NeoHookean" call GetKmat(obj,mdl,obj%FEMDomain%ShapeFunction,Kmat_e,gvec_e,dim_num,elemnod_num,i,j) else ! infinitesimal strain theory ! ---Compute stress/strain measures of the Finite Elasto-Plast- call F_tensor_ICU(obj,i,j,mdl%F_iJ_n,mdl%F_iJ) ! -----B,W�}�g���N�X�̌v�Z---------------------------------------- call B_mat(dim_num,obj%FEMDomain%ShapeFunction%dNdgzi, & obj%FEMDomain%ShapeFunction%Jmat, obj%FEMDomain%ShapeFunction%detJ, mdl%Bmat,m) call GetKmat(obj,mdl,obj%FEMDomain%ShapeFunction,Kmat_e,gvec_e,dim_num,elemnod_num,i,j) endif !do k=1,size(Kmat_e,1) ! write(200,*) Kmat_e(k,:),"|",gvec_e(k) !enddo !stop "debug" !call K_mat_e(j,obj%FEMDomain%ShapeFunction%GaussIntegWei, & ! BTmat, Ce_neoHK, Bmat, obj%FEMDomain%ShapeFunction%detJ, Kmat_e,F_iJ) !call GetGvec(obj,mdl,obj%FEMDomain%ShapeFunction,gvec_e,dim_num,elemnod_num,i) !call g_vector_e(i,j,obj%FEMDomain%ShapeFunction%GaussIntegWei,& ! BTmat,obj%DeformStress, obj%FEMDomain%ShapeFunction%detJ, gvec_e) ! -----�v�f�����}�g���N�X�̑�������-------------------------------- enddo ! ------�S�̍����}�g���N�X�ւ̏d�ˍ��킹------------------------------ call K_mat_ICU(obj%DeformStressMat, obj%FEMDomain%Mesh%ElemNod, i, Kmat_e) call g_vector_ICU(i,obj%FEMDomain%Mesh%ElemNod,gvec_e,obj%DeformStressRHS) enddo end subroutine !================================================================================ !================================================================================ subroutine GetKmat(obj,mdl,sf,Kmat_e,gvec_e,dim_num,elemnod_num,elem,gp) type(FiniteDeform_),intent(inout) :: obj type(ConstModel_),intent(inout)::mdl type(ShapeFunction_),intent(in)::sf real(real64),intent(inout) :: Kmat_e(:,:),gvec_e(:) integer(int32),intent(in)::dim_num,elemnod_num,elem,gp real(real64):: a0(3,3),dXdgzi(3,3) real(real64):: Xmat(elemnod_num,3) real(real64),allocatable:: a0_inv(:,:), dgzidX(:,:),Jgzimat_inv(:,:),Dmat(:,:),Sigma(:),Sigma_i(:,:),& dumat_i(:,:),Sigma_tot(:,:) real(real64):: a1(elemnod_num,3) real(real64):: K1(elemnod_num*dim_num,elemnod_num*dim_num) real(real64):: a2(elemnod_num,3) real(real64):: A3(elemnod_num,elemnod_num) real(real64)::dumat(elemnod_num,3) real(real64):: a4(elemnod_num,3,elemnod_num,3) real(real64)::Jgzimat(3,3),dxxdgzi(3,3),dNdX(elemnod_num,3) integer(int32) ::alpha,beta,ganma,ii,I,J,jj,K,kk,L,ll,n,m,o,p,q,rr,R,s,loc_1,loc_2 real(real64) :: detJgzi real(real64) :: PorePressure character*70 :: DerType !DerType="F_iJ" DerType="c_current" allocate(dumat_i(elemnod_num*3,1) ) !mdl%infinitesimal = obj%infinitesimal n=size(sf%ElemCoord,1) m=size(sf%ElemCoord,2) do i=1,n !num of node per elems do j=1,m ! num of dim dumat(i,j)=obj%DeformVecEBEInc(elem, m*(i-1) + j ) dumat_i( (i-1)*3 + j ,1) =dumat(i,j) enddo enddo !Xmat(:,:) =sf%ElemCoord(:,:) !- dumat(:,:) dXdgzi(:,:)=matmul(sf%dNdgzi,sf%ElemCoord) !dXdgzi(:,:)=matmul(sf%dNdgzi,Xmat) call inverse_rank_2(dXdgzi,dgzidX) dNdX(:,:) =matmul(transpose(sf%dNdgzi),dgzidX ) !a0(:,:)=matmul(sf%dNdgzi,sf%ElemCoord) !dxdgzi Jgzimat(:,:)=matmul(mdl%F_iJ,dXdgzi) call inverse_rank_2(Jgzimat,Jgzimat_inv) detJgzi=det_mat(Jgzimat,size(Jgzimat,1) ) if(obj%infinitesimal .eqv. .false.)then !call inverse_rank_2(a0,a0_inv) call HyperElasticStress(mdl) call HyperElasticDer(mdl,DerType) !get cc mat call GetDmat(Dmat,mdl%StressDer,dim_num) call GetSigmaVec(Sigma,mdl%sigma,dim_num) K1(:,:)=0.0d0 A3(:,:)=matmul( matmul( dNdX , mdl%sigma ) , transpose(dNdX)) do i=1,elemnod_num do j=1,elemnod_num K1( (I-1)*dim_num+1,(J-1)*dim_num+1 ) =K1( (I-1)*dim_num+1,(J-1)*dim_num+1 ) + A3(I,J) K1( (I-1)*dim_num+2,(J-1)*dim_num+2 ) =K1( (I-1)*dim_num+2,(J-1)*dim_num+2 ) + A3(I,J) K1( (I-1)*dim_num+3,(J-1)*dim_num+3 ) =K1( (I-1)*dim_num+3,(J-1)*dim_num+3 ) + A3(I,J) enddo enddo ! debug !print *, "size(obj%PorePressure)" ,size(obj%PorePressure) !write(3001,*) maxval(obj%PorePressure),minval(obj%PorePressure) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! For water-soil coupling analysis if(.not. allocated(obj%PorePressure) )then PorePressure=0.0d0 elseif(size(obj%PorePressure) ==0 )then PorePressure=0.0d0 else PorePressure=obj%PorePressure(elem) endif Sigma(1)=Sigma(1)-PorePressure/dble(obj%FEMDomain%ShapeFunction%NumOfGp) Sigma(2)=Sigma(2)-PorePressure/dble(obj%FEMDomain%ShapeFunction%NumOfGp) Sigma(3)=Sigma(3)-PorePressure/dble(obj%FEMDomain%ShapeFunction%NumOfGp) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Kmat_e(:,:)=Kmat_e(:,:)+matmul(matmul(transpose(mdl%Bmat),Dmat),mdl%Bmat )/det_mat(mdl%F_ij,size(mdl%F_ij,1) )*sf%detJ & +K1(:,:)*sf%detJ !stop "debug l3" gvec_e(:)=gvec_e(:)+matmul(transpose(mdl%Bmat),Sigma)*sf%detJ else !call inverse_rank_2(a0,a0_inv) !call HyperElasticStress(mdl) !call HyperElasticDer(mdl,DerType) !get cc mat !call GetDmat(Dmat,mdl%StressDer,dim_num) !call GetSigmaVec(Sigma,mdl%sigma,dim_num) ! clear Dmat and import St. Venant stiffness matrix allocate(Dmat(6,6) ) Dmat(:,:)=0.0d0 Dmat(1,1)= 2.0d0*mdl%mu + mdl%lamda Dmat(1,2)= mdl%lamda Dmat(1,3)= mdl%lamda Dmat(2,1)= mdl%lamda Dmat(2,2)= 2.0d0*mdl%mu + mdl%lamda Dmat(2,3)= mdl%lamda Dmat(3,1)= mdl%lamda Dmat(3,2)= mdl%lamda Dmat(3,3)= 2.0d0*mdl%mu + mdl%lamda Dmat(4,4)= mdl%mu Dmat(5,5)= mdl%mu Dmat(6,6)= mdl%mu ! stress dsigma = [D][B]{du} allocate(Sigma_i(6,1)) allocate(Sigma_tot(6,1)) Sigma_i(:,:) = matmul(Dmat, matmul(mdl%Bmat, dumat_i) ) !\partial \sigma / \artial t * dt Sigma_tot(1:6,1) =obj%DeformStress(elem,gp,1:6) ! debug !print *, "size(obj%PorePressure)" ,size(obj%PorePressure) !write(3001,*) maxval(obj%PorePressure),minval(obj%PorePressure) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! For water-soil coupling analysis if(.not. allocated(obj%PorePressure) )then PorePressure=0.0d0 elseif(size(obj%PorePressure) ==0 )then PorePressure=0.0d0 else PorePressure=obj%PorePressure(elem) endif Sigma_tot(1,1)=Sigma_tot(1,1)-PorePressure!/dble(obj%FEMDomain%ShapeFunction%NumOfGp) Sigma_tot(2,1)=Sigma_tot(2,1)-PorePressure!/dble(obj%FEMDomain%ShapeFunction%NumOfGp) Sigma_tot(3,1)=Sigma_tot(3,1)-PorePressure!/dble(obj%FEMDomain%ShapeFunction%NumOfGp) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Kmat_e(:,:)=Kmat_e(:,:)+matmul(matmul(transpose(mdl%Bmat),Dmat),mdl%Bmat )*sf%detJ !stop "debug l3" gvec_e(:)=gvec_e(:)+matmul(transpose(mdl%Bmat),Sigma_tot(:,1) )*sf%detJ& +matmul(transpose(mdl%Bmat),Sigma_i(:,1) )*sf%detJ ! update d \sigma obj%DeformStressinc(elem,gp,1:6)=Sigma_i(1:6,1) endif ! !do I=1,elemnod_num ! do ii=1,3 ! do K=1,elemnod_num ! do kk=1,3 ! do R=1,3 ! do j=1,3 ! a4(II,ii,K,kk)=a4(II,ii,K,kk)& ! +a1(I,j)*mdl%StressDer(ii,j,kk,R)*dNdX(K,R)*detJgzi ! enddo ! enddo ! enddo ! enddo ! enddo !enddo ! !do I=1,elemnod_num ! do ii=1,3 ! do K=1,elemnod_num ! do kk=1,3 ! loc_1=(I - 1)*3+ii ! loc_2=(K - 1)*3+kk ! a0(:,:)=matmul(mdl%F_inv_iJ,mdl%sigma) ! a1(:,:)=matmul(dNdX,a0) ! a2(:,:)=matmul(dNdX,mdl%F_inv_iJ) ! Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)& ! -a1(K,ii)*a2(I,kk)*detJgzi ! ! a0(:,:)=matmul(mdl%F_inv_iJ,mdl%sigma) ! a1(:,:)=matmul(dNdX,a0) ! a2(:,:)=matmul(dNdX,dXdgzi) ! a3(:,:)=matmul(a2,Jgzimat_inv) ! Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)& ! +a1(I,ii)*a3(K,kk)*detJgzi ! ! Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)& ! +a4(II,ii,K,kk)*detJgzi ! enddo ! enddo ! enddo !enddo ! !do i=1,size(Kmat_e,1) ! write(10,*) Kmat_e(i,:) !enddo !write (10,*) " " ! do alpha=1,elemnod_num ! do ganma=1,3 ! do I=1,elemnod_num ! do ii=1,3 ! do m=1,3 ! do j=1,3 ! a4(alpha,ganma,I,ii)=a4(alpha,ganma,I,ii)& ! +a1(alpha,j)*a3(I,m)*mdl%StressDer(ii,j,ganma,m)*sf%detJ ! enddo ! enddo ! enddo ! enddo ! enddo ! enddo ! ! do alpha=1,elemnod_num ! do ganma=1,3 ! do I=1,elemnod_num ! do ii=1,3 ! loc_1=(alpha-1)*3+ii ! loc_2=(I - 1)*3+ganma ! Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)& ! -a1(alpha,ganma)*a2(I,ii)*sf%detJ& ! +a2(alpha,ii)*a1(I,ganma)*sf%detJ& ! -a4(alpha,ganma,I,ii) ! enddo ! enddo ! enddo ! enddo ! end subroutine !================================================================================= !================================================================================ subroutine GetGvec(obj,mdl,sf,gvec_e,dim_num,elemnod_num,elem) type(FiniteDeform_),intent(in) :: obj type(ConstModel_),intent(inout)::mdl type(ShapeFunction_),intent(in)::sf real(real64),intent(inout) :: gvec_e(:) integer(int32),intent(in)::dim_num,elemnod_num,elem real(real64):: a0(3,3),Jgzimat(3,3),dXdgzi(3,3) real(real64),allocatable:: a0_inv(:,:),dgzidX(:,:) real(real64):: a1(elemnod_num,3) real(real64):: a2(elemnod_num,3) real(real64):: dumat(elemnod_num,3) real(real64):: Xmat(elemnod_num,3) real(real64):: a3(elemnod_num,3) real(real64):: dNdX(elemnod_num,3) real(real64):: a4(elemnod_num,3,elemnod_num,3),detJgzi integer(int32) ::alpha,beta,ganma,ii,I,j,k,l,m,o,p,q,r,s,loc_1,n character*70 :: DerType DerType="F_iJ" n=size(sf%ElemCoord,1) m=size(sf%ElemCoord,2) do i=1,n !num of node per elems do j=1,m ! num of dim dumat(i,j)=obj%DeformVecEBEInc(elem, m*(i-1) + j ) enddo enddo Xmat(:,:) =sf%ElemCoord(:,:) - dumat(:,:) dXdgzi(:,:)=matmul(sf%dNdgzi,Xmat) a0(:,:)=matmul(sf%dNdgzi,sf%ElemCoord) call inverse_rank_2(a0,a0_inv) call HyperElasticStress(mdl) call HyperElasticDer(mdl,DerType) a1(:,:)=matmul(transpose(sf%dNdgzi),a0) a2(:,:)=matmul(a1,mdl%sigma) a3(:,:)=matmul(a1,mdl%F_ij) a4(:,:,:,:)=0.0d0 Jgzimat(:,:)=matmul(mdl%F_iJ,dXdgzi) detJgzi=det_mat(Jgzimat,size(Jgzimat,1) ) dXdgzi(:,:)=matmul(sf%dNdgzi,Xmat) call inverse_rank_2(dXdgzi,dgzidX) dNdX(:,:) =matmul(transpose(sf%dNdgzi),dgzidX ) a1(:,:)=matmul(dNdX,mdl%F_inv_iJ) do I=1,elemnod_num do ii=1,3 do j=1,3 loc_1=(I-1)*3+ii gvec_e(loc_1)=gvec_e(loc_1)& +a1(I,j)*mdl%sigma(ii,j)*detJgzi enddo enddo enddo end subroutine !================================================================================= !================================================================================= ! K�}�g���N�X�ւ̏d�ˍ��킹 subroutine K_mat_ICU(Kmat, elem_nod, i, Kemat) integer(int32), intent(in) :: i, elem_nod(:,:) real(real64), intent(in) :: Kemat(:,:) real(real64), intent(inout) :: Kmat(:,:,:) Kmat(i,:,:)=Kmat(i,:,:) + Kemat(:,:) end subroutine !================================================================================= !================================================================================= subroutine g_vector_ICU(elem,elem_nod,gvec_e,gvec) integer(int32), intent(in) :: elem_nod(:,:),elem real(real64),intent(in):: gvec_e(:) real(real64),intent(inout)::gvec(:,:) gvec(elem,:)=gvec_e(:) end subroutine !================================================================================= !================================================================================= subroutine F_tensor_ICU(obj,elem,gauss,F_iJ_n,F_iJ) class(FiniteDeform_),intent(inout)::obj integer(int32), intent(in)::elem,gauss real(real64),allocatable:: F_iJ_n(:,:),F_iJ(:,:),dNdx(:,:),x_dNdxi_inv(:,:),x_dNdxi(:,:),& dumat_t(:,:),f_n_n_1(:,:),dumat(:,:),x_u(:,:),tr_dNdgzi(:,:),dF_iJ(:,:) real(real64)::LungeKutta(3,3),LungeKutta_1(3,3),LungeKutta_2(3,3),LungeKutta_3(3,3),& LungeKutta_4(3,3),delta(3,3),dt integer(int32) n,m,i,j n=size(obj%FEMDomain%ShapeFunction%ElemCoord,1) m=size(obj%FEMDomain%ShapeFunction%ElemCoord,2) if(.not.allocated( F_iJ_n )) allocate( F_iJ_n(3,3) ) if(.not.allocated( f_n_n_1 )) allocate( f_n_n_1(3,3) ) if(.not.allocated( F_iJ )) allocate( F_iJ(3,3) ) if(.not.allocated( dF_iJ )) then allocate( dF_iJ(3,3) ) dF_iJ(:,:)=0.0d0 endif if(.not.allocated( dNdx )) allocate( dNdx(n,m) ) if(.not.allocated( x_dNdxi )) allocate( x_dNdxi(m,m) ) if(.not.allocated( x_dNdxi_inv )) allocate( x_dNdxi_inv(m,m) ) if(.not.allocated( dumat_t )) allocate( dumat_t(m,n) ) if(.not.allocated( dumat )) allocate( dumat(n,m) ) if(.not.allocated( x_u )) allocate( x_u(n,m) ) if(.not.allocated( tr_dNdgzi )) allocate( tr_dNdgzi(n,m) ) n=size(obj%FEMDomain%ShapeFunction%ElemCoord,1) m=size(obj%FEMDomain%ShapeFunction%ElemCoord,2) do i=1,n !num of node per elems do j=1,m ! num of dim dumat(i,j)=obj%DeformVecEBEInc(elem, m*(i-1) + j ) enddo enddo dumat_t(:,:)=transpose(dumat) delta(:,:)=0.0d0 delta(1,1)=1.0d0 delta(2,2)=1.0d0 delta(3,3)=1.0d0 x_dNdxi(:,:)=0.0d0 tr_dNdgzi(:,:)=transpose(obj%FEMDomain%ShapeFunction%dNdgzi) x_u(:,:)=obj%FEMDomain%ShapeFunction%ElemCoord(:,:) -dumat(:,:) x_dNdxi(:,:)=transpose(matmul(obj%FEMDomain%ShapeFunction%dNdgzi,x_u)) !dxndgzi at tn call inverse_rank_2(x_dNdxi,x_dNdxi_inv)!dgzidxn at tn dNdx(:,:)=matmul(tr_dNdgzi,x_dNdxi_inv) if(m==2)then f_n_n_1(:,:)=0.0d0 f_n_n_1(3,3)=1.0d0 F_iJ(:,:)=0.0d0 F_iJ(3,3)=1.0d0 F_iJ(1,1)=1.0d0 F_iJ(1,2)=0.0d0 F_iJ(2,1)=0.0d0 F_iJ(2,2)=1.0d0 F_iJ_n(:,:)=0.0d0 F_iJ_n(1,1)=obj%DeformStrain(elem,gauss,7) F_iJ_n(1,2)=obj%DeformStrain(elem,gauss,9) F_iJ_n(1,3)=0.0d0 F_iJ_n(2,1)=obj%DeformStrain(elem,gauss,10) F_iJ_n(2,2)=obj%DeformStrain(elem,gauss,8) F_iJ_n(2,3)=0.0d0 F_iJ_n(3,1)=0.0d0 F_iJ_n(3,2)=0.0d0 F_iJ_n(3,3)=1.0d0 f_n_n_1(1:2,1:2)=F_iJ(1:2,1:2) + matmul(transpose(dumat),dNdx) F_iJ(:,:)=matmul(f_n_n_1,F_iJ_n) obj%DeformStrain(elem,gauss,11)=F_iJ(1,1) obj%DeformStrain(elem,gauss,12)=F_iJ(2,2) obj%DeformStrain(elem,gauss,13)=F_iJ(1,2) obj%DeformStrain(elem,gauss,14)=F_iJ(2,1) elseif(m==3)then f_n_n_1(:,:)=0.0d0 F_iJ(:,:)=0.0d0 F_iJ(1,1)=1.0d0 F_iJ(2,2)=1.0d0 F_iJ(3,3)=1.0d0 F_iJ_n(:,:)=0.0d0 F_iJ_n(1,1)=obj%DeformStrain(elem,gauss,7) F_iJ_n(1,2)=obj%DeformStrain(elem,gauss,9) F_iJ_n(1,3)=obj%DeformStrain(elem,gauss,23) F_iJ_n(2,1)=obj%DeformStrain(elem,gauss,10) F_iJ_n(2,2)=obj%DeformStrain(elem,gauss,8) F_iJ_n(2,3)=obj%DeformStrain(elem,gauss,24) F_iJ_n(3,1)=obj%DeformStrain(elem,gauss,25) F_iJ_n(3,2)=obj%DeformStrain(elem,gauss,26) F_iJ_n(3,3)=obj%DeformStrain(elem,gauss,22) dF_iJ(1,1)=obj%DeformStrain(elem,gauss,32) dF_iJ(1,2)=obj%DeformStrain(elem,gauss,33) dF_iJ(1,3)=obj%DeformStrain(elem,gauss,34) dF_iJ(2,1)=obj%DeformStrain(elem,gauss,35) dF_iJ(2,2)=obj%DeformStrain(elem,gauss,36) dF_iJ(2,3)=obj%DeformStrain(elem,gauss,37) dF_iJ(3,1)=obj%DeformStrain(elem,gauss,38) dF_iJ(3,2)=obj%DeformStrain(elem,gauss,39) dF_iJ(3,3)=obj%DeformStrain(elem,gauss,40) ! Clank-Nicolson f_n_n_1(:,:)=F_iJ(:,:) + 0.50d0*matmul(transpose(dumat),dNdx) + 0.50d0*dF_iJ(:,:) ! 4th order Lunge-Kutta !dt=1.0d0 !f_n_n_1(:,:)=matmul(transpose(dumat),dNdx) !LungeKutta_1(:,:)=f_n_n_1(:,:) !LungeKutta_2(:,:)=f_n_n_1(:,:)+dt/2.0d0*matmul(f_n_n_1,LungeKutta_1) !LungeKutta_3(:,:)=f_n_n_1(:,:)+dt/2.0d0*matmul(f_n_n_1,LungeKutta_2) !LungeKutta_4(:,:)=f_n_n_1(:,:)+dt*matmul(f_n_n_1,LungeKutta_3) !LungeKutta(:,:)=dt/6.0d0*( LungeKutta_1(:,:)+2.0d0*LungeKutta_2(:,:)+2.0d0*LungeKutta_3(:,:)+LungeKutta_4(:,:) ) !F_iJ(:,:)=F_iJ_n(:,:) + matmul(LungeKutta,F_iJ_n) F_iJ(:,:)=matmul(f_n_n_1,F_iJ_n) obj%DeformStrain(elem,gauss,11)= F_iJ(1,1) obj%DeformStrain(elem,gauss,12)= F_iJ(2,2) obj%DeformStrain(elem,gauss,13)= F_iJ(1,2) obj%DeformStrain(elem,gauss,14)= F_iJ(2,1) obj%DeformStrain(elem,gauss,27)= F_iJ(3,3) obj%DeformStrain(elem,gauss,28)= F_iJ(1,3) obj%DeformStrain(elem,gauss,29)= F_iJ(2,3) obj%DeformStrain(elem,gauss,30)= F_iJ(3,1) obj%DeformStrain(elem,gauss,31)= F_iJ(3,2) dF_iJ(:,:)=matmul(transpose(dumat),dNdx) obj%DeformStrain(elem,gauss,41) = dF_iJ(1,1) obj%DeformStrain(elem,gauss,42) = dF_iJ(2,2) obj%DeformStrain(elem,gauss,43) = dF_iJ(1,2) obj%DeformStrain(elem,gauss,44) = dF_iJ(2,1) obj%DeformStrain(elem,gauss,45) = dF_iJ(3,3) obj%DeformStrain(elem,gauss,46) = dF_iJ(1,3) obj%DeformStrain(elem,gauss,47) = dF_iJ(2,3) obj%DeformStrain(elem,gauss,48) = dF_iJ(3,1) obj%DeformStrain(elem,gauss,49) = dF_iJ(3,2) else print*, "FiniteDeformationClass/F_tensor_ICU :: >> Dimension",m,"is not Supported." stop endif end subroutine !================================================================================= !================================================================================= subroutine C_tensor(F,C_IJ,b_ij,itr,dim_num) real(real64),allocatable::C_IJ(:,:),b_ij(:,:),F_T(:,:) real(real64),intent(in)::F(:,:) integer(int32),intent(in)::itr,dim_num integer(int32) i if(.not.allocated( C_IJ )) allocate( C_IJ(3,3) ) if(.not.allocated( b_ij )) allocate( b_ij(3,3) ) call trans_rank_2(F,F_T) if(dim_num==2)then C_IJ(:,:)=0.0d0 b_ij(:,:)=0.0d0 C_IJ(3,3)=1.0d0 b_ij(3,3)=1.0d0 C_IJ(1:2,1:2)=matmul(F_T(1:2,1:2) ,F(1:2,1:2) ) b_ij(1:2,1:2)=matmul(F(1:2,1:2) ,F_T(1:2,1:2) ) elseif(dim_num==3)then C_IJ(:,:)=matmul(F_T ,F ) b_ij(:,:)=matmul(F ,transpose(F) ) else print *, "C_tensor :: dimension : ",dim_num,"is not supported" stop endif end subroutine C_tensor !================================================================================= subroutine Cp_tensor(elem,gauss,strain_measure,Cp_IJ_n,Cp_IJ,Cp_IJ_inv,dim_num) real(real64),allocatable:: Cp_iJ_n(:,:),Cp_iJ(:,:),Cp_IJ_inv(:,:) real(real64),intent(in)::strain_measure(:,:,:) integer(int32), intent(in)::elem,gauss,dim_num integer(int32) i if(.not.allocated( Cp_iJ_n )) allocate( Cp_iJ_n(3,3) ) if(.not.allocated( Cp_iJ )) allocate( Cp_iJ(3,3) ) Cp_iJ(:,:) =0.0d0 Cp_iJ_n(:,:)=0.0d0 Cp_iJ(3,3) =1.0d0 Cp_iJ_n(3,3)=1.0d0 Cp_iJ(1,1)=strain_measure(elem,gauss,4) Cp_iJ(1,2)=strain_measure(elem,gauss,6) Cp_iJ(2,1)=strain_measure(elem,gauss,6) Cp_iJ(2,2)=strain_measure(elem,gauss,5) Cp_iJ_n(1,1)=strain_measure(elem,gauss,1) Cp_iJ_n(1,2)=strain_measure(elem,gauss,3) Cp_iJ_n(2,1)=strain_measure(elem,gauss,3) Cp_iJ_n(2,2)=strain_measure(elem,gauss,2) if(dim_num==3)then Cp_iJ_n(3,3) = strain_measure(elem,gauss,16) Cp_iJ_n(1,3) = strain_measure(elem,gauss,17) Cp_iJ_n(2,3) = strain_measure(elem,gauss,18) Cp_iJ_n(3,1) = strain_measure(elem,gauss,17) Cp_iJ_n(3,2) = strain_measure(elem,gauss,18) Cp_iJ(3,3) = strain_measure(elem,gauss,19) Cp_iJ(1,3) = strain_measure(elem,gauss,20) Cp_iJ(2,3) = strain_measure(elem,gauss,21) Cp_iJ(3,1) = strain_measure(elem,gauss,20) Cp_iJ(3,2) = strain_measure(elem,gauss,21) endif call inverse_rank_2(Cp_IJ,Cp_IJ_inv) end subroutine Cp_tensor !================================================================================= subroutine M_neo_Hookean(C_IJ,Cp_IJ,Cp_IJ_inv,M_IJ,Lamda,mu,elem,gauss) real(real64),intent(in)::C_IJ(:,:),Lamda,mu,Cp_IJ(:,:),Cp_IJ_inv(:,:) real(real64),allocatable::M_IJ(:,:),G_IJ(:,:),C_Cp_1(:,:) integer(int32), intent(in):: elem,gauss integer(int32) i,j,n if(.not. allocated(M_IJ) )allocate(M_IJ(3,3)) allocate(C_Cp_1(3,3),G_IJ(3,3)) !3-D for constitutive equations n=size(C_IJ,1) !Metric tensor G_IJ(:,:)=0.0d0 G_IJ(1,1) =1.0d0 G_IJ(2,2) =1.0d0 G_IJ(3,3) =1.0d0 C_Cp_1(:,:)=matmul(C_IJ,Cp_IJ_inv) M_IJ(:,:)=Lamda/2.0d0*( det_mat(C_IJ,n)/det_mat(Cp_IJ,n)-1.0d0)*G_IJ(:,:) + mu*(C_Cp_1(:,:)-G_IJ(:,:)) end subroutine M_neo_Hookean !================================================================================= subroutine Return_Mapping_MCDP(dim_num,elem,gauss,C_IJ,Cp_IJ,Cp_IJ_n,Cp_IJ_inv,M_IJ,MatPara,& itr_rm,tol,sigma,F_T,F_T_inv,itr,itr_contact,strain_measure,step) real(real64),intent(in)::C_IJ(:,:),Cp_IJ_n(:,:),F_T(:,:),F_T_inv(:,:),MatPara(:) real(real64),intent(inout)::Cp_IJ(:,:),sigma(:,:,:),strain_measure(:,:,:) real(real64),allocatable,intent(inout)::M_IJ(:,:),Cp_IJ_inv(:,:) integer(int32), intent (in)::elem,gauss,itr,itr_rm,step,itr_contact,dim_num real(real64),allocatable::G_IJ(:,:),Jmat(:,:),Xvec(:),Yvec(:),dXvec(:),Zmat(:,:),sigm(:,:),M_FT(:,:),C_IJ_inv(:,:),& M_1(:,:),M_2(:,:),M_2_T(:,:),M_3(:,:),M_4(:,:),M_5(:,:),C_1(:,:),C_2(:,:),C_3(:,:),B_6(:,:),B_7(:,:),B_8(:,:),B_9(:,:),B_10(:,:) real(real64) detF,I1_M,J2_M,J3_M,Theta_M,BI_MC,BI_DP,fc_MC,er,tol,residual_0,yield_function_mc,detC,detCp,alpha,beta,gunma,omega,& a11,a12,a13,a14,a15,a16,a17,a18,a21,a22,a23,a24,a25,a26,MM,xx,c1,c2,c3,c4,c5,c6,c7,val,E,v,Lamda,mu,c,phy,psy integer(int32) n,A,B,I,J,K,L,R,M,P,Q,nn,itrmax,retmap_itr,mesh allocate(G_IJ(3,3),Jmat(4,4),Xvec(4),dXvec(4),Yvec(4),Zmat(3,3),sigm(3,3),M_FT(3,3) ) allocate(M_1(3,3),M_2(3,3),M_2_T(3,3),M_3(3,3),M_4(3,3),M_5(3,3),C_1(3,3),C_2(3,3),C_3(3,3),& B_6(3,3),B_7(3,3),B_8(3,3),B_9(3,3),B_10(3,3)) E=MatPara(1) v=MatPara(2) Lamda=MatPara(3) mu=MatPara(4) c =MatPara(5) phy=MatPara(6) psy=MatPara(7) !compute vriables for R-M algorithm G_IJ(:,:)=0.0d0 G_IJ(1,1)=1.0d0 G_IJ(2,2)=1.0d0 G_IJ(3,3)=1.0d0 n =size(C_IJ,1) detF =sqrt(det_mat(C_IJ,n)) I1_M =M_iJ(1,1) + M_iJ(2,2) + M_iJ(3,3) J2_M =0.0d0 do i=1, 3 do j=1,3 J2_M=J2_M + 1.0d0/2.0d0*( M_IJ(i,j)- I1_M/3.0d0*G_IJ(i,j) )*( M_IJ(i,j)- I1_M/3.0d0*G_IJ(i,j) ) enddo enddo J3_M =0.0d0 do I=1, 3 do J=1,3 do K=1,3 J3_M=J3_M + 1.0d0/3.0d0*( M_IJ(i,j)- I1_M/3.0d0*G_IJ(i,j) )*( M_IJ(j,k)- I1_M/3.0d0*G_IJ(j,k) )& *( M_IJ(k,i)- I1_M/3.0d0*G_IJ(k,i) ) enddo enddo enddo if(J2_M==0.0d0)then Theta_M=0.0d0 else val=1.50d0*sqrt(3.0d0)*J3_M/(J2_M*sqrt(J2_M) ) if(val<-1.0d0)then val=-1.0d0 elseif(val>1.0d0)then val=1.0d0 endif Theta_M=1.0d0/3.0d0*asin(val) endif BI_MC=( 1.0d0 + sin(phy) )/( 1.0d0-sin(phy) ) BI_DP=( 1.0d0 + sin(psy) )/( 1.0d0-sin(psy) ) fc_MC=(2.0d0*c*cos(phy))/( 1.0d0-sin(phy) ) if( ( 1.0d0-sin(phy) ) ==0.0d0) stop " ( 1.0d0-sin(phy) ) =0" Yield_function_MC=1.0d0/detF*( sqrt(J2_M) + ( (BI_MC-1.0d0)*I1_M-3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) ) !write(54,*)& ! sqrt(2.0d0/3.0d0)*sqrt(2.0d0*J2_M)*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& ! sqrt(2.0d0/3.0d0)*sqrt(2.0d0*J2_M)*sin(Theta_M ) + I1_M/3.0d0,& ! sqrt(2.0d0/3.0d0)*sqrt(2.0d0*J2_M)*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0 !Plastic corrector algorithm !print *, sqrt(J2_M),Yield_function_MC,(BI_MC-1.0d0),I1_M,3.0d0*fc_MC !print *, Yield_function_MC !write(109,*) step,itr,Yield_function_MC,detF,I1_M,fc_Mc if(Yield_function_MC>0.0d0 .and. itr_contact >= 3 )then !initialization retmap_itr=1 Xvec(1)=Cp_IJ_n(1,1) Xvec(2)=Cp_IJ_n(2,2) Xvec(3)=Cp_IJ_n(1,2) Xvec(4)=0.0d0 Yvec(1)=0.0d0 Yvec(2)=0.0d0 Yvec(3)=0.0d0 Yvec(4)=yield_function_mc call inverse_rank_2(C_IJ,C_IJ_inv) !Loop for Return-Mapping algorithm do dXvec(:) =0.0d0 Zmat(:,:)=0.0d0 Jmat(:,:)=0.0d0 I1_M =M_iJ(1,1) + M_iJ(2,2) + M_iJ(3,3) J2_M =0.0d0 do i=1, 3 do j=1,3 J2_M=J2_M + 1.0d0/2.0d0*( M_IJ(i,j)- I1_M/3.0d0*G_IJ(i,j) )*( M_IJ(i,j)- I1_M/3.0d0*G_IJ(i,j) ) enddo enddo J3_M =0.0d0 do I=1, 3 do J=1,3 do K=1,3 J3_M=J3_M + 1.0d0/3.0d0*( M_IJ(i,j)- I1_M/3.0d0*G_IJ(i,j) )*( M_IJ(j,k)- I1_M/3.0d0*G_IJ(j,k) )& *( M_IJ(k,i)- I1_M/3.0d0*G_IJ(k,i) ) enddo enddo enddo MM=0.0d0 do i=1,3 do j=1, 3 MM=MM + M_IJ(I,J)*M_IJ(J,I) enddo enddo if(J2_M==0.0d0)then Theta_M=0.0d0 else Theta_M=1.0d0/3.0d0*asin(1.50d0*sqrt(3.0d0)*J3_M/(J2_M*sqrt(J2_M) )) if(J2_M==0.0d0) stop "J2_M=0" endif Yield_function_MC=1.0d0/detF*( sqrt(J2_M) + ( (BI_MC-1.0d0)*I1_M-3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) ) Yvec(4)=yield_function_mc if( detF ==0.0d0) stop " detF =0" !debug print *, "YC", yield_function_mc if(retmap_itr==1)then residual_0=sqrt(dot_product(Yvec,Yvec)) endif print *, "itr, residual R-M",retmap_itr,sqrt(dot_product(Yvec,Yvec)) detC=det_mat(C_IJ,n) detCp=det_mat(Cp_IJ,n) !update dYij/dCpkl a18=-Lamda/2*detC/detCp c1 = -4.0d0*(1.0d0/J2_M*(MM-I1_M*I1_M/3.0d0) + ((BI_DP-1.0d0)/(BI_DP + 2.0d0))**2.0d0 )**(-3.0d0/2.0d0) c2 = 2.0d0*(1.0d0/J2_M*(MM-I1_M*I1_M/3.0d0) + ((BI_DP-1.0d0)/(BI_DP + 2.0d0))**2.0d0 )**(-1.0d0/2.0d0) c3 = -I1_M*I1_M/3.0d0*a18*xvec(4)*c1/J2_M/sqrt(J2_M) c4 = -mu/2.0d0*xvec(4)*c1/J2_M/sqrt(J2_M) c5 = xvec(4)*c1/sqrt(3.0d0)*(BI_DP-1.0d0)/(BI_DP + 2.0d0)/J2_M c6 = c3 + a18*c5*(3.0d0-I1_M) M_1(:,:)=M_IJ(:,:)-I1_M/3.0d0*G_IJ(:,:) !M_2(:,:)=matmul(matmul(C_IJ,M_IJ),Cp_IJ_inv) !M_3(:,:)=matmul(Cp_IJ_inv,M_2) M_4(:,:)=matmul(M_1,Cp_IJ_n) M_5(:,:)=matmul(matmul(Cp_IJ_inv,C_IJ ),matmul(M_1,Cp_IJ_inv ) ) C_1(:,:)=matmul(C_IJ,Cp_IJ_inv ) C_2(:,:)=matmul(Cp_IJ_inv,Cp_IJ_n ) C_3(:,:)=matmul(Cp_IJ_inv,matmul(C_IJ,Cp_IJ_inv) ) alpha=-0.50d0/J2_M/sqrt(J2_M) beta=1.0d0/sqrt(J2_M) gunma=-xvec(4)*mu*c2 c7=c4 + alpha*gunma !in case of the edge if((3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) )==0.0d0)then stop "singular stress of Mohr-Coulomb" endif a11=0.50d0/detF/sqrt(J2_M) a12=-(BI_MC-1.0d0)/detF/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) a13=-1.0d0/detF*( (BI_MC-1.0d0)*I1_M + 3.0d0*detF*fc_MC )/& (3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) )**2.0d0& *(3.0d0*(BI_MC + 1.0d0)*sin(Theta_M)-sqrt(3.0d0)*(BI_MC-1.0d0)*cos(Theta_M) ) a14=sqrt(3.0d0)/2.0d0/sqrt(J2_M)/J2_M/& sqrt(1.0d0-( 27.0d0/4.0d0*J3_M*J3_M/J2_M/J2_M/J2_M )**2.0d0 ) a15=-3.0d0*sqrt(3.0d0)/4.0d0/sqrt(J2_M)/J2_M/J2_M*J3_M/& sqrt(1.0d0-( 27.0d0/4.0d0*J3_M*J3_M/J2_M/J2_M/J2_M )**2.0d0 ) a16=-2.0d0*I1_M/3.0d0 a17=2.0d0*I1_M*I1_M/9.0d0 a18=-Lamda/2*detC/detCp if( (3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) ==0.0d0)& stop " (3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) =0" a21=-a13*a14 a22=a11-a13*a15-a13*a14*a16 a23=-I1_M/3.0d0*a11-a12 + a13*a15*I1_M/3.0d0 + a13*a14/3.0d0*MM-a13*a14*a17 a25=a18*(a21*MM + a22*I1_M + 3.0d0*a23) B_6(:,:)=matmul(C_IJ,matmul(M_IJ,Cp_IJ_inv )) B_7(:,:)=matmul(Cp_IJ_inv,B_6) B_8(:,:)=matmul(Cp_IJ_inv,matmul(C_IJ,Cp_IJ_inv) ) B_6(:,:)=matmul(Cp_IJ_inv,matmul(M_IJ,C_IJ )) B_9(:,:)=matmul(B_6,Cp_IJ_inv) B_6(:,:)=matmul(Cp_IJ_inv,matmul(M_IJ,M_IJ )) B_10(:,:)=matmul(matmul(B_6,C_IJ),Cp_IJ_inv ) Jmat(:,:)=0.0d0 I=1 J=1 K=1 L=1 Jmat(1,1)=1.0d0-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=1 J=1 K=2 L=2 Jmat(1,2)=-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=1 J=1 K=1 L=2 Jmat(1,3)=-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=2 J=2 K=1 L=1 Jmat(2,1)=-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=2 J=2 K=2 L=2 Jmat(2,2)=1.0d0-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=2 J=2 K=1 L=2 Jmat(2,3)=-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=1 J=2 K=1 L=1 Jmat(3,1)=-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=1 J=2 K=2 L=2 Jmat(3,2)=-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) I=1 J=2 K=1 L=2 Jmat(3,3)=1.0d0-c6*Cp_IJ_n(I,J)*Cp_IJ_inv(L,K)-c7*M_4(I,J)*M_5(K,L) + mu*c5*Cp_IJ_n(I,J)*M_5(K,L)& -beta*gunma*C_1(I,K)*C_2(L,J) + 1.0d0/3.0d0*beta*gunma*C_3(L,K)*Cp_IJ_n(I,J) do i=1,size(zmat,1) do j=1, size(zmat,1) zmat(i,j)=zmat(i,j) + (1.0d0/sqrt(J2_M)*( M_IJ(i,j)-I1_M/3.0d0*G_IJ(i,j) ) + & 2.0d0*(BI_DP-1.0d0)/sqrt(3.0d0)/(BI_DP + 2.0d0)*G_IJ(i,j) ) if( (BI_DP + 2.0d0) ==0.0d0) stop " (BI_DP + 2.0d0) =0" enddo enddo zmat(:,:)=2.0d0*( (MM-I1_M*I1_M/3.0d0)/J2_M + & ((BI_DP-1.0d0)/(BI_DP + 2.0d0))**2.0d0 )**(-1.0d0/2.0d0)*zmat(:,:) do M=1,3 Jmat(1,4)=Jmat(1,4)-Zmat(1,M)*Cp_IJ_n(M,1) Jmat(2,4)=Jmat(2,4)-Zmat(2,M)*Cp_IJ_n(M,2) Jmat(3,4)=Jmat(3,4)-Zmat(1,M)*Cp_IJ_n(M,2) enddo I=1 J=1 Jmat(4,1)=a25*Cp_IJ_inv(J,I)-mu*(a21*B_10(J,I) + a22*B_7(I,J) + a23*B_8(J,I)) I=2 J=2 Jmat(4,2)=a25*Cp_IJ_inv(J,I)-mu*(a21*B_10(J,I) + a22*B_7(I,J) + a23*B_8(J,I)) I=1 J=2 Jmat(4,3)=a25*Cp_IJ_inv(J,I)-mu*(a21*B_10(J,I) + a22*B_7(I,J) + a23*B_8(J,I)) Jmat(4,4)=0.0d0 !control parameters !write(53,*)"I,yield_function_mc,norm",retmap_itr,yield_function_mc,sqrt(dot_product(Yvec,Yvec)) ! stop "debug r-m" itrmax=itr_rm er=1.0e-14 nn=size(dXvec) print *, "detJmat=",det_mat(Jmat,4),det_mat(zmat,2) !do i=1,size(Jmat,1) ! write(52,*)Jmat(i,:) !enddo Yvec(1)=Cp_IJ(1,1) Yvec(2)=Cp_IJ(2,2) Yvec(3)=Cp_IJ(1,2) Yvec(4)=yield_function_mc do K=1,3 Yvec(1)=Yvec(1)-(G_IJ(1,K) + Xvec(4)*Zmat(1,K) )*Cp_IJ_n(K,1) Yvec(2)=Yvec(2)-(G_IJ(2,K) + Xvec(4)*Zmat(2,K) )*Cp_IJ_n(K,2) Yvec(3)=Yvec(3)-(G_IJ(1,K) + Xvec(4)*Zmat(1,K) )*Cp_IJ_n(K,2) enddo !solve call gauss_jordan_pv(Jmat, dXvec, Yvec, nn) Xvec(:)=Xvec(:)-dXvec(:) !update variables Cp_IJ(:,:)=0.0d0 Cp_IJ(3,3)=1.0d0 Cp_IJ(1,1)=Xvec(1) Cp_IJ(2,2)=Xvec(2) Cp_IJ(1,2)=Xvec(3) Cp_IJ(2,1)=Xvec(3) call inverse_rank_2(Cp_IJ,Cp_IJ_inv) call M_neo_Hookean(C_IJ,Cp_IJ,Cp_IJ_inv,M_IJ,Lamda,mu,elem,gauss) !Judgement of convergence if( sqrt(dot_product(Yvec,Yvec)) <=TOL)then !converged if( residual_0 ==0.0d0) stop " residual_0 =0" print *, "R-M converge itr=",itr exit elseif(retmap_itr>=itrmax)then print *, "Return-Mapping (E-P) did not converge" !debug: print yield surface xx=0.40d0 mesh=20 do i=1, mesh do j=1, 12 !1 I1_M=xx*(1.0d0-2.0d0/dble(mesh)*dble(i)) Theta_M=3.14159d0/6.0d0-3.14159d0/3.0d0/12.0d0*dble(j) J2_M=( -(BI_MC-1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) write(55,*)& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M ) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0 enddo do j=1, 12 !6 I1_M=xx*(1.0d0-2.0d0/dble(mesh)*dble(i)) Theta_M=-3.14159d0/6.0d0 + 3.14159d0/3.0d0/12.0d0*dble(j) J2_M=( -(BI_MC-1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) write(55,*)& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M ) + I1_M/3.0d0 enddo do j=1, 12 !3 I1_M=xx*(1.0d0-2.0d0/dble(mesh)*dble(i)) Theta_M=3.14159d0/6.0d0-3.14159d0/3.0d0/12.0d0*dble(j) J2_M=( -(BI_MC-1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) write(55,*)& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M ) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0 enddo do j=1, 12 !5 I1_M=xx*(1.0d0-2.0d0/dble(mesh)*dble(i)) Theta_M=-3.14159d0/6.0d0 + 3.14159d0/3.0d0/12.0d0*dble(j) J2_M=( -(BI_MC-1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) write(55,*)& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M ) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0 enddo do j=1, 12 !4 I1_M=xx*(1.0d0-2.0d0/dble(mesh)*dble(i)) Theta_M=3.14159d0/6.0d0-3.14159d0/3.0d0/12.0d0*dble(j) J2_M=( -(BI_MC-1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) write(55,*)& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M ) + I1_M/3.0d0 enddo do j=1, 12 !2 I1_M=xx*(1.0d0-2.0d0/dble(mesh)*dble(i)) Theta_M=-3.14159d0/6.0d0 + 3.14159d0/3.0d0/12.0d0*dble(j) J2_M=( -(BI_MC-1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + & sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) ) write(55,*)& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M ) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,& 2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0 enddo enddo stop else retmap_itr=retmap_itr + 1 cycle endif enddo else !print *, "elastic" endif M_FT(:,:)=matmul(M_IJ,F_T) sigm(:,:)=1.0d0/detF*matmul( F_T_inv , M_FT ) if(dim_num==2)then sigma(elem,gauss,1)=sigm(1,1) sigma(elem,gauss,2)=sigm(2,2) sigma(elem,gauss,3)=sigm(1,2) sigma(elem,gauss,4)=sigm(3,3) elseif(dim_num==3)then sigma(elem,gauss,1)=sigm(1,1) sigma(elem,gauss,2)=sigm(2,2) sigma(elem,gauss,3)=sigm(3,3) sigma(elem,gauss,4)=sigm(1,2) sigma(elem,gauss,5)=sigm(2,3) sigma(elem,gauss,6)=sigm(3,1) else stop "Return_Mapping_MCDP >> dimension is to be 2 or 3" endif strain_measure(elem,gauss,4)=Cp_IJ(1,1) strain_measure(elem,gauss,6)=Cp_IJ(1,2) strain_measure(elem,gauss,6)=Cp_IJ(2,1) strain_measure(elem,gauss,5)=Cp_IJ(2,2) if(dim_num==3)then strain_measure(elem,gauss,19) = Cp_iJ(3,3) strain_measure(elem,gauss,20) = Cp_iJ(1,3) strain_measure(elem,gauss,21) = Cp_iJ(2,3) endif end subroutine Return_Mapping_MCDP !================================================================================= ! D�}�g���N�X�̌v�Z subroutine Ce_neoHK_current(dim_num, elem, gauss,Lame1,Lame2,C_IJ,Cp_IJ,b_ij,M_IJ,Ce_neoHK,F_T,F_T_inv,ij) integer(int32), intent(in) :: dim_num,elem, gauss real(real64), intent(in) :: Lame1,Lame2,C_IJ(:,:),Cp_IJ(:,:),b_ij(:,:),F_T(:,:),F_T_inv(:,:),M_IJ(:,:) real(real64), allocatable, intent(out) :: Ce_neoHK(:,:) integer(int32),allocatable,intent(out)::ij(:,:) real(real64), allocatable::t_ij(:,:),G_IJ(:,:) integer(int32) n,m,p,q,i,j,k,l real(real64) detF,detCp,detC allocate(G_IJ(3,3)) G_IJ(:,:)=0.0d0 G_IJ(1,1) =1.0d0 G_IJ(2,2) =1.0d0 G_IJ(3,3) =1.0d0 if(allocated(ij) ) deallocate(ij) if(dim_num==2)then n = 3 ! �Ђ��݂���11,��22,��12��3���� allocate(ij(3,1) ) ij(1,1)=1 ; ij(1,2)=1; ij(2,1)=2 ; ij(2,2)=2; ij(3,1)=1 ; ij(3,2)=2; m=size(C_IJ,1) if(.not.allocated(Ce_neoHK) ) allocate(Ce_neoHK(n,n)) allocate (t_ij(n,n)) detC=det_mat(C_IJ,m) detCp=det_mat(Cp_IJ,m) detF=sqrt(detC) t_ij=matmul( F_T_inv ,matmul(M_IJ,F_T) ) do p=1,3 do q=1,3 i=ij(p,1) j=ij(p,2) k=ij(q,1) l=ij(q,2) Ce_neoHK(p,q)=Lame1*(detC/detCp*b_ij(i,j)*b_ij(l,k)& - Lame1*(detC/detCp-1)*b_ij(i,k)*b_ij(l,j) )& + 2.0d0*Lame2*b_ij(i,k)*b_ij(l,j)& + G_IJ(i,k)*t_ij(l,j) enddo enddo Ce_neoHK(:,:)=1.0d0/detF*Ce_neoHK(:,:) elseif(dim_num==3)then n = 6 ! �Ђ��݂���11,��22,��12��3���� allocate(ij(6,2) ) ij(1,1)=1 ; ij(1,2)=1; ij(2,1)=2 ; ij(2,2)=2; ij(3,1)=3 ; ij(3,2)=3; ij(4,1)=1 ; ij(4,2)=2; ij(5,1)=2 ; ij(5,2)=3; ij(6,1)=3 ; ij(6,2)=1; m=size(C_IJ,1) if(.not.allocated(Ce_neoHK) ) allocate(Ce_neoHK(n,n)) allocate (t_ij(n,n)) detC=det_mat(C_IJ,m) detCp=det_mat(Cp_IJ,m) detF=sqrt(detC) t_ij=matmul( F_T_inv ,matmul(M_IJ,F_T) ) do p=1,6 do q=1,6 i=ij(p,1) j=ij(p,2) k=ij(q,1) l=ij(q,2) Ce_neoHK(p,q)=Lame1*(detC/detCp*b_ij(i,j)*b_ij(l,k)& - Lame1*(detC/detCp-1)*b_ij(i,k)*b_ij(l,j) )& + 2.0d0*Lame2*b_ij(i,k)*b_ij(l,j) !+ G_IJ(i,k)*t_ij(l,j) enddo enddo Ce_neoHK(:,:)=1.0d0/detF*Ce_neoHK(:,:) else print *, "Ce_neoHK_current >> Dimension, ",dim_num," is not supported" stop endif end subroutine Ce_neoHK_current !================================================================================= !================================================================================= ! D�}�g���N�X�̌v�Z subroutine GetSigmaVec(Sigma,Sigma_ij,dim_num) real(real64), intent(in) :: Sigma_ij(:,:) real(real64), allocatable, intent(inout) :: Sigma(:) integer(int32),allocatable::ij(:,:) integer(int32),intent(in)::dim_num integer(int32) n,m,p,q,i,j,k,l real(real64) detF,detCp,detC if(allocated(ij) ) deallocate(ij) if(dim_num==2)then n = 3 ! �Ђ��݂���11,��22,��12��3���� allocate(ij(3,1) ) ij(1,1)=1 ; ij(1,2)=1; ij(2,1)=2 ; ij(2,2)=2; ij(3,1)=1 ; ij(3,2)=2; if(.not.allocated(Sigma) ) allocate(Sigma(n)) do p=1,3 i=ij(p,1) j=ij(p,2) Sigma(p)=Sigma_ij(i,j) enddo elseif(dim_num==3)then n = 6 ! �Ђ��݂���11,��22,��12��3���� allocate(ij(6,2) ) ij(1,1)=1 ; ij(1,2)=1; ij(2,1)=2 ; ij(2,2)=2; ij(3,1)=3 ; ij(3,2)=3; ij(4,1)=1 ; ij(4,2)=2; ij(5,1)=2 ; ij(5,2)=3; ij(6,1)=3 ; ij(6,2)=1; if(.not.allocated(Sigma) ) allocate(Sigma(n)) do p=1,6 Sigma(p)=Sigma_ij(ij(p,1),ij(p,2)) enddo else print *, "Dmat >> Dimension, ",dim_num," is not supported" stop endif end subroutine !================================================================================= !================================================================================= ! D�}�g���N�X�̌v�Z subroutine GetDmat(Dmat,c_ijkl,dim_num) real(real64), intent(in) :: c_ijkl(:,:,:,:) real(real64), allocatable, intent(inout) :: Dmat(:,:) integer(int32),allocatable::ij(:,:) integer(int32),intent(in)::dim_num integer(int32) n,m,p,q,i,j,k,l real(real64) detF,detCp,detC if(allocated(ij) ) deallocate(ij) if(dim_num==2)then n = 3 ! �Ђ��݂���11,��22,��12��3���� allocate(ij(3,1) ) ij(1,1)=1 ; ij(1,2)=1; ij(2,1)=2 ; ij(2,2)=2; ij(3,1)=1 ; ij(3,2)=2; if(.not.allocated(Dmat) ) allocate(Dmat(n,n)) do p=1,3 do q=1,3 i=ij(p,1) j=ij(p,2) k=ij(q,1) l=ij(q,2) Dmat(p,q)=c_ijkl(i,j,k,l) enddo enddo elseif(dim_num==3)then n = 6 ! �Ђ��݂���11,��22,��12��3���� allocate(ij(6,2) ) ij(1,1)=1 ; ij(1,2)=1; ij(2,1)=2 ; ij(2,2)=2; ij(3,1)=3 ; ij(3,2)=3; ij(4,1)=1 ; ij(4,2)=2; ij(5,1)=2 ; ij(5,2)=3; ij(6,1)=3 ; ij(6,2)=1; if(.not.allocated(Dmat) ) allocate(Dmat(n,n)) do p=1,6 do q=1,6 i=ij(p,1) j=ij(p,2) k=ij(q,1) l=ij(q,2) Dmat(p,q)=c_ijkl(i,j,k,l) enddo enddo else print *, "Dmat >> Dimension, ",dim_num," is not supported" stop endif end subroutine !================================================================================= !B�}�g���N�X�̌v�Z subroutine B_mat(dim_num,Psymat, Jmat, detJ, Bmat,mm) real(real64), intent(in) :: Psymat(:,:), Jmat(:,:), detJ ! J�̋t�s�� real(real64), allocatable, intent(inout) :: Bmat(:,:) integer(int32),intent(in)::dim_num real(real64), allocatable :: JPsy(:,:), Jin(:,:) integer(int32) k, l,m, n, a, b, p,mm,i,j,q if(dim_num==2)then k=3 elseif(dim_num==3)then k=6 else stop "B_mat >> dim_num = tobe 2 or 3 " endif !k = size(ij,1) ! �Ђ��݂���11,��22,��12��3���� l = mm ! 8�ړ_�v�f*2���� m = size(Psymat, 1) n = size(Psymat, 2) if(allocated(Bmat )) deallocate(Bmat) allocate (Bmat(k, l)) allocate(JPsy(m,n)) allocate(Jin(m,m)) ! J:Psymat�̌v�Z if(mm/2==3)then stop 'Now Constructing' elseif(mm/2==4)then if(detJ==0.0d0) stop "Bmat,detJ=0" Jin(1,1) = (1.0d0 / detJ) * Jmat(2,2) Jin(2,2) = (1.0d0 / detJ) * Jmat(1,1) Jin(1,2) = (-1.0d0 / detJ) * Jmat(1,2) Jin(2,1) = (-1.0d0 / detJ) * Jmat(2,1) JPsy(:,:) = matmul(Jin, Psymat) Bmat(1,1) = JPsy(1,1) Bmat(1,2) = 0.0d0 Bmat(1,3) = 1.0d0 * JPsy(1,2) Bmat(1,4) = 0.0d0 Bmat(1,5) = 1.0d0 * JPsy(1,3) Bmat(1,6) = 0.0d0 Bmat(1,7) = 1.0d0 * JPsy(1,4) Bmat(1,8) = 0.0d0 Bmat(2,1) = 0.0d0 Bmat(2,2) = 1.0d0 * JPsy(2,1) Bmat(2,3) = 0.0d0 Bmat(2,4) = 1.0d0 * JPsy(2,2) Bmat(2,5) = 0.0d0 Bmat(2,6) = 1.0d0 * JPsy(2,3) Bmat(2,7) = 0.0d0 Bmat(2,8) = 1.0d0 * JPsy(2,4) Bmat(3,1) = Bmat(2,2) Bmat(3,2) = Bmat(1,1) Bmat(3,3) = Bmat(2,4) Bmat(3,4) = Bmat(1,3) Bmat(3,5) = Bmat(2,6) Bmat(3,6) = Bmat(1,5) Bmat(3,7) = Bmat(2,8) Bmat(3,8) = Bmat(1,7) elseif(mm/2==8)then Jin(1,1) = (1.0d0 / detJ) * Jmat(2,2) Jin(2,2) = (1.0d0 / detJ) * Jmat(1,1) Jin(1,2) = (-1.0d0 / detJ) * Jmat(2,1) Jin(2,1) = (-1.0d0 / detJ) * Jmat(1,2) JPsy(:,:) = matmul(Jin, Psymat) Bmat(1,1) = -JPsy(1,1) Bmat(1,2) = 0.0d0 Bmat(1,3) = -1.0d0 * JPsy(1,2) Bmat(1,4) = 0.0d0 Bmat(1,5) = -1.0d0 * JPsy(1,3) Bmat(1,6) = 0.0d0 Bmat(1,7) = -1.0d0 * JPsy(1,4) Bmat(1,8) = 0.0d0 Bmat(1,9) = -1.0d0 * JPsy(1,5) Bmat(1,10) = 0.0d0 Bmat(1,11) = -1.0d0 * JPsy(1,6) Bmat(1,12) = 0.0d0 Bmat(1,13) = -1.0d0 * JPsy(1,7) Bmat(1,14) = 0.0d0 Bmat(1,15) = -1.0d0 * JPsy(1,8) Bmat(1,16) = 0.0d0 Bmat(2,1) = 0.0d0 Bmat(2,2) = -1.0d0 * JPsy(2,1) Bmat(2,3) = 0.0d0 Bmat(2,4) = -1.0d0 * JPsy(2,2) Bmat(2,5) = 0.0d0 Bmat(2,6) = -1.0d0 * JPsy(2,3) Bmat(2,7) = 0.0d0 Bmat(2,8) = -1.0d0 * JPsy(2,4) Bmat(2,9) = 0.0d0 Bmat(2,10) = -1.0d0 * JPsy(2,5) Bmat(2,11) = 0.0d0 Bmat(2,12) = -1.0d0 * JPsy(2,6) Bmat(2,13) = 0.0d0 Bmat(2,14) = -1.0d0 * JPsy(2,7) Bmat(2,15) = 0.0d0 Bmat(2,16) = -1.0d0 * JPsy(2,8) Bmat(3,1) = Bmat(2,2) Bmat(3,2) = Bmat(1,1) Bmat(3,3) = Bmat(2,4) Bmat(3,4) = Bmat(1,3) Bmat(3,5) = Bmat(2,6) Bmat(3,6) = Bmat(1,5) Bmat(3,7) = Bmat(2,8) Bmat(3,8) = Bmat(1,7) Bmat(3,9) = Bmat(2,10) Bmat(3,10) = Bmat(1,9) Bmat(3,11) = Bmat(2,12) Bmat(3,12) = Bmat(1,11) Bmat(3,13) = Bmat(2,14) Bmat(3,14) = Bmat(1,13) Bmat(3,15) = Bmat(2,16) Bmat(3,16) = Bmat(1,15) elseif(k==6 )then if(detJ==0.0d0) stop "Bmat,detJ=0" call inverse_rank_2(Jmat,Jin) JPsy(:,:) = transpose(matmul(transpose(Psymat),Jin)) !dNdgzi* dgzidx Bmat(:,:)=0.0d0 do q=1,size(JPsy,2) do p=1,dim_num Bmat(p,dim_num*(q-1) + p )=JPsy(p,q) enddo Bmat(4,dim_num*(q-1) + 1 )=JPsy(2,q); Bmat(4, dim_num*(q-1) + 2 )=JPsy(1,q);Bmat(4, dim_num*(q-1) + 3 )=0.0d0 ; Bmat(5,dim_num*(q-1) + 1 )=0.0d0 ; Bmat(5, dim_num*(q-1) + 2 )=JPsy(3,q);Bmat(5, dim_num*(q-1) + 3 )=JPsy(2,q); Bmat(6,dim_num*(q-1) + 1 )=JPsy(3,q); Bmat(6, dim_num*(q-1) + 2 )=0.0d0 ;Bmat(6, dim_num*(q-1) + 3 )=JPsy(1,q); enddo Bmat(4:6,:)=0.50d0*Bmat(4:6,:) else stop "Bmat >> The element is not supported." endif end subroutine B_mat !================================================================================= ! �K�E�X�ϕ� subroutine K_mat_e(j,s, BTmat, Ce_neoHK, Bmat, detJ, Kmat_e,F_iJ) integer(int32), intent(in) :: j real(real64), intent(in) :: BTmat(:,:), Ce_neoHK(:,:), Bmat(:,:), detJ, s(:),F_iJ(:,:) real(real64), intent(out) :: Kmat_e(:,:) real(real64), allocatable :: DBmat(:,:) integer(int32) nm, e,n n= size(F_iJ,1) nm = size(Bmat, 2) e = size(Bmat,1) allocate ( DBmat(e, nm)) DBmat = matmul(Ce_neoHK, Bmat) Kmat_e(:,:) = Kmat_e(:,:) + s(j) * detJ *matmul(BTmat, DBmat) !/det_mat(F_iJ,n) end subroutine K_mat_e !================================================================================= subroutine g_vector_e(elem,gauss,s, BTmat,sigma, detJ, gvec_e) integer(int32), intent(in) :: elem,gauss real(real64), intent(in) :: BTmat(:,:), sigma(:,:,:), detJ, s(:) real(real64), intent(inout) :: gvec_e(:) real(real64), allocatable :: sigm(:) integer(int32) nm, i if(size(BTmat,2)==3)then allocate(sigm(3)) nm = size(BTmat, 1) do i=1,size(BTmat,2) sigm(i)=sigma(elem,gauss,i) enddo elseif(size(BTmat,2)==6)then allocate(sigm(6)) nm = size(BTmat, 1) do i=1,size(BTmat,2) sigm(i)=sigma(elem,gauss,i) enddo else stop "Error :: g_vec_e :: Invalid size of Bmat" endif gvec_e(:)=gvec_e(:) + s(gauss)*detJ*matmul(BTmat,sigm) end subroutine g_vector_e !================================================================================= !######################## Solve DiffusionEq ######################## subroutine SolveFiniteDeform(obj,OptionItr,Solvertype,nr_tol) class(FiniteDeform_),intent(inout)::obj integer(int32),optional,intent(in)::OptionItr character(*),optional,intent(in)::Solvertype real(real64) ,optional,intent(in) :: nr_tol character*70 ::solver,defaultsolver type(IO_) :: f real(real64),allocatable::Amat(:,:),bvec(:),xvec(:) real(real64)::val,er integer(int32) ::i,j,n,m,k,l,dim1,dim2,nodeid1,nodeid2,localid,itrmax,SetBC,int1,int2 integer(int32) :: dim_num,node_num,elem_num,node_num_elmtl !obj%nr_tol=1.0e-08 if(present(nr_tol) )then obj%nr_tol = nr_tol endif node_num=size(obj%FEMDomain%Mesh%NodCoord,1) dim_num=size(obj%FEMDomain%Mesh%NodCoord,2) elem_num=size(obj%FEMDomain%Mesh%ElemNod,1) node_num_elmtl=size(obj%FEMDomain%Mesh%ElemNod,2) if(present(OptionItr) )then ! Interation in Netwon's Loop >> Set all BCs for zero. if(OptionItr<=0)then SetBC=1 elseif(OptionItr>=1)then SetBC=0 else print *, "ERROR :: FiniteDeformationClass.f90 >> SolveDeformStress :: Please check OptionalItr" endif else SetBC=1 endif !if sorving initial value if(SetBC==1)then obj%DeformVecEBEInc(:,:)=0.0d0 obj%DeformVecGloInc(:)=0.0d0 endif if(present(SolverType) )then solver=SolverType else solver="BiCGSTAB" endif n=node_num m=dim_num allocate(Amat(n*m,n*m) ,bvec(n*m),xvec(n*m) ) Amat(:,:)=0.0d0 bvec(:) =0.0d0 !print *, "stop debug" !call showArray(Amat, Name="Amat.txt") !stop !=================================== ! assemble matrix do i=1,elem_num do j=1, node_num_elmtl nodeid1=obj%FEMDomain%Mesh%ElemNod(i,j) do k=1,node_num_elmtl nodeid2=obj%FEMDomain%Mesh%ElemNod(i,k) do dim1=1,dim_num do dim2=1,dim_num int1=(dim_num*(j-1)) + dim1 int2=(dim_num*(k-1)) + dim2 Amat(dim_num*(nodeid1-1) + dim1, dim_num*(nodeid2-1) + dim2 )=& Amat(dim_num*(nodeid1-1) + dim1, dim_num*(nodeid2-1) + dim2 ) + & (obj%DeformStressMat(i, int1, int2 )) enddo enddo enddo enddo enddo !=================================== call GetTractionVector(obj) call GetInternalVector(obj) call GetResidualVector(obj) bvec(:)=obj%ResidualVecGlo(:) !=================================== ! introduce D.B.C do k=1,size(obj%FEMDomain%Boundary%DBoundNum) do i=1,size(obj%FEMDomain%Boundary%DBoundNodID,1) nodeid1=obj%FEMDomain%Boundary%DBoundNodID(i,k) if(nodeid1<1)then cycle else val=obj%FEMDomain%Boundary%DBoundValInc(i,k) do j=1,size(bvec) bvec(j)=bvec(j)-Amat(j,dim_num*(nodeid1-1) + k)*val*dble(SetBC) enddo Amat(dim_num*(nodeid1-1) + k,:) = 0.0d0 Amat(:,dim_num*(nodeid1-1) + k) = 0.0d0 Amat(dim_num*(nodeid1-1) + k,dim_num*(nodeid1-1) + k) = 1.0d0 bvec(dim_num*(nodeid1-1) + k) = val*dble(SetBC) endif enddo enddo !=================================== obj%ResidualVecGlo(:)=bvec(:) itrmax=1000 er=1.0e-15 n=size(bvec) xvec(:)=0.0d0 if(trim(solver(1:11) )=="GaussJordan")then print *, "Solver type :: GaussJordan" call gauss_jordan_pv(Amat, xvec, bvec, n) elseif(trim(solver(1:8) )=="BiCGSTAB")then print *, "Solver type :: BiCGSTAB" call bicgstab_dirichlet(Amat, bvec, xvec, size(xvec), itrmax, er,& obj%FEMDomain%Boundary%DBoundNodID,obj%FEMDomain%Boundary%DBoundValInc,SetBC) else print *, "Solver Name :",trim(solver) print *, "Critical ERROR :: No solvers are selected in FiniteDeform_" stop endif ! print *, size(bvec) ! call showArray(Amat,name="Amat.txt") ! call f%open("Bvec.txt") ! do i=1,size(bvec) ! writE(f%fh,*) bvec(:) ! enddo ! call f%close() ! ! stop !===================================== ! Export Values to Element-by-Element form do i=1,elem_num do j=1, node_num_elmtl nodeid1=obj%FEMDomain%Mesh%ElemNod(i,j) do dim1=1,dim_num obj%DeformVecEBEInc(i, dim_num*(j-1) + dim1)=& obj%DeformVecEBEInc(i, dim_num*(j-1) + dim1) + & (xvec( (dim_num*(nodeid1-1)) + dim1 )) enddo enddo enddo obj%DeformVecGloInc(:)=obj%DeformVecGloInc(:) + xvec(:) if(dot_product(obj%DeformVecGloInc,obj%DeformVecGloInc) ==0.0d0)then obj%error= dot_product(xvec,xvec) else obj%error = abs(1.0d0- dot_product(obj%DeformVecGloInc-xvec,obj%DeformVecGloInc-xvec)/& dot_product(obj%DeformVecGloInc,obj%DeformVecGloInc)) endif !===================================== end subroutine !######################## Solve DiffusionEq ######################## !######################## Display Finite Deformation ######################## subroutine resultFiniteDeform(obj,path,name,step) class(FiniteDeform_),intent(inout)::obj character(*),intent(in) :: Name,path integer(int32),optional,intent(in)::step real(real64),allocatable::DBCvec(:,:),DispVector(:,:) integer(int32) :: i,j,n,m,fstep,dim_num if(size(obj%FEMDomain%Mesh%NodCoord,2)==2)then ! 2-D condition if(present(step) )then fstep=step else fstep=1 endif call GmshExportStress(obj%FEMDomain,obj%DeformVecGloTot,obj%DeformStress,& obj%DeformStrain,fstep,Name=trim(path)//"/"//trim(adjustl(name))) call obj%getDispVector(DispVector) call GmshPlotVector(obj%FEMDomain,Vector=DispVector,Name=trim(path)//"/"& //trim(adjustl(name)),FieldName="Displacement",& Step=step,withMsh=.true.,NodeWize=.true.,onlyDirichlet=.true.) elseif(size(obj%FEMDomain%Mesh%NodCoord,2)==3)then ! 3-D condition if(present(step) )then fstep=step else fstep=1 endif call GmshExportStress(obj%FEMDomain,obj%DeformVecGloTot,obj%DeformStress,& obj%DeformStrain,step,Name=trim(path)//"/"//trim(adjustl(name)) ) call obj%getDBCVector(DBCvec) call GmshPlotVector(obj=obj%FEMDomain,Vector=DBCvec,Name=trim(path)//"/"& //trim(adjustl(name)),Step=fstep,FieldName="DispBoundary",& NodeWize=.true. ) !call obj%exportAsPly() endif end subroutine !######################## Display Finite Deformation ######################## subroutine exportAsPlyFiniteDeform(obj) class(FiniteDeform_),intent(inout) :: obj real(real64),allocatable :: scalar(:) integer(int32) :: n,m n=size(obj%FEMDomain%Mesh%NodCoord,1) allocate(scalar(n) ) end subroutine !######################## Display Finite Deformation ######################## subroutine DisplayDeformStress(obj,DisplayMode,OptionalStep,Name,withDirichlet) class(FiniteDeform_),intent(inout)::obj character(*),optional,intent(in) :: Name,DisplayMode logical,optional,intent(in) :: withDirichlet integer(int32),optional,intent(in)::OptionalStep real(real64),allocatable::DBCvec(:,:),DispVector(:,:) integer(int32) :: i,j,n,m,step,dim_num if(.not.associated(obj%FEMDomain) )then return endif if(obj%FEMDomain%Mesh%empty() .eqv. .true. )then return endif if(size(obj%FEMDomain%Mesh%NodCoord,2)==2)then ! 2-D condition if(present(OptionalStep) )then step=OptionalStep else step=1 endif if(present(DisplayMode) )then if(trim(DisplayMode)=="Terminal" .or. trim(DisplayMode)=="terminal")then do i=1,size(obj%DeformVecEBEInc,1) print *, obj%DeformVecEBEInc(i,:) enddo elseif(trim(DisplayMode)=="gmsh" .or. trim(DisplayMode)=="Gmsh" )then call GmshExportStress(obj%FEMDomain,obj%DeformVecGloTot,obj%DeformStress,obj%DeformStrain,step,Name=Name) call obj%getDispVector(DispVector) call GmshPlotVector(obj%FEMDomain,Vector=DispVector,Name=Name,FieldName="Displacement",& Step=step,withMsh=.true.,NodeWize=.true.,onlyDirichlet=.true.) elseif(trim(DisplayMode)=="gnuplot" .or. trim(DisplayMode)=="Gnuplot" )then call GnuplotExportStress(obj%FEMDomain,obj%DeformVecGloTot,obj%DeformStress,obj%DeformStrain,step ) else print *, "Invalid DisplayMode >> DisplayDiffusionEq " print *, "DisplayMode '",trim(DisplayMode),"' is not defined" endif return endif elseif(size(obj%FEMDomain%Mesh%NodCoord,2)==3)then ! 3-D condition if(present(OptionalStep) )then step=OptionalStep else step=1 endif if(present(DisplayMode) )then if(trim(DisplayMode)=="Terminal" .or. trim(DisplayMode)=="terminal")then do i=1,size(obj%DeformVecEBEInc,1) print *, obj%DeformVecEBEInc(i,:) enddo elseif(trim(DisplayMode)=="gmsh" .or. trim(DisplayMode)=="Gmsh" )then call GmshExportStress(obj%FEMDomain,obj%DeformVecGloTot,obj%DeformStress,obj%DeformStrain,step,Name=Name ) if(present(withDirichlet) )then if(withDirichlet .eqv. .true.)then ! Export dirichlet (Deformation) boundary conditions call obj%getDBCVector(DBCvec) call GmshPlotVector(obj=obj%FEMDomain,Vector=DBCvec,Name=Name,Step=step,FieldName="DispBoundary",NodeWize=.true. ) endif endif elseif(trim(DisplayMode)=="gnuplot" .or. trim(DisplayMode)=="Gnuplot" )then call GnuplotExportStress(obj%FEMDomain,obj%DeformVecGloTot,obj%DeformStress,obj%DeformStrain,step ) else print *, "Invalid DisplayMode >> DisplayDiffusionEq " print *, "DisplayMode '",trim(DisplayMode),"' is not defined" endif return endif else stop "only 2D and 3D can be displayed" endif end subroutine !######################## Display Finite Deformation ######################## !############# Get Traction Vector ###################### subroutine GetTractionVector(obj) class(FiniteDeform_),intent(inout)::obj integer(int32) :: i,j,dim_num,node if(allocated(obj%TractionVecGlo) .and. size(obj%DeformVecGloTot )/=size(obj%TractionVecGlo) )then deallocate(obj%TractionVecGlo) endif if(.not.allocated(obj%TractionVecGlo))then allocate(obj%TractionVecGlo(size(obj%DeformVecGloTot ) ) ) endif obj%TractionVecGlo(:)=0.0d0 if(allocated(obj%FEMDomain%Boundary%NBoundNodID) )then dim_num=size(obj%FEMDomain%Boundary%NBoundNodID,2 ) if(dim_num==0)then return endif do i=1, size(obj%FEMDomain%Boundary%NBoundNodID,1 ) do j=1, size(obj%FEMDomain%Boundary%NBoundNodID,2 ) if(obj%FEMDomain%Boundary%NBoundNodID(i,j)>0 )then node=obj%FEMDomain%Boundary%NBoundNodID(i,j) obj%TractionVecGlo(dim_num*(node-1) + j )=obj%TractionVecGlo(dim_num*(node-1) + j ) + & obj%FEMDomain%Boundary%NBoundVal(i,j) endif enddo enddo endif end subroutine !############# Get Traction Vector ###################### !############# Get Internal Vector ###################### subroutine GetInternalVector(obj) class(FiniteDeform_),intent(inout)::obj integer(int32) :: i,j,dim_num,node_num,nodeid1,node_num_elmtl,elem_num,dim1 node_num=size(obj%FEMDomain%Mesh%NodCoord,1) dim_num=size(obj%FEMDomain%Mesh%NodCoord,2) elem_num=size(obj%FEMDomain%Mesh%ElemNod,1) node_num_elmtl=size(obj%FEMDomain%Mesh%ElemNod,2) if(allocated(obj%InternalVecGlo) .and. size(obj%DeformVecGloTot )/=size(obj%InternalVecGlo) )then deallocate(obj%InternalVecGlo) endif if(.not.allocated(obj%InternalVecGlo))then allocate(obj%InternalVecGlo(size(obj%DeformVecGloTot ) ) ) endif obj%InternalVecGlo(:)=0.0d0 do i=1,elem_num do j=1, node_num_elmtl nodeid1=obj%FEMDomain%Mesh%ElemNod(i,j) do dim1=1,dim_num obj%InternalVecGlo(dim_num*(nodeid1-1) + dim1 )=& obj%InternalVecGlo(dim_num*(nodeid1-1) + dim1 ) + & (obj%DeformStressRHS(i, (dim_num*(j-1)) + dim1) ) enddo enddo enddo end subroutine !############# Get Internal Vector ###################### !############# Get Residual Vector ###################### subroutine GetResidualVector(obj) class(FiniteDeform_),intent(inout)::obj integer(int32) :: i,j,dim_num,node if(allocated(obj%ResidualVecGlo) .and. size(obj%DeformVecGloTot )/=size(obj%ResidualVecGlo) )then deallocate(obj%ResidualVecGlo) endif if(.not.allocated(obj%ResidualVecGlo))then allocate(obj%ResidualVecGlo(size(obj%DeformVecGloTot ) ) ) endif obj%ResidualVecGlo(:)=obj%TractionVecGlo(:)-obj%InternalVecGlo(:) end subroutine !############# Get Residual Vector ###################### subroutine UpdateStressMeasure(obj) class(FiniteDeform_),intent(inout)::obj if(obj%infinitesimal .eqv. .true.)then obj%DeformStress(:,:,:)=obj%DeformStress(:,:,:)+obj%DeformStressinc(:,:,:) obj%DeformStressinc(:,:,:)=0.0d0 endif end subroutine !############# Uodate strain variables ###################### subroutine UpdateStrainMeasure(obj) class(FiniteDeform_),intent(inout)::obj integer(int32) :: i,j,dim_num dim_num=size(obj%FEMDomain%Mesh%NodCoord,2) obj%VolInitCurrEBE(:,3)=obj%VolInitCurrEBE(:,2)-obj%VolInitCurrEBE(:,1) obj%VolInitCurrEBE(:,1)=obj%VolInitCurrEBE(:,2) do i=1,size(obj%DeformStrain,1) do j=1,size(obj%DeformStrain,2) obj%DeformStrain(i,j, 1)=obj%DeformStrain(i,j, 4) obj%DeformStrain(i,j, 2)=obj%DeformStrain(i,j, 5) obj%DeformStrain(i,j, 3)=obj%DeformStrain(i,j, 6) obj%DeformStrain(i,j, 7)=obj%DeformStrain(i,j, 11) obj%DeformStrain(i,j, 8)=obj%DeformStrain(i,j, 12) obj%DeformStrain(i,j, 9)=obj%DeformStrain(i,j, 13) obj%DeformStrain(i,j, 10)=obj%DeformStrain(i,j, 14) if(dim_num==3)then obj%DeformStrain(:,:,16) = obj%DeformStrain(:,:,19) obj%DeformStrain(:,:,17) = obj%DeformStrain(:,:,20) obj%DeformStrain(:,:,18) = obj%DeformStrain(:,:,21) obj%DeformStrain(:,:,22) = obj%DeformStrain(:,:,27) obj%DeformStrain(:,:,23) = obj%DeformStrain(:,:,28) obj%DeformStrain(:,:,24) = obj%DeformStrain(:,:,29) obj%DeformStrain(:,:,25) = obj%DeformStrain(:,:,30) obj%DeformStrain(:,:,26) = obj%DeformStrain(:,:,31) obj%DeformStrain(:,:,32) = obj%DeformStrain(:,:,41) obj%DeformStrain(:,:,33) = obj%DeformStrain(:,:,42) obj%DeformStrain(:,:,34) = obj%DeformStrain(:,:,43) obj%DeformStrain(:,:,35) = obj%DeformStrain(:,:,44) obj%DeformStrain(:,:,36) = obj%DeformStrain(:,:,45) obj%DeformStrain(:,:,37) = obj%DeformStrain(:,:,46) obj%DeformStrain(:,:,38) = obj%DeformStrain(:,:,47) obj%DeformStrain(:,:,39) = obj%DeformStrain(:,:,48) obj%DeformStrain(:,:,40) = obj%DeformStrain(:,:,49) endif enddo enddo end subroutine !############# Uodate strain variables ###################### !############# Reaction Force at Loading Dirichlet Boundary ###################### subroutine DisplayReactionForce(obj,id) class(FiniteDeform_),intent(in)::obj integer(int32),optional,intent(in) :: id integer(int32) :: i,j,k,dim_num,dbc_num character(200) :: filename real(real64),allocatable :: ReactionForce(:),ReactionForce_wall(:) real(real64) :: val dim_num=size(obj%FEMDomain%Boundary%DBoundNodID,2) dbc_num=size(obj%FEMDomain%Boundary%DBoundNodID,1) allocate(ReactionForce(dim_num),ReactionForce_wall(dim_num) ) ReactionForce(:)=0.0d0 ReactionForce_wall(:)=0.0d0 do i=1,dim_num do j=1,dbc_num k=obj%FEMDomain%Boundary%DBoundNodID(j,i) val=obj%FEMDomain%Boundary%DBoundValInc(j,i) if(k<1)then cycle elseif(k>=1 .and. Val/=0.0d0 )then ReactionForce(i)=ReactionForce(i) + obj%InternalVecGlo( dim_num*(k-1) + i ) elseif(k>=1 .and. Val==0.0d0 )then ReactionForce_wall(i)=ReactionForce_wall(i) + abs(obj%InternalVecGlo( dim_num*(k-1) + i )) else cycle endif enddo enddo k=input(default=0,option=id) filename="ReactionForce"//trim(str(k))//"_wall.txt" if(obj%Step==1)then open(101,file="all_"//trim(filename),status="replace") else open(101,file="all_"//trim(filename),position="append") endif write(101,*) obj%Step,ReactionForce(:) close(101) if(obj%Step==1)then open(101,file=trim(filename),status="replace") else open(101,file=trim(filename),position="append") endif write(101,*) obj%Step,ReactionForce_wall(:) close(101) end subroutine !############# Reaction Force at Loading Dirichlet Boundary ###################### ! ################################################## subroutine getDBCVectorDeform(obj,DBCvec) class(FiniteDeform_),intent(in)::obj real(real64),allocatable,intent(inout)::DBCvec(:,:) integer(int32) :: i,j,n,m,k,l n=size(obj%FEMDomain%Mesh%NodCoord,1) m=size(obj%FEMDomain%Mesh%NodCoord,2) if(.not. allocated(DBCvec ) )then allocate(DBCvec(n,m) ) DBCvec(:,:)=0.0d0 endif ! check number of DBC do i=1,size(obj%FEMDomain%Boundary%DBoundNum) k=countif(Array=obj%FEMDomain%Boundary%DBoundNodID(:,i),Value=-1,notEqual=.true.) l=obj%FEMDomain%Boundary%DBoundNum(i) if(k /= l)then print *, "Caution :: FiniteDeformationClass::getDBCVector :: check number of DBC :: k /= l" endif enddo do i=1,size(obj%FEMDomain%Boundary%DBoundNodID,1) do j=1,size(obj%FEMDomain%Boundary%DBoundNodID,2) if(obj%FEMDomain%Boundary%DBoundNodID(i,j) <=0)then cycle endif DBCvec(obj%FEMDomain%Boundary%DBoundNodID(i,j),j )=obj%FEMDomain%Boundary%DBoundVal(i,j) enddo enddo end subroutine ! ################################################## ! ################################################## subroutine getDispVectorDeform(obj,Vector) class(FiniteDeform_),intent(in)::obj real(real64),allocatable,intent(inout)::Vector(:,:) integer(int32) :: i,j,n,m n=size(obj%FEMDomain%Mesh%NodCoord,1) m=size(obj%FEMDomain%Mesh%NodCoord,2) if(.not.allocated(Vector) )then allocate(Vector(n,m) ) Vector(:,:)=0.0d0 endif do i=1,n do j=1,m Vector(i,j)=obj%DeformVecGloTot( m*(i-1)+j) enddo enddo end subroutine ! ################################################## ! ################################################## subroutine exportFiniteDeform(obj,restart,path) class(FiniteDeform_),intent(inout) :: obj logical,optional, intent(in) :: restart character(*),intent(in) :: path type(IO_) :: f, ff integer(int32) :: fh_ if(present(restart) )then ! make new dir call execute_command_line("mkdir -p "//trim(path)//"/FiniteDeform") if(associated(obj%FEMDomain) )then call obj%FEMDomain%export(path=path//"/FiniteDeform",restart=restart) endif call f%open("./",trim(path)//"/FiniteDeform","/FiniteDeform.prop") write(f%fh, *) obj%DeformStress(:,:,:) write(f%fh, *) obj%DeformStrain(:,:,:) write(f%fh, *) obj%DeformStressInit(:,:,:) write(f%fh, *) obj%DeformStressinc(:,:,:) write(f%fh, *) obj%DeformStressMat(:,:,:) write(f%fh, *) obj%DeformStressRHS(:,:) write(f%fh, *) obj%DeformVecEBETot(:,:) write(f%fh, *) obj%DeformVecEBEInc(:,:) write(f%fh, *) obj%DeformVecGloTot(:) write(f%fh, *) obj%DeformVecGloInc(:) write(f%fh, *) obj%TractionVecGlo(:) write(f%fh, *) obj%ResidualVecGlo(:) write(f%fh, *) obj%InternalVecGlo(:) write(f%fh, *) obj%VolInitCurrEBE(:,:) write(f%fh, *) obj%YoungsModulus(:) write(f%fh, *) obj%PoissonsRatio(:) write(f%fh, *) obj%PorePressure(:) write(f%fh, *) obj%dt,obj%error,obj%reactionforce write(f%fh, *) obj%nr_tol write(f%fh, *) obj%ReducedIntegration write(f%fh, *) obj%infinitesimal write(f%fh, *) obj%itr write(f%fh, *) obj%Step call f%close() endif end subroutine ! ################################################## ! ################################################## function getVolumeDeform(obj) result(volume) class(FiniteDeform_),intent(inout) :: obj real(real64),allocatable :: volume(:) integer(int32) :: numnode, numelem,i,j,gp_num obj%FEMDomain%ShapeFunction%ElemType=obj%FEMDomain%Mesh%GetElemType() call SetShapeFuncType(obj%FEMDomain%ShapeFunction) gp_num=obj%FEMDomain%ShapeFunction%NumOfGp allocate(volume(size(obj%FEMDomain%Mesh%ElemNod,1)) ) do i = 1, size(obj%FEMDomain%Mesh%ElemNod,1) do j = 1, obj%FEMDomain%ShapeFunction%NumOfGp !�K�E�X�ϕ����ƃ��[�v ! -----J�}�g���N�X�̌v�Z----------------------------------------- call GetAllShapeFunc(obj%FEMDomain%ShapeFunction,elem_id=i,nod_coord=obj%FEMDomain%Mesh%NodCoord,& elem_nod=obj%FEMDomain%Mesh%ElemNod,OptionalGpID=j) volume(i) = obj%FEMDomain%ShapeFunction%detJ enddo enddo end function getVolumeDeform ! ################################################## end module FiniteDeformationClass