solvcon.kerpak.elastic

1 # -*- coding: UTF-8 -*- 2 # 3 # Copyright (C) 2010 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 from solvcon.gendata import SingleAssignDict, AttributeDict 20 from solvcon.anchor import Anchor 21 from solvcon.hook import BlockHook 22 from .cese import CeseSolver 23 from .cese import CeseCase 24 from .cese import CeseBC

25 26 ############################################################################### 27 # Metadata for materials. 28 ############################################################################### 29 30 -class MtrlTypeRegistry(SingleAssignDict, AttributeDict):

31 """ 32 BC type registry class, and its instance holds BC type classes, which can 33 be indexed by BC type name and BC type number. 34 35 In current design, there should exist only one registry singleton in 36 package. 37 38 BC classes in registry should not be altered, in any circumstances. 39 """

40 - def register(self, bctype):

41 name = bctype.__name__ 42 self[name] = bctype 43 return bctype

44 mltregy = MtrlTypeRegistry() # registry singleton.

45 46 -class MaterialMeta(type):

47 """ 48 Meta class for material class. 49 """

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

51 newcls = super(MaterialMeta, cls).__new__(cls, name, bases, namespace) 52 # register. 53 mltregy.register(newcls) 54 return newcls

55

56 ############################################################################### 57 # Solver. 58 ############################################################################### 59 60 -class ElasticSolver(CeseSolver):

61 """ 62 Basic elastic solver. 63 64 @ivar cfldt: the time_increment for CFL calculation at boundcond. 65 @itype cfldt: float 66 @ivar cflmax: the maximum CFL number. 67 @itype cflmax: float 68 @ivar mtrldict: map from names to material objects. 69 @itype mtrldict: dict 70 @ivar mtrllist: list of all material objects. 71 @itype mtrllist: list 72 """ 73 from solvcon.dependency import getcdll 74 __clib_elastic = { 75 2: getcdll('elastic2d'), 76 3: getcdll('elastic3d'), 77 } 78 del getcdll 79 @property

80 - def _clib_elastic(self):

81 return self.__clib_elastic[self.ndim]

82 _gdlen_ = 9 * 9 * 2 83 @property

84 - def _jacofunc_(self):

85 return self._clib_elastic.calc_jaco

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

87 self.cfldt = kw.pop('cfldt', None) 88 self.cflmax = 0.0 89 self.mtrldict = kw.pop('mtrldict', {}) 90 self.mtrllist = None 91 super(ElasticSolver, self).__init__(*args, **kw)

92 @staticmethod

93 - def _build_mtrllist(grpnames, mtrldict):

94 """ 95 Build the material list out of the mapping dict. 96 97 @param grpnames: sequence of group names. 98 @type grpnames: list 99 @param mtrldict: the map from names to material objects. 100 @type mtrldict: dict 101 @return: the list of material object. 102 @rtype: Material 103 """ 104 mtrllist = list() 105 default_mtuple = mtrldict.get(None, None) 106 for grpname in grpnames: 107 try: 108 mtrl = mtrldict.get(grpname, default_mtuple) 109 except KeyError, e: 110 args = e.args[:] 111 args.append('no material named %s in mtrldict'%grpname) 112 e.args = args 113 raise 114 mtrllist.append(mtrl) 115 return mtrllist

116 - def provide(self):

117 from ctypes import byref, c_int 118 # build material list. 119 self.mtrllist = self._build_mtrllist(self.grpnames, self.mtrldict) 120 for igrp in range(len(self.grpnames)): 121 mtrl = self.mtrllist[igrp] 122 jaco = self.grpda[igrp].reshape(self.neq, self.neq, self.ndim) 123 jaco[:,:,0] = mtrl.jacox 124 jaco[:,:,1] = mtrl.jacoy 125 if self.ndim == 3: 126 jaco[:,:,2] = mtrl.jacoz 127 # pre-calculate CFL. 128 self._set_time(self.time, self.cfldt) 129 self._clib_elastic.calc_cfl( 130 byref(self.exd), c_int(0), c_int(self.ncell)) 131 self.cflmax = self.cfl.max() 132 # super method. 133 super(ElasticSolver, self).provide()

134 - def calccfl(self, worker=None):

135 if self.marchret is None: 136 self.marchret = [0.0, 0.0, 0, 0] 137 self.marchret[0] = self.cflmax 138 self.marchret[1] = self.cflmax 139 self.marchret[2] = 0 140 self.marchret[3] = 0 141 return self.marchret

