sora/epfl-archive
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Disabled external gits

Browse Source
This commit is contained in:
choelzl
2022-04-07 18:46:57 +02:00
parent 88cb3426ad
commit 15e7120d6d
5316 changed files with 4563444 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

403
cs457-gc/assignment_2_3/src/Utils.py Normal file
View File

@ -0,0 +1,403 @@
import numpy as np
import time
import math
def conjugate_gradient(L_method, b):
'''
Finds an inexact Newton descent direction using Conjugate Gradient (CG)
Solves partially L(x) = b where A is positive definite, using CG.
The method should be take care of checking whether the added direction is
an ascent direction or not, and whether the residuals are small enough.
Details can be found in the handout.
Input:
- L_method : a method that computes the Hessian vector product. It should
take an array of shape (n,) and return an array of shape (n,)
- b : right hand side of the linear system (n,)
Output:
- p_star : np array of shape (n,) solving the linear system approximately
'''
N = b.shape[0]
fshape = (3*N,N)
z = np.zeros(fshape)
r = np.zeros(fshape)
d = np.zeros(fshape)
alpha = np.zeros(3*N)
beta = np.zeros(3*N)
r[0] = -b
d[0] = -r[0]
ck = 0
for k in range(0,3*N - 1):
Ldk = L_method(d[k])
if (d[k].T @ Ldk) <= 0:
if (k==0):
return b
else:
return z[k]
alpha[k] = (r[k].T @ r[k])/(d[k] @ Ldk)
z[k+1] = z[k] + alpha[k] * d[k]
r[k+1] = r[k] + alpha[k] * Ldk
va = np.linalg.norm(r[k+1])
vb = math.sqrt(np.linalg.norm(b))
if va < min(0.5,vb) * np.linalg.norm(b):
break
beta[k+1] = (r[k].T @ r[k+1])/(r[k].T @ r[k])
d[k+1] = -r[k+1] + beta[k+1]*d[k]
ck = k
return z[ck+1]
def compute_inverse_approximate_hessian_matrix(sk, yk, invB_prev):
'''
Input:
- sk : previous step x_{k+1} - x_k, shape (n, 1)
- yk : grad(f)_{k+1} - grad(f)_{k}, shape (n, 1)
- invB_prev : previous Hessian estimate Bk, shape (n, n)
Output:
- invB_new : previous Hessian estimate Bk, shape (n, n)
'''
invB_new = invB_prev.copy()
invB_new += (sk.T @ yk + yk.T @ invB_prev @ yk) / ((sk.T @ yk) ** 2) * (sk @ sk.T)
prod = (invB_prev @ yk) @ sk.T
invB_new -= (prod + prod.T) / (sk.T @ yk)
return invB_new
def equilibrium_convergence_report_GD(solid, v_init, n_steps, step_size, thresh=1e-3):
'''
Finds the equilibrium by minimizing the total energy using gradient descent.
Input:
- solid : an elastic solid to optimize
- v_init : the initial guess for the equilibrium position
- n_step : number of optimization steps
- step_size : scaling factor of the gradient when taking the step
- thresh : threshold to stop the optimization process on the gradient's magnitude
Ouput:
- report : a dictionary containing various quantities of interest
'''
solid.update_def_shape(v_init)
energies_el = np.zeros(shape=(n_steps+1,))
energies_ext = np.zeros(shape=(n_steps+1,))
residuals = np.zeros(shape=(n_steps+1,))
times = np.zeros(shape=(n_steps+1,))
step_sizes = np.zeros(shape=(n_steps,))
energies_el[0] = solid.energy_el
energies_ext[0] = solid.energy_ext
residuals[0] = np.linalg.norm((solid.f + solid.f_ext)[solid.free_idx, :])
idx_stop = n_steps
energy_tot_prev = energies_el[0] + energies_ext[0]
t_start = time.time()
for i in range(n_steps):
## TODO: Find the descent direction
## descent_dir has shape (#v, 3)
descent_dir = (solid.f + solid.f_ext)
step_size_tmp = step_size
max_l_iter = 20
for l_iter in range(max_l_iter):
step_size_tmp *= 0.5
solid.displace(step_size_tmp * descent_dir)
## TODO: Check if the armijo rule is satisfied
## energy_tot_tmp is the current total energy
## armijo is a boolean that says whether the condition is satisfied
energy_tot_tmp = solid.energy_ext + solid.energy_el
grad_tmp = -descent_dir
armijo = energy_tot_tmp < energy_tot_prev + thresh * step_size_tmp * ((descent_dir*grad_tmp).sum())
if armijo or l_iter == max_l_iter-1:
break
else:
solid.displace(-step_size_tmp * descent_dir)
step_sizes[i] = step_size_tmp
# Measure the force residuals
energies_el[i+1] = solid.energy_el
energies_ext[i+1] = solid.energy_ext
residuals[i+1] = np.linalg.norm((solid.f + solid.f_ext)[solid.free_idx, :])
energy_tot_prev = energy_tot_tmp
if residuals[i+1] < thresh:
energies_el[i+1:] = energies_el[i+1]
energies_ext[i+1:] = energies_ext[i+1]
residuals[i+1:] = residuals[i+1]
idx_stop = i
break
times[i+1] = time.time() - t_start
report = {}
report['energies_el'] = energies_el
report['energies_ext'] = energies_ext
report['residuals'] = residuals
report['times'] = times
report['idx_stop'] = idx_stop
report['step_sizes'] = step_sizes
return report
def equilibrium_convergence_report_BFGS(solid, v_init, n_steps, step_size, thresh=1e-3):
'''
Finds the equilibrium by minimizing the total energy using BFGS.
Input:
- solid : an elastic solid to optimize
- v_init : the initial guess for the equilibrium position
- n_step : number of optimization steps
- step_size : scaling factor of the direction when taking the step
- thresh : threshold to stop the optimization process on the gradient's magnitude
Ouput:
- report : a dictionary containing various quantities of interest
'''
solid.update_def_shape(v_init)
energies_el = np.zeros(shape=(n_steps+1,))
energies_ext = np.zeros(shape=(n_steps+1,))
residuals = np.zeros(shape=(n_steps+1,))
times = np.zeros(shape=(n_steps+1,))
energies_el[0] = solid.energy_el
energies_ext[0] = solid.energy_ext
## TODO: Collect free vertex positions
## grad_tmp is the current flattened gradient of the total energy
## with respect to the free vertices
grad_tmp = - (solid.f + solid.f_ext)[solid.free_idx,:].reshape(-1,1)
residuals[0] = np.linalg.norm(grad_tmp)
idx_stop = n_steps
energy_tot_prev = energies_el[0] + energies_ext[0]
## TODO: Collect free vertex positions
## v_tmp are the current flattened free vertices
v_tmp = solid.v_def[solid.free_idx,:].reshape(-1,1)
dir_zeros = np.zeros_like(solid.v_def)
invB_prev = np.eye(v_tmp.shape[0])
t_start = time.time()
for i in range(n_steps):
dir_tmp = - invB_prev @ grad_tmp
dir_zeros[solid.free_idx, :] = dir_tmp.reshape(-1, 3)
step_size_tmp = step_size
max_l_iter = 20
for l_iter in range(max_l_iter):
step_size_tmp *= 0.5
solid.displace(step_size_tmp * dir_zeros)
## TODO: Check if the armijo rule is satisfied
## energy_tot_tmp is the current total energy
## armijo is a boolean that says whether the condition is satisfied
energy_tot_tmp = solid.energy_el + solid.energy_ext
armijo = energy_tot_tmp < energy_tot_prev + thresh * step_size_tmp*((dir_tmp*grad_tmp).sum())
if armijo or l_iter == max_l_iter-1:
break
else:
solid.displace(-step_size_tmp * dir_zeros)
## TODO: Update all quantities
## v_new are the new flattened free vertices
## grad_new is the new flattened gradient of the total energy
## with respect to the free vertices
v_new = solid.v_def[solid.free_idx,:].reshape(-1,1)
grad_new = - (solid.f+solid.f_ext)[solid.free_idx,:].reshape(-1,1)
invB_prev = compute_inverse_approximate_hessian_matrix(v_new - v_tmp,
grad_new - grad_tmp,
invB_prev)
v_tmp = v_new.copy()
grad_tmp = grad_new.copy()
energies_el[i+1] = solid.energy_el
energies_ext[i+1] = solid.energy_ext
residuals[i+1] = np.linalg.norm(grad_tmp)
if residuals[i+1] < thresh:
residuals[i+1:] = residuals[i+1]
energies_el[i+1:] = energies_el[i+1]
energies_ext[i+1:] = energies_ext[i+1]
idx_stop = i
break
times[i+1] = time.time() - t_start
report = {}
report['energies_el'] = energies_el
report['energies_ext'] = energies_ext
report['residuals'] = residuals
report['times'] = times
report['idx_stop'] = idx_stop
return report
def equilibrium_convergence_report_NCG(solid, v_init, n_steps, thresh=1e-3):
'''
Finds the equilibrium by minimizing the total energy using Newton CG.
Input:
- solid : an elastic solid to optimize
- v_init : the initial guess for the equilibrium position
- n_step : number of optimization steps
- thresh : threshold to stop the optimization process on the gradient's magnitude
Ouput:
- report : a dictionary containing various quantities of interest
'''
solid.update_def_shape(v_init)
energies_el = np.zeros(shape=(n_steps+1,))
energies_ext = np.zeros(shape=(n_steps+1,))
residuals = np.zeros(shape=(n_steps+1,))
times = np.zeros(shape=(n_steps+1,))
energies_el[0] = solid.energy_el
energies_ext[0] = solid.energy_ext
residuals[0] = np.linalg.norm((solid.f + solid.f_ext)[solid.free_idx, :])
idx_stop = n_steps
t_start = time.time()
for i in range(n_steps):
# Take a Newton step
solid.equilibrium_step()
# Measure the force residuals
energies_el[i+1] = solid.energy_el
energies_ext[i+1] = solid.energy_ext
residuals[i+1] = np.linalg.norm((solid.f + solid.f_ext)[solid.free_idx, :])
if residuals[i+1] < thresh:
residuals[i+1:] = residuals[i+1]
energies_el[i+1:] = energies_el[i+1]
energies_ext[i+1:] = energies_ext[i+1]
idx_stop = i
break
times[i+1] = time.time() - t_start
report = {}
report['energies_el'] = energies_el
report['energies_ext'] = energies_ext
report['residuals'] = residuals
report['times'] = times
report['idx_stop'] = idx_stop
return report
import matplotlib.pyplot as plt
def fd_validation_ext(solid):
epsilons = np.logspace(-9,-3,100)
perturb_global = np.random.uniform(-1e-3, 1e-3, size=solid.v_def.shape)
solid.displace(perturb_global)
v_def = solid.v_def.copy()
perturb = np.random.uniform(-1, 1, size=solid.v_def.shape)
errors = []
for eps in epsilons:
# Back to original
solid.update_def_shape(v_def)
grad = np.zeros(solid.f_ext.shape)
grad[solid.free_idx] = -solid.f_ext.copy()[solid.free_idx]
an_delta_E = (grad*perturb).sum()
# One step forward
solid.displace(perturb * eps)
E1 = solid.energy_ext.copy()
# Two steps backward
solid.displace(-2*perturb * eps)
E2 = solid.energy_ext.copy()
# Compute error
fd_delta_E = (E1 - E2)/(2*eps)
errors.append(abs(fd_delta_E - an_delta_E)/abs(an_delta_E))
plt.loglog(epsilons, errors)
plt.grid()
plt.show()
def fd_validation_elastic(solid):
epsilons = np.logspace(-9,-3,100)
perturb_global = np.random.uniform(-1e-3, 1e-3, size=solid.v_def.shape)
solid.displace(perturb_global)
v_def = solid.v_def.copy()
solid.make_elastic_forces()
perturb = np.random.uniform(-1, 1, size=solid.v_def.shape)
errors = []
for eps in epsilons:
# Back to original
solid.update_def_shape(v_def)
solid.make_elastic_forces()
grad = np.zeros(solid.f.shape)
grad[solid.free_idx] = -solid.f.copy()[solid.free_idx]
an_delta_E = (grad*perturb).sum()
# One step forward
solid.displace(perturb * eps)
E1 = solid.energy_el.copy()
# Two steps backward
solid.displace(-2*perturb * eps)
E2 = solid.energy_el.copy()
# Compute error
fd_delta_E = (E1 - E2)/(2*eps)
errors.append(abs(fd_delta_E - an_delta_E)/abs(an_delta_E))
solid.displace(perturb * eps)
plt.loglog(epsilons, errors)
plt.grid()
plt.show()
def fd_validation_elastic_differentials(solid):
epsilons = np.logspace(-9, 3,500)
perturb_global = 1e-3*np.random.uniform(-1., 1., size=solid.v_def.shape)
solid.displace(perturb_global)
v_def = solid.v_def.copy()
perturb = np.random.uniform(-1, 1, size=solid.v_def.shape)
errors = []
for eps in epsilons:
# Back to original
solid.update_def_shape(v_def)
perturb_0s = np.zeros_like(perturb)
perturb_0s[solid.free_idx] = perturb[solid.free_idx]
an_df = solid.compute_force_differentials(perturb_0s)[solid.free_idx, :]
an_df_full = np.zeros(solid.f.shape)
an_df_full[solid.free_idx] = an_df.copy()
# One step forward
solid.displace(perturb * eps)
f1 = solid.f[solid.free_idx, :]
f1_full = np.zeros(solid.f.shape)
f1_full[solid.free_idx] = f1
# Two steps backward
solid.displace(-2*perturb * eps)
f2 = solid.f[solid.free_idx, :]
f2_full = np.zeros(solid.f.shape)
f2_full[solid.free_idx] = f2
# Compute error
fd_delta_f = (f1_full - f2_full)/(2*eps)
norm_an_df = np.linalg.norm(an_df_full)
norm_error = np.linalg.norm(an_df_full - fd_delta_f)
errors.append(norm_error/norm_an_df)
plt.loglog(epsilons, errors)
plt.grid()
plt.show()

163
cs457-gc/assignment_2_3/src/elasticenergy.py Normal file
View File

@ -0,0 +1,163 @@
import numpy as np
from numpy.core.einsumfunc import einsum
from numpy.core.fromnumeric import swapaxes
class ElasticEnergy:
def __init__(self, young, poisson):
'''
Input:
- young : Young's modulus [Pa]
- poisson : Poisson ratio
'''
self.young = young
self.poisson = poisson
self.lbda = young * poisson / ((1 + poisson) * (1 - 2 * poisson))
self.mu = young / (2 * (1 + poisson))
self.psi = None
self.E = None
self.P = None
self.dE = None
self.dP = None
def make_energy_density(self, jac):
'''
This method computes the energy density at each tetrahedron (#t,),
and stores the result in self.psi
Input:
- jac : jacobian of the deformation (#t, 3, 3)
Updated attributes:
- psi : energy density per tet (#t,)
'''
print("Please specify the kind of elasticity model.")
raise NotImplementedError
def make_strain_tensor(self, jac):
'''
This method computes the strain tensor (#t, 3, 3), and stores it in self.E
Input:
- jac : jacobian of the deformation (#t, 3, 3)
Updated attributes:
- E : strain induced by the deformation (#t, 3, 3)
'''
print("Please specify the kind of elasticity model.")
raise NotImplementedError
def make_piola_kirchhoff_stress_tensor(self, jac):
'''
This method computes the stress tensor (#t, 3, 3), and stores it in self.P
Input:
- jac : jacobian of the deformation (#t, 3, 3)
Updated attributes:
- P : stress tensor induced by the deformation (#t, 3, 3)
'''
print("Please specify the kind of elasticity model.")
raise NotImplementedError
def make_differential_strain_tensor(self, jac, dJac):
'''
This method computes the differential of strain tensor (#t, 3, 3),
and stores it in self.dE
Input:
- jac : jacobian of the deformation (#t, 3, 3)
- dJac : differential of the jacobian of the deformation (#t, 3, 3)
Updated attributes:
- dE : differential of the strain tensor (#t, 3, 3)
'''
print("Please specify the kind of elasticity model.")
raise NotImplementedError
def make_differential_piola_kirchhoff_stress_tensor(self, jac, dJac):
'''
This method computes the differential of the stress tensor (#t, 3, 3), and stores it in self.dP
Input:
- jac : jacobian of the deformation (#t, 3, 3)
- dJac : differential of the jacobian of the deformation (#t, 3, 3)
Updated attributes:
- dP : differential of the stress tensor (#t, 3, 3)
'''
print("Please specify the kind of elasticity model.")
raise NotImplementedError
class LinearElasticEnergy(ElasticEnergy):
def __init__(self, young, poisson):
super().__init__(young, poisson)
def make_energy_density(self, jac):
fa = self.mu * np.einsum("ijk,ijk->i",self.E, self.E)
fb = self.lbda/2.0 * np.einsum('ijj->i' ,self.E)**2
self.psi = fa + fb
def make_strain_tensor(self, jac):
self.E = (jac+np.transpose(jac,axes=(0,2,1)))/2.0 - np.identity(3)
def make_piola_kirchhoff_stress_tensor(self, jac):
fa = 2.0*self.mu*self.E
fb = self.lbda* np.einsum('i,jk->ijk',np.einsum('ijj->i',self.E), np.identity(3))
self.P = fa + fb
def make_differential_strain_tensor(self, jac, dJac):
self.dE = (dJac+np.transpose(dJac,axes=(0,2,1)))/2.0
def make_differential_piola_kirchhoff_stress_tensor(self, jac, dJac):
fa = 2.0*self.mu*self.dE
fb = self.lbda* np.einsum("i,jk->ijk",np.trace(self.dE, axis1=1,axis2=2), np.identity(3))
self.dP = fa + fb
class NeoHookeanElasticEnergy(ElasticEnergy):
def __init__(self, young, poisson):
super().__init__(young, poisson)
self.logJ = None
self.Finv = None
def make_energy_density(self, jac):
i1 = np.einsum('ijk,ijk->i',jac,jac)
j = np.log(np.linalg.det(jac))
fa = self.mu/2.0 * (i1-3 - 2*j)
fb = self.lbda/2.0 * j**2
self.psi = fa + fb
def make_strain_tensor(self, jac):
pass
def make_piola_kirchhoff_stress_tensor(self, jac):
'''
Additional updated attributes:
- logJ ; log of the determinant of the jacobians (#t,)
- Finv : inverse of the jacobians (#t, 3, 3)
'''
self.logJ = np.log(np.linalg.det(jac))
self.Finv = np.linalg.inv(jac)
FinvT = np.transpose(self.Finv,axes=(0,2,1))
fa = self.mu * (jac - FinvT)
fb = self.lbda * np.einsum("i,ijk->ijk",self.logJ, FinvT)
self.P = fa +fb
def make_differential_strain_tensor(self, jac, dJac):
pass
def make_differential_piola_kirchhoff_stress_tensor(self, jac, dJac):
self.logJ = np.log(np.linalg.det(jac))
self.Finv = np.linalg.inv(jac)
FinvT = np.transpose(self.Finv,axes=(0,2,1))
dFT = np.transpose(dJac,axes=(0,2,1))
fa = self.mu * dJac
fb = np.einsum("i,ijk->ijk",(self.mu - self.lbda * self.logJ),FinvT @ dFT @ FinvT)
fc = self.lbda * np.einsum("i,ijk->ijk",np.trace(self.Finv @ dJac, axis1=1,axis2=2),FinvT)
self.dP = fa + fb + fc

