Skip to content
Snippets Groups Projects
Select Git revision
  • 745baaa7981342c3d4036f09e829324ba25cc850
  • master default protected
  • ujwal_21_08_2024
  • dev_mm_pf
  • cyrielle
  • vinayak
  • dev_mm_torchSCRU
  • debug_mm_pf
  • newStructureNonLocal
  • Mohamed_stochasticDMN
  • dev_mm_bench
  • stochdmn
  • revert-351ff7aa
  • ujwal_29April2024
  • dev_mm_ann
  • mohamed_vevp
  • ujwal_debug
  • ujwal_2ndApril2024
  • ujwal_October_2023
  • gabriel
  • SFEM
  • v4.0
  • v3.2.3_multiplePhase
  • v3.5
  • v3.3.2
  • v3.4
  • v3.3
  • ver3.2
  • verJulienWork
  • ver3.1
  • ver2
  • ver1.1.2
  • ver1.1.1
  • ver1.1
34 results

composite.py

Blame
  • ludovic noels's avatar
    Ludovic Noels authored
    745baaa7
    History
    composite.py 3.57 KiB
    #coding-Utf-8-*-
    from gmshpy import *
    from dG3Dpy import*
    
    #script to launch beam problem with a python script
    
    # material law
    lawnumEP = 2 
    rhoEP   = 1000.
    youngEP = 3.2e9
    nuEP    = 0.3
    sy0EP   = 25.e6
    hEP     = 7.10e9
    
    lawnumF = 3 
    rhoF   = 1000.
    youngF = 320e9
    nuF    = 0.3
    sy0F   = 25.e16
    hF     = 7.10e9 
    
    lawcnumEP = 4
    GcEP   =  78.
    sigmacEP = 83.e6
    deltacEP = 2*GcEP/sigmacEP
    KcEP = sigmacEP/(0.005*deltacEP)
    
    
    lawcnumInt = 5
    GcInt   = 60.
    sigmaInt = 30.e6
    deltacInt= 2*GcInt/sigmaInt
    KcInt = sigmaInt/(0.005*deltacInt)
    
    beta =0.87 # ratio KII/KI
    mu = 0. # friction coefficient ??
    fsmin = 1.
    fsmax = 1.
    
    # geometry
    meshfile="composite.msh" # name of mesh file
    
    # solver
    sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
    soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
    nstep = 500   # number of step (used only if soltype=1)
    ftime =1 # Final time (used only if soltype=1)
    tol=1.e-4   # relative tolerance for NR scheme (used only if soltype=1)
    nstepArch=10 # Number of step between 2 archiving (used only if soltype=1)
    fullDg =1 #O = CG, 1 = DG
    space1 = 0 # function space (Lagrange=0)
    beta1  = 100
    
    #  compute solution and BC (given directly to the solver
    # creation of law
    law1EP = J2LinearDG3DMaterialLaw(lawnumEP,rhoEP,youngEP,nuEP,sy0EP,hEP)
    law2EP = LinearCohesive3DLaw(lawcnumEP,GcEP,sigmacEP,beta,mu,fsmin,fsmax,KcEP)
    law3EP = FractureByCohesive3DLaw(6,lawnumEP,lawcnumEP)
    
    lawF = J2LinearDG3DMaterialLaw(lawnumF,rhoF,youngF,nuF,sy0F,hF)
    law2Int = LinearCohesive3DLaw(lawcnumInt,GcInt,sigmaInt,beta,mu,fsmin,fsmax,KcInt)
    
    # creation of ElasticField
    
    myfieldEP = dG3DDomain(10,78,space1,6,fullDg)
    #myfieldEP.matrixByPerturbation(1,1,1,1e-8)
    myfieldEP.stabilityParameters(beta1)
    
    myfieldF = dG3DDomain(11,77,space1,lawnumF,fullDg)
    #myfieldF.matrixByPerturbation(1,1,1,1e-8)
    myfieldF.stabilityParameters(beta1)
    
    #in // we need to create the ininterdomain after adding the domains
    myinterfield = interDomainBetween3D(12,myfieldEP,myfieldF,lawcnumInt)
    myinterfield.stabilityParameters(beta1)
    #myinterfield.matrixByPerturbation(1,1e-10)
    
    
    # creation of Solver
    mysolver = nonLinearMechSolver(1000)
    mysolver.loadModel(meshfile)
    mysolver.addDomain(myfieldEP)
    mysolver.addDomain(myfieldF)
    mysolver.addMaterialLaw(law1EP)
    mysolver.addMaterialLaw(law2EP)
    mysolver.addMaterialLaw(law3EP)
    mysolver.addMaterialLaw(lawF)
    
    mysolver.addDomain(myinterfield)
    mysolver.addMaterialLaw(law2Int)
    mysolver.Scheme(soltype)
    mysolver.Solver(sol)
    mysolver.snlData(nstep,ftime,tol,1e-12)
    mysolver.stepBetweenArchiving(nstepArch)
    
    
    #mysolver.explicitSpectralRadius(ftime,0.5,0.)
    #mysolver.dynamicRelaxation(0.1, ftime, 1.e-3,2)
    #mysolver.explicitTimeStepEvaluation(nstep)
    
    mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package petsc")
    
    # BC
    mysolver.displacementBC("Volume",77,2,0.)
    mysolver.displacementBC("Volume",78,2,0.)
    mysolver.displacementBC("Face",82,0,0.)
    mysolver.displacementBC("Face",81,1,0.)
    mysolver.displacementBC("Face",80,0,5e-6)
    
    
    
    mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
    mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
    mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
    mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
    mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
    mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
    mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
    mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
    
    mysolver.archivingForceOnPhysicalGroup("Face", 82, 0, 1)
    mysolver.archivingNodeDisplacement(86,0,1)
    
    mysolver.solve()
    
    check = TestCheck()
    check.equal(1.056398e-08,mysolver.getEnergy(3),5.e-3)