Set boundary condition for a tensor field

variational formulation
1

Hi everyone. I am working with a problem where I have to apply a no-flux boundary condition on a tensor-valued field. More specifically, I would like to set \boldsymbol{C} \cdot \boldsymbol{n} = \boldsymbol{0} on the boundary of my domain, where C is a second-rank tensor (3 x 3 components). To do so, I have tried the following:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from __future__ import print_function
from ufl import *  # imports everything from ufl
from dolfin import *  # imports everything from dolfin
import numpy as np
import os
import sys

set_log_level(LogLevel.ERROR)

Lbox = 20.0  # Length of the simulation box in x-direction
Hbox = 4.0  # Length of the simulation box in y-direction
Wbox = 100000.0  # Length of the simulation box in z-direction

N_el_x, N_el_y, N_el_z = 100, 20, 1
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Lbox, Hbox, Wbox), N_el_x, N_el_y, N_el_z)

# Define function space and functions
V_C   = TensorFunctionSpace(mesh, "CG", 2)

# Define function C_bound shuch that C_bound \cdot n = g
class BoundarySource(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = 0.
        values[0] = g*n[0]
        values[1] = g*n[1]
        values[2] = g*n[2]
    def value_shape(self):
        return (3,)

C_bound = BoundarySource(mesh, degree=0)

# Define the boundaries of the domain
def boundary(x, on_boundary):
    return (near(x[0], 0.) or near(x[0], Lbox) or \
            near(x[1], 0.) or near(x[1], Hbox)) and on_boundary

# Define the essential boundary condition on chi
bc = DirichletBC(V_C, C_bound, boundary)

This throws me the following error:

Traceback (most recent call last):
  File "MWE.py", line 46, in <module>
    bc = DirichletBC(V_C, C_bound, boundary)
  File "/home/lima/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/dirichletbc.py", line 131, in __init__
    super().__init__(*args)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  Illegal value rank (1), expecting (2).
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  149cf7bbfffe157a560d262289a63cfe2d1e7244
*** -------------------------------------------------------------------------

Could you help me solving this issue?

2

First of all, you need to use the correct shape for your expression:

should be:

    def value_shape(self):
        return (3, 3)

If you want to apply this condition on boundaries aligning with the coordinate system, i.e. n=(1,0,0), (0,1,0) or (0,0,1) you should set the boundary condition on the appropriate sub space instead of the full space.
If your boundaries (normals) are not aligning with the coordinate system, you need to enforce the condition weakly (through the variational form), or with a multi point constraint method.

2 Likes
3

Thanks for your answer! It worked fine, but I also changed eval_cell method to the following:

    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = 0.
        values[0] = g*n[0]
        values[1] = g*n[1]
        values[2] = g*n[2]
        values[3] = g*n[0]
        values[4] = g*n[1]
        values[5] = g*n[2]
        values[6] = g*n[0]
        values[7] = g*n[1]
        values[8] = g*n[2]

In this specific case, as my domain is rectangular, the face normals are aligned with the coordinate systems, but if it were not the case, I think that this implementation should also work (tell me if I’m wrong).

4

Actually, by doing this, I’m setting all components C_{ij} = 0 on the boundary, which may not be true. Would it be possible to set C_{i1}n_1 + C_{i2}n_2 + C_{i3}n_3 = 0 for i = 1, 2, 3 ?

5

If C Is your unknown, then dot(C,n)=0 has to be enforced by either a Multi point constraint or weakly through your variational form.

1 Like