Subject: [FFC-dev] square root of a sum

From: Jake Ostien <[email-masked]>
Date: Tue, 10 Jul 2007 14:47:54 -0400

Hi,

Given the following form

    DG09 = VectorElement("Discontinuous Lagrange", "tetrahedron", 0, 9)

    V  = TestFunction(DG09)
    Hp = TrialFunction(DG09)

    chi = 0.1

    # definitions
    def Mat(A):
        return ( [A[0], A[3], A[8]],
                 [A[6], A[1], A[4]],
                 [A[5], A[7], A[2]] )

    def sym(A):
        return 0.5*(A + transp(A))

    def skw(A):
        return 0.5*(A - transp(A))

    def denom(A):
         return 1/sqrt(dot(2*sym(A),2*sym(A)) + chi*dot(2*skw(A),2*skw(A)))

    a = denom(Mat(Hp))*dot(Mat(V),Mat(Hp))*dx

I get an error stating:

        a = denom(Mat(Hp))*dot(Mat(V),Mat(Hp))*dx
      File "invsqrtsum.py", line 30, in denom
        return 1/sqrt(dot(2*sym(A),2*sym(A)) + chi*dot(2*skw(A),2*skw(A)))
      File
    "/home/jtostie/builds/lib/python/ffc/compiler/language/operators.py",
    line 231, in sqrt
        return Form(v).sqrt()
      File
    "/home/jtostie/builds/lib/python/ffc/compiler/language/algebra.py",
    line 628, in sqrt
        raise FormError, (self, "Cannot take square root of sum.")

I can see that there is a check for len(self.monomials) > 1.  Is there a 
fundamental reason why this can't be relaxed given a (guaranteed) 
positive sum?  Is there another way around this restriction? 

Thanks,
Jake

---

From: Anders Logg <[email-masked]>
Date: Tue, 10 Jul 2007 21:49:46 +0200

Yes, it's intentional. You can only take the square root of a Function, 
not a sum of two or more Functions. This is because the square root is 
applied to the *coefficients* of the basis function expansion of the 
Function.

A couple of suggestions:

1. Do you really want to integrate the square root of something? I'm not 
sure what you are integrating, but if you are computing something like 
an L2 norm then the square root is applied after integration.

2. First compute the projection of the expression you want to take the 
square root of, say

    (v, g) = (v, dot(f, f)) for all v

Then integrate the square root of 1/g:

    sqrt(1/g)*dx

/Anders

Jake Ostien wrote:
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev