Homec4science
Diffusion MeshDRefinement (master) Browse

R12728/FEM.py
1e34ae800aa7master

Last Change
Raw File

FEM.py
View Options

import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.linalg import lu_factor, lu_solve
def solve_DD(coords, conn, areas, E, isDD, dataset, applyBC, max_iter = 1000, seed = 42, initial_guess=None, parallel=False):
"""
Coupled DD-FEM solver.
Input :
coords : [nxdim] float array of the points coordinates
conn : [ex2] integer array of connectivity
areas : [ex1] float array or bars cross-sections
E : stiffness of FE bars
isDD : [ex1] boolean array of whether the element is to be solved with DDCM
dataset : Dataset object
applyBC : functor who gives blocked DoF, force vector with BC included and displacement vector with BC included
max_iter : integer
seed : random generator seed
initial_guess : [ex2] or [eDDx2] array of strain-stress states for either each element or each DD element
"""
random.seed(seed)
isFEM = (1-isDD).astype('bool')
Ndd = sum(isDD)
Nfem = sum(isFEM)
Ndata = len(dataset.DP)
Nel = len(conn)
Ce = dataset.Ce
# FE stiffness matrix
K = compute_K(coords,conn[isFEM],areas,E)
# DD equivalent stiffness matrix
Keq = compute_K(coords,conn[isDD],areas,Ce)
# Matrix of shape function derivative for DD and element weights
B_DD, weights_DD = compute_B_and_weights(coords,conn,areas,isDD)
# Initial and boundary conditions
blocked_DOF, free_DOF, u, F = applyBC(coords, conn)
Nfree = np.sum(free_DOF)
F_free = F[free_DOF]
Keq_free = Keq[free_DOF,:][:,free_DOF]
K_free = K[free_DOF,:][:,free_DOF]
B_DD_free = B_DD[:,free_DOF]
# Assembly of coupling matrix
Coupled = compute_Coupled_Wuhan(K_free, Keq_free)
# Pre-factorization of coupling matrix
lu, piv = lu_factor(Coupled)
# List of guess index in dataset. Initialized to -1 as initial guess not from dataset
guess_list = np.array([-1 for _ in range(Ndd)]).astype('int')
if initial_guess is None:
guess = np.zeros((Ndd,2)) # Default initial guess
elif len(initial_guess) == Ndd:
guess = initial_guess
elif len(initial_guess) == Nel:
guess = initial_guess[isDD]
else:
raise Exception("Initial guess is in an unknown format")
converged = False
for i in range(max_iter):
rhs_top =( Ce*np.einsum('e,e,ei->i',weights_DD[isDD],guess[:,0],B_DD_free[isDD,:])
- (Keq[free_DOF,:][:,blocked_DOF]).dot(u[blocked_DOF]))
rhs_bottom = (F_free.ravel()
- np.einsum('e,e,ei->i',weights_DD[isDD],guess[:,1],B_DD_free[isDD,:])
- (K[free_DOF,:][:,blocked_DOF]).dot(u[blocked_DOF]))
rhs = np.vstack((rhs_top.reshape((-1,1)), rhs_bottom.reshape((-1,1)))).ravel()
u_eta = lu_solve((lu, piv), rhs)
u[free_DOF]=u_eta[:Nfree]
eta = u_eta[Nfree:]
epsi_comp = (B_DD).dot(u)[isDD]
sigma_comp = guess[:,1] + Ce*(B_DD_free).dot(eta)[isDD]
closest, distance = dataset.find_closest(epsi_comp, sigma_comp, weights_DD[isDD], parallel)
### ### ###
#print("Changed DP : ", sum([0 if a==b else 1 for a,b in zip(closest, guess_list)]))
if (closest == guess_list).all():
converged = True
break
guess_list=closest
guess=dataset.get_points(guess_list)
if converged:
print("Succesfull convergence in ", i, " steps")
else:
print("Not converged after ", i, " steps")
F = K.dot(u) + np.einsum('e,e,ei->i',weights_DD[isDD],sigma_comp,B_DD[isDD,:])
return u, F, np.vstack([epsi_comp, sigma_comp]).T, guess, eta
def solve_DD_drefine(coords, conn, areas, E, dataset, limit, applyBC, isDD=None, max_iter_DD = 1000, max_iter_r = 100):
"""
Solve the problem with iterative D-refinement. E is the linear stiffness modulus, limit is the "yield stress"
"""
# We need B to check the strains
B, _ = compute_B_and_weights(coords,conn,areas)
# Full FEM at first unless specified differently
if isDD is None:
isDD = np.zeros((len(conn),1)).ravel()
isDD = isDD.astype('bool')
# Solve a first time
u, F, comp, guess, eta = solve_DD(coords, conn, areas, E, isDD, dataset, applyBC, max_iter_DD)
converged = False
for i in range(max_iter_r):
#Check which elements are over the limit
above = np.abs(B.dot(u))>limit/E
new_isDD = isDD + above
if (new_isDD == isDD).all():
converged = True
print("Refinement complete after ", i, " steps.")
print(np.sum(isDD), "DD driven elements for a total of ", len(conn), " elements")
break
print(np.sum(np.logical_xor(new_isDD,isDD)), " element(s) have been refined")
isDD = new_isDD
u, F, comp, guess, eta = solve_DD(coords, conn, areas, E, isDD, dataset, applyBC, max_iter_DD)
if not converged:
print("Refinement complete after ", i, " steps.")
return u, F, comp, guess, eta, isDD
def NR_compute_K(coords, conn, areas, Elaw, epsi):
"""
Compute K given a strain-dependent stiffness for bar elements, in 2D or 3D
Ideally, this function should be merged with compute_K
"""
Nnodes, dim =coords.shape
Ktot = np.zeros((Nnodes*dim, Nnodes*dim))
for i,e in enumerate(conn):
Xa = coords[e[0]]
Xb = coords[e[1]]
dX = Xb-Xa
L =np.linalg.norm(dX)
if dim == 2:
C = dX[0]/L
S = dX[1]/L
K = areas[i]*Elaw(epsi[i])/L*np.array([[1,0,-1,0],
[0,0, 0,0],
[-1,0,1,0],
[0,0, 0,0]])
R = np.array([[C,S,0,0],
[-S,C,0,0],
[0,0,C,S],
[0,0,-S,C]])
K = R.T.dot(K).dot(R)
dof = np.array([e[0]*2,e[0]*2+1,e[1]*2,e[1]*2+1])
elif dim == 3:
Cx = dX[0]/L
Cy = dX[1]/L
Cz = dX[2]/L
inter = np.array([[Cx**2,Cx*Cy,Cx*Cz],
[Cx*Cy,Cy**2, Cy*Cz],
[Cx*Cz,Cy*Cz,Cz**2]])
K = areas[i]*Elaw(epsi[i]) / L * np.vstack([np.hstack([inter,-inter]), np.hstack([-inter,inter])])
dof = np.array([e[0]*3,e[0]*3+1,e[0]*3+2, e[1]*3,e[1]*3+1,e[1]*3+2])
else:
print("Unsupported dimension ", dim)
for m, d in enumerate(dof):
for n, p in enumerate(dof):
Ktot[d,p]+= K[m,n]
return Ktot
def NR(coords,conn, areas, sigma_law, Elaw, u, F, free_DOF,
tol=1e-10, max_iter = 100,
constantK=False, evolvingK=False, alpha=0.1):
"""
Solve FEM problem with Newton-Raphson. sigma_law and Elaw should be functor of strains
"""
B, weights = compute_B_and_weights(coords,conn,areas)
if constantK:
K = compute_K(coords, conn, areas, Elaw(0))
if evolvingK:
K0 = compute_K(coords, conn, areas, Elaw(0))
Nfree=np.sum(free_DOF)
for i in range(max_iter):
epsi = B.dot(u)
if constantK:
pass
elif evolvingK:
Kn = NR_compute_K(coords, conn, areas, Elaw, epsi)
K = K0*(1-alpha)+Kn*alpha
else:
K = NR_compute_K(coords, conn, areas, Elaw, epsi)
Fint = B.T.dot(np.multiply(weights,sigma_law(epsi)))
r = Fint-F
x = np.linspace(-0.1,0.1,10000)
error = np.linalg.norm(r[free_DOF])/Nfree
if error<tol:
print("Succesfull converge in", i, " steps")
return u, Fint
rhs = r[free_DOF]
du = np.linalg.solve(K[free_DOF,:][:,free_DOF], rhs)
u[free_DOF] -= du
print("Not converged in ", i, "steps\nError is ", np.linalg.norm(r))
return u, Fint
def compute_B_and_weights(coords,conn,areas,isDD=None):
"""
Compute B (shape derivative) matrix for bar elements, in 2D or 3D
Also return element weigths
"""
Nnodes, dim =coords.shape
Nel=len(conn)
B = np.zeros((Nel, Nnodes*dim))
weights = np.zeros((Nel,1))
if isDD is None:
isDD = np.ones_like(areas).astype('bool')
for i,e in enumerate(conn):
if not isDD[i]:
continue
Xa = coords[e[0]]
Xb = coords[e[1]]
dX = Xb-Xa
L =np.linalg.norm(dX)
if dX[0] != 0:
B[i,e[0]*dim]=-dX[0]/L**2
B[i,e[1]*dim]=dX[0]/L**2
if dX[1] != 0:
B[i,e[0]*dim+1]=-dX[1]/L**2
B[i,e[1]*dim+1]=dX[1]/L**2
if dim==3:
if dX[2] != 0:
B[i,e[0]*dim+2]=-dX[2]/L**2
B[i,e[1]*dim+2]= dX[2]/L**2
weights[i]=L*areas[i]
weights=weights.ravel()
return B, weights
def compute_M(coords, conn, areas, E, isDD=None):
"""
Compute mass matrix for bar elements, in 2D or 3D
"""
Nnodes, dim =coords.shape
M = np.zeros((Nnodes*dim, Nnodes*dim))
if isDD is None:
isDD = np.ones_like(areas).astype('bool')
for i,e in enumerate(conn):
if not isDD[i]:
continue
Xa = coords[e[0]]
Xb = coords[e[1]]
dX = Xb-Xa
L =np.linalg.norm(dX)
weights=L*areas[i]
M[e[0]*dim, e[0]*dim] += weights/2
M[e[0]*dim+1,e[0]*dim+1] += weights/2
M[e[1]*dim, e[1]*dim] += weights/2
M[e[1]*dim+1,e[1]*dim+1] += weights/2
if dim == 3:
M[e[0]*dim+2,e[0]*dim+2] += weights/2
M[e[1]*dim+2,e[1]*dim+2] += weights/2
for m, d in enumerate(dof):
for n, p in enumerate(dof):
Ktot[d,p]+= K[m,n]
return M
def compute_K(coords, conn, areas, E, isDD=None):
"""
Compute stiffness matrix for bar elements, in 2D or 3D
"""
Nnodes, dim =coords.shape
Ktot = np.zeros((Nnodes*dim, Nnodes*dim))
if isDD is None:
isDD = np.ones_like(areas).astype('bool')
for i,e in enumerate(conn):
if not isDD[i]:
continue
Xa = coords[e[0]]
Xb = coords[e[1]]
dX = Xb-Xa
L =np.linalg.norm(dX)
if dim == 2:
C = dX[0]/L
S = dX[1]/L
K = areas[i]*E/L*np.array([[1,0,-1,0],
[0,0, 0,0],
[-1,0,1,0],
[0,0, 0,0]])
R = np.array([[C,S,0,0],
[-S,C,0,0],
[0,0,C,S],
[0,0,-S,C]])
K = R.T.dot(K).dot(R)
dof = np.array([e[0]*2,e[0]*2+1,e[1]*2,e[1]*2+1])
elif dim == 3:
Cx = dX[0]/L
Cy = dX[1]/L
Cz = dX[2]/L
inter = np.array([[Cx**2,Cx*Cy,Cx*Cz],
[Cx*Cy,Cy**2, Cy*Cz],
[Cx*Cz,Cy*Cz,Cz**2]])
K = areas[i]*E / L * np.vstack([np.hstack([inter,-inter]), np.hstack([-inter,inter])])
dof = np.array([e[0]*3,e[0]*3+1,e[0]*3+2, e[1]*3,e[1]*3+1,e[1]*3+2])
else:
print("Unsupported dimension ", dim)
for m, d in enumerate(dof):
for n, p in enumerate(dof):
Ktot[d,p]+= K[m,n]
return Ktot
def compute_Coupled_Wuhan(K, Keq):
return np.hstack((np.vstack((Keq,K)),np.vstack((-K,Keq))))
def get_free(blocked_DOF):
return (1-blocked_DOF).astype('bool')
def solve_static(K,u,F,blocked_DOF):
"""
Solve the static FEM problem
"""
free_DOF = (1-blocked_DOF).astype('bool')
F_free = F[free_DOF]
K_free = K[free_DOF,:][:,free_DOF]
rhs = F_free - K[free_DOF,:][:,blocked_DOF].dot(u[blocked_DOF])
u[free_DOF]=np.linalg.solve(K_free,rhs)
return u
def plot_truss(coords, conn, disp=None, duplicate = False, isDD=None, nodes=False, aspect = None):
"""
Plot a truss in 2 or 3D.
Disp : displacements field, same format as coords
duplicate : plot both the deformed and undeformed case
isDD : boolean array, if True for a bar, it is plotted in blue
nodes : to plot the nodes as discs
aspect : aspect ration
"""
dim = coords.shape[1]
if aspect is None:
aspect = [1]*dim
base_string='-'
if nodes:
base_string='.-'
if disp is None:
disp=np.zeros_like(coords)
if isDD is None:
color = ['k'+base_string]*len(conn)
else:
color = ['r'+base_string if i else 'k' for i in isDD]
if dim == 2:
plot_truss_2d(coords, conn, disp, duplicate, color, aspect)
elif dim == 3:
plot_truss_3d(coords, conn, disp, duplicate, color, aspect)
else:
print("Dimension not supported : ", dim)
def plot_truss_3d(coords, conn, disp, duplicate = False, color=None, aspect = (1,1,1)):
ax = plt.figure().add_subplot(projection='3d')
ax.set_box_aspect(aspect)
for e,bar in enumerate(conn):
barcoords = coords[bar]
bardisp = disp[bar]
if duplicate:
plt.plot(barcoords[:,0],barcoords[:,1],barcoords[:,2],color[e], alpha=0.5)
plt.plot(barcoords[:,0]+bardisp[:,0],barcoords[:,1]+bardisp[:,1],barcoords[:,2]+bardisp[:,2],color[e])
def plot_truss_2d(coords, conn, disp, duplicate = False, color=None, aspect = (1,1)):
fig = plt.figure()
fax = fig.add_subplot(111, aspect='equal')
for e,bar in enumerate(conn):
barcoords = coords[bar]
bardisp = disp[bar]
if duplicate:
plt.plot(barcoords[:,0],barcoords[:,1],color[e], alpha=0.5)
plt.plot(barcoords[:,0]+bardisp[:,0],barcoords[:,1]+bardisp[:,1],color[e])
c4science · Help