-
Notifications
You must be signed in to change notification settings - Fork 0
/
get_plate.py
148 lines (122 loc) · 3.91 KB
/
get_plate.py
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
# -*- coding: utf-8 -*-
# Description:
# Generate a mesh triangulation of the reference square (a,b)^2.
#
# Arguments:
# n Number of nodes in each spatial direction (n^2 total nodes).
# a Lower limit for x and y
# b Upper limit for x and y
#
# Returns:
# p Nodal points, (x,y)-coordinates for point i given in row i.
# tri Elements. Index to the three corners of element i given in row i.
# edge Index list of all nodal points on the outer edge.
#
# Written by Olav M. S. Gran using the
# Boiler code by: Kjetil a. Johannessen, Abdullah Abdulhaque (October 2019)
# from https://wiki.math.ntnu.no/tma4220/2020h/project : getplate.py
import numpy as np
import scipy.spatial as spsa
def _make_edge(n):
"""
Make the edge
Parameters
----------
n : int
Number of nodes in each spatial direction (n^2 total nodes).
Returns
-------
edge : np.array
Index list of all nodal points on the outer edge.
"""
n2 = n * n
# Generating nodal points on outer _edge.
south_edge = np.array([np.arange(1, n), np.arange(2, n + 1)]).T
east_edge = np.array([np.arange(n, n2 - n + 1, n), np.arange(2 * n, n2 + 1, n)]).T
north_edge = np.array([np.arange(n2, n2 - n + 1, -1), np.arange(n2 - 1, n2 - n, -1)]).T
west_edge = np.array([np.arange(n2 - n + 1, n - 1, -n), np.arange(n2 - 2 * n + 1, 0, -n)]).T
edge = np.vstack((south_edge, east_edge, north_edge, west_edge))
# Added this to get this script too work.
edge -= 1
return edge
def getPlatev2(n, a=0, b=1):
"""
Get the plate (a,b)^2, version 2
Parameters
----------
n : int
Number of nodes in each spatial direction (n^2 total nodes).
a : float, optional
Lower limit for x and y. The default is 0.
b : float, optional
Upper limit for x and y. The default is 1.
Returns
-------
p : np.array
Nodal points, (x,y)-coordinates for point i given in row i.
tri : np.array
Elements. Index to the three corners of element i given in row i.
edge : np.array
Index list of all nodal points on the outer edge.
"""
# Defining auxiliary variables.
l = np.linspace(a, b, n)
y, x = np.meshgrid(l, l)
# Generating nodal points.
n2 = n * n
p = np.zeros((n2, 2))
p[:, 0] = x.T.ravel()
p[:, 1] = y.T.ravel()
# Generating delaunay elements.
mesh = spsa.Delaunay(p)
tri = mesh.simplices
edge = _make_edge(n)
return p, tri, edge
def getPlatev3(n, a=0, b=1):
"""
Get the plate (a,b)^2, version 3
Parameters
----------
n : int
Number of nodes in each spatial direction (n^2 total nodes).
a : float, optional
Lower limit for x and y. The default is 0.
b : float, optional
Upper limit for x and y. The default is 1.
Returns
-------
p : np.array
Nodal points, (x,y)-coordinates for point i given in row i.
tri : np.array
Elements. Index to the three corners of element i given in row i.
edge : np.array
Index list of all nodal points on the outer edge.
"""
# Defining auxiliary variables.
l = np.linspace(a, b, n)
y, x = np.meshgrid(l, l)
# Generating nodal points.
n2 = n * n
p = np.zeros((n2, 2))
p[:, 0] = x.T.ravel()
p[:, 1] = y.T.ravel()
# Generating elements.
n12 = (n - 1) * (n - 1) * 2
tri = np.zeros((n12, 3), dtype=int)
def index_map(i, j):
return i + n * j
k = 0
for i in range(n - 1):
for j in range(n - 1):
tri[k, 0] = index_map(i, j)
tri[k, 1] = index_map(i + 1, j)
tri[k, 2] = index_map(i + 1, j + 1)
k += 1
tri[k, 0] = index_map(i, j)
tri[k, 1] = index_map(i + 1, j + 1)
tri[k, 2] = index_map(i, j + 1)
k += 1
arg = np.argsort(tri[:, 0])
tri = tri[arg]
edge = _make_edge(n)
return p, tri, edge