forked from eX-Mech/pipeMeshNek
-
Notifications
You must be signed in to change notification settings - Fork 0
/
plotMesh.py
116 lines (108 loc) · 2.54 KB
/
plotMesh.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
#!/usr/bin/env python
#
# load and plot 2D mesh created by pipeMeshNek
#------------------------------------------------------------------------------
# import pymech
#
import os
try:
assert os.path.exists('./pymech')
except AssertionError:
import subprocess
subprocess.call('git clone https://github.com/jcanton/pymech.git', shell=True)
import sys
sys.path.append('./pymech/src/')
import neksuite as ns
import numpy as np
import matplotlib.pyplot as plt
#------------------------------------------------------------------------------
# plotting function
#
def arc_2points(x, y, r):
#
# plot arc given 2 points and radius
#
import numpy as np
#
x1 = x[0]; x2 = x[1]
y1 = y[0]; y2 = y[1]
#
alpha = (x1 - x2)/(y2 - y1)
beta = (x2**2 - x1**2 + y2**2 - y1**2)/(2*(y2 - y1))
a = 1 + alpha**2
b = -2*x1 - 2*y1*alpha + 2*alpha*beta
c = x1**2 + y1**2 - 2*y1*beta + beta**2 - r**2
#
delta = b**2 - 4*a*c
#
xc1 = (-b + np.sqrt(delta))/(2*a)
xc2 = (-b - np.sqrt(delta))/(2*a)
#
yc1 = alpha*xc1 + beta
yc2 = alpha*xc2 + beta
#
d1 = np.sqrt(xc1**2 + yc1**2)
d2 = np.sqrt(xc2**2 + yc2**2)
#
if ( d1<d2 ):
xc = xc1; yc = yc1
else:
xc = xc2; yc = yc2
#
tt1 = np.arctan2( (y1-yc), (x1-xc) )
tt2 = np.arctan2( (y2-yc), (x2-xc) )
#
if (abs(tt1)<1e-3):
tt1 = 0
if (abs(tt2)<1e-3):
tt2 = 0
#
if (tt1 < 0):
tt1 = tt1 + 2*np.pi
if (tt2 < 0):
tt2 = tt2 + 2*np.pi
#
start_t = min(tt1,tt2)
end_t = max(tt1,tt2)
#
if (start_t==0 and end_t>3./2.*np.pi):
start_t = end_t
end_t = 2*np.pi
#
tt = np.linspace(start_t, end_t, 100)
xx = xc + np.abs(r)*np.cos(tt)
yy = yc + np.abs(r)*np.sin(tt)
#
plt.plot(xx, yy, '-k')
#------------------------------------------------------------------------------
# load mesh
#
fname = 'base2d.rea'
field = ns.readrea(fname)
#------------------------------------------------------------------------------
# plot
#
plt.figure(1)
plt.clf()
nedges = 4
for iel in range(field.nel):
xv = np.reshape(field.elem[iel].pos[0, 0,:,:], (4,1) )
yv = np.reshape(field.elem[iel].pos[1, 0,:,:], (4,1) )
xl = (np.min(xv)+np.max(xv))/2
yl = (np.min(yv)+np.max(yv))/2
plt.text(xl, yl, '%d' % (iel+1), horizontalalignment='center', verticalalignment='center')
for iedge in range(nedges):
xe = np.roll(xv, -iedge)[0:2]
ye = np.roll(yv, -iedge)[0:2]
if (field.elem[iel].curv[iedge] == 0):
plt.plot(xe, ye, '-k')
else:
arc_2points(xe, ye, field.elem[iel].curv[iedge])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'2D section of the mesh')
plt.grid(True)
plt.axis('equal')
plt.draw()
plt.show(block=True)