Page MenuHomec4science
Files F85658681

test_hertz.py
No OneTemporary
Actions

File Metadata

Created
Mon, Sep 30, 20:42

test_hertz.py
View Options

from rough_contact import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import math
#generate surface
n = 32
L = 1.
a = SurfaceReal(n,L)
r = 10.
applied_load = .001
N2 = a.size()*a.size()
ss = a.getThis()
print type(ss)
for i in range(0,n):
for j in range(0,n):
x = 1.*i - n/2;
y = 1.*j - n/2;
d = (x*x + y*y)*1./n/n*L*L;
if d < r*r: ss[i,j] = - r + math.sqrt(r*r-d);
else : ss[i,j] = - r;
#plot surface function
fig_by_titles = {}
axe_by_titles = {}
img_by_titles = {}
cbar_by_titles = {}
def plotCurve(curve,title):
if not title in fig_by_titles:
fig = plt.figure()
axe = fig.add_subplot(1,1,1)
fig_by_titles[title] = fig
axe_by_titles[title] = axe
fig = fig_by_titles[title]
axe = axe_by_titles[title]
axe.clear()
axe.set_title(title)
axe.plot(curve[:,0],np.log(curve[:,1]),'-')
# print curve
fig.show()
plt.draw()
def plotSurface(surf,title):
try: s = surf.getWrappedNumpy()
except: s = surf
print "{0} range: {1} - {2}".format(title,s.min(),s.max())
if len(s.shape) == 1:
n = np.sqrt(s.shape[0])
s = np.reshape(s,(n,n))
if not title in fig_by_titles:
fig = plt.figure()
axe = fig.add_subplot(1,1,1)
img = axe.imshow(s.real)
cbar = plt.colorbar(img)
fig_by_titles[title] = fig
axe_by_titles[title] = axe
img_by_titles[title] = img
cbar_by_titles[title] = cbar
fig = fig_by_titles[title]
axe = axe_by_titles[title]
img = img_by_titles[title]
img.set_data(s.real)
img.autoscale()
axe.set_title(title)
fig.show()
# plt.show()
plt.draw()
#plot the surface
plotSurface(a,'surface')
#construct the bem
bem = BemFFTBase(a)
#make a bem calculation
bem.setEffectiveModulus(1.)
bem.computeDisplacement()
disp = bem.getDisplacements()
################################################################
class ConvergenceMonitor():
def __init__(self):
self.cpt = 1
self.curves = {}
def addValue(self,title,value):
if title not in self.curves:
self.curves[title] = np.zeros((0,2))
a = np.array([[self.cpt,value]])
self.curves[title] = np.concatenate((self.curves[title],a))
def incCpt(self):
self.cpt += 1
def plot(self):
for title in self.curves.keys():
plotCurve(self.curves[title],title)
convergence = ConvergenceMonitor()
################################################################
#compute the objective function
def display(x):
F = computeF(x,bem,terms,display=True)
convergence.addValue('F',F)
_kwargs = computeDisplacement(x,bem)
for f in ['disp','trac','gap']:
plotSurface(_kwargs[f],f)
grad = computeGradientF(x,bem,terms)
convergence.addValue('gradientNorm',np.linalg.norm(grad))
plotSurface(grad[:n*n],'gradient')
convergence.incCpt()
convergence.plot()
pass
################################################################
class EnergyTerm:
def __init__(self,need_multiplier=False):
self.need_multiplier = need_multiplier
self.name = self.__class__.__name__
def setMultiplierRank(self,m_rank):
self.m_rank = m_rank
################################################################
class DeformationEnergy(EnergyTerm):
def __init__(self):
EnergyTerm.__init__(self)
def compute(self,trac=None,disp=None,**kwargs):
return 0.5*(trac*disp).sum()
def computeGradient(self,disp=None,**kwargs):
return 0.5*disp
################################################################
class OrthogonalityCondition(EnergyTerm):
def __init__(self):
EnergyTerm.__init__(self)
def compute(self,trac=None,gap=None,**kwargs):
return (trac*gap).sum()
def computeGradient(self,gap=None,**kwargs):
return gap
################################################################
class LoadBalance(EnergyTerm):
def __init__(self,applied_pressure,beta):
EnergyTerm.__init__(self)
self.applied_pressure = applied_pressure
self.beta = beta
def compute(self,trac=None,**kwargs):
avg_pressure = np.average(trac)
return self.beta*(avg_pressure-self.applied_pressure)**2
def computeGradient(self,trac=None,multipliers=None,N=None,**kwargs):
avg_pressure = np.average(trac)
grad = self.beta*(avg_pressure-self.applied_pressure)*1./N*np.ones(trac.shape)
return grad
################################################################
class LoadBalanceMultiplier(EnergyTerm):
def __init__(self,applied_pressure):
EnergyTerm.__init__(self,need_multiplier=True)
self.applied_pressure = applied_pressure
def compute(self,trac=None,multipliers=None,**kwargs):
mult = multipliers[self.m_rank]
avg_pressure = np.average(trac)
return mult*(avg_pressure-self.applied_pressure)**2
def computeGradient(self,trac=None,multipliers=None,N=None,**kwargs):
mult = multipliers[self.m_rank]
avg_pressure = np.average(trac)
grad1 = mult*(avg_pressure-self.applied_pressure)*1./N*np.ones(trac.shape)
grad2 = (avg_pressure-self.applied_pressure)**2
return grad1,grad2
################################################################
class LoadBalanceQuadraticMultiplier(EnergyTerm):
def __init__(self,applied_pressure):
EnergyTerm.__init__(self,need_multiplier=True)
self.applied_pressure = applied_pressure
def compute(self,trac=None,multipliers=None,**kwargs):
mult = multipliers[self.m_rank]
avg_pressure = np.average(trac)
return mult**2*(avg_pressure-self.applied_pressure)**2
def computeGradient(self,trac=None,multipliers=None,N=None,**kwargs):
mult = multipliers[self.m_rank]
avg_pressure = np.average(trac)
grad1 = mult**2*(avg_pressure-self.applied_pressure)*1./N*np.ones(trac.shape)
grad2 = 2*mult*(avg_pressure-self.applied_pressure)**2
return grad1,grad2
################################################################
class PositivePenalty(EnergyTerm):
def __init__(self,beta):
EnergyTerm.__init__(self)
self.beta = beta
def compute(self,gradient_flag=False,beta=None,trac=None,gap=None,**kwargs):
penalty1 = np.copy(trac)
penalty1[penalty1 > 0] = 0
penalty2 = np.copy(gap)
penalty2[penalty2 > 0] = 0
return self.beta*((penalty1**2).sum()+(penalty2**2).sum())
def computeGradient(self,gradient_flag=False,beta=None,trac=None,gap=None,**kwargs):
penalty1 = np.copy(trac)
penalty1[penalty1 > 0] = 0
penalty2 = np.copy(gap)
penalty2[penalty2 > 0] = 0
return 2.*(penalty1+penalty2)*self.beta
################################################################
def computeDisplacement(unknowns,bemModel):
surface = bemModel.getSurface().getWrappedNumpy().real
N = surface.shape[0]*surface.shape[1]
trac = unknowns[0:N]
trac = np.reshape(trac,surface.shape)
multipliers = unknowns[N:]
np.copyto(bemModel.getTractions().getWrappedNumpy(),trac)
bemModel.computeDisplacement()
disp = bemModel.getDisplacements().getWrappedNumpy().real
bemModel.computeTrueDisplacement()
true_disp = bemModel.getTrueDisplacements().getWrappedNumpy().real
bemModel.computeGap()
gap = bemModel.getGap().getWrappedNumpy().real
# plotSurface(trac,'trac')
# print disp
# print true_disp
# print surface
# print gap
# plotSurface(disp,'disp')
# plotSurface(true_disp,'truedisp')
# plotSurface(surface,'surface')
# plotSurface(gap,'gap')
# print trac
return {'trac':trac,'multipliers':multipliers,'disp':true_disp,'surface':surface,'gap':gap,'N':N}
################################################################
def computeGradF(unknowns,bemModel,terms):
_kwargs = computeDisplacement(unknowns,bemModel)
avg_pressure = np.average(trac)
# grad1 = disp - surface + 2.*multipliers*np.ones(surface.shape)*(sum_pressure-applied_pressure)
N = surface.shape[0]*surface.shape[1]
grad1 = disp - surface + multipliers**2*np.ones(surface.shape)*(avg_pressure-applied_pressure)*1./N
# grad1 += positivePenaltyGradient(trac,1.)
grad1 = np.reshape(grad1,N)
grad = np.empty(N+1)
grad[:-1] = grad1
grad[-1] = 2.*multipliers*(avg_pressure-applied_pressure)**2
# grad[-1] = 0
# print "current unknowns(grad)"
# print unknowns
# print "current gradient"
# print grad
return grad
################################################################
def computeGradientF(unknowns,bemModel,terms):
_kwargs = computeDisplacement(unknowns,bemModel)
grad = np.zeros(unknowns.shape)
N = _kwargs['N']
for t in terms:
if t.need_multiplier:
g,m_var = t.computeGradient(**_kwargs)
else: g = t.computeGradient(**_kwargs)
# plotSurface(g,'gradient-{0}'.format(t.name))
if not len(g.shape) == 1:
g = np.reshape(g,g.shape[0]*g.shape[1])
grad[:N] += g[:N]
if t.need_multiplier:
grad[N+t.m_rank] = m_var
# plotSurface(grad[:N],'gradient-total')
return grad
################################################################
def computeF(unknowns,bemModel,terms,display=False):
_kwargs = computeDisplacement(unknowns,bemModel)
res = [t.compute(**_kwargs) for t in terms]
F = sum(res)
if display:
avg_pressure = np.average(_kwargs['trac'])
formated = ["{0:.5e}".format(t) for t in [avg_pressure]+res+[F]]
print "\t".join(formated)
return F
################################################################
def printHead(terms):
titles = ['avg_pressure']
titles += [i.name for i in terms]
titles += 'F'
formated = ["{0}".format(t) for t in titles]
print "\t".join(formated)
################################################################
terms = [DeformationEnergy(),OrthogonalityCondition(),LoadBalance(applied_load,10000.),PositivePenalty(10000.)]
#terms = [LoadBalance(applied_load)]
multipliers_number = 0
for t in terms:
if t.need_multiplier:
t.setMultiplierRank(multipliers_number)
multipliers_number += 1
x0 = np.zeros((n*n+multipliers_number,))
#x0[:n*n] = applied_load
x0[n*n:] = 1.
printHead(terms)
x = scipy.optimize.fmin_cg(computeF,x0,args=(bem,terms),callback=display,gtol=1e-3,fprime=computeGradientF)
x = scipy.optimize.fmin_cg(computeF,x,args=(bem,terms),callback=display,gtol=1e-8)#,fprime=computeGradientF)
plt.show()
xx = np.reshape(x[:n*n],a.getWrappedNumpy().shape)
print xx
plotSurface(xx,'pressure solution')
plt.show()
F = computeF(x,bem,terms)
bem.computeDisplacement()
plotSurface(bem.getTrueDisplacements(),'displacements')
plt.show()

Event Timeline

c4science · Help