ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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:
> 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
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/ffc-dev