<jupyter_start><jupyter_text>生成网格<jupyter_code>with pygmsh.geo.Geometry() as geom:
    geom.add_polygon(
        [
            [0.0, 0.0],
            [1.0, 0.0],
            [1.0, 1.0],
            [0.0, 1.0],
        ],
        mesh_size=0.5,
    )
    mesh = geom.generate_mesh()<jupyter_output><empty_output><jupyter_text>展示节点数据<jupyter_code>mesh.point_data
nodes=mesh.points
nodes
plt.scatter(mesh.points.T[0],mesh.points.T[1])<jupyter_output><empty_output><jupyter_text>展示网格<jupyter_code>mesh.get_cell_data
mesh.cells_dict['triangle']
elements=mesh.cells_dict['triangle']
E=10.0
nv=0.3
D=E/(1.0-nv**2)*np.array([
    [1.0,nv,0.0], 
    [nv,1.0,0.0], 
    [0.0,0.0,0.5*(1.0-nv)]])
len_global=2*len(nodes)
K=np.zeros((len_global,len_global))
for element in elements:
    x0=nodes[element[0]][0]
    y0=nodes[element[0]][1]
    x1=nodes[element[1]][0]
    y1=nodes[element[1]][1]
    x2=nodes[element[2]][0]
    y2=nodes[element[2]][1]
    area_2=np.linalg.det(np.array([
        [1.0,x0,y0], 
        [1.0,x1,y1], 
        [1.0,x2,y2]]))
    B=0.5/area_2*np.array([
        [y1-y2,0.0,y2-y0,0.0,y0-y1,0.0], 
        [0.0,x2-x1,0.0,x0-x2,0.0,x1-x0], 
        [x2-x1,y1-y2,x0-x2,y2-y0,x1-x0,y0-y1]])
    K_element=area_2*B.T@D@B
    deform_global_index=np.array([[2*element[i],2*element[i]+1] for i in range(3)],dtype=np.uint64).reshape(-1)
    for i_local,i_global in enumerate(deform_global_index):
        for j_local,j_global in enumerate(deform_global_index):
            K[i_global,j_global]+=K_element[i_local,j_local]
def get_position_local(element):
    x0=nodes[element[0]][0]
    y0=nodes[element[0]][1]
    x1=nodes[element[1]][0]
    y1=nodes[element[1]][1]
    x2=nodes[element[2]][0]
    y2=nodes[element[2]][1]
    return x0,y0,x1,y1,x2,y2
def det_B(x0,y0,x1,y1,x2,y2):
    matrix=np.array([\
        [1.0,x0,y0], 
        [1.0,x1,y1], 
        [1.0,x2,y2]])
    det=np.linalg.det(matrix)
    return abs(det)
def b_matrix(x0,y0,x1,y1,x2,y2,area_2):
    return 1/area_2*np.array([\
        [y1-y2,0.0,y2-y0,0.0,y0-y1,0.0], 
        [0.0,x2-x1,0.0,x0-x2,0.0,x1-x0], 
        [x2-x1,y1-y2,x0-x2,y2-y0,x1-x0,y0-y1]])
def element_integrate(B,D,area_2):
    return area_2*B.T@D@B
def Ke2K(element,K_element,K):
    deform_global_index=np.array([[2*element[i],2*element[i]+1] for i in range(3)],dtype=np.uint64).reshape(-1)
    for i_local,i_global in enumerate(deform_global_index):
        for j_local,j_global in enumerate(deform_global_index):
            K[i_global,j_global]+=K_element[i_local,j_local]
for element in elements:
    # calculate stiffness of element
    x0,y0,x1,y1,x2,y2=get_position_local(element)
    area_2=det_B(x0,y0,x1,y1,x2,y2)
    B=b_matrix(x0,y0,x1,y1,x2,y2,area_2)
    K_element=element_integrate(B,D,area_2)
    # add local Ke into global K
    Ke2K(element,K_element,K)
