Pas de problème, mais c'est un peu long et complexe
Voici les classes définissant les objets géométriques
class node:
"""
node class
"""
def __init__(self,num,x,y,z):
self._globNum = num
self._x = x
self._y = y
self._z = z
def getNum(self):
return self._globNum
def getXYZ(self):
return (self._x,self._y,self._z)
def numberOfNodes(self):
return 1
def getInterpolation(self,U,u,v,w):
return u(1)
########################################################################
class line(node):
"""
line class
"""
def __init__(self,num,tag,partition,node1,node2):
self._globNum = num
self._tag = tag
self._partition = partition
self._node1 = node1
self._node2 = node2
def getNum(self):
return self._globNum
def getXYZ(self):
return (self._x,self._y,self._z)
def getTag(self):
return self._tag
def getPartition(self):
return self._partition
def getNodes(self):
return [self._node1,self._node2]
def numberOfNodes(self):
return 2
def getShapesFunctions(self,u,v,w):
shf = [ 0.5 * (1-u), 0.5 * (1+u) ]
return shf
def getGradShapesFunctions(self,u,v,w):
gx = [-0.5 , 0.5]
gy = [ 0.0 , 0.0]
gz = [ 0.0 , 0.0]
gdsf = [ gx , gy , gz]
return gdsf
def getInterpolation(self,U,u,v,w):
shf = self.getShapesFunctions(u,v,w)
value = shf(0) * U(0) + shf(1) * U(1)
return value
########################################################################
class triangle(node):
"""
triangle class
"""
def __init__(self,num,tag,partition,node1,node2,node3):
self._globNum = num
self._tag = tag
self._partition = partition
self._node1 = node1
self._node2 = node2
self._node3 = node3
def getNum(self):
return self._globNum
def getXYZ(self):
return (self._x,self._y,self._z)
def getTag(self):
return self._tag
def getPartition(self):
return self._partition
def getNodes(self):
return [self._node1,self._node2,self._node3]
def numberOfNodes(self):
return 3
def getShapesFunctions(self,u,v,w):
shf = [ (1 - u - v) , u , v ]
return shf
def getGradShapesFunctions(self,u,v,w):
gx = [ -1.0 , 1.0 , 0.0]
gy = [ -1.0 , 0.0 , 1.0]
gz = [ 0.0 , 0.0 , 0.0]
gdsf = [ gx , gy , gz]
return gdsf
def getInterpolation(self,U,u,v,w):
shf = self.getShapesFunctions(u,v,w)
value = shf(0) * U(0) + shf(1) * U(1) + shf(2) * U(2)
return value
#########################################################################
class tetrahedron(node):
"""
tetrahedron class
"""
def __init__(self,num,tag,partition,node1,node2,node3,node4):
self._globNum = num
self._tag = tag
self._partition = partition
self._node1 = node1
self._node2 = node2
self._node3 = node3
self._node4 = node4
def getNum(self):
return self._globNum
def getXYZ(self):
return (self._x,self._y,self._z)
def getTag(self):
return self._tag
def getPartition(self):
return self._partition
def getNodes(self):
return [self._node1,self._node2,self._node3,self._node4]
def numberOfNodes(self):
return 4
def getShapesFunctions(self,u,v,w):
shf = [ (1 - u - v - w) , u , v , w]
return shf
def getGradShapesFunctions(self,u,v,w):
gx = [ -1.0 , 1.0 , 0.0 , 0.0 ]
gy = [ -1.0 , 0.0 , 1.0 , 0.0 ]
gz = [ -1.0 , 0.0 , 0.0 , 1.0 ]
gdsf = [ gx , gy , gz]
return gdsf
def getInterpolation(self,U,u,v,w):
shf = self.getShapesFunctions(u,v,w)
value = shf(0) * U(0) + shf(1) * U(1) + shf(2) * U(2) + shf(3) * U(3)
return value
########################################################################################
La classe définissant le maillage
from geom import *
class mesh :
"""
mesh class
"""
def __init__(self,file_name):
"""
constructor of the mesh class, from a gmsh file
"""
self._NODES = []
self._LINES = []
self._TRIANGLES = []
self._TETRAHEDRONS = []
#open the file
f = open(file_name,"r")
# nodes lecture
str = "a"
while str != "$Nodes":
str = f.readline().strip()
self._numNodes = int(f.readline())
for i in xrange(self._numNodes):
str = f.readline().strip().split()
n = node(int(str[0]) , float(str[1]) , float(str[2]) , float(str[3]))
self._NODES.append(n)
str = f.readline()
# elements lecture
str = f.readline()
self._numElements = int(f.readline())
for i in xrange(self._numElements):
str = f.readline().strip().split()
if str[1] == "1":
num = int(str[0])
tag = int(str[3])
partition = int(str[5])
numNode1 = int(str[6]) - 1
numNode2 = int(str[7]) - 1
node1 = self._NODES[numNode1]
node2 = self._NODES[numNode2]
n = line(num,tag,partition,node1,node2)
self._LINES.append(n)
elif str[1] == "2":
num = int(str[0])
tag = int(str[3])
partition = int(str[5])
numNode1 = int(str[6]) - 1
numNode2 = int(str[7]) - 1
numNode3 = int(str[8]) - 1
node1 = self._NODES[numNode1]
node2 = self._NODES[numNode2]
node3 = self._NODES[numNode3]
n = triangle(num,tag,partition,node1,node2,node3)
self._TRIANGLES.append(n)
elif str[1] == "4":
num = int(str[0])
tag = int(str[3])
partition = int(str[5])
numNode1 = int(str[6]) - 1
numNode2 = int(str[7]) - 1
numNode3 = int(str[8]) - 1
numNode4 = int(str[9]) - 1
node1 = self._NODES[numNode1]
node2 = self._NODES[numNode2]
node3 = self._NODES[numNode3]
node4 = self._NODES[numNode4]
n = tetrahedron(num,tag,partition,node1,node2,node3,node4)
self._TETRAHEDRONS.append(n)
else:
if str[1] != "15":
print "There are non implemented type of element in the mesh"
print "-Mesh reading"
print "Mesh file name :", file_name
print "Number of nodes :", self._numNodes
print "Number of lines :", len(self._LINES)
print "Number of triangles :", len(self._TRIANGLES)
print "Number of tetrahedrons :", len(self._TETRAHEDRONS)
print "-End of mesh reading "
# close the file
f.close()
def getNumNodes(self):
return self._numNodes
def getElementsByTag(self,dim,tag):
"""
function returning a list of elements of dimension 'dim' and tag 'tag'
"""
elements = []
if dim == 1:
for line in self._LINES:
if line.getTag() == tag:
elements.append(line)
elif dim == 2:
for triangle in self._TRIANGLES:
if triangle.getTag() == tag:
elements.append(triangle)
elif dim == 3:
for tetrahedron in self._TETRAHEDRONS:
if tetrahedron.getTag() == tag:
elements.append(tetrahedron)
return elements
def getVerticesByTag(self,dim,tag):
"""
function returning a list a elements if dimension 'dim' and tag 'tag'
"""
elements = self.getElementsByTag(dim,tag)
vertices = set([])
for elem in elements :
a = elem.getNodes()
vertices.update(a)
return vertices
##############################################################################################
La classe gérant les degrés de libertés
from geom import *
from collections import *
from linear_system import *
class dofkey:
def __init__(self,node,name):
self._node = node
self._name = name
def __eq__(self,other):
node = other.getNode()
name = other.getName()
if (self._node == node) and (self._name == name):
return True
else:
return False
def __hash__(self):
value = hash(self._name)
return value
def getNode(self):
return self._node
def getName(self):
return self._name
class dof_manager:
def __init__(self):
self._numbering = defaultdict()
self._fixed = defaultdict()
def fixVertex(self,node,name,val):
dof = dofkey(node,name)
self._fixed[dof]=val
def numberVertex(self,node,name):
dof = dofkey(node,name)
if self._fixed.get(dof,"_") == "_":
self._numbering[dof]=(len(self._numbering))
def fixVertices(self,nodes,name,val):
fix = self.fixVertex
for node in nodes:
fix(node,name,val)
def numberVertices(self,nodes,name):
number = self.numberVertex
for node in nodes:
number(node,name)
def getNumberingSize(self):
return len(self._numbering)
def getFixedSize(self):
return len(self._fixed)
def linkToLinearSystem(self,lsys):
self._lsys = lsys
def fixInSolution(self,node,name,val):
dof = dofkey(node,name)
num = self._numbering.get(dof,-1)
if not num == -1:
self._lsys.fixInSolution(num,val)
def addToMatrix(self,dofi,dofj,val):
if (dofi in self._numbering.keys()):
if (dofj in self._numbering.keys()):
R = self._numbering[dofi]
C = self._numbering[dofj]
self._lsys.addToMatrix(R,C,val)
else:
R = self._numbering[dofi]
F = self._fixed[dofj]
self._lsys.addToRHS(R,-val*F)
def addToRHS(self,dofi,val):
if (dofi in self._numbering.keys()):
R = self._numbering[dofi]
self._lsys.addToRHS(R,val)
def assembleInMatrix(self,node1,node2,n,local):
unknows = ["RHO","P","UX","UY","UZ","Y","e","eps","T","k","h","YEQ","DEFF"]
for i in xrange(n-1):
f = "X" + str(i+1)
unknows.append(f)
unknows.append("XN")
for i in xrange(n-1):
f = "T" + str(i+1)
unknows.append(f)
unknows.append("TN")
for i in xrange(13+2*n):
dofi = dofkey(node1,unknows[i])
for j in xrange(13+2*n):
dofj = dofkey(node2,unknows[j])
self.addToMatrix(dofi,dofj,local[i,j])
def assembleInVector(self,node1,n,local):
unknows = ["RHO","P","UX","UY","UZ","Y","e","eps","T","k","h","DEFF"]
for i in xrange(n-1):
f = "X" + str(i+1)
unknows.append(f)
unknows.append("XN")
for i in xrange(n-1):
f = "T" + str(i+1)
unknows.append(f)
unknows.append("TN")
for i in xrange(13+2*n):
dofi = dofkey(node1,unknows[i])
self.addToRHS(dofi,local[i])
def getNewDofValue(self,node,name):
dof = dofkey(node,name)
if (dof in self._numbering.keys()):
R = self._numbering[dof]
val = self._lsys.getFromNewSolution(R)
else:
val = self._fixed[dof]
return val
def getOldDofValue(self,node,name):
dof = dofkey(node,name)
if (dof in self._numbering.keys()):
R = self._numbering[dof]
val = self._lsys.getFromOldSolution(R)
else:
val = self._fixed[dof]
return val
def systemSolve(self):
self._lsys.solve()
def getSolution(self):
return self._lsys.getSolution()
