Subject: [FIAT-dev] Re: [FFC-dev] Element area

From: Anders Logg <[email-masked]>
Date: Tue, 12 Jul 2005 09:59:16 -0500

This is the problem that Garth spotted, but it looks like you are
using an old version of FIAT since

  1. you have shapes_new in there, and

  2. the debug message "[[((-1.0, -1.0, -1.0),), " etc is still there

/Anders

On Tue, Jul 12, 2005 at 01:13:08PM +0200, [email-masked] wrote:
> 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
> > [email-masked]
> > http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
> >
> 
> 
> 
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
> 

-- 
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: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

On Tue, Jul 12, 2005 at 01:24:25PM +0200, Garth N. Wells wrote:
> 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 [email-masked]:
> 
> > Hi Anders,
> >