#!/usr/bin/env python
from nutils import *
import scipy.sparse.linalg
import time
def precon(M,tol_orth=.99,tol_rem=1e-12):
islands = []
p = numpy.ones(M.shape[0],dtype=bool)
D = scipy.sparse.lil_matrix(scipy.sparse.diags(1./numpy.sqrt(M.diagonal()),0))
i_v,j_v,dummy = scipy.sparse.find( numpy.abs( scipy.sparse.tril( D.dot(M.dot(D.T)) - scipy.sparse.diags( D.dot(M.dot(D.T)).diagonal(),0 ) ) ) > tol_orth )
while len(i_v): # while the matrix requires orthonormalisation
## sort indices
for index in zip( i_v, j_v ):
index = set( index )
for island in list(islands):
if island & index:
islands.remove( island )
index |= island
islands.append( index )
## orthonormalise
for island in islands:
island_unsort = list( island )
other_indices = list( set(numpy.arange(M.shape[0])) - island )
intersections = numpy.ravel((numpy.abs(D.dot(M.dot(D.T))[:,other_indices][island_unsort,:])>tol_rem).sum(1))
island = []
for island_index in intersections.argsort():
island.append( island_unsort[island_index] )
M_p = M[:,island][island,:]
D_p = numpy.diagflat(1./numpy.sqrt(M_p.diagonal()),0)
for ind_1 in numpy.arange( len( island ) ):
for ind_2 in numpy.arange( ind_1 ):
D_p[ind_1,:] -= D_p[ind_2,:]*D_p[ind_2,:].dot(M_p.dot(D_p[ind_1,:]))
if D_p[ind_1,:].dot(M_p.dot(D_p[ind_1,:])) <= tol_rem:
D_p[ind_1,:] = 0
p[island[ind_1]] = False
else:
D_p[ind_1,:] = D_p[ind_1,:]/numpy.sqrt(D_p[ind_1,:].dot(M_p.dot(D_p[ind_1,:])))
D[:,island][island,:] = D_p
i_v,j_v,dummy = scipy.sparse.find( numpy.abs( scipy.sparse.tril( D.dot(M.dot(D.T)) - scipy.sparse.diags( D.dot(M.dot(D.T)).diagonal(),0 ) ) ) > tol_orth )
return D[p,:]
def plotmesh ( domain, geom, disp=0., norm_stress=0., mag_stress=0., name='curved' ):
points, disp, norm_stress, mag_stress = domain.boundary.simplex.elem_eval( [geom,disp,norm_stress,mag_stress], ischeme='vtk', separate=True )
with plot.VTKFile(name+'_boundary') as vtkfile:
vtkfile.unstructuredgrid( points )
vtkfile.pointdataarray( 'disp', disp )
vtkfile.pointdataarray( 'norm_stress', norm_stress )
vtkfile.pointdataarray( 'mag_stress', mag_stress )
# points, displacement, variable = domain .simplex.elem_eval( [geom,disp,var], ischeme='vtk', separate=True )
# with plot.VTKFile(name+'_volume') as vtkfile:
# vtkfile.unstructuredgrid( points )
# vtkfile.pointdataarray( 'disp', disp )
# vtkfile.pointdataarray( 'norm_stress', norm_stress )
# vtkfile.pointdataarray( 'mag_stress', mag_stress )
return
def main( R = 2. , #Radius curve
r = 1. , #Radius rod
l = 1. , #Length rod
lmbda = 1. , #First Lame parameter
mu = 1. , #Second Lame parameter
local = False , #Stabilisation
mmr = 0 , #Minmaxrefine
degree = 1 , #Degree
solvetol = 1e-10 , #Solvetol
h_v = numpy.array([2.,1.,.5]) #Gridsize
):
stress = library.Hooke(lmbda=lmbda,mu=mu)
energy_norm_v = numpy.zeros(h_v.shape)
energy_tang_v = numpy.zeros(h_v.shape)
for h_n in numpy.arange(h_v.shape[0]):
h = h_v[h_n]
log.user('starting with h_n =',h_n)
t_old = time.time()
ni = numpy.ceil( l /h) #Number of inner elements
no = numpy.ceil((R+r)/h) #Number of outer elements
nz = numpy.ceil( r /h) #Number of z elements
domain, geom = mesh.rectilinear([ numpy.linspace( -h*(ni+degree), h*(no+degree), ni+no+2*degree+1 ),
numpy.linspace( -h*(ni+degree), h*(no+degree), ni+no+2*degree+1 ),
numpy.linspace( -h*(nz+degree), h*(nz+degree), 2*nz+2*degree+1 ) ])
x, y, z = geom
bc_norm = function.asarray([ 1., 0., 0. ])
bc_tang = function.asarray([ 0., -z , y-R ])
dxy = function.piecewise( function.arctan2(y,x), [0,numpy.pi/2], (x-R)**2 ,
x**2 + y**2 + R**2 - 2*R*function.sqrt(x**2+y**2),
(y-R)**2 )
lvl_rad = r**2 - dxy - z**2
lvl_bottom = l + y
lvl_left = l + x
domain = domain.trim( lvl_rad , maxrefine = mmr+h_v.shape[0]-h_n-1, name = 'trim_rad' )
domain = domain.trim( lvl_bottom, maxrefine = mmr+h_v.shape[0]-h_n-1, name = 'trim_bottom' )
domain = domain.trim( lvl_left , maxrefine = mmr+h_v.shape[0]-h_n-1, name = 'trim_left' )
tang_tensor = function.eye( 3 ) - geom.normal()[:,_]*geom.normal()[_,:]
funcsp = domain.basis( 'spline', degree=degree ).vector( 3 )
symgrad_funcsp = funcsp.symgrad(geom)
stresssp = stress(symgrad_funcsp)
log.user('created geometry in',time.time()-t_old)
##############################
# Volume integration #
##############################
t_old = time.time()
elasticity_matrix = function.outer( symgrad_funcsp, stresssp ).sum([-1,-2])
log.user('created volumetric operators in',time.time()-t_old)
t_old = time.time()
volume_matrix, locs, vols = domain.integrate( [ elasticity_matrix,
(symgrad_funcsp**2).sum([-1,-2])[:,_]*geom[_,:],
(symgrad_funcsp**2).sum([-1,-2])
], geometry=geom, ischeme='gauss'+str(2*degree) )
volume_matrix = volume_matrix.toscipy()
locs = locs.toarray()
xlocs = locs[2*funcsp.shape[0]/3:: ,0]/vols[2*funcsp.shape[0]/3:: ]
ylocs = locs[ : funcsp.shape[0]/3,1]/vols[ : funcsp.shape[0]/3]
zlocs = locs[ funcsp.shape[0]/3:2*funcsp.shape[0]/3,2]/vols[ funcsp.shape[0]/3:2*funcsp.shape[0]/3]
xfacts = (ylocs-numpy.average(ylocs)) * volume_matrix.diagonal()[ : funcsp.shape[0]/3]
yfacts = (zlocs-numpy.average(zlocs)) * volume_matrix.diagonal()[ funcsp.shape[0]/3:2*funcsp.shape[0]/3]
zfacts = (xlocs-numpy.average(xlocs)) * volume_matrix.diagonal()[2*funcsp.shape[0]/3:: ]
x_index = numpy.argmax(volume_matrix.diagonal()[ : funcsp.shape[0]/3])
y_index = funcsp.shape[0]/3 + numpy.argmax(volume_matrix.diagonal()[ funcsp.shape[0]/3 : 2*funcsp.shape[0]/3])
z_index = 2*funcsp.shape[0]/3 + numpy.argmax(volume_matrix.diagonal()[2*funcsp.shape[0]/3 : : ])
x_index1 = numpy.argmin(xfacts); x_index2 = numpy.argmax(xfacts)
y_index1 = funcsp.shape[0]/3 + numpy.argmin(yfacts); y_index2 = funcsp.shape[0]/3 + numpy.argmax(yfacts)
z_index1 = 2*funcsp.shape[0]/3 + numpy.argmin(zfacts); z_index2 = 2*funcsp.shape[0]/3 + numpy.argmax(zfacts)
log.user('finished volume integration in',time.time()-t_old)
##############################
# Setting stability constants#
##############################
if local:
if h_n == 0:
#Create local bases indices
x_power = numpy.reshape(numpy.tile(numpy.arange(degree+1),(degree+1,degree+1,1)).transpose((0,1,2)),-1)
y_power = numpy.reshape(numpy.tile(numpy.arange(degree+1),(degree+1,degree+1,1)).transpose((1,2,0)),-1)
z_power = numpy.reshape(numpy.tile(numpy.arange(degree+1),(degree+1,degree+1,1)).transpose((2,0,1)),-1)
basis_indices_mu = []; basis_indices_lmbda = []
for i in numpy.arange((degree+1)**3):
basis_indices_mu .append((x_power[i],y_power[i],z_power[i]))
basis_indices_lmbda.append((x_power[i],y_power[i],z_power[i]))
basis_indices_mu .remove(( 0 , 0 , 0 ))
basis_indices_mu .remove(( 1 , 0 , 0 ))
basis_indices_mu .remove(( 0 , 1 , 0 ))
basis_indices_mu .remove(( 0 , 0 , 1 ))
basis_indices_lmbda.remove(( degree, degree, degree ))
#Initialise
beta_mu_norm_dict = {}; beta_mu_tang_dict = {}; beta_lmbda_dict = {}
#Index Dirichlet boundary elements
elem_belems = {}
for belem in domain.boundary['trim_bottom,trim_left']:
etrans = belem.transform.lookup( domain.edict )
ielem = domain.edict[ etrans ]
elem_belems.setdefault(ielem,[]).append( belem )
#Loop over Dirichlet boundary elements
for ielem, belems in log.iter( 'boundary element', elem_belems.items() ):
velem = domain.elements[ielem]
#Create local topology
vdomain = topology.Topology( [velem] )
bdomain = topology.Topology( belems )
#Calculate centre of mass
volume, geom_sep = vdomain.integrate( [ 1., geom ], geometry=geom, ischeme='gauss1' )
geom_local = geom - geom_sep / volume
#Create local bases
basis_blocks = [ f * function.eye(3) for f in function.product( function.power( geom_local[_,:], basis_indices_mu ), axis=-1 ) ]
basis_blocks.append( function.diagonalize(geom_local) )
basis_blocks.append([ [ geom_local[1], geom_local[0], 0 ],
[ 0 , geom_local[2], geom_local[1] ],
[ geom_local[2], 0 , geom_local[0] ] ])
basis_stab_mu = function.concatenate( basis_blocks, axis=0 )
basis_stab_lmbda = function.product( function.power( geom_local[_,:], basis_indices_lmbda ), axis=-1 )
#Operators
symgrad_basis_mu = basis_stab_mu.symgrad(geom)
E_mu_norm_operator = function.outer( symgrad_basis_mu.dotnorm(geom).dotnorm(geom) )
E_mu_tang_operator = function.outer( (symgrad_basis_mu.dotnorm(geom)[:,_,:] * tang_tensor[_,:,:]).sum(-1) ).sum(-1)
E_lmbda_operator = function.outer( basis_stab_lmbda )
A_mu_operator = function.outer( symgrad_basis_mu ).sum([-1,-2])
A_lmbda_operator = E_lmbda_operator
#Integrate
E_mu_norm_elem, E_mu_tang_elem, E_lmbda_elem = bdomain.integrate( [ E_mu_norm_operator, E_mu_tang_operator, E_lmbda_operator ], geometry=geom, ischeme='gauss'+str(2*degree) )
A_mu_elem, A_lmbda_elem = vdomain.integrate( [ A_mu_operator, A_lmbda_operator ], geometry=geom, ischeme='gauss'+str(2*degree) )
E_mu_norm_elem = E_mu_norm_elem.toscipy(); E_mu_tang_elem = E_mu_tang_elem.toscipy(); E_lmbda_elem = E_lmbda_elem.toscipy()
A_mu_elem = A_mu_elem.toscipy(); A_lmbda_elem = A_lmbda_elem.toscipy()
#Eigenvalue problem
# D = precon(A_mu_elem ); D = D.toarray()
# beta_mu_norm_dict[velem.transform] = 4*mu *max(scipy.linalg.eigh(D.dot(E_mu_norm_elem.dot(D.T)),D.dot(A_mu_elem .dot(D.T)),eigvals_only=True))
# beta_mu_tang_dict[velem.transform] = 4*mu *max(scipy.linalg.eigh(D.dot(E_mu_tang_elem.dot(D.T)),D.dot(A_mu_elem .dot(D.T)),eigvals_only=True))
# D = precon(A_lmbda_elem); D = D.toarray()
# beta_lmbda_dict [velem.transform] = 2*lmbda*max(scipy.linalg.eigh(D.dot(E_lmbda_elem .dot(D.T)),D.dot(A_lmbda_elem.dot(D.T)),eigvals_only=True))
D = precon(A_mu_elem ); D = D.toarray()
beta_h_mu_norm = 4*mu *max(scipy.linalg.eigh(D.dot(E_mu_norm_elem.dot(D.T)),D.dot(A_mu_elem .dot(D.T)),eigvals_only=True))*h
beta_h_mu_tang = 4*mu *max(scipy.linalg.eigh(D.dot(E_mu_tang_elem.dot(D.T)),D.dot(A_mu_elem .dot(D.T)),eigvals_only=True))*h
D = precon(A_lmbda_elem); D = D.toarray()
beta_h_lmbda = 2*lmbda*max(scipy.linalg.eigh(D.dot(E_lmbda_elem .dot(D.T)),D.dot(A_lmbda_elem.dot(D.T)),eigvals_only=True))*h
# D = scipy.sparse.csr_matrix(precon(A_mu_elem ))
# (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_norm_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem .dot(D.T)), which='LM' )
# beta_mu_norm_dict[velem.transform] = 4*mu * eigenvalue
# (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_tang_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem .dot(D.T)), which='LM' )
# beta_mu_tang_dict[velem.transform] = 4*mu * eigenvalue
# D = scipy.sparse.csr_matrix(precon(A_lmbda_elem))
# (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_lmbda_elem .dot(D.T)), k=1, M=D.dot(A_lmbda_elem.dot(D.T)), which='LM' )
# beta_lmbda_dict [velem.transform] = 2*lmbda * eigenvalue
# D = scipy.sparse.csr_matrix(precon(A_mu_elem ))
# (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_norm_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem .dot(D.T)), which='LM' )
# beta_h_mu_norm = 4*mu * eigenvalue * h
# (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_tang_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem .dot(D.T)), which='LM' )
# beta_h_mu_tang = 4*mu * eigenvalue * h
# D = scipy.sparse.csr_matrix(precon(A_lmbda_elem))
# (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_lmbda_elem .dot(D.T)), k=1, M=D.dot(A_lmbda_elem.dot(D.T)), which='LM' )
# beta_h_lmbda = 2*lmbda * eigenvalue * h
break
# #Assign
# beta_mu_norm = function.Elemwise( beta_mu_norm_dict, () )
# beta_mu_tang = function.Elemwise( beta_mu_tang_dict, () )
# beta_lmbda = function.Elemwise( beta_lmbda_dict , () )
beta_mu_norm = beta_h_mu_norm/h
beta_mu_tang = beta_h_mu_tang/h
beta_lmbda = beta_h_lmbda /h
##############################
# Boundary integration #
##############################
t_old = time.time()
#Operators
flux_norm_matrix = -function.outer( stresssp.dotnorm(geom).dotnorm(geom) , funcsp.dotnorm(geom) )
flux_norm_matrix = flux_norm_matrix + flux_norm_matrix.T
flux_norm_vector = -stresssp.dotnorm(geom).dotnorm(geom) * bc_norm.dotnorm(geom)
penalty_norm_matrix = (beta_mu_norm+beta_lmbda) * function.outer( funcsp.dotnorm(geom) )
penalty_norm_vector = (beta_mu_norm+beta_lmbda) * funcsp.dotnorm(geom) * bc_norm.dotnorm(geom)
flux_tang_matrix = -function.outer( tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_], funcsp[:,_,:] ).sum([-1,-2])
flux_tang_matrix = flux_tang_matrix + flux_tang_matrix.T
flux_tang_vector = -(tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
penalty_tang_matrix = beta_mu_tang * function.outer( ( tang_tensor[_,:,:] * funcsp[:,_,:] ).sum(-1) ).sum(-1)
penalty_tang_vector = beta_mu_tang * ( tang_tensor[_,:,:] * funcsp[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
log.user('created boundary operators in',time.time()-t_old)
t_old = time.time()
#Integrate
boundary_norm_matrix_bottom, boundary_tang_matrix_bottom = domain.boundary['trim_bottom'].integrate( [ flux_norm_matrix + penalty_norm_matrix,
flux_tang_matrix + penalty_tang_matrix,
], geometry=geom, ischeme='gauss'+str(2*degree) )
log.user('finished bottom integration in',time.time()-t_old)
t_old = time.time()
boundary_norm_matrix_left, boundary_norm_vector, boundary_tang_matrix_left, boundary_tang_vector = domain.boundary['trim_left' ].integrate( [ flux_norm_matrix + penalty_norm_matrix,
flux_norm_vector + penalty_norm_vector,
flux_tang_matrix + penalty_tang_matrix,
flux_tang_vector + penalty_tang_vector
], geometry=geom, ischeme='gauss'+str(2*degree) )
log.user('finished left integration in',time.time()-t_old)
boundary_norm_matrix = boundary_norm_matrix_bottom.toscipy() + boundary_norm_matrix_left.toscipy() + boundary_tang_matrix_bottom.toscipy()
boundary_tang_matrix = boundary_tang_matrix_bottom.toscipy() + boundary_tang_matrix_left.toscipy()
else:
t_old = time.time()
#Operators
flux_norm_matrix = -function.outer( stresssp.dotnorm(geom).dotnorm(geom) , funcsp.dotnorm(geom) )
flux_norm_matrix = flux_norm_matrix + flux_norm_matrix.T
flux_norm_vector = -stresssp.dotnorm(geom).dotnorm(geom) * bc_norm.dotnorm(geom)
penalty_norm_matrix = function.outer( funcsp.dotnorm(geom) )
penalty_norm_vector = funcsp.dotnorm(geom) * bc_norm.dotnorm(geom)
stability_norm_matrix = function.outer( stresssp.dotnorm(geom).dotnorm(geom) )
flux_tang_matrix = -function.outer( tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_], funcsp[:,_,:] ).sum([-1,-2])
flux_tang_matrix = flux_tang_matrix + flux_tang_matrix.T
flux_tang_vector = -(tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
penalty_tang_matrix = function.outer( ( tang_tensor[_,:,:] * funcsp[:,_,:] ).sum(-1) ).sum(-1)
penalty_tang_vector = ( tang_tensor[_,:,:] * funcsp[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
stability_tang_matrix = function.outer( ( tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,_,:] ).sum(-1) ).sum(-1)
log.user('created boundary operators in',time.time()-t_old)
t_old = time.time()
#Integrate
flux_norm_matrix_bottom, \
penalty_norm_matrix_bottom, \
stability_norm_matrix_bottom, \
flux_tang_matrix_bottom, \
penalty_tang_matrix_bottom, \
stability_tang_matrix_bottom = domain.boundary['trim_bottom'].integrate( [ flux_norm_matrix,
penalty_norm_matrix,
stability_norm_matrix,
flux_tang_matrix,
penalty_tang_matrix,
stability_tang_matrix
], geometry=geom, ischeme='gauss'+str(2*degree) )
log.user('finished bottom integration in',time.time()-t_old)
t_old = time.time()
flux_norm_matrix_left, \
flux_norm_vector, \
penalty_norm_matrix_left, \
penalty_norm_vector, \
stability_norm_matrix_left, \
flux_tang_matrix_left, \
flux_tang_vector, \
penalty_tang_matrix_left, \
penalty_tang_vector, \
stability_tang_matrix_left = domain.boundary['trim_left' ].integrate( [ flux_norm_matrix,
flux_norm_vector,
penalty_norm_matrix,
penalty_norm_vector,
stability_norm_matrix,
flux_tang_matrix,
flux_tang_vector,
penalty_tang_matrix,
penalty_tang_vector,
stability_tang_matrix
], geometry=geom, ischeme='gauss'+str(2*degree) )
log.user('finished left integration in',time.time()-t_old)
flux_norm_matrix = flux_norm_matrix_bottom .toscipy() + flux_norm_matrix_left .toscipy() + flux_tang_matrix_bottom .toscipy()
penalty_norm_matrix = penalty_norm_matrix_bottom .toscipy() + penalty_norm_matrix_left .toscipy() + penalty_tang_matrix_bottom .toscipy()
stability_norm_matrix = stability_norm_matrix_bottom.toscipy() + stability_norm_matrix_left.toscipy() + stability_tang_matrix_bottom.toscipy()
flux_tang_matrix = flux_tang_matrix_bottom .toscipy() + flux_tang_matrix_left .toscipy()
penalty_tang_matrix = penalty_tang_matrix_bottom .toscipy() + penalty_tang_matrix_left .toscipy()
stability_tang_matrix = stability_tang_matrix_bottom.toscipy() + stability_tang_matrix_left.toscipy()
#Remove singularities
p = numpy.ones(funcsp.shape[0],dtype=bool)
p[ x_index1 ] = False
p[ x_index2 ] = False
p[ y_index1 ] = False
p[ y_index2 ] = False
p[ z_index1 ] = False
p[ z_index2 ] = False
t_old = time.time()
#Eigenvalue problem
A_matrix = volume_matrix [:,p][p,:]
E_norm_matrix = stability_norm_matrix[:,p][p,:]
E_tang_matrix = stability_tang_matrix[:,p][p,:]
D = precon(A_matrix).toarray()
log.user('assembled eigenvalue problems in',time.time()-t_old)
t_old = time.time()
beta_norm = 2*max(scipy.linalg.eigh(D.dot(E_norm_matrix.dot(D.T)),D.dot(A_matrix.dot(D.T)),eigvals_only=True))
beta_tang = 2*max(scipy.linalg.eigh(D.dot(E_tang_matrix.dot(D.T)),D.dot(A_matrix.dot(D.T)),eigvals_only=True))
# (beta_norm,),dummy = 2*scipy.sparse.linalg.eigsh( A=D.dot(E_norm_matrix.dot(D.T)), k=1, M=D.dot(A_matrix.dot(D.T)), which='LM' )
# (beta_tang,),dummy = 2*scipy.sparse.linalg.eigsh( A=D.dot(E_tang_matrix.dot(D.T)), k=1, M=D.dot(A_matrix.dot(D.T)), which='LM' )
log.user('solved eigenvalue problems in',time.time()-t_old)
#Assemble
boundary_norm_matrix = flux_norm_matrix + beta_norm * penalty_norm_matrix
boundary_norm_vector = flux_norm_vector + beta_norm * penalty_norm_vector
boundary_tang_matrix = flux_tang_matrix + beta_tang * penalty_tang_matrix
boundary_tang_vector = flux_tang_vector + beta_tang * penalty_tang_vector
##############################
# Solve #
##############################
#Remove singularities
p_norm = numpy.ones(funcsp.shape[0],dtype=bool)
p_tang = numpy.ones(funcsp.shape[0],dtype=bool)
p_tang[ x_index2 ] = False
t_old = time.time()
#Assemble
full_norm_matrix = (volume_matrix + boundary_norm_matrix)[:,p_norm][p_norm,:]
full_norm_vector = boundary_norm_vector[p_norm]
full_tang_matrix = (volume_matrix + boundary_tang_matrix)[:,p_tang][p_tang,:]
full_tang_vector = boundary_tang_vector[p_tang]
log.user('assembled full systems in',time.time()-t_old)
#Solve
t_old = time.time()
D = scipy.sparse.csr_matrix(precon(full_norm_matrix))
lhs_free = D.T.dot(scipy.sparse.linalg.cg(D.dot(full_norm_matrix.dot(D.T)),D.dot(full_norm_vector),tol=solvetol,maxiter=10000)[0])
lhs_norm = numpy.zeros(funcsp.shape[0]); lhs_norm[p_norm] = lhs_free
log.user('solved norm in',time.time()-t_old)
t_old = time.time()
D = scipy.sparse.csr_matrix(precon(full_tang_matrix))
lhs_free = D.T.dot(scipy.sparse.linalg.cg(D.dot(full_tang_matrix.dot(D.T)),D.dot(full_tang_vector),tol=solvetol,maxiter=10000)[0])
lhs_tang = numpy.zeros(funcsp.shape[0]); lhs_tang[p_tang] = lhs_free
log.user('solved tang in',time.time()-t_old)
##############################
# Postprocessing #
##############################
#Convergence check
energy_norm_v[h_n] = lhs_norm.dot(volume_matrix.dot(lhs_norm))
energy_tang_v[h_n] = lhs_tang.dot(volume_matrix.dot(lhs_tang))
if h_n>0:
with plot.PyPlot('energy_normal_' ) as plt:
plt.loglog(h_v[:h_n],numpy.abs(energy_norm_v[:h_n]-energy_norm_v[1:h_n+1]),c='r',marker='o',markersize=10,markeredgecolor='r',linewidth=0,fillstyle='full')
plt.loglog(h_v[:h_n],numpy.abs(energy_norm_v[h_n-1]-energy_norm_v[h_n])/h_v[h_n-1]**(2*degree)*h_v[:h_n]**(2*degree),c='k',linewidth=1,)
plt.grid(True)
#plt.legend(fontsize=20)
plt.xlabel(r'$h$',fontsize=20)
plt.ylabel(r'$E(v^h)$',fontsize=20)
log.user('normal' ,numpy.abs(energy_norm_v[:h_n]-energy_norm_v[1:h_n+1]))
with plot.PyPlot('energy_tangent_') as plt:
plt.loglog(h_v[:h_n],numpy.abs(energy_tang_v[:h_n]-energy_tang_v[1:h_n+1]),c='r',marker='o',markersize=10,markeredgecolor='r',linewidth=0,fillstyle='full')
plt.loglog(h_v[:h_n],numpy.abs(energy_tang_v[h_n-1]-energy_tang_v[h_n])/h_v[h_n-1]**(2*degree)*h_v[:h_n]**(2*degree),c='k',linewidth=1,)
plt.grid(True)
#plt.legend(fontsize=20)
plt.xlabel(r'$h$',fontsize=20)
plt.ylabel(r'$E(v^h)$',fontsize=20)
log.user('tangent',numpy.abs(energy_tang_v[:h_n]-energy_tang_v[1:h_n+1]))
#Plotting
plotmesh( domain, geom, funcsp.dot(lhs_norm), stresssp.dot(lhs_norm).dotnorm(geom), (stresssp.dot(lhs_norm)**2).sum([-1,-2]), 'h_'+str(int(h_n))+'_normal' )
plotmesh( domain, geom, funcsp.dot(lhs_tang), stresssp.dot(lhs_tang).dotnorm(geom), (stresssp.dot(lhs_tang)**2).sum([-1,-2]), 'h_'+str(int(h_n))+'_tangent' )
return
util.run( main )