<jupyter_start><jupyter_text>Energy-Based Neural Network  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
import time<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 derivative(x, Net, func):
    
    w = Net(x)*func(x).view(-1,1)
    
    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):
    
    
    _, du_x, du_y = derivative(x, Net_u, func_u) 
    _, dv_x, dv_y = derivative(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, Nf):
    
    x0 = 1
    y0 = 1
    
    #X2 = np.hstack([ x0*np.ones([Nx,1]), np.random.rand(Nx,1)*y0])
    Xb = np.hstack([ x0*np.ones([Nx,1]), y0*np.linspace(0,1,Nx).reshape(-1,1)])
    Xb = torch.tensor(Xb, dtype=torch.float32, requires_grad=True)

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

    
    Xf = np.random.rand(Nf,2)
    Xf = torch.tensor(Xf, dtype=torch.float32, requires_grad=True)
    
    return Xb, Yb, Xf
def err(X, Y):
    
    return torch.mean(torch.mean((X-Y)**2))
import pandas as pd
data = pd.read_csv('rec-2.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
# Prepare training data
Nx = 1000
Ny = 1000
Nf = 10000
Xb1, Xb2, Xf = train_data(Nx, Ny, Nf)

#print(Xb1, Xb2, Xb3, Xf)
# Construct neural network
Net_u = Net(2, 1, 5, 5)
Net_v = Net(2, 1, 5, 5)
func_u = lambda x:  x[:,0]
func_v = lambda x:  x[:,1]
E = 70
mu = 0.3
a = 10
b = 10
q0 = 1.0
C = E/(1 - mu**2)
# Construct neural network
# optimizer
nepoches = 6000
learning_rate = 1.0e-3
optimizer = torch.optim.Adam(list(Net_u.parameters())+list(Net_v.parameters()), lr=learning_rate)

training_loss = []
sig_loss = []

for epoch in range(nepoches):
    
    ## Calculate loss
    Eint = PDE(Xf, Net_u, func_u, Net_v, func_v)

    q = q0 * torch.cos(3.1415/2*Xb1[:,1]).view(-1,1)
    Eext = torch.mean(q *Net_u(Xb1)*func_u(Xb1).view(-1,1))

    loss = Eint - Eext
    
    loss.backward()

    if (epoch+1) % 200 == 0:
        Xb1, Xb2, Xf = train_data(Nx, Ny, Nf)
    
    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} ')
        
    if (epoch+1) % 10 ==0: 
        with torch.no_grad():
            U_pred = Net_u(X_train)*func_u(X_train).view(-1,1)
            V_pred = Net_v(X_train)*func_v(X_train).view(-1,1)
            loss = err(torch.cat((U_pred, V_pred), 1), U)
            training_loss.append(loss)
        
        _, du_x, du_y,= derivative(X_train, Net_u, func_u)
        _, dv_x, dv_y,= derivative(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.
 
        sig11_loss = err(sig_x, torch.tensor(sig11))
        sig22_loss = err(sig_y, torch.tensor(sig22))
        sig12_loss = err(sig_xy, torch.tensor(sig12))
        sig_loss.append([sig11_loss, sig22_loss, sig12_loss]) 
torch.save({'Net_u':Net_u.state_dict(), 'Net_v': Net_v.state_dict() }, 'energy-5-5.pt')
xx = [[xx[0].detach().numpy (),xx[1].detach().numpy(),xx[2].detach().numpy()] for xx in sig_loss1]
xx = [xx.detach().numpy() for xx in training_loss1]
pd.DataFrame(xx).to_clipboard()
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))
# 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.14, 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('Flat-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.

sig_x = sig_x.detach().numpy()
sig_y = sig_y.detach().numpy() 
sig_xy = sig_xy.detach().numpy()
np.corrcoef(sig12.reshape(1,-1), sig_xy.reshape(1,-1))
max(sig_x)<jupyter_output><empty_output><jupyter_text>Post-processing<jupyter_code>plt.rcParams['font.family'] = 'Times New Roman'
fig, ax = plt.subplots(figsize=(5.8, 4.8)) 
surf = ax.scatter(X, Y, c = sig_x, vmin=0, vmax=1.0, 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('Flat-S-Energy-10-10.png', dpi=600, transparent=True)
plt.show()
data_out = np.hstack([X.reshape(-1,1),Y.reshape(-1,1),u1,u2,sig_x,sig_y,sig_xy])
data_out = np.hstack([X.reshape(-1,1),Y.reshape(-1,1),U1,U2,sig11,sig22,sig12])
df_out = pd.DataFrame(data_out, columns=['X', 'Y', 'U', 'V', 'Sig_x', 'Sig_y', 'Sig_xy'])
df_out.to_csv('Flat-Energy-10-10-1.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17633674630261661, 'fraction_numerical': 0.04175199089874858, 'mean_word_length': 3.754461871281774, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '386253', 'license': 'MIT', 'n_tokens_mistral': 3910, 'n_tokens_neox': 3756, 'n_words': 686, 'path': 'Ex1 - Uniaxial tension/Flat_tension_Energy_based_NN.ipynb', 'repo_name': 'weili101/Deep_Plates', 'size': '8790'}