Strange behaviour of assemble

m.valdez · September 20, 2019, 4:28pm

Hi All,

I’m getting very strange results when calculating integrals using assemble over subdomains, consider the MWE:

import dolfin as fem

mesh = fem.UnitSquareMesh(100, 100)

class SubSquare(fem.SubDomain):
    def __init__(self, xi, xf, yi, yf, **kwargs):
        self.xi, self.xf = xi, xf
        self.yi, self.yf = yi, yf
        super().__init__(**kwargs)
    def inside(self, x, on_boundary):
        return self.xi <= x[0] <= self.xf and self.yi <= x[1] <= self.yf

cellfn = fem.MeshFunction('size_t', mesh, mesh.topology().dim())
cellfn.set_all(0)

bl = SubSquare(0.0, 0.5, 0.0, 0.5) # bottom left
br = SubSquare(0.5, 1.0, 0.0, 0.5) # bottom right
tr = SubSquare(0.5, 1.0, 0.5, 1.0) # top right
tl = SubSquare(0.0, 0.5, 0.5, 1.0) # top left

bl.mark(cellfn, 1) # bottom left
br.mark(cellfn, 2) # bottom right
tr.mark(cellfn, 3) # top right
tl.mark(cellfn, 4) # top left

V = fem.FunctionSpace(mesh, 'CG', 1)
u = fem.TrialFunction(V); v = fem.TestFunction(V)

dx = fem.Measure('dx', domain= mesh, subdomain_data= cellfn)
ds = fem.Measure('ds', domain= mesh)

a = (fem.dot(fem.grad(u), fem.grad(v)) + u*v)*dx + u*v*ds
L =  v*dx

m = fem.Function(V)
fem.solve(a == L, m) # so far, everything works

# but strange results when I try to calculate some integrals, e.g.:
print(fem.assemble(m*dx)) # gives ~0.22474021790739562 (works)

# integrating over the subdomains should give roughly 1/4 the integral over
# the domain (~0.056185054476848906) due to symmetry, however:
print(fem.assemble(m*dx(1))) # gives 0.0
print(fem.assemble(m*dx(2))) # gives 0.0
print(fem.assemble(m*dx(3))) # gives 0.0
print(fem.assemble(m*dx(4))) # gives 0.0

# to make things stranger, consider
print(fem.assemble((m*1.000000000000001)*dx(1))) # gives 0.05618539993070857
print(fem.assemble((m*1.000000000000001)*dx(2))) # gives 0.056184709022991684
print(fem.assemble((m*1.000000000000001)*dx(3))) # gives 0.05618539993070186
print(fem.assemble((m*1.000000000000001)*dx(4))) # gives 0.05618470902299123

Indeed, multiplying m by any number other than one works as expected. I find this behaviour quite strange. Anybody got an idea?

nate · September 20, 2019, 8:38pm

I can’t reproduce your error. I have nonzero values of \approx 0.56185 each time. Which version of FEniCS are you working with?

m.valdez · September 20, 2019, 8:57pm

I’m using 2019.1.0 but I should add that this is all on the remote environment of colab.google with fenics installed via apt get install
Also, this behaviour is quite inconsistent as this was not happening yesterday and multiplying by a constant stopped working in later integrations in my script

m.valdez · September 23, 2019, 2:52pm

As this seems to be irreproducible and more likely related to the colab.google environment, I’ll consider this settled. Thank you, @nate

I’m posting here, however, a solution for those working in the cloud environment of colab.google:
Since what seems to be failing is using dx(i) a coefficient which takes the value 1 in a specified domain and 0 elsewhere can be used. Here’s the MWE:

import dolfin as fem

mesh = fem.UnitSquareMesh(100, 100)

class SubSquare(fem.SubDomain):
    def __init__(self, xi, xf, yi, yf, **kwargs):
        self.xi, self.xf = xi, xf
        self.yi, self.yf = yi, yf
        super().__init__(**kwargs)
    def inside(self, x, on_boundary):
        return self.xi <= x[0] <= self.xf and self.yi <= x[1] <= self.yf

class DomainSelector(fem.UserExpression):
    def __init__(self, cellfn, i, **kwargs):
        self.cellfn = cellfn.array()
        self.i = i
        super().__init__(**kwargs)
    def eval_cell(self, values, x, cell):
        values[0] = 1. if self.cellfn[cell.index] == self.i else 0.
    def value_shape(self):
        return ()

