Applying a spatially variable Dirichlet b.c

Hello, I have a code that doesn’t work (returns an error) when I try to apply a spatially variable Dirichlet boundary condition, that I am unable to resolve.

I tried to create a MWE, except that this MWE works, i.e., does not display my error. The only major difference between the MWE and my code, aside a totally different paradigm (OO vs straight sequential code) is that I use a mesh file in my code and a dolfinx generated code in the MWE.

The relevant part in my code that yields an error is:

        # Define the boundary conditions
        left_facets = facet_markers.find(inj_current_curve)
        right_facets = facet_markers.find(out_current_curve)
        left_dofs = locate_dofs_topological(ME.sub(1), mesh.topology.dim-1, left_facets)

        left_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, left_facets)
        right_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, right_facets)

        #bc_temp_left = dirichletbc(ScalarType(in_params.T_cold), left_dofs_temp, ME.sub(0))
        #bc_temp_right = dirichletbc(ScalarType(in_params.T_cold + in_params.ΔT), right_dofs_temp, ME.sub(0))
        
        def T_left_func(x): 
            return np.full(x.shape[1], 250*x[1])
        T0 = ME.sub(0)
        Q, _ = T0.collapse()
        T_left = Function(Q)
        T_left.interpolate(T_left_func)
        
        T_right=T_left

        # Apply Dirichlet boundary conditions with position-dependent functions
        bc_temp_left = dirichletbc(T_left, left_dofs_temp, T0)
        bc_temp_right = dirichletbc(T_right, right_dofs_temp, T0)        
        
        bcs = [bc_temp_left, bc_temp_right]

with traceback:

Info    : Reading 'meshes/quarter_annulus.msh'...
Info    : 10 entities
Info    : 10234 nodes
Info    : 20466 elements
Info    : Done reading 'meshes/quarter_annulus.msh'
Started simulation
debug curves: 12, 13
Traceback (most recent call last):
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py", line 175, in dirichletbc
    bc = bctype(_value, dofs, V)
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace<double>)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace<double>)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double, double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double, double>, dofs: Annotated[List[numpy.ndarray[numpy.int32]], FixedSize(2)], V: dolfinx::fem::FunctionSpace<double>)

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7a351c43d930>, array([    1,    20,    23,    24,    74,    77,    78,   152,   155,
         156,   248,   251,   252,   374,   377,   378,   524,   527,
         528,   710,   713,   714,   926,   929,   930,  1178,  1181,
        1182,  1436,  1439,  1440,  1724,  1727,  1728,  2036,  2039,
        2040,  2396,  2399,  2400,  2780,  2783,  2784,  3182,  3185,
        3186,  3614,  3617,  3618,  4070,  4073,  4074,  4556,  4559,
        4560,  5060,  5063,  5064,  5606,  5609,  5610,  6158,  6161,
        6162,  6746,  6749,  6750,  7358,  7361,  7362,  7988,  7991,
        7992,  8642,  8645,  8646,  9326,  9329,  9330, 10028, 10031,
       10032, 10760, 10763, 10764, 11516, 11519, 11520, 12290, 12293,
       12294, 13094, 13097, 13098, 13922, 13925, 13926, 14456, 14459,
       14460, 14468, 14474, 14475, 14834, 14837, 14838, 14846, 14852,
       14853, 15416, 15419, 15420, 15428, 15434, 15435, 15824, 15827,
       15828, 16418, 16421, 16422, 16844, 16847, 16848, 17450, 17453,
       17454, 17894, 17897, 17898, 17906, 17912, 17913, 18536, 18539,
       18540, 19016, 19019, 19020, 19664, 19667, 19668], dtype=int32), FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, equispaced, unset, False), (2,)), 0), FiniteElement('Lagrange', triangle, 3))

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/root/shared/debug_quarter_annulus_renormalize_PF/Bridgman_Nye_example.py", line 477, in <module>
    post_proc = PostProcessor()#thermoelec_sim, display_output=True, save_J_e=True)
  File "/root/shared/debug_quarter_annulus_renormalize_PF/Bridgman_Nye_example.py", line 252, in __init__
    self.thermoelec_sim_res = thermoelec_sim.run_sim()
  File "/root/shared/debug_quarter_annulus_renormalize_PF/Bridgman_Nye_example.py", line 163, in run_sim
    bc_temp_left = dirichletbc(T_left, left_dofs_temp, T0)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py", line 177, in dirichletbc
    bc = bctype(_value, dofs, V._cpp_object)
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace<double>)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace<double>)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double, double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double, double>, dofs: Annotated[List[numpy.ndarray[numpy.int32]], FixedSize(2)], V: dolfinx::fem::FunctionSpace<double>)

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7a351c43d930>, array([    1,    20,    23,    24,    74,    77,    78,   152,   155,
         156,   248,   251,   252,   374,   377,   378,   524,   527,
         528,   710,   713,   714,   926,   929,   930,  1178,  1181,
        1182,  1436,  1439,  1440,  1724,  1727,  1728,  2036,  2039,
        2040,  2396,  2399,  2400,  2780,  2783,  2784,  3182,  3185,
        3186,  3614,  3617,  3618,  4070,  4073,  4074,  4556,  4559,
        4560,  5060,  5063,  5064,  5606,  5609,  5610,  6158,  6161,
        6162,  6746,  6749,  6750,  7358,  7361,  7362,  7988,  7991,
        7992,  8642,  8645,  8646,  9326,  9329,  9330, 10028, 10031,
       10032, 10760, 10763, 10764, 11516, 11519, 11520, 12290, 12293,
       12294, 13094, 13097, 13098, 13922, 13925, 13926, 14456, 14459,
       14460, 14468, 14474, 14475, 14834, 14837, 14838, 14846, 14852,
       14853, 15416, 15419, 15420, 15428, 15434, 15435, 15824, 15827,
       15828, 16418, 16421, 16422, 16844, 16847, 16848, 17450, 17453,
       17454, 17894, 17897, 17898, 17906, 17912, 17913, 18536, 18539,
       18540, 19016, 19019, 19020, 19664, 19667, 19668], dtype=int32), <dolfinx.cpp.fem.FunctionSpace_float64 object at 0x7a351c760330>

The MWE that I attempt to follow (because it works!) is:

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
                         locate_dofs_topological, Expression, create_matrix, assemble_matrix,
                         assemble_vector, apply_lifting)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags, compute_incident_entities
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8    
    
domain = mesh.create_unit_square(MPI.COMM_WORLD, 40, 40, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
    return np.isclose(x[0], 0)

def boundary_D_bottom(x):
    return np.isclose(x[1], 0)

def T_cold_func(x): 
    return np.full(x.shape[1], 250*x[1])
def T_hot_func(x): 
    return np.full(x.shape[1], 200.0)

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])
the_current = 0.0

# Define the boundary conditions
T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_bottom)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T1 = ME.sub(0)
Q1, _ = T1.collapse()
T_cold = Function(Q1)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T1, Q1), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T1)

bcs = [bc_temp_left, bc_temp_right]

Any idea what’s wrong with my code? Compared to the MWE, I do not see what’s the difference.

This is not quite right when you have a collapsed sup space.
You need to do the following

then

left_dofs_temp = locate_dofs_topological((T0,Q), mesh.topology.dim-1, left_facets)

As this yields a map between the dofs in the collapsed space and the dofs in the sub space (T0).

Then this should work

1 Like