# interface
x_fix={}
y_fix={}
f_given={}
for index,position in enumerate(nodes):
    x=position[0].tolist()
    if x<1e-6:
        x_fix.update({index:0.0})
        y_fix.update({index:0.0})
    if 1.0-x<1e-6:
        f_given.update({index:(1.0,0.0)})
deform_free_index=[]
for node,_ in enumerate(nodes):
    if node not in x_fix:
        deform_free_index.append(2*node)
    if node not in y_fix:
        deform_free_index.append(2*node+1)
len_reduce=len(deform_free_index)
deform=np.zeros(len_global)
force=np.zeros(len_global)
for node in x_fix:
    deform[2*node]=x_fix[node]
for node in y_fix:
    deform[2*node+1]=y_fix[node]
for node in f_given:
    force[2*node]=f_given[node][0]
    force[2*node+1]=f_given[node][1]
deform_reduce=np.zeros(len_reduce)
force_reduce=np.zeros(len_reduce)
K_reduce=np.zeros((len_reduce,len_reduce))
for i_reduce,i_global in enumerate(deform_free_index):
    for j_reduce,j_global in enumerate(deform_free_index):
        K_reduce[i_reduce,j_reduce]=K[i_global,j_global]
for i_reduce,i_global in enumerate(deform_free_index):
    force_reduce[i_reduce]=force[i_global]
deform_reduce=np.linalg.solve(K_reduce,force_reduce)
for i_reduce,i_global in enumerate(deform_free_index):
    deform[i_global]=deform_reduce[i_reduce]
force=K@deform
deform_free_index
deform.reshape(len_global//2,2)
def get_strain(deform_local):
    return B@deform_local
def get_stress(deform_local):
    return D@B@deform_local
def get_deform(element):
    n0=element[0].tolist()
    n1=element[1].tolist()
    n2=element[2].tolist()
    u0=deform[2*n0]
    v0=deform[2*n0+1]
    u1=deform[2*n1]
    v1=deform[2*n1+1]
    u2=deform[2*n2]
    v2=deform[2*n2+1]
    return u0,v0,u1,v1,u2,v2
def center_point(x0,y0,x1,y1,x2,y2):
    x=(x0+x1+x2)/3.0
    y=(y0+y1+y2)/3.0
    return [x,y]
show_data={'position':{'x':[],'y':[]},'strain':{'e11':[],'e22':[],'e12':[]},'stress':{'S11':[],'S22':[],'S12':[]}}
for element in elements:
    x0,y0,x1,y1,x2,y2=get_position_local(element)
    position_element=center_point(x0,y0,x1,y1,x2,y2)
    deform_local=get_deform(element)
    strain=get_strain(deform_local)
    stress=get_stress(deform_local)
    show_data['position']['x'].append(position_element[0])
    show_data['position']['y'].append(position_element[1])
    for i,name in enumerate(('e11','e22','e12')):
        show_data['strain'][name].append(strain[i])
    for i,name in enumerate(('S11','S22','S12')):
        show_data['stress'][name].append(stress[i])
points=np.array([show_data['position']['x'],show_data['position']['y']]).T
values=show_data['stress']['S11']
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:100j]
z=griddata(points, values, (grid_x, grid_y), method='nearest')

fig, ax = plt.subplots()
ax.pcolormesh(grid_x,grid_y,z)<jupyter_output>C:\Users\MuTong\AppData\Local\Temp/ipykernel_14592/3417120346.py:7: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  ax.pcolormesh(grid_x,grid_y,z)
{'fraction_non_alphanumeric': 0.18169672678690715, 'fraction_numerical': 0.06462925851703406, 'mean_word_length': 4.623474178403756, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '638594', 'license': 'MIT', 'n_tokens_mistral': 2815, 'n_tokens_neox': 2744, 'n_words': 315, 'path': 'scratch_po.ipynb', 'repo_name': 'MuTong320/NaiveFeaDocument', 'size': '5988'}