solvcon.kerpak.vslin

1 # -*- coding: UTF-8 -*- 2 # 3 # Copyright (C) 2010-2011 Yung-Yu Chen <yyc@solvcon.net>. 4 # 5 # This program is free software; you can redistribute it and/or modify 6 # it under the terms of the GNU General Public License as published by 7 # the Free Software Foundation; either version 2 of the License, or 8 # (at your option) any later version. 9 # 10 # This program is distributed in the hope that it will be useful, 11 # but WITHOUT ANY WARRANTY; without even the implied warranty of 12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 # GNU General Public License for more details. 14 # 15 # You should have received a copy of the GNU General Public License along 16 # with this program; if not, write to the Free Software Foundation, Inc., 17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 18 19 """ 20 Anisotropic velocity-stress equations solver using the linear CESE method with 21 CUDA support. 22 """ 23 24 from solvcon.gendata import TypeNameRegistry 25 from solvcon.anchor import Anchor 26 from solvcon.hook import BlockHook 27 from solvcon.kerpak.cuse import CUDA_RAISE_ON_FAIL, CuseBC 28 from solvcon.kerpak.lincuse import (LincuseSolver, LincuseCase, 29 PlaneWaveSolution)

30 31 ############################################################################### 32 # Solver. 33 ############################################################################### 34 35 -class VslinSolver(LincuseSolver):

36 """ 37 Basic elastic solver. 38 39 @ivar mtrldict: map from names to material objects. 40 @itype mtrldict: dict 41 @ivar mtrllist: list of all material objects. 42 @itype mtrllist: list 43 """ 44 from solvcon.dependency import getcdll 45 __clib_vslin_c = { 46 2: getcdll('vslin2d_c'), 47 3: getcdll('vslin3d_c'), 48 } 49 __clib_vslin_cu = { 50 2: getcdll('vslin2d_cu', raise_on_fail=CUDA_RAISE_ON_FAIL), 51 3: getcdll('vslin3d_cu', raise_on_fail=CUDA_RAISE_ON_FAIL), 52 } 53 del getcdll 54 @property

55 - def _clib_vslin_c(self):

56 return self.__clib_vslin_c[self.ndim]

57 @property

58 - def _clib_vslin_cu(self):

59 return self.__clib_vslin_cu[self.ndim]

60 @property

61 - def _gdlen_(self):

62 return 9 * 9 * self.ndim

63 - def __init__(self, *args, **kw):

64 super(VslinSolver, self).__init__(*args, **kw) 65 self.mtrldict = kw.pop('mtrldict', {}) 66 self.mtrllist = None

67 - def make_grpda(self):

68 self.mtrllist = self._build_mtrllist(self.grpnames, self.mtrldict) 69 for igrp in range(len(self.grpnames)): 70 mtrl = self.mtrllist[igrp] 71 jaco = self.grpda[igrp].reshape(self.neq, self.neq, self.ndim) 72 mjacos = mtrl.get_jacos() 73 for idm in range(self.ndim): 74 jaco[:,:,idm] = mjacos[idm,:,:]

75 @staticmethod

76 - def _build_mtrllist(grpnames, mtrldict):

77 """ 78 Build the material list out of the mapping dict. 79 80 @type grpnames: list 81 @param mtrldict: the map from names to material objects. 82 @type mtrldict: dict 83 @return: the list of material object. 84 @rtype: Material 85 """ 86 mtrllist = list() 87 default_mtuple = mtrldict.get(None, None) 88 for grpname in grpnames: 89 try: 90 mtrl = mtrldict.get(grpname, default_mtuple) 91 except KeyError, e: 92 args = e.args[:] 93 args.append('no material named %s in mtrldict'%grpname) 94 e.args = args 95 raise 96 mtrllist.append(mtrl) 97 return mtrllist

98

99 ############################################################################### 100 # Case. 101 ############################################################################### 102 103 -class VslinCase(LincuseCase):

104 """ 105 Case for anisotropic elastic solids. 106 """ 107 defdict = { 108 'execution.neq': 9, 109 'solver.solvertype': VslinSolver, 110 'solver.mtrldict': dict, 111 }

112 - def make_solver_keywords(self):

113 kw = super(VslinCase, self).make_solver_keywords() 114 # setup material mapper. 115 kw['mtrldict'] = self.solver.mtrldict 116 return kw

117

118 ############################################################################### 119 # Boundary conditions. 120 ############################################################################### 121 122 -class VslinBC(CuseBC):

123 """ 124 Basic BC class for elastic problems. 125 """ 126 from solvcon.dependency import getcdll 127 __clib_vslinb_c = { 128 2: getcdll('vslinb2d_c'), 129 3: getcdll('vslinb3d_c'), 130 } 131 __clib_vslinb_cu = { 132 2: getcdll('vslinb2d_cu', raise_on_fail=CUDA_RAISE_ON_FAIL), 133 3: getcdll('vslinb3d_cu', raise_on_fail=CUDA_RAISE_ON_FAIL), 134 } 135 del getcdll 136 @property

137 - def _clib_vslinb_c(self):

138 return self.__clib_vslinb_c[self.svr.ndim]

139 @property

140 - def _clib_vslinb_cu(self):

141 return self.__clib_vslinb_cu[self.svr.ndim]

142

143 -class VslinTraction(VslinBC):

144 vnames = [ 145 'bfcsys', 'tau1', 'tau2', 'tau3', 'freq', 'phase', 146 ] 147 vdefaults = { 148 'bfcsys': 0.0, 149 'tau1': 0.0, 'tau2': 0.0, 'tau3': 0.0, 'freq': 0.0, 'phase': 0.0, 150 } 151 _ghostgeom_ = 'compress'

152 - def soln(self):

153 from ctypes import byref, c_int 154 self._clib_boundcond.bound_traction_soln( 155 byref(self.svr.exd), 156 c_int(self.facn.shape[0]), 157 self.facn.ctypes._as_parameter_, 158 c_int(self.value.shape[1]), 159 self.value.ctypes._as_parameter_, 160 )

161 - def dsoln(self):

162 from ctypes import byref, c_int 163 self._clib_boundcond.bound_traction_dsoln( 164 byref(self.svr.exd), 165 c_int(self.facn.shape[0]), 166 self.facn.ctypes._as_parameter_, 167 )

168

169 -class VslinTractionFree(VslinBC):

170 _ghostgeom_ = 'mirror'

171 - def soln(self):

172 from ctypes import byref, c_int 173 self._clib_boundcond.bound_traction_free_soln( 174 byref(self.svr.exd), 175 c_int(self.facn.shape[0]), 176 self.facn.ctypes._as_parameter_, 177 )

178 - def dsoln(self):

179 from ctypes import byref, c_int 180 self._clib_boundcond.bound_traction_free_dsoln( 181 byref(self.svr.exd), 182 c_int(self.facn.shape[0]), 183 self.facn.ctypes._as_parameter_, 184 )

185

186 -class VslinTractionFree2(VslinBC):

187 _ghostgeom_ = 'mirror'

188 - def soln(self):

189 from ctypes import byref, c_int 190 self._clib_boundcond.bound_traction_free2_soln( 191 byref(self.svr.exd), 192 c_int(self.facn.shape[0]), 193 self.facn.ctypes._as_parameter_, 194 )

195 - def dsoln(self):

196 from ctypes import byref, c_int 197 self._clib_boundcond.bound_traction_free2_dsoln( 198 byref(self.svr.exd), 199 c_int(self.facn.shape[0]), 200 self.facn.ctypes._as_parameter_, 201 )

202

203 ################################################################################ 204 # Plane wave solution. 205 ################################################################################ 206 207 -class VslinPWSolution(PlaneWaveSolution):

208 - def _calc_eigen(self, **kw):

209 from numpy import sqrt 210 from numpy.linalg import eig 211 wvec = kw['wvec'] 212 mtrl = kw['mtrl'] 213 idx = kw['idx'] 214 nml = wvec/sqrt((wvec**2).sum()) 215 jacos = mtrl.get_jacos() 216 jaco = jacos[0] * nml[0] 217 for idm in range(1, len(nml)): 218 jaco += jacos[idm] * nml[idm] 219 evl, evc = eig(jaco) 220 srt = evl.argsort() 221 evl = evl[srt[idx]].real 222 evc = evc[:,srt[idx]].real 223 evc *= evc[0]/abs(evc[0]+1.e-200) 224 return evl, evc

225

226 ################################################################################ 227 # Anchor. 228 ################################################################################ 229 230 -class VslinOAnchor(Anchor):

