Rectangle area element split in Python

Program overview

To help create test data for FEM programs that use quadrilateral elements, we created a rectangular area element division program. By inputting the dimensions and number of divisions of the rectangular area you want to model, the number of nodes, the number of elements, the element-node relationship, and the node coordinates are output. In order to confirm the created data, the model diagram is also output by matplotlib.

The article below is very helpful on how to use matplotlib patches. http://matthiaseisen.com/matplotlib/

Rectangle area element division program

py_rectmesh.py


# Meshing of rectangular domain
import numpy as np
import sys
import matplotlib.pyplot as plt
import matplotlib.patches as patches
#(Reference site) http://matthiaseisen.com/matplotlib/

param=sys.argv
aa=float(param[1]) # length in x-direction
bb=float(param[2]) # length in y-direction
nn=int(param[3])   # division number in x-direction
mm=int(param[4])   # division number in y-direction
x0=float(param[5]) # x-coordinate at left bottom of the domain
y0=float(param[6]) # y-coordinate at left bottom of the domain

npoin=(nn+1)*(mm+1)
nele=nn*mm
node=np.zeros([4,nele],dtype=np.int)
x=np.zeros([2,npoin],dtype=np.float64)

k=-1
for i in range(0,mm+1):
    for j in range(0,nn+1):
        k=k+1
        x[0,k]=aa/float(nn)*float(j)+x0
        x[1,k]=bb/float(mm)*float(i)+y0
ne=-1
for k in range(0,mm):
    i0=1+(nn+1)*k
    for j in range(0,nn):
        i=i0+j
        ne=ne+1
        node[0,ne]=i
        node[1,ne]=i+1
        node[2,ne]=node[1,ne]+nn+1
        node[3,ne]=node[2,ne]-1

print('{0:5d} {1:5d}'.format(npoin,nele))
for ne in range(0,nele):
    ss=''
    for i in range(0,4):
        ss=ss+'{0:5d}'.format(node[i,ne])
    ss=ss+'{0:5d}'.format(ne+1)
    print(ss)
for nd in range(0,npoin):
    ss=''
    for i in range(0,2):
        ss=ss+'{0:10.3f}'.format(x[i,nd])
    ss=ss+'{0:5d}'.format(nd+1)
    print(ss)

# Drawing
axmin=np.min(x[0,:])
axmax=np.max(x[0,:])
aymin=np.min(x[1,:])
aymax=np.max(x[1,:])
ra=np.min([axmax-axmin,aymax-aymin])
xmin=axmin-0.2*ra
xmax=axmax+0.2*ra
ymin=aymin-0.2*ra
ymax=aymax+0.2*ra

fnameF='fig_mesh.png'
fig = plt.figure()
ax1=plt.subplot(111)
ax1.set_xlim([xmin,xmax])
ax1.set_ylim([ymin,ymax])
ax1.set_xlabel('x-direction (m)')
ax1.set_ylabel('y-direction (m)')
ax1.spines['right'].set_visible(False)
ax1.spines['top'].set_visible(False)
ax1.yaxis.set_ticks_position('left')
ax1.xaxis.set_ticks_position('bottom')
ax1.set_aspect('equal', 'datalim')

# mesh
for ne in range(0,nele):
    n1=node[0,ne]-1; x1=x[0,n1]; y1=x[1,n1]
    n2=node[1,ne]-1; x2=x[0,n2]; y2=x[1,n2]
    n3=node[2,ne]-1; x3=x[0,n3]; y3=x[1,n3]
    n4=node[3,ne]-1; x4=x[0,n4]; y4=x[1,n4]
    ax1.fill([x1,x2,x3,x4],[y1,y2,y3,y4],facecolor='#ffffcc',edgecolor='#777777',lw=0.5)
# node number
for i in range(0,npoin):
    ax1.add_patch(patches.Circle((x[0,i],x[1,i]),0.05*ra,facecolor='#ffffcc',edgecolor='#777777',lw=0.5))
    ax1.text(x[0,i],x[1,i],str(i+1),ha='center',va='center',fontsize=12,color='#0000ff')
