# Demo entry 6357579

py

Submitted by anonymous on Apr 23, 2017 at 11:01
Language: Python 3. Code size: 13.1 kB.

```#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
define the points and meshes
Created on Sun Apr 16 14:24:08 2017

@author: jf
"""
import numpy as np
import math
import random

class Element:
'''
this is a class that contain
(point)the point coordinates data and
(element)serial number data of element and
(elementType)material data of element
'''
def __init__(self):
print ('call Element')
# memery the coordinates of all points into a dic{pointBLB:[x y z]}
# in point {} , pointBLB is the serial number of point
self.point = {}
self.element = {}#memery the points of each elements  into a tuple,{pointBLB:(pointBLB,........)}
self.elementType = {}#memery the element what material it is{pointBLB:1, , ,}
self.point_element()
def point_element(self):
length=800 # the length of the RVE cube
coorStart = -length/2 # the min coordinate
edgePointNum = 101 #points number of each edge of RVE cube
elementLength = length/(edgePointNum-1) # the length of the edge of element
'''
the order of the points in the an element
8 -------- 7
/|         /|
5----------6 |
| |        | |
| |        | |
| 4----------3
|/         |/
1----------2

and the serial number of the points in the whole RVE cube
n*n*n-n+1----------------- n*n*n
/|                         /|
n*n(n-1)+1--------n*n(n-1)+n |
| |                        | |
| |                        | |
| |                        | |
| |                        | |
| n*(n-1)+1------------------n*n
/                         | /
1---------------------------n

the serial number in an element,
p is the serial number of the bottom-left-behind point
p+n+n*n--- p+n*n+n+1
/|         /|
p+n*n------p+n*n+1
| |        | |
| |        | |
| p+n-------p+1+n
|/         |/
p----------P+1

the coordinate system
z
| y
|/
O---x
'''

'''calculate the coordinates of every point'''
'''the serials number of blow-left-behind point in an element'''
pointBLB=1
for k in range(1,edgePointNum+1):# the z dimension
for j in range(1,edgePointNum+1):# the y dimension
for i in range(1,edgePointNum+1):# the x dimension
self.point[pointBLB]=[((i-1)*elementLength + coorStart),#x
((j-1)*elementLength + coorStart),#y
((k-1)*elementLength + coorStart)]#z
'''caculate the point number of elements '''
if not (i == edgePointNum or j == edgePointNum or\
k == edgePointNum):
self.element[pointBLB]=(pointBLB,
pointBLB+1,
pointBLB+1+edgePointNum,
pointBLB+edgePointNum,
pointBLB+edgePointNum*edgePointNum,
pointBLB+edgePointNum*edgePointNum+1,
pointBLB+edgePointNum*edgePointNum+1+edgePointNum,
pointBLB+edgePointNum*edgePointNum+edgePointNum)
self.elementType[pointBLB]=1
pointBLB= pointBLB+1

def trans_matrix(a,b,c,alfa,beta,gama):
'''return a translation matrix
a,b,c are the translation at x,y,z axis, respectively.
alfa,beta,gama are the rotation
around x,y,z axis, respectively
'''
translation = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[a,b,c,1]])#translation
rotation_x = np.array([[1,0,0,0],[0,math.cos(alfa),math.sin(alfa),0],
[0,-math.sin(alfa),math.cos(alfa),0],[0,0,0,1]])
rotation_y = np.array([[math.cos(alfa),0,math.sin(alfa),0],[0,1,0,0],
[-math.sin(alfa),0,math.cos(alfa),0],[0,0,0,1]])
rotation_z = np.array([[math.cos(alfa),math.sin(alfa),0,0],
[-math.sin(alfa),math.cos(alfa),0,0],
[0,0,1,0],[0,0,0,1]])
matrix = translation.dot(rotation_x.dot(rotation_y.dot(rotation_z)))
return matrix

def tube_function(RVElength,tubeR,volumefraction=0.03,tubelength_x=100):
'''
input: RVElength:length of the RVE cube,
volumefraction: the volume fraction of the composit.
tubelength_x: the length of cnt in x axis.
output: a list of points coordinates of the axis of cnt
array([x1 y1 z1][x2 y2 z2]...[xn yn zn])

call: trans_matrix().
'''
xdistance = 1 #the distance in x axis between adjacent two points
x = np.arange(-tubelength_x / 2, tubelength_x / 2 + 1, xdistance)
y = 10*np.sin(x / (tubelength_x/2) * np.pi)
z = np.zeros(len(x))
one = np.ones(len(x))
# combine the x,y,z and one into an array, and perform transpose
#array([[x1 y1 z1 1],
#       [x2 y2 z2 1],
#       ...........
#       [xn yn zn 1]])
xyz = np.array([np.reshape(x,len(x)),np.reshape(y,len(y)),
np.reshape(z,len(z)),np.reshape(one,len(one))]).transpose()
cntnum = 1 # the number or cnt
while True:
# process transform, and delete the last column of transformed xyz
pointdata = np.delete(xyz.dot(trans_matrix(a=random.uniform(-RVElength/2,RVElength/2),
b=random.uniform(-RVElength/2,RVElength/2),
c=random.uniform(-RVElength/2,RVElength/2),
alfa = random.uniform(0,360),
beta = random.uniform(0,360),
gama = random.uniform(0,360))),[3],axis=1)
if cntnum == 1:# the first cnt
# request the length of the axis of cnt (cntlength)
i = 0 # fist point in cube
cntlength = 0 #the length of the axis of cnt
# judge if the point in the cube
for coordinate in pointdata:
incube = (coordinate[0] <= RVElength/2 and coordinate[0] >= -RVElength/2
and coordinate[1] <= RVElength/2 and coordinate[1] >= -RVElength/2
and coordinate[2] <= RVElength/2 and coordinate[2] >= -RVElength/2)
if incube:
if i==0:# first point
old_coordinate = coordinate
i= i+1
continue
else:
cntlength = cntlength +np.linalg.norm(coordinate-old_coordinate)
old_corrdinate = coordinate
i = i + 1
else:
i = 0 # next point will be firs point
oldpointdata = pointdata
temp_volumefraction = cntlength*(np.pi*tubeR**2)/RVElength**3
if temp_volumefraction >= volumefraction:
break
cntnum = cntnum+1
continue
#if not the firs cnt, compare the distance between new and old cnt points
# coordinates
if min_distance(oldpointdata,pointdata) >= tubeR*2:
# request the length of the axis of cnt (cntlength)
i = 0 # fist point in cube
# judge if the point in the cube
for coordinate in pointdata:
incube = (coordinate[0] <= RVElength/2 and coordinate[0] >= -RVElength/2
and coordinate[1] <= RVElength/2 and coordinate[1] >= -RVElength/2
and coordinate[2] <= RVElength/2 and coordinate[2] >= -RVElength/2)
if incube:
if i==0:# first point
old_coordinate = coordinate
i= i+1
continue
else:
cntlength = cntlength +np.linalg.norm(coordinate-old_coordinate)
old_coordinate = coordinate
i = i + 1
else:
i = 0 # next point will be firs point
temp_volumefraction = cntlength*(np.pi*tubeR**2)/RVElength**3
if temp_volumefraction >= volumefraction:
break
oldpointdata = np.concatenate((oldpointdata,pointdata))
cntnum = cntnum + 1 # next cnt will not be the first
return oldpointdata

def min_distance(pointset1,pointset2):
'''
input: two sets of points coordinates.
array([x1 y1 z1],...[xn yn zn]) or [x y z]
output: the min distance between the two point sets.
'''
listdistance=[] # memory the distance data
# turn pointset data into np.array
pointset1 = np.array(pointset1)
pointset2 = np.array(pointset2)
if pointset1.ndim !=1:# not one dimension
for point1 in pointset1:
if pointset2.ndim !=1:
for point2 in pointset2:
listdistance.append(np.linalg.norm(point1-point2))
else:
point2 = pointset2 # pointset2 is one dimension
listdistance.append(np.linalg.norm(point1-point2))
else:
point1 = pointset1
if pointset2.ndim !=1:
for point2 in pointset2:
listdistance.append(np.linalg.norm(point1-point2))
else:
point2 = pointset2
listdistance.append(np.linalg.norm(point1-point2))
return min(listdistance)

def writeinp(mesh_data):
'''write data into inp file'''
inp = open ('inp.inp','w')
inp.write ("*Heading\n** Job name: 1 Model name: MMC\n")
inp.write ("** Generated by: Abaqus/CAE 6.12-3\n")
inp.write ("*Preprint,echo=NO, model=NO, history=NO, contact=NO\n")
inp.write ("**\n**PART\n**\n")

'''change the data structure of mesh_data'''
matrix_point_num = set()
matrix_element = []
reinforcement_point_num = set()
reinforcement_element = []
interface_point_num = set()
interface_element = []
for (e,t) in mesh_data.elementType.items():
if t == 1:
for point_num in mesh_data.element[e]:
matrix_element.append(mesh_data.element[e])
if t == 2:
for point_num in mesh_data.element[e]:
reinforcement_element.append(mesh_data.element[e])
if t == 3:
for point_num in mesh_data.element[e]:
interface_element.append(mesh_data.element[e])

'''write part matrix'''
inp.write("*Part, name=MATRIX"'\n')
inp.write("*Node"'\n')
for point_num in matrix_point_num:
inp.write(str(point_num)+",    "+str(mesh_data.point[point_num])[1:-1]+"\n")
inp.write("*Element, type=C3D8\n")
element_num=0
for element in matrix_element:
element_num = element_num + 1
inp.write(str(element_num)+",    "+str(element)[1:-1]+"\n")
inp.write("*End Part\n")

'''write part reinforcement'''
inp.write("*Part, name=REINFORCEMENT"'\n')
inp.write("*Node"'\n')
for point_num in reinforcement_point_num:
inp.write(str(point_num)+",    "+str(mesh_data.point[point_num])[1:-1]+"\n")
inp.write("*Element, type=C3D8\n")
element_num=0
for element in reinforcement_element:
element_num = element_num + 1
inp.write(str(element_num)+",    "+str(element)[1:-1]+"\n")
inp.write("*End Part\n")

'''
update Element()
'''
tubeR = 10 # the radius of cnt is 20
mesh_data = Element()
print ('Have got Element()')
pointdataset = tube_function(RVElength = 800, tubeR = tubeR)
print ('Have got pointdataset')
for pointdata in pointdataset:
# to judge whether the coordinate is in a cube whose center is
# pointdate and length is diameter of cnt
xmax = pointdata[0]+tubeR
xmin = pointdata[0]-tubeR
ymax = pointdata[1]+tubeR
ymin = pointdata[1]-tubeR
zmax = pointdata[2]+tubeR
zmin = pointdata[2]-tubeR
for (pointBLB, nouse) in mesh_data.elementType.items():
coordinate = mesh_data.point[pointBLB]# [x,y,z]
aroundpointdata=(coordinate[0] <= xmax and
coordinate[0] >= xmin and
coordinate[1] <= ymax and
coordinate[1] >= ymin and
coordinate[2] <= zmax and
coordinate[2] >= zmin)
if aroundpointdata:
if min_distance(coordinate,pointdata)<=tubeR:
mesh_data.elementType[pointBLB]=2
writeinp(mesh_data)
```

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.