modellib/py_src/mechanics.py


##############################################################################
#
# Copyright (c) 2003-2020 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
# Development from 2019 by School of Earth and Environmental Sciences
#
##############################################################################

from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2020 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

#from esys.escript import *
from esys.escript.modelframe import Model,IterationDivergenceError
from esys.escript.linearPDEs import LinearPDE

class Mechanics(Model):
      """
      base class for mechanics models in updated lagrangean framework

      :note: Instance variable domain - domain (in)
      :note: Instance variable internal_force - =Data()
      :note: Instance variable external_force - =Data()
      :note: Instance variable prescribed_velocity - =Data()
      :note: Instance variable location_prescribed_velocity - =Data()
      :note: Instance variable temperature -  = None
      :note: Instance variable expansion_coefficient -  = 0.
      :note: Instance variable bulk_modulus - =1.
      :note: Instance variable shear_modulus - =1.
      :note: Instance variable rel_tol - =1.e-3
      :note: Instance variable abs_tol - =1.e-15
      :note: Instance variable max_iter - =10
      :note: Instance variable displacement - =None
      :note: Instance variable stress - =None
      """
      SAFTY_FACTOR_ITERATION=1./100.
      def __init__(self,**kwargs):
         """
         set up the model
         
         :keyword debug: debug flag
         :type debug: ``bool``
         """
         super(Mechanics, self).__init__(self,**kwargs)
         self.declareParameter(domain=None, \
                               displacement=None, \
                               stress=None, \
                               velocity=None, \
                               internal_force=None, \
                               external_force=None, \
                               prescribed_velocity=None, \
                               location_prescribed_velocity=None, \
                               temperature = None, \
                               expansion_coefficient = 0., \
                               bulk_modulus=2., \
                               shear_modulus=1., \
                               rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
         self.__iter=0

      def doInitialization(self):
           """
           initialize model
           """
           if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
           if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
           if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
           if not self.internal_force: self.internal_force = Data()
           if not self.external_force: self.external_force = Data()
           if not self.prescribed_velocity: self.prescribed_velocity = Data()
           if not self.location_prescribed_velocity: self.location_prescribed_velocity =Data()
           # save the old values:
           self.__stress_safe=self.stress
           self.__temperature_safe=self.temperature
           self.__displacement_safe=self.displacement
           self.__velocity_safe=self.velocity
           self.__velocity_old=None
           self.__old_dt=None
           self.__very_old_dt=None
           # get node cooridnates and apply initial displacement
           self.__x=self.domain.getX()
           self.domain.setX(self.__x+self.displacement)
           # open PDE:
           self.__pde=LinearPDE(self.domain)
           self.__pde.setSolverMethod(self.__pde.DIRECT)
           self.__solver_options=self.__pde.getSolverOptions()
           self.__solver_options.setSolverMethod(self.__solver_options.DIRECT)
           self.__solver_options.setVerbosity(self.debug)

           # self.__pde.setSymmetryOn()

      def doStepPreprocessing(self,dt):
            """
            step up pressure iteration

            if run within a time dependend problem extrapolation of pressure from previous time steps is used to
            get an initial guess (that needs some work!!!!!!!)
            """
            # reset iteration counters:
            self.__iter=0
            self.__diff=self.UNDEF_DT
            # set initial guesses for the iteration:
            self.displacement=self.__displacement_safe
            self.stress=self.__stress_safe
            self.velocity=self.__velocity_safe
            # update geometry
            self.domain.setX(self.__x+self.displacement)

      def doStep(self,dt):
          """
          """
          self.__iter+=1
          k3=kronecker(self.domain)
          # set new thermal stress increment
          if self.temperature == None:
             self.deps_th=0.
          else:
             self.deps_th=self.expansion_coefficient*(self.temperature-self.__temperature_safe)
          # set PDE coefficients:
          self.__pde.setValue(A=self.S)
          self.__pde.setValue(X=-self.stress-self.bulk_modulus*self.deps_th*k3)
          if self.internal_force: self.__pde.setValue(Y=self.internal_force)
          if self.external_force: self.__pde.setValue(y=self.external_force)
          self.__pde.setValue(q=self.location_prescribed_velocity, \
                              r=Data())
          if not self.prescribed_velocity.isEmpty() and self.__iter==1:
               self.__pde.setValue(r=dt*self.prescribed_velocity)
          # solve the PDE:
          self.__solver_options.setTolerance(self.rel_tol**2)
          self.du=self.__pde.getSolution()
          # update geometry
          self.displacement=self.displacement+self.du
          self.domain.setX(self.__x+self.displacement)
          self.velocity=(self.displacement-self.__displacement_safe)/dt

          if self.debug:
             for i in range(self.domain.getDim()):
                self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))
             for i in range(self.domain.getDim()):
                self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
             for i in range(self.domain.getDim()):
                self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
          self.__stress_last=self.stress

      def terminateIteration(self):
          """iteration is terminateIterationd if relative pressure change is less than rel_tol"""
          if self.__iter>self.max_iter:
              raise IterationDivergenceError("Maximum number of iterations steps reached")
          if self.__iter==0:
             self.__diff=self.UNDEF_DT
          else:
             self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff
             s_sup=Lsup(self.stress)
             self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
             if self.__iter>2 and diff_safe<self.__diff:
                 raise IterationDivergenceError("no improvement in stress iteration")
             return self.__diff<=self.rel_tol*self.SAFTY_FACTOR_ITERATION*s_sup+self.abs_tol

      def doStepPostprocessing(self,dt):
           """
           accept all the values:
           """
           self.__displacement_safe=self.displacement
           self.__temperature_safe=self.temperature
           self.__stress_safe=self.stress
           self.__velocity_safe=self.velocity

      def getSafeTimeStepSize(self,dt):
           """
           returns new step size
           """
           a=sup(length(self.velocity)/self.domain.getSize())
           if a>0:
              return 1./a
           else:
              return self.UNDEF_DT


