Hi all,
I am testing dolfinx and try to convert some codes from dolfin to new version. I modified some of the tutorials (e.g. A nonlinear Poisson equation — FEniCSx tutorial) I want to solve the stationary heat equation - \nabla \cdot (\kappa \cdot \nabla T) = 0 with \kappa the conduction tensor and T the temperature. \kappa varies depending on the subdomain.
The error occurs in the line
n, converged = solver.solve(T)
with Did you forget to '#include <pybind11/stl.h>'? Or <pybind11/complex.h>, <pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic conversions are optional and require extra headers to be included when compiling your pybind11 module.
Here is my MWE:
''' Thermal Quasistatic Problem -> Test for implementing in FenicsX'''
# Packages
from dolfinx import fem, mesh, io, nls, log
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, TensorFunctionSpace, Function, dirichletbc, locate_dofs_geometrical, form
from dolfinx.io import XDMFFile
import numpy as np
from mpi4py import MPI
from ufl import grad, div, dx, ds, dot, inner, FacetNormal, Measure, TrialFunction, TestFunction, FiniteElement, lhs, rhs, derivative
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
## Mesh
meshFile = "RVE_2D_f03_Circle_fenics" # mesh file name
with XDMFFile(MPI.COMM_WORLD, meshFile+"_mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(domain, name="Grid") # cell tags
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, meshFile+"_facets.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(domain, name="Grid") # facet tags
# Different surfaces
dx = Measure('dx', domain=domain, subdomain_data=ct)
# number of cells
num_cells = domain.topology.index_map(domain.topology.dim).size_local
## Parameters
S = TensorFunctionSpace(domain, ("DG", 0))
kappa = Function(S)
block_size_s = S.dofmap.index_map_bs
for i in range(num_cells):
if ct.values[i] == 1: # matrix
kappa.x.array[[i*block_size_s, i*block_size_s+1, i*block_size_s+2, i*block_size_s+3]] = [2, 0, 0, 1]
elif ct.values[i] == 2: # inclusions
kappa.x.array[[i*block_size_s, i*block_size_s+1, i*block_size_s+2, i*block_size_s+3]] = [5, 0, 0, 5]
## Element, FunctionSpace and Functions
# Element
P1 = FiniteElement("CG", domain.ufl_cell(), 1)
element = P1
# FunctionSpace
V = FunctionSpace(domain, element)
#Test- und Trial-Functions
dfunc = TrialFunction(V)
T = Function(V)
dT = TestFunction(V)
## Boundary Conditions
T_left = 273.0 # K
T_right = 310.0 # K
def left(x):
return np.isclose(x[0], domain.geometry.x.min())
def right(x):
return np.isclose(x[0], domain.geometry.x.max())
fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 3), np.full_like(right_facets, 4)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(3))
right_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(4))
bc_left = [dirichletbc(ScalarType(T_left), left_dofs, V)]
bc_right = [dirichletbc(ScalarType(T_right), right_dofs, V)]
bcs = [bc_left, bc_right]
## variational problem
eqn = dot(dot(kappa, grad(T)), grad(dT)) * dx
Jac = derivative(eqn, T, dfunc)
problem = fem.petsc.NonlinearProblem(eqn, T, bcs=bcs, J = Jac)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(T)
assert(converged)
print(f"Number of interations: {n:d}")
I know that this is actually a linear problem. However, I just wanted to do it like that to have a working example for nonlinear cases. In addition, it also does not work with a linear solver.
Here iare the mesh files: Cloudstore
Thank you for your help.