Custom integrals and source terms

General
1

Hi, I would like advice from the community regarding non trivial source terms.
I used info from this topic for the integrations.

The problem is the following: consider the 2D problem of the known FEniCSx tutorial:
- {\nabla ^2}u(x,y) = f(x,y) where x \in [0.1] and y \in [0.1] with \Omega = [0,1] \times [0,1], f(x,y)=-6 and Dirichlet BCs: u_{D}(x,y)=1+x^2+2y^2. With the known definitions for the trial and test functions the variational problem is \int_\Omega {{\rm{d}}x} \nabla u \cdot \nabla v = \int_\Omega {{\rm{d}}x} fv and the analytical solution is u(x,y)=1+x^2+2y^2.

Consider now that the problem changes with the addition of various non trivial source terms.

Firstly, say f(x,y)=-6+\theta(0.5-x)\theta(0.5-y). The new variational problem is:

\int_\Omega {{\rm{d}}x} \nabla u \cdot \nabla v = \int_\Omega {{\rm{d}}x} fv + \int_\Omega {{\rm{d}}x} \theta (0.5 - x)\theta (0.5 - y)v

How can I define the proper measure? Since this is a integration of the same spatial dimension as the problem I tried (full MWE in the end)

def facet_eq_2D(x):
   tol=1e-6
   less1 =  x[0] < 0.5 + tol
   less2 = x[1] < 0.5 + tol
   return np.logical_and(less1, less2)
cells = dolfinx.mesh.locate_entities(domain, domain.topology.dim, facet_eq_2D)
tags = np.full_like(cells, 1)
cell_tags = dolfinx.mesh.meshtags(domain, domain.topology.dim, cells, tags)
dx1 = ufl.Measure('dx', domain, subdomain_data=cell_tags, subdomain_id=1 )
print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(uh * dx1(1))))
L1 = source * v * ufl.dx + v * dx1(1)
problem = dolfinx.fem.petsc.LinearProblem(a, L1, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh1 = problem. Solve()

but it doesn’t work.
Still though the line print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(uh * dx1(1)))) very accurately calculates the integral.

Same for the case f(x,y)=-6+\delta(x-0.5)\theta(0.5-y) for the source the following code fails. This is the case of a line integral so I tried

def facet_eq_1D(x):
    tol=1e-6
    close =  np.isclose(x[0], 0.5)
    less = x[1] < 0.5 + tol
    return np.logical_and(close, less)
cells = dolfinx.mesh.locate_entities(domain, domain.topology.dim-1, facet_eq_1D)
tags = np.full_like(cells, 1)
cell_tags = dolfinx.mesh.meshtags(domain, domain.topology.dim-1, cells, tags)
dx2 = ufl.Measure("dS", subdomain_data=cell_tags, subdomain_id=1 )
print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(uh * dx2(1))))
L2 = source * v * ufl.dx + v * dx2(1)
problem = dolfinx.fem.petsc.LinearProblem(a, L1, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

with no success.

Any help on how to define and solve those variational problem with FEniCSx is greatly appreciated.

MWE

from mpi4py import MPI
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
import dolfinx

domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.FunctionSpace(domain, ("CG", 2))
uD = dolfinx.fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, domain.topology.dim-1, boundary_facets)
bc = dolfinx.fem.dirichletbc(uD, boundary_dofs)
source = dolfinx.fem.Constant(domain, ScalarType(-6))
#Solution of initial problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = source * v * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

def facet_eq_2D(x):
   tol=1e-6
   less1 =  x[0] < 0.5 + tol
   less2 = x[1] < 0.5 + tol
   return np.logical_and(less1, less2)
cells = dolfinx.mesh.locate_entities(domain, domain.topology.dim, facet_eq_2D)
tags = np.full_like(cells, 1)
cell_tags = dolfinx.mesh.meshtags(domain, domain.topology.dim, cells, tags)
dx1 = ufl.Measure('dx', domain, subdomain_data=cell_tags, subdomain_id=1 )
print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(uh * dx1(1))))
L1 = source * v * ufl.dx + v * dx1(1)
problem = dolfinx.fem.petsc.LinearProblem(a, L1, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh1 = problem.solve()


def facet_eq_1D(x):
    tol=1e-6
    close =  np.isclose(x[0], 0.5)
    less = x[1] < 0.5 + tol
    return np.logical_and(close, less)
cells = dolfinx.mesh.locate_entities(domain, domain.topology.dim-1, facet_eq_1D)
tags = np.full_like(cells, 1)
cell_tags = dolfinx.mesh.meshtags(domain, domain.topology.dim-1, cells, tags)
dx2 = ufl.Measure("dS", subdomain_data=cell_tags, subdomain_id=1 )
print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(uh * dx2(1))))
L2 = source * v * ufl.dx + v * dx2(1)
problem = dolfinx.fem.petsc.LinearProblem(a, L1, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh2 = problem.solve()
2

I don’t understand what you mean by measure here.
As far as I can tell, you want to integrate the function \theta(0.5-x)\theta(0.5-y) over the whole domain \Omega. This does not require a new measure.

Your definition of dx1 generates an integral measure over all cells whose vertices satisfies

which would be an integral over all cells where x<0.5+tol and y<0.5+tol in your domain.

As you are stating that this is a 3D problem, you are now searching for all facets where both vertices has an x coordinate close to 0.5 (which cannot be satisfied). If you want to implement the point-source at a given degree of freedom, you should locate that degree of freedom and add to the correct position in the rhs. If you want to use a smoothed version of delta, see for instance; Dirac delta distribution (dolfinx) - #3 by dokken

3

Thanks for taking the time to write this.

In my original question I forgot to write that \theta is the Heaviside step function (usually defined as H(x) and less commonly as \theta(x)). So in my first problem the product of the step functions restricts the extra source in the area where x,y \in (0,0.5) so I was wondering if there is a more elegant way to write the variational problem by changing the integration measure thus avoiding analytical approximations for the step function.

The idea was to try to do something similar with the delta function. I have actually tried the smooth analytical approximation of delta with a gaussian and it works in this simple problem but in more complex ones the results are highly dependent on the mesh and the standard deviation of the gaussian which is understandable.