Hi all,
I am very new to Fenics and I am hoping to find some help to setup my first 3D vector problem.
For simplicity, I simplified the problem as much as possible and designed the following toy example:
\begin{aligned}
\boldsymbol{u}\left( \boldsymbol{x} \right) = \boldsymbol{f}\left( \boldsymbol{x} \right) \\
\boldsymbol{f}\left( \boldsymbol{x} \right) = \hat{\phi}
\end{aligned}
where $$ \hat{\phi} $$ is the azimuthal unit vector in a cylindrical coordinate system.
My geometry is a torus, centered in the origin: my forcing term is then a constant function flowing orthogonal to the cross section of the torus.
What am I doing wrong in the code below? I suspect I am not assembling correctly the boundary conditions, because if I change problem = LinearProblem(a, L, bcs=[bcs]) to problem = LinearProblem(a, L) the program runs. In this case the b.c. are superfluous, but the example is didactic to understand how to set them (e.g., if I enclose the torus inside a sphere and I want to apply Dirichlet b.c. on the boundary of the sphere).
I apologize if the code is a bit long, but the first part is only to generate the mesh (included for reproducibility), and the interesting portion is exclusively in the bottom half.
I am running the snipped below inside the docker image dokken92/dolfinx_custom:15072022 (DOLFINx version: 0.4.2.0).
Thanks in advance!
import gmsh
import numpy as np
from mpi4py import MPI
import os
from dolfinx.io import gmshio
import pyvista
from dolfinx.plot import create_vtk_mesh
_OUTPUT_PATH = "/root/shared/out"
geometry_name = "torus"
a = 0.1 # Radius of torus cross section [m]
r_torus = 1 # Radius of torus [m]
gdim = 3 # Geometric dimension of the mesh
mesh_density = 50E-3 # Mesh density [m]
rank = MPI.COMM_WORLD.rank
gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 0)
model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:
# Define torus geometry
top_cross_section_surface = gmsh.model.occ.addDisk(r_torus, 0, 0, a, a, zAxis = [0, 1, 0])
gmsh.model.occ.synchronize()
c1_plus = gmsh.model.occ.revolve([(2, top_cross_section_surface)], 0,0,0, 0,0,1, np.pi)
c1_minus = gmsh.model.occ.revolve([(2, top_cross_section_surface)], 0,0,0, 0,0,1, -np.pi)
c1_vol = [v for v in c1_plus + c1_minus if v[0] == 3]
# glue volumes, set mesh size and sync to model
c1 = gmsh.model.occ.fragment(c1_vol, [])
gmsh.model.occ.synchronize()
groups = []
groups.append(gmsh.model.addPhysicalGroup(gdim, [v[1] for v in c1[0]], name = f"Group {len(groups)+1}"))
gmsh.model.geo.synchronize()
vGroups = gmsh.model.getPhysicalGroups() # (dim, group_id)
gmsh.model.occ.mesh.setSize(gmsh.model.occ.getEntities(0), mesh_density) # Only entities of dimension 0 (points) are handled by OCC setSize.
gmsh.model.occ.synchronize()
gmsh.model.geo.synchronize()
# Generate mesh
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.option.setNumber("Mesh.Algorithm3D", 4)
gmsh.model.mesh.generate(gdim)
meshfile = os.path.join(_OUTPUT_PATH, f'{geometry_name}.msh')
gmsh.write(meshfile)
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()
from dolfinx.fem import (dirichletbc, Expression, Function, FunctionSpace,
VectorFunctionSpace, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary
from ufl import FiniteElement, SpatialCoordinate, TestFunction, TrialFunction, as_vector, dot, dx, grad, inner, Identity
from petsc4py.PETSc import ScalarType
V = FunctionSpace(mesh, ("N1curl", 1))
x = SpatialCoordinate(mesh)
# Compute phi_hat
def phi_hat_expression(x):
vals = np.zeros((mesh.geometry.dim, x.shape[1]))
phi = np.arctan2(x[1], x[0])
vals[0] = - np.sin(phi)
vals[1] = np.cos(phi)
return vals
phi_hat = Function(V)
phi_hat.interpolate(phi_hat_expression)
v = TestFunction(V)
L = inner(phi_hat, v) * dx
u = TrialFunction(V)
a = inner(u, v) * dx
boundary_facets = locate_entities_boundary(
mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True, dtype=bool))
boundary_dofs = locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)
zero_u = Function(V)
zero_u.x.array[:] = 0.0
bcs = [dirichletbc(zero_u, boundary_dofs)]
problem = LinearProblem(a, L, bcs=[bcs])
uh = problem.solve()
This is the error I am getting, which I can’t understand:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [51], in <cell line: 41>()
37 # A = Function(V)
38 # problem = LinearProblem(a, L, u=A, bcs=[bcs])
39 # problem.solve()
40 problem = LinearProblem(a, L, bcs=[bcs])
---> 41 uh = problem.solve()
File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:604, in LinearProblem.solve(self)
602 # Assemble lhs
603 self._A.zeroEntries()
--> 604 _assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
605 self._A.assemble()
607 # Assemble rhs
File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:362, in _assemble_matrix_mat(A, a, bcs, diagonal, constants, coeffs)
360 constants = _pack_constants(a) if constants is None else constants
361 coeffs = _pack_coefficients(a) if coeffs is None else coeffs
--> 362 _cpp.fem.petsc.assemble_matrix(A, a, constants, coeffs, bcs)
363 if a.function_spaces[0] is a.function_spaces[1]:
364 A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
TypeError: assemble_matrix(): incompatible function arguments. The following argument types are supported:
1. (A: mat, a: dolfinx::fem::Form<double>, constants: numpy.ndarray[numpy.float64], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]], bcs: List[dolfinx::fem::DirichletBC<double>], unrolled: bool = False) -> None
2. (A: mat, a: dolfinx::fem::Form<double>, constants: numpy.ndarray[numpy.float64], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]], rows0: numpy.ndarray[numpy.int8], rows1: numpy.ndarray[numpy.int8], unrolled: bool = False) -> None
Invoked with: <petsc4py.PETSc.Mat object at 0x7ff75f8c3290>, <dolfinx.fem.forms.Form object at 0x7ff77599dda0>, array([], dtype=float64), {(<IntegralType.cell: 0>, -1): array([], shape=(0, 0), dtype=float64)}, [[<dolfinx.fem.bcs.DirichletBC object at 0x7ff75f8c3880>]]