Page MenuHomec4science
Files F90471819

test_hertz_spectral.py
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 2, 00:02

test_hertz_spectral.py
View Options

from python_surface import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import scipy
import math
import sys
#generate surface
n = 32
L = 1.
a = Surface(n,L)
r = 2.
applied_load = .001
N2 = a.size()*a.size()
ss = a.getWrappedNumpy()
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)
fig.canvas.set_window_title(title)
axe.plot(curve[:,0],curve[:,1],'-')
# print curve
fig.show()
plt.draw()
def plotSurface(surf,title,complex_part='real'):
try: s = surf.getWrappedNumpy()
except: s = surf
if complex_part == 'real': s = s.real
else: s = s.imag
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)
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)
img.autoscale()
axe.set_title(title)
fig.canvas.set_window_title(title)
fig.show()
# plt.show()
plt.draw()
#plot the surface
plotSurface(a,'surface')
################################################################
class ConvergenceMonitor():
def __init__(self):
self.cpt = 1
self.curves = {}
def addValue(self,title,value,log_flag=False):
if title not in self.curves:
self.curves[title] = np.zeros((0,2))
if log_flag: a = np.array([[self.cpt,np.log10(np.abs(value))]])
else: 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 = computeInternals(x,bem)
for f in ['disp','trac']:
plotSurface(_kwargs[f],f)
grad = computeGradientF(x,bem,terms)
convergence.addValue('gradientNorm',np.linalg.norm(grad),log_flag=True)
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_spectral=None,N=None,**kwargs):
grad = 1./N*disp_spectral.conj()
grad[0,0] = 0.
return grad
################################################################
class OrthogonalityCondition(EnergyTerm):
def __init__(self):
EnergyTerm.__init__(self)
def compute(self,trac=None,disp=None,surface=None,**kwargs):
return (trac*(disp-surface)).sum()
def computeGradient(self,disp_spectral=None,surface_spectral=None,N=None,**kwargs):
grad = 1./N*(2.*disp_spectral - surface_spectral).conj()
grad[0,0] = 0
# plotSurface(grad,'orthogonality gradient')
return grad
################################################################
class OrthogonalityCondition2(EnergyTerm):
def __init__(self):
EnergyTerm.__init__(self)
def compute(self,trac=None,disp=None,surface=None,**kwargs):
return (trac*(disp-surface)).sum()
def computeGradient(self,disp_spectral=None,surface_spectral=None,N=None,**kwargs):
grad = 1./N*(disp_spectral - surface_spectral).conj()
grad[0,0] = 0
# plotSurface(grad,'orthogonality gradient')
return grad
################################################################
class NegativePenalty(EnergyTerm):
def __init__(self,beta):
EnergyTerm.__init__(self)
self.beta = beta
def compute(self,trac=None,disp=None,surface=None,**kwargs):
trac = np.copy(trac)
trac[trac>0] = 0
plotSurface(trac,'negative traction')
return self.beta*(trac**2).sum()
def computeGradient(self,trac=None,N=None,**kwargs):
trac = np.copy(trac)
trac[trac>0] = 0
trac = 2.*trac
trac = np.fft.fft2(trac)
grad = self.beta/N*trac.conj()
grad[0,0] = 0.
plotSurface(grad,'penalty gradient')
return grad
################################################################
def computeInternals(unknowns,bemModel):
surface = bemModel.getSurface().getWrappedNumpy().real
N = surface.shape[0]*surface.shape[1]
spectral_trac = unknowns[:N]
spectral_trac = np.reshape(spectral_trac,surface.shape)
np.copyto(bemModel.getSpectralTractions().getWrappedNumpy(),spectral_trac)
bemModel.computeSpectralDisplacements()
bemModel.computeRealDisplacementsFromSpectral()
bemModel.computeRealTractionsFromSpectral()
disp = bemModel.getDisplacements().getWrappedNumpy().real
sdisp = bemModel.getSpectralDisplacements().getWrappedNumpy()
trac = bemModel.getTractions().getWrappedNumpy().real
ssurf = bemModel.getSpectralSurface().getWrappedNumpy()
return {'trac':trac,'disp':disp,'surface':surface,'N':N,'disp_spectral':sdisp,'surface_spectral':ssurf}
################################################################
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 = computeInternals(unknowns,bemModel)
grad = np.zeros(unknowns.shape,dtype=complex)
N = _kwargs['N']
for t in terms:
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]
# plotSurface(grad[:N],'gradient-total')
return grad
################################################################
def computeF(unknowns,bemModel,terms,display=False):
_kwargs = computeInternals(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)
################################################################
def stepFunction(trac, Delta):
trac = np.copy(trac.getWrappedNumpy())
trac[trac>0] = 0
res = 0.5*(trac**2)/(Delta**2)
return res
def stepFunctionVariation(trac, Delta):
trac = trac.getWrappedNumpy()
trac[trac>0] = 0
res = 2.*trac/(Delta**2)
return res
################################################################
##make a bem calculation (old method first)
#old_bem = BemFFT(a)
#F = old_bem.computeEquilibrium(1e-9,applied_load,0);
#
#trac = old_bem.getTractions()
#plotSurface(trac,'initial solution trac')
#plotSurface(old_bem.getDisplacements(),'initial solution displacement')
#sdisp = old_bem.getSpectralDisplacements().getWrappedNumpy()
#plotSurface(sdisp,'initial solution spectral displacement',complex_part='imag')
#strac = old_bem.getSpectralTractions().getWrappedNumpy()
#plotSurface(strac,'initial solution spectral tractions',complex_part='imag')
#plotSurface(stepFunction(trac,1./n),'initial solution step function')
#plotSurface(stepFunctionVariation(trac,1./n),'initial solution step function variation')
#
##old_bem.computeSpectralDisplacement
#plt.show()
#sys.exit(0)
#make a bem calculation
bem = BemFFTBase(a)
bem.setEffectiveModulus(1.)
x = np.zeros((n*n,),dtype=complex)
x[:n*n] = 0
x[0] = applied_load*n*n
for beta in [1e4, 2e4, 3e4]:
terms = [DeformationEnergy(),OrthogonalityCondition(),NegativePenalty(beta)]
#terms = [LoadBalance(applied_load)]
printHead(terms)
x = scipy.optimize.fmin_cg(computeF,x,args=(bem,terms),callback=display,gtol=1e-4,fprime=computeGradientF)
xx = np.reshape(x[:n*n],a.getWrappedNumpy().shape)
print xx
plotSurface(xx,'spectral pressure solution')
F = computeF(x,bem,terms)
trac = bem.getTractions()
plotSurface(trac,'pressure solution')
plotSurface(bem.getDisplacements(),'displacements')
plt.show()

Event Timeline

c4science · Help