cellfn = fem.MeshFunction('size_t', mesh, mesh.topology().dim())
cellfn.set_all(0)

bl = SubSquare(0.0, 0.5, 0.0, 0.5) # bottom left
br = SubSquare(0.5, 1.0, 0.0, 0.5) # bottom right
tr = SubSquare(0.5, 1.0, 0.5, 1.0) # top right
tl = SubSquare(0.0, 0.5, 0.5, 1.0) # top left

bl.mark(cellfn, 1) # bottom left
br.mark(cellfn, 2) # bottom right
tr.mark(cellfn, 3) # top right
tl.mark(cellfn, 4) # top left

V = fem.FunctionSpace(mesh, 'CG', 1)
u = fem.TrialFunction(V); v = fem.TestFunction(V)

dx = fem.Measure('dx', domain= mesh, subdomain_data= cellfn)
ds = fem.Measure('ds', domain= mesh)

a = (fem.dot(fem.grad(u), fem.grad(v)) + u*v)*dx + u*v*ds
L =  v*dx

m = fem.Function(V)
fem.solve(a == L, m)

# the integral over 'everywhere'
print(fem.assemble(m*dx)) # gives ~0.22474021790739562 (works)

# the integrals over subdomains
print(fem.assemble(DomainSelector(cellfn, 1)*m*dx)) # gives 0.05618539993070849
print(fem.assemble(DomainSelector(cellfn, 2)*m*dx)) # gives 0.056184709022991594
print(fem.assemble(DomainSelector(cellfn, 3)*m*dx)) # gives 0.056185399930701777
print(fem.assemble(DomainSelector(cellfn, 4)*m*dx)) # gives 0.05618470902299114

# the sum of integrals over subdomains, coincides with the integral over 'everywhere'
print(fem.assemble(DomainSelector(cellfn, 1)*m*dx) +
      fem.assemble(DomainSelector(cellfn, 2)*m*dx) +
      fem.assemble(DomainSelector(cellfn, 3)*m*dx) +
      fem.assemble(DomainSelector(cellfn, 4)*m*dx)) # gives 0.224740217907393

Daniel_Garcia · October 27, 2020, 11:45am

Hi all,

I want to revive this topic since I have the exact same problem.

In my case, I am using subdomain integration to build a set of multiple RHS. Each of the RHS involves a single subdomain while the remaining ones are inactive. The solution posted on the topic is ok, but in my case, I need to compute many RHS in various iterations. Thus, I prefer to integrate just in the subdomains and not need to integrate into all the domain, since the integration outside the active single subdomain will be cero.

I am adding a simple example of what I am observing

from dolfin import *
import numpy as np
##
##
# Clase to set the subdomains
class Omega(SubDomain): 
	def __init__(self, xL, xU): 
		self.xL = xL 
		self.xU = xU 
		SubDomain.__init__(self) 
	def inside(self, x, on_boundary): 
		dim = len(x) 
		L, U = self.xL, self.xU 
		return all(between(x[i], (L[i], U[i])) for i in range(dim))
##
##
# Main code
mesh = RectangleMesh(Point(0,0),Point(1,1),6,6) 
mesh2 = RectangleMesh.create([Point(0,0),Point(1,1)],[3,3],CellType.Type.quadrilateral)
#
ncell = mesh2.num_cells() 
dim   = mesh2.topology().dim()
#
idx = np.zeros(ncell,dtype="int32") 
val = np.zeros((ncell,4)) 
for cell in cells(mesh2): 
    idx[cell.index()] = cell.global_index() 
    for i in range(dim): 
        val[cell.index(),    i] = cell.get_coordinate_dofs()[     i] 
        val[cell.index(),dim+i] = cell.get_coordinate_dofs()[-dim+i] 
# Set global coordinates of the pixels 
full_idx = MPI.comm_world.allgather(idx) 
full_val = MPI.comm_world.allgather(val) 
#
domains = MeshFunction("size_t", mesh, mesh.topology().dim()) 
domains.set_all(9999) 
#
xL_px = np.zeros(dim) 
xU_px = np.zeros(dim) 
for i in range(MPI.comm_world.Get_size()): 
    for j in range(len(full_idx[i])): 
        for k in range(dim): 
            xL_px[k] = full_val[i][j][    k] 
            xU_px[k] = full_val[i][j][dim+k] 
        Omg = Omega(xL_px, xU_px) 
        Omg.mark(domains,full_idx[i][j]) 
