Skip to content
Snippets Groups Projects
Select Git revision
  • f62a888ab88477a8bbeb04a458fe06eb2f706591
  • master default protected
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

polynomialBasis.h

Blame
  • cube.py 7.93 KiB
    #coding-Utf-8-*-
    from gmshpy import*
    from dG3Dpy import*
    from math import*
    
    # Script for plane notched specimen with Gurson
    
    
    # Material law creation
    # ===================================================================================
    BulkLawNum1 = 11
    
    rho   = 7850
    young = 205.4e3 #MPa
    nu    = 0.3 
    sy0   = 856.4 #MPa
    h = 3399.
    p1 = 0.012
    n1 = 0.084
    p2 = 0.056
    n2 = 0.058
    
    harden_Law1 = LinearFollowedByMultiplePowerLawJ2IsotropicHardening(BulkLawNum1,sy0,h,p1,n1,p2,n2)
    
    # - non-local law
    l = 1e3
    cl = l*l
    Cl_Law1 = IsotropicCLengthLaw(BulkLawNum1, cl)
    
    # - gurson
    q1    = 1.5
    q2    = 1.
    q3    = 1.5 
    fVinitial = 5e-3 # initial porosity
    ffstar = (q1- sqrt(q1*q1-q3*q3))/(q3*q3) # critical porosity 0.66667
    print("critical porosity = ", ffstar)
    
    # material law creation
    BulkLaw1 = NonLocalPorousCoupledDG3DMaterialLaw(BulkLawNum1, young, nu, rho, q1,q2,q3, fVinitial, 3., harden_Law1, Cl_Law1)
    BulkLaw1.setOrderForLogExp(9)
    BulkLaw1.setShearPorosityGrowthFactor(10.)
    
    # - nucleation
    function1 = LinearNucleationFunctionLaw(0, 0.5)
    criterion1 = PorosityBasedNucleationCriterionLaw(0, 0.007, 1.)
    BulkLaw1.setNucleationLaw(function1, criterion1)
    
    BulkLaw1.setCrackTransition(True)
    BulkLaw1.setPostBlockedDissipationBehavior(1)
    
    interfaceLawNum1 = 2
    lBand = 1. # m
    Kp = 1e12 # N/m
    cohLaw1 = PorosityCohesiveBand3DLaw(interfaceLawNum1,lBand,Kp,True)
    
    fracLawNum1 = 6
    fracLaw1 = FractureByCohesive3DLaw(fracLawNum1,BulkLawNum1,interfaceLawNum1)
    
    
    # Solver parameters - implicit
    # ===================================================================================
    soltype = 1 		# StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, # Multi=3, Implicit=4, Eigen=5)
    nstep = 1000 		# Number of step
    ftime =1.0 		# Final time
    tol=1.e-6 		# Relative tolerance for NR scheme
    tolAbs = 1.e-8		# Absolute tolerance for NR scheme
    nstepArchIP=1		# Number of step between 2 archiving
    nstepArchForce = 1	# Number of step between 2 force archiving
    nstepArchEnergy = nstepArchForce # Number of step between 2 energy computation and archiving
    
    MaxIter = 12		# Maximum number of iterations
    StepIncrease = 2	# Number of successfull timestep before reincreasing it
    StepReducFactor = 2 	# Timestep reduction factor
    NumberReduction = 100	# Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)
    
    fullDg = True
    eqRatio = 1.e2		# Ratio between "elastic" and non-local equations 
    space1 = 0 		# Function space (Lagrange=0)
    beta1  = 100.0		# Penality parameter for DG
    
    # Domain creation
    
    meshfile= "cube.msh" 	# name of mesh file
    
    ## ===================================================================================
    field1 = dG3DDomain(11,11,0,fracLawNum1,fullDg,3,1)
    field1.stabilityParameters(beta1) 			# Adding stability parameters (for DG)
    field1.setNonLocalStabilityParameters(beta1,fullDg) 		# Adding stability parameters (for DG)
    field1.setNonLocalEqRatio(eqRatio)
    field1.setBulkDamageBlockedMethod(0)
    field1.averageStrainBased(False)
    
    # Solver creation
    # ===================================================================================
    mysolver = nonLinearMechSolver(1000) 		# Solver associated with numSolver1
    
    mysolver.loadModel(meshfile)		# add mesh
    mysolver.addDomain(field1) 		# add domain
    mysolver.addMaterialLaw(BulkLaw1) 	# add material law
    mysolver.addMaterialLaw(cohLaw1) 	# add material law
    mysolver.addMaterialLaw(fracLaw1) 	# add material law
    
    # solver parametrisation
    mysolver.Scheme(soltype) 		# solver scheme
    mysolver.Solver(2) 			# Library solver: Gmm=0 (default) Taucs=1 PETsc=2
    mysolver.snlData(nstep,ftime,tol,tolAbs) # solver parameters for imp. solving (numb. of step, final time and tolerance)
    mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction)
    
    # solver archiving
    mysolver.stepBetweenArchiving(nstepArchIP) 	# archiving frequency
    mysolver.energyComputation(nstepArchEnergy)	# archiving frequency for energy
    
    
    # Boundary conditions
    # ===============================
    mysolver.initialBC("Volume","Position",11,3,fVinitial)
    
    mysolver.displacementBC("Face",1,0,0.)
    mysolver.displacementBC("Face",3,1,0.)
    mysolver.displacementBC("Face",5,2,0.)
    
    mysolver.sameDisplacementBC("Face",2,7,0)
    mysolver.sameDisplacementBC("Face",4,7,1)
    mysolver.sameDisplacementBC("Face",6,7,2)
    
    method=0
    mysolver.pathFollowing(True,method)
    if method==0:
    	mysolver.setPathFollowingIncrementAdaptation(True,6,3)
    	mysolver.setPathFollowingControlType(0)
    	mysolver.setPathFollowingCorrectionMethod(0)
    	mysolver.setPathFollowingArcLengthStep(2e-1)
    	mysolver.setBoundsOfPathFollowingArcLengthSteps(0,5e-1);
    elif method==1:
    	mysolver.setPathFollowingIncrementAdaptation(False,5)
    	mysolver.setPathFollowingLocalSteps(1e-1,1e-2)
    	mysolver.setPathFollowingLocalIncrementType(1); 
    	mysolver.setPathFollowingSwitchCriterion(0)
    
    T = 2
    fact1 = 0.3; # (T-1./3.)/(T+2./3)
    fact2 = 0.3; #
    mysolver.pressureOnPhysicalGroupBC("Face",2,1e3,0.)
    mysolver.pressureOnPhysicalGroupBC("Face",4,fact1*1e3,0.)
    mysolver.pressureOnPhysicalGroupBC("Face",6,fact2*1e3,0.)
    
    # Variable storage
    # ===============================
    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("damage",IPField.DAMAGE,1,1)
    mysolver.internalPointBuildView("local_fV",IPField.LOCAL_POROSITY,1,1)
    mysolver.internalPointBuildView("local_fV_max",IPField.LOCAL_POROSITY,1,IPField.MAX_VALUE)
    mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
    mysolver.internalPointBuildView("triaxiality",IPField.STRESS_TRIAXIALITY, 1, 1)
    mysolver.internalPointBuildView("LODE_PARAMETER",IPField.LODE_PARAMETER, 1, 1)
    mysolver.internalPointBuildView("FIRST_PRINCIPAL_STRESS",IPField.FIRST_PRINCIPAL_STRESS, 1, 1)
    mysolver.internalPointBuildView("SECOND_PRINCIPAL_STRESS",IPField.SECOND_PRINCIPAL_STRESS, 1, 1)
    mysolver.internalPointBuildView("THIRD_PRINCIPAL_STRESS",IPField.THIRD_PRINCIPAL_STRESS, 1, 1)
    mysolver.internalPointBuildView("COALESCENCE",IPField.COALESCENCE, 1, 1)
    mysolver.internalPointBuildView("NONLOCAL_POROSITY",IPField.NONLOCAL_POROSITY, 1, 1)
    mysolver.internalPointBuildView("CORRECTED_POROSITY",IPField.CORRECTED_POROSITY, 1, 1)
    mysolver.internalPointBuildView("CORRECTED_POROSITY_MAX",IPField.CORRECTED_POROSITY, 1, IPField.MAX_VALUE)
    mysolver.internalPointBuildView("POROSITY_AT_COALESCENCE",IPField.POROSITY_COALESCENCE_ONSET, 1, 1)
    mysolver.internalPointBuildView("FAILED",IPField.FAILED, 1, 1)
    
    mysolver.OneUnknownBuildView("nonlocalVar",3,1)
    
    mysolver.archivingForceOnPhysicalGroup("Face", 1, 0, nstepArchForce)
    mysolver.archivingForceOnPhysicalGroup("Face", 3, 1, nstepArchForce)
    mysolver.archivingForceOnPhysicalGroup("Face", 5, 2, nstepArchForce)
    mysolver.archivingNodeDisplacement(7,0, nstepArchForce)
    mysolver.archivingNodeDisplacement(7,1, nstepArchForce)
    mysolver.archivingNodeDisplacement(7,2, nstepArchForce)
    
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.NONLOCAL_POROSITY,IPField.MAX_VALUE,1)
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.LOCAL_POROSITY,IPField.MAX_VALUE,1)
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.CORRECTED_POROSITY,IPField.MAX_VALUE,1)
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.YIELD_POROSITY,IPField.MAX_VALUE,1)
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.VOLUMETRIC_PLASTIC_STRAIN,IPField.MAX_VALUE,1)
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.VOLUMETRIC_PLASTIC_STRAIN,IPField.MAX_VALUE,1)
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.SVM,IPField.MAX_VALUE,1)
    mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.LIGAMENT_RATIO,IPField.MAX_VALUE,1)
    
    mysolver.solve()
    
    check = TestCheck()
    import csv
    data = csv.reader(open('IPVolume11val_LOCAL_POROSITYMax.csv'), delimiter=';')
    porosity = list(data)
    check.equal(2.660465e-02,float(porosity[-1][1]),1e-6)