# How to define a function defined on boundary? **URL:** \ **Category:** General\ **Created:** [July 22, 2024, 2:26am UTC](https://fenicsproject.discourse.group/t/how-to-define-a-function-defined-on-boundary/15306 "2024-07-22T02:26:16Z")\ **Posts on this page:** 6\ **Page:** 1
**Author:** ![OsbertWang](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/osbertwang/32/2138_2.png) [@OsbertWang](https://fenicsproject.discourse.group/u/OsbertWang)\ **Post date:** [July 22, 2024, 2:26am UTC](https://fenicsproject.discourse.group/t/how-to-define-a-function-defined-on-boundary/15306/1 "2024-07-22T02:26:16Z")
Hello anyone, I once again do a numerical experiment based on [My previous question](https://fenicsproject.discourse.group/t/how-to-solve-a-heat-equation-whose-analysis-solution-is-known-and-expressed-by-spatialcoordinate/15222). The problem is: find u = u(\boldsymbol{x}, t) obeying \begin{alignat}{2} \frac{\partial u}{\partial t} - \nabla \cdot (\kappa \nabla u) &= f, && (\boldsymbol{x}, t) \in \Omega \times [0,T],\\ -\kappa \frac{\partial u}{\partial \boldsymbol{n}} &= g, &\quad& (\boldsymbol{x}, t) \in \partial \Omega \times [0,T],\\ u(0) &= u\_0, && \boldsymbol{x} \in \Omega, \end{alignat} where \Omega is unit square, T =2, and the diffusion coefficient \begin{equation} \kappa = \kappa (t). \end{equation} The sources \begin{align} f &= 2t - \kappa (6x + 12y),\\ g &= \begin{cases} 0, & y=0,\\ -6 \kappa, & y=1,\\ 0, & x=0,\\ -3 \kappa. & x=1, \end{cases}\\ u\_0 &= 1 + x^3 + 2y^3 \end{align} are derived from exact solution \begin{equation} u\_{ex} = 1 + x^3 + 2y^3 + t^2. \end{equation} The diffusion coefficient \begin{equation} \kappa = \sqrt{t} \end{equation} in this demo. I follow [this tutorial](https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html) to code as: ```auto from petsc4py import PETSc from mpi4py import MPI import ufl from dolfinx import mesh, fem, default_scalar_type from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_matrix, create_vector, set_bc import numpy as np import time start = time.time() t = 0.0 T = 2.0 num_steps = 20 dt = (T - t) / num_steps nx, ny = 5, 5 domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle) V = fem.functionspace(domain, ("Lagrange", 1)) class exact_solution(): def __init__ (self, t): self.t = t def __call__ (self, x): return 1 + x[0] **3 + 2 * x[1]** 3 + self.t ** 2 u_exact = exact_solution(t) u_n = fem.Function(V) u_n.interpolate(u_exact) kappa_expression = np.sqrt(t) kappa = fem.Constant(domain, kappa_expression) class source_term(): def __init__ (self, t, kappa): self.t = t self.kappa = kappa def __call__ (self, x): return 2 * self.t - self.kappa * (6 * x[0] + 12 * x[1]) f_expression = source_term(t, kappa_expression) f = fem.Function(V) f.interpolate(f_expression) g = fem.Function(V) boundaries = [ (1, lambda x: np.isclose(x[1], 0)), (2, lambda x: np.isclose(x[1], 1)), (3, lambda x: np.isclose(x[0], 0)), (4, lambda x: np.isclose(x[0], 1)) ] facet_indices, facet_markers = [], [] fdim = domain.topology.dim - 1 for (marker, locator) in boundaries: facets = mesh.locate_entities(domain, fdim, locator) if locator == 1: g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type) elif locator == 2: g.x.array[facets] = np.full_like(facets, -6.0 * kappa_expression, dtype = default_scalar_type) elif locator == 3: g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type) elif locator == 4: g.x.array[facets] = np.full_like(facets, -3.0 * kappa_expression, dtype = default_scalar_type) else: g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type) facet_indices.append(facets) facet_markers.append(np.full_like(facets, marker)) facet_indices = np.hstack(facet_indices).astype(np.int32) facet_markers = np.hstack(facet_markers).astype(np.int32) sorted_facets = np.argsort(facet_indices) facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]) ds = ufl.Measure("ds", domain = domain, subdomain_data = facet_tag) u, v = ufl.TrialFunction(V), ufl.TestFunction(V) F = u * v * ufl.dx + dt * kappa * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - u_n * v * ufl.dx - dt * f * v * ufl.dx + dt * g * v * ds(1) + dt * g * v * ds(2) + dt * g * v * ds(3) + dt * g * v * ds(4) a = fem.form(ufl.lhs(F)) L = fem.form(ufl.rhs(F)) A = create_matrix(a) b = create_vector(L) uh = fem.Function(V) solver = PETSc.KSP().create(domain.comm) solver.setOperators(A) solver.setType(PETSc.KSP.Type.PREONLY) solver.getPC().setType(PETSc.PC.Type.LU) errors_L2 = [] for n in range(num_steps): t += dt kappa_expression = np.sqrt(t) kappa.value = kappa_expression u_exact.t = t f_expression.t = t f_expression.kappa = kappa_expression f.interpolate(f_expression) A.zeroEntries() assemble_matrix(A, a) A.assemble() with b.localForm() as loc_b: loc_b.set(0) assemble_vector(b, L) apply_lifting(b, [a], [[]]) b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) set_bc(b, []) solver.solve(b, uh.vector) uh.x.scatter_forward() u_n.x.array[:] = uh.x.array V_ex = fem.functionspace(domain, ("Lagrange", 2)) u_ex = fem.Function(V_ex) u_ex.interpolate(u_exact) error_L2 = domain.comm.allreduce(fem.assemble_scalar(fem.form((uh - u_ex)**2 * ufl.dx)), op = MPI.SUM) errors_L2.append(error_L2) L2_error = np.sqrt(np.trapz(errors_L2, dx = dt)) print(f"L2-error: {L2_error:.2e}") end = time.time() print(f"Execution time: {(end - start) * 10 ** 3:.6f} ms") ``` then the system tells me ```auto Traceback (most recent call last): File "/media/math/个人文件空间/WangRan/work-fenics/SIR/demo4/ex.py", line 69, in g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type) IndexError: index 37 is out of bounds for axis 0 with size 36 ``` I guess the error results from the definition of `g` since it is defined on the square, but how to define a function only on facet? Moreover, is there any simple way to define the function `g`? Thanks very much. ---
**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\ **Post date:** [July 22, 2024, 7:59am UTC](https://fenicsproject.discourse.group/t/how-to-define-a-function-defined-on-boundary/15306/2 "2024-07-22T07:59:45Z")
`g.x.array` is indexed over the degrees of freedom, not the facets. Once you have `facets`, find the corresponding dofs using [dolfinx.fem — DOLFINx 0.9.0.0 documentation](https://docs.fenicsproject.org/dolfinx/main/python/generated/dolfinx.fem.html#dolfinx.fem.locate_dofs_topological) ---
**Author:** ![OsbertWang](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/osbertwang/32/2138_2.png) [@OsbertWang](https://fenicsproject.discourse.group/u/OsbertWang)\ **Post date:** [July 23, 2024, 1:46am UTC](https://fenicsproject.discourse.group/t/how-to-define-a-function-defined-on-boundary/15306/3 "2024-07-23T01:46:57Z")
Thanks. Now I define `g` as ```auto g = fem.Function(V) domain.topology.create_connectivity(fdim, domain.topology.dim) bottom_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(1)) top_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(2)) left_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(3)) right_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(4)) g.x.array[bottom_dofs] = np.full_like(bottom_dofs, 0.0, dtype = default_scalar_type) g.x.array[top_dofs] = np.full_like(top_dofs, -6.0 * kappa_expression, dtype = default_scalar_type) g.x.array[left_dofs] = np.full_like(left_dofs, 0.0, dtype = default_scalar_type) g.x.array[right_dofs] = np.full_like(right_dofs, -3.0 * kappa_expression, dtype = default_scalar_type) ``` Also, I update `g.x.array` in the iteration, and I obtain the error ```auto L2-error: 6.86e-01 ``` I don’t know why the error is so large. ---
**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\ **Post date:** [July 23, 2024, 7:26am UTC](https://fenicsproject.discourse.group/t/how-to-define-a-function-defined-on-boundary/15306/4 "2024-07-23T07:26:48Z")
You should carry out a convergence analysis, i.e. progressively double `nx` and `ny` and see if the FE solution converges at the expected rate. ---
**Author:** ![OsbertWang](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/osbertwang/32/2138_2.png) [@OsbertWang](https://fenicsproject.discourse.group/u/OsbertWang)\ **Post date:** [July 23, 2024, 8:43am UTC](https://fenicsproject.discourse.group/t/how-to-define-a-function-defined-on-boundary/15306/5 "2024-07-23T08:43:53Z")
Thank you for your advice. I change the `nx` and `ny` from 2^1 to 2^7, then obtain the errors: ```auto L2-error: 1.81e+00 L2-error: 8.83e-01 L2-error: 3.83e-01 L2-error: 1.28e-01 L2-error: 5.67e-02 L2-error: 1.01e-01 L2-error: 1.32e-01 ``` I think these results show that my code is wrong… However, I haven’t found the error. ---
**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\ **Post date:** [July 23, 2024, 9:55am UTC](https://fenicsproject.discourse.group/t/how-to-define-a-function-defined-on-boundary/15306/6 "2024-07-23T09:55:15Z")
The suggestion would be to start from/go back to a stationary problem, and not progress to an unsteady one until your steady code has the right convergence rates.