Skip to content
Snippets Groups Projects
Select Git revision
  • 3f73249a8770896b03a093b3fe10bb7cc62f42e0
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

PViewVertexArrays.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    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)