DOLFIN

Computing numerical integral of a function on a cell

Asked by James Avery

Hi all,

I have a very simple question that, for some reason, I'm having a hell
of a time figuring out. I need to evaluate the integral of some
function on each cell -- in order to determine whether or not to
subdivide the cell in question. With Deal.II, I would do something like
the following:

Let 'density' be the function that should be integrated on the cell.
Then

  fe_values.reinit (cell);
  quadrature_points = fe_values.get_quadrature_points();

  // Compute the function value in the cell's quadrature points
  density.value_list(quadrature_points,function_values);

  // Integral is function values dotted with the Jacobian-times-weights.
  double cell_density = 0;
  for(size_t q=0;q<N_q;q++)
    cell_density += function_values[q]*fe_values.JxW(q);

computes the cell-integral, where the types above are as follows:

  Function<dim> density; // This is a continuous function, like the FEniCS 'Expression' class.
  FEValues<dim> fe_values(fe, quadrature_formula, update_quadrature_points | update_JxW_values);
  vector< Point<dim> > quadrature_points(N_q);

I've been searching both the documentation and the source code for a
long time, and I still haven't figured out how to do this with
Dolfin. It seems I'm completely blind; a bit of help crossing the
street would be much appreciated!

Many thanks in advance.

Cheers,
James

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Marie Rognes
Solved:
Last query:
Last reply:
Revision history for this message
Marie Rognes (meg-simula) said : 		#1

Define a DG0 function space on your mesh, and define a test function v on this function space. (The DG0 basis functions will be 1 in one of the cells and 0 on all the rest). Assemble the form density*v*dx to yield a vector, where each entry corresponds to the integral of 'density' over a cell.

In other words

density = ...
DG0 = FunctionSpace(mesh, "DG", 0)
v = TestFunction(DG0)
L = density*v*dx
b = assemble(L)

Could this be what you are looking for?

Revision history for this message
James Avery (avery-diku) said : 		#2

Hi Marie, thanks so much for your swift reply!

I'm not quite sure I understand, though: Once this bit of UFL is FFC'ed, how do I then compute the cell integrals of the density from Dolfin/C++?

Just to make sure we're on the same page:

1. 'density' is a spatial scalar function given at runtime (I suppose the UFL data type would be Coefficient(mesh)?).
2. For each cell, I need to evaluate the cell integral of 'density'. If there is more than a given amount in a cell, I mark it for subdivision. This proceeds iteratively until each cell contains no more than a given amount of density.

Now, in your above example, I understand that 'b' defines a vector of length (number of cells), the entries being the cell integrals of 'density'.

1. This is UFL-code, which is compiled with FFC, is that correct? From the C++-program, how do I then execute the cell-integration and access the results?
2. How are the integrals calculated, then?
3. Can't I just get the quadrature points and quadrature weights in an easy way? Then the operation is essentially trivial.

I think I may still be harboring some confusion about how UFL/FFC and Dolfin interoperate -- about what is done offf-line by FFC, and what can be done on-line by Dolfin. I hope the confusion doesn't make it too difficult to answer my question.

Thanks very much for your help!

Cheers,
James

Revision history for this message
Best Marie Rognes (meg-simula) said : 		#3

Oh, you are in c++, ok. Define the form in a UFL file, compile it with FFC, and then assemble the resulting form in DOLFIN.
Sample code included below. (Not compiled/tested, all of DOLFIN installations are currently broken)

Some more comments:

(1) For more info on the UFL/FFC/DOLFIN interplay, consider looking at some of the cpp demos or reading for instance the DOLFIN chapter in the FEniCS book.
(2) This way of doing it should be more trivial (and perhaps more generally useful) that evaluating the quadrature yourself.

--- DensityIntegral.ufl ---
# Compile using FFC with
# >> ffc -l dolfin -O DensityIntegral.ufl

# Define some element and the to-be given coefficient for density
cell = "tetrahedron"
V = FiniteElement("CG", cell, 2)
density = Coefficient(V)

# Define DG element
Q = FiniteElement("DG", cell, 0)
v = TestFunction(Q)

