ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

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

/Anders