agd.Metrics.Seismic.hooke API documentation

Hooke(hooke, *args, **kwargs) View Source

37	def __init__(self,hooke,*args,**kwargs):
38		super(Hooke,self).__init__(*args,**kwargs)
39		self.hooke = hooke
40		self._to_common_field()

hooke

def is_definite(self): View Source

42	def is_definite(self):
43		return Riemann(self.hooke).is_definite()

Attempts to check wether the data defines a mathematically valid norm.

vdim View Source

58	@property
59	def vdim(self):
60		return self._vdim(len(self.hooke))

The ambient vector space dimension, often denoted $d$ in mathematical formulas.

shape View Source

62	@property
63	def shape(self): return self.hooke.shape[2:]

Dimensions of the underlying domain. Expected to be the empty tuple, or a tuple of length vdim.

def model_HFM(self): View Source

65	def model_HFM(self):
66		d = self.vdim
67		suffix = "" if self.inverse_transformation is None else "Topographic"
68		return f"Seismic{suffix}{d}"

The name of the 'model' for parameter, as input to the HFM library.

def flatten(self): View Source

70	def flatten(self):
71		hooke = misc.flatten_symmetric_matrix(self.hooke)
72		if self.inverse_transformation is None: 
73			return hooke
74		else: 
75			inv_trans= self.inverse_transformation.reshape((self.vdim**2,)+self.shape)
76			return np.concatenate((hooke,inv_trans),axis=0)

Flattens and concatenate the member fields into a single array.

@classmethod

def expand(cls, arr): View Source

78	@classmethod
79	def expand(cls,arr):
80		return cls(misc.expand_symmetric_matrix(arr))

Inverse of the flatten member function. Turns a suitable array into a metric.

def with_cost(self, cost): View Source

87	def with_cost(self,cost):
88		other = copy.copy(self)
89		hooke,cost = fd.common_field((self.hooke,cost),depths=(2,0))
90		other.hooke = hooke / cost**2
91		return other

Produces a norm $N'$ obeying $N'(x) = N(cost*x)$.

def extract_xz(self): View Source

118	def extract_xz(self):
119		"""
120		Extract a two dimensional Hooke tensor from a three dimensional one, 
121		corresponding to a slice through the X and Z axes.
122		"""
123		assert self.vdim==3
124		h=self.hooke
125		return Hooke(ad.asarray([ 
126			[h[0,0], h[0,2], h[0,4] ],
127			[h[2,0], h[2,2], h[2,4] ],
128			[h[4,0], h[4,2], h[4,4] ]
129			]))

Extract a two dimensional Hooke tensor from a three dimensional one, corresponding to a slice through the X and Z axes.

@classmethod

def from_VTI_2(cls, Vp, Vs, ε, δ): View Source

131	@classmethod
132	def from_VTI_2(cls,Vp,Vs,ε,δ):
133		"""
134		X,Z slice of a Vertical Transverse Isotropic medium
135		based on Thomsen parameters
136		"""
137		c33=Vp**2
138		c44=Vs**2
139		c11=c33*(1+2*ε)
140		c13=-c44+np.sqrt( (c33-c44)**2+2*δ*c33*(c33-c44) )
141		return cls.from_hexagonal(c11,None,c13,c33,c44,vdim=2)

X,Z slice of a Vertical Transverse Isotropic medium based on Thomsen parameters

@classmethod

def from_Ellipse(cls, m): View Source

143	@classmethod
144	def from_Ellipse(cls,m):
145		"""
146		Rank deficient Hooke tensor,
147		equivalent, for pressure waves, to the Riemannian metric defined by $m ^ {-2}$.
148		Shear waves are infinitely slow.
149		"""
150		assert(len(m)==2)
151		a,b,c=m[0,0],m[1,1],m[0,1]
152		return Hooke(ad.asarray( [ [a*a, a*b,a*c], [a*b, b*b, b*c], [a*c, b*c, c*c] ] ))

Rank deficient Hooke tensor, equivalent, for pressure waves, to the Riemannian metric defined by $m ^ {-2}$. Shear waves are infinitely slow.

@classmethod

def from_cast(cls, metric): View Source

154	@classmethod
155	def from_cast(cls,metric): 
156		if isinstance(metric,cls):	return metric
157		riemann = Riemann.from_cast(metric)
158		
159		m = riemann.dual().m
160		assert not ad.is_ad(m)
161		from scipy.linalg import sqrtm
162		return cls.from_Ellipse(sqrtm(m))

Produces a metric by casting another metric of a compatible type.

def to_depth4(self): View Source

180	def to_depth4(self):
181		"""
182		Produces the full Hooke tensor, of shape
183		(vdim,vdim,vdim,vdim, n1,...,nk)
184		where vdim is the ambient space dimension.
185		"""
186		Voigt = self._Voigt
187		d = self.vdim
188		return ad.array([ [ [ [ self.hooke[Voigt[i,j],Voigt[k,l]]
189			for i in range(d)] for j in range(d)] for k in range(d)] for l in range(d)])

Produces the full Hooke tensor, of shape (vdim,vdim,vdim,vdim, n1,...,nk) where vdim is the ambient space dimension.

def rotate(self, r): View Source

191	def rotate(self,r):
192		other = copy.copy(self)
193		hooke,r = common_field((self.hooke,r),(2,2))
194		Voigti = self._Voigti
195		R = ad.array([ [ 
196			r[i0,i1]*r[j0,j1] if i1==j1 else (r[i0,i1]*r[j0,j1]+r[j0,i1]*r[i0,j1])
197			for (i0,j0) in Voigti] for (i1,j1) in Voigti])
198		other.hooke = lp.dot_AA(lp.transpose(R),lp.dot_AA(hooke,R))
199		return other

Rotation of the norm, by a given rotation matrix. The new unit ball is the direct image of the previous one.

def to_Mandel(self, a=1.4142135623730951): View Source

208	def to_Mandel(self,a=np.sqrt(2)):
209		r"""Introduces the $\sqrt 2$ and $2$ factors involved in Mandel's notation"""
210		return self.hooke*self._Mandel_factors(self.vdim,self.shape)

Introduces the $\sqrt 2$ and $2$ factors involved in Mandel's notation

@classmethod

def from_Mandel(cls, mandel, a=1.4142135623730951): View Source

212	@classmethod
213	def from_Mandel(cls,mandel,a=np.sqrt(2)):
214		r"""Removes the $\sqrt 2$ and $2$ factors involved in Mandel's notation"""
215		vdim = cls._vdim(len(mandel))
216		return Hooke(mandel/cls._Mandel_factors(vdim,mandel.shape[2:],a))

Removes the $\sqrt 2$ and $2$ factors involved in Mandel's notation

@classmethod

def from_orthorombic(cls, c11, c12, c13, c22, c23, c33, c44, c55, c66): View Source

218	@classmethod
219	def from_orthorombic(cls,c11,c12,c13,c22,c23,c33,c44,c55,c66):
220		z=0*c11 # np.zeros_like(c11) raises issue in combination with cupy
221		return cls(ad.array([
222		[c11,c12,c13,  z,  z,  z],
223		[c12,c22,c23,  z,  z,  z],
224		[c13,c23,c33,  z,  z,  z],
225		[  z,  z,  z,c44,  z,  z],
226		[  z,  z,  z,  z,c55,  z],
227		[  z,  z,  z,  z,  z,c66]
228		]))

@classmethod

def from_orthorombic2(cls, c11, c21, c22, c31, c32, c33, c44, c55, c66): View Source

230	@classmethod
231	def from_orthorombic2(cls,c11,c21,c22,c31,c32,c33,c44,c55,c66): 
232		"""Orthorombic medium with a different ordering of the first block coefficients"""
233		c12,c13,c23 = c21,c31,c32
234		return cls.from_orthorombic(c11,c12,c13,c22,c23,c33,c44,c55,c66)

Orthorombic medium with a different ordering of the first block coefficients

@classmethod

def from_tetragonal(cls, c11, c12, c13, c33, c44, c66): View Source

236	@classmethod
237	def from_tetragonal(cls,c11,c12,c13,c33,c44,c66):
238		c22,c23,c55 = c11,c13,c44
239		return cls.from_orthorombic(c11,c12,c13,c22,c23,c33,c44,c55,c66)

@classmethod

def from_hexagonal(cls, c11, c12, c13, c33, c44, vdim=3): View Source

241	@classmethod
242	def from_hexagonal(cls,c11,c12,c13,c33,c44,vdim=3):
243		if vdim==2:
244			z = 0*c11
245			return cls(ad.asarray( [ [c11,c13,z], [c13,c33,z], [z,z,c44] ] ))
246		c66 = (c11-c12)/2.
247		return cls.from_tetragonal(c11,c12,c13,c33,c44,c66)

@classmethod

def from_ThomsenElastic(cls, tem): View Source

249	@classmethod 
250	def from_ThomsenElastic(cls,tem):
251		"""Hooke tensor (m/s)^2 and density (g/cm^3)"""
252		if not isinstance(tem,ThomsenElasticMaterial): tem = ThomsenElasticMaterial(*tem)
253		hex,ρ = tem.to_hexagonal()
254		return cls.from_hexagonal(*hex),ρ

Hooke tensor (m/s)^2 and density (g/cm^3)

def to_orthorombic(self): View Source

