Subject: [FFC-dev] Nonlinear stokes flow, computing the form (again)

From: Jed Brown <[email-masked]>
Date: Thu, 02 Aug 2007 18:54:30 +0200

(It seems that the inline attachments were washed, sorry for the repost.)

I am having trouble getting FFC to generate code for nonlinear stokes
flow.  I realize that the fractional powers code is just not written,
but I'm not sure why the other cases don't work.  Here is my form file:

# Compile this form with FFC: ffc -l dolfin GlenStokes.form
From numpy import square as sqr

P2 = VectorElement("Lagrange", "triangle", 2)
P1 = FiniteElement("Lagrange", "triangle", 1)
TH = P2 + P1

(v, q) = TestFunctions(TH)
(u, p) = TrialFunctions(TH)

(U, P) = Functions(TH)
g = Function(P2)

D2 = (sqr(U[0].dx(0)) + sqr(U[1].dx(1))
      + 0.5 * sqr(U[0].dx(1) + U[1].dx(0)))
dD2 = (2 * U[0].dx(0) * u[0].dx(0) + 2 * U[1].dx(1) * u[1].dx(1)
       + (U[0].dx(1) + U[1].dx(0)) * (u[0].dx(1) + u[1].dx(0)))

## The *real* Glen Stokes form: powers are not yet implemented
# n = 3
# eta = D2 ** ((1 - n) / (2 * n))
# deta = ((1 - n) / (2 * n)) * D2 ** ((1 - 3 * n) / (2 * n)) * dD2
## Only use sqrt and inverse: error `Cannot take square root of sum.'
# eta = 1 / sqrt(D2)
# deta = - 1 / (2 * sqrt(D2) * D2) * dD2
## Something simpler: error `Index does not appear exactly twice in term.'
eta = D2
deta = dD2

## Linearized form to be used with fixed point method.
a = (eta * dot(grad(v), grad(u)) - div(v) * p + q * div(u)) * dx
L = dot(v, g) * dx

## Newton's method form
# a = (eta * dot(grad(v), grad(u)) - div(v) * p + q * div(u)
#      + deta * dot(grad(v), grad(U))) * dx
# L = (dot(v, g) - (eta * dot(grad(v), grad(U))
#                   - div(v) * P + q * div(U)))
#### EOF

This generates

$ ffc -l dolfin -d1 GlenStokes.form >& GlenStokes.out
This is FFC, the FEniCS Form Compiler, version 0.4.1.
For further information, go to http://www/fenics.org/ffc/.

Preprocessing form file: GlenStokes.form --> GlenStokes.py

Phase 1: Analyzing form
-----------------------

  Checking validity of form... ok
  Reassigning form indices... done
  Checking validity of form... ok
  Simplifying form... Reassigning form indices...
  done
  Reassigning form indices... done
  Reassigning form indices... done
  Reassigning form indices... done
  Reassigning form indices... done
  Reassigning form indices... done
  Reassigning form indices... done
  Reassigning form indices... done
  Reassigning form indices... done
  done
  Checking validity of form...
Traceback (most recent call last):
  File "/home/jed/usr/bin/ffc", line 180, in ?
    sys.exit(main(sys.argv[1:]))
  File "/home/jed/usr/bin/ffc", line 107, in main
    execfile(outname, ns)
  File "GlenStokes.py", line 48, in ?
    compile([a, L, M, element], "GlenStokes", "tensor", "dolfin", {'blas': False, 'precision=': '15', 'optimize': False})
  File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/compiler.py", line 66, in compile
    form_data = __compile_forms(forms, prefix, representation, language, options)
  File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/compiler.py", line 92, in __compile_forms
    form_data = analyze_form(form)
  File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/compiler.py", line 183, in analyze_form
    form_data = analyze(form)
  File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/analysis/analyze.py", line 42, in analyze
    check_form(form)
  File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/analysis/checks.py", line 24, in check_form
    check_completeness(form)
  File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/analysis/checks.py", line 62, in check_completeness
    raise FormError, (m, "Index does not appear exactly twice in term.")
ffc.common.exceptions.FormError

Am I doing something wrong or is there work to be done so that FFC can
compile these forms?

Jed

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Aug 2007 19:59:15 +0200

This seems to be a bug in FFC. I simplified your example and get the 
same error with the following example:

   from numpy import square as sqr

   P1 = VectorElement("Lagrange", "triangle", 1)

   v = TestFunction(P1)
   f = Function(P1)

   L = (sqr(f[0].dx(0)) + sqr(f[1].dx(1))) * v[0] * dx

I'll see what I can do. I know who to ask... ;-)

/Anders

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

---

From: Anders Logg <[email-masked]>
Date: Wed, 22 Aug 2007 12:53:52 +0200

This should be fixed now (thanks to Marie).

/Anders