ffc team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxxxxxxx>]

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


Quoting jhoffman@xxxxxxxxxxx:

> 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'
> 
> 
> 
> 
> 
> 
> > On Fri, Jul 08, 2005 at 11:02:04PM +0200, Garth N. Wells wrote:
> >> 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.
> >
> > 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).
> >
> >> 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
> >
> > 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
> >
> > _______________________________________________
> > FFC-dev mailing list
> > FFC-dev@xxxxxxxxxx
> > http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
> >
> 
> 
> 
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
>