256	def to_orthorombic(self):
257		"""Inverse function of from_orthorombic. No reconstruction check."""
258		assert self.vdim==3
259		return tuple(self.hooke[i,j] for i,j in 
260			((0,0),(0,1),(0,2),(1,1),(1,2),(2,2),(3,3),(4,4),(5,5)))

Inverse function of from_orthorombic. No reconstruction check.

def to_orthorombic2(self): View Source

261	def to_orthorombic2(self):
262		a,b,c,d,e,f,g,h,i = self.to_orthorombic()
263		return (a,b,d,c,e,f,g,h,i)

def to_tetragonal(self): View Source

264	def to_tetragonal(self):
265		a,b,c,_a,_c,d,e,_e,f = self.to_orthorombic()
266		return (a,b,c,d,e,f)

def to_hexagonal(self): View Source

267	def to_hexagonal(self):
268		a,b,c,d,e,_a_b_2 = self.to_tetragonal()
269		return (a,b,c,d,e)

def is_TTI(self, tol=None): View Source

271	def is_TTI(self,tol=None):
272		"""
273		Determine if the metric is in a TTI form.
274		"""
275		# Original code by F. Desquilbet, 2020
276		if tol is None: # small value (acts as zero for floats)
277			tol = max(1e-9, MA(hooke)*1e-12)
278
279		def small(arr): return np.max(np.abs(arr))<tol
280		is_Sym = small(hooke-lp.transpose(hooke)) # symmetrical
281
282		if metric.vdim==2:
283			return is_Sym and small(hooke[2,0]) and small(hooke[2,1])
284		if metric.vdim==3:
285			return (is_Sym 
286				and small((hooke[0,0]-hooke[0,1])/2-hooke[5,5])
287				and all(small(hooke[i,j]-hooke[k,l]) 
288					for ((i,j),(k,l)) in [((0,0),(1,1)), ((2,0),(2,1)), ((3,3),(4,4))])
289				and all(small(hooke[i,j]) 
290					for (i,j) in [(3,0),(4,0),(5,0),(3,1),(4,1),(5,1),
291					(3,2),(4,2),(5,2),(4,3),(5,3),(5,4)]) )

Determine if the metric is in a TTI form.

def dot_A(self, m, sym=True): View Source

316	def dot_A(self,m,sym=True):
317		"""
318		Dot product associated with a Hooke tensor, which turns a strain tensor epsilon
319		into a stress tensor sigma.
320
321		Input:
322		- m : the strain tensor.
323		"""
324		v,hooke = fd.common_field((self._Voigt_m2v(m,sym),self.hooke),(1,2))
325		w = lp.dot_AV(hooke,v)
326		return ad.array( ((w[0],w[2]),(w[2],w[1])) )

Dot product associated with a Hooke tensor, which turns a strain tensor epsilon into a stress tensor sigma.

Input:

m : the strain tensor.

def dot_AA(self, m1, m2=None, sym=True): View Source

328	def dot_AA(self,m1,m2=None,sym=True):
329		"""
330		Inner product associated with a Hooke tensor, on the space of symmetric matrices.
331
332		Inputs:
333		- m1 : first symmetric matrix
334		- m2 : second symmetric matrix. Defaults to m1.
335		"""
336		if m2 is None: 
337			v1,hooke = fd.common_field((self._Voigt_m2v(m1,sym),self.hooke),(1,2))
338			v2=v1
339		else: 
340			v1,v2,hooke = fd.common_field(
341				(self._Voigt_m2v(m1,sym),self._Voigt_m2v(m2,sym),self.hooke),(1,1,2))
342		return lp.dot_VV(v1,lp.dot_AV(hooke,v2))

Inner product associated with a Hooke tensor, on the space of symmetric matrices.

Inputs:

m1 : first symmetric matrix
m2 : second symmetric matrix. Defaults to m1.

def Selling(self): View Source

344	def Selling(self):
345		r"""
346		Returns a decomposition of the hooke tensor in the mathematical form
347		$$
348		hooke = \sum_i \rho_i m_i  m_i^\top,
349		$$
350		where $\rho_i$ is a non-negative coefficient, $m_i$ is symmetric nonzero and has 
351		integer entries, and $\sum_i \rho_i$ is maximal. 
352		"""
353		assert(self.inverse_transformation is None)
354		if self.vdim<=2: coefs,offsets = Selling.Decomposition(self.hooke)
355		else: 
356			from ... import Eikonal
357			coefs,offsets = Eikonal.VoronoiDecomposition(self.hooke)
358		if self.vdim==1: 
359			moffsets = np.expand_dims(offsets,axis=0)
360		elif self.vdim==2:
361			moffsets = ad.array(((offsets[0],offsets[2]),(offsets[2],offsets[1])))
362		elif self.vdim==3:
363			moffsets = ad.array(( #Voigt notation
364				(offsets[0],offsets[5],offsets[4]),
365				(offsets[5],offsets[1],offsets[3]),
366				(offsets[4],offsets[3],offsets[2])))
367		else :
368			raise ValueError("Unsupported dimension")
369		return coefs,moffsets

