ffc team mailing list archive

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

On Mon, 2006-05-29 at 20:43 +0200, Anders Logg wrote:
> On Mon, May 29, 2006 at 08:34:13PM +0200, Garth N. Wells wrote:
> > On Mon, 2006-05-29 at 20:20 +0200, Anders Logg wrote:
> > > On Mon, May 29, 2006 at 07:57:35PM +0200, Garth N. Wells wrote:
> > > > On Mon, 2006-05-29 at 19:32 +0200, Anders Logg wrote:
> > > > > On Mon, May 29, 2006 at 07:10:23PM +0200, Garth N. Wells wrote:
> > > > > > For evaluating functions at vertxes, FFC writes a function
> > > > > > 
> > > > > > void vertexeval(real values[], unsigned int vertex, const real x[],
> > > > > > const Mesh& mesh) const
> > > > > > 
> > > > > > Would there be a problem in changing "const real x[]" to "const
> > > > > > GenericVector x"? This would make it easier to represent a function
> > > > > > using different vector data structures.
> > > > > > 
> > > > > > Garth
> > > > > 
> > > > > The problem would be that vertexeval() would then need to call
> > > > > GenericMatrix::operator() (uint i) which is slow for PETScSparseMatrix.
> > > > > 
> > > > > But perhaps calling PETScSparseMatrix::array() is just as slow?
> > > > > 
> > > > 
> > > > When manipulating a lot of term, PETScSparseMatrix::array() is much
> > > > faster. I think that this has to with the calls
> > > > 
> > > >   VecAssemblyBegin(x);
> > > >   VecAssemblyEnd(x);
> > > > 
> > > > after accessing each entry.
> > > > 
> > > > > vertexeval() is just a temporary fix for Lagrange elements. Somehow,
> > > > > we need to compute the values at all vertices so we can write the
> > > > > solution to file. What we should really do would be to compute the
> > > > > projection onto piecewise linears and then just pick the values.
> > > > > 
> > > > 
> > > > OK. A simple temporary fix would be if vextexeval returned an array of
> > > > integers containing the positions in the big vector. Then the values
> > > > could be pulled as a block out of the vector. 
> > > 
> > > Even simpler, maybe we can just stick the code in DiscreteFunction and
> > > remove vertexeval()? It's special for Lagrange anyway. We just need to
> > > pick the value at the vertices and we know how the dofs are ordered.
> > > 
> > 
> > Sounds like the simplest and best approach.
> 
> ok, I'll take a look tomorrow.
>  
> > > > > Maybe we could always generate code for the projection onto piecewise
> > > > > linears when a form is compiled in FFC? If U is in V' and V is the
> > > > > space of piecewise linears, then we need to compute the projection
> > > > > W of U onto V:
> > > > > 
> > > > >     int v W dx = int v U dx for all v in V
> > > > > 
> > > > > On matrix form, we have
> > > > > 
> > > > >     MW = b
> > > > > 
> > > > > where b represents the right-hand side int U v dx.
> > > > > 
> > > > > The mass matrix M is already computed by default in DOLFIN so we don't
> > > > > need to generate it. We only need to generate the linear form v*U*dx.
> > > > > 
> > > > 
> > > > . . . and solve the linear system.
> > > 
> > > Yes, but presumably it's an easy system to solve and we do it only
> > > once at postprocessing.
> > > 
> > 
> > This can add up for time-dependent problems and problems parametrised by
> > quasi-time. I'd like to have both approaches (projection and
> > interpolation). 
> > 
> > Garth
> 
> The problem with interpolation is that its not always possible. For
> say Crouzeix-Raviart, we can't just pick the values at the vertices.
> Maybe we can treat Lagrange as a special case.
> 

It could be worth having it as a special case since an arbitrary
function will generally be projected onto a Lagrange basis for
post-processing.

Garth

> /Anders
> 
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev