Improve integral accuracy with assemble

1

Hi,

I have found related posts but none of those solves my problem.

After solving a == L, I want to compute an integral over a subdomain of the function space, but the integral is very inaccurate.

The code is like

a = assemble(u*dx(2))

Are there any ways we can use to improve the accuracy? I tried a mesh with higher resolution but still inaccurate.

Update:

Here is a sample code. Is there a way to improve the accuracy of integral inside the circle?

from fenics import *
import numpy as np
from mshr import *
import matplotlib.pyplot as plt


# Create mesh and define function space
domain = Rectangle(Point(-0.5,-0.5), Point(0.5,0.5))
domain_setter = domain
domain.set_subdomain(1, domain_setter)

cir = Circle(Point(0, 0), 0.05)
domain.set_subdomain(2, cir)
mesh = generate_mesh(domain, 8)#UnitSquareMesh(nx, ny)

markers = MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())

non_cir_mesh = MeshView.create(markers, 1)
cir_mesh = MeshView.create(markers, 2)

V1 = FunctionSpace(cir_mesh, 'CG', 1)
V2 = FunctionSpace(non_cir_mesh, 'CG', 1)

total_vol = 1
u1_bar_true = 0.001*(1-np.pi*0.1**2) / total_vol
u2_bar_true = 0.0002*(np.pi*0.1**2) / total_vol
u1 = project(0.001, V1)
u2 = project(0.0002, V2)

dx_all = dx_halfeps_packing = Measure('dx', domain=mesh)
dx_non_cir = dx_halfeps_packing = Measure('dx', domain=non_cir_mesh)
dx_cir = dx_halfeps_packing = Measure('dx', domain=cir_mesh)


u1_bar_cal = assemble(u1*dx_non_cir)/total_vol
u2_bar_cal = assemble(u2*dx_cir)/total_vol
print(u1_bar_cal, u1_bar_true) 
print(u2_bar_cal,  u2_bar_true)
2

You can always increase the quadrature degree:

dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': degree})
a = assemble(u*dx(2))

As you have not provided any example illustrating what inaccuracy you are observing, I can’t give you any further suggestions.

3

I’ve added a sample script. Thank you.

4

First of all your code exact solution is wrong, as you are computing

but assuming radius 0.1

Changing this improves your convergence.
Following is a minimal code, also highlighting the area estimation error:

from fenics import *
import numpy as np
from mshr import *


def compute_volume(N, r=0.05):
    # Create mesh and define function space
    domain = Rectangle(Point(-0.5, -0.5), Point(0.5, 0.5))
    domain_setter = domain
    domain.set_subdomain(1, domain_setter)
    cir = Circle(Point(0, 0), r)
    domain.set_subdomain(2, cir)
    mesh = generate_mesh(domain, N)

    markers = MeshFunction(
        'size_t', mesh, mesh.topology().dim(), mesh.domains())
    File("markers.pvd") << markers
    V = FunctionSpace(mesh, 'CG', 1)

    total_vol = 1
    u1_bar_true = 0.001*(1-np.pi*r**2) / total_vol
    u2_bar_true = 0.0002*(np.pi*r**2) / total_vol
    u1 = project(0.001, V)
    u2 = project(0.0002, V)

    dx_all = Measure('dx', domain=mesh, subdomain_data=markers)
    dx_non_cir = dx_all(1)
    dx_cir = dx_all(2)

    tot_vol_cal = assemble(1*dx_all)
    u1_bar_cal = assemble(u1*dx_non_cir)/tot_vol_cal
    u2_bar_cal = assemble(u2*dx_cir)/tot_vol_cal

    volume1_true = 1-np.pi*r**2
    volume1_estimated = assemble(1*dx_non_cir)
    print(
        f"N:{N} Area square (True): {volume1_true}, (Estimated) {volume1_estimated} (Error) {(volume1_true-volume1_estimated)/volume1_true}")
    volume2_true = np.pi*r**2
    volume2_estimated = assemble(1*dx_cir)
    print(
        f"N:{N} Area circle (True): {volume2_true}, (Estimated) {volume2_estimated} (Error) {(volume2_true-volume2_estimated)/volume2_true}")

    print(N, u1_bar_cal, u1_bar_true, (u1_bar_cal-u1_bar_true)/u1_bar_true)
    print(N, u2_bar_cal,  u2_bar_true, (u2_bar_cal-u2_bar_true)/u2_bar_true)


for N in [8, 16, 32, 128]:
    compute_volume(N)


N:8 Area square (True): 0.9921460183660256, (Estimated) 0.9940558967731556 (Error) -0.0019249973005742438
N:8 Area circle (True): 0.007853981633974483, (Estimated) 0.005944103226844711 (Error) 0.24317327135934297
8 0.0009940558967731555 0.0009921460183660257 0.0019249973005739866
8 1.188820645368942e-06 1.5707963267948969e-06 -0.2431732713593432
N:16 Area square (True): 0.9921460183660256, (Estimated) 0.9940558967731558 (Error) -0.0019249973005743557
N:16 Area circle (True): 0.007853981633974483, (Estimated) 0.00594410322684471 (Error) 0.24317327135934308
16 0.000994055896773155 0.0009921460183660257 0.0019249973005735495
16 1.1888206453689421e-06 1.5707963267948969e-06 -0.24317327135934305
N:32 Area square (True): 0.9921460183660256, (Estimated) 0.9929289321881329 (Error) -0.0007891114892511149
N:32 Area circle (True): 0.007853981633974483, (Estimated) 0.0070710678118654745 (Error) 0.09968368384289407
32 0.000992928932188136 0.0009921460183660257 0.0007891114892541431
32 1.4142135623730988e-06 1.5707963267948969e-06 -0.09968368384289172
N:128 Area square (True): 0.9921460183660256, (Estimated) 0.9921963871193602 (Error) -5.076748019166827e-05
N:128 Area circle (True): 0.007853981633974483, (Estimated) 0.007803612880645131 (Error) 0.006413148855794224
128 0.0009921963871193467 0.0009921460183660257 5.076748017790442e-05
128 1.5607225761290172e-06 1.5707963267948969e-06 -0.006413148855800095
5

Is it possible to improve the accuracy but keep the resolution smaller?

6

Create a second order mesh. This is supported in DOLFINx
using for instance Gmsh for generating the mesh, as shown in dolfinx/demo_gmsh.py at main · FEniCS/dolfinx · GitHub
and in various tutorials at:
The FEniCSx tutorial — FEniCSx tutorial