Hello.
In the examples kindly provided (electro-magnetism and sub-domains), the material properties (constants) are set by sub-domain. I would like to know if it is possible to split the weak form in subdomains.
In the tutorial of subdomains, a variable \kappa{} is defined such that
\kappa=\begin{cases} 1 &\text{if } x\in\Omega_0\\ 0.1& \text{if } x\in\Omega_1\\ \end{cases}
kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=ScalarType) kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=ScalarType)
or
def eval_kappa(x): values = np.zeros(x.shape[1], dtype=ScalarType) # Create a boolean array indicating which dofs (corresponding to cell centers) # that are in each domain top_coords = x[1]>0.5 bottom_coords = x[1]<0.5 values[top_coords] = np.full(sum(top_coords), 0.1) values[bottom_coords] = np.full(sum(bottom_coords), 1) return values kappa2 = Function(Q) kappa2.interpolate(eval_kappa)
Would this scale for larger domains?
Does it mean that having kappa for a larger domain results in a much bigger memory allocation?, and that we would have to multiply that amount by the number of coefficients (kappa_1, kappa_2, etc.)?
Is it possible and convenient to split the formulation?
Instead of a constant, I would like to have a mathematical function which depends on the solution variable (u below)
-
The original idea is that
The whole domain will be \Omega=[0,1]\times[0,1], which consists of two subdomains \Omega_0=[0,1]\times [0,1/2] and \Omega_1=[0,1]\times[1/2, 1].
-
When defining the weak form in the example, it is still done in the full domain as
a = inner(kappa*grad(u), grad(v)) * dx -
Other examples (hyperelasticity, ns_code2–navier stokes, robin_neumann_dirichlet) apply different forms with subdomains, but they only do it on \partial{\Omega{}}
(from hyperelasticity)
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)(from ns~code2~)
dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle_marker)(from robin~neumanndirichlet~)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
Is something like this also convenient? is there loss of efficiency due to data access?
::: captioned-content
::: caption
Alternative definition of the weak form (with subdomains?)
:::
kappa = 1
dx0 = ufl.Measure("dx", domain=msh, subdomain_data=cells_0)
dx1 = ufl.Measure("dx", domain=msh, subdomain_data=cells_1)
a0 = inner(kappa * grad(u), grad(v)) * dx0
a1 = inner(0.1 * kappa * grad(u), grad(v)) * dx1
a = a0 + a1
:::
Minimal non-working example
# This code is based on the FEniCSx' subdomains tutorial
# from Jørgen S. Dokken. It was modified to portray the idea
# of subdomains in a weak form
from dolfinx import fem
from dolfinx import mesh
from dolfinx import plot as dlfplt
import dolfinx.fem.petsc as petsc
from dolfinx import default_scalar_type
import numpy as np
from mpi4py import MPI
import ufl
import pyvista
def Omega_0(x):
return x[1] <= 0.5
def Omega_1(x):
return x[1] >= 0.5
msh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
# # -:- skipped: for definition of kappa -:-
# Q = fem.functionspace(msh, ("DG", 0))
cells_0 = mesh.locate_entities(msh, msh.topology.dim, Omega_0)
cells_1 = mesh.locate_entities(msh, msh.topology.dim, Omega_1)
# # -:- skipped: kappa definition -:-
# kappa = fem.Function(Q)
# kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
# kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=default_scalar_type)
kappa = 1
V = fem.functionspace(msh, ("Lagrange", 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
# # -:- skipped: redifinition of weak form -:-
# a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
# L = fem.Constant(msh, default_scalar_type(1)) * v * ufl.dx
# Alternative definition of the weak form
kappa = default_scalar_type(0.1)
dx0 = ufl.Measure("dx", domain=msh, subdomain_data=cells_0)
dx1 = ufl.Measure("dx", domain=msh, subdomain_data=cells_1)
a0 = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx0
a1 = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx1
a = a0 + kappa * a1
L = (fem.Constant(msh, default_scalar_type(1)) * v * dx0
+ fem.Constant(msh, default_scalar_type(1)) * v * dx1)
# x = SpatialCoordinate(msh)
dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
bcs = [fem.dirichletbc(default_scalar_type(1), dofs, V)]
problem = petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
[EDIT]
My system
-
FEniCSx software
dolfinx: 0.8.0.dev0_r27624.095ef96-1
basix: 0.8.0.dev0_r975.15b915b-1
ufl: 2023.3.0.dev0_r3583.e8d3077-1
ffcx: 0.8.0.dev0_r7098.833967c-1 -
Dependencies
python: 3.11.5-2
python-numpy: 1.26.0-1
petsc: 3.19.5.1172.g92f94a14170-1
hdf5-openmpi: 1.14.2-1
boost: 1.83.0-2
adios2: 2.8.3-5
scotch: 7.0.4-1
pybind11: 2.11.1-1
python-build: 1.0.1-1
python-cffi: 1.15.1-4
python-cppimport: 22.08.02.r6.g0849d17-1
python-imageio: 2.33.0-2
python-imageio-ffmpeg: 0.4.9-1
python-installer: 0.7.0-3
python-meshio: 5.3.4-2
python-mpi4py: 3.1.4-3
python-pytest: 7.4.2-1
python-scikit-build: 0.17.6-1
python-setuptools: 1:68.0.0-1
python-wheel: 0.40.0-3
xtensor: 0.24.0-2
xtensor-blas: 0.20.0-1
python-cattrs: 23.1.2-1
python-scikit-build: 0.17.6-1
nanobind: 1.8.0-1 -
Operating system
Arch Linux: 6.5.4-arch2-1
Disclaimer
I have been working on this for a while, and missed these which look similar:
- Defining subdomains for different materials, but with coefficients that aren't constant (without MWE)
- Split a 2d solution to add different constants to each component and write to file (legacy DOLFIN)
- Difference between split() and sub() in mixed formulation (for mixed formulation; related?)
- How to assign different values to different elements and How to extract values of a variable defined as a function of the solution u? (for mixed formulation; related?)
