# ExampleAxialStressedBeam: AxialStressedBeam-v1.1.py

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

Line | |
---|---|

1 | # Main OpenTURNS module |

2 | from 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: |

8 | from openturns.viewer import View |

9 | |

10 | # Constant definition: |

11 | stochasticDimension = 2 |

12 | |

13 | # Random generator initialization: |

14 | RandomGenerator().SetSeed(0) |

15 | |

16 | #Analytical model definition: |

17 | |

18 | # Analytical construction : Input |

19 | inputFunction = Description(stochasticDimension) |

20 | inputFunction[0] = "R" |

21 | inputFunction[1] = "F" |

22 | |

23 | # Analytical construction : Output |

24 | outputFunction = Description(1) |

25 | outputFunction[0] = "G" |

26 | |

27 | formulas = Description(outputFunction.getSize()) |

28 | # Here, _pi is the built-in constant _pi=3.14159265.... |

29 | formulas[0] = "R-F/(_pi*100.0)" |

30 | |

31 | LimitState = NumericalMathFunction(inputFunction, outputFunction, formulas) |

32 | |

33 | # Test of the limit state function: |

34 | x = NumericalPoint(stochasticDimension, 0) |

35 | x[0]=300.0 |

36 | x[1]=75000.0 |

37 | print "x=" , x |

38 | print "G(x)=" , LimitState(x) |

39 | |

40 | # Stochastic model definition: |

41 | |

42 | # Mean |

43 | mean = NumericalPoint(stochasticDimension, 0.0) |

44 | mean[0] = 300.0 |

45 | mean[1] = 75000.0 |

46 | |

47 | # Standard deviation |

48 | sigma = NumericalPoint(stochasticDimension, 0.0) |

49 | sigma[0] = 30.0 |

50 | sigma[1] = 5000.0 |

51 | |

52 | # Additional parameters for the lognormal distribution: |

53 | component = Description(1) |

54 | BorneInf = 0.0 |

55 | |

56 | # Initialization of the distribution collection: |

57 | aCollection = DistributionCollection() |

58 | |

59 | # Create a first marginal : LogNormal distribution 1D, parameterized by its mean and standard deviation |

60 | marginal = LogNormal(mean[0], sigma[0], BorneInf, LogNormal.MUSIGMA) |

61 | marginal.setName("Yield strength") |

62 | component[0] = "R" |

63 | marginal.setDescription(component) |

64 | # Graphical output of the PDF |

65 | pdfgraph1=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 |

68 | pdfgraph1.draw("pdf_R", 640, 480) |

69 | # Visualize the graph |

70 | # Up to OT v1.0 : |

71 | # ViewImage(pdfgraph1.getBitmap()) |

72 | # With OT v1.1 : |

73 | View(pdfgraph1) |

74 | # Fill the first marginal of aCollection |

75 | aCollection.add( Distribution(marginal, "Yield strength") ) |

76 | |

77 | # Create a second marginal : Normal distribution 1D |

78 | marginal = Normal(mean[1], sigma[1]) |

79 | marginal.setName("Traction_load") |

80 | component[0] = "F" |

81 | marginal.setDescription(component) |

82 | # Graphical output of the PDF |

83 | pdfgraph2=marginal.drawPDF() |

84 | pdfgraph2.draw("pdf_F", 640, 480) |

85 | View(pdfgraph2) |

86 | # Fill the second marginal of aCollection |

87 | aCollection.add( Distribution(marginal, "Traction_load") ) |

88 | |

89 | # Create a copula : IndependentCopula (no correlation) |

90 | aCopula = IndependentCopula(aCollection.getSize()) |

91 | aCopula.setName("Independent copula") |

92 | |

93 | # Instanciate one distribution object |

94 | myDistribution = ComposedDistribution(aCollection, Copula(aCopula)) |

95 | myDistribution.setName("myDist") |

96 | |

97 | # We create a 'usual' RandomVector from the Distribution |

98 | vect = RandomVector(Distribution(myDistribution)) |

99 | |

100 | # We create a composite random vector |

101 | G = RandomVector(LimitState, vect) |

