Finally have a good model for material axes.

Implemented the three basic material axis types available in FEBio, and created
tests for them.  Now to create materials that use these axes.

febabel/materials.py

@@ -7,6 +7,7 @@ class Material(object):
 
     def _store(self, params):
         del params['self']
+        self.parameters = params.copy()
         for k,v in params.iteritems():
             setattr(self, k, v)
 
@@ -36,3 +37,91 @@ class Rigid(Material, Constrainable):
     def __init__(self, COM=None):
         Constrainable.__init__(self)
         self._store(locals())
+
+
+class TransIsoElastic(Material):
+    "Adds fibers to a (hyper)elastic base material."
+    def __init__(self, c3, c4, c5, lam_max, fiber_dist, base):
+        self._store(locals())
+
+
+
+class AxisOrientation(object):
+    """Base class for special material axis orientation objects.
+
+    Axes are defined by callable objects.  The call should take a single
+    Element object, and return a set of three mutually-orthogonal unit vectors
+    describing the orientation of the material axes within the given Element.
+
+    Arbitrary callable objects can be used for axis orientation; special
+    objects are used to identify special cases that solvers may have built-in."""
+
+    def __init__(self, pos1, pos2):
+        p = iter(pos1)
+        self.pos1 = [p.next(), p.next(), p.next()]
+        p = iter(pos2)
+        self.pos2 = [p.next(), p.next(), p.next()]
+
+    @staticmethod
+    def _normalize(v1, v2):
+        """Takes two vectors and returns three normalized unit vectors
+        describing the corresponding axes.
+        e1 = v1 / abs(v1), e2 = e3 x e1, e3 = (v1 x v2) / abs(v1 x v2)"""
+        # All to avoid a numpy dependency...
+        from math import sqrt
+        absv1 = sqrt( sum(i*i for i in v1) )
+        e1 = [i/absv1 for i in v1]
+
+        v1xv2 = [
+            v1[1]*v2[2] - v1[2]*v2[1],
+            v1[2]*v2[0] - v1[0]*v2[2],
+            v1[0]*v2[1] - v1[1]*v2[0],
+        ]
+        absv1xv2 = sqrt( sum(i*i for i in v1xv2) )
+        e3 = [i/absv1xv2 for i in v1xv2]
+
+        e2 = [
+            e3[1]*e1[2] - e3[2]*e1[1],
+            e3[2]*e1[0] - e3[0]*e1[2],
+            e3[0]*e1[1] - e3[1]*e1[0],
+        ]
+
+        return (e1, e2, e3)
+
+
+class VectorOrientation(AxisOrientation):
+    """Gives constant material axes throughout the entire domain.
+    pos1 gives the primary axis, while pos2 indicates the direction of the
+    secondary axis in the form:
+    e1 = pos1 / abs(pos1), e2 = e3 x e1, e3 = (pos1 x pos2) / abs(pos1 x pos2)"""
+    def __call__(self, element):
+        return self._normalize(self.pos1, self.pos2)
+
+class SphericalOrientation(AxisOrientation):
+    """Gives primary material axes radiating outward from a central point.
+    pos1 is the location of the central point.  pos2 is a vector indicating the
+    direction of the secondary axis."""
+    def __call__(self, element):
+        v1 = [ e - p for e,p in zip(element.get_vertex_avg(), self.pos1) ]
+        return self._normalize(v1, self.pos2)
+
+class NodalOrientation(AxisOrientation):
+    """Gives axes based on the positions of an Element's Nodes.
+    edge1 and edge2 are each a pair of integers corresponding to 0-indexed node
+    numbers.  Each element is given a primary axis in the direction of edge1,
+    and a secondary axis based on the direction of edge2."""
+    def __init__(self, edge1, edge2):
+        # Where edge1 and edge2 are each a pair of integers corresponding to
+        # element node indices.
+        p = iter(edge1)
+        self.edge1 = [p.next(), p.next()]
+        p = iter(edge2)
+        self.edge2 = [p.next(), p.next()]
+        # NOTE: Orthotropic materials in FEBio expect edge1[0] == edge2[0].
+
+    def __call__(self, element):
+        v1 = [ element[ self.edge1[1] ][i] - element[ self.edge1[0] ][i]
+            for i in xrange(3) ]
+        v2 = [ element[ self.edge2[1] ][i] - element[ self.edge2[0] ][i]
+            for i in xrange(3) ]
+        return self._normalize(v1, v2)