142

143 ############################################################################### 144 # Case. 145 ############################################################################### 146 147 -class ElasticCase(CeseCase):

148 """ 149 Case for anisotropic elastic solids. 150 """ 151 from solvcon.domain import Domain 152 defdict = { 153 'execution.neq': 9, 154 'solver.solvertype': ElasticSolver, 155 'solver.domaintype': Domain, 156 'solver.alpha': 0, 157 'solver.cfldt': None, 158 'solver.mtrldict': dict, 159 } 160 del Domain

161 - def make_solver_keywords(self):

162 kw = super(ElasticCase, self).make_solver_keywords() 163 # setup delta t for CFL calculation. 164 cfldt = self.solver.cfldt 165 cfldt = self.execution.time_increment if cfldt is None else cfldt 166 kw['cfldt'] = cfldt 167 # setup material mapper. 168 kw['mtrldict'] = self.solver.mtrldict 169 return kw

170

171 ############################################################################### 172 # Boundary conditions. 173 ############################################################################### 174 175 -class ElasticBC(CeseBC):

176 """ 177 Basic BC class for elastic problems. 178 """ 179 typn = -10200 180 from solvcon.dependency import getcdll 181 __clib_elasticb = { 182 2: getcdll('elasticb2d'), 183 3: getcdll('elasticb3d'), 184 } 185 del getcdll 186 @property

187 - def _clib_elasticb(self):

188 return self.__clib_elasticb[self.svr.ndim]

189

190 -class ElasticTraction(ElasticBC):

191 typn = 10201 192 vnames = [ 193 'bfcsys', 'tau1', 'tau2', 'tau3', 'freq', 'phase', 194 ] 195 vdefaults = { 196 'bfcsys': 0.0, 197 'tau1': 0.0, 'tau2': 0.0, 'tau3': 0.0, 'freq': 0.0, 'phase': 0.0, 198 } 199 _ghostgeom_ = 'compress'

200 - def sol(self):

201 from solvcon.dependency import intptr 202 from ctypes import byref, c_int 203 self._clib_boundcond.bound_traction_soln( 204 byref(self.svr.exd), 205 c_int(self.facn.shape[0]), 206 self.facn.ctypes.data_as(intptr), 207 c_int(self.value.shape[1]), 208 self.value.ctypes.data_as(self.fpptr), 209 )

210 - def dsol(self):

211 from solvcon.dependency import intptr 212 from ctypes import byref, c_int 213 self._clib_boundcond.bound_traction_dsoln( 214 byref(self.svr.exd), 215 c_int(self.facn.shape[0]), 216 self.facn.ctypes.data_as(intptr), 217 )

218

219 -class ElasticTractionFree(ElasticBC):

220 typn = 10202 221 _ghostgeom_ = 'mirror'

222 - def sol(self):

223 from solvcon.dependency import intptr 224 from ctypes import byref, c_int 225 self._clib_boundcond.bound_traction_free_soln( 226 byref(self.svr.exd), 227 c_int(self.facn.shape[0]), 228 self.facn.ctypes.data_as(intptr), 229 )

230 - def dsol(self):

231 from solvcon.dependency import intptr 232 from ctypes import byref, c_int 233 self._clib_boundcond.bound_traction_free_dsoln( 234 byref(self.svr.exd), 235 c_int(self.facn.shape[0]), 236 self.facn.ctypes.data_as(intptr), 237 )

238

239 -class ElasticTractionFree2(ElasticBC):

240 typn = 10203 241 _ghostgeom_ = 'mirror'

242 - def sol(self):

243 from solvcon.dependency import intptr 244 from ctypes import byref, c_int 245 self._clib_boundcond.bound_traction_free2_soln( 246 byref(self.svr.exd), 247 c_int(self.facn.shape[0]), 248 self.facn.ctypes.data_as(intptr), 249 )

250 - def dsol(self):

251 from solvcon.dependency import intptr 252 from ctypes import byref, c_int 253 self._clib_boundcond.bound_traction_free2_dsoln( 254 byref(self.svr.exd), 255 c_int(self.facn.shape[0]), 256 self.facn.ctypes.data_as(intptr), 257 )

258

259 ################################################################################ 260 # Anchor. 261 ################################################################################ 262 263 -class ElasticOAnchor(Anchor):