231 """ 232 Calculate total energy, i.e., the summation of kinetic energy and strain 233 energy. 234 """

235 - def _calculate_physics(self):

236 from ctypes import byref 237 from numpy import empty 238 from numpy.linalg import inv 239 svr = self.svr 240 # input arrays. 241 rhos = empty(svr.ngroup, dtype=svr.fpdtype) 242 comps = empty((svr.ngroup, 6, 6), dtype=svr.fpdtype) 243 for igp in range(svr.ngroup): 244 mtrl = svr.mtrllist[igp] 245 rhos[igp] = mtrl.rho 246 comps[igp,:,:] = inv(mtrl.stiff).T 247 # output arrays. 248 svr._clib_vslin_c.calc_energy( 249 byref(svr.exd), 250 rhos.ctypes._as_parameter_, 251 comps.ctypes._as_parameter_, 252 svr.der['energy'].ctypes._as_parameter_, 253 )

254 - def provide(self):

255 from numpy import empty 256 svr = self.svr 257 svr.der['energy'] = empty(svr.ngstcell+svr.ncell, dtype=svr.fpdtype) 258 self._calculate_physics()

259 - def postfull(self):

260 self._calculate_physics()

261 262 ################################################################################ 263 # Material definition. 264 ################################################################################ 265 266 mltregy = TypeNameRegistry() # registry singleton.

267 -class MaterialMeta(type):

268 """ 269 Meta class for material class. 270 """

271 - def __new__(cls, name, bases, namespace):

272 newcls = super(MaterialMeta, cls).__new__(cls, name, bases, namespace) 273 # register. 274 mltregy.register(newcls) 275 return newcls

276

277 -class Material(object):

278 """ 279 Material properties. The constitutive relation needs not be symmetric. 280 281 @cvar _zeropoints_: list of tuples for indices where the content should be 282 zero. 283 @ctype _zeropoints_: list 284 @ivar rho: density 285 @ivar al: alpha angle. 286 @ivar be: beta angle. 287 @ivar ga: gamma angle. 288 @ivar origstiff: stiffness matrix in the crystal coordinate. 289 @ivar stiff: stiffness matrix in the transformed global coordinate. 290 """ 291 292 __metaclass__ = MaterialMeta 293 294 _zeropoints_ = [] 295 296 from numpy import array 297 K = array([ [ 298 [1, 0, 0, 0, 0, 0], 299 [0, 0, 0, 0, 0, 1], 300 [0, 0, 0, 0, 1, 0], 301 ], [ 302 [0, 0, 0, 0, 0, 1], 303 [0, 1, 0, 0, 0, 0], 304 [0, 0, 0, 1, 0, 0], 305 ], [ 306 [0, 0, 0, 0, 1, 0], 307 [0, 0, 0, 1, 0, 0], 308 [0, 0, 1, 0, 0, 0], 309 ], ], dtype='float64') 310 del array 311

312 - def __init__(self, *args, **kw):

313 from numpy import empty, dot 314 self.rho = kw.pop('rho') 315 self.al = kw.pop('al') 316 self.be = kw.pop('be') 317 self.ga = kw.pop('ga') 318 # set stiffness matrix. 319 origstiff = empty((6,6), dtype='float64') 320 origstiff.fill(0.0) 321 for key in kw.keys(): # becaues I pop out the key. 322 if len(key) == 4 and key[:2] == 'co': 323 try: 324 i = int(key[2])-1 325 j = int(key[3])-1 326 except: 327 continue 328 assert i < origstiff.shape[0] 329 assert j < origstiff.shape[1] 330 val = kw.pop(key) 331 origstiff[i,j] = val 332 self.origstiff = origstiff 333 # check for zeros. 334 self._check_origstiffzero(self.origstiff) 335 # compute the stiffness matrix in the transformed global coordinate 336 # system. 337 bondmat = self.get_bondmat() 338 self.stiff = dot(bondmat, dot(self.origstiff, bondmat.T)) 339 super(Material, self).__init__(*args, **kw)

340

341 - def __getattr__(self, key):

342 if len(key) == 4 and key[:2] == 'co': 343 i = int(key[2]) 344 j = int(key[3]) 345 if 1 <= i <= 6 and 1 <= j <= 6: 346 return self.origstiff[i-1,j-1] 347 elif len(key) == 3 and key[0] == 'c': 348 i = int(key[1]) 349 j = int(key[2]) 350 if 1 <= i <= 6 and 1 <= j <= 6: 351 return self.stiff[i-1,j-1] 352 else: 353 raise AttributeError

354

355 - def __str__(self):

356 from math import pi 357 return '[%s: al=%.2f be=%.2f ga=%.2f (deg)]' % (self.__class__.__name__, 358 self.al/(pi/180), self.be/(pi/180), self.ga/(pi/180))

359 360 @classmethod

361 - def _check_origstiffzero(cls, origstiff):

362 """ 363 Check for zero in original stiffness matrix. 364 365 @note: no assumed symmetry. 366 """ 367 for i, j in cls._zeropoints_: 368 assert origstiff[i,j] == 0.0

369

370 - def get_rotmat(self):

371 """ 372 Coordinate transformation matrix for three successive rotations through 373 the Euler angles. 374 375 @return: the transformation matrix. 376 @rtype: numpy.ndarray 377 """ 378 from numpy import array, cos, sin, dot 379 al = self.al; be = self.be; ga = self.ga 380 almat = array([ 381 [cos(al), sin(al), 0], 382 [-sin(al), cos(al), 0], 383 [0, 0, 1], 384 ], dtype='float64') 385 bemat = array([ 386 [1, 0, 0], 387 [0, cos(be), sin(be)], 388 [0, -sin(be), cos(be)], 389 ], dtype='float64') 390 gamat = array([ 391 [cos(ga), sin(ga), 0], 392 [-sin(ga), cos(ga), 0], 393 [0, 0, 1], 394 ], dtype='float64') 395 return dot(gamat, dot(bemat, almat))

396

397 - def get_bondmat(self):

398 """ 399 The Bond's matrix M as a shorthand of coordinate transformation for the 400 6-component stress vector. 401 402 @return: the Bond's matrix. 403 @rtype: numpy.ndarray 404 """ 405 from numpy import empty 406 rotmat = self.get_rotmat() 407 bond = empty((6,6), dtype='float64') 408 # upper left. 409 bond[:3,:3] = rotmat[:,:]**2 410 # upper right. 411 bond[0,3] = 2*rotmat[0,1]*rotmat[0,2] 412 bond[0,4] = 2*rotmat[0,2]*rotmat[0,0] 413 bond[0,5] = 2*rotmat[0,0]*rotmat[0,1] 414 bond[1,3] = 2*rotmat[1,1]*rotmat[1,2] 415 bond[1,4] = 2*rotmat[1,2]*rotmat[1,0] 416 bond[1,5] = 2*rotmat[1,0]*rotmat[1,1] 417 bond[2,3] = 2*rotmat[2,1]*rotmat[2,2] 418 bond[2,4] = 2*rotmat[2,2]*rotmat[2,0] 419 bond[2,5] = 2*rotmat[2,0]*rotmat[2,1] 420 # lower left. 421 bond[3,0] = rotmat[1,0]*rotmat[2,0] 422 bond[3,1] = rotmat[1,1]*rotmat[2,1] 423 bond[3,2] = rotmat[1,2]*rotmat[2,2] 424 bond[4,0] = rotmat[2,0]*rotmat[0,0] 425 bond[4,1] = rotmat[2,1]*rotmat[0,1] 426 bond[4,2] = rotmat[2,2]*rotmat[0,2] 427 bond[5,0] = rotmat[0,0]*rotmat[1,0] 428 bond[5,1] = rotmat[0,1]*rotmat[1,1] 429 bond[5,2] = rotmat[0,2]*rotmat[1,2] 430 # lower right. 431 bond[3,3] = rotmat[1,1]*rotmat[2,2] + rotmat[1,2]*rotmat[2,1] 432 bond[3,4] = rotmat[1,0]*rotmat[2,2] + rotmat[1,2]*rotmat[2,0] 433 bond[3,5] = rotmat[1,1]*rotmat[2,0] + rotmat[1,0]*rotmat[2,1] 434 bond[4,3] = rotmat[0,1]*rotmat[2,2] + rotmat[0,2]*rotmat[2,1] 435 bond[4,4] = rotmat[0,0]*rotmat[2,2] + rotmat[0,2]*rotmat[2,0] 436 bond[4,5] = rotmat[0,1]*rotmat[2,0] + rotmat[0,0]*rotmat[2,1] 437 bond[5,3] = rotmat[0,1]*rotmat[1,2] + rotmat[0,2]*rotmat[1,1] 438 bond[5,4] = rotmat[0,0]*rotmat[1,2] + rotmat[0,2]*rotmat[1,0] 439 bond[5,5] = rotmat[0,1]*rotmat[1,0] + rotmat[0,0]*rotmat[1,1] 440 return bond

441

442 - def get_jacos(self):

443 """ 444 Obtain the Jacobian matrices for the solid. 445 446 @param K: the K matrix. 447 @type K: numpy.ndarray 448 @return: the Jacobian matrices 449 @rtype: numpy.ndarray 450 """ 451 from numpy import zeros, dot 452 rho = self.rho 453 sf = self.stiff 454 jacos = zeros((3,9,9), dtype='float64') 455 for idm in range(3): 456 K = self.K[idm] 457 jaco = jacos[idm] 458 jaco[:3,3:] = K/(-rho) # the upper right submatrix. 459 jaco[3:,:3] = -dot(sf, K.T) # the lower left submatrix. 460 return jacos

461

462 ################################################################################ 463 # Symmetry. 464 ################################################################################ 465 466 -class Triclinic(Material):

467 """ 468 The stiffness matrix has to be symmetric. 469 """ 470 _zeropoints_ = []

471 - def __init__(self, *args, **kw):

472 for key in kw.keys(): # becaues I modify the key. 473 if len(key) == 4 and key[:2] == 'co': 474 try: 475 i = int(key[2]) 476 j = int(key[3]) 477 except: 478 continue 479 symkey = 'co%d%d' % (j, i) 480 if i != j: 481 assert symkey not in kw 482 kw[symkey] = kw[key] 483 super(Triclinic, self).__init__(*args, **kw)

484 @classmethod

485 - def _check_origstiffzero(cls, origstiff):

486 for i, j in cls._zeropoints_: 487 assert origstiff[i,j] == origstiff[j,i] == 0.0

488

489 -class Monoclinic(Triclinic):

490 _zeropoints_ = [ 491 (0,3), (0,5), 492 (1,3), (1,5), 493 (2,3), (2,5), 494 (3,4), (4,5), 495 ]

496

497 -class Orthorhombic(Triclinic):

498 _zeropoints_ = [ 499 (0,3), (0,4), (0,5), 500 (1,3), (1,4), (1,5), 501 (2,3), (2,4), (2,5), 502 (3,4), (3,5), (4,5), 503 ]

504

505 -class Tetragonal(Triclinic):

506 _zeropoints_ = [ 507 (0,3), (0,4), 508 (1,3), (1,4), 509 (2,3), (2,4), (2,5), 510 (3,4), (3,5), (4,5), 511 ]

512 - def __init__(self, *args, **kw):

513 kw['co22'] = kw['co11'] 514 kw['co23'] = kw['co13'] 515 kw['co26'] = -kw.get('co16', 0.0) 516 kw['co55'] = kw['co44'] 517 super(Tetragonal, self).__init__(*args, **kw)

518

519 -class Trigonal(Triclinic):

520 _zeropoints_ = [ 521 (0,5), (1,5), 522 (2,3), (2,4), (2,5), 523 (3,4), 524 ]

525 - def __init__(self, *args, **kw):

526 kw['co15'] = -kw.get('co25', 0.0) 527 kw['co22'] = kw['co11'] 528 kw['co23'] = kw['co13'] 529 kw['co24'] = -kw.get('co14', 0.0) 530 kw['co46'] = kw.get('co25', 0.0) 531 kw['co55'] = kw['co44'] 532 kw['co56'] = kw.get('co14', 0.0) 533 kw['co66'] = (kw['co11'] - kw['co12'])/2 534 super(Trigonal, self).__init__(*args, **kw)

535

536 -class Hexagonal(Trigonal):

537 _zeropoints_ = [ 538 (0,3), (0,4), (0,5), 539 (1,3), (1,4), (1,5), 540 (2,3), (2,4), (2,5), 541 (3,4), (3,5), (4,5), 542 ]

543

544 -class Cubic(Triclinic):

545 _zeropoints_ = [ 546 (0,3), (0,4), (0,5), 547 (1,3), (1,4), (1,5), 548 (2,3), (2,4), (2,5), 549 (3,4), (3,5), (4,5), 550 ]

551 - def __init__(self, *args, **kw):

