ExampleAxialStressedBeam: AxialStressedBeam.py

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