dolfin team mailing list archive

["Dag Lindbo" <dag@xxxxxxxxxx>]

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.
>
> --
> Anders
>