-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathdefine_mapping.py
112 lines (101 loc) · 4.33 KB
/
define_mapping.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
###############################################################################
import numpy as np
import math as math
###############################################################################
# Q1 nodes for mapping are corners of Q2 basis functions
# Q2 nodes for mapping are same as Q2 basis functions
# Q3,Q4,Q5,Q6 nodes for mapping are built
# note that python uses row-major storage for 2S arrays and since we need
# to do dot products with xmapping and zmapping it makes more sense
# to use column-major, i.e. F_CONTIGUOUS in python jargon.
###############################################################################
def define_mapping(mapping,mmapping,xV,zV,iconV,nel,axisymmetric,rad,theta,nelt):
xmapping=np.zeros((mmapping,nel),dtype=np.float64,order='F')
zmapping=np.zeros((mmapping,nel),dtype=np.float64,order='F')
if mapping=='Q1':
for iel in range(0,nel):
xmapping[0,iel]=xV[iconV[0,iel]] ; zmapping[0,iel]=zV[iconV[0,iel]]
xmapping[1,iel]=xV[iconV[1,iel]] ; zmapping[1,iel]=zV[iconV[1,iel]]
xmapping[2,iel]=xV[iconV[2,iel]] ; zmapping[2,iel]=zV[iconV[2,iel]]
xmapping[3,iel]=xV[iconV[3,iel]] ; zmapping[3,iel]=zV[iconV[3,iel]]
if mapping=='Q2':
for iel in range(0,nel):
xmapping[0,iel]=xV[iconV[0,iel]] ; zmapping[0,iel]=zV[iconV[0,iel]]
xmapping[1,iel]=xV[iconV[1,iel]] ; zmapping[1,iel]=zV[iconV[1,iel]]
xmapping[2,iel]=xV[iconV[2,iel]] ; zmapping[2,iel]=zV[iconV[2,iel]]
xmapping[3,iel]=xV[iconV[3,iel]] ; zmapping[3,iel]=zV[iconV[3,iel]]
xmapping[4,iel]=xV[iconV[4,iel]] ; zmapping[4,iel]=zV[iconV[4,iel]]
xmapping[5,iel]=xV[iconV[5,iel]] ; zmapping[5,iel]=zV[iconV[5,iel]]
xmapping[6,iel]=xV[iconV[6,iel]] ; zmapping[6,iel]=zV[iconV[6,iel]]
xmapping[7,iel]=xV[iconV[7,iel]] ; zmapping[7,iel]=zV[iconV[7,iel]]
xmapping[8,iel]=xV[iconV[8,iel]] ; zmapping[8,iel]=zV[iconV[8,iel]]
if mapping=='Q3':
if not axisymmetric:
dtheta=2*np.pi/nelt/3
else:
dtheta=np.pi/nelt/3
for iel in range(0,nel):
thetamin=theta[iconV[0,iel]]
rmin=rad[iconV[0,iel]]
rmax=rad[iconV[2,iel]]
counter=0
for j in range(0,4):
for i in range(0,4):
ttt=thetamin+i*dtheta
rrr=rmin+j*(rmax-rmin)/3
xmapping[counter,iel]=np.sin(ttt)*rrr
zmapping[counter,iel]=np.cos(ttt)*rrr
counter+=1
if mapping=='Q4':
if not axisymmetric:
dtheta=2*np.pi/nelt/4
else:
dtheta=np.pi/nelt/4
for iel in range(0,nel):
thetamin=theta[iconV[0,iel]]
rmin=rad[iconV[0,iel]]
rmax=rad[iconV[2,iel]]
counter=0
for j in range(0,5):
for i in range(0,5):
ttt=thetamin+i*dtheta
rrr=rmin+j*(rmax-rmin)/4
xmapping[counter,iel]=np.sin(ttt)*rrr
zmapping[counter,iel]=np.cos(ttt)*rrr
counter+=1
if mapping=='Q5':
if not axisymmetric:
dtheta=2*np.pi/nelt/5
else:
dtheta=np.pi/nelt/5
for iel in range(0,nel):
thetamin=theta[iconV[0,iel]]
rmin=rad[iconV[0,iel]]
rmax=rad[iconV[2,iel]]
counter=0
for j in range(0,6):
for i in range(0,6):
ttt=thetamin+i*dtheta
rrr=rmin+j*(rmax-rmin)/5
xmapping[counter,iel]=np.sin(ttt)*rrr
zmapping[counter,iel]=np.cos(ttt)*rrr
counter+=1
if mapping=='Q6':
if not axisymmetric:
dtheta=2*np.pi/nelt/6
else:
dtheta=np.pi/nelt/6
for iel in range(0,nel):
thetamin=theta[iconV[0,iel]]
rmin=rad[iconV[0,iel]]
rmax=rad[iconV[2,iel]]
counter=0
for j in range(0,7):
for i in range(0,7):
ttt=thetamin+i*dtheta
rrr=rmin+j*(rmax-rmin)/6
xmapping[counter,iel]=np.sin(ttt)*rrr
zmapping[counter,iel]=np.cos(ttt)*rrr
counter+=1
return xmapping,zmapping
###############################################################################