forked from cedrict/fieldstone
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathapp_useful_python.tex
More file actions
180 lines (126 loc) · 4.72 KB
/
Copy pathapp_useful_python.tex
File metadata and controls
180 lines (126 loc) · 4.72 KB
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
\begin{flushright} {\tiny {\color{gray} app\_useful\_python.tex}} \end{flushright}
%--------------------------
\subsection{Sparse matrices}
So far, the best way I have found to deal with sparse matrices is to
declare the matrix as a 'lil\_matrix' (linked list).
\begin{lstlisting}
from scipy.sparse import csr_matrix, lil_matrix
A_mat = lil_matrix((Nfem,Nfem),dtype=np.float64)
\end{lstlisting}
One then adds terms to it as if it was a full/dense matrix.
Once the assembly is done, the conversion to CSR format is trivial:
\begin{lstlisting}
A_mat=A_mat.tocsr()
\end{lstlisting}
Finally the solver can be called:
\begin{lstlisting}
sol=sps.linalg.spsolve(A_mat,rhs)
\end{lstlisting}
%--------------------------
\subsection{condition number}
if the matrix has been declared as lil\_matrix, first convert it to a dense matrix:
\begin{lstlisting}
A_mat=A_mat.dense()
\end{lstlisting}
The condition number of the matrix is simply obtained as follows:
\begin{lstlisting}
from numpy import linalg as LA
print(LA.cond(A_mat))
\end{lstlisting}
%--------------------------
\subsection{Weird stuff}
Python is touted as the one language students should learn and master. However it is a language which allows *way* too much liberty in its syntax and encourages students to be sloppy.
For instance the following code runs just fine:
\begin{lstlisting}
for k in range(0,5):
for k in range(0,5):
for k in range(0,5):
print (k)
\end{lstlisting}
This alone should disqualify this language. It is easy to see the obvious problem with this code, but adding a few lines of code in between each 'for' line hides the problem and the absence of any warning makes this code a nightmare to debug.
\begin{verbatim}
https://www.w3schools.com/python/default.asp
https://www.codecademy.com/learn/learn-python
https://learnpythonthehardway.org/book/
\end{verbatim}
%----------------------------------
\subsection{Making simple 2D plots}
\begin{lstlisting}
import matplotlib.pyplot as plt
# number of points
N=100
# a despicable way of filling two arrays
x_data=[]
y_data=[]
for i in range(0,N):
x=i
y=i**2+2*i+1.
x_data.append(x)
y_data.append(y)
# generating a 2D figure with the data
plt.figure()
plt.plot(x_data,y_data, label = 'name of data')
plt.xlabel('x-axis label')
plt.ylabel('y-axis label')
plt.legend()
plt.savefig('myplot.pdf', bbox_inches='tight')
plt.show()
\end{lstlisting}
\begin{center}
\includegraphics[width=7cm]{images/python/myplot}
\end{center}
%----------------------------------
\subsection{Making simple 3D plots of scatter}
\begin{lstlisting}
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_title("insert here text for title")
size = ..some value..
ax.scatter3D(x, y, z, s = size)
\end{lstlisting}
%----------------------------------
\subsection{How to debug your code}
Debugging a FE code is by no means trivial. There is (at least) a grid, a connectivity array, basis functions and their derivatives, the elemental matrices and rhs, the assembly, boundary conditions, and a call to a solver before the solution (if the solver returns one!) can be visualised.
\begin{itemize}
\item First and foremost, make sure that your grid of points is correct.
For instance, you can resort to exporting it to an ascii file as follows:
\begin{lstlisting}
np.savetxt('velocity.ascii',np.array([x,y,u,v]).T,header='# x,y,u,v')
\end{lstlisting}
In two dimensions, you should set for example nelx=3 and nely=2, so that
for $Q_1$ elements the grid counts 12 points. Then make sure the coordinates and the order of the points makes sense.
Repeat the process for pressure nodes, temperature nodes, etc ...
\item Then it is time to check the connectivity array(s).
\begin{lstlisting}
for iel in range (0,nel):
print ("iel=",iel)
for k in range(0,m)
print ("node ",icon[0,iel],"at pos.",x[icon[0,iel]], y[icon[0,iel]])
\end{lstlisting}
This displays the list of nodes and their positions making each element.
Repeat the process for every connectivity array.
\item We can go on with testing that the all basis functions are 1 on their node and zero elsewhere:
\begin{lstlisting}
for i in range(0,m):
print ('node',i,':',NNV(rnodes[i],snodes[i]))
\end{lstlisting}
here the arrays rnodes and snodes contain the (r,s) coordinates of the m nodes
\item test jacobian, compute volume of domain
\item sum(dNNNdx)=0
\item print nodes where bc
\end{itemize}
%----------------------------------
\subsection{Optional arguments}
Courtesy of Henry Brett.
\begin{lstlisting}
def myfunc(a,b, *optional_arguments, **keyword_arguments):
print(a)
print(b)
for ar in optional_arguments:
print(ar)
d=keyword_arguments.get("d", None)
print(d)
a="dog"
b="cat"
myfunc(a,b,"shrek","fiona",d="donkey")
\end{lstlisting}