#
FEinfo = FiniteElement("CG",mesh.ufl_cell(),1) 
W = FunctionSpace(mesh,FEinfo) 
#
v = TestFunction(W)
#
dx_s = Measure('dx', domain=mesh, subdomain_data=domains)
#
# ----------------------------------
# Full integration
ops = v*dx_s() 
vec = assemble(ops) 
print("Full domain integration \n",vec[:])
# Results
#[0.00462963, 0.01388889, 0.01388889, 0.01388889, 0.02777778,
# 0.01388889, 0.01388889, 0.02777778, 0.02777778, 0.01388889,
# 0.01388889, 0.02777778, 0.02777778, 0.02777778, 0.01388889,
# 0.01388889, 0.02777778, 0.02777778, 0.02777778, 0.02777778,
# 0.01388889, 0.00925926, 0.02777778, 0.02777778, 0.02777778,
# 0.02777778, 0.02777778, 0.00925926, 0.01388889, 0.02777778,
# 0.02777778, 0.02777778, 0.02777778, 0.01388889, 0.01388889,
# 0.02777778, 0.02777778, 0.02777778, 0.01388889, 0.01388889,
# 0.02777778, 0.02777778, 0.01388889, 0.01388889, 0.02777778,
# 0.01388889, 0.01388889, 0.01388889, 0.00462963])
# ----------------------------------
# 4 Integrations one in the i-th subdomain only
for i in range(4): 
    ops = v*dx_s(i) 
    vec = assemble(ops) 
    print(str(i)+"-th subdomain integration \n",vec[:]) 
# Results
#[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0.]
#[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0.]
#[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0.]
#[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
# 0.]

I am working on a laptop and I installed FEniCs using conda.

volkerk · October 27, 2020, 4:06pm

I get the following output when i run your code:

Full domain integration
[0.00462963 0.01388889 0.01388889 0.01388889 0.02777778 0.01388889
0.01388889 0.02777778 0.02777778 0.01388889 0.01388889 0.02777778
0.02777778 0.02777778 0.01388889 0.01388889 0.02777778 0.02777778
0.02777778 0.02777778 0.01388889 0.00925926 0.02777778 0.02777778
0.02777778 0.02777778 0.02777778 0.00925926 0.01388889 0.02777778
0.02777778 0.02777778 0.02777778 0.01388889 0.01388889 0.02777778
0.02777778 0.02777778 0.01388889 0.01388889 0.02777778 0.02777778
0.01388889 0.01388889 0.02777778 0.01388889 0.01388889 0.01388889
0.00462963]
0-th subdomain integration
[0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0.00462963 0.
0. 0. 0. 0.01388889 0.01388889 0.
0. 0. 0. 0.00925926 0.02777778 0.00925926
0. 0. 0. 0. 0.01388889 0.01388889
0. 0. 0. 0. 0.00462963 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. ]
1-th subdomain integration
[0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.00462963
0. 0. 0. 0. 0. 0.01388889
0.01388889 0. 0. 0. 0.00925926 0.02777778
0.00925926 0. 0. 0.01388889 0.01388889 0.
0. 0.00462963 0. 0. 0. 0.
0. ]
2-th subdomain integration
[0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0.00462963 0. 0. 0. 0.01388889 0.01388889
0. 0.00925926 0.02777778 0.00925926 0.01388889 0.01388889
0.00462963]
3-th subdomain integration
[0. 0. 0. 0.00462963 0. 0.
0.01388889 0.01388889 0. 0. 0.00925926 0.02777778
0.00925926 0. 0. 0. 0.01388889 0.01388889
0. 0. 0. 0. 0. 0.00462963
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. ]

Daniel_Garcia · October 27, 2020, 7:48pm

Hi Volkerk,

It is nice that it works for you. The code manages to understand that you are asking him to integrate just on the subdomains.

In my case, the results is an array of only zeros. Big problem is that is not showing any error, it is just not integrating, or integrating but returning ceros. This is kind of dangerous, this makes me ask myself how I will be secure that FeniCS will integrate what I am asking him.

Anyway, I will be grateful If someone help me to find the reason of this problem.

dokken · October 27, 2020, 8:22pm

As @volkerk said, he cannot reproduce your issue. How Did you install fenics, and Which version are you running ?

Daniel_Garcia · October 27, 2020, 8:36pm

Hi Dokken,

I installed FEniCs using anaconda cloud and the version I have in my computer is 2019.1.0.

