<jupyter_start><jupyter_text>Example of central tension  Import packages<jupyter_code>import numpy as np
import torch
import torch.nn as nn
import matplotlib.pyplot as plt 
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D<jupyter_output><empty_output><jupyter_text>User defined neural network* A fully-connected feed-forward network    * **n_input** - dimension of input, 2 in this case(x and y lcation)    * **n_output** - dimension of output, 1 in this case (u - horizontal or v - vertical displacemnt)    * **n_layer** - number of hidden layers    * **n_nodes** - number of nodes of each hidden layer* **two networks are defined seperately for u and v**<jupyter_code>class Net(nn.Module):

    def __init__(self, n_input, n_output, n_layer, n_nodes):
        super(Net, self).__init__()
        self.n_layer = n_layer
        
        self.Input = nn.Linear(n_input, n_nodes)   # linear layer
        nn.init.xavier_uniform_(self.Input.weight) # wigths and bias initiation
        nn.init.normal_(self.Input.bias)

        self.Output = nn.Linear(n_nodes, n_output)
        nn.init.xavier_uniform_(self.Output.weight)
        nn.init.normal_(self.Output.bias)
        
        self.Hidden = nn.ModuleList() # hidden layer list
        for i in range(n_layer):
            self.Hidden.append(nn.Linear(n_nodes, n_nodes))
        for layer in self.Hidden:
            nn.init.xavier_uniform_(layer.weight)
            nn.init.normal_(layer.bias)
        

    def forward(self, x):
        y = torch.tanh(self.Input(x)) # tanh activation function
        for layer in self.Hidden:
            y = torch.tanh(layer(y))
        y = self.Output(y)
        return y<jupyter_output><empty_output><jupyter_text>Functions for first and sencod partial derivtives<jupyter_code>def derivativs(x, Net, func):
    
    w = Net(x)*func(x).view(-1,1)
    #print(Net(x),func(x).view(-1,1),w)
    dw_xy = torch.autograd.grad(w, x, torch.ones_like(w), retain_graph=True,
    create_graph=True, allow_unused=True)
    dw_x = dw_xy[0][:,0].view(-1,1)
    dw_y = dw_xy[0][:,1].view(-1,1)
    # second partial derivative
    dw_xxy = torch.autograd.grad(dw_x, x, torch.ones_like(dw_x), retain_graph=True,
    create_graph=True,allow_unused=True)
    dw_xx = dw_xxy[0][:,0].view(-1,1)
    dw_xy = dw_xxy[0][:,1].view(-1,1)
    dw_yy = torch.autograd.grad(dw_y, x, torch.ones_like(dw_y), retain_graph=True,
    create_graph=True,allow_unused=True)[0][:,1].view(-1,1)

    return w, dw_x, dw_y, dw_xx, dw_yy, dw_xy
def derivativs2(x, Net, func):
    
    w = Net(x)*func(x).view(-1,1)
    #print(Net(x),func(x).view(-1,1),w)
    dw_xy = torch.autograd.grad(w, x, torch.ones_like(w), retain_graph=True,
    create_graph=True, allow_unused=True)
    dw_x = dw_xy[0][:,0].view(-1,1)
    dw_y = dw_xy[0][:,1].view(-1,1)

    return w, dw_x, dw_y
def derivativs2u(x, Net, func):
    
    w = Net(x)*func(x).view(-1,1) + 0.02*(x[:,0]+10).view(-1,1)
    #print(Net(x),func(x).view(-1,1),w)
    dw_xy = torch.autograd.grad(w, x, torch.ones_like(w), retain_graph=True,
    create_graph=True, allow_unused=True)
    dw_x = dw_xy[0][:,0].view(-1,1)
    dw_y = dw_xy[0][:,1].view(-1,1)

    return w, dw_x, dw_y<jupyter_output><empty_output><jupyter_text>Construct the loss of PDE and boundary condition<jupyter_code>def PDE(x, Net_u, func_u, Net_v, func_v):
    
    E = 70
    mu = 0.3
    a = 10
    b = 10
    C = E/(1 - mu**2)
    
    _, du_x, du_y = derivativs2(x, Net_u, func_u) 
    _, dv_x, dv_y = derivativs2(x, Net_v, func_v)
    
    du_y, dv_x = du_y*a/b, dv_x*b/a
    
    eps_xx = du_x 
    eps_yy = dv_y 
    eps_xy = 0.5*(du_y + dv_x)
    
    sig_xx = C*(eps_xx + mu*eps_yy)
    sig_yy = C*(eps_yy + mu*eps_xx)
    sig_xy = C*(1 - mu)*eps_xy

    res = 0.5*(eps_xx*sig_xx + eps_yy*sig_yy + 2*eps_xy*sig_xy)

    return torch.mean(res)<jupyter_output><empty_output><jupyter_text>Function for preparing training data* **uniform sampling at the boundary and radom sampling within the solution domain**<jupyter_code>def train_data(Nx, Ny, Nc, Nf1, Nf2):
    
    a = 10
    b = 10
    r = 2
    
    x0 = 1
    y0 = 1
    #r0 = r
    ra = 2/a
    rb = 2/b
    
    c = 3./a

    #X2 = np.hstack([ x0*np.ones([Nx,1]), np.random.rand(Nx,1)*y0])
    X2 = np.hstack([ x0*np.ones([Nx,1]), np.linspace(0,y0,Nx).reshape([Nx,1])])
    
    Xb = torch.tensor(X2, dtype=torch.float32, requires_grad=True)

    Y2 = np.hstack([np.random.rand(Nx,1)*x0,  y0*np.ones([Ny,1])])
    Yb = torch.tensor(Y2, dtype=torch.float32, requires_grad=True)
    
    Xf = np.random.rand(Nf1, 2)
    ind1 = []
    #Xf = np.copy(Xf0)
    for i, (x, y) in enumerate(Xf):

        #x, y = x*x0, y*y0
        #if (x >= c) or (y >= c):
        if (x**2/c**2 + y**2/c**2) >= 1:
            ind1.append(i)  
    Xf1 = Xf[ind1]
    
    Xf = np.random.rand(Nf2, 2)* c
    ind2 = []
    for i, (x, y) in enumerate(Xf):
        if (x**2/ra**2 + y**2/rb**2) >= 1 and (x**2/c**2 + y**2/c**2) < 1 and (x <= y):
            ind2.append(i)
    Xf2 = Xf[ind2]
    
    Xf3 = np.vstack((Xf2[:,1], Xf2[:,0])).T
    
    Xf1 = torch.tensor(Xf1, dtype=torch.float32, requires_grad=True)
    Xf2 = torch.tensor(Xf2, dtype=torch.float32, requires_grad=True)
    Xf3 = torch.tensor(Xf3, dtype=torch.float32, requires_grad=True)
    
    return Xb, Yb, Xf1, Xf2, Xf3
def train_data2(Nx, Ny, Nc, Nf0, Nf1, Nf2):
    
    a = 10
    b = 10
    r = 2
    
    x0 = 1
    y0 = 1
    #r0 = r
    ra = 2/a
    rb = 2/b
    
    c = 3.0/a
    d = 5.0/a

    X2 = np.hstack([ x0*np.ones([Nx,1]), np.random.rand(Nx,1)*y0])
    
    Xb = torch.tensor(X2, dtype=torch.float32, requires_grad=True)

    Y2 = np.hstack([np.random.rand(Nx,1)*x0,  y0*np.ones([Ny,1])])
    Yb = torch.tensor(Y2, dtype=torch.float32, requires_grad=True)
    
    Xf = np.random.rand(Nf0, 2)
    ind0 = []
    #Xf = np.copy(Xf0)
    for i, (x, y) in enumerate(Xf):

        #x, y = x*x0, y*y0
        #if (x >= c) or (y >= c):
        if (x**2/d**2 + y**2/d**2) >= 1:
            ind0.append(i)  
    Xf0 = Xf[ind0]
    
    Xf = np.random.rand(Nf1, 2)
    ind1 = []
    #Xf = np.copy(Xf0)
    for i, (x, y) in enumerate(Xf):

        #x, y = x*x0, y*y0
        #if (x >= c) or (y >= c):
        if (x**2/c**2 + y**2/c**2) >= 1 and (x**2/d**2 + y**2/d**2) <= 1:
            ind1.append(i)  
    Xf1 = Xf[ind1]
    
    Xf = np.random.rand(Nf2, 2)* c
    ind2 = []
    for i, (x, y) in enumerate(Xf):
        if (x**2/ra**2 + y**2/rb**2) >= 1 and (x <= y):
            ind2.append(i)
    Xf2 = Xf[ind2]
    
    Xf3 = np.vstack((Xf2[:,1], Xf2[:,0])).T
   
    Xf0 = torch.tensor(Xf0, dtype=torch.float32, requires_grad=True)
    Xf1 = torch.tensor(Xf1, dtype=torch.float32, requires_grad=True)
    Xf2 = torch.tensor(Xf2, dtype=torch.float32, requires_grad=True)
    Xf3 = torch.tensor(Xf3, dtype=torch.float32, requires_grad=True)
    
    return Xb, Yb, Xf0, Xf1, Xf2, Xf3
Xf2
# Prepare training data
Nx = 1000
Ny = 1000
Nc = 250

Nf1 = 8000
Nf2 = 4000
Xb1, Xb2, Xf1, Xf2, Xf3 = train_data(Nx, Ny, Nc, Nf1, Nf2)

#print(Xb1, Xb2, Xb3, Xf)
# Construct neural network
Net_u = Net(2, 1, 10, 20)
Net_v = Net(2, 1, 10, 20)
func_u = lambda x:  x[:,0]
func_v = lambda x:  x[:,1]
Xf2.size()
def plot_pos(Xf0, Xf1, Xf2, Xb1):
    
    Xf0 = Xf0.detach().numpy()
    Xf1 = Xf1.detach().numpy()
    Xf2 = Xf2.detach().numpy()
    #Xf3 = Xf3.detach().numpy()
    
    Xb1 = Xb1.detach().numpy()
    #Xb2 = Xb2.detach().numpy()
    plt.rcParams['font.family'] = 'Times New Roman'
    fig, ax = plt.subplots(figsize=(5.8, 4.8)) 
    surf = ax.scatter(Xf0[:,0], Xf0[:,1], s=0.5, c='b', marker='x')
    plt.scatter(Xf1[:,0], Xf1[:,1], s=0.5, c='r', marker='x')
    plt.scatter(Xf2[:,0], Xf2[:,1], s=0.5, c='r', marker='x')
    #plt.scatter(Xf3[:,0], Xf3[:,1], s=1, c='b', marker='x')
    
    #plt.plot([0,0], [0,1], c='k', alpha=0.5)
    #plt.plot([0,1], [0,0], c='k', alpha=0.5)
    #plt.plot([1,1], [0,1], c='k', alpha=0.5)
    #plt.plot([0,1], [1,1], c='k', alpha=0.5)
    plt.scatter(Xb1[:,0], Xb1[:,1], s=0.1, c='k', marker='x')
    #plt.scatter(Xb2[:,0], Xb2[:,1], s=0.1, c='k', marker='x')
    #cbar = fig.colorbar(ax)
    #cb = fig.colorbar(surf, ticks=[0,0.03,0.06,0.09,0.12,0.15,0.18])
    #cb.ax.locator_params(nbins=7)
    #cb.ax.tick_params(labelsize=16)
    #cb.set_label(label =r'$\sigma_xx (MPa)$', fontsize=16)
    #cb.set_label(fontsize=16)
    ax.axis('equal')
    ax.set_xlabel('X Position (mm)', fontsize=18)
    ax.set_ylabel('Y Position (mm)', fontsize=18)
    for tick in ax.get_xticklabels():
        #tick.set_fontname('Times New Roman')
        tick.set_fontsize(16)
    for tick in ax.get_yticklabels():
        #tick.set_fontname('Times New Roman')
        tick.set_fontsize(16)
    #plt.savefig('Ch-R2-Sample_Uni.png', dpi=600, transparent=True)
    plt.show()
