dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Jun 12, 2006 at 09:45:00AM +0200, Garth N. Wells wrote:
> 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

Sounds good. Does the projection to linears give reasonable results?
Or does it introduce any artifacts?

/Anders