264 """ 265 Calculate total energy, i.e., the summation of kinetic energy and strain 266 energy. 267 """

268 - def _calculate_physics(self):

269 from ctypes import byref 270 from numpy import empty 271 from numpy.linalg import inv 272 svr = self.svr 273 # input arrays. 274 rhos = empty(svr.ngroup, dtype=svr.fpdtype) 275 comps = empty((svr.ngroup, 6, 6), dtype=svr.fpdtype) 276 for igp in range(svr.ngroup): 277 mtrl = svr.mtrllist[igp] 278 rhos[igp] = mtrl.rho 279 comps[igp,:,:] = inv(mtrl.stiff).T 280 # output arrays. 281 svr._clib_elastic.calc_energy( 282 byref(svr.exd), 283 rhos.ctypes.data_as(svr.fpptr), 284 comps.ctypes.data_as(svr.fpptr), 285 svr.der['energy'].ctypes.data_as(svr.fpptr), 286 )

287 - def provide(self):

288 from numpy import empty 289 svr = self.svr 290 svr.der['energy'] = empty(svr.ngstcell+svr.ncell, dtype=svr.fpdtype) 291 self._calculate_physics()

292 - def postfull(self):

293 self._calculate_physics()

294

295 ################################################################################ 296 # Material definition. 297 ################################################################################ 298 299 -class Material(object):

300 """ 301 Material properties. The constitutive relation needs not be symmetric. 302 303 @cvar _zeropoints_: list of tuples for indices where the content should be 304 zero. 305 @ctype _zeropoints_: list 306 @ivar rho: density 307 @ivar al: alpha angle. 308 @ivar be: beta angle. 309 @ivar ga: gamma angle. 310 @ivar origstiff: stiffness matrix in the crystal coordinate. 311 @ivar stiff: stiffness matrix in the transformed global coordinate. 312 """ 313 314 __metaclass__ = MaterialMeta 315 316 _zeropoints_ = [] 317 318 from numpy import array 319 K = array([ 320 [ 321 [1, 0, 0, 0, 0, 0], 322 [0, 0, 0, 0, 0, 1], 323 [0, 0, 0, 0, 1, 0], 324 ], 325 [ 326 [0, 0, 0, 0, 0, 1], 327 [0, 1, 0, 0, 0, 0], 328 [0, 0, 0, 1, 0, 0], 329 ], 330 [ 331 [0, 0, 0, 0, 1, 0], 332 [0, 0, 0, 1, 0, 0], 333 [0, 0, 1, 0, 0, 0], 334 ], 335 ], dtype='float64') 336 del array 337

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

339 from numpy import empty, dot 340 self.rho = kw.pop('rho') 341 self.al = kw.pop('al') 342 self.be = kw.pop('be') 343 self.ga = kw.pop('ga') 344 # set stiffness matrix. 345 origstiff = empty((6,6), dtype='float64') 346 origstiff.fill(0.0) 347 for key in kw.keys(): # becaues I pop out the key. 348 if len(key) == 4 and key[:2] == 'co': 349 try: 350 i = int(key[2])-1 351 j = int(key[3])-1 352 except: 353 continue 354 assert i < origstiff.shape[0] 355 assert j < origstiff.shape[1] 356 val = kw.pop(key) 357 origstiff[i,j] = val 358 self.origstiff = origstiff 359 # check for zeros. 360 self._check_origstiffzero(self.origstiff) 361 # compute the stiffness matrix in the transformed global coordinate 362 # system. 363 bondmat = self.bondmat 364 self.stiff = dot(bondmat, dot(self.origstiff, bondmat.T)) 365 super(Material, self).__init__(*args, **kw)

366

367 - def __getattr__(self, key):

368 if len(key) == 4 and key[:2] == 'co': 369 i = int(key[2]) 370 j = int(key[3]) 371 if 1 <= i <= 6 and 1 <= j <= 6: 372 return self.origstiff[i-1,j-1] 373 elif len(key) == 3 and key[0] == 'c': 374 i = int(key[1]) 375 j = int(key[2]) 376 if 1 <= i <= 6 and 1 <= j <= 6: 377 return self.stiff[i-1,j-1] 378 else: 379 raise AttributeError

380

381 - def __str__(self):

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

385 386 @classmethod

387 - def _check_origstiffzero(cls, origstiff):

388 """ 389 Check for zero in original stiffness matrix. 390 391 @note: no assumed symmetry. 392 """ 393 for i, j in cls._zeropoints_: 394 assert origstiff[i,j] == 0.0

395 396 @property

397 - def rotmat(self):

398 """ 399 Coordinate transformation matrix for three successive rotations through 400 the Euler angles. 401 402 @return: the transformation matrix. 403 @rtype: numpy.ndarray 404 """ 405 from numpy import array, cos, sin, dot 406 al = self.al; be = self.be; ga = self.ga 407 almat = array([ 408 [cos(al), sin(al), 0], 409 [-sin(al), cos(al), 0], 410 [0, 0, 1], 411 ], dtype='float64') 412 bemat = array([ 413 [1, 0, 0], 414 [0, cos(be), sin(be)], 415 [0, -sin(be), cos(be)], 416 ], dtype='float64') 417 gamat = array([ 418 [cos(ga), sin(ga), 0], 419 [-sin(ga), cos(ga), 0], 420 [0, 0, 1], 421 ], dtype='float64') 422 return dot(gamat, dot(bemat, almat))

423 424 @property

425 - def bondmat(self):

426 """ 427 The Bond's matrix M as a shorthand of coordinate transformation for the 428 6-component stress vector. 429 430 @return: the Bond's matrix. 431 @rtype: numpy.ndarray 432 """ 433 from numpy import empty 434 rotmat = self.rotmat 435 bond = empty((6,6), dtype='float64') 436 # upper left. 437 bond[:3,:3] = rotmat[:,:]**2 438 # upper right. 439 bond[0,3] = 2*rotmat[0,1]*rotmat[0,2] 440 bond[0,4] = 2*rotmat[0,2]*rotmat[0,0] 441 bond[0,5] = 2*rotmat[0,0]*rotmat[0,1] 442 bond[1,3] = 2*rotmat[1,1]*rotmat[1,2] 443 bond[1,4] = 2*rotmat[1,2]*rotmat[1,0] 444 bond[1,5] = 2*rotmat[1,0]*rotmat[1,1] 445 bond[2,3] = 2*rotmat[2,1]*rotmat[2,2] 446 bond[2,4] = 2*rotmat[2,2]*rotmat[2,0] 447 bond[2,5] = 2*rotmat[2,0]*rotmat[2,1] 448 # lower left. 449 bond[3,0] = rotmat[1,0]*rotmat[2,0] 450 bond[3,1] = rotmat[1,1]*rotmat[2,1] 451 bond[3,2] = rotmat[1,2]*rotmat[2,2] 452 bond[4,0] = rotmat[2,0]*rotmat[0,0] 453 bond[4,1] = rotmat[2,1]*rotmat[0,1] 454 bond[4,2] = rotmat[2,2]*rotmat[0,2] 455 bond[5,0] = rotmat[0,0]*rotmat[1,0] 456 bond[5,1] = rotmat[0,1]*rotmat[1,1] 457 bond[5,2] = rotmat[0,2]*rotmat[1,2] 458 # lower right. 459 bond[3,3] = rotmat[1,1]*rotmat[2,2] + rotmat[1,2]*rotmat[2,1] 460 bond[3,4] = rotmat[1,0]*rotmat[2,2] + rotmat[1,2]*rotmat[2,0] 461 bond[3,5] = rotmat[1,1]*rotmat[2,0] + rotmat[1,0]*rotmat[2,1] 462 bond[4,3] = rotmat[0,1]*rotmat[2,2] + rotmat[0,2]*rotmat[2,1] 463 bond[4,4] = rotmat[0,0]*rotmat[2,2] + rotmat[0,2]*rotmat[2,0] 464 bond[4,5] = rotmat[0,1]*rotmat[2,0] + rotmat[0,0]*rotmat[2,1] 465 bond[5,3] = rotmat[0,1]*rotmat[1,2] + rotmat[0,2]*rotmat[1,1] 466 bond[5,4] = rotmat[0,0]*rotmat[1,2] + rotmat[0,2]*rotmat[1,0] 467 bond[5,5] = rotmat[0,1]*rotmat[1,0] + rotmat[0,0]*rotmat[1,1] 468 return bond

469

470 - def _jaco(self, K):

