Manually compute normal vector on inner boundary/subdomain

Hi, I am trying to get the normal vector on the inner subdomain. When I use
n = FacetNormal(mesh)
I only get the normal vector outer boundary. I use the unit square mesh and have two subdomains; one circle and the other is the whole domain - the circle. Is there a way to do this?

The facetnormal Also works for inner cells. However, you have to restrict it to either + or minus side, i.e. n("+”) or n("-").
For instance consider assemble(inner(avg(grad(u)),jump(v,n))*ds

Thank you for your reply. I am not able to use n(‘-’) and n(‘+’)… I tried the following working example:

V = VectorFunctionSpace(mesh, “CG”, 1)
n = FacetNormal(mesh)

f = Expression((‘1.0’, ‘1.0’), degree=2)
f = interpolate(f, V)

g1 = dot(f, n(‘+’))*ds(2)
flux_1 = assemble(g1)

And what is your error message?

I tried with both py2.7 and 3.7.

Error for 3.7 is:
AssertionError Traceback (most recent call last)
in
6
7 g1 = dot(f, n(’-’))*ds(2)
----> 8 flux_1 = assemble(g1)
9
10 ### Second method ###

~/anaconda3/envs/ft3/lib/python3.7/site-packages/dolfin/fem/assembling.py in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
196 dolfin_form = form
197 else:
–> 198 dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
199
200 # Create tensor

~/anaconda3/envs/ft3/lib/python3.7/site-packages/dolfin/fem/assembling.py in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
56 return Form(form,
57 form_compiler_parameters=form_compiler_parameters,
—> 58 function_spaces=function_spaces)
59 else:
60 raise TypeError(“Invalid form type %s” % (type(form),))

~/anaconda3/envs/ft3/lib/python3.7/site-packages/dolfin/fem/form.py in init(self, form, **kwargs)
42
43 ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
—> 44 mpi_comm=mesh.mpi_comm())
45 ufc_form = cpp.fem.make_ufc_form(ufc_form[0])
46

~/anaconda3/envs/ft3/lib/python3.7/site-packages/dolfin/jit/jit.py in mpi_jit(*args, **kwargs)
45 # Just call JIT compiler when running in serial
46 if MPI.size(mpi_comm) == 1:
—> 47 return local_jit(*args, **kwargs)
48
49 # Default status (0 == ok, 1 == fail)

~/anaconda3/envs/ft3/lib/python3.7/site-packages/dolfin/jit/jit.py in ffc_jit(ufl_form, form_compiler_parameters)
95 p.update(dict(parameters[“form_compiler”]))
96 p.update(form_compiler_parameters or {})
—> 97 return ffc.jit(ufl_form, parameters=p)
98
99

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/jitcompiler.py in jit(ufl_object, parameters, indirect)
215
216 # Inspect cache and generate+build if necessary
–> 217 module = jit_build(ufl_object, module_name, parameters)
218
219 # Raise exception on failure to build or import module

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/jitcompiler.py in jit_build(ufl_object, module_name, parameters)
131 name=module_name,
132 params=params,
–> 133 generate=jit_generate)
134 return module
135

~/anaconda3/envs/ft3/lib/python3.7/site-packages/dijitso/jit.py in jit(jitable, name, params, generate, send, receive, wait)
163 elif generate:
164 # 1) Generate source code
–> 165 header, source, dependencies = generate(jitable, name, signature, params[“generator”])
166 # Ensure we got unicode from generate
167 header = as_unicode(header)

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/jitcompiler.py in jit_generate(ufl_object, module_name, signature, parameters)
64
65 code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
—> 66 prefix=module_name, parameters=parameters, jit=True)
67
68 # Jit compile dependent objects separately,

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/compiler.py in compile_form(forms, object_names, prefix, parameters, jit)
141 “”“This function generates UFC code for a given UFL form or list of UFL forms.”""
142 return compile_ufl_objects(forms, “form”, object_names,
–> 143 prefix, parameters, jit)
144
145

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/compiler.py in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
188 # Stage 2: intermediate representation
189 cpu_time = time()
–> 190 ir = compute_ir(analysis, prefix, parameters, jit)
191 _print_timing(2, time() - cpu_time)
192

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/representation.py in compute_ir(analysis, prefix, parameters, jit)
181 info(“Computing representation of integrals”)
182 irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
–> 183 for (form_id, fd) in enumerate(form_datas)]
184 ir_integrals = list(chain(*irs))
185

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/representation.py in (.0)
181 info(“Computing representation of integrals”)
182 irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
–> 183 for (form_id, fd) in enumerate(form_datas)]
184 ir_integrals = list(chain(*irs))
185

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/representation.py in _compute_integral_ir(form_data, form_id, prefix, element_numbers, classnames, parameters, jit)
458 element_numbers,
459 classnames,
–> 460 parameters)
461
462 # Remember jit status

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/uflacs/uflacsrepresentation.py in compute_integral_ir(itg_data, form_data, form_id, element_numbers, classnames, parameters)
137 coefficient_numbering,
138 quadrature_rules,
–> 139 parameters)
140 ir.update(uflacs_ir)
141

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/uflacs/build_uflacs_ir.py in build_uflacs_ir(cell, integral_type, entitytype, integrands, tensor_shape, coefficient_numbering, quadrature_rules, parameters)
333 for num_points, expressions in cases:
334 # Rebalance order of nested terminal modifiers
–> 335 expressions = [balance_modifiers(expr) for expr in expressions]
336
337 # Build initial scalar list-based graph representation

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/uflacs/build_uflacs_ir.py in (.0)
333 for num_points, expressions in cases:
334 # Rebalance order of nested terminal modifiers
–> 335 expressions = [balance_modifiers(expr) for expr in expressions]
336
337 # Build initial scalar list-based graph representation

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/uflacs/analysis/balancing.py in balance_modifiers(expr)
90 def balance_modifiers(expr):
91 mf = BalanceModifiers()
—> 92 return map_expr_dag(mf, expr)

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ufl/corealg/map_dag.py in map_expr_dag(function, expression, compress)
35 Return the result of the final function call.
36 “”"
—> 37 result, = map_expr_dags(function, [expression], compress=compress)
38 return result
39

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ufl/corealg/map_dag.py in map_expr_dags(function, expressions, compress)
84 r = handlersv.ufl_typecode
85 else:
—> 86 r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
87
88 # Optionally check if r is in rcache, a memory optimization

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/uflacs/analysis/balancing.py in _modifier(self, expr, *ops)
77
78 def _modifier(self, expr, *ops):
—> 79 return balance_modified_terminal(expr)
80
81 reference_value = _modifier

