Thank you very much @dokken
However, I have one question, when i add the extra term Constant(0)*dx, i get an error saying,
This integral is missing an integration domain.
Traceback (most recent call last):
L2 = dot(grad(f)('-'), n('-'))*dS(10) + Constant(0)*dx
~~~~~~~~~~~^~~
File "/home/dputluru/.conda/envs/cent7/2020.11-py38/fenicsproject/lib/python3.12/site-packages/ufl/measure.py", line 441, in __rmul__
error("This integral is missing an integration domain.")
File "/home/dputluru/.conda/envs/cent7/2020.11-py38/fenicsproject/lib/python3.12/site-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: This integral is missing an integration domain.
No clue why it occured.
I used to define my dx as,
domains = MeshFunction("size_t", mesh, nsd)
dx = dx(metadata={'quadrature_degree': 2}, subdomain_data=domains)
hdf.read(domains, "/domains")
Now when i add
dx = Measure('dx', domain=mesh, subdomain_data=domains), it runs fine with sensible results.
Here is minimal working code if required,
from dolfin import *
import numpy as np
mesh = Mesh()
h5_file = "dS_mesh.h5"
hdf = HDF5File(mesh.mpi_comm(), h5_file, "r")
hdf.read(mesh, "/mesh", False)
# define mesh parameters
nsd = mesh.topology().dim()
n = FacetNormal(mesh)
I = Identity(nsd)
h = CellDiameter(mesh)
x_i = SpatialCoordinate(mesh)
# read mesh functions
boundaries = MeshFunction("size_t", mesh, nsd - 1)
ds = ds(metadata={'quadrature_degree': 2}, subdomain_data=boundaries)
dS = dS(metadata={'quadrature_degree': 2}, subdomain_data=boundaries)
hdf.read(boundaries, "/boundaries")
domains = MeshFunction("size_t", mesh, nsd)
dx = dx(metadata={'quadrature_degree': 2}, subdomain_data=domains)
hdf.read(domains, "/domains")
# start ----- new lines that made it work
dx = Measure('dx', domain=mesh, subdomain_data=domains)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries, subdomain_id=10)
# end ----- new lines that made it work
# output paraview file for domains
output_file = File("domains.pvd")
output_file << domains
V = FunctionSpace(mesh, 'CG', 1)
# Make up something with discontinuous gradient
# Define the level set function for the sphere (x^2 + y^2 + z^2 - 25 = 0)
class SphereLevelSet(UserExpression):
def eval(self, value, x):
r_sqr = x[0]**2 + x[1]**2 + x[2]**2
value[0] = r_sqr - 25
def value_shape(self):
return ()
# Create the level set function
phi = SphereLevelSet(degree=1)
# Define the expression with different values inside and outside the sphere
class MyExpression(UserExpression):
def __init__(self, level_set, **kwargs):
super().__init__(**kwargs)
self.level_set = level_set
def eval(self, value, x):
if self.level_set(x) <= 0: # Inside the sphere
value[0] = 0
else: # Outside the sphere
# Linear increase proportional to the distance from the sphere's surface
distance_from_surface = np.sqrt(self.level_set(x) + 25) - 5
value[0] = distance_from_surface
def value_shape(self):
return ()
# Create and interpolate the expression
my_expression = MyExpression(phi, degree=1)
f = interpolate(my_expression, V)
n = FacetNormal(mesh)
# +, - sides are by default assigned arbitrarily so assembling below
# will not yield meaningful result
L0 = dot(grad(f)('+'), n('+'))*dS
print(assemble(L0))
L1 = dot(grad(f)('-'), n('-'))*dS(10)
print(assemble(L1))
# Add dummy terms that will allow for controlling restrictions
# This is (0, 0) * (1, 0) * 2
L2 = dot(grad(f)('-'), n('-'))*dS(10) + Constant(0)*dx
print(assemble(L2))
# This is (1, 0) . (-1, 0) * 2
L3 = dot(grad(f)('+'), n('+'))*dS(10) + Constant(0)*dx
print(assemble(L3))
“I think the issue is with the way i import my mesh. Hence i cannot use mshr package.”
I am unable to upload the mesh file here 