ExampleAxialStressedBeam: AxialStressedBeam-v1.1.py

File AxialStressedBeam-v1.1.py, 5.8 KB (added by michael.baudin@…, 4 years ago)

Updated graphics commands for OT v1.1

Line 
1# Main OpenTURNS module
2from openturns import *
3# OpenTURNS viewer capability (since release 0.9.2)
4# from openturns_viewer import ViewImage
5# Up to v1.0:
6#from openturns.viewer import ViewImage
7# Since v1.1:
8from openturns.viewer import View
9
10# Constant definition:
11stochasticDimension = 2
12
13# Random generator initialization:
14RandomGenerator().SetSeed(0)
15
16#Analytical model definition:
17
18# Analytical construction : Input
19inputFunction = Description(stochasticDimension)
20inputFunction[0] = "R"
21inputFunction[1] = "F"
22
23# Analytical construction : Output
24outputFunction = Description(1)
25outputFunction[0] = "G"
26
27formulas = Description(outputFunction.getSize())
28# Here, _pi is the built-in constant _pi=3.14159265....
29formulas[0] = "R-F/(_pi*100.0)"
30
31LimitState = NumericalMathFunction(inputFunction, outputFunction, formulas)
32
33# Test of the limit state function:
34x = NumericalPoint(stochasticDimension, 0)
35x[0]=300.0
36x[1]=75000.0
37print "x=" , x
38print "G(x)=" , LimitState(x)
39
40# Stochastic model definition:
41
42# Mean
43mean = NumericalPoint(stochasticDimension, 0.0)
44mean[0] = 300.0
45mean[1] = 75000.0
46
47# Standard deviation
48sigma = NumericalPoint(stochasticDimension, 0.0)
49sigma[0] = 30.0
50sigma[1] = 5000.0
51
52# Additional parameters for the lognormal distribution:
53component = Description(1)
54BorneInf = 0.0
55
56# Initialization of the distribution collection:
57aCollection = DistributionCollection()
58
59# Create a first marginal : LogNormal distribution 1D, parameterized by its mean and standard deviation
60marginal = LogNormal(mean[0], sigma[0], BorneInf, LogNormal.MUSIGMA)
61marginal.setName("Yield strength")
62component[0] = "R"
63marginal.setDescription(component)
64# Graphical output of the PDF
65pdfgraph1=marginal.drawPDF()
66# If the directory name is omitted, the graph will be produced in the current directory.
67# If the dimensions are omitted, the default is 400x300 up to release 0.10.0, and 640x480 for the later releases
68pdfgraph1.draw("pdf_R", 640, 480)
69# Visualize the graph
70# Up to OT v1.0 :
71# ViewImage(pdfgraph1.getBitmap())
72# With OT v1.1 :
73View(pdfgraph1)
74# Fill the first marginal of aCollection
75aCollection.add( Distribution(marginal, "Yield strength") )
76
77# Create a second marginal : Normal distribution 1D
78marginal = Normal(mean[1], sigma[1])
79marginal.setName("Traction_load")
80component[0] = "F"
81marginal.setDescription(component)
82# Graphical output of the PDF
83pdfgraph2=marginal.drawPDF()
84pdfgraph2.draw("pdf_F", 640, 480)
85View(pdfgraph2)
86# Fill the second marginal of aCollection
87aCollection.add( Distribution(marginal, "Traction_load") )
88
89# Create a copula : IndependentCopula (no correlation)
90aCopula = IndependentCopula(aCollection.getSize())
91aCopula.setName("Independent copula")
92
93# Instanciate one distribution object
94myDistribution = ComposedDistribution(aCollection, Copula(aCopula))
95myDistribution.setName("myDist")
96
97# We create a 'usual' RandomVector from the Distribution
98vect = RandomVector(Distribution(myDistribution))
99
100# We create a composite random vector
101G = RandomVector(LimitState, vect)
102
103# We create an Event from this RandomVector
104myEvent = Event(G, ComparisonOperator(Less()), 0.0, "Event 1")
105
106# Using Monte Carlo simulations
107
108# Resolution options:
109cv = 0.05
110NbSim = 100000
111
112algoMC = MonteCarlo(myEvent)
113algoMC.setMaximumOuterSampling(NbSim)
114algoMC.setBlockSize(1)
115algoMC.setMaximumCoefficientOfVariation(cv)
116# For statistics about the algorithm
117initialNumberOfCall = LimitState.getEvaluationCallsNumber()
118
119# Perform the analysis:
120algoMC.run()
121
122# Results:
123result = algoMC.getResult()
124probability = result.getProbabilityEstimate()
125print "MonteCarlo result=" , result
126print "Number of executed iterations =" , result.getOuterSampling()
127print "Number of calls to the limit state =" , LimitState.getEvaluationCallsNumber() - initialNumberOfCall
128print "Pf = " , probability
129print "CV =" , result.getCoefficientOfVariation()
130MCcvgraph = algoMC.drawProbabilityConvergence()
131MCcvgraph.draw("montecarlo_convergence", 640, 480)
132View(MCcvgraph)
133
134# Using FORM analysis
135
136# Resolution options:
137eps=1E-3
138
139# We create a NearestPoint algorithm
140myCobyla = Cobyla()
141myCobyla.setSpecificParameters(CobylaSpecificParameters())
142myCobyla.setMaximumIterationsNumber(100)
143myCobyla.setMaximumAbsoluteError(eps)
144myCobyla.setMaximumRelativeError(eps)
145myCobyla.setMaximumResidualError(eps)
146myCobyla.setMaximumConstraintError(eps)
147
148# For statistics about the algorithm
149initialNumberOfCall = LimitState.getEvaluationCallsNumber()
150
151# We create a FORM algorithm
152# The first parameter is a NearestPointAlgorithm
153# The second parameter is an event
154# The third parameter is a starting point for the design point research
155
156algoFORM = FORM(NearestPointAlgorithm(myCobyla), myEvent, mean)
157
158# Perform the analysis:
159algoFORM.run()
160
161# Results:
162result = algoFORM.getResult()
163print "Number of calls to the limit state =" , LimitState.getEvaluationCallsNumber() - initialNumberOfCall
164print "Pf =" , result.getEventProbability()
165
166# Graphical result output
167importanceFactorsGraph = result.drawImportanceFactors()
168importanceFactorsGraph.draw("importance_factors", 640, 480)
169View(importanceFactorsGraph)
170
171# Using Directional sampling
172
173# Resolution options:
174cv = 0.05
175NbSim = 100000
176
177algoDS = DirectionalSampling(myEvent)
178algoDS.setMaximumOuterSampling(NbSim)
179algoDS.setBlockSize(1)
180algoDS.setMaximumCoefficientOfVariation(cv)
181# For statistics about the algorithm
182initialNumberOfCall = LimitState.getEvaluationCallsNumber()
183
184# Perform the analysis:
185algoDS.run()
186
187# Results:
188result = algoDS.getResult()
189probability = result.getProbabilityEstimate()
190print "Number of executed iterations =" , result.getOuterSampling()
191print "Number of calls to the limit state =" , LimitState.getEvaluationCallsNumber() - initialNumberOfCall
192print "Pf = " , probability
193print "CV =" , result.getCoefficientOfVariation()
194DScvgraph = algoDS.drawProbabilityConvergence()
195DScvgraph.draw("directionalsampling_convergence", 640, 480)
196View(DScvgraph)
197