plot_pos(Xf1, Xf2, Xf3, Xb1)
import time
start_time = time.time()
# Construct neural network
# optimizer
nepoches = 6000
learning_rate = 1.0e-5
optimizer = torch.optim.Adam(list(Net_u.parameters())+list(Net_v.parameters()), lr=learning_rate)

training_loss = []
for epoch in range(nepoches):
    ## Calculate loss
    a = 10
    b = 10
    q = 1.0
    ra = 2
    rb = 2
    c = 3
    
    Eint1 = PDE(Xf1, Net_u, func_u, Net_v, func_v)*(a*b - 0.25*3.1415926*c*c)
    #Eint1 = PDE(Xf1, Net_u, func_u, Net_v, func_v)*0.25*3.1415926*(d*d - c*c)
    Eint2 = PDE(Xf2, Net_u, func_u, Net_v, func_v)*0.25*3.1415926*(c*c - ra*rb)/2
    Eint3 = PDE(Xf3, Net_u, func_u, Net_v, func_v)*0.25*3.1415926*(c*c - ra*rb)/2
    Eint =  Eint1 + Eint2 + Eint3
    Eext = torch.mean(q *Net_u(Xb1)*func_u(Xb1).view(-1,1))*a*b
    #Eint, Eext = Eint, Eext
    loss = Eint - Eext
    #loss = torch.abs(loss)
    loss.backward()
    #ind = 0
    #loss_list = [loss1, loss2 + loss3 + loss4, loss5 + loss6, loss7 + loss8]
    loss_temp = loss
    training_loss.append(loss_temp.detach().numpy())
    if (epoch+1) % 200 == 0:
        Xb1, Xb2, Xf1, Xf2, Xf3 = train_data(Nx, Ny, Nc, Nf1, Nf2)
    
    optimizer.step()
    optimizer.zero_grad()
    if (epoch+1) % 100 == 0:
        print(f'epoch:{epoch+1}: total loss:{loss:.4e}, Eint:{Eint:.4e}, Eext:{Eext:.4e} ')
elapsed_time = time.time() - start_time
print(elapsed_time)
import pandas as pd
data = pd.read_csv('CH-R2.csv')
X_train = data.iloc[:, 5:7].to_numpy()
U = data.iloc[:,11:13].to_numpy()
LE = data.iloc[:,13:16].to_numpy()
sig = data.iloc[:,16:19].to_numpy()
U1 = U[:,0].reshape(-1, 1)
U2 = U[:,1].reshape(-1, 1)
eps11 = LE[:,0].reshape(-1, 1)
eps22 = LE[:,1].reshape(-1, 1)
eps12 = LE[:,2].reshape(-1, 1)
sig11 = sig[:,0].reshape(-1, 1)
sig22 = sig[:,1].reshape(-1, 1)
sig12 = sig[:,2].reshape(-1, 1)
X_train = torch.tensor(X_train, dtype=torch.float32, requires_grad=True)
X_train = X_train/10
U = torch.tensor(U, dtype=torch.float32)
U = U / 10
u1 = Net_u(X_train)*func_u(X_train).view(-1,1)
u2 = Net_v(X_train)*func_v(X_train).view(-1,1)

u1 = u1.detach().numpy().reshape(-1,1)*10
u2 = u2.detach().numpy().reshape(-1,1)*10

