# ExampleAxialStressedBeam: AxialStressedBeam.py

File AxialStressedBeam.py, 5.7 KB (added by regis.lebrun@…, 8 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 | # New syntax |

6 | from openturns.viewer import ViewImage |

7 | |

8 | # Constant definition: |

9 | stochasticDimension = 2 |

10 | |

11 | # Random generator initialization: |

12 | RandomGenerator().SetSeed(0) |

13 | |

14 | #Analytical model definition: |

15 | |

16 | # Analytical construction : Input |

17 | inputFunction = Description(stochasticDimension) |

18 | inputFunction[0] = "R" |

19 | inputFunction[1] = "F" |

20 | |

21 | # Analytical construction : Output |

22 | outputFunction = Description(1) |

23 | outputFunction[0] = "G" |

24 | |

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

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

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

28 | |

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

30 | |

31 | # Test of the limit state function: |

32 | x = NumericalPoint(stochasticDimension, 0) |

33 | x[0]=300.0 |

34 | x[1]=75000.0 |

35 | print "x=" , x |

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

37 | |

38 | # Stochastic model definition: |

39 | |

40 | # Mean |

41 | mean = NumericalPoint(stochasticDimension, 0.0) |

42 | mean[0] = 300.0 |

43 | mean[1] = 75000.0 |

44 | |

45 | # Standard deviation |

46 | sigma = NumericalPoint(stochasticDimension, 0.0) |

47 | sigma[0] = 30.0 |

48 | sigma[1] = 5000.0 |

49 | |

50 | # Additional parameters for the lognormal distribution: |

51 | component = Description(1) |

52 | BorneInf = 0.0 |

53 | |

54 | # Initialization of the distribution collection: |

55 | aCollection = DistributionCollection() |

56 | |

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

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

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

60 | component[0] = "R" |

61 | marginal.setDescription(component) |

62 | # Graphical output of the PDF |

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

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

67 | # Visualize the graph |

68 | ViewImage(pdfgraph1.getBitmap()) |

69 | # Fill the first marginal of aCollection |

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

71 | |

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

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

74 | marginal.setName("Traction_load") |

75 | component[0] = "F" |

76 | marginal.setDescription(component) |

77 | # Graphical output of the PDF |

78 | pdfgraph2=marginal.drawPDF() |

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

80 | ViewImage(pdfgraph2.getBitmap()) |

81 | # Fill the second marginal of aCollection |

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

83 | |

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

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

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

87 | |

88 | # Instanciate one distribution object |

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

90 | myDistribution.setName("myDist") |

91 | |

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

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

94 | |

95 | # We create a composite random vector |

96 | G = RandomVector(LimitState, vect) |

97 | |

98 | # We create an Event from this RandomVector |

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

100 | |

101 | # Using Monte Carlo simulations |

102 | |

103 | # Resolution options: |

104 | cv = 0.05 |

105 | NbSim = 100000 |

106 | |

107 | algoMC = MonteCarlo(myEvent) |

108 | algoMC.setMaximumOuterSampling(NbSim) |

109 | algoMC.setBlockSize(1) |

110 | algoMC.setMaximumCoefficientOfVariation(cv) |

111 | # For statistics about the algorithm |

112 | initialNumberOfCall = LimitState.getEvaluationCallsNumber() |

113 | |

114 | # Perform the analysis: |

115 | algoMC.run() |

116 | |

117 | # Results: |

118 | result = algoMC.getResult() |

119 | probability = result.getProbabilityEstimate() |

120 | print "MonteCarlo result=" , result |

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

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

123 | print "Pf = " , probability |

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

125 | MCcvgraph = algoMC.drawProbabilityConvergence() |

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

127 | ViewImage(MCcvgraph.getBitmap()) |

128 | |

129 | # Using FORM analysis |

130 | |

131 | # Resolution options: |

132 | eps=1E-3 |

133 | |

134 | # We create a NearestPoint algorithm |

135 | myCobyla = Cobyla() |

136 | myCobyla.setSpecificParameters(CobylaSpecificParameters()) |

137 | myCobyla.setMaximumIterationsNumber(100) |

138 | myCobyla.setMaximumAbsoluteError(eps) |

139 | myCobyla.setMaximumRelativeError(eps) |

140 | myCobyla.setMaximumResidualError(eps) |

141 | myCobyla.setMaximumConstraintError(eps) |

142 | |

143 | # For statistics about the algorithm |

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

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

152 | |

153 | # Perform the analysis: |

154 | algoFORM.run() |

155 | |

156 | # Results: |

157 | result = algoFORM.getResult() |

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

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

160 | |

161 | # Graphical result output |

162 | importanceFactorsGraph = result.drawImportanceFactors() |

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

164 | ViewImage(importanceFactorsGraph.getBitmap()) |

165 | |

166 | # Using Directional sampling |

167 | |

168 | # Resolution options: |

169 | cv = 0.05 |

170 | NbSim = 100000 |

171 | |

172 | algoDS = DirectionalSampling(myEvent) |

173 | algoDS.setMaximumOuterSampling(NbSim) |

174 | algoDS.setBlockSize(1) |

175 | algoDS.setMaximumCoefficientOfVariation(cv) |

176 | # For statistics about the algorithm |

177 | initialNumberOfCall = LimitState.getEvaluationCallsNumber() |

178 | |

179 | # Perform the analysis: |

180 | algoDS.run() |

181 | |

182 | # Results: |

183 | result = algoDS.getResult() |

184 | probability = result.getProbabilityEstimate() |

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

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

187 | print "Pf = " , probability |

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

189 | DScvgraph = algoDS.drawProbabilityConvergence() |

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

191 | ViewImage(DScvgraph.getBitmap()) |

192 |