1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
"""
Create a mesh
=============
"""
# %%
import openturns as ot
import openturns.viewer as otv
import math as m
# %%
# Creation of a regular grid
# --------------------------
#
# In this first part we demonstrate how to create a regular grid. Note that a regular grid is a particular mesh of :math:`\mathcal{D}=[0,T] \in \mathbb{R}`.
#
# Here we assume it represents the time :math:`t` as it is often the case, but it can represent any physical quantity.
#
# A regular time grid is a regular discretization of the interval :math:`[0, T] \in \mathbb{R}` into :math:`N` points, noted :math:`(t_0, \dots, t_{N-1})`.
#
# The time grid can be defined using :math:`(t_{Min}, \Delta t, N)` where :math:`N` is the number of points in the time grid.
# :math:`\Delta t` the time step between two consecutive instants and :math:`t_0 = t_{Min}`.
# Then, :math:`t_k = t_{Min} + k \Delta t` and :math:`t_{Max} = t_{Min} + (N-1) \Delta t`.
#
#
# Consider :math:`X: \Omega \times \mathcal{D} \rightarrow \mathbb{R}^d` a multivariate stochastic process of dimension :math:`d`,
# where :math:`n=1`, :math:`\mathcal{D}=[0,T]` and :math:`t\in \mathcal{D}` is interpreted as a time stamp.
# Then the mesh associated to the process :math:`X` is a (regular) time grid.
# %%
# We define a time grid from a starting time :math:`tMin`, a time step :math:`tStep` and a number of time steps :math:`n`.
tMin = 0.0
tStep = 0.1
n = 10
time_grid = ot.RegularGrid(tMin, tStep, n)
# %%
# We get the first and the last instants, the step and the number of points :
start = time_grid.getStart()
step = time_grid.getStep()
grid_size = time_grid.getN()
end = time_grid.getEnd()
print("start=", start, "step=", step, "grid_size=", grid_size, "end=", end)
# %%
# We draw the grid.
time_grid.setName("time")
graph = time_grid.draw()
graph.setTitle("Time grid")
graph.setXTitle("t")
graph.setYTitle("")
view = otv.View(graph)
# %%
# Creation of a mesh
# ------------------
# %%
# In this paragraph we create a mesh :math:`\mathcal{M}` associated to a domain :math:`\mathcal{D} \in \mathbb{R}^n`.
#
# A mesh is defined from vertices in :math:`\mathbb{R}^n` and a topology that connects the vertices: the simplices.
# The simplex :math:`Indices([i_1,\dots, i_{n+1}])` relies the vertices of index :math:`(i_1,\dots, i_{n+1})` in :math:`\mathbb{N}^n`.
# In dimension 1, a simplex is an interval :math:`Indices([i_1,i_2])`; in dimension 2, it is a triangle :math:`Indices([i_1,i_2, i_3])`.
#
# The library enables to easily create a mesh which is a box of dimension :math:`d=1` or :math:`d=2` regularly meshed in all its directions, thanks to the object IntervalMesher.
#
# Consider :math:`X: \Omega \times \mathcal{D} \rightarrow \mathbb{R}^d` a multivariate stochastic process of dimension :math:`d`, where :math:`\mathcal{D} \in \mathbb{R}^n`.
# The mesh :math:`\mathcal{M}` is a discretization of the domain :math:`\mathcal{D}`.
# %%
# A one dimensional mesh is created and represented by :
vertices = [[0.5], [1.5], [2.1], [2.7]]
simplicies = [[0, 1], [1, 2], [2, 3]]
mesh1D = ot.Mesh(vertices, simplicies)
graph1 = mesh1D.draw()
graph1.setTitle("One dimensional mesh")
view = otv.View(graph1)
# %%
# We define a bidimensional mesh :
vertices = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [1.5, 1.0], [2.0, 1.5], [0.5, 1.5]]
simplicies = [[0, 1, 2], [1, 2, 3], [2, 3, 4], [2, 4, 5], [0, 2, 5]]
mesh2D = ot.Mesh(vertices, simplicies)
graph2 = mesh2D.draw()
graph2.setTitle("Bidimensional mesh")
graph2.setLegendPosition("lower right")
view = otv.View(graph2)
# %%
# We can also define a mesh which is a regularly meshed box in dimension 1 or 2.
# We define the number of intervals in each direction of the box :
myIndices = [5, 10]
myMesher = ot.IntervalMesher(myIndices)
# %%
# We then create the mesh of the box :math:`[0, 2] \times [0, 4]` :
lowerBound = [0.0, 0.0]
upperBound = [2.0, 4.0]
myInterval = ot.Interval(lowerBound, upperBound)
myMeshBox = myMesher.build(myInterval)
mygraph3 = myMeshBox.draw()
mygraph3.setTitle("Bidimensional mesh on a box")
view = otv.View(mygraph3)
# %%
# It is possible to perform a transformation on a regularly meshed box.
myIndices = [20, 20]
mesher = ot.IntervalMesher(myIndices)
# r in [1., 2.] and theta in (0., pi]
lowerBound2 = [1.0, 0.0]
upperBound2 = [2.0, m.pi]
myInterval = ot.Interval(lowerBound2, upperBound2)
meshBox2 = mesher.build(myInterval)
# %%
# We define the mapping function and draw the transformation :
f = ot.SymbolicFunction(["r", "theta"], ["r*cos(theta)", "r*sin(theta)"])
oldVertices = meshBox2.getVertices()
newVertices = f(oldVertices)
meshBox2.setVertices(newVertices)
graphMappedBox = meshBox2.draw()
graphMappedBox.setTitle("Mapped box mesh")
view = otv.View(graphMappedBox)
# %%
# Finally we create a mesh of a heart in dimension 2.
# sphinx_gallery_thumbnail_number = 6
def meshHeart(ntheta, nr):
# First, build the nodes
nodes = ot.Sample(0, 2)
nodes.add([0.0, 0.0])
for j in range(ntheta):
theta = (m.pi * j) / ntheta
if abs(theta - 0.5 * m.pi) < 1e-10:
rho = 2.0
elif (abs(theta) < 1e-10) or (abs(theta - m.pi) < 1e-10):
rho = 0.0
else:
absTanTheta = abs(m.tan(theta))
rho = absTanTheta ** (1.0 / absTanTheta) + m.sin(theta)
cosTheta = m.cos(theta)
sinTheta = m.sin(theta)
for k in range(nr):
tau = (k + 1.0) / nr
r = rho * tau
nodes.add([r * cosTheta, r * sinTheta - tau])
# Second, build the triangles
triangles = []
# First heart
for j in range(ntheta):
triangles.append([0, 1 + j * nr, 1 + ((j + 1) % ntheta) * nr])
# Other hearts
for j in range(ntheta):
for k in range(nr - 1):
i0 = k + 1 + j * nr
i1 = k + 2 + j * nr
i2 = k + 2 + ((j + 1) % ntheta) * nr
i3 = k + 1 + ((j + 1) % ntheta) * nr
triangles.append([i0, i1, i2 % (nr * ntheta)])
triangles.append([i0, i2, i3 % (nr * ntheta)])
return ot.Mesh(nodes, triangles)
mesh4 = meshHeart(48, 16)
graphMesh = mesh4.draw()
graphMesh.setTitle("Bidimensional mesh")
graphMesh.setLegendPosition("")
view = otv.View(graphMesh)
# %%
# Display all figures
otv.View.ShowAll()
|
