dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Nice. The obvious thing would be to implement this in DiscreteFunction
and map that to the function call

  virtual void Function::eval(real* values, const real* x) const;

so that any Function (discrete, constant or user-defined) can be
evaluated at an arbitrary point.

It should be possible to implement this for any kind of element, and
the code will look about the same as the code you have done for simple
elements.

We might add some kind of caching so that evaluation at multiple
points that lie close to each other is efficient. (But maybe GTS is
smart and handles this already.)

-- 
Anders


On Sat, Feb 16, 2008 at 01:13:16PM +0100, Dag Lindbo wrote:
> This has been a fairly long discussion, but I have code in place now that
> does the evaluation at an arbitrary point (2d only, but extension is easy)
> - so thanks!
> 
> I settled with a solution where the evaluation works on simple elements
> (single valued, non-mixed etc). The calling context is responsible for
> takning the form and stripping it down to simple finite elements and
> dof-maps, and calling the evaluation for each sub-function as required
> (example below).
> 
> As for interfaces in the future, I don't think that this kind of
> evaluation belongs in the DiscreteFunction (since it probably shouln't be
> encouraged due to the computational complexity). However, it would be nice
> to be able to extract the element and DofMap from the discrete function
> (as in pre-ufc dolfin), instead of having to invoke the form.
> 
> For future reference if anyone is interested this is the code I use
> (hoping that it is correct :) )
> 
> // Evaluate func at points (xx,uu), put result into uu
> void evaluate_FE_func(const Function& func,
> 		      const ufc::finite_element* fe,
> 		      ufc::dof_map* ufc_dof_map,
> 		      const double* xx, const double* yy,
> 		      double* uu, int N)
> {
>   int i,j;
>   real element_sum;
> 
>   IntersectionDetector idt;
>   idt.init(mesh);
> 
>   DofMap* dof_map = new DofMap(*ufc_dof_map, mesh);
>   unsigned int* dofs = new unsigned int[dof_map->local_dimension()];
>   real* dof_values = new real[fe->space_dimension()];
>   real* basis_val = new real[fe->value_dimension(0)];
>   real* coord = new real[2];
> 
>   Vector vec = func.vector();
> 
>   for(i=0;i<N;i++)
>     {
>       Point p(xx[i],yy[i]);
> 
>       Array<unsigned int> cells;
>       idt.overlap(p,cells);
>       int cell_nr = cells[0];
>       Cell cell(mesh,cell_nr);
>       UFCCell ufc_cell(cell);
> 
>       dof_map->tabulate_dofs(dofs,ufc_cell);
>       vec.get(dof_values, dof_map->local_dimension(), dofs);
> 
>       coord[0] = p.x();
>       coord[1] = p.y();
> 
>       element_sum=0;
>       for(j=0;j<fe->space_dimension();j++)
> 	{
> 	  fe->evaluate_basis(j,basis_val,coord,ufc_cell);
> 	  element_sum += dof_values[j]*basis_val[0];
> 	}
> 
>       uu[i] = element_sum;
>     }
> }
> 
> Example of calling function (Taylor-Hood element, velocity and pressure):
> 
> // Evaluate velocity components at all points (xx,yy)
> // a and u are members of the class
> void some_class::eval_velocity_TH(const double* xx, const double* yy,
> 		                  double* u_mx, int N)
> {
>   // Strip down to simple elements
>   ufc::finite_element* fe_up = a->form().create_finite_element(0);
>   ufc::finite_element* fe_u = fe_up->create_sub_element(0);
>   ufc::finite_element* fe_ux = fe_u->create_sub_element(0);
>   ufc::finite_element* fe_uy = fe_u->create_sub_element(1);
> 
>   // Strip down to simple dof-maps
>   ufc::dof_map* dof_map_up = a->form().create_dof_map(0);
>   ufc::dof_map* dof_map_u = dof_map_up->create_sub_dof_map(0);
>   ufc::dof_map* dof_map_ux = dof_map_u->create_sub_dof_map(0);
>   ufc::dof_map* dof_map_uy = dof_map_u->create_sub_dof_map(1);
> 
>   // x-component
>   Function ux(u[0]);
>   evaluate_FE_func(ux,fe_ux,dof_map_ux,xx,yy,u_mx,N);
> 
>   // y-component
>   Function uy(u[1]);
>   evaluate_FE_func(uy,fe_uy,dof_map_uy,xx,yy,u_mx+N,N);
> 
> }
> 
>    /Dag
> 
> >
> > This is what the functions create_sub_element and create_sub_dof_map
> > are for. If you have an element (or dof_map), you can pick out the
> > element for any sub systems. You can have arbitrary depths of sub
> > system nesting. For example in your case, you have
> >
> >   w = (u, p)
> >
> > so you can ask for the finite element for w and then ask that element
> > for the elements for u and p. Then again, u is a system of size d
> > (space dimension) so then you can again ask for the sub systems for
> > the individual components.
> >
> > Look in the UFC manual for details on the create_sub_foo functions.
> >
> >
> 
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev