• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Apr 25th, 2019 92 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
2. import numpy as np
3. import mpmath
4. import sympy
5. from sympy import Point3D, Line3D, Plane, Segment3D, Symbol
6. from sympy.geometry import intersection, Segment3D
7. import re
8. from sympy.solvers import solve
10. import h5py
11. import sys
12. import math
13.
14. #CONVERT FUNCTION
15. #this handles turning values returned as "a/b" to decimal format
16. #this occurs for some intersection points
17. #takes a value as input
18. def convert(z):
19.     try:
20.             return float(z)
21.     except ValueError:
22.             num, denom = z.split('/')
23.             return float(num) / float(denom)
24.
25. #INTERSECTION FUNCTION
26. #finds intersections between plane and line entities
27. #takes a plane, line, and storage list as input
28. def find_intersection(a,b,XYZ):
29. #loop over the geometry entities and find intersections
30.     for i in range(6):
31.         for j in range(12):
32.             s=intersection(a[i],b[j])#get intersection
33.             #sort through the intersection data and retrieve numbers from the output string
34.             if s and (b[j].contains(s[0])==True):#if intersection exists and point is on edge
35.                 d=str(s[0])
36.                 l=len(d)
37.                 r1=[]
38.                 w1=[]
39.                 if(type(s[0])==sympy.geometry.line.Segment3D):#if intersection is line segment type
40.                     e=d[18:(l-1)]
41.                     q=re.search('P(.*?)\)', e)
42.                     l=len(q.group(1))
43.                     u=q.group(1)[7:l]
44.                     if(u not in XYZ):
45.                         XYZ.append(u)
46.
47.                     w=u.split(",")
48.                     f=re.search('(.*?)\)', e)
49.                     g=f.group(1)
50.                     if(g not in XYZ):
51.                         XYZ.append(g)
52.                 if(type(s[0])==sympy.geometry.point.Point3D):#if intersection is point type
53.                     e=d[8:(l-1)]
54.                     if(e not in XYZ):
55.                         XYZ.append(e)
56.     return XYZ
57.
58. #FUNCTION SOLVING FOR ISOPARAMENTRIC COORDINATES
59. #takes x,y,z coefficients, point coordinates, and symbols xi, eta, zeta as inputs
60. def isoparametric(a,b,c,r1,xi,eta,zeta):
61.     Q=[] #list storage for output
62.     #set up a system of linear equations for x,y,z in terms of xi,eta,zeta
63.     #the equation is solved and the solutions are returned
64.     iso_point=solve((
65.     a[0]+a[1]*xi+a[2]*eta+a[3]*zeta+a[4]*xi*eta+a[5]*eta*zeta+a[6]*xi*zeta+a[7]*xi*eta*zeta-r1[0],
66.     b[0]+b[1]*xi+b[2]*eta+b[3]*zeta+b[4]*xi*eta+b[5]*eta*zeta+b[6]*xi*zeta+b[7]*xi*eta*zeta-r1[1],
67.     c[0]+c[1]*xi+c[2]*eta+c[3]*zeta+c[4]*xi*eta+c[5]*eta*zeta+c[6]*xi*zeta+c[7]*xi*eta*zeta-r1[2]),
68.     xi, eta, zeta)
69.
70.     #perform tests for the output type
71.     #the data type will vary depending upon if there are multiple solutions or not
72.     test='good'
73.     try:
74.         Q.append(iso_point[0][0])
75.     except:
78.         Q.append(iso_point[xi])
79.         Q.append(iso_point[eta])
80.         Q.append(iso_point[zeta])
81.     if(test=='good'):
82.     #check whether each solution is real or not
86.             #fill Q will dummy values that will get rejected
87.             Q[0]=10000000
88.             Q.append(10000000)
89.             Q.append(10000000)
90.
91.
92.         if(iso_point[0][0].is_real ==True and
93.             iso_point[0][1].is_real ==True and
94.             iso_point[0][2].is_real ==True):
95.             #check if the point is within the isoparametric bounds
96.             if(round(min(iso_point[0]),7)<-1 or round(max(iso_point[0]),7)>1):
97.                 Q[0]=iso_point[1][0]
98.                 Q.append(iso_point[1][1])
99.                 Q.append(iso_point[1][2])
100.             if(round(min(iso_point[1]),7)<-1 or round(max(iso_point[1]),7)>1):
101.                 Q[0]=iso_point[0][0]
102.                 Q.append(iso_point[0][1])
103.                 Q.append(iso_point[0][2])
104.     return Q
105.
106. #JACOBIAN FUNCTION
107. #takes x,y,z coefficients, and lists of coordinates as input
108. def J1(a,b,c,x,y,z):
109.     #partial derivatives
110.     dx_1=a[1]+(a[4]*y)+(a[6]*z)+(a[7]*y*z)
111.     dx_2=a[2]+(a[4]*x)+(a[5]*z)+(a[7]*x*z)
112.     dx_3=a[3]+(a[5]*y)+(a[6]*x)+(a[7]*x*y)
113.
114.     dy_1=b[1]+(b[4]*y)+(b[6]*z)+(b[7]*y*z)
115.     dy_2=b[2]+(b[4]*x)+(b[5]*z)+(b[7]*x*z)
116.     dy_3=b[3]+(b[5]*y)+(b[6]*x)+(b[7]*x*y)
117.
118.     dz_1=c[1]+(c[4]*y)+(c[6]*z)+(c[7]*y*z)
119.     dz_2=c[2]+(c[4]*x)+(c[5]*z)+(c[7]*x*z)
120.     dz_3=c[3]+(c[5]*y)+(c[6]*x)+(c[7]*x*y)
121.
122.     #jacobian matrix
123.     B=np.zeros([3,3])
124.     B[0,:]=[dx_1,dx_2,dx_3]
125.     B[1,:]=[dy_1,dy_2,dy_3]
126.     B[2,:]=[dz_1,dz_2,dz_3]
127.     j=np.linalg.det(B)
128.
129.     return abs(j)
130.
131. #FUNCTION FOR FINDING INTEGRAL VALUES USING SHAPE FUNCTIONS
132. #takes x,y,z coefficients, point coordinates, and point values as input
133. def shape(a,b,c,x,y,z,u):
134. #matrix of coefficients for xi,et,zeta terms in shape functions
135.     A=np.ones([NEN,3])
136.     A[0,:]=[-1,-1,-1]
137.     A[1,:]=[1,-1,-1]
138.     A[2,:]=[1,1,-1]
139.     A[3,:]=[-1,1,-1]
140.     A[4,:]=[-1,-1,1]
141.     A[5,:]=[1,-1,1]
142.     A[6,:]=[1,1,1]
143.     A[7,:]=[-1,1,1]
144.
145.     #get jacobian
146.     j=J1(a,b,c,x,y,z)
147.
148.
149. #sum contributions from each shape function
150.     #called as part of gauss quadrature
151.     NSUM=0
152.     for i in range(NEN):
153.         NSUM=NSUM+((abs(j)*u[i]*(1+(x*A[i,0]))*(1+(y*A[i,1]))*(1+(z*A[i,2])))/8)
154.     return NSUM
155.
156. #FUNCTION FOR FINDING POINT VALUES USING SHAPE FUNCTIONS
157. #takes x,y,z coefficients, point coordinates, and point values as input
158. def shapeVal(a,b,c,x,y,z,u):
159.     A=np.ones([NEN,3])
160.     A[0,:]=[-1,-1,-1]
161.     A[1,:]=[1,-1,-1]
162.     A[2,:]=[1,1,-1]
163.     A[3,:]=[-1,1,-1]
164.     A[4,:]=[-1,-1,1]
165.     A[5,:]=[1,-1,1]
166.     A[6,:]=[1,1,1]
167.     A[7,:]=[-1,1,1]
168.
169.     #sum contributions from each shape function
170.     NSUM=0
171.     for i in range(NEN):
172.         NSUM=NSUM+(u[i]*(1+(x*A[i,0]))*(1+(y*A[i,1]))*(1+(z*A[i,2])))/8
173.     return NSUM
174.
175. #FUNCTION FOR TET SHAPE FUNCTION INTEGRATION
176. #takes array of x,y,z coordinates and values, point coordinates as input
177. def tet_iso(T,x,y,z):
178.     a=np.zeros([4,3]) #isoparametric coefficients for tets
179.     J=np.zeros([3,3]) #jacobian matrix
180.     p=[]
181.     for i in range(3):
182.         a[0,i]=T[0,i]
183.         a[1,i]=T[1,i]-T[0,i]
184.         a[2,i]=T[2,i]-T[0,i]
185.         a[3,i]=T[3,i]-T[0,i]
186.     for i in range(3):
187.         for j in range(3):
188.             J[i,j]=a[j+1,i]
189.     jacobian=np.linalg.det(J)
190.     U=abs(jacobian)*T[0,3]+(x*(T[1,3]-T[0,3]))+(y*(T[2,3]-T[0,3]))+(z*(T[3,3]-T[0,3]))
191.     return U
192.
193. #FUNCTION FOR HEX SHAPE FUNCTION INTEGRATION
194. #takes array of x,y,z coordinates and values,x,y,z coefficients, and  point coordinates as input
195. def hex_iso(T,a,b,c,x,y,z):
196.     A=np.ones([NEN,3])
197.     A[0,:]=[-1,-1,-1]
198.     A[1,:]=[1,-1,-1]
199.     A[2,:]=[1,1,-1]
200.     A[3,:]=[-1,1,-1]
201.     A[4,:]=[-1,-1,1]
202.     A[5,:]=[1,-1,1]
203.     A[6,:]=[1,1,1]
204.     A[7,:]=[-1,1,1]
205.
206.     #get jacobian
207.     j=J1(a,b,c,x,y,z)
208.
209.     NSUM=0
210.     for i in range(NEN):
211.         NSUM=NSUM+((abs(j)*T[i,3]*(1+(x*A[i,0]))*(1+(y*A[i,1]))*(1+(z*A[i,2])))/8)
212.     return NSUM
213.
215. #takes takes array of x,y,z coordinates and values, quadrature points and weights,
216. #and sets of x,y,z coefficients as input
218.     sumg=0
219.     L=len(T[:,0])
220.     for i in range(3):
221.         for j in range(3):
222.             for k in range(3):
223.                 if(L==4):
224.                     sumg=sumg+(tet_iso(T,xg[i],xg[j],xg[k])*wg[i]*wg[j]*wg[k]*J1(a,b,c,xg[i],xg[j],xg[k]))
225.                 if(L==8):
226.                     sumg=sumg+(hex_iso(T,a1,b1,c1,xg[i],xg[j],xg[k])*wg[i]*wg[j]*wg[k]*J1(a,b,c,xg[i],xg[j],xg[k]))
227.     return sumg
228.
229. #FUNCTION FOR FINDING HEX VOLUME
230. #takes x,y,z coordinate lists, and NEN as input
231. def hexVol(x,y,z,NEN):
232.     #fill h with node coordinates
233.     h=np.zeros([NEN,3])
234.     for i in range(NEN):
235.         h[i,0]=x[i]
236.         h[i,1]=y[i]
237.         h[i,2]=z[i]
238.     A=np.ones([4,4])
239.     B=np.ones([4,4])
240.     C=np.ones([4,4])
241.     D=np.ones([4,4])
242.     E=np.ones([4,4])
243.
244.     #fill matrices A,B,C,D,E using handles
245.     #matrices define the 5 tets needed to get hex volume
246.     for i in range(3):
247.         A[0,i]=h[0,i]
248.         A[1,i]=h[1,i]
249.         A[2,i]=h[3,i]
250.         A[3,i]=h[4,i]
251.
252.         B[0,i]=h[1,i]
253.         B[1,i]=h[2,i]
254.         B[2,i]=h[3,i]
255.         B[3,i]=h[6,i]
256.
257.         C[0,i]=h[1,i]
258.         C[1,i]=h[4,i]
259.         C[2,i]=h[5,i]
260.         C[3,i]=h[6,i]
261.
262.         D[0,i]=h[3,i]
263.         D[1,i]=h[4,i]
264.         D[2,i]=h[6,i]
265.         D[3,i]=h[7,i]
266.
267.         E[0,i]=h[1,i]
268.         E[1,i]=h[3,i]
269.         E[2,i]=h[4,i]
270.         E[3,i]=h[6,i]
271.
272.     VA=abs(np.linalg.det(A))
273.     VB=abs(np.linalg.det(B))
274.     VC=abs(np.linalg.det(C))
275.     VD=abs(np.linalg.det(D))
276.     VE=abs(np.linalg.det(E))
277.     V=(VA+VB+VC+VD+VE)/6
278.     return V
279.
280. #FUNCTION FOR CHECKING IF A POINT IS INSIDE A HEX VOLUME
281. #takes x,y,z coordinate lists,NEN, and point coordinates as input
282. def NewHexVol(x,y,z,NEN,p):
283.     #fill h with node coordinates
284.     h=np.zeros([NEN,3])
285.     for i in range(NEN):
286.         h[i,0]=x[i]
287.         h[i,1]=y[i]
288.         h[i,2]=z[i]
289.     A=np.ones([4,4])
290.     VA=0.0
291.     #find volumes of 12 tets cre4ated using the test point
292.     #A is iteratively set up as each of the 12 required tets
293.     for j in range(12):
294.         for i in range(3):
295.             if(j==0):
296.                 A[0,i]=h[0,i]
297.                 A[1,i]=h[1,i]
298.                 A[2,i]=h[2,i]
299.                 A[3,i]=p[i]
300.             if(j==1):
301.                 A[0,i]=h[0,i]
302.                 A[1,i]=h[3,i]
303.                 A[2,i]=h[2,i]
304.                 A[3,i]=p[i]
305.             if(j==2):
306.                 A[0,i]=h[0,i]
307.                 A[1,i]=h[3,i]
308.                 A[2,i]=h[4,i]
309.                 A[3,i]=p[i]
310.             if(j==3):
311.                 A[0,i]=h[0,i]
312.                 A[1,i]=h[1,i]
313.                 A[2,i]=h[4,i]
314.                 A[3,i]=p[i]
315.             if(j==4):
316.                 A[0,i]=h[1,i]
317.                 A[1,i]=h[2,i]
318.                 A[2,i]=h[5,i]
319.                 A[3,i]=p[i]
320.             if(j==5):
321.                 A[0,i]=h[1,i]
322.                 A[1,i]=h[4,i]
323.                 A[2,i]=h[5,i]
324.                 A[3,i]=p[i]
325.             if(j==6):
326.                 A[0,i]=h[2,i]
327.                 A[1,i]=h[5,i]
328.                 A[2,i]=h[6,i]
329.                 A[3,i]=p[i]
330.             if(j==7):
331.                 A[0,i]=h[2,i]
332.                 A[1,i]=h[3,i]
333.                 A[2,i]=h[7,i]
334.                 A[3,i]=p[i]
335.             if(j==8):
336.                 A[0,i]=h[2,i]
337.                 A[1,i]=h[7,i]
338.                 A[2,i]=h[6,i]
339.                 A[3,i]=p[i]
340.             if(j==9):
341.                 A[0,i]=h[5,i]
342.                 A[1,i]=h[6,i]
343.                 A[2,i]=h[4,i]
344.                 A[3,i]=p[i]
345.             if(j==10):
346.                 A[0,i]=h[4,i]
347.                 A[1,i]=h[6,i]
348.                 A[2,i]=h[7,i]
349.                 A[3,i]=p[i]
350.             if(j==11):
351.                 A[0,i]=h[4,i]
352.                 A[1,i]=h[3,i]
353.                 A[2,i]=h[7,i]
354.                 A[3,i]=p[i]
355.         #sum volumes
356.         VA=VA+abs(np.linalg.det(A))
357.     V=VA/6
358.     return V
359. #ELEMENT FUNCTION - perform operations on each PROTEUS element
360. def element(ee,nEl,MCNP_el_val,M,m,xyz,NEN,dim,vDat,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,center,max_dim):
361.     #storage lists
362.     points=[]
363.     x=[]
364.     y=[]
365.     z=[]
366.     val=[]
367.     #load the PROTEUS node coordinates and values
368.     #rounding is due to the python float data type appending digits
369.     #the source file provides 7 decimals
370.     for n in range(NEN):#loop over nodes
371.         x.append(round((xyz[n+(NEN*ee*dim)]),7))
372.         y.append(round((xyz[n+(NEN)+(NEN*ee*dim)]),7))
373.         z.append(round((xyz[n+(NEN*2)+(NEN*ee*dim)]),7))
374.         val.append(round((vDat[3,n+(NEN*ee)]),7))
375.
376.     #perform checks to see if PROTEUS elements are near the target MCNP element
377.     neighbor=False
378.
379.     #limiting search to a single z slice of 1 pin
380.     if(round(z[0],7)>=0.05 and round(z[0],7)<=0.15):
381.         #check the distance of each node from the MCNP element center
382.         for i in range(NEN):
383.             d=(((center[0]-x[i])**2)+((center[1]-y[i])**2)+((center[2]-z[i])**2))**0.5
384.             if(d<=((2**.5)*max_dim/2)):
385.                 neighbor=True
386.                 break
387.     #only run for neighboring nodes
388.     if(neighbor==True):
389.         A=np.ones([8,8])
390.         A[0,:]=[1,-1,-1,-1,1,1,1,-1]
391.         A[1,:]=[1,1,-1,-1,-1,1,-1,1]
392.         A[2,:]=[1,1,1,-1,1,-1,-1,-1]
393.         A[3,:]=[1,-1,1,-1,-1,-1,1,1]
394.         A[4,:]=[1,-1,-1,1,1,-1,-1,1]
395.         A[5,:]=[1,1,-1,1,-1,-1,1,-1]
396.         A[7,:]=[1,-1,1,1,-1,1,-1,-1]
397.         A_inv=np.linalg.matrix_power(A,-1)
398.
399.         #coefficients for isoparametric coordinate functions
400.         a=np.matmul(A_inv,x)
401.         b=np.matmul(A_inv,y)
402.         c=np.matmul(A_inv,z)
403.
404.         #array to store intersection points - assume maximum of 10 points
405.         N=np.zeros([10,3])
406.         pnts=0
407.         V=hexVol(x,y,z,NEN) #full PROTEUS element volume
408.         inside=0
409.         for i in range(NEN):
410.             #check if any PROTEUS nodes are inside the MCNP element
411.             if(x[i]>=XMIN and x[i]<=XMAX and y[i]>=YMIN and y[i]<=YMAX and z[i]>=ZMIN and z[i]<=ZMAX):
412.                 rr1=[x[i],y[i],z[i]]
413.                 Q=isoparametric(a,b,c,rr1,xi,eta,zeta)
414.                 qw.write(str(x[i])+','+str(y[i])+','+str(z[i])+'\n')
415.                 N[pnts,0]=round(Q[0],7)
416.                 N[pnts,1]=round(Q[1],7)
417.                 N[pnts,2]=round(Q[2],7)
418.                 pnts=pnts+1 #counter for number of intersection points
419.
420.                 inside=inside+1 #count the number of PROTEUS nodes inside the MCNP element
421.             r1=[xm[i],ym[i],zm[i]] #MCNP node coordinates
422.             V2=NewHexVol(x,y,z,NEN,r1)
423.             #check if any MCNP nodes are inside the PROTEUS element
424.             if(V2<=V):
425.                 Q=isoparametric(a,b,c,r1,xi,eta,zeta)
426.                 N[pnts,0]=round(Q[0],7)
427.                 N[pnts,1]=round(Q[1],7)
428.                 N[pnts,2]=round(Q[2],7)
429.                 pnts=pnts+1
430.         #create 3D point entities for finding intersections
431.         for i in range(len(x)):
432.             points.append(Point3D(x[i],y[i],z[i]))
433.
434.         if(inside!=NEN): #only execute the rest if all PROTEUS nodes are not inside MCNP element
435.             P=[] #PROTEUS faces
436.             p=[] #PROTEUS edges
437.
438.             #create lists of planes and edges for each element
439.             P.append(Plane(points[0],points[1],points[2]))
440.             P.append(Plane(points[0],points[1],points[4]))
441.             P.append(Plane(points[0],points[3],points[4]))
442.             P.append(Plane(points[5],points[6],points[2]))
443.             P.append(Plane(points[5],points[6],points[7]))
444.             P.append(Plane(points[2],points[3],points[6]))
445.
446.             p.append(Segment3D(points[0],points[1]))
447.             p.append(Segment3D(points[1],points[2]))
448.             p.append(Segment3D(points[2],points[3]))
449.             p.append(Segment3D(points[3],points[0]))
450.             p.append(Segment3D(points[4],points[5]))
451.             p.append(Segment3D(points[5],points[6]))
452.             p.append(Segment3D(points[6],points[7]))
453.             p.append(Segment3D(points[7],points[4]))
454.             p.append(Segment3D(points[1],points[5]))
455.             p.append(Segment3D(points[2],points[6]))
456.             p.append(Segment3D(points[3],points[7]))
457.             p.append(Segment3D(points[0],points[4]))
458.
459.
460.             XYZ=[] #storage for intersection coordinates
461.
462.             #check for PROTEUS edges intersecting MCNP planes
463.             #check for MCNP edges intersecting PROTEUS planes
464.             find_intersection(M,p,XYZ)
465.             find_intersection(P,m,XYZ)
466.             l=len(XYZ)
467.
468.             #iterate over all possible intersections
469.             for i in range(l):
470.                 r1=[]
471.                 Q=[]
472.                 #convert intersections to numeric format
473.                 r=XYZ[i].split(",")
474.                 for k in range(3):
475.                     r1.append(convert(r[k]))
476.                 Q=isoparametric(a,b,c,r1,xi,eta,zeta)
477.                 #check is Q is within PROTEUS element limits
478.                 if(round(max(Q),7)<=1 and round(min(Q),7)>=-1
479.                    and r1[0]>=XMIN and r1[0]<=XMAX and r1[1]>=YMIN
480.                    and r1[1]<=YMAX and r1[2]>=ZMIN and r1[2]<=ZMAX):
481.
482.                     qw.write(str(r1[0])+','+str(r1[1])+','+str(r1[2])+'\n')
483.                     #store valid points in N
484.                     N[pnts,0]=round((Q[0]),7)
485.                     N[pnts,1]=round((Q[1]),7)
486.                     N[pnts,2]=round((Q[2]),7)
487.                     pnts=pnts+1
488.             #at this point all the intersections are found
489.
490.             if(pnts>3):
491.                     #reorder the intersection points
492.                     #nodes are ordered in a counterclockwise fashion in each z plane
493.                     k=pnts-1
494.                     while (k>0):
495.                             l=0
496.                             for j in range(k):
497.                                 if((N[j,2]>N[j+1,2]) or (N[j,2] == N[j+1,2] and N[j,0] >  N[j+1,0]) or
498.                                 (N[j,2] == N[j+1,2] and N[j,0] == N[j+1,0] and N[j,1] > N[j+1,1])):
499.                                     for mm in range(3):
500.                                         tmp=N[j,mm]
501.                                         N[j,mm]=N[j+1,mm]
502.                                         N[j+1,mm]=tmp
503.                                     l=j
504.                                 k=l
505.                     int_vals=np.zeros([pnts])
506.                     #calculate value for each intersection point after reordering
507.                     for i in range(pnts):
508.                         int_vals[i]=shapeVal(a,b,c,N[i,0],N[i,1],N[i,2],val)
509.                     #assign nodes to each piece of the intersection
510.                     #assign nodal values
511.
512.                     #this section assigns the intersection nodes to matrices
513.                     #procedure varies by number of nodes
514.                     #some shapes are broken into simpler shapes
515.                     #integration is then done over the intersection regions
516.                         #the intersection is transformed to its own isoparametric coordinates
517.                         #jacobians from the original and new mapping are used in the calculations
518.
519.                     #hex - not split
520.                     if(pnts==8):
521.                     #lists storing x,y,z values
522.                         ax=[]
523.                         ay=[]
524.                         az=[]
525.                         Nh=np.zeros([8,4])
526.                         #fill Nh and do slight reordering
527.                         #order is modified for matrix solving
528.                         for i in range(3):
529.                             for j in range(8):
530.                                 Nh[j,i]=N[j,i]
531.                             Nh[1,i]=N[3,i]
532.                             Nh[3,i]=N[1,i]
533.                             Nh[5,i]=N[7,i]
534.                             Nh[7,i]=N[5,i]
535.                         for i in range(8):
536.                             Nh[i,3]=int_vals[i]
537.                             ax.append(Nh[i,0])
538.                             ay.append(Nh[i,1])
539.                             az.append(Nh[i,2])
540.                         Nh[1,3]=int_vals[3]
541.                         Nh[3,3]=int_vals[1]
542.                         Nh[5,3]=int_vals[7]
543.
544.                         #get isoparametric coefficients
545.                         a1=np.matmul(A_inv,ax)
546.                         b1=np.matmul(A_inv,ay)
547.                         c1=np.matmul(A_inv,az)
548.
550.
551.                         #tet - not split
552.                     if(pnts==4):
553.                         Nt1=np.zeros([4,4])
554.                         for i in range(3):
555.                             Nt1[0,i]=N[0,i]
556.                             Nt1[1,i]=N[1,i]
557.                             Nt1[2,i]=N[2,i]
558.                             Nt1[3,i]=N[3,i]
559.                         Nt1[0,3]=int_vals[0]
560.                         Nt1[1,3]=int_vals[1]
561.                         Nt1[2,3]=int_vals[2]
562.                         Nt1[3,3]=int_vals[3]
563.                         print(Nt1)
565.
566.                     #pentahedral - split into 1 hex and 3 tets
567.                     if(pnts==10):
568.                         ax=[]
569.                         ay=[]
570.                         az=[]
571.                         Nh=np.zeros([8,4])
572.                         Nt1=np.zeros([4,4])
573.                         Nt2=np.zeros([4,4])
574.                         Nt3=np.zeros([4,4])
575.                         for i in range(3):
576.                             Nh[0,i]=N[0,i]
577.                             Nh[1,i]=N[4,i]
578.                             Nh[2,i]=N[2,i]
579.                             Nh[3,i]=N[1,i]
580.                             Nh[4,i]=N[5,i]
581.                             Nh[5,i]=N[9,i]
582.                             Nh[6,i]=N[7,i]
583.                             Nh[7,i]=N[6,i]
584.
585.                             Nt1[0,i]=N[7,i]
586.                             Nt1[1,i]=N[8,i]
587.                             Nt1[2,i]=N[9,i]
588.                             Nt1[3,i]=N[2,i]
589.
590.                             Nt2[0,i]=N[8,i]
591.                             Nt2[1,i]=N[9,i]
592.                             Nt2[2,i]=N[2,i]
593.                             Nt2[3,i]=N[3,i]
594.
595.                             Nt3[0,i]=N[9,i]
596.                             Nt3[1,i]=N[2,i]
597.                             Nt3[2,i]=N[3,i]
598.                             Nt3[3,i]=N[4,i]
599.
600.                         for i in range(8):
601.                             ax.append(Nh[i,0])
602.                             ay.append(Nh[i,1])
603.                             az.append(Nh[i,2])
604.                         Nh[0,3]=int_vals[0]
605.                         Nh[1,3]=int_vals[4]
606.                         Nh[2,3]=int_vals[2]
607.                         Nh[3,3]=int_vals[1]
608.                         Nh[4,3]=int_vals[5]
609.                         Nh[5,3]=int_vals[9]
610.                         Nh[6,3]=int_vals[7]
611.                         Nh[7,3]=int_vals[6]
612.
613.                         a1=np.matmul(A_inv,ax)
614.                         b1=np.matmul(A_inv,ay)
615.                         c1=np.matmul(A_inv,az)
616.
617.                         Nt1[0,3]=int_vals[7]
618.                         Nt1[1,3]=int_vals[8]
619.                         Nt1[2,3]=int_vals[9]
620.                         Nt1[3,3]=int_vals[2]
621.
622.                         Nt2[0,3]=int_vals[8]
623.                         Nt2[1,3]=int_vals[9]
624.                         Nt2[2,3]=int_vals[2]
625.                         Nt2[3,3]=int_vals[3]
626.
627.                         Nt3[0,3]=int_vals[9]
628.                         Nt3[1,3]=int_vals[2]
629.                         Nt3[2,3]=int_vals[3]
630.                         Nt3[3,3]=int_vals[4]
631.
635.
636.                     #prism - split into tets
637.                     if(pnts==6):
638.                         Nt1=np.zeros([4,4])
639.                         Nt2=np.zeros([4,4])
640.                         Nt3=np.zeros([4,4])
641.                         for i in range(3):
642.
643.                             Nt1[0,i]=N[3,i]
644.                             Nt1[1,i]=N[4,i]
645.                             Nt1[2,i]=N[5,i]
646.                             Nt1[3,i]=N[1,i]
647.
648.                             Nt2[0,i]=N[5,i]
649.                             Nt2[1,i]=N[3,i]
650.                             Nt2[2,i]=N[1,i]
651.                             Nt2[3,i]=N[2,i]
652.
653.                             Nt3[0,i]=N[3,i]
654.                             Nt3[1,i]=N[1,i]
655.                             Nt3[2,i]=N[2,i]
656.                             Nt3[3,i]=N[0,i]
657.
658.                         Nt1[0,3]=int_vals[3]
659.                         Nt1[1,3]=int_vals[4]
660.                         Nt1[2,3]=int_vals[5]
661.                         Nt1[3,3]=int_vals[1]
662.
663.                         Nt2[0,3]=int_vals[5]
664.                         Nt2[1,3]=int_vals[3]
665.                         Nt2[2,3]=int_vals[1]
666.                         Nt2[3,3]=int_vals[2]
667.
668.                         Nt3[0,3]=int_vals[3]
669.                         Nt3[1,3]=int_vals[1]
670.                         Nt3[2,3]=int_vals[2]
671.                         Nt3[3,3]=int_vals[0]
672.
675.         #this part of the code is called is all PROTEUS nodes are inside the MCNP element
676.         if(inside==NEN):
677.             sumg=0.0 #initialize summation
679.             for i in range(3):
680.                 for j in range(3):
681.                     for k in range(3):
682.                         sumg=sumg+(shape(a,b,c,xg[i],xg[j],xg[k],val)*wg[i]*wg[j]*wg[k])
683.             MCNP_el_val=MCNP_el_val+sumg # add to MCNP element total
684.
685.     return MCNP_el_val
686.
687. #MAIN FUNCTION
688.
689. #symbols for equation solving
690. xi=Symbol('xi')
691. eta=Symbol('eta')
692. zeta=Symbol('zeta')
693.
694. #get points and weights for gauss quadrature
695. xg, wg = gauss_legendre(3, 10)
696.
698. f=h5py.File(r'C:\Users\Peter\Documents\FEP\Project\3D Mesh\RCF_2x2.h5')
699. #get the hdf5 keys
700. tier1=list(f.keys())
701. con=f[('CONTROL')]
702. blocks=con[0] #number of material regions
703. dim=con[1] #dimension
704. #open output files
705. qw=open('intersections.csv','w')
706. windows=open('windows.csv','w')
707.
708.
709. #MCNP mesh parameters
710. #outer limits and number of divisions
711. iints=2
712. jints=2
713. kints=1
714. Xhigh=0.9
715. Yhigh=0.9
716. Zhigh=.1
717. Xlow=-0.9
718. Ylow=-0.9
719. Zlow=0
720. del_x=(Xhigh-Xlow)/iints
721. del_y=(Yhigh-Ylow)/jints
722. del_z=(Zhigh-Zlow)/kints
723.
724. #loop iterating over the MCNPm elements
725. for mk in range(kints):
726.     for mj in range(jints):
727.         for mi in range(iints):
728.
729.             MCNP_el_val=0.0 #initialize the value for the element
730.
731.             #MCNP element limits
732.             XMIN=Xlow+(mi*del_x)
733.             YMIN=Ylow+(mj*del_y)
734.             ZMIN=Zlow+(mk*del_y)
735.             XMAX=XMIN+del_x
736.             YMAX=YMIN+del_y
737.             ZMAX=ZMIN+del_z
738.             #find the element center
739.             center=[]
740.             center.append(XMIN+((XMAX-XMIN)/2))
741.             center.append(YMIN+((YMAX-YMIN)/2))
742.             center.append(ZMIN+((ZMAX-ZMIN)/2))
743.             #find the maximim distance to a face of the element
744.             Dist=[]
745.             Dist.append(XMAX-XMIN)
746.             Dist.append(YMAX-YMIN)
747.             Dist.append(ZMAX-ZMIN)
748.             max_dim=max(Dist)
749.
750.             #MCNP element node coordinates
751.             xm=[XMIN,XMAX,XMAX,XMIN,XMIN,XMAX,XMAX,XMIN]
752.             ym=[YMIN,YMIN,YMAX,YMAX,YMIN,YMIN,YMAX,YMAX]
753.             zm=[ZMIN,ZMIN,ZMIN,ZMIN,ZMAX,ZMAX,ZMAX,ZMAX]
754.
755.
756.             #setup SymPy style points, edges, and planes for finding intersections
757.             #create Point3D entities using MCNP node coordinates
758.             #the Point3D entities are used for constructing planes and edges
759.             pointsM=[]
760.             for i in range(8):
761.                     pointsM.append(Point3D(xm[i],ym[i],zm[i]))
762.
763.             M=[] #MCNP faces
764.             m=[] #MCNP edges
765.
766.             #fill lists with planes and edges
767.             M.append(Plane(pointsM[0],pointsM[1],pointsM[2]))
768.             M.append(Plane(pointsM[0],pointsM[1],pointsM[4]))
769.             M.append(Plane(pointsM[0],pointsM[3],pointsM[4]))
770.             M.append(Plane(pointsM[5],pointsM[6],pointsM[2]))
771.             M.append(Plane(pointsM[5],pointsM[6],pointsM[7]))
772.             M.append(Plane(pointsM[2],pointsM[3],pointsM[6]))
773.
774.             m.append(Segment3D(pointsM[0],pointsM[1]))
775.             m.append(Segment3D(pointsM[1],pointsM[2]))
776.             m.append(Segment3D(pointsM[2],pointsM[3]))
777.             m.append(Segment3D(pointsM[3],pointsM[0]))
778.             m.append(Segment3D(pointsM[4],pointsM[5]))
779.             m.append(Segment3D(pointsM[5],pointsM[6]))
780.             m.append(Segment3D(pointsM[6],pointsM[7]))
781.             m.append(Segment3D(pointsM[7],pointsM[4]))
782.             m.append(Segment3D(pointsM[1],pointsM[5]))
783.             m.append(Segment3D(pointsM[2],pointsM[6]))
784.             m.append(Segment3D(pointsM[3],pointsM[7]))
785.             m.append(Segment3D(pointsM[0],pointsM[4]))
786.
787.             total=0 #initialize contributions from PROTEUS elements
788.
789.             #loop over blocks
790.             for bl in range(blocks):
791.                 dset=f[tier1[bl]]
792.                 subD=list(dset)
793.                 info=dset[(subD[2])]
794.                 vDat=dset[(subD[3])] #output data
795.                 xyz=dset[(subD[4])] #node coordinates
796.                 NEN=info[1] #nodes per element
797.                 nEl=info[0] #elements per block
798.
799.                 #loop over elements
800.                 for ee in range(nEl):
801.                     #call element function to get contribution from each element
802.                     output=element(ee,nEl,MCNP_el_val,M,m,xyz,NEN,dim,vDat,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,center,max_dim)
803.                     total=total+output #add contributions to the total
804.
805.             print('total: '+str(total))
806.             #write output to file once all element contributions from all blocks are summed
807.             windows.write(str(total)+',\n')
808.             #move on to next MCNP element
809.
810. qw.close()
811. windows.close()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top