Hello everyone,
I am trying to create a unit circular domain with subdomains of radius 0.6 and radius 0.8.
I have tried to create a circular mesh and I want to solve
\int_\Omega\kappa \nabla u \nabla v + \int_{\partial \Omega} \epsilon uv = \int_{\Omega_{0}} f v
\kappa = \left\{ \begin{array}{} 1, x \in \Omega_{0} \cup \Omega_{2} \\ 1/50, x \in \Omega_{1} \end{array} \right\}
My code is as follows:
import gmsh
import numpy as np
import pyvista
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, compute_midpoints
from dolfinx.plot import create_vtk_mesh
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
dx, grad, inner)
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import dolfinx
import matplotlib.pyplot as plt
import ufl
from dolfinx import mesh
gmsh.initialize()
c2 = gmsh.model.occ.addDisk(0, 0, 0, 1,1)
gmsh.model.occ.synchronize()
gdim = 2
status2 = gmsh.model.addPhysicalGroup(gdim, [c2], 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.001)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.1)
gmsh.model.mesh.generate(gdim)
from dolfinx import io
if MPI.COMM_WORLD.rank == 0:
# Get mesh geometry
geometry_data = io.extract_gmsh_geometry(gmsh.model)
# Get mesh topology for each element
topology_data = io.extract_gmsh_topology_and_markers(gmsh.model)
if MPI.COMM_WORLD.rank == 0:
# Extract the cell type and number of nodes per cell and broadcast
# it to the other processors
gmsh_cell_type = list(topology_data.keys())[0]
properties = gmsh.model.mesh.getElementProperties(gmsh_cell_type)
name, dim, order, num_nodes, local_coords, _ = properties
cells = topology_data[gmsh_cell_type]["topology"]
cell_id, num_nodes = MPI.COMM_WORLD.bcast([gmsh_cell_type, num_nodes], root=0)
else:
cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
cells, geometry_data = np.empty([0, num_nodes]), np.empty([0, gdim])
from dolfinx import mesh, cpp
# Permute topology data from MSH-ordering to dolfinx-ordering
ufl_domain = io.ufl_mesh_from_gmsh(cell_id, gdim)
gmsh_cell_perm = io.cell_perm_gmsh(cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]
# Create distributed mesh
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, geometry_data[:, :gdim], ufl_domain)
from dolfinx import fem
Q = fem.FunctionSpace(domain, ("DG", 0))
V = fem.FunctionSpace(domain, ("CG", 1))
def Omega_0(x):
return (x[0]*x[0] + x[1]*x[1] <= 0.6*0.6)
def Omega_1(x):
return [x and y for (x,y) in zip(x[0]*x[0] + x[1]*x[1] >= 0.6*0.6, x[0]*x[0] + x[1]*x[1] <= 0.8*0.8)]
def Omega_2(x):
return (x[0]*x[0] + x[1]*x[1] >= 0.8*0.8)
kappa = Function(Q)
cells_0 = locate_entities(domain, domain.topology.dim, Omega_0)
cells_1 = locate_entities(domain, domain.topology.dim, Omega_1)
cells_2 = locate_entities(domain, domain.topology.dim, Omega_2)
kappa.x.array[:] = 1
kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=ScalarType)
kappa.x.array[cells_1] = np.full_like(cells_1, 1/50, dtype=ScalarType)
kappa.x.array[cells_2] = np.full_like(cells_2, 1, dtype=ScalarType)
from petsc4py.PETSc import ScalarType
x = ufl.SpatialCoordinate(domain)
# cell tag
tdim = domain.topology.dim
fdim = tdim - 1
num_cells = domain.topology.index_map(tdim).size_local + \
domain.topology.index_map(tdim).num_ghosts
marker1 = lambda x: x[0]*x[0] + x[1]*x[1] <= 0.6*0.6
indices1 = mesh.locate_entities(domain, 2, marker1)
cell_marker = np.zeros(num_cells, dtype=np.int32)
cell_marker[indices1] = 1
cell_tag = dolfinx.mesh.meshtags(domain, dim=2, indices=np.arange(num_cells),values=cell_marker)
dx_sub = ufl.Measure('dx', subdomain_data=cell_tag, domain=domain)
Vs = dolfinx.fem.FunctionSpace(domain,("CG",1))
uDs = dolfinx.fem.Function(Vs)
uDs.interpolate(lambda x:x[0]**2 + x[1]**2)
us,vs = ufl.TrialFunction(Vs), ufl.TestFunction(Vs)
domain.topology.create_connectivity(tdim-1, tdim)
boundary_facets = np.flatnonzero(dolfinx.mesh.compute_boundary_facets(domain.topology)) # substitutive: facet_marker + np.where
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, tdim-1, boundary_facets)
bcsd = dolfinx.fem.dirichletbc(uDs, boundary_dofs)
bcs = [bcsd]
import numpy
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = numpy.flatnonzero(mesh.compute_boundary_facets(domain.topology))
def on_boundary(x):
return np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1)
boundary_dofs = fem.locate_dofs_geometrical(V, on_boundary)
bc = fem.dirichletbc(ScalarType(0), boundary_dofs, V)
import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
from petsc4py.PETSc import ScalarType
#f = fem.Constant(domain, ScalarType(-6))
f = fem.Function(Vs)
f.interpolate(lambda x: x[0]**2-x[1]**2)
eps = 1e-3
ae = ufl.inner(kappa*ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.inner(eps*u, v) * ufl.ds
Le = ufl.inner(f, vs)* dx_sub
problem = fem.petsc.LinearProblem(ae, Le, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
ue = problem.solve()
!pip3 install pyvista pythreejs ipygany --upgrade
!jupyter nbextension enable --py --sys-prefix ipygany
from dolfinx import plot
import pyvista
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
pyvista.set_jupyter_backend("pythreejs")
p = pyvista.Plotter(window_size=[800, 800], shape=(1,2))
# Filter out ghosted cells
num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
marker = np.zeros(num_cells_local, dtype=np.int32)
cells_0 = cells_0[cells_0<=num_cells_local]
cells_1 = cells_1[cells_1<=num_cells_local]
cells_2 = cells_2[cells_2<=num_cells_local]
marker[cells_0] = 1
marker[cells_1] = 2
marker[cells_2] = 3
topology, cell_types, x = create_vtk_mesh(domain, domain.topology.dim, np.arange(num_cells_local, dtype=np.int32))
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
grid.cell_data["Marker"] = marker
grid.set_active_scalars("Marker")
p.subplot(0,0)
actor0 = p.add_mesh(grid, show_edges=True)
p.subplot(0,1)
grid_uh = pyvista.UnstructuredGrid(*create_vtk_mesh(V))
grid_uh.point_data["u"] = ue.x.array.real
grid_uh.set_active_scalars("u")
actor1 = p.add_mesh(grid_uh, show_edges=True)
if not pyvista.OFF_SCREEN:
p.show()
else:
pyvista.start_xvfb()
figure = p.screenshot("subdomains_structured.png")
But I got the output as follows where I am getting a gap between the subdomains :
Kindly let me know why am I getting the gap and the mistake with this code.
I am new to dolfinx and I am learning it slowly, It would be great if anyone could help me with this. Thank you in advance!