dolfin team mailing list archive

[Kristian Oelgaard <K.Oelgaard@xxxxxxxxxxxxxxx>]

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

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