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

I need to think about this (and not right now). Will try to think
about it tomorrow...

/Anders