Coverage for /usr/local/lib/python3.8/site-packages/spam/mesh/tetrahedra.py: 76%
88 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1# Library of SPAM functions for dealing with assembly maxtrix of unstructured 3D meshes made of tetrahedra
2# Copyright (C) 2020 SPAM Contributors
3#
4# This program is free software: you can redistribute it and/or modify it
5# under the terms of the GNU General Public License as published by the Free
6# Software Foundation, either version 3 of the License, or (at your option)
7# any later version.
8#
9# This program is distributed in the hope that it will be useful, but WITHOUT
10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
12# more details.
13#
14# You should have received a copy of the GNU General Public License along with
15# this program. If not, see <http://www.gnu.org/licenses/>.
16import numpy
19def tetCentroid(points):
20 """Compute coordinates of the centroid of a tetrahedron
22 Parameters
23 ----------
24 points : 4x3 array
25 Array of the 3D coordinates of the 4 points
27 Returns
28 -------
29 3x1 array
30 the coordinates of the centroid
31 """
33 centroid = numpy.zeros(3)
34 for i in range(3):
35 a = points[0][i]
36 b = points[1][i]
37 c = points[2][i]
38 d = points[3][i]
40 # mid point of AB (1/2 of A -> B)
41 mid_ab = a + (b - a) / 2.0
43 # centroid of ABC (1/3 of AB -> C median)
44 face_centroid = mid_ab + (c - mid_ab) / 3.0
46 # centroid of ABCD (1/4 of ABC -> D median)
47 centroid[i] = face_centroid + (d - face_centroid) / 4.0
49 return centroid
52def tetVolume(points):
53 """Compute the volume of a tetrahedron
55 Parameters
56 ----------
57 points : 4x3 array
58 Array of the 3D coordinates of the 4 points
60 Returns
61 -------
62 float
63 the volume of the tetrahedron
64 """
66 # compute the jacobian matrix by padding the points with a line of 1
67 pad = numpy.array([[1], [1], [1], [1]])
68 jacobian = numpy.hstack([pad, points])
70 # the volume is 1/6 of the determinant of the jacobian matrix
71 return numpy.linalg.det(jacobian) / 6.0
74def shapeFunctions(points):
75 """
76 This function computes the shape coefficients matrux from the coordinates
77 of the 4 nodes of a linear tetrahedron (see 4.11 in Zienkiewicz).
79 .. code-block:: text
81 coefficients_matrix = [
82 a1 b1 c1 d1
83 a2 b2 c2 d2
84 a3 b3 c3 d3
85 a4 b4 c4 d4
86 ]
87 with
88 N1 = a1 + b1x + c1y + d1z
89 N2 = a2 + b2x + c2y + d2z
90 N3 = a3 + b3x + c3y + d3z
91 N4 = a4 + b4x + c4y + d4z
93 Parameters
94 ----------
95 points : 4x3 array
96 Array of the 3D coordinates of the 4 points
98 Returns
99 -------
100 volume : float
101 The volume of the tetrahedron
102 coefficients_matrix : 4x4 array
103 The coefficients 4 coefficients of the 4 shape functions
105 Note
106 -----
107 Pure python function.
108 """
109 # cast into numpy array and check shape
110 points = numpy.array(points)
111 if points.shape != (4, 3):
112 raise ValueError("Points coordinates in bad format. Can I have a 4x3 please?")
114 # jacobian matrix
115 pad = numpy.array([[1], [1], [1], [1]])
116 jacobian = numpy.hstack([pad, points])
117 six_v = numpy.linalg.det(jacobian)
118 volume = six_v / 6.0
120 coefficients_matrix = numpy.zeros((4, 4))
121 for l_number in range(4): # loop over the 4 nodes
122 for c_number in range(4): # loop over the 4 coefficients for 1 node
123 # create mask for the sub matrix (ie: [1, 2, 3], [0, 2, 3], ...)
124 lines = [_ for _ in range(4) if _ != l_number]
125 columns = [_ for _ in range(4) if _ != c_number]
126 # creates the sub_jacobian (needs two runs for some numpy reason)
127 sub_jacobian = jacobian[:, columns]
128 sub_jacobian = sub_jacobian[lines, :]
130 # compute the determinant and fill coefficients matrix
131 det = numpy.linalg.det(sub_jacobian)
132 sign = (-1.0) ** (l_number + c_number)
133 coefficients_matrix[l_number, c_number] = sign * det / six_v
135 return volume, coefficients_matrix
138def elementaryStiffnessMatrix(points, young, poisson):
139 """
140 This function computes elementary stiffness matrix from the coordinates
141 of the 4 nodes of a linear tetrahedron.
143 .. code-block:: text
145 D = [something with E and mu] (6x6)
146 B = [B1, B2, B3, B4] (4 times 6x3 -> 6x12)
147 ke = V.BT.D.B (12x12)
149 Parameters
150 ----------
151 points : 4x3 array
152 Array of the 3D coordinates of the 4 points
153 young: float
154 Young modulus
155 poisson: float
156 Poisson ratio
158 Returns
159 -------
160 stiffness_matrix : 12x12 array
161 The stiffness matrix of the tetrahedron
163 Note
164 -----
165 Pure python function.
166 """
168 # get volume and shape functions coefficients
169 volume, N = shapeFunctions(points)
171 # compute the full B matrix
172 def Ba(N, a):
173 return numpy.array(
174 [
175 [N[a, 1], 0, 0],
176 [0, N[a, 2], 0],
177 [0, 0, N[a, 3]],
178 [0, N[a, 3], N[a, 2]],
179 [N[a, 3], 0, N[a, 1]],
180 [N[a, 2], N[a, 1], 0],
181 ]
182 )
184 # initialise B with first Ba
185 B = Ba(N, 0)
186 for a in [1, 2, 3]:
187 B = numpy.hstack([B, Ba(N, a)])
189 # compute D matrix
190 D = (1.0 - poisson) * numpy.eye(6)
191 D[0, [1, 2]] = poisson
192 D[1, [0, 2]] = poisson
193 D[2, [0, 1]] = poisson
194 D[3:, 3:] -= 0.5 * numpy.eye(3)
195 D *= young / ((1 + poisson) * (1 - 2 * poisson))
197 # compute stiffness matrix
198 ke = volume * numpy.matmul(numpy.transpose(B), numpy.matmul(D, B))
200 return ke
203def globalStiffnessMatrix(points, connectivity, young, poisson):
204 """
205 This function assembles elementary stiffness matrices to computes an elastic
206 (3 dof) global stiffness matrix for a 4 noded tetrahedra mesh
208 Notations
209 ---------
210 neq: number of global equations
211 nel: number of elements
212 npo: number of points
213 neq = 3 x npo
215 Parameters
216 ----------
217 points : npo x 3 array
218 Array of the 3D coordinates of the all nodes points
219 connectivity : nel x 4 array
220 Connectivity matrix
221 young: float
222 Young modulus
223 poisson: float
224 Poisson ratio
226 Returns
227 -------
228 stiffness_matrix : neq x neq array
229 The stiffness matrix of the tetrahedron
231 Note
232 -----
233 Pure python function.
234 """
235 # cast to numpy array
236 points = numpy.array(points)
237 connectivity = numpy.array(connectivity)
239 neq = 3 * points.shape[0]
240 # nel = connectivity.shape[0]
241 K = numpy.zeros((neq, neq), dtype=float)
243 for e, element_nodes_number in enumerate(connectivity):
244 # print(f'Element number {e}: {element_nodes_number}')
246 # built local stiffness matrix (12x12)
247 element_points = [points[A] for A in element_nodes_number]
248 ke = elementaryStiffnessMatrix(element_points, young, poisson)
250 # loop over stiffness matrix (lines and rows)
251 for p1 in range(12):
252 i1 = p1 % 3 # deduce degree of freedom
253 a1 = (p1 - i1) // 3 # deduce local node number
254 P1 = 3 * int(connectivity[e, a1]) + i1 # get global equation number
255 for p2 in range(p1, 12):
256 i2 = p2 % 3
257 a2 = (p2 - i2) // 3
258 P2 = 3 * int(connectivity[e, a2]) + i2
259 K[P1, P2] += ke[p1, p2]
261 K = K + K.transpose()
262 K[numpy.diag_indices(neq)] *= 0.5
264 return K
267if __name__ == "__main__":
268 import matplotlib.pyplot as plt
269 import spam.mesh
271 # points = [
272 # [0, 0, 0],
273 # [0, 1, 0],
274 # [0, 0, 1],
275 # [1, 0, 0]
276 # ]
277 #
278 # coefficients_matrix = shapeFunctions(points)
279 # ke = elementaryStiffnessMatrix(points, 1, 0.3)
280 # print(f'Elementary connectivity matrix: {ke.shape}')
281 # # print(ke)
283 lengths = [1, 1, 1]
284 lc = 0.2
285 points, connectivity = spam.mesh.createCuboid(
286 lengths,
287 lc,
288 origin=[0.0, 0.0, 0.0],
289 periodicity=False,
290 verbosity=1,
291 gmshFile=None,
292 vtkFile=None,
293 binary=False,
294 skipOutput=False,
295 )
297 K = globalStiffnessMatrix(points, connectivity, 1, 0.3)
299 plt.imshow(numpy.log(K[:100, :100]))
300 plt.show()
302 # ID function
303 def ID(i: int, A: int):
304 """
305 Defined as in Hughes: The Finite Element Method eq. 2.8.3
306 WARNING: without the condition of dirichlet boundary conditions
307 Takes
308 - the degree of freedom: i = 0, 1 or 2
309 - the node global number: A = 0, 1, 2, ... n_nodes
310 Returns
311 - the global equation number: P = 0, 1, 2 ... 3 * n_nodes
312 """
313 P = 3 * A + i
314 return int(P)
316 # for A, point in enumerate(points):
317 # for i in range(3):
318 # print(f'Node number {A} dof {i}: {ID(i, A)}')
320 # IEN function
321 def IEN(a: int, e: int, connectivity):
322 """
323 Defined as in Hughes: The Finite Element Method eq. 2.6.1
324 Takes
325 - the local node number: a = 0, 1, 2 or 3
326 - the element number: e = 0, 1, 2, ... n_elem
327 Returns
328 - the node global number: A = 0, 1, 2, ... n_nodes
329 Uses the connectivity matrix
330 """
331 # connectivity matrix shape: n_elem x 4
332 A = connectivity[e, a]
333 return int(A)
335 # for e, _ in enumerate(connectivity):
336 # for a in range(4):
337 # print(f'Element number {e} node number {a}: {IEN(a, e, connectivity)}')
339 # LM function
340 def LM(i: int, a: int, e: int, connectivity):
341 """
342 Defined as in Hughes: The Finite Element Method eq. 2.10.1
343 Takes
344 - the degree of freedom: i = 0, 1 or 2
345 - the local node number: a = 0, 1, 2 or 3
346 - the element number: e = 0, 1, 2, ... n_elem
347 Returns
348 - the global equation number: P = 0, 1, 2 ... 3 * n_nodes
349 Uses the connectivity matrix
350 """
351 P = ID(i, IEN(a, e, connectivity))
352 return P
354 # for e, _ in enumerate(connectivity):
355 # for a in range(4):
356 # for i in range(3):
357 # P = 3 * int(connectivity[e, a]) + i
358 # print(f'Element number {e} node number {a} dof {i}: {LM(i, a, e, connectivity)}')
359 #
360 #
361 # if e == 3:
362 # break