dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 06, 2006 at 03:09:58PM +0200, Kristian Oelgaard wrote:
> Quoting Anders Logg <logg@xxxxxxxxx>:
> 
> > On Tue, Jun 06, 2006 at 12:05:10PM +0200, Kristian Oelgaard wrote:
> > > 
> > > 
> > > Hi,
> > > I have some difficulties with the subfunction facility in Dolfin. I'm
> > interested
> > > in saving the normal and shear components of a stress function (in 3D) for
> > post
> > > processing purposes.
> > > 
> > > the stress function is initialised in this way:
> > > 
> > > 	stress.init(mesh, a_strain.trial());
> > > 
> > > and the trial element 'a_strain.trial()' is defined as:
> > > 
> > > 	elementB = FiniteElement("Discontinuous vector Lagrange", "tetrahedron",
> > 0, 3)
> > > 	elementC = FiniteElement("Discontinuous vector Lagrange", "tetrahedron",
> > 0, 3)
> > > 
> > > 	strain_element = elementB + elementC
> > > 
> > > 	(wn,ws) = TrialFunctions(strain_element)
> > > 
> > > Now, when doing:
> > > 
> > > 	Function n_stress = stress[0];
> > > 	Function s_stress = stress[1];
> > > 
> > >  	File file_normal_stress("normal_stress.pvd");
> > >  	File file_shear_stress("shear_stress.pvd");
> > > 
> > > 	while( t < T)
> > > 	{
> > > 
> > >     nonlinear_solver.solve(nonlinear_problem, x);
> > > 
> > > 		n_stress = stress[0];
> > > 		s_stress = stress[1];
> > > 
> > > 	  file_normal_stress << n_stress;
> > > 	  file_shear_stress  << s_stress;
> > > 
> > >   }
> > > 
> > > the values of the two files are exactly the same as if only the first half
> > of
> > > stress.vector() is used as a basis for both subfunctions.
> > > 
> > > Funny thing is that a similar approach works fine in the elasticity demo.
> > Only
> > > difference is that in this case the strains are computed in FFC, where my
> > > stresses are computed in the non-linear problem, on element level, and
> > then
> > > written into the stress.vector(). So, what am I doing wrong?
> > > 
> > > Regards,
> > > 
> > > Kristian
> > > 
> > > _______________________________________________
> > > DOLFIN-dev mailing list
> > > DOLFIN-dev@xxxxxxxxxx
> > > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> > 
> > Which versions of DOLFIN and FFC are you using? We recently updated
> > the evaluation at vertices (that is called when saving to file), so we
> > might have broken something.
> 
> The versions i used were a couple of days old, but to be sure I just grapped new
> versions of FFC and Dolfin, and I still get problems.
> 
> > 
> > Check the code generated by FFC, what does the function(s) vertexeval
> > say?
> 
> In class trial function:
> 
> 	public:
> 
> void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const
>     {
>       // FIXME: Temporary fix for Lagrange elements
>       vertex_nodes[0] = vertex;
>       int offset = mesh.numCells();
>       vertex_nodes[1] = offset + vertex;
>       offset = offset + mesh.numCells();
>       vertex_nodes[2] = offset + vertex;
>       offset = offset + mesh.numCells();
>       vertex_nodes[3] = offset + vertex;
>       offset = offset + mesh.numCells();
>       vertex_nodes[4] = offset + vertex;
>       offset = offset + mesh.numCells();
>       vertex_nodes[5] = offset + vertex;
>     }
> 
> in private:
> 
>     class SubElement_0 : public dolfin::FiniteElement
>     {
>     public:
> 
> 			void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const
>       {
>         // FIXME: Temporary fix for Lagrange elements
>         vertex_nodes[0] = vertex;
>         int offset = mesh.numCells();
>         vertex_nodes[1] = offset + vertex;
>         offset = offset + mesh.numCells();
>         vertex_nodes[2] = offset + vertex;
>       }
> 		}
> 
>     class SubElement_1 : public dolfin::FiniteElement
>     {
>     public:
> 
>       void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh&
> mesh) const
>       {
>         // FIXME: Temporary fix for Lagrange elements
>         vertex_nodes[0] = vertex;
>         int offset = mesh.numCells();
>         vertex_nodes[1] = offset + vertex;
>         offset = offset + mesh.numCells();
>         vertex_nodes[2] = offset + vertex;
>       }
> 		}
> 
> I guess the two functions look too much alike?
> 
> > 
> > You're using discontinuous Lagrange which I don't think is supported
> > correctly by vertexeval().
> 
> Could very well be, if I try to project the strains onto a discontinuous basis
> in elasticity demo I get the same error.
> 
> /Kristian

Here's the problem...

First of all, vertexeval() is just a quick hack that works for
Lagrange basis functions. It assumes that degrees of freedom
associated with vertices are always stored first for each
component. So if you have a vector x for a vector-valued function that
has degrees of freedom at each vertex and each edge, then

    x[0] = value of component 0 at vertex 0
    x[1] = value of component 0 at vertex 1
    ...
    x[num_vertices] = value of component 0 at edge 0
    x[num_vertices + 1] = value of component 0 at edge 1
    ...
    x[num_vertices + num_edges] = value of component 1 at vertex 0
    x[num_vertices + num_edges] = value of component 1 at vertex 1
    ...
    x[2*num_vertices + num_edges] = value of component 1 at edge 0
    x[2*num_vertices + num_edges + 1] = value of component 1 at edge 1
    ...

For a mixed element, we use the same pattern but all dofs associated
with the first sub function are stored first (like above), then the
second sub function etc.

So, vertexeval() only works for Lagrange.

Now, as you point out, Garth has been clever and worked out a special
solution for piecewise constant elements in VTKFile. This seems to be
a correct implementation for vector-valued elements, but it does not
work for mixed elements since the mixed offset is not known. I'm a
little reluctant to export mixed_offset in the Function interface, so
the best solution would be to add something like

    real operator() (const Cell& cell, uint i = 0);

to the Function class and put the implementation in DiscreteFunction
where mixed_offset is known.

Then, we would have

    real operator() (const Vertex& vertex, uint i = 0);
    real operator() (const Cell& cell, uint i = 0);

in the Function interface. If in addition, we implement projection to
piecewise linears and constants, this should have postprocessing
covered for a fairly wide range of problems. Right now, one can do
manual projection by defining the variational problem with FFC, but we
should maybe consider generating the projection for every element that
is used in an FFC form so it's always available.

/Anders