# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: fenicsker
#     language: python
#     name: fenicsker
# ---

#import the required libraries
import numpy as np
import matplotlib.pyplot as pl
from dolfin import *
import meshio

# +
## Import the mesh file
msh = meshio.read("meshings/slot_rbf2.msh")

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
                         cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write("plate.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)

from dolfin import *
set_log_level(LogLevel.ERROR)

mesh = Mesh()
with XDMFFile("plate.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
# define the facet meshes
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("plate.xdmf") as infile:
    infile.read(mvc2, "name_to_read")
#define the cells meshes    
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)

# +
from dolfin import *
# Scaled variables for cantilever beam
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
V = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
tol = 1E-14


bc = DirichletBC(V, Constant((0, 0, 0)),mf, 1)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)


ds = Measure("ds", domain=mesh, subdomain_data=mf)
dS = Measure("dS", domain=mesh, subdomain_data=mf)

# Compute solution
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds(2)

u = Function(V)
solve(a == L, u, bc)

file1= XDMFFile('results/u1.xdmf')
file1.parameters["flush_output"] = True
file1.write(u)
# # Plot solution
# plot(u, title='Displacement', mode='displacement')


# +
# define a degree of freedom map
v2d = vertex_to_dof_map(V)
dofs = []
for facet in facets(mesh):
    if mf[facet.index()] == 1:
        vertices = facet.entities(0)
        for vertex in vertices:
            dofs.append(v2d[vertex])

unique_dofs = np.array(list(set(dofs)), dtype=np.int32)
boundary_coords = V.tabulate_dof_coordinates()[unique_dofs]
N=len(boundary_coords[:,1])
g=np.zeros((N,3))
for i, dof in enumerate(unique_dofs):
    g[i,0]=u.vector()[dof]
#     print(boundary_coords[i], v.vector()[dof])

# -

g.shape

# define coordinates and control points 
coor=V.tabulate_dof_coordinates()
num_coor=len(coor[:,1])
ctr_pts=np.array(boundary_coords)

# +
# nx, ny, nz = (10, 10, 10)
# mesh = np.zeros((nx * ny * nz, 3))

# xv = np.linspace(0, 1, nx)
# yv = np.linspace(0, 1, ny)
# zv = np.linspace(0, 1, nz)
# z_1, y_1, x_1 = np.meshgrid(zv, yv, xv)

# # mesh = np.array([x.ravel(), y.ravel(), z.ravel()])
# # mesh = mesh.T

# +
# num_coor=len(x_1.flatten())
# coor=np.zeros((num_coor,3))
# coor[:,0]=x_1.flatten()
# coor[:,1]=y_1.flatten()
# coor[:,2]=z_1.flatten()
# -

#Plot in matplotlib
fig = pl.figure(1)
ax = fig.add_subplot(111, projection='3d')
ax.scatter(coor[:, 0], coor[:, 1], coor[:, 2], c='blue', marker='o')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
pl.show()

# +
# N=len(ctr_pts)
# g=np.zeros((N,3))
# for i in range(N):
#     g[i,0]=
# -

# define the P matrix from control points 
P=np.zeros((N,4))
for i in range(N):
    for j in range(4):
        if (j==0): P[i,j]=1
        else: P[i,j]=ctr_pts[i,j-1]
# known displacements


P.shape


def rd(c1, c2):
    return np.sqrt((c1[0]-c2[0])**2+(c1[1]-c2[1])**2+(c1[2]-c2[2])**2)
#rbf as global support spline type
def rbf(r):
    return np.exp(-r**2)


# +
# # define lhs RBF kernel 
def M(coor,ctr_pts):
    M_m=np.zeros((N,N))
    for i in range(N):
        for j in range(N):
            M_m[i,j]=rbf(rd(ctr_pts[i],ctr_pts[j]))
    return M_m
# formulate the lhs part of the systems of linear equations to be dolved for coefficients
def lhs_(P,M_m):
    lhs=np.zeros((N+4,N+4))
    P_t=np.transpose(P)
    for i in range(N+4):
        for j in range(N+4):
            if (i<N and j<N): lhs[i,j]=M_m[i,j]
            elif (i<N and j>(N-1)): lhs[i, j]=P[i,j-N]
            elif (i>(N-1) and j<N):lhs[i,j]= P_t[i-N,j]
    return lhs
# Solve for the coefficents 
def sol_yb(lhs, g):
    rhs=np.zeros((N+4,3))
    for i in range(N):
        for j in range(3):    
            rhs[i,j]=g[i,j]
    
    inv_lhs=np.linalg.inv(lhs)
    yb= inv_lhs @ rhs
    gmma=np.zeros((N,3))
    beta=np.zeros((4,3))
    for i in range(N+4):
        if (i<N): gmma[i,:]=yb[i,:]
        else : beta[i-N,:]=yb[i,:]
    return (gmma, beta)
                

# -

# calculating the new coordinates from the previous ones
def s_x(ctr_pts, coor, beta, gmma):
    s_x=np.zeros(3)
    for i in range(3):
        for j in range(N):
            s_x[i]=s_x[i]+gmma[j,i]*rbf(rd(coor, ctr_pts[j]))
        s_x[i]=s_x[i]+beta[0,i]*1 + beta[1,i]*coor[0] + coor[1]*beta[2,i] + coor[2]*beta[3,i]
    return s_x


# ### Calling the required functions and solving in order for individual points 

M_m=M(coor,ctr_pts)
# M_m

lhs=lhs_(P,M_m)
# lhs

gmma, beta=sol_yb(lhs, g)

coor_new=np.zeros((num_coor,3))

cor_t=(1.0,1.0,1.0)
s_x(ctr_pts, cor_t, beta, gmma)

for i in range(num_coor):
    coor_new[i,:]=coor[i,:]+s_x(ctr_pts, coor[i], beta, gmma)

coor_new.shape

# plotting the deformed mesh in  
fig = pl.figure(1)
ax = fig.add_subplot(111, projection='3d')
ax.scatter(coor_new[:, 0], coor_new[:, 1], coor_new[:, 2], c='blue', marker='o')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
pl.show()

# Export to VTK
from pyevtk.hl import pointsToVTK
import numpy as np
# Example 1
npoints = 216
pressure = np.random.rand(npoints)
temp = np.linspace(1,npoints,npoints)
pointsToVTK("./rnd_points", np.array(coor_new[:, 0]), np.array(coor_new[:, 1]), np.array(coor_new[:, 2]), data={"temp": temp, "pressure": pressure})



StarcoderdataJupyterScriptsDedupFiltered