~/anaconda3/envs/ft3/lib/python3.7/site-packages/ffc/uflacs/analysis/balancing.py in balance_modified_terminal(expr)
52 layers = [expr]
53 while not expr.ufl_is_terminal:
—> 54 assert expr.ufl_is_terminal_modifier
55 expr = expr.ufl_operands[0]
56 layers.append(expr)

AssertionError:

and with 2.7:

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 perform just-in-time compilation of form.
  • Reason: ffc.jit failed with message:
    Traceback (most recent call last):
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/dolfin/compilemodules/jit.py”, line 142, in jit
    result = ffc.jit(ufl_object, parameters=p)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/jitcompiler.py”, line 218, in jit
    module = jit_build(ufl_object, module_name, parameters)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/jitcompiler.py”, line 134, in jit_build
    generate=jit_generate)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/dijitso/jit.py”, line 167, in jit
    header, source, dependencies = generate(jitable, name, signature, params[“generator”])
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/jitcompiler.py”, line 67, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/compiler.py”, line 143, in compile_form
    prefix, parameters, jit)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/compiler.py”, line 190, in compile_ufl_objects
    ir = compute_ir(analysis, prefix, parameters, jit)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/representation.py”, line 186, in compute_ir
    for (form_id, fd) in enumerate(form_datas)]
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/representation.py”, line 463, in _compute_integral_ir
    parameters)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/uflacs/uflacsrepresentation.py”, line 139, in compute_integral_ir
    parameters)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/uflacs/build_uflacs_ir.py”, line 335, in build_uflacs_ir
    expressions = [balance_modifiers(expr) for expr in expressions]
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/uflacs/analysis/balancing.py”, line 92, in balance_modifiers
    return map_expr_dag(mf, expr)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ufl/corealg/map_dag.py”, line 37, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ufl/corealg/map_dag.py”, line 86, in map_expr_dags
    r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/uflacs/analysis/balancing.py”, line 79, in _modifier
    return balance_modified_terminal(expr)
    File “/home/usr/anaconda3/envs/ft2/lib/python2.7/site-packages/ffc/uflacs/analysis/balancing.py”, line 54, in balance_modified_terminal
    assert expr.ufl_is_terminal_modifier
    AssertionError
    .
  • Where: This error was encountered inside jit.py.
  • Process: 0
  • DOLFIN version: 2017.2.0
  • Git changeset: 774aa9b05f5a21fcf3d1bd632e722933a05fdb45

You are using a really old version of FEniCS, I would suggest updating to latest release.

Is the one on py3.7 also an old release? I downloaded it today

I am using the 2019.1.0 version on py3.7, is there a newer one?

Is it possible to define three normal vectors? If I have three domains?

Could you supply the full code (that you ran in python3.7). It’s impossible to tell what you are doing wrong without the full code.

I was a bit too quick further up. You should use the dS measure. Here is an example:

from dolfin import *
mesh = UnitSquareMesh(3,3)
V = VectorFunctionSpace(mesh, "DG", 1)

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
domains.set_all(0)
class InteriorMarker(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1/3)

marker = InteriorMarker()
marker.mark(domains, 1)
u = Function(V)
u.vector()[:] = 1
n = FacetNormal(mesh)
dInterior = Measure("dS", subdomain_data=domains, subdomain_id=1)
print(assemble(inner(avg(u), n("+"))*dInterior))
1 Like

Thank you! I think that works, but I get an error when I try to use this for my whole system. Do you know what to do when getting this error message: “Discontinuous type Jacobian must be restricted.”?

If your function is discontinuous, i.e. if you use grad(u) on a CG function, you need to restrict it to either side of the interface. You have to really consider the interface conditions properly, as you will end up with jumps and averages that has to be handled correctly in an variational setting. See for instance: https://onlinelibrary.wiley.com/doi/epdf/10.1002/nme.4823

1 Like

Thanks a lot! I will have a closer look at the article. Is it possible to plot the normal vectors?

You would have to project it to an appropriate function space. However, https://bitbucket.org/Epoxid/femorph/src/c7317791c8f00d70fe16d593344cb164a53cad9b/femorph/SAD_geometry.py?at=dokken%2Frestructuring
Has several different functions for finding appropriate representations of the normal function.

1 Like