class DruckerPrager(Mechanics):
      """

      """

      def __init__(self,**kwargs):
           """
           set up model
           """
           super(DruckerPrager, self).__init__(**kwargs)
           self.declareParameter(plastic_stress=0.,
                                 hardening=0.,
                                 friction_parameter=0.,
                                 dilatancy_parameter=0.,
                                 shear_length=1.e15)
      def doInitialization(self):
          """
          """
          super(DruckerPrager, self).doInitialization()
          self.__plastic_stress_safe=self.plastic_stress
          self.__shear_length_safe=self.shear_length
          self.__hardening_safe=self.hardening
          self.__chi_safe=0
          self.__tau_safe=0

      def doStepPreprocessing(self,dt):
          """
          """
          super(DruckerPrager, self).doStepPreprocessing(dt)
          # set initial guess for iteration:
          self.shear_length=self.__shear_length_safe
          self.plastic_stress=self.__plastic_stress_safe
          self.hardening=self.__hardening_safe
          self.__chi=self.__chi_safe
          self.__tau=self.__tau_safe

      def doStep(self,dt):
          # set new tangential operator:
          self.setTangentialTensor()
          # do the update step:
          super(DruckerPrager, self).doStep(dt)
          # update stresses:
          self.setStress()

      def doStepPostprocessing(self,dt):
          super(DruckerPrager, self).doStepPostprocessing(dt)
          self.__plastic_stress_safe=self.plastic_stress
          self.__shear_length_safe=self.shear_length
          self.__hardening_safe=self.hardening
          self.__chi_safe=self.__chi
          self.__tau_safe=self.__tau

      def setStress(self):
           d=self.domain.getDim()
           G=self.shear_modulus
           K=self.bulk_modulus
           alpha=self.friction_parameter
           beta=self.dilatancy_parameter
           h=self.hardening
           k3=kronecker(self.domain)
           # elastic trial stress:
           g=grad(self.du)
           D=symmetric(g)
           W=nonsymmetric(g)
           s_e=self.stress+K*self.deps_th*k3+ \
                      2*G*D+(K-2./3*G)*trace(D)*k3 \
                      +2*symmetric(matrix_mult(W,self.stress))
           p_e=-1./d*trace(s_e)
           s_e_dev=s_e+p_e*k3
           tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
           # yield conditon for elastic trial stress:
           F=tau_e-alpha*p_e-self.shear_length
           self.__chi=whereNonNegative(F+(self.rel_tol*(self.SAFTY_FACTOR_ITERATION)**2)*self.shear_length)
           # plastic stress increment:
           l=self.__chi*F/(h+G+beta*K)
           self.__tau=tau_e-G*l
           self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*k3
           self.plastic_stress=self.plastic_stress+l
           # update hardening
           self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))

      def setTangentialTensor(self):
           d=self.domain.getDim()
           G=self.shear_modulus
           K=self.bulk_modulus
           alpha=self.friction_parameter
           beta=self.dilatancy_parameter
           tau_Y=self.shear_length
           chi=self.__chi
           tau=self.__tau
           h=self.hardening
           k3=kronecker(Function(self.domain))

           sXk3=outer(self.stress,k3)
           k3Xk3=outer(k3,k3)
           s_dev=self.stress-trace(self.stress)*(k3/d)
           tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))

           self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
                 + (K-2./3*G)*k3Xk3 \
                 + (sXk3-swap_axes(swap_axes(sXk3,1,2),2,3))   \
                 + 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3)   \
                        -swap_axes(swap_axes(sXk3,0,3),2,3)   \
                        -swap_axes(sXk3,1,2)                  \
                        +swap_axes(sXk3,1,3)                ) \
                 - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)