404
cs457-gc/assignment_2_3/src/elasticsolid.py Normal file
View File

@ -0,0 +1,404 @@
import numpy as np
from numpy import linalg
from scipy import sparse
from Utils import *
import igl
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# ELASTIC SOLID CLASS
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
class ElasticSolid(object):
def __init__(self, v_rest, t, ee, rho=1, pin_idx=[], f_mass=None):
'''
Input:
- v_rest : position of the vertices of the mesh (#v, 3)
- t : indices of the element's vertices (#t, 4)
- ee : elastic energy object that can be found in elasticenergy.py
- rho : mass per unit volume [kg.m-3]
- pin_idx : list or np array of vertex indices to pin
- f_mass : external force per unit mass (3,) [N.kg-1]
'''
self.v_rest = v_rest.copy()
self.v_def = v_rest.copy()
self.t = t
self.ee = ee
self.rho = rho
self.pin_idx = np.array(pin_idx)
self.f_mass = f_mass.copy()
self.free_idx = None
self.pin_mask = None
self.W0 = None
self.Dm = None
self.Bm = None
self.rest_barycenters = None
self.W = None
self.Ds = None
self.F = None
self.def_barycenters = None
self.energy_el = None
self.energy_ext = None
self.f = None
self.f_vol = None
self.f_ext = None
self.make_free_indices_and_pin_mask()
self.update_rest_shape(self.v_rest)
self.update_def_shape(self.v_def)
## Utils ##
def vertex_tet_sum(self, data):
'''
Distributes data specified at each tetrahedron to the neighboring vertices.
All neighboring vertices will receive the value indicated at the corresponding tet position in data.
Input:
- data : np array of shape (#t,) or (4*#t,)
Output:
- data_sum : np array of shape (#v,), containing the summed data
'''
i = self.t.flatten('F') # (4*#t,)
j = np.arange(len(self.t)) # (#t,)
j = np.tile(j, 4) # (4*#t,)
if len(data) == len(self.t):
data = data[j]
# Has shape (#v, #t)
m = sparse.coo_matrix((data, (i, j)), (len(self.v_rest), len(self.t)))
return np.array(m.sum(axis=1)).flatten()
## Precomputation ##
def make_free_indices_and_pin_mask(self):
'''
Should list all the free indices and the pin mask.
Updated attributes:
- free_index : np array of shape (#free_vertices,) containing the list of unpinned vertices
- pin_mask : np array of shape (#v, 1) containing 1 at free vertex indices and 0 at pinned vertex indices
'''
self.pin_mask = np.ones((self.v_rest.shape[0],1))
self.free_idx = np.arange(0,self.v_rest.shape[0])
if(self.pin_idx.shape[0]<=0):
return
self.pin_mask[self.pin_idx] = 0
self.free_idx = np.where(self.pin_mask==1)[0]
## Methods related to rest quantities ##
def make_rest_barycenters(self):
'''
Construct the barycenters of the undeformed configuration
Updated attributes:
- rest_barycenters : np array of shape (#t, 3) containing the position of each tet's barycenter
'''
self.rest_barycenters = np.einsum("ijk->ik",self.v_rest[self.t])/4
def make_rest_shape_matrices(self):
'''
Construct Ds that has shape (#t, 3, 3), and its inverse Bm
Updated attributes:
- Dm : np array of shape (#t, 3, 3) containing the shape matrix of each tet
- Bm : np array of shape (#t, 3, 3) containing the inverse shape matrix of each tet
'''
tv = self.v_rest[self.t]
tva = tv[:,:3,:]
tvb = tv[:,3,:].reshape((self.t.shape[0],1,3))
tvb= np.repeat(tvb, 3, axis=1)
tvs = tva - tvb
tvsf = np.transpose(tvs,(0,2,1))
self.Dm = tvsf
self.Bm = np.linalg.inv(tvsf)
def update_rest_shape(self, v_rest):
'''
Updates the vertex position, the shape matrices Dm and Bm, the volumes W0,
and the mass matrix at rest
Input:
- v_rest : position of the vertices of the mesh at rest state (#v, 3)
Updated attributes:
- v_rest : np array of shape (#v, 3) containing the position of each vertex at rest
- W0 : np array of shape (#t,) containing the signed volume of each tet
'''
self.v_rest = v_rest
self.make_rest_barycenters()
self.make_rest_shape_matrices()
self.W0 = -np.linalg.det(self.Dm)/6
self.update_def_shape(self.v_def)
self.make_volumetric_and_external_forces()
## Methods related to deformed quantities ##
def make_def_barycenters(self):
'''
Construct the barycenters of the deformed configuration
Updated attributes:
- def_barycenters : np array of shape (#t, 3) containing the position of each tet's barycenter
'''
self.def_barycenters = np.einsum("ijk->ik",self.v_def[self.t])/4
def make_def_shape_matrices(self):
'''
Construct Ds that has shape (#t, 3, 3)
Updated attributes:
- Ds : np array of shape (#t, 3, 3) containing the shape matrix of each tet
'''
tv = self.v_def[self.t]
tva = tv[:,:3,:]
tvb = tv[:,3,:].reshape((self.t.shape[0],1,3))
tvb= np.repeat(tvb, 3, axis=1)
tvs = tva - tvb
tvsf = np.transpose(tvs,(0,2,1))
self.Ds = tvsf
def make_jacobians(self):
'''
Compute the current Jacobian of the deformation
Updated attributes:
- F : np array of shape (#t, 3, 3) containing Jacobian of the deformation in each tet
'''
self.F = self.Ds @ self.Bm # <=> np.einsum("lij,ljk->lik",self.Ds,self.Bm)
def update_def_shape(self, v_def):
'''
Updates the vertex position, the Jacobian of the deformation, and the
resulting elastic forces.
Input:
- v_def : position of the vertices of the mesh (#v, 3)
Updated attributes:
- v_def : np array of shape (#v, 3) containing the position of each vertex after deforming the solid
- W : np array of shape (#t,) containing the signed volume of each tet
'''
self.v_def = self.v_rest
self.v_def[self.free_idx] = v_def[self.free_idx]
self.make_def_barycenters()
self.make_def_shape_matrices()
self.W = -np.linalg.det(self.Ds)/6
self.make_jacobians()
if(self.f_vol is None):
self.make_volumetric_and_external_forces()
self.make_elastic_forces()
self.make_elastic_energy()
self.make_external_energy()
def displace(self, v_disp):
'''
Displace the whole mesh so that v_def += v_disp
Input:
- v_disp : displacement of the vertices of the mesh (#v, 3)
'''
self.update_def_shape(self.v_def + v_disp)
## Energies ##
def make_elastic_energy(self):
'''
This updates the elastic energy
Updated attributes:
- energy_el : elastic energy of the system [J]
'''
self.ee.make_strain_tensor(self.F)
self.ee.make_energy_density(self.F)
self.energy_el = np.sum(self.W0 * self.ee.psi)
def make_external_energy(self):
'''
This computes the external energy potential
Updated attributes:
- energy_ext : postential energy due to external forces [J]
'''
dx = self.def_barycenters - self.rest_barycenters
ifv = np.einsum('i,ij -> ij', self.W0, self.f_vol)
self.energy_ext = - np.einsum("i,ij,ij->",self.W0, self.f_vol, dx)
## Forces ##
def make_elastic_forces(self):
'''
This method updates the elastic forces stored in self.f (#v, 3)
Updated attributes:
- f : elastic forces per vertex (#v, 3)
- ee : elastic energy, some attributes should be updated
'''
self.ee.make_strain_tensor(self.F)
self.ee.make_piola_kirchhoff_stress_tensor(self.F)
tdm = np.transpose(self.Bm,(0,2,1))
ptdm = self.ee.P @ tdm
f123 = -np.einsum("i,ijk->ikj",self.W0,ptdm )
f4 = -f123[:,0] - f123[:,1] - f123[:,2]
H = np.hstack((f123,f4[:,None]))
x = self.vertex_tet_sum(H[:,:,0].flatten('F'))
y = self.vertex_tet_sum(H[:,:,1].flatten('F'))
z = self.vertex_tet_sum(H[:,:,2].flatten('F'))
self.f = np.hstack((x[:,None],y[:,None],z[:,None]))
def make_volumetric_and_external_forces(self):
'''
Convert force per unit mass to volumetric forces, then distribute
the forces to the vertices of the mesh.
Updated attributes:
- f_vol : np array of shape (#t, 3) external force per unit volume acting on the tets
- f_ext : np array of shape (#v, 3) external force acting on the vertices
'''
mpv = self.f_mass*self.rho
self.f_vol = np.ones((self.W0.shape[0],1)) * mpv
ifv = np.einsum('i,ij -> ij', self.W0, self.f_vol)
f_ext_x = self.vertex_tet_sum(ifv[:,0])
f_ext_y = self.vertex_tet_sum(ifv[:,1])
f_ext_z = self.vertex_tet_sum(ifv[:,2])
self.f_ext = np.stack((f_ext_x,f_ext_y,f_ext_z),axis=1) / 4.0
## Force Differentials
def compute_force_differentials(self, v_disp):
'''
This computes the differential of the force given a displacement dx,
where df = df/dx|x . dx = - K(x).dx. Where K(x) is the stiffness matrix (or Hessian)
of the solid. Note that the implementation doesn't need to construct the stiffness matrix explicitly.
Input:
- v_disp : displacement of the vertices of the mesh (#v, 3)
Output:
- df : force differentials at the vertices of the mesh (#v, 3)
Updated attributes:
- ee : elastic energy, some attributes should be updated
'''
# First compute the differential of the Jacobian
tv = v_disp[self.t]
tva = tv[:,:3,:]
tvb = tv[:,3,:].reshape((-1,1,3))
tvb= np.repeat(tvb, 3, axis=1)
tvs = tva - tvb
tvsf = np.transpose(tvs,(0,2,1))
dDs = tvsf
dF = dDs @ self.Bm
# Then update differential quantities in self.ee
self.ee.make_differential_strain_tensor(self.F,dF)
self.ee.make_differential_piola_kirchhoff_stress_tensor(self.F,dF)
# Compute the differential of the forces
tdm = np.transpose(self.Bm,(0,2,1))
ptdm = self.ee.dP @ tdm
f123 = -np.einsum("i,ijk->ikj",self.W0,ptdm)
f4 = -f123[:,0] - f123[:,1] - f123[:,2]
H = np.hstack((f123,f4[:,None]))
x = self.vertex_tet_sum(H[:,:,0].flatten('F'))
y = self.vertex_tet_sum(H[:,:,1].flatten('F'))
z = self.vertex_tet_sum(H[:,:,2].flatten('F'))
return np.hstack((x[:,None],y[:,None],z[:,None]))
def equilibrium_step(self, step_size_init = 2, max_l_iter = 20, c1 = 1e-4):
'''
This function displaces the whole solid to the next deformed configuration
using a Newton-CG step.
Updated attributes:
- LHS : The hessian vector product
- RHS : Right hand side for the conjugate gradient linear solve
Other than them, only attributes updated by displace(self, v_disp) should be changed
'''
# Define LHS
def LHS(dx):
'''
Should implement the Hessian-Vector Product L(dx), and take care of pinning constraints
as described in the handout.
'''
if len(dx.shape)==1:
dx = dx.reshape((-1,3))
if (dx.shape[0] < self.v_rest.shape[0]):
dxf = np.zeros((self.v_rest.shape[0],3))
dxf[self.free_idx] = dx
dx = dxf
return -self.compute_force_differentials(dx)[self.free_idx,:].flatten()
self.LHS = LHS # Save to class for testing
# Define RHS
RHS = (self.f + self.f_ext)[self.free_idx].flatten()
self.RHS = RHS # Save to class for testing
# Use conjugate gradient to find a descent direction
# (see conjugate_gradient in Utils.py)
dx_CG = conjugate_gradient(LHS,RHS).reshape(-1,3) # shape (#v, 3)
if (dx_CG.shape[0] < self.v_rest.shape[0]):
dxf = np.zeros((self.v_rest.shape[0],3))
dxf[self.free_idx] = dx_CG
dx_CG = dxf
# Run line search on the direction
step_size = step_size_init
ft = self.f + self.f_ext
energy_tot_prev = (self.energy_el + self.energy_ext)
for l_iter in range(max_l_iter):
step_size *= 0.5 # DO NOT CHANGE
dx_search = dx_CG * step_size
self.displace(dx_search)
# Current force norm for unpinned vertices
f_tmp = (self.f + self.f_ext)[self.free_idx]
g = np.linalg.norm(f_tmp)
energy_tot_tmp = self.energy_el + self.energy_ext
#Make sure you only use the gradient of the unpinned vertices
# USE parameter c1 = 1e-4 (from the function arguments) when you submit
armijo = energy_tot_tmp < energy_tot_prev + c1 * step_size*(-dx_CG*dx_CG).sum()
if armijo or l_iter == max_l_iter-1:
print("Energy: " + str(energy_tot_tmp) + " Force norm: " + str(g) + " Line search Iters: " + str(l_iter))
break
else:
self.displace(-dx_search)