# Define form to be used to evaluate \int_T density for each cell T
L = density*v*dx

--- main.cpp ---

#include <dolfin.h>
#include "DensityIntegral.h"

using namespace dolfin;

// Density term (right-hand side)
class Density : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;
    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
  }
};

int main()
{
  // Create mesh
  UnitCube mesh(4, 4, 4);

  // Create density function (just an example)
  Density density;

  // Create test space
  DensityIntegral::FunctionSpace Q(mesh);

  // Create linear form
  DensityIntegral::LinearForm L(Q);
  L.density = density;

  // Assemble form into vector
  Vector b;
  assemble(b, L);
  info(b);

  // Next steps ...

  // 1. Filter vector to yield MeshFunction of bools over cells
  // ('markers') indicating which to refine or not

  // 2. Call
  // new_mesh = refine(mesh, markers)
  // to refine mesh locally

  return 0;
}

Revision history for this message
Kristian B. Ølgaard (k.b.oelgaard) said : 		#4

Put this in a file say Foo.ufl:

DG0 = FiniteElement("DG", triangle, 0)
element = FiniteElement("CG", triangle, 2)

v = TestFunction(DG0)
f = Coefficient(element)
L = f*v*dx

Compile it with FFC, don't forget the '-l dolfin' argument:

$ ffc -l dolfin Foo.ufl

which will produce Foo.h

Then put the following in a main.cpp file:

#include <dolfin.h>
#include "Foo.h"

using namespace dolfin;

class Source : public Expression
{
  public:

    void eval(Array<double>& values, const Array<double>& x) const
    {
      values[0] = sin(x[0]*x[1]);
    }
};

int main()
{
  UnitSquare mesh(3,3);

  Avg::FunctionSpace V(mesh);
  Avg::LinearForm L(V);

  Source f;
  L.f = f;

  Vector b;
  assemble(b, L);
  info(b, true);

  return 0;
}

Now you have your values in the vector b. What you'll do with it is up to you.

Kristian

On 6 May 2011 16:31, James Avery <email address hidden> wrote:
> Question #156272 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/156272
>
>    Status: Answered => Open
>
> James Avery is still having a problem:
> Hi Marie, thanks so much for your swift reply!
>
> I'm not quite sure I understand, though: Once this bit of UFL is FFC'ed,
> how do I then compute the cell integrals of the density from Dolfin/C++?
>
> Just to make sure we're on the same page:
>
> 1. 'density' is a spatial scalar function given at runtime (I suppose the UFL data type would be Coefficient(mesh)?).
> 2. For each cell, I need to evaluate the cell integral of 'density'. If there is more than a given amount in a cell, I mark it for subdivision. This proceeds iteratively until each cell contains no more than a given amount of density.
>
> Now, in your above example, I understand that 'b' defines a vector of
> length (number of cells), the entries being the cell integrals of
> 'density'.
>
> 1. This is UFL-code, which is compiled with FFC, is that correct? From the C++-program, how do I then execute the cell-integration and access the results?
> 2. How are the integrals calculated, then?
> 3. Can't I just get the quadrature points and quadrature weights in an easy way? Then the operation is essentially trivial.
>
> I think I may still be harboring some confusion about how UFL/FFC and
> Dolfin interoperate -- about what is done offf-line by FFC, and what can
> be done on-line by Dolfin. I hope the confusion doesn't make it too
> difficult to answer my question.
>
> Thanks very much for your help!
>
> Cheers,
> James
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>

Revision history for this message
James Avery (avery-diku) said : 		#5

Hi Marie,

Thanks a lot! It seems I had mistaken your python code for UFL, causing some of my confusion. The code in the second reply makes a lot of sense to me. I'll have a go at getting it to work -- but it looks like it shouldn't be much trouble. :)

Thanks again for your swift and informative reply!

Cheers,
James

Revision history for this message
James Avery (avery-diku) said : 		#6

...and thank you Kristian as well. :)

Revision history for this message
James Avery (avery-diku) said : 		#7

Thanks Marie Rognes, that solved my question.