If you need more information just let me know.

dokken · October 27, 2020, 8:42pm

I can also not reproduce this using docker:

docker run -ti -v $(pwd):/home/fenics/shared -w /home/fenics/shared/ --rm quay.io/fenicsproject/stable:latest

python3 -c "import dolfin;print(dolfin.__version__)

2019.1.0

fenics@ee647a3cdb2c:~/shared$ python3 script.py 
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Full domain integration 
 [ 0.00462963  0.01388889  0.01388889  0.01388889  0.02777778  0.01388889
  0.01388889  0.02777778  0.02777778  0.01388889  0.01388889  0.02777778
  0.02777778  0.02777778  0.01388889  0.01388889  0.02777778  0.02777778
  0.02777778  0.02777778  0.01388889  0.00925926  0.02777778  0.02777778
  0.02777778  0.02777778  0.02777778  0.00925926  0.01388889  0.02777778
  0.02777778  0.02777778  0.02777778  0.01388889  0.01388889  0.02777778
  0.02777778  0.02777778  0.01388889  0.01388889  0.02777778  0.02777778
  0.01388889  0.01388889  0.02777778  0.01388889  0.01388889  0.01388889
  0.00462963]
Calling FFC just-in-time (JIT) compiler, this may take some time.
0-th subdomain integration 
 [ 0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.00462963  0.          0.          0.
  0.          0.01388889  0.01388889  0.          0.          0.          0.
  0.00925926  0.02777778  0.00925926  0.          0.          0.          0.
  0.01388889  0.01388889  0.          0.          0.          0.
  0.00462963  0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.        ]
Calling FFC just-in-time (JIT) compiler, this may take some time.
1-th subdomain integration 
 [ 0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.00462963  0.          0.          0.          0.
  0.          0.01388889  0.01388889  0.          0.          0.
  0.00925926  0.02777778  0.00925926  0.          0.          0.01388889
  0.01388889  0.          0.          0.00462963  0.          0.          0.
  0.          0.        ]
Calling FFC just-in-time (JIT) compiler, this may take some time.
2-th subdomain integration 
 [ 0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.00462963  0.          0.          0.          0.01388889
  0.01388889  0.          0.00925926  0.02777778  0.00925926  0.01388889
  0.01388889  0.00462963]
Calling FFC just-in-time (JIT) compiler, this may take some time.
3-th subdomain integration 
 [ 0.          0.          0.          0.00462963  0.          0.
  0.01388889  0.01388889  0.          0.          0.00925926  0.02777778
  0.00925926  0.          0.          0.          0.01388889  0.01388889
  0.          0.          0.          0.          0.          0.00462963
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.          0.          0.          0.
  0.          0.          0.          0.        ]

Daniel_Garcia · October 28, 2020, 8:56am

Hi all,

Thanks for the answers. I understand that is really hard to help me with this issue if you can not reproduce it, maybe impossible.

But my idea to revive this topic was to find someone that has similar problems and that could give me a feedback, or someone that maybe have some clue of the reason that the integration on subdomains behaves like that to me.

Good day

Daniel_Garcia · October 28, 2020, 10:59am

Hi again,

For the code I posted, if I multiply the variational form by an Expression or a Constant

e.g., Constant(1.0)vdx_s(i)

the code integrates correctly on the subdomains. It is very interesting, but I do not understand the reason. It is like the integration lacks some information that the expression or the Constant have.

Now, in the real code that I am working with and that it is a little bit more complex to include in an easy example, I multiplied the variational form by an Expression or Constant but it doesn’t work correctly.

However, using the class DomainSelector(UserExpression) that posted m.valez in this topic (second post), the code integrates correctly using both dx or dx_s(i). I will use this for now.

I am still interested in the strange behavior of the integration in my nonreproducible problem. I hope to find a reason for this and post it in the future

Thanks to all

Topic		Replies	Views
Assemble is not working with fenics_adjoint and imported mesh with MESHIO dolfin-adjoint	0	224	October 8, 2023
Question about an issue with assemble() and project() for equations defined as UFL Forms Legacy DOLFIN legacy-dolfin , ufl	2	35	August 5, 2024
Incorrect result for assembly over boundary part Errors	7	819	October 14, 2021
Problem with assemble_scalar in thermoelastic problem dolfinx	0	301	September 12, 2023
Improve integral accuracy with assemble	5	661	November 9, 2021

Strange behaviour of assemble

Related topics