Returns a decomposition of the hooke tensor in the mathematical form $$ hooke = \sum_i \rho_i m_i m_i^\top, $$ where $\rho_i$ is a non-negative coefficient, $m_i$ is symmetric nonzero and has integer entries, and $\sum_i \rho_i$ is maximal.

def apply_transform(self): View Source

371	def apply_transform(self):
372		"""
373		Applies the transformation, if any stored, to the hooke tensor. 
374
375		CAUTION : this feature is required for some applications to elasticity, 
376		but is incompatible with the eikonal equation solver.
377		"""
378		r = self.inverse_transformation
379		if r is None: return self
380		other = copy.copy(self)
381		other.inverse_transformation = None
382		return other.rotate(lp.transpose(r))

Applies the transformation, if any stored, to the hooke tensor.

CAUTION : this feature is required for some applications to elasticity, but is incompatible with the eikonal equation solver.

@classmethod

def from_Lame(cls, λ, μ, vdim=2): View Source

384	@classmethod
385	def from_Lame(cls,λ,μ,vdim=2):
386		"""
387		Constructs a Hooke tensor from the Lame coefficients, in dimension 2 or 3.
388		Positive definite provided μ>0 and 2μ+dλ>0
389		"""
390		z,s = 0*λ,λ+2*μ
391		if vdim==2: return cls(ad.array([
392			[s,λ,z],
393			[λ,s,z],
394			[z,z,μ]]))
395		elif vdim==3: return cls(ad.array([
396			[s,λ,λ,z,z,z],
397			[λ,s,λ,z,z,z],
398			[λ,λ,s,z,z,z],
399			[z,z,z,μ,z,z],
400			[z,z,z,z,μ,z],
401			[z,z,z,z,z,μ]]))
402		else: raise ValueError("Unsupported dimension")

Constructs a Hooke tensor from the Lame coefficients, in dimension 2 or 3. Positive definite provided μ>0 and 2μ+dλ>0

def contract(self, w): View Source

404	def contract(self,w):
405		r"""Returns the contracted tensor $\sum_{j,l}c_{ijkl} w_j w_l$."""
406		voi = self._Voigt
407		hooke,w = fd.common_field((self.hooke,w),depths=(2,1))
408		def c(i,j,k,l): return hooke[voi[i,j],voi[k,l]]
409		d = self.vdim; assert len(w)==d
410		return ad.array([[
411			sum(c(i,j,k,l)*w[j]*w[l] for j in range(d) for l in range(d))
412			for i in range(d)] for k in range(d)])

Returns the contracted tensor $\sum_{j,l}c_{ijkl} w_j w_l$.

def waves(self, k, ρ): View Source

414	def waves(self,k,ρ):
415		"""Returns the pulsation and direction of the waves with the given wave vector."""
416		m = np.moveaxis(self.contract(k),(0,1),(-2,-1))
417		eVal,eVec = np.linalg.eigh(m)
418		eVec = np.moveaxis(eVec,(-2,-1),(0,1))
419		return np.sqrt(eVal/ρ),eVec

Returns the pulsation and direction of the waves with the given wave vector.

mica = (Hooke(array([[178. , 42.4, 14.5, 0. , 0. , 0. ], [ 42.4, 178. , 14.5, 0. , 0. , 0. ], [ 14.5, 14.5, 54.9, 0. , 0. , 0. ], [ 0. , 0. , 0. , 12.2, 0. , 0. ], [ 0. , 0. , 0. , 0. , 12.2, 0. ], [ 0. , 0. , 0. , 0. , 0. , 67.8]]), None, 6, (), 0.0), 2.79)

stishovite = (Hooke(array([[453, 211, 203, 0, 0, 0], [211, 453, 203, 0, 0, 0], [203, 203, 776, 0, 0, 0], [ 0, 0, 0, 252, 0, 0], [ 0, 0, 0, 0, 252, 0], [ 0, 0, 0, 0, 0, 302]]), None, 6, (), 0.0), 4.29)

olivine = (Hooke(array([[323.7, 66.4, 71.6, 0. , 0. , 0. ], [ 66.4, 197.6, 75.6, 0. , 0. , 0. ], [ 71.6, 75.6, 235.1, 0. , 0. , 0. ], [ 0. , 0. , 0. , 64.6, 0. , 0. ], [ 0. , 0. , 0. , 0. , 78.7, 0. ], [ 0. , 0. , 0. , 0. , 0. , 79. ]]), None, 6, (), 0.0), 3.311)

agd.Metrics.Seismic.hooke

Inherited Members