# element number
for ne in range(0,nele):
    n1=node[0,ne]-1; x1=x[0,n1]; y1=x[1,n1]
    n2=node[1,ne]-1; x2=x[0,n2]; y2=x[1,n2]
    n3=node[2,ne]-1; x3=x[0,n3]; y3=x[1,n3]
    n4=node[3,ne]-1; x4=x[0,n4]; y4=x[1,n4]
    x0=0.25*(x1+x2+x3+x4)
    y0=0.25*(y1+y2+y3+y4)
    ax1.text(x0,y0,str(ne+1),ha='center',va='center',fontsize=12,color='#ff0000')

ss='npoin='+str(npoin)+', nele='+str(nele)
ax1.set_title(ss)
plt.savefig(fnameF, bbox_inches="tight", pad_inches=0.2)
plt.show()

Output example

Script format for execution

python py_rectmesh.py aa bb nn mm x0 y0 > out.txt
aa : Length of the area in the x direction
bb : The length of the area in the y direction
nn : Number of divisions in the x direction of the area
mm : Number of divisions in the y direction of the area
x0 : x coordinate of the lower left corner of the area
y0 : y coordinate in the lower left corner of the area
out.txt : Output file name

Execution script case

python py_rectmesh.py 5 3 5 3 0 0 > _test.txt

Output data example

The output is the number of nodes, the number of elements, the element-node relationship, and the node coordinates. The element number is output at the end of the element-node relationship output column, which is reference data. Also, the node number is output at the end of the output column of node coordinates, which is also reference data.

   24    15                 # number of nodes, number of elements
    1    2    8    7    1   # element-node relationship
    2    3    9    8    2   # node-1, node-2, node-3, node-4, element_number(reference)
    3    4   10    9    3
    4    5   11   10    4
    5    6   12   11    5
    7    8   14   13    6
    8    9   15   14    7
    9   10   16   15    8
   10   11   17   16    9
   11   12   18   17   10
   13   14   20   19   11
   14   15   21   20   12
   15   16   22   21   13
   16   17   23   22   14
   17   18   24   23   15
     0.000     0.000    1  # x-coordinate, y-coordinate, node_number(reference)
     1.000     0.000    2
     2.000     0.000    3
     3.000     0.000    4
     4.000     0.000    5
     5.000     0.000    6
     0.000     1.000    7
     1.000     1.000    8
     2.000     1.000    9
     3.000     1.000   10
     4.000     1.000   11
     5.000     1.000   12
     0.000     2.000   13
     1.000     2.000   14
     2.000     2.000   15
     3.000     2.000   16
     4.000     2.000   17
     5.000     2.000   18
     0.000     3.000   19
     1.000     3.000   20
     2.000     3.000   21
     3.000     3.000   22
     4.000     3.000   23
     5.000     3.000   24

Output diagram example

fig_mesh.png

that's all

Recommended Posts

Rectangle area element split in Python
Split iterator into chunks in python
How to extract polygon area in Python
Manipulate namespaced XML in Python (Element Tree)
Quadtree in Python --2
Python in optimization
Metaprogramming in Python
Python 3.3 in Anaconda
Geocoding in python
SendKeys in Python
Meta-analysis in Python
Unittest in python
Epoch in Python
Discord in Python
Sudoku in Python
nCr in python
N-Gram in Python
Programming in python
Plink in Python
Lifegame in Python.
FizzBuzz in Python
Sqlite in python
N-gram in python
LINE-Bot [0] in Python
Csv in python
Disassemble in Python
Constant in python
nCr in Python.
format in python
Scons in Python3
Puyo Puyo in python
python in virtualenv
PPAP in Python
Quad-tree in Python
Reflection in Python
Chemistry in Python
Hashable in python
DirectLiNGAM in Python
LiNGAM in Python
Flatten in python
flatten in python
[Python] Display the Altair legend in the plot area
Split files when writing vim plugin in python
Split camel case string word by word in Python
Sorted list in Python
Daily AtCoder # 36 in Python
Daily AtCoder # 2 in Python
Implement Enigma in python
Daily AtCoder # 32 in Python
Daily AtCoder # 18 in Python
Singleton pattern in Python
File operations in Python
Key input in Python
Daily AtCoder # 33 in Python
Logistic distribution in Python
Daily AtCoder # 7 in Python
LU decomposition in Python
One liner in Python
Daily AtCoder # 24 in Python
case class in python
RNN implementation in python