Skip to content
Snippets Groups Projects
Select Git revision
  • 49041a7f5059255b87ecc09d49873bd034e3c32a
  • master default protected
  • dev_mm_pf
  • cyrielle
  • ujwal_05_08_2025
  • vinayak
  • complexpetsc
  • ujwal_21_08_2024
  • 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
  • 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

cubeTransverseAnisotropy.py

Blame
    • cubeTransverseAnisotropy.py 3.43 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 
    

    

    

    #coding-Utf-8-*-
from gmshpy import *
from dG3Dpy import*

#script to launch beam problem with a python script

# material law
lawnum   = 1 # unique number of law
lawcoh   = 2 # unique number of law
lawfrac  = 3 # unique number of law
rho      = 7850
young    = 40.e9
youngA   = 230.e9;
nu       = 0.2 #
nu_minor = 0.256*young/youngA 
GA       = 24.e9

GcL      = 2.e3;
sigmacL  = 500.e6;
GcT      = 0.5e3;
sigmacT  = 50.e6;
beta     = 0.87;
mu       = -1.;
fsmin    = 0.999;
fsmax    = 1.001;


#longitudinal loading
Ax       = 1.; #direction of anisotropy
Ay       = 0.;
Az       = 1.;
d1=0.75e-5



# geometry
geofile="cube.geo"
meshfile="cube.msh" # name of mesh file

# solver
sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 7500   # number of step (used only if soltype=1)
ftime =1.   # Final time (used only if soltype=1)
tol=1.e-6   # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=50 # Number of step between 2 archiving (used only if soltype=1)
fullDg = 1 #O = CG, 1 = DG
space1 = 0 # function space (Lagrange=0)


#  compute solution and BC (given directly to the solver
# creation of law
law1 = TransverseIsotropicDG3DMaterialLaw(lawnum,rho,young,nu,youngA,GA,nu_minor,Ax,Ay,Az)
lawcoh1  = TransverseIsotropicLinearCohesive3D(lawcoh,GcL,sigmacL,GcT,sigmacT,beta,mu,Ax,Ay,Az,fsmin,fsmax)
lawfrac1 = FractureByCohesive3DLaw(lawfrac,lawnum,lawcoh)


# creation of ElasticField
nfield = 10 # number of the field (physical number of surface)
myfield1 = dG3DDomain(1000,nfield,space1,lawfrac,fullDg)
myfield1.stabilityParameters(30.)
myfield1.matrixByPerturbation(0,0,0,1e-8)
#myfield1.forceCohesiveInsertionAllIPs(True,0.)
# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.createModel(geofile,meshfile,3,2)
#mysolver.loadModel(meshfile)
mysolver.addDomain(myfield1)
mysolver.addMaterialLaw(law1)
mysolver.addMaterialLaw(lawcoh1)
mysolver.addMaterialLaw(lawfrac1)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.stepBetweenArchiving(nstepArch)
# BC
#shearing
#mysolver.displacementBC("Face",1234,0,0.)
#mysolver.displacementBC("Face",1234,1,0.)
#mysolver.displacementBC("Face",1234,2,0.)
#mysolver.displacementBC("Face",5678,1,0.00)
#mysolver.displacementBC("Face",5678,0,0.00001)
#mysolver.forceBC("Face",2376,2,100000000)
#mysolver.forceBC("Face",4158,2,-100000000)
#mysolver.forceBC("Face",5678,0,100000000)


#tension along z
cyclicFunction1=cycleFunctionTime(0.,0.,ftime/4., d1/2., ftime/2., d1/4., 3.*ftime/4., d1/2., ftime, d1);

mysolver.displacementBC("Face",1234,2,0.)
mysolver.displacementBC("Face",5678,2,cyclicFunction1)
mysolver.displacementBC("Face",2376,0,0.)
mysolver.displacementBC("Face",1265,1,0.)


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", 1234, 2)
mysolver.archivingForceOnPhysicalGroup("Face", 5678, 2)

mysolver.solve()

check = TestCheck()
check.equal(2.614308e-01,mysolver.getArchivedForceOnPhysicalGroup("Face", 5678, 2),1.e-3)