471 """ 472 @param K: the K matrix. 473 @type K: numpy.ndarray 474 @return: the Jacobian matrix in x-dir. 475 @rtype: numpy.ndarray 476 """ 477 from numpy import zeros, dot 478 rho = self.rho 479 sf = self.stiff 480 jaco = zeros((9,9), dtype='float64') 481 jaco[:3,3:] = K/(-rho) # the upper right submatrix. 482 jaco[3:,:3] = -dot(sf, K.T) # the lower left submatrix. 483 return jaco

484 @property

485 - def jacox(self):

486 """ 487 @return: the Jacobian matrix in x-dir. 488 @rtype: numpy.ndarray 489 """ 490 return self._jaco(self.K[0])

491 @property

492 - def jacoy(self):

493 """ 494 @return: the Jacobian matrix in y-dir. 495 @rtype: numpy.ndarray 496 """ 497 return self._jaco(self.K[1])

498 @property

499 - def jacoz(self):

500 """ 501 @return: the Jacobian matrix in z-dir. 502 @rtype: numpy.ndarray 503 """ 504 return self._jaco(self.K[2])

505

506 - def _eig(self, jaco):

507 from numpy.linalg import eig 508 evl, evc = eig(jaco) 509 srt = evl.argsort() 510 evl = evl[srt] 511 evc = evc[:,srt] 512 return evl, evc

513 514 @property

515 - def eigx(self):

516 return self._eig(self.jacox)

517 @property

518 - def eigy(self):

519 return self._eig(self.jacoy)

520 @property

521 - def eigz(self):

522 return self._eig(self.jacoz)

523 524 @property

525 - def vpx(self):

526 return self.eigx[0][-3:]

527 @property

528 - def vpy(self):

529 return self.eigy[0][-3:]

530 @property

531 - def vpz(self):

532 return self.eigz[0][-3:]

533

534 -class Triclinic(Material):

535 """ 536 The stiffness matrix has to be symmetric. 537 """ 538 _zeropoints_ = [] 539

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

541 for key in kw.keys(): # becaues I modify the key. 542 if len(key) == 4 and key[:2] == 'co': 543 try: 544 i = int(key[2]) 545 j = int(key[3]) 546 except: 547 continue 548 symkey = 'co%d%d' % (j, i) 549 if i != j: 550 assert symkey not in kw 551 kw[symkey] = kw[key] 552 super(Triclinic, self).__init__(*args, **kw)

553 554 @classmethod

555 - def _check_origstiffzero(cls, origstiff):

556 """ 557 Check for zero in original stiffness matrix. 558 559 @note: assumed symmetric. 560 """ 561 for i, j in cls._zeropoints_: 562 assert origstiff[i,j] == origstiff[j,i] == 0.0

563

564 -class Monoclinic(Triclinic):

565 _zeropoints_ = [ 566 (0,3), (0,5), 567 (1,3), (1,5), 568 (2,3), (2,5), 569 (3,4), (4,5), 570 ]

571

572 -class Orthorhombic(Triclinic):

573 _zeropoints_ = [ 574 (0,3), (0,4), (0,5), 575 (1,3), (1,4), (1,5), 576 (2,3), (2,4), (2,5), 577 (3,4), (3,5), (4,5), 578 ]

579

580 -class Tetragonal(Triclinic):

581 _zeropoints_ = [ 582 (0,3), (0,4), 583 (1,3), (1,4), 584 (2,3), (2,4), (2,5), 585 (3,4), (3,5), (4,5), 586 ] 587

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

589 kw['co22'] = kw['co11'] 590 kw['co23'] = kw['co13'] 591 kw['co26'] = -kw.get('co16', 0.0) 592 kw['co55'] = kw['co44'] 593 super(Tetragonal, self).__init__(*args, **kw)

594

595 -class Trigonal(Triclinic):

596 _zeropoints_ = [ 597 (0,5), (1,5), 598 (2,3), (2,4), (2,5), 599 (3,4), 600 ] 601

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

603 kw['co15'] = -kw.get('co25', 0.0) 604 kw['co22'] = kw['co11'] 605 kw['co23'] = kw['co13'] 606 kw['co24'] = -kw.get('co14', 0.0) 607 kw['co46'] = kw.get('co25', 0.0) 608 kw['co55'] = kw['co44'] 609 kw['co56'] = kw.get('co14', 0.0) 610 kw['co66'] = (kw['co11'] - kw['co12'])/2 611 super(Trigonal, self).__init__(*args, **kw)

612

613 -class Hexagonal(Trigonal):

614 _zeropoints_ = [ 615 (0,3), (0,4), (0,5), 616 (1,3), (1,4), (1,5), 617 (2,3), (2,4), (2,5), 618 (3,4), (3,5), (4,5), 619 ]

620

621 -class Cubic(Triclinic):

622 _zeropoints_ = [ 623 (0,3), (0,4), (0,5), 624 (1,3), (1,4), (1,5), 625 (2,3), (2,4), (2,5), 626 (3,4), (3,5), (4,5), 627 ] 628

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

630 kw['co13'] = kw['co12'] 631 kw['co22'] = kw['co11'] 632 kw['co23'] = kw['co12'] 633 kw['co33'] = kw['co11'] 634 kw['co55'] = kw['co44'] 635 kw['co66'] = kw['co44'] 636 super(Cubic, self).__init__(*args, **kw)

637

638 -class Isotropic(Triclinic):

639 _zeropoints_ = [ 640 (0,3), (0,4), (0,5), 641 (1,3), (1,4), (1,5), 642 (2,3), (2,4), (2,5), 643 (3,4), (3,5), (4,5), 644 ] 645

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

647 kw['co12'] = kw['co11']-2*kw['co44'] 648 kw['co13'] = kw['co11']-2*kw['co44'] 649 kw['co22'] = kw['co11'] 650 kw['co23'] = kw['co11']-2*kw['co44'] 651 kw['co33'] = kw['co11'] 652 kw['co55'] = kw['co44'] 653 kw['co66'] = kw['co44'] 654 super(Isotropic, self).__init__(*args, **kw)

655

656 -class GaAs(Cubic):

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

658 kw.setdefault('rho', 5307.0) 659 kw.setdefault('co11', 11.88e10) 660 kw.setdefault('co12', 5.38e10) 661 kw.setdefault('co44', 5.94e10) 662 super(GaAs, self).__init__(*args, **kw)

663

664 -class ZnO(Hexagonal):

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

666 kw.setdefault('rho', 5680.0) 667 kw.setdefault('co11', 20.97e10) 668 kw.setdefault('co12', 12.11e10) 669 kw.setdefault('co13', 10.51e10) 670 kw.setdefault('co33', 21.09e10) 671 kw.setdefault('co44', 4.247e10) 672 super(ZnO, self).__init__(*args, **kw)

673

674 -class CdS(Hexagonal):

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

676 kw.setdefault('rho', 4820.0) 677 kw.setdefault('co11', 9.07e10) 678 kw.setdefault('co12', 5.81e10) 679 kw.setdefault('co13', 5.1e10) 680 kw.setdefault('co33', 9.38e10) 681 kw.setdefault('co44', 1.504e10) 682 super(CdS, self).__init__(*args, **kw)

683

684 -class Zinc(Hexagonal):

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

686 kw.setdefault('rho', 7.1*1.e-3/(1.e-2**3)) 687 kw.setdefault('co11', 14.3e11*1.e-5/(1.e-2**2)) 688 kw.setdefault('co12', 1.7e11*1.e-5/(1.e-2**2)) 689 kw.setdefault('co13', 3.3e11*1.e-5/(1.e-2**2)) 690 kw.setdefault('co33', 5.0e11*1.e-5/(1.e-2**2)) 691 kw.setdefault('co44', 4.0e11*1.e-5/(1.e-2**2)) 692 super(Zinc, self).__init__(*args, **kw)

693

694 -class Beryl(Hexagonal):

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

696 kw.setdefault('rho', 2.7*1.e-3/(1.e-2**3)) 697 kw.setdefault('co11', 26.94e11*1.e-5/(1.e-2**2)) 698 kw.setdefault('co12', 9.61e11*1.e-5/(1.e-2**2)) 699 kw.setdefault('co13', 6.61e11*1.e-5/(1.e-2**2)) 700 kw.setdefault('co33', 23.63e11*1.e-5/(1.e-2**2)) 701 kw.setdefault('co44', 6.53e11*1.e-5/(1.e-2**2)) 702 super(Beryl, self).__init__(*args, **kw)

703

704 -class Albite(Triclinic):

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

