Subject: Re: [FFC-dev] vertexeval function

From: Anders Logg <[email-masked]>
Date: Mon, 29 May 2006 20:20:05 +0200

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

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 29 May 2006 20:34:13 +0200

> pick the value at the vertices and we know how the dofs are ordered.

Sounds like the simplest and best approach.

> Yes, but presumably it's an easy system to solve and we do it only

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

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

---

From: Anders Logg <[email-masked]>
Date: Mon, 29 May 2006 20:43:38 +0200

> Sounds like the simplest and best approach.

ok, I'll take a look tomorrow.
 
> quasi-time. I'd like to have both approaches (projection and

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.

/Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 29 May 2006 20:56:59 +0200

> 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