552 kw['co13'] = kw['co12'] 553 kw['co22'] = kw['co11'] 554 kw['co23'] = kw['co12'] 555 kw['co33'] = kw['co11'] 556 kw['co55'] = kw['co44'] 557 kw['co66'] = kw['co44'] 558 super(Cubic, self).__init__(*args, **kw)

559

560 -class Isotropic(Triclinic):

561 _zeropoints_ = [ 562 (0,3), (0,4), (0,5), 563 (1,3), (1,4), (1,5), 564 (2,3), (2,4), (2,5), 565 (3,4), (3,5), (4,5), 566 ]

567 - def __init__(self, *args, **kw):

568 kw['co12'] = kw['co11']-2*kw['co44'] 569 kw['co13'] = kw['co11']-2*kw['co44'] 570 kw['co22'] = kw['co11'] 571 kw['co23'] = kw['co11']-2*kw['co44'] 572 kw['co33'] = kw['co11'] 573 kw['co55'] = kw['co44'] 574 kw['co66'] = kw['co44'] 575 super(Isotropic, self).__init__(*args, **kw)

576

577 ################################################################################ 578 # Material properties. 579 ################################################################################ 580 581 -class GaAs(Cubic):

582 - def __init__(self, *args, **kw):

583 kw.setdefault('rho', 5307.0) 584 kw.setdefault('co11', 11.88e10) 585 kw.setdefault('co12', 5.38e10) 586 kw.setdefault('co44', 5.94e10) 587 super(GaAs, self).__init__(*args, **kw)

588

589 -class ZnO(Hexagonal):

590 - def __init__(self, *args, **kw):

591 kw.setdefault('rho', 5680.0) 592 kw.setdefault('co11', 20.97e10) 593 kw.setdefault('co12', 12.11e10) 594 kw.setdefault('co13', 10.51e10) 595 kw.setdefault('co33', 21.09e10) 596 kw.setdefault('co44', 4.247e10) 597 super(ZnO, self).__init__(*args, **kw)

598

599 -class CdS(Hexagonal):

600 - def __init__(self, *args, **kw):

601 kw.setdefault('rho', 4820.0) 602 kw.setdefault('co11', 9.07e10) 603 kw.setdefault('co12', 5.81e10) 604 kw.setdefault('co13', 5.1e10) 605 kw.setdefault('co33', 9.38e10) 606 kw.setdefault('co44', 1.504e10) 607 super(CdS, self).__init__(*args, **kw)

608

609 -class Zinc(Hexagonal):

610 - def __init__(self, *args, **kw):

611 kw.setdefault('rho', 7.1*1.e-3/(1.e-2**3)) 612 kw.setdefault('co11', 14.3e11*1.e-5/(1.e-2**2)) 613 kw.setdefault('co12', 1.7e11*1.e-5/(1.e-2**2)) 614 kw.setdefault('co13', 3.3e11*1.e-5/(1.e-2**2)) 615 kw.setdefault('co33', 5.0e11*1.e-5/(1.e-2**2)) 616 kw.setdefault('co44', 4.0e11*1.e-5/(1.e-2**2)) 617 super(Zinc, self).__init__(*args, **kw)

618

619 -class Beryl(Hexagonal):

620 - def __init__(self, *args, **kw):

621 kw.setdefault('rho', 2.7*1.e-3/(1.e-2**3)) 622 kw.setdefault('co11', 26.94e11*1.e-5/(1.e-2**2)) 623 kw.setdefault('co12', 9.61e11*1.e-5/(1.e-2**2)) 624 kw.setdefault('co13', 6.61e11*1.e-5/(1.e-2**2)) 625 kw.setdefault('co33', 23.63e11*1.e-5/(1.e-2**2)) 626 kw.setdefault('co44', 6.53e11*1.e-5/(1.e-2**2)) 627 super(Beryl, self).__init__(*args, **kw)

628

629 -class Albite(Triclinic):

630 - def __init__(self, *args, **kw):

