Approximate area in a 2D field

1

Hello,

I am trying to figure out a way to calculate the (approximated) area of a 2D field that meets certain criteria, for example, area in the field that has field value greater than zero. I was thinking that solving the problem with constant elements, and if that value in that element is greater than zero, add the area of that element into the calculation, yet not sure how to perform it.

in the following code I create an example field, simply a circle with radius of 0.2 in a unit square, the field value inside the circle is 1 and outside of it is 0, if the area approximation if correct, it should be about 0.1257

import dolfin
mesh = UnitSquareMesh(100, 100)
#V = FunctionSpace(mesh, 'R', 0)  ## not sure why the constant element couldn't get the expression interpolated correctly
V = FunctionSpace(mesh, 'Lagrange', 1)
u0 = Expression('(pow(x[0]-0.5,2) + pow(x[1]-0.5,2) > pow(R,2)) ? m : n',
                 degree=2,R=0.2,m=0.0,n=1.0)
u = interpolate(u0, V)

Any help would be appreciated, thank you!

2

The R space is a constant single constant for the whole mesh, not per element. A per element function space is the DG, 0 space:

from dolfin import *
mesh = UnitSquareMesh(100, 100)

R = 0.2
V = FunctionSpace(mesh, 'Lagrange', 1)
u0 = Expression('(pow(x[0]-0.5,2) + pow(x[1]-0.5,2) > pow(R,2)) ? m : n',
                 degree=2,R=R,m=0.0,n=1.0)
u = interpolate(u0, V)
print("CG", assemble(u*dx))

Q = FunctionSpace(mesh, "DG", 0)
u_DG = interpolate(u0, Q)
print("DG", assemble(u_DG*dx))
print("Exact", pi*R**2)
3

Hello Dokken

Thanks for correcting my misunderstanding of the R space! I like your approach to calculate the area! However, this approach only works with this particular test case where the field of interest is a binary field. If the parameter n is changed to 2.0 in u0, or m is changed to -1.0, assembling u*dx would give an incorrect answer. For a more realistic application, n would not be a constant value. My end goal is to obtain the area of the field where there is positive values. I wonder if there is a way to convert all positive value in u to one and all negative value in u to zero.

4

Please provide a minimal example of such a field, as I do not understand what you want to do.

The resulting DG 0 function is a binary function, that is one or zero on each element (or 2 and -1 if you change m and n).

5

One example of such field would be:

R=0.2
u0 = Expression('(pow(x[0]-0.5,2) + pow(x[1]-0.5,2) > pow(R,2)) ? m*x[0] : n*x[0]',
                 degree=2,R=R,m=-1.5,n=2.0)
u = interpolate(u0, V)

I just tried something with Conditional(), it works to convert the field to binary, but the calculation accuracy is very bad

from dolfin import *
mesh = UnitSquareMesh(100, 100)

R = 0.2
u0 = Expression('(pow(x[0]-0.5,2) + pow(x[1]-0.5,2) > pow(R,2)) ? m : n',
                 degree=2,R=R,m=-1.0,n=3.0)

V = FunctionSpace(mesh, 'Lagrange', 1)
u = interpolate(u0, V)

F = conditional( gt(u,Constant(0.0)), Constant(1.0), Constant(0.0) )
print("CG-1 conditional ", assemble(F*dx))

This gives me 0.1326999999999935 while the exact value should be 0.12566370614359174 computed from pi*R**2

6

I just did another test using Conditional() on the DG,0 space, the accuracy seems to improve a lot, but I ran into plotting issue:

Object cannot be plotted directly, projecting to piecewise linears.
Traceback (most recent call last):
  File "foo.py", line 56, in <module>
    h = plot(u_DG, mode='color',cmap=cm2,vmin=0.0,vmax=1.0)
......
AttributeError: 'PolyCollection' object has no property 'mode'

I am not very familiar with DG, is there any better way to work out the area calculation using CG,1? If not, I will use DG for now and close this discussion until I run into further complications.

7

You can visualize DG in paraview using XDMFFile.write_checkpoint: