Select Git revision
Scal2Vec.cpp
-
Axel Modave authoredAxel Modave authored
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)