1	ksteube	1809
2	jfenwick	3981	##############################################################################
3	ksteube	1312	#
4	uqaeller	6939	# Copyright (c) 2003-2020 by The University of Queensland
5	jfenwick	3981	# http://www.uq.edu.au
6	ksteube	1312	#
7	ksteube	1809	# Primary Business: Queensland, Australia
8	jfenwick	6112	# Licensed under the Apache License, version 2.0
9			# http://www.apache.org/licenses/LICENSE-2.0
10	ksteube	1312	#
11	jfenwick	3981	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12	jfenwick	4657	# Development 2012-2013 by School of Earth Sciences
13			# Development from 2014 by Centre for Geoscience Computing (GeoComp)
14	uqaeller	6939	# Development from 2019 by School of Earth and Environmental Sciences
15	jfenwick	3981	#
16			##############################################################################
17	gross	814
18	sshaw	5706	from __future__ import print_function, division
19
20	uqaeller	6939	__copyright__="""Copyright (c) 2003-2020 by The University of Queensland
21	jfenwick	3981	http://www.uq.edu.au
22	ksteube	1809	Primary Business: Queensland, Australia"""
23	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
24			http://www.apache.org/licenses/LICENSE-2.0"""
25	jfenwick	2344	__url__="https://launchpad.net/escript-finley"
26	gross	814
27	jfenwick	3432	#from esys.escript import *
28	gross	814	from esys.escript.modelframe import Model,IterationDivergenceError
29			from esys.escript.linearPDEs import LinearPDE
30
31			class Mechanics(Model):
32			"""
33			base class for mechanics models in updated lagrangean framework
34
35	jfenwick	6647	:note: Instance variable domain - domain (in)
36			:note: Instance variable internal_force - =Data()
37			:note: Instance variable external_force - =Data()
38			:note: Instance variable prescribed_velocity - =Data()
39			:note: Instance variable location_prescribed_velocity - =Data()
40			:note: Instance variable temperature - = None
41			:note: Instance variable expansion_coefficient - = 0.
42			:note: Instance variable bulk_modulus - =1.
43			:note: Instance variable shear_modulus - =1.
44			:note: Instance variable rel_tol - =1.e-3
45			:note: Instance variable abs_tol - =1.e-15
46			:note: Instance variable max_iter - =10
47			:note: Instance variable displacement - =None
48			:note: Instance variable stress - =None
49	gross	814	"""
50	gross	831	SAFTY_FACTOR_ITERATION=1./100.
51	gross	918	def __init__(self,**kwargs):
52	gross	814	"""
53			set up the model
54
55	jfenwick	2625	:keyword debug: debug flag
56			:type debug: ``bool``
57	gross	814	"""
58	gross	918	super(Mechanics, self).__init__(self,**kwargs)
59	gross	814	self.declareParameter(domain=None, \
60			displacement=None, \
61			stress=None, \
62			velocity=None, \
63	gross	918	internal_force=None, \
64			external_force=None, \
65			prescribed_velocity=None, \
66			location_prescribed_velocity=None, \
67	gross	814	temperature = None, \
68			expansion_coefficient = 0., \
69	gross	821	bulk_modulus=2., \
70	gross	814	shear_modulus=1., \
71			rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
72			self.__iter=0
73
74			def doInitialization(self):
75			"""
76			initialize model
77			"""
78			if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
79			if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
80			if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
81	gross	918	if not self.internal_force: self.internal_force = Data()
82			if not self.external_force: self.external_force = Data()
83			if not self.prescribed_velocity: self.prescribed_velocity = Data()
84			if not self.location_prescribed_velocity: self.location_prescribed_velocity =Data()
85	gross	831	# save the old values:
86	gross	832	self.__stress_safe=self.stress
87			self.__temperature_safe=self.temperature
88			self.__displacement_safe=self.displacement
89			self.__velocity_safe=self.velocity
90	gross	838	self.__velocity_old=None
91			self.__old_dt=None
92			self.__very_old_dt=None
93	gross	831	# get node cooridnates and apply initial displacement
94	gross	829	self.__x=self.domain.getX()
95			self.domain.setX(self.__x+self.displacement)
96			# open PDE:
97			self.__pde=LinearPDE(self.domain)
98	gross	831	self.__pde.setSolverMethod(self.__pde.DIRECT)
99	gross	3014	self.__solver_options=self.__pde.getSolverOptions()
100			self.__solver_options.setSolverMethod(self.__solver_options.DIRECT)
101			self.__solver_options.setVerbosity(self.debug)
102
103	gross	841	# self.__pde.setSymmetryOn()
104	gross	831
105	gross	814	def doStepPreprocessing(self,dt):
106			"""
107			step up pressure iteration
108
109			if run within a time dependend problem extrapolation of pressure from previous time steps is used to
110			get an initial guess (that needs some work!!!!!!!)
111			"""
112	gross	831	# reset iteration counters:
113	gross	814	self.__iter=0
114			self.__diff=self.UNDEF_DT
115	gross	834	# set initial guesses for the iteration:
116	gross	832	self.displacement=self.__displacement_safe
117			self.stress=self.__stress_safe
118			self.velocity=self.__velocity_safe
119	gross	829	# update geometry
120			self.domain.setX(self.__x+self.displacement)
121
122	gross	814	def doStep(self,dt):
123			"""
124			"""
125	gross	831	self.__iter+=1
126	gross	814	k3=kronecker(self.domain)
127			# set new thermal stress increment
128	gross	918	if self.temperature == None:
129			self.deps_th=0.
130	gross	814	else:
131	gross	918	self.deps_th=self.expansion_coefficient*(self.temperature-self.__temperature_safe)
132	gross	814	# set PDE coefficients:
133			self.__pde.setValue(A=self.S)
134	gross	831	self.__pde.setValue(X=-self.stress-self.bulk_modulusself.deps_thk3)
135	gross	814	if self.internal_force: self.__pde.setValue(Y=self.internal_force)
136			if self.external_force: self.__pde.setValue(y=self.external_force)
137	gross	831	self.__pde.setValue(q=self.location_prescribed_velocity, \
138			r=Data())
139			if not self.prescribed_velocity.isEmpty() and self.__iter==1:
140			self.__pde.setValue(r=dt*self.prescribed_velocity)
141	gross	814	# solve the PDE:
142	gross	3014	self.__solver_options.setTolerance(self.rel_tol**2)
143			self.du=self.__pde.getSolution()
144	gross	829	# update geometry
145	gross	831	self.displacement=self.displacement+self.du
146	gross	834	self.domain.setX(self.__x+self.displacement)
147	gross	832	self.velocity=(self.displacement-self.__displacement_safe)/dt
148	gross	829
149	gross	821	if self.debug:
150			for i in range(self.domain.getDim()):
151	gross	831	self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))
152	gross	829	for i in range(self.domain.getDim()):
153			self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
154	gross	832	for i in range(self.domain.getDim()):
155			self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
156	gross	831	self.__stress_last=self.stress
157	gross	829
158			def terminateIteration(self):
159	caltinay	4285	"""iteration is terminateIterationd if relative pressure change is less than rel_tol"""
160	gross	814	if self.__iter>self.max_iter:
161	jfenwick	3892	raise IterationDivergenceError("Maximum number of iterations steps reached")
162	gross	831	if self.__iter==0:
163			self.__diff=self.UNDEF_DT
164			else:
165	gross	832	self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff
166	gross	831	s_sup=Lsup(self.stress)
167			self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
168	gross	834	if self.__iter>2 and diff_safe<self.__diff:
169	jfenwick	3892	raise IterationDivergenceError("no improvement in stress iteration")
170	gross	831	return self.__diff<=self.rel_tolself.SAFTY_FACTOR_ITERATIONs_sup+self.abs_tol
171	gross	814
172			def doStepPostprocessing(self,dt):
173			"""
174			accept all the values:
175			"""
176	gross	832	self.__displacement_safe=self.displacement
177			self.__temperature_safe=self.temperature
178			self.__stress_safe=self.stress
179			self.__velocity_safe=self.velocity
180	gross	814
181			def getSafeTimeStepSize(self,dt):
182			"""
183			returns new step size
184	gross	832	"""
185	gross	841	a=sup(length(self.velocity)/self.domain.getSize())
186			if a>0:
187			return 1./a
188	gross	814	else:
189	gross	832	return self.UNDEF_DT
190	gross	814
191
192
193			class DruckerPrager(Mechanics):
194			"""
195
196			"""
197
198	gross	918	def __init__(self,**kwargs):
199	gross	814	"""
200			set up model
201			"""
202	gross	918	super(DruckerPrager, self).__init__(**kwargs)
203	gross	814	self.declareParameter(plastic_stress=0.,
204	gross	831	hardening=0.,
205	gross	814	friction_parameter=0.,
206			dilatancy_parameter=0.,
207	gross	829	shear_length=1.e15)
208	gross	814	def doInitialization(self):
209	gross	829	"""
210			"""
211			super(DruckerPrager, self).doInitialization()
212	gross	832	self.__plastic_stress_safe=self.plastic_stress
213			self.__shear_length_safe=self.shear_length
214			self.__hardening_safe=self.hardening
215			self.__chi_safe=0
216			self.__tau_safe=0
217	gross	814
218	gross	831	def doStepPreprocessing(self,dt):
219	gross	829	"""
220			"""
221			super(DruckerPrager, self).doStepPreprocessing(dt)
222			# set initial guess for iteration:
223	gross	832	self.shear_length=self.__shear_length_safe
224			self.plastic_stress=self.__plastic_stress_safe
225			self.hardening=self.__hardening_safe
226			self.__chi=self.__chi_safe
227			self.__tau=self.__tau_safe
228	gross	814
229	gross	831	def doStep(self,dt):
230	gross	829	# set new tangential operator:
231	gross	831	self.setTangentialTensor()
232	gross	829	# do the update step:
233			super(DruckerPrager, self).doStep(dt)
234			# update stresses:
235	gross	831	self.setStress()
236	gross	829
237			def doStepPostprocessing(self,dt):
238			super(DruckerPrager, self).doStepPostprocessing(dt)
239	gross	832	self.__plastic_stress_safe=self.plastic_stress
240			self.__shear_length_safe=self.shear_length
241			self.__hardening_safe=self.hardening
242			self.__chi_safe=self.__chi
243			self.__tau_safe=self.__tau
244	gross	829
245			def setStress(self):
246	gross	831	d=self.domain.getDim()
247	gross	814	G=self.shear_modulus
248			K=self.bulk_modulus
249			alpha=self.friction_parameter
250			beta=self.dilatancy_parameter
251	gross	829	h=self.hardening
252			k3=kronecker(self.domain)
253			# elastic trial stress:
254	gross	831	g=grad(self.du)
255	gross	829	D=symmetric(g)
256			W=nonsymmetric(g)
257	gross	831	s_e=self.stress+Kself.deps_thk3+ \
258			2GD+(K-2./3G)trace(D)*k3 \
259			+2*symmetric(matrix_mult(W,self.stress))
260			p_e=-1./d*trace(s_e)
261	gross	829	s_e_dev=s_e+p_e*k3
262			tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
263			# yield conditon for elastic trial stress:
264			F=tau_e-alpha*p_e-self.shear_length
265	gross	838	self.__chi=whereNonNegative(F+(self.rel_tol(self.SAFTY_FACTOR_ITERATION)2)self.shear_length)
266	gross	829	# plastic stress increment:
267			l=self.__chiF/(h+G+betaK)
268			self.__tau=tau_e-G*l
269	gross	832	self.stress=self.__tau/(tau_e+self.abs_tolwhereZero(tau_e,self.abs_tol))s_e_dev-(p_e+lbetaK)*k3
270	gross	831	self.plastic_stress=self.plastic_stress+l
271	gross	829	# update hardening
272	gross	832	self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))
273	gross	829
274			def setTangentialTensor(self):
275	gross	831	d=self.domain.getDim()
276	gross	829	G=self.shear_modulus
277			K=self.bulk_modulus
278			alpha=self.friction_parameter
279			beta=self.dilatancy_parameter
280	gross	814	tau_Y=self.shear_length
281	gross	829	chi=self.__chi
282	gross	831	tau=self.__tau
283	gross	825	h=self.hardening
284	gross	832	k3=kronecker(Function(self.domain))
285	gross	814
286	gross	831	sXk3=outer(self.stress,k3)
287	gross	814	k3Xk3=outer(k3,k3)
288	gross	831	s_dev=self.stress-trace(self.stress)*(k3/d)
289			tmp=Gs_dev/(tau+self.abs_tolwhereZero(tau,self.abs_tol))
290	gross	836
291	gross	834	self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
292	gross	831	+ (K-2./3G)k3Xk3 \
293	gross	841	+ (sXk3-swap_axes(swap_axes(sXk3,1,2),2,3)) \
294	gross	836	+ 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3) \
295			-swap_axes(swap_axes(sXk3,0,3),2,3) \
296			-swap_axes(sXk3,1,2) \
297			+swap_axes(sXk3,1,3) ) \
298	gross	831	- outer(chi/(h+G+alphabetaK)(tmp+betaKk3),tmp+alphaK*k3)