algebra.py

finite element library for mathematic majored research

    """

    def __init__(self, other = None):
        "Create Monomial."
        if other == None:
            # Create default Monomial (unity)
            self.numeric = 1.0
            self.constants = []
            self.coefficients = []
            self.transforms = []
            self.basisfunctions = []
            self.determinants = []
            self.integral = None
        elif isinstance(other, int) or isinstance(other, float):
            # Create Monomial from scalar
            self.numeric = float(other)
            self.constants = []
            self.coefficients = []
            self.transforms = []
            self.basisfunctions = []
            self.determinants = []
            self.integral = None
        elif isinstance(other, Function):
            # Create Monomial from Function
            index = Index()
            self.numeric = 1.0
            self.constants = []
            self.coefficients = [Coefficient(other, index)]
            self.transforms = []
            self.basisfunctions = [BasisFunction(other.e1, index)]
            self.determinants = []
            self.integral = None
        elif isinstance(other, BasisFunction):
            # Create Monomial from BasisFunction
            self.numeric = 1.0
            self.constants = []
            self.coefficients = []
            self.transforms = []
            self.basisfunctions = [BasisFunction(other)]
            self.determinants = []
            self.integral = None
        elif isinstance(other, Monomial):
            # Create Monomial from Monomial (copy constructor)
            self.numeric = float(other.numeric)
            self.constants = listcopy(other.constants)
            self.coefficients = listcopy(other.coefficients)
            self.transforms = listcopy(other.transforms)
            self.basisfunctions = listcopy(other.basisfunctions)
            self.determinants = listcopy(other.determinants)
            self.integral = other.integral
        else:
            raise FormError, (other, "Unable to create Monomial from given expression.")

        return
    
    def __mul__(self, other):
        "Operator: Monomial * Element"
        if isinstance(other, Integral):
            if not self.integral == None:
                raise FormError, (self, "Integrand can only be integrated once.")
            w = Monomial(self)
            w.integral = Integral(other)
            return w
        elif isinstance(other, Form):
            return Form(self) * Form(other)
        else:
            # Create two copies
            w0 = Monomial(self)
            w1 = Monomial(other)
            # Check ranks
            if not w0.value_rank() == w1.value_rank() == 0:
                raise FormError, (self, "Operands for monomial must be scalar.")
            # Reassign all complete Indices to avoid collisions
            reassign_complete(w0, Index.SECONDARY)
            reassign_complete(w0, Index.AUXILIARY)
            reassign_complete(w1, Index.SECONDARY)
            reassign_complete(w1, Index.AUXILIARY)
            # Compute monomial
            w = Monomial()
            w.numeric = float(w0.numeric * w1.numeric)
            w.constants = listcopy(w0.constants + w1.constants)
            w.coefficients = listcopy(w0.coefficients + w1.coefficients)
            w.transforms = listcopy(w0.transforms + w1.transforms)
            w.basisfunctions = listcopy(w0.basisfunctions + w1.basisfunctions)
            w.determinants = listcopy(w0.determinants + w1.determinants)
            if w0.integral and w1.integral:
                raise FormError, (self, "Integrand can only be integrated once.")
            elif w0.integral:
                w.integral = Integral(w0.integral)
            elif w1.integral:
                w.integral = Integral(w1.integral)
            else:
                w.integral = None;

            return w

    def __neg__(self):
        "Operator: -Monomial"
        w = Monomial(self)
        w.numeric = -w.numeric
        return w

    def __invert__(self):
        "Operator: ~Monomial"
        w = Monomial(self)
        w.numeric = 1.0/w.numeric
        for i in range(len(w.coefficients)):
            w.coefficients[i].ops = [Operators.INVERSE] + w.coefficients[i].ops
        return w

    def modulus(self):
        "Take absolute value of monomial"
        w = Monomial(self)
        w.numeric = abs(w.numeric)
        for i in range(len(w.coefficients)):
            w.coefficients[i].ops = [Operators.MODULUS] + w.coefficients[i].ops
        return w

    def sqrt(self):
        "Take square root of monomial"
        w = Monomial(self)
        w.numeric = math.sqrt(w.numeric)
        for i in range(len(w.coefficients)):
            w.coefficients[i].ops = [Operators.SQRT] + w.coefficients[i].ops
        return w

    def  __getitem__(self, component):
        "Operator: Monomial[component], pick given component."
        # Always scalar if monomial of more than one basis function
        if not len(self.basisfunctions) == 1:
            raise FormError, (self, "Cannot pick component of scalar expression.")
        # Otherwise, return component of first and only BasisFunction
        p = Monomial(self)
        p.basisfunctions = [] 
        w = Monomial(self.basisfunctions[0][component]) 
        return w*p

    def __len__(self):
        "Operator: len(Monomial)"
        # Always scalar if monomial of more than one basis function
        if not len(self.basisfunctions) == 1:
            raise FormError, (self, "Vector length of scalar expression is undefined.")
        # Otherwise, return length of first and only BasisFunction
        return len(self.basisfunctions[0])

    def __call__(self, r):
        v = Monomial(self)
        v.basisfunctions = ([w(r) for w in v.basisfunctions])
        # Same restriction for all basis functions so pick first
        restriction = v.basisfunctions[0].restriction
        for i in range(len(v.transforms)):
            v.transforms[i].restriction = restriction
        for i in range(len(v.determinants)):
            v.determinants[i].restriction = restriction
        return v

    def __repr__(self):
        "Print nicely formatted representation of Monomial."
        if not (self.coefficients or self.transforms or self.basisfunctions):
            return str(self.numeric)
        if self.numeric == -1.0:
            s = "-"
        elif not self.numeric == 1.0:
            s = str(self.numeric)
        else:
            s = ""
        d = "".join([d.__repr__() for d in self.determinants])
        c = "".join([w.__repr__() for w in self.constants])
        w = "".join([w.__repr__() for w in self.coefficients])
        t = "".join([t.__repr__() for t in self.transforms])
        v = "*".join([v.__repr__() for v in self.basisfunctions])
        if not self.integral == None:
            i = "*" + str(self.integral)
        else:
            i = ""
        return s + c + d + w + t + " | " + v + i

    def dx(self, index = None):
        "Operator: (d/dx)Monomial in given coordinate direction."
        w = Form()
        for i in range(len(self.basisfunctions)):
            p = Monomial(self)
            p.basisfunctions = []
            for j in range(len(self.basisfunctions)):
                if i == j:
                    p = p * self.basisfunctions[i].dx(index)
                else:
                    p = p * self.basisfunctions[j]
            w = w + p
        return w

    def value_rank(self):
        "Return value rank of Monomial."
        if not self.basisfunctions:
            return 0
        if len(self.basisfunctions) > 1:
            for v in self.basisfunctions:
                if not v.value_rank() == 0:
                    raise FormError, (self, "Illegal rank for BasisFunction of Monomial (non-scalar).")
        return self.basisfunctions[0].value_rank()
            
class Form(Element):
    """A Form represents a sum of Monomials. Each Monomial will be
    compiled separately, since different Monomials are probably of
    different rank.

    Attributes:

        monomials - a list of Monomials (terms) in the Form
    """
    
    def __init__(self, other = None):
        "Create Form."
        if other == None:
            # Create default Form (zero)
            self.monomials = []
        elif isinstance(other, int) or isinstance(other, float):
            # Create Form from float
            self.monomials = [Monomial(other)]
        elif isinstance(other, BasisFunction):
            # Create Form from BasisFunction
            self.monomials = [Monomial(other)]
        elif isinstance(other, Function):
            # Create Form from Function
            self.monomials = [Monomial(other)]
        elif isinstance(other, Monomial):
            # Create Form from Monomial
            self.monomials = [Monomial(other)]
        elif isinstance(other, Form):
            # Create Form from Form (copy constructor)
            self.monomials = listcopy(other.monomials)
        else:
            raise FormError, (other, "Unable to create Form from given expression.")
        return

    def __add__(self, other):
        "Operator: Form + Element"
        w0 = Form(self)
        w1 = Form(other)
        if not w0.value_rank() == w1.value_rank():
            raise FormError, (self, "Operands for addition have non-matching ranks.")
        w = Form()
        w.monomials = w0.monomials + w1.monomials
        return w
    
    def __mul__(self, other):
        "Operator: Form * Element"
        if isinstance(other, Form):
            w = Form()
            w.monomials = [p*q for p in self.monomials for q in other.monomials]
            return w
        else:
            w = Form()
            w.monomials = [p*other for p in self.monomials]
            return w

    def __neg__(self):
        "Operator: -Form"
        w = Form()
        w.monomials = [-p for p in self.monomials]
        return w

    def __invert__(self):
        "Operator: ~Form"
        if len(self.monomials) > 1:
            raise FormError, (self, "Cannot take inverse of sum.")
        w = Form()
        w.monomials = [~p for p in self.monomials]
        return w

    def modulus(self):
        "Take absolute value of form"
        if len(self.monomials) > 1:
            raise FormError, (self, "Cannot take absolute value of sum.")
        w = Form()
        w.monomials = [p.modulus() for p in self.monomials]
        return w

    def sqrt(self):
        "Take square root of form"
        if len(self.monomials) > 1:
            raise FormError, (self, "Cannot take square root of sum.")
        w = Form()
        w.monomials = [p.sqrt() for p in self.monomials]
        return w

    def  __getitem__(self, component):
        "Operator: indexing, pick given component."
        w = Form()
        w.monomials = [p[component] for p in self.monomials]
        return w

    def __len__(self):
        "Operator: len(Form)"
        # Check that all terms have the same length
        for i in range(len(self.monomials) - 1):
            if not len(self.monomials[i]) == len(self.monomials[i + 1]):
                raise FormError, (self, "Terms have different vector length.")
        # Return length of first term
        return len(self.monomials[0])

    def __call__(self, r):
        v = Form(self)
        v.monomials = ([w(r) for w in v.monomials])
        return v

    def __repr__(self):
        "Print nicely formatted representation of Form."
        return " + ".join([p.__repr__() for p in self.monomials])

    def dx(self, index = None):
        "Operator: (d/dx)Form in given coordinate direction."
        w = Form()
        for p in self.monomials:
            w = w + p.dx(index)
        return w

    def value_rank(self):
        "Return value rank of Form."
        if not self.monomials:
            return 0
        for j in range(len(self.monomials) - 1):
            if not self.monomials[j].value_rank() == self.monomials[j + 1].value_rank():
                raise FormError, (self, "Terms have different rank.")
        return self.monomials[0].value_rank()
  
class TestFunction(BasisFunction):
    """A TestFunction is the BasisFunction with the lowest primary
    index. We simply pick an index lower than all others (-2)."""

    def __init__(self, element):
        index = Index("primary")
        index.index = -2
        BasisFunction.__init__(self, element, index)
        return

class TrialFunction(BasisFunction):
    """A TrialFunction is the BasisFunction with the next lowest primary
    index. We simply pick an index lower than almost all others (-1)."""

    def __init__(self, element):
        index = Index("primary")
        index.index = -1
        BasisFunction.__init__(self, element, index)
        return