102 | |

103 | # We create an Event from this RandomVector |

104 | myEvent = Event(G, ComparisonOperator(Less()), 0.0, "Event 1") |

105 | |

106 | # Using Monte Carlo simulations |

107 | |

108 | # Resolution options: |

109 | cv = 0.05 |

110 | NbSim = 100000 |

111 | |

112 | algoMC = MonteCarlo(myEvent) |

113 | algoMC.setMaximumOuterSampling(NbSim) |

114 | algoMC.setBlockSize(1) |

115 | algoMC.setMaximumCoefficientOfVariation(cv) |

116 | # For statistics about the algorithm |

117 | initialNumberOfCall = LimitState.getEvaluationCallsNumber() |

118 | |

119 | # Perform the analysis: |

120 | algoMC.run() |

121 | |

122 | # Results: |

123 | result = algoMC.getResult() |

124 | probability = result.getProbabilityEstimate() |

125 | print "MonteCarlo result=" , result |

126 | print "Number of executed iterations =" , result.getOuterSampling() |

127 | print "Number of calls to the limit state =" , LimitState.getEvaluationCallsNumber() - initialNumberOfCall |

128 | print "Pf = " , probability |

129 | print "CV =" , result.getCoefficientOfVariation() |

130 | MCcvgraph = algoMC.drawProbabilityConvergence() |

131 | MCcvgraph.draw("montecarlo_convergence", 640, 480) |

132 | View(MCcvgraph) |

133 | |

134 | # Using FORM analysis |

135 | |

136 | # Resolution options: |

137 | eps=1E-3 |

138 | |

139 | # We create a NearestPoint algorithm |

140 | myCobyla = Cobyla() |

141 | myCobyla.setSpecificParameters(CobylaSpecificParameters()) |

142 | myCobyla.setMaximumIterationsNumber(100) |

143 | myCobyla.setMaximumAbsoluteError(eps) |

144 | myCobyla.setMaximumRelativeError(eps) |

145 | myCobyla.setMaximumResidualError(eps) |

146 | myCobyla.setMaximumConstraintError(eps) |

147 | |

148 | # For statistics about the algorithm |

149 | initialNumberOfCall = 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 | |

156 | algoFORM = FORM(NearestPointAlgorithm(myCobyla), myEvent, mean) |

157 | |

158 | # Perform the analysis: |

159 | algoFORM.run() |

160 | |

161 | # Results: |

162 | result = algoFORM.getResult() |

163 | print "Number of calls to the limit state =" , LimitState.getEvaluationCallsNumber() - initialNumberOfCall |

164 | print "Pf =" , result.getEventProbability() |

165 | |

166 | # Graphical result output |

167 | importanceFactorsGraph = result.drawImportanceFactors() |

168 | importanceFactorsGraph.draw("importance_factors", 640, 480) |

169 | View(importanceFactorsGraph) |

170 | |

171 | # Using Directional sampling |

172 | |

173 | # Resolution options: |

174 | cv = 0.05 |

175 | NbSim = 100000 |

176 | |

177 | algoDS = DirectionalSampling(myEvent) |

178 | algoDS.setMaximumOuterSampling(NbSim) |

179 | algoDS.setBlockSize(1) |

180 | algoDS.setMaximumCoefficientOfVariation(cv) |

181 | # For statistics about the algorithm |

182 | initialNumberOfCall = LimitState.getEvaluationCallsNumber() |

183 | |

184 | # Perform the analysis: |

185 | algoDS.run() |

186 | |

187 | # Results: |

188 | result = algoDS.getResult() |

189 | probability = result.getProbabilityEstimate() |

190 | print "Number of executed iterations =" , result.getOuterSampling() |

191 | print "Number of calls to the limit state =" , LimitState.getEvaluationCallsNumber() - initialNumberOfCall |

192 | print "Pf = " , probability |

193 | print "CV =" , result.getCoefficientOfVariation() |

194 | DScvgraph = algoDS.drawProbabilityConvergence() |

195 | DScvgraph.draw("directionalsampling_convergence", 640, 480) |

196 | View(DScvgraph) |

197 |