|
2 | 2 | import unittest |
3 | 3 | from util import BurnManTest |
4 | 4 | import numpy as np |
| 5 | +from sympy import Matrix |
| 6 | +from fractions import Fraction |
| 7 | +import importlib |
5 | 8 |
|
6 | 9 | from burnman import Composite |
7 | | -from burnman.minerals import SLB_2011 |
| 10 | +from burnman.minerals import SLB_2011, JH_2015 |
8 | 11 | from burnman.tools.polytope import solution_polytope_from_charge_balance |
9 | 12 | from burnman.tools.polytope import solution_polytope_from_endmember_occupancies |
10 | 13 | from burnman.tools.polytope import simplify_composite_with_composition |
| 14 | +from burnman import MaterialPolytope |
11 | 15 |
|
12 | 16 |
|
13 | 17 | class polytope(BurnManTest): |
@@ -62,6 +66,44 @@ def test_simplify_composite_and_composition(self): |
62 | 66 | self.assertEqual(strings[1], "[Mg]3[Mg][Si]") |
63 | 67 | self.assertArraysAlmostEqual([0.1, 0.9], new_gt.molar_fractions) |
64 | 68 |
|
| 69 | + def test_cddlib_versions(self): |
| 70 | + gt = JH_2015.garnet() |
| 71 | + endmember_occupancies = gt.solution_model.endmember_occupancies |
| 72 | + |
| 73 | + n_sites = sum(endmember_occupancies[0]) |
| 74 | + n_occs = endmember_occupancies.shape[1] |
| 75 | + |
| 76 | + nullspace = np.array(Matrix(endmember_occupancies).nullspace(), dtype=float) |
| 77 | + |
| 78 | + equalities = np.zeros((len(nullspace) + 1, n_occs + 1)) |
| 79 | + equalities[0, 0] = -n_sites |
| 80 | + equalities[0, 1:] = 1 |
| 81 | + if len(nullspace) > 0: |
| 82 | + try: |
| 83 | + equalities[1:, 1:] = nullspace |
| 84 | + except ValueError: |
| 85 | + equalities[1:, 1:] = nullspace[:, :, 0] |
| 86 | + |
| 87 | + pos_constraints = np.concatenate( |
| 88 | + ( |
| 89 | + np.zeros((len(equalities[0]) - 1, 1)), |
| 90 | + np.identity(len(equalities[0]) - 1), |
| 91 | + ), |
| 92 | + axis=1, |
| 93 | + ) |
| 94 | + |
| 95 | + poly = MaterialPolytope( |
| 96 | + equalities, |
| 97 | + pos_constraints, |
| 98 | + independent_endmember_occupancies=endmember_occupancies, |
| 99 | + ) |
| 100 | + |
| 101 | + try: |
| 102 | + _ = importlib.import_module("cdd.gmp") |
| 103 | + self.assertTrue(type(poly.raw_vertices[0][0]) is Fraction) |
| 104 | + except ImportError: |
| 105 | + self.assertTrue(type(poly.raw_vertices[0][0]) is float) |
| 106 | + |
65 | 107 |
|
66 | 108 | if __name__ == "__main__": |
67 | 109 | unittest.main() |
0 commit comments