API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 25th, 2019
113
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 29.83 KB | None | 0 0
raw download clone embed print report
  1. #load python modules used
  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
  9. from sympy.integrals.quadrature import gauss_legendre
  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:
  76.         test='bad'
  77.     if(test=='bad'):
  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
  83.         if(iso_point[0][0].is_Add ==True or
  84.             iso_point[0][1].is_Add ==True or
  85.             iso_point[0][2].is_Add ==True):
  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.  
  214. #GAUSS QUADRATURE FUNCTION
  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    
  217. def Gquad(T,xg,wg,a1,b1,c1,a,b,c):
  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.                        
  549.                         MCNP_el_val=MCNP_el_val+(Gquad(Nh,xg,wg,a1,b1,c1,a,b,c))
  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)
  564.                         MCNP_el_val=MCNP_el_val+(Gquad(Nt1,xg,wg,0,0,0,a,b,c))
  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.                        
  632.                         MCNP_el_val=MCNP_el_val+(Gquad(Nt1,xg,wg,0,0,0,a,b,c)+
  633.                         Gquad(Nt2,xg,wg,0,0,0)+Gquad(Nt3,xg,wg,0,0,0,a,b,c)+
  634.                         Gquad(Nh,xg,wg,a1,b1,c1,a,b,c))
  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.                        
  673.                         MCNP_el_val=MCNP_el_val+(Gquad(Nt1,xg,wg,0,0,0,a,b,c)+
  674.                         Gquad(Nt2,xg,wg,0,0,0,a,b,c)+Gquad(Nt3,xg,wg,0,0,0,a,b,c))
  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
  678.             #perform gauss quadrature
  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.  
  697. #load PROTEUS output file
  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()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!