Subject: [FFC-dev] Element area

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 08 Jul 2005 23:02:04 +0200

Hi Anders,

I'm trying to stabilise a problem and for this I need a measure of the element
area. Do you have any plans to make this sort of information available through
ffc? I have seen that Johan H. hasn't yet defined the element length in the
Navier-Stokes module yet. Johan, do you have a plan how to do this? It could be
done easily by defining a constant in the .form file and then adding a line c4 =
sqrt(det), but this would destroy the beauty of the ffc generated file. 

A more general question - do you have a plan on how more complicated nonlinear
problems, such as plasticity, will fit in with ffc + DOLFIN? I mean here solving
equations in which the coefficients of an equation themselves come from the
solution of a non-trivial nonlinear problem, such a plasticity return mapping
algorithm.

Garth

---

From: Anders Logg <[email-masked]>
Date: Fri, 8 Jul 2005 16:45:35 -0500

> sqrt(det), but this would destroy the beauty of the ffc generated file. 

I will probably add a default Function h(x) that gives you the size of
the element. Do you need the side length (diameter) or the area?

This will probably not be until after the summer though, so what you
can do now is to define a piecewise constant function that is defined
as the mesh size or area or whatever, like so:

   P0 = FiniteElement("Discontinuous Lagrange", "triangle", 0)
   h  = Function(P0)

and then define a corresponding Function in DOLFIN from a Vector and
set the values of the Vector to the local mesh size (x(i) area of
triangle i).

> solution of a non-trivial nonlinear problem, such a plasticity return mapping

This is already implemented. When you define a Function in FFC, then
this corresponds to a Function in DOLFIN, and the DOLFIN Function
could be either defined by an expression or it can be the solution of
another problem.

The Function concept is really central to both FFC and DOLFIN, and it
is the key to working with FFC and DOLFIN. Each Function in FFC
corresponds to a Function in DOLFIN, and it should not matter how the
Function is defined in DOLFIN. There are a couple of things missing,
like reading Functions from file, projection between different
function spaces, but the basic stuff is implemented. (Take a look at
Function.cpp, it's pretty short.)

Ask again if something is unclear.

/Anders

---

From: [email-masked]
Date: Tue, 12 Jul 2005 13:13:08 +0200 (MEST)

Hi Anders,

I have started to look at how to solve the stabilization problem in the
best way. I first tried the approach with a piecewise constant function
for the stabilization parameter, but I seem to have missed something? Do
you find anything obvious wrong with the form-file below?

/Johan

name = "NSEContinuity"
scalar_element = FiniteElement("Lagrange", "tetrahedron", 1)
vector_element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
disc_scalar_element = FiniteElement("Discontinuous Lagrange",
"tetrahedron", 0)

q = BasisFunction(scalar_element)
p = BasisFunction(scalar_element)
uc = Function(vector_element)
f = Function(vector_element)
d1 = Function(disc_scalar_element)

i0 = Index()
i1 = Index()
i2 = Index()

a = q.dx(i0)*p.dx(i0)*dx
L = q.dx(i0)*f[i0]*dx - q.dx(i0)*uc[i1]*uc[i0].dx(i1)*dx -
(1/d1)*q*uc[i0].dx(i0)*dx

skalman:~/local> ffc NSEContinuity.form
This is FFC, the FEniCS Form Compiler, version 0.1.8.
For further information, go to http://www/fenics.org/ffc/.

Parsing NSEContinuity.form
Output written to NSEContinuity.py
[[((-1.0, -1.0, -1.0),), ((1.0, -1.0, -1.0),), ((-1.0, 1.0, -1.0),),
((-1.0, -1.0, 1.0),)], [(), (), (), (), (), ()], [(), (), (), ()], [()]]
((-1.0, -1.0, -1.0), (1.0, -1.0, -1.0), (-1.0, 1.0, -1.0), (-1.0, -1.0, 1.0))
Warning: element untested
Traceback (most recent call last):
  File "/usr/bin/ffc", line 81, in ?
    main(sys.argv[1:])
  File "/usr/bin/ffc", line 59, in main
    execfile(outname)
  File "NSEContinuity.py", line 32, in ?
    disc_scalar_element = FiniteElement("Discontinuous Lagrange",
"tetrahedron", 0)
  File "/home/hoffman/local/lib/python/ffc/compiler/finiteelement.py",
line 66, in __init__
    self.element = DiscontinuousLagrange(fiat_shape, degree)
  File "/home/hoffman/local/lib/python/FIAT/DiscontinuousLagrange.py",
line 87, in DiscontinuousLagrange
    if n == 0: return P0.P0( shape )
  File "/home/hoffman/local/lib/python/FIAT/P0.py", line 35, in __init__
    Udual = P0Dual( shape , U )
  File "/home/hoffman/local/lib/python/FIAT/P0.py", line 21, in __init__
    d = shapes_new.dims[ shape ]
AttributeError: 'module' object has no attribute 'dims'

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

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 12 Jul 2005 13:24:25 +0200

Hi Johan,

I ran into this problem with FIAT when defining a discontinuous constant scalar
basis function in ffc (the vector version works). I sent an email with a patch
to the FIAT-dev list about this. Below is what I have, which multiplies the
masss matrix by the element length (just for testing) and seems to work well now.

scalar       = FiniteElement("Lagrange", "triangle", 1)
vector       = FiniteElement("Vector Lagrange", "triangle", 1)
scalar_const = FiniteElement("Discontinuous Lagrange", "triangle", 0)

v  = BasisFunction(scalar)     # test function
u  = BasisFunction(scalar)     # trial function
f  = Function(scalar)          # source term
h  = Function(scalar_const)    # element size h

a = h*v*u*dx 
L = v*f*dx

Then, following Anders's advice, I set up a function

Vector hvector;
hvector.init(mesh.noCells());
for (CellIterator cell(mesh); !cell.end(); ++cell)
{
  hvector((*cell).id()) = sqrt(2*(*cell).volume());
}
Function h(hvector);

Garth

---

From: Anders Logg <[email-masked]>
Date: Tue, 12 Jul 2005 10:02:10 -0500

Rob is at a conference, so I will go ahead and apply your patch (the
one-line fix) to the CVS version of FIAT. (And if I do it within the next
minute, it will make the 10 o'clock tar-ball...)

/Anders

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/

---

From: Anders Logg <[email-masked]>
Date: Tue, 12 Jul 2005 10:06:40 -0500

Garth, there is also a function called Cell::diameter() that you can
use directly to get the mesh size.

/Anders