706 #kw.setdefault('rho', ) 707 kw.setdefault('co11', 69.9e9) 708 kw.setdefault('co22', 183.5e9) 709 kw.setdefault('co33', 179.5e9) 710 kw.setdefault('co44', 24.9e9) 711 kw.setdefault('co55', 26.8e9) 712 kw.setdefault('co66', 33.5e9) 713 kw.setdefault('co12', 34.0e9) 714 kw.setdefault('co13', 30.8e9) 715 kw.setdefault('co14', 5.1e9) 716 kw.setdefault('co15', -2.4e9) 717 kw.setdefault('co16', -0.9e9) 718 kw.setdefault('co23', 5.5e9) 719 kw.setdefault('co24', -3.9e9) 720 kw.setdefault('co25', -7.7e9) 721 kw.setdefault('co26', -5.8e9) 722 kw.setdefault('co34', -8.7e9) 723 kw.setdefault('co35', 7.1e9) 724 kw.setdefault('co36', -9.8e9) 725 kw.setdefault('co45', -2.4e9) 726 kw.setdefault('co46', -7.2e9) 727 kw.setdefault('co56', 0.5e9) 728 super(Albite, self).__init__(*args, **kw)

729

730 -class Acmite(Monoclinic):

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

732 kw.setdefault('rho', 3.5e3) 733 kw.setdefault('co11', 185.8e9) 734 kw.setdefault('co22', 181.3e9) 735 kw.setdefault('co33', 234.4e9) 736 kw.setdefault('co44', 62.9e9) 737 kw.setdefault('co55', 51.0e9) 738 kw.setdefault('co66', 47.4e9) 739 kw.setdefault('co12', 68.5e9) 740 kw.setdefault('co13', 70.7e9) 741 kw.setdefault('co15', 9.8e9) 742 kw.setdefault('co23', 62.9e9) 743 kw.setdefault('co25', 9.4e9) 744 kw.setdefault('co35', 21.4e9) 745 kw.setdefault('co46', 7.7e9) 746 super(Acmite, self).__init__(*args, **kw)

747

748 -class AlphaUranium(Orthorhombic):

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

750 #kw.setdefault('rho', ) 751 kw.setdefault('rho', 8.2e3) # a false value. 752 kw.setdefault('co11', 215.e9) 753 kw.setdefault('co22', 199.e9) 754 kw.setdefault('co33', 267.e9) 755 kw.setdefault('co44', 124.e9) 756 kw.setdefault('co55', 73.e9) 757 kw.setdefault('co66', 74.e9) 758 kw.setdefault('co12', 46.e9) 759 kw.setdefault('co13', 22.e9) 760 kw.setdefault('co23', 107.e9) 761 super(AlphaUranium, self).__init__(*args, **kw)

762

763 -class BariumTitanate(Tetragonal):

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

765 kw.setdefault('rho', 6.2e3) 766 kw.setdefault('co11', 275.0e9) 767 kw.setdefault('co33', 165.0e9) 768 kw.setdefault('co44', 54.3e9) 769 kw.setdefault('co66', 113.0e9) 770 kw.setdefault('co12', 179.0e9) 771 kw.setdefault('co13', 151.0e9) 772 super(BariumTitanate, self).__init__(*args, **kw)

773

774 -class AlphaQuartz(Trigonal):

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

776 kw.setdefault('rho', 2.651e3) 777 kw.setdefault('co11', 87.6e9) 778 kw.setdefault('co33', 106.8e9) 779 kw.setdefault('co44', 57.2e9) 780 kw.setdefault('co12', 6.1e9) 781 kw.setdefault('co13', 13.3e9) 782 kw.setdefault('co14', 17.3e9) 783 super(AlphaQuartz, self).__init__(*args, **kw)

784

785 -class RickerSample(Isotropic):

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

787 kw.setdefault('rho', 2200.e0) 788 kw.setdefault('co11', 3200.e0**2*2200.e0) 789 kw.setdefault('co44', 1847.5e0**2*2200.e0) 790 super(RickerSample, self).__init__(*args, **kw)

791 -class RickerSampleLight(Isotropic):

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

793 scale = 1.e-3 794 kw.setdefault('rho', 2200.e0*scale) 795 kw.setdefault('co11', 3200.e0**2*2200.e0*scale) 796 kw.setdefault('co44', 1847.5e0**2*2200.e0*scale) 797 super(RickerSampleLight, self).__init__(*args, **kw)

798

Source Code for Module solvcon.kerpak.elastic