X = X_train[:,0].detach().numpy().reshape(-1,1)*10
Y = X_train[:,1].detach().numpy().reshape(-1,1)*10
np.corrcoef(U1.reshape(1,-1), u1.reshape(1,-1))
np.corrcoef(U2.reshape(1,-1), u2.reshape(1,-1))
# Then, "ALWAYS use sans-serif fonts"
plt.rcParams['font.family'] = 'Times New Roman'
fig, ax = plt.subplots(figsize=(5.8, 4.8)) 
surf = ax.scatter(X, Y, c = u1, vmin=0, vmax=0.18, cmap=cm.rainbow)
#cbar = fig.colorbar(ax)
cb = fig.colorbar(surf, ticks=[0,0.03,0.06,0.09,0.12,0.15,0.18])
#cb.ax.locator_params(nbins=7)
cb.ax.tick_params(labelsize=16)
#cb.set_label(label =r'$\sigma_xx (MPa)$', fontsize=16)
#cb.set_label(fontsize=16)
ax.axis('equal')
ax.set_xlabel('X Position (mm)', fontsize=18)
ax.set_ylabel('Y Position (mm)', fontsize=18)
for tick in ax.get_xticklabels():
    #tick.set_fontname('Times New Roman')
    tick.set_fontsize(16)
for tick in ax.get_yticklabels():
    #tick.set_fontname('Times New Roman')
    tick.set_fontsize(16)
#plt.savefig('CH-R2-U-Energy-10-10.png', dpi=600, transparent=True)
plt.show()
_, du_x, du_y,= derivativs2(X_train, Net_u, func_u)
_, dv_x, dv_y,= derivativs2(X_train, Net_v, func_v)

E, mu = 70, 0.3
sig_x = (du_x + mu*dv_y)*E/(1 - mu**2) 
sig_y = (dv_y + mu*du_x)*E/(1 - mu**2)
sig_xy = (dv_x + du_y)*E/(1 + mu)/2.

#X = Xf[:,0].detach().numpy().reshape(-1,1)
#Y = Xf[:,1].detach().numpy().reshape(-1,1)

sig_x = sig_x.detach().numpy()
sig_y = sig_y.detach().numpy() 
sig_xy = sig_xy.detach().numpy()

#err = abs(sig_x - sig11)
max(sig11)
plt.rcParams['font.family'] = 'Times New Roman'
fig, ax = plt.subplots(figsize=(5.8, 4.8)) 
#surf = ax.scatter(X, Y, c = abs(sig_x-sig11), vmin=0, vmax=0.1, cmap=cm.rainbow)
surf = ax.scatter(X, Y, c = sig_x, vmin=0, vmax=3, cmap=cm.rainbow)
#cbar = fig.colorbar(ax)
cb = fig.colorbar(surf)
cb.ax.locator_params(nbins=7)
cb.ax.tick_params(labelsize=16)
#cb.set_label(label =r'$\sigma_xx (MPa)$', fontsize=16)
#cb.set_label(fontsize=16)
ax.axis('equal')
ax.set_xlabel('X Position (mm)', fontsize=18)
ax.set_ylabel('Y Position (mm)', fontsize=18)
for tick in ax.get_xticklabels():
    #tick.set_fontname('Times New Roman')
    tick.set_fontsize(16)
for tick in ax.get_yticklabels():
    #tick.set_fontname('Times New Roman')
    tick.set_fontsize(16)
#plt.savefig('CH-R5-S-Energy-5-5-2.png', dpi=600, transparent=True)
plt.show()
np.corrcoef(sig22.reshape(1,-1), sig_y.reshape(1,-1))
np.corrcoef(sig11.reshape(1,-1), sig_x.reshape(1,-1))
np.corrcoef(sig12.reshape(1,-1), sig_xy.reshape(1,-1))
data_out = np.hstack([X.reshape(-1,1),Y.reshape(-1,1),u1,u2,sig_x,sig_y,sig_xy])
df_out = pd.DataFrame(data_out, columns=['X', 'Y', 'U', 'V', 'Sig_x', 'Sig_y', 'Sig_xy'])
df_out.to_csv('CH-R2-Energy-5-5-4.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19478998975559783, 'fraction_numerical': 0.05722230352700132, 'mean_word_length': 3.2883589582679638, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 32, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '550100', 'license': 'MIT', 'n_tokens_mistral': 6479, 'n_tokens_neox': 6249, 'n_words': 1120, 'path': 'Ex2 - Centrol-hole tension/CH_tension_Energy.ipynb', 'repo_name': 'weili101/Deep_Plates', 'size': '13666'}