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 
    ---------------radian--------------------------------
    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, 
           tubeR: the radius of cnt,  
           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')
    # write inp Head
    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_point_num.add(point_num)
            matrix_element.append(mesh_data.element[e])           
        if t == 2:
            for point_num in mesh_data.element[e]:
                reinforcement_point_num.add(point_num)
            reinforcement_element.append(mesh_data.element[e])
        if t == 3:
            for point_num in mesh_data.element[e]:
                interface_point_num.add(point_num)
            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.

Delete this entry (admin only).