## module m_prepBIE ''' Module containing procedures, which compile the *.msh mesh file of the studied domain. The *.msh file should be generated by the GMSH 2.8.5 as 2D mesh. The *.msh file can generate mesh such the way that the studied domain is separated to the subdomains in the form of 'chessboard'. The two subdomains is a minimum. The boundary of each subdomain has to be oriented counterclockwise, holes clocwise. Moreover, the material properties of 'black' and 'white' subdomains can mismatch. For better understanding look at *.geo examples. Classes: matdomchar,mesh ''' from math import sqrt class matdomchar(object): ''' Generates the lists of elastic constants and lists of nodes, boundary elements and triangular elements of both subdomains. Inputs are plane stress/strain conditions, elastic constants and name of the msh file. There are also the methods searching node coordinates, nodes of the boundary and domain elements. Input(example): (iplane,el,pn,file_name) #plane stress/strain - 1/0: iplane=1 #Young's modulus of material 1 and 2: el=(2.1E+05,2.1E+05) #Poissons' ratio for material 1 and 2: pn=(0.3E+00,0.3E+00) #MSH file name file_name='file.msh' Parameters: iplane - [int] plane stress/strain el1 - [float] Young's modulus of 1st subdomain el2 - [float] Young's modulus of 2nd subdomain pn1 - [float] Poissons'ratio of 1st subdomain pn2 - [float] Poissons'ratio of 2nd subdomain file_name - [string] msh file name Methods: load() - generates dictionaries Nodes - No - [int] node ID number - x - [float] node x-coordinate - y - [float] node y-coordinate - z - [float] node z-coordinate (should be zero) Elements - No - [int] boundary element ID number - line - [int] line ID number - node1 - [int] node 1 of the boundary element - node2 - [int] node 2 of the boundary element TriElements - No - [int] subdomain element ID number - dom - [int] subdomain ID number - node1 - [int] ID of the subdomain element node 1 - node2 - [int] ID of the subdomain element node 2 - node3 - [int] ID of the subdomain element node 3 nNode(n) - returns tuple (x,y) of x [float] and y [float] coordinate of the node with id ID number n [int] lElements(n) - returns a list of tuples (x1,y1,x2,y2) of x [float] and y [float] coordinates of node1 and node2 of the line boundary elements of the boundary line with ID number n [int] NodeElementPlot() - returns the tuple (xynodes,dom,elements) of the subdomain mesh informations dedicated to the matplotlib plotting, where - xynodes is list of the tuples (x,y), x [float] and y [float] coordinates of the triangle nodes of the domain, - dom is list of ID number [int] of subdomain in which the coresponding node lies - elements is list of tuples (index1 [int],index2 [int],index3 [int]) of node indexes corresponding to the list xynodes ''' def __init__(self,plane,E,nu,name): #default values of input parameters self.iplane=plane self.el1=E[0] self.el2=E[1] self.pn1=nu[0] self.pn2=nu[1] self.file_name=name self.load() def load(self): self.Nodes=[] self.Elements=[] self.TriElements=[] with open(self.file_name,'r') as f: while True: try: line=f.next() if line[:6]=='$Nodes': line=f.next() for ii in range(int(line)): No,x,y,z=f.next().split() self.Nodes.append({'No':int(No),'x':float(x),'y':float(y),'z':float(z)}) if line[:9]=='$Elements': line=f.next() for ii in range(int(line)): a=[int(jj) for jj in f.next().split()] if a[1]==1: self.Elements.append({'No':a[0],'line':a[4],'node1':a[5],'node2':a[6]}) if a[1]==2: self.TriElements.append({'No':a[0],'dom':a[4],'node1':a[5],'node2':a[6],'node3':a[7]}) except StopIteration: break def nNode(self,n): ii=0 while self.Nodes[ii]['No']!=n: ii+=1 return (self.Nodes[ii]['x'],self.Nodes[ii]['y']) def lElements(self,nline): elements=[] for ii in self.Elements: if ii['line']==nline: xy1=self.nNode(ii['node1']) xy2=self.nNode(ii['node2']) elements.append((xy1[0],xy1[1],xy2[0],xy2[1])) return elements def NodeElementPlot(self): dom=[0 for ii in self.Nodes] idnodes=[ii['No'] for ii in self.Nodes] xynodes=[(ii['x'],ii['y']) for ii in self.Nodes] elements=[(idnodes.index(ii['node1']),idnodes.index(ii['node2']),idnodes.index(ii['node3'])) for ii in self.TriElements] for ii in self.TriElements: dom[idnodes.index(ii['node1'])]=ii['dom'] dom[idnodes.index(ii['node2'])]=ii['dom'] dom[idnodes.index(ii['node3'])]=ii['dom'] return xynodes,dom,elements class mesh(matdomchar): ''' Generates the lists of dictionaries elements1 and elements2 in the form, which the class >solution< of the module m_procBIE.py requirein input Input: see class matdomchar doms - tuple of two lists (dom1,dom2) prescribing boundary conditions along the boundaries of both subdomains in the form - dom1 (example) is composed of four lines with ID 1,2,4,6 and prescribed displacements, stress free boundaries and mutual subdomain2 boundary, dom1=[(1,'uu',(0.,0.)),(2,'tt',(0.,0.)),(4,'tt',(0.,0.)),(6,'12')] - dom2 (example) is composed of mutual subdomain1 line and line with ID 3 and 5 with prescribed boundary conditions combining displacements and tractions, the mutual boundary is omitted, dom2=[(3,'ut',(0.,100.)),(5,'tu',(0.,0.))] Parameters: dom1 - [list] list of boundary conditions of subdomain1 dom2 - [list] list of boundary conditions of subdomain2 elements1/elements2 - [list]/[list] of dictionaries in which 0: [float tuple] coordinates of the midpoint of the element 1: [float tuple] coordinates of the first endpoint of the element 2: [float tuple] coordinates of he second endpoint of the element '0.5*l': [float] halflength of the element '$b': [string] kind of the boundary condition - tt,uu,tu,ut,12 'bxy': [float tuple] coordinates of the boundary condition 'n': [float tuple] coordinates of the external normal to the element Methods: make_mesh() - generates the lists elements1 and elements2 ''' # def __init__(self,plane,E,nu,name,doms): super(mesh,self).__init__(plane,E,nu,name) self.dom1=doms[0] self.dom2=doms[1] self.make_mesh() def make_mesh(self): self.elements1,self.elements2=[],[] elements12,elements34=[],[] for ii in self.dom1: if ii[1]=='12': elements=self.lElements(ii[0]) for jj in elements: elements12.append({1:(jj[0],jj[1]), \ 0:(0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])), \ 2:(jj[2],jj[3]), \ '0.5*l':0.5*sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ 'n':((jj[3]-jj[1])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ -(jj[2]-jj[0])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2)), \ 'b$':'12'}) elements34.append({1:(jj[2],jj[3]), \ 0:(0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])), \ 2:(jj[0],jj[1]), \ '0.5*l':0.5*sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ 'n':(-(jj[3]-jj[1])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ +(jj[2]-jj[0])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2)), \ 'b$':'12'}) else: elements=self.lElements(ii[0]) for jj in elements: self.elements1.append({1:(jj[0],jj[1]), \ 0:(0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])), \ 2:(jj[2],jj[3]), \ '0.5*l':0.5*sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ 'n':((jj[3]-jj[1])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ -(jj[2]-jj[0])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2)), \ 'b$':ii[1], \ 'bxy':(ii[2][0](0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])), \ ii[2][1](0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])))}) self.elements1.extend(elements12) for ii in self.dom2: elements=self.lElements(ii[0]) for jj in elements: self.elements2.append({1:(jj[0],jj[1]), \ 0:(0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])), \ 2:(jj[2],jj[3]), \ '0.5*l':0.5*sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ 'n':((jj[3]-jj[1])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2), \ -(jj[2]-jj[0])/sqrt((jj[2]-jj[0])**2+(jj[3]-jj[1])**2)), \ 'b$':ii[1], \ 'bxy':(ii[2][0](0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])), \ ii[2][1](0.5*(jj[0]+jj[2]),0.5*(jj[1]+jj[3])))}) self.elements2.extend(elements34[::-1])