dolfin team mailing list archive

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

On Mon, 2006-06-12 at 09:15 +0200, Anders Logg wrote:
> 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.
> 

I spoke with Kristian about this last week, and our conclusion was that
the best approach (at this point) is to remove the export of constant
functions to a VTK file, and just project onto a linear basis. We can
reintroduce the output of constant functions later when 

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

is available.

Garth

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