13
cs457-gc/assignment_2_3/src/objectives.py Normal file
View File

@ -0,0 +1,13 @@
import numpy as np
def objective_target_BV(v, vt, bv):
'''
Input:
- v : array of shape (#v, 3), containing the current vertices position
- vt : array of shape (#bv, 3), containing the target surface
- bv : boundary vertices index (#bv,)
Output:
- objective : single scalar measuring the deviation from the target shape
'''
return 0.

51
cs457-gc/assignment_2_3/src/stretch_utils.py Normal file
View File

@ -0,0 +1,51 @@
import numpy as np
from objectives import *
from Utils import *
def stretch_from_point(v, stretches, v_idx):
'''
Input:
- v : vertices to be deformed (#v, 3)
- stretches : np array of shape (3,) giving the stretching coefficients along x, y, z
- v_idx : vertex index around which the mesh is stretched
Output:
- v_stretched : deformed vertices (#v, 3)
'''
v_stretched = v
return v_stretched
def report_stretches(solid, v_rest_init, bv, v_target, stretches, lowest_pin_idx):
'''
Input:
- solid : an elastic solid to deform
- v_rest_init : reference vertex position to compress/stretch (#v, 3)
- bv : boundary vertices index (#bv,)
- v_target : target boundary vertices position
- stretches : np array of shape (n_stretches,) containing a stretch factors to try out
- lowest_pin_idx : index of the pinned index that has the lowest z coordinate
Output:
- list_v_rest : list of n_stretches rest vertex positions that have been tried out
- list_v_eq : list of the corresponding n_stretches equilibrium vertex positions
- target_closeness : np array of shape (n_stretches,) containing the objective values for each stretch factor
'''
list_v_rest = []
list_v_eq = []
target_closeness = np.zeros(shape=(stretches.shape[0],))
v_rest_tmp = v_rest_init.copy()
for i, stretch in enumerate(stretches):
# Compute and update the rest shape of the solid
# Compute new equilibrium (use the previous equilibrium state as initial guess if available)
# You may use equilibrium_convergence_report_NCG to find the equilibrium (20 steps and thresh=1 should do)
# Update guess for next stretch factor
# Fill in the
list_v_rest.append(solid.v_rest.copy())
list_v_eq.append(solid.v_def.copy())
target_closeness[i] = objective_target_BV(solid.v_def, v_target, bv)
return list_v_rest, list_v_eq, target_closeness