Skip to content
Snippets Groups Projects
Select Git revision
  • 30965a9700675dd7cda0947e442b598385818ee0
  • master default protected
  • hierarchical-basis
  • alphashapes
  • bl
  • 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
  • tmp_jcjc24
  • fixedMeshIF
  • save_edges
  • 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

Eigenvalues.cpp

Blame
  • beam.py 3.60 KiB
    #-*-coding:Utf-8-*-
    from gmshpy import *
    from dgshellpy import *
    # bending beam
    
    # material law
    lawnum = 1 #  unique number of the law
    E = 100.e9 # Young's modulus
    nu = 0.3   #  Poisson's ratio
    rho =7850. 
    # geometry
    h = 0.01  # thickness
    geofile="beam.geo"
    meshfile="beam.msh" # name of mesh file
    # integration
    nsimp = 3 #  number of Simpson's points (odd)
    
    # solver
    sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
    beta1 = 10. # value of stabilization parameter
    beta2 = 50.
    beta3 = 0.1
    soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
    nstep = 1   # number of step (used only if soltype=1)
    ftime =1.   #  Final time (used only if soltype=1)
    tol=1.e-3  #  relative tolerance for NR scheme (used only if soltype=1)
    nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
    
     # compute solution and BC (given directly to the solver
    
    #  creation of law
    law1 = linearElasticShellLaw(lawnum,E,nu,h,nsimp,rho)
    
    #  creation of ElasticField
    nfield =99 # number of the field (physical number of surface)
    fullDg = 0 #  formulation CgDg=0 fullDg =1
    space1 = 0 # function space Lagrange=0
    myfield1 = dgLinearShellDomain(1000,nfield,space1,lawnum,fullDg)
    myfield1.stabilityParameters(beta1,beta2,beta3)
    #myfield1.gaussIntegration(dgLinearShellDomain.Lobatto,5,3)
    #  creation of Solver
    mysolver = nonLinearMechSolver(1000)
    mysolver.createModel(geofile,meshfile,2,2,False)
    mysolver.addDomain(myfield1)
    mysolver.addMaterialLaw(law1)
    mysolver.Scheme(soltype)
    mysolver.Solver(sol)
    mysolver.snlData(nstep,ftime,tol)
    mysolver.stepBetweenArchiving(nstepArch)
    # BC
    mysolver.displacementBC("Edge",41,0,0.)
    mysolver.displacementBC("Edge",41,1,0.)
    mysolver.displacementBC("Edge",41,2,0.)
    #mysolver.displacementBC("Edge",21,0,0.)
    #mysolver.displacementBC("Edge",21,1,0.)
    #mysolver.displacementBC("Edge",21,2,0.)
    #mysolver.displacementBC("Edge",21,2,0.4)
    #mysolver.forceBC("Face",99,2,1000.)
    fct1= simpleFunctionTimeDouble(-1.e6)
    mysolver.forceBC("Edge",21,2,fct1)
    #mysolver.displacementBC("Face",99,1,0.)
    #mysolver.displacementBC("Face",99,0,0.)
    # mysolver.pressureOnPhysicalGroupBC(99,1000.)
    mysolver.thetaBC(41)
    #mysolver:AddThetaConstraint(21)
    
    mysolver.archivingForceOnPhysicalGroup("Edge",41,2)
    mysolver.archivingNodeDisplacement(22,2)
    mysolver.archivingNodeDisplacement(22,1)
    mysolver.archivingNodeDisplacement(22,0)
    mysolver.solve()
    
    # save value for later check
    fprevious = mysolver.getArchivedForceOnPhysicalGroup("Edge",41,2)
    eprevious = mysolver.getEnergy(energeticField.deformation)
    
    # explicit analysis
    ftime = 0.000001
    gamma_s = 0.06666666
    rhoinfty = 0.99
    nstep = 1
    nstepArch = 1
    mysolver.Scheme(2)
    #mysolver.Solver(0)
    mysolver.explicitSpectralRadius(ftime,gamma_s,rhoinfty)
    mysolver.explicitTimeStepEvaluation(nstep)
    mysolver.stepBetweenArchiving(nstepArch)
    
    # new BC
    mysolver.resetBoundaryConditions()
    mysolver.displacementBC("Edge",41,0,0.)
    mysolver.displacementBC("Edge",41,1,0.)
    mysolver.displacementBC("Edge",41,2,0.)
    
    mysolver.archivingNodeDisplacement(22,2)
    mysolver.archivingNodeDisplacement(22,1)
    mysolver.archivingNodeDisplacement(22,0)
    
    mysolver.solve()
    
    # check 
    check = TestCheck()
    check.equal(9.999950e+03,mysolver.getArchivedForceOnPhysicalGroup('Edge',41,2))
    check.equal(-3.943707e-02,mysolver.getArchivedNodalValue(22,2,mysolver.displacement))
    # previous value
    try:
        import linecache
        lforce = linecache.getline('previousScheme/force41comp2.csv',1)
        lenergy= linecache.getline('previousScheme/energy.csv',3)
    except:
        print('Cannot get values in the files')
        import os
        os._exit(1)
    else:
        force = float(lforce.split(';')[1])
        energy = float(lenergy.split(';')[2])
        check.equal(fprevious,force,1.e-6)
        check.equal(eprevious,energy,1.e-6)