test/test_materials.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python2
+import unittest
+from math import sqrt
+
+import sys, os
+# For Python 3, use the translated version of the library.
+# For Python 2, find the library one directory up.
+if sys.version < '3':
+    sys.path.append(os.path.dirname(sys.path[0]))
+import febabel as f
+m = f.materials
+
+
+class TestMaterials(unittest.TestCase):
+
+
+    def test_parameter_storing(self):
+        l = m.LinearIsotropic(12.5, 15)
+        self.assertEqual(l.E, 12.5)
+        self.assertEqual(l.v, 15)
+        for k,v in l.parameters.items():
+            self.assertEqual(v, getattr(l, k))
+        mr = m.MooneyRivlin(1,4,9)
+        self.assertEqual(mr.c1, 1)
+        self.assertEqual(mr.c2, 4)
+        self.assertEqual(mr.k, 9)
+        for k,v in mr.parameters.items():
+            self.assertEqual(v, getattr(mr, k))
+
+
+
+class TestAxes(unittest.TestCase):
+    def setUp(self):
+        h = f.geometry.Hex8
+        n = f.geometry.Node
+        self.elements = [
+            h( [n((1,1,1)), n((2,1,1)), n((2,2,1)), n((1,2,1)),
+                n((1,1,2)), n((2,1,2)), n((2,2,2)), n((1,2,2))] ),
+            h( [n((0,4,3))] * 8 )
+        ]
+
+    def assertArrayAlmostEqual(self, a, b):
+        '''Convenience function for testing arrays of inexact floats.'''
+        self.assertArrayEqual(a, b, almost=True)
+
+    def assertArrayEqual(self, a, b, n=[], almost=False):
+        '''Convenience function for testing arrays exactly.'''
+        if hasattr(a, '__getitem__') and hasattr(b, '__getitem__'):
+            n = n[:]; n.append(0)
+            for m,(i,j) in enumerate(zip(a,b)):
+                n[-1] = m
+                self.assertArrayEqual(i,j, n, almost)
+        else:
+            f = self.assertEqual if not almost else self.assertAlmostEqual
+            f(a,b, msg='Elements differ at %s' % n)
+
+
+    def test_normalization(self):
+        self.assertArrayAlmostEqual( m.AxisOrientation._normalize([1,0,0],[0,1,0]),
+            ([1,0,0],[0,1,0],[0,0,1]) )
+        self.assertArrayAlmostEqual( m.AxisOrientation._normalize([1,0,0],[0.234,1,0]),
+            ([1,0,0],[0,1,0],[0,0,1]) )
+        self.assertArrayAlmostEqual( m.AxisOrientation._normalize([3,4,0],[3,4,1]),
+            ([0.6, 0.8, 0], [0,0,1], [0.8, -0.6, 0]) )
+
+        vs = m.AxisOrientation._normalize([2.312211, 745.3442, 3.505], [543, 1.1, 1.2])
+        for v in vs:
+            length = sqrt( sum( i*i for i in v) )
+            self.assertAlmostEqual(length, 1,
+                msg='Basis vector is not unit length!')
+        for va in vs:
+            for vb in vs:
+                if not (va is vb):
+                    dot = sum( i*j for i,j in zip(va, vb) )
+                    self.assertAlmostEqual(dot, 0,
+                        msg='Basis vectors are not orthogonal!')
+
+
+    def test_vector_orientation(self):
+        self.assertEqual( m.VectorOrientation([1,0,0],[0,1,0])(self.elements[0]),
+            ([1,0,0],[0,1,0],[0,0,1]) )
+        self.assertEqual( m.VectorOrientation([1,0,0],[0.234,1,0])(self.elements[0]),
+            ([1,0,0],[0,1,0],[0,0,1]) )
+        self.assertEqual( m.VectorOrientation([3,4,0],[3,4,1])(self.elements[0]),
+            ([0.6, 0.8, 0], [0,0,1], [0.8, -0.6, 0]) )
+
+
+    def test_spherical_orientation(self):
+        self.assertEqual( m.SphericalOrientation([0,0,0],[0,0,1])(self.elements[1]),
+            ([0, 0.8, 0.6], [0, -0.6, 0.8], [1,0,0]) )
+        # Imprecision in the floats makes standard equality fail here.
+        self.assertArrayAlmostEqual(
+            m.SphericalOrientation([0,0,0],[0,0,1])(self.elements[0]),
+            ([1/sqrt(3), 1/sqrt(3), 1/sqrt(3)], [-1/sqrt(6), -1/sqrt(6), sqrt(2.0/3)],
+                [1/sqrt(2), -1/sqrt(2), 0]) )
+
+
+    def test_nodal_orientation(self):
+        self.assertEqual( m.NodalOrientation((0,1),(0,3))(self.elements[0]),
+            ([1,0,0],[0,1,0],[0,0,1]) )
+        self.assertEqual( m.NodalOrientation((0,1),(0,2))(self.elements[0]),
+            ([1,0,0],[0,1,0],[0,0,1]) )
+        self.assertEqual( m.NodalOrientation((0,1),(1,2))(self.elements[0]),
+            ([1,0,0],[0,1,0],[0,0,1]) )
+        self.assertEqual( m.NodalOrientation((0,1),(0,4))(self.elements[0]),
+            ([1,0,0],[0,0,1],[0,-1,0]) )
+
+
+
+
+if __name__ == '__main__':
+    unittest.main()