631 #kw.setdefault('rho', ) 632 kw.setdefault('co11', 69.9e9) 633 kw.setdefault('co22', 183.5e9) 634 kw.setdefault('co33', 179.5e9) 635 kw.setdefault('co44', 24.9e9) 636 kw.setdefault('co55', 26.8e9) 637 kw.setdefault('co66', 33.5e9) 638 kw.setdefault('co12', 34.0e9) 639 kw.setdefault('co13', 30.8e9) 640 kw.setdefault('co14', 5.1e9) 641 kw.setdefault('co15', -2.4e9) 642 kw.setdefault('co16', -0.9e9) 643 kw.setdefault('co23', 5.5e9) 644 kw.setdefault('co24', -3.9e9) 645 kw.setdefault('co25', -7.7e9) 646 kw.setdefault('co26', -5.8e9) 647 kw.setdefault('co34', -8.7e9) 648 kw.setdefault('co35', 7.1e9) 649 kw.setdefault('co36', -9.8e9) 650 kw.setdefault('co45', -2.4e9) 651 kw.setdefault('co46', -7.2e9) 652 kw.setdefault('co56', 0.5e9) 653 super(Albite, self).__init__(*args, **kw)

654

655 -class Acmite(Monoclinic):

656 - def __init__(self, *args, **kw):

657 kw.setdefault('rho', 3.5e3) 658 kw.setdefault('co11', 185.8e9) 659 kw.setdefault('co22', 181.3e9) 660 kw.setdefault('co33', 234.4e9) 661 kw.setdefault('co44', 62.9e9) 662 kw.setdefault('co55', 51.0e9) 663 kw.setdefault('co66', 47.4e9) 664 kw.setdefault('co12', 68.5e9) 665 kw.setdefault('co13', 70.7e9) 666 kw.setdefault('co15', 9.8e9) 667 kw.setdefault('co23', 62.9e9) 668 kw.setdefault('co25', 9.4e9) 669 kw.setdefault('co35', 21.4e9) 670 kw.setdefault('co46', 7.7e9) 671 super(Acmite, self).__init__(*args, **kw)

672

673 -class AlphaUranium(Orthorhombic):

674 - def __init__(self, *args, **kw):

675 #kw.setdefault('rho', ) 676 kw.setdefault('rho', 8.2e3) # a false value. 677 kw.setdefault('co11', 215.e9) 678 kw.setdefault('co22', 199.e9) 679 kw.setdefault('co33', 267.e9) 680 kw.setdefault('co44', 124.e9) 681 kw.setdefault('co55', 73.e9) 682 kw.setdefault('co66', 74.e9) 683 kw.setdefault('co12', 46.e9) 684 kw.setdefault('co13', 22.e9) 685 kw.setdefault('co23', 107.e9) 686 super(AlphaUranium, self).__init__(*args, **kw)

687

688 -class BariumTitanate(Tetragonal):

689 - def __init__(self, *args, **kw):

690 kw.setdefault('rho', 6.2e3) 691 kw.setdefault('co11', 275.0e9) 692 kw.setdefault('co33', 165.0e9) 693 kw.setdefault('co44', 54.3e9) 694 kw.setdefault('co66', 113.0e9) 695 kw.setdefault('co12', 179.0e9) 696 kw.setdefault('co13', 151.0e9) 697 super(BariumTitanate, self).__init__(*args, **kw)

698

699 -class AlphaQuartz(Trigonal):

700 - def __init__(self, *args, **kw):

701 kw.setdefault('rho', 2.651e3) 702 kw.setdefault('co11', 87.6e9) 703 kw.setdefault('co33', 106.8e9) 704 kw.setdefault('co44', 57.2e9) 705 kw.setdefault('co12', 6.1e9) 706 kw.setdefault('co13', 13.3e9) 707 kw.setdefault('co14', 17.3e9) 708 super(AlphaQuartz, self).__init__(*args, **kw)

709

710 -class RickerSample(Isotropic):

711 - def __init__(self, *args, **kw):

712 kw.setdefault('rho', 2200.e0) 713 kw.setdefault('co11', 3200.e0**2*2200.e0) 714 kw.setdefault('co44', 1847.5e0**2*2200.e0) 715 super(RickerSample, self).__init__(*args, **kw)

716 -class RickerSampleLight(Isotropic):

717 - def __init__(self, *args, **kw):

718 scale = 1.e-3 719 kw.setdefault('rho', 2200.e0*scale) 720 kw.setdefault('co11', 3200.e0**2*2200.e0*scale) 721 kw.setdefault('co44', 1847.5e0**2*2200.e0*scale) 722 super(RickerSampleLight, self).__init__(*args, **kw)

723

Source Code for Module solvcon.kerpak.vslin