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
    ref value
    Ludovic Noels authored
    745baaa7
    History
    ludovic noels's avatar 745baaa7
    History
    composite.py 3.57 KiB
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128 
    

    

    

    #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)