Discontinuous type Coefficient must be restricted type error

Ashini · May 30, 2024, 4:02pm

Hello Everyone

I am Aadarshini, working in Multiphenics and need some help. I define the Lagrange Multiplier (LM) in DG1 space defined by “P0” and lam_ex is the exact solution for LM. then e_la is the error in the L^2 norm but it shows error.

P2v = VectorFunctionSpace(mesh, "CG", 2)
P2 = VectorFunctionSpace(mesh, "CG", 1)
P1  = FunctionSpace(mesh, "CG", 1)
P0  = FunctionSpace(mesh, "DG", 1)
Hh = BlockFunctionSpace([P2v, P2v, P2v, P1, P1, P0], restrict = [OmS,  OmH, OmH, OmS,  OmH, Sig])
trial = BlockTrialFunction(Hh)
uf, ur,  ys, pf, ph, la= block_split(trial)
test  = BlockTestFunction(Hh)
vf, vr,  ws, qf, qh, ze= block_split(test)

lam_ex = dot(sigma_exS('+')*n('+'),n('+'))
 e_la = e_la + dt*abs(assemble(la_h - lam_ex)**2*dS(interf))

error:

ufl_legacy.log.UFLException: Discontinuous type Coefficient must be restricted.

I think this happens because I am using DG function and lam_ex is also DG type and I computed the error norm for other variables but they didn’t show any error. Please help me in this regard.

Thanks

francesco-ballarin · May 30, 2024, 4:08pm

I guess your are missing a restriction in la_h, but it’s impossible to confirm or deny the guess but you copied a very narrow snippet of your code.

Please read Read before posting: How do I get my question answered?

Ashini · May 30, 2024, 6:15pm

Below I attach minimal example. The problem is that if I consider Dg1 element for LM then lam_ex (exact LM) create problem. I don’t know how to resolve it.

from __future__ import print_function
from dolfin import *
from multiphenics import *

def epsilon(u):
    return 0.5*(grad(u)+grad(u).T) 


nu    = Constant(1)
alpha = Constant(1.)
kappa = Constant(1.)


# ******* construct mesh and define normal, tangent ****** #
darcy = 10; stokes = 13; outlet = 14; inlet = 15;
interf = 16; wallS = 17; wallD = 18
nps = 2
mesh = RectangleMesh(Point(0,-1.0),Point(1.0,1.0),nps,2*nps)
subdomains = MeshFunction("size_t", mesh, 2)
subdomains.set_all(0)
boundaries = MeshFunction("size_t", mesh, 1)
boundaries.set_all(0)

class Top(SubDomain):
      def inside(self, x, on_boundary):
          return (near(x[1], 1.0) and on_boundary)

class SRight(SubDomain):
      def inside(self, x, on_boundary):
          return (near(x[0], 1.0) and between(x[1], (0.0, 1.0)) and on_boundary)

class SLeft(SubDomain):
      def inside(self, x, on_boundary):
          return (near(x[0], 0.0) and between(x[1], (0.0, 1.0)) and on_boundary)

class DRight(SubDomain):
      def inside(self, x, on_boundary):
          return (near(x[0], 1.0) and between(x[1], (-1.0, 0.0)) and on_boundary)

class DLeft(SubDomain):
      def inside(self, x, on_boundary):
          return (near(x[0], 0.0) and between(x[1], (-1.0, 0.0)) and on_boundary)

class Bot(SubDomain):
      def inside(self, x, on_boundary):
          return (near(x[1], -1.0) and on_boundary)

class MStokes(SubDomain):
     def inside(self, x, on_boundary):
         return x[1]>=0.

class MDarcy(SubDomain):
      def inside(self, x, on_boundary):
          return x[1]<=0.  

class Interface(SubDomain):
      def inside(self, x, on_boundary):
          return near(x[1], 0.0)

MDarcy().mark(subdomains, darcy)
MStokes().mark(subdomains, stokes)
Interface().mark(boundaries, interf)

Top().mark(boundaries,inlet)
SRight().mark(boundaries,wallS)
SLeft().mark(boundaries,wallS)
DRight().mark(boundaries,wallD)
DLeft().mark(boundaries,wallD)
Bot().mark(boundaries,outlet)


n = FacetNormal(mesh); t = as_vector((-n[1], n[0]))

# ******* set subdomains, boundaries, and interface ****** #

OmS = MeshRestriction(mesh, MStokes())
OmD = MeshRestriction(mesh, MDarcy())
Sig = MeshRestriction(mesh, Interface())

dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dS = Measure("dS", domain=mesh, subdomain_data=boundaries)

p_exD = Expression("x[0]*x[1]*(1-x[0])*(1-x[1])", degree= 3, domain =mesh)

#extra terms due to manufactured solutions
lam_ex = -p_exD('-')


# ***** Global FE spaces and their restrictions ****** #

P2v = VectorFunctionSpace(mesh, "CG", 2)
P1  = FunctionSpace(mesh, "CG", 1)
RT0 = FunctionSpace(mesh, "RT", 1)
P0  = FunctionSpace(mesh, "DG", 0)

#                         uS,   uD, pS, pD, la, beta0   
Hh = BlockFunctionSpace([P2v, RT0, P1, P0, P0], restrict = [OmS, OmD, OmS, OmD, Sig])

trial = BlockTrialFunction(Hh)
uS, uD, pS, pD, la = block_split(trial)
test  = BlockTestFunction(Hh)
vS, vD, qS, qD, ze = block_split(test)

print("DoFs = ", Hh.dim())

# ******** Other parameters and BCs ************* #

bcUin  = DirichletBC(Hh.sub(0), Constant((0,0)), boundaries, inlet)
bcUS   = DirichletBC(Hh.sub(0),Constant((0,0)), boundaries, wallS)
bcUD   = DirichletBC(Hh.sub(3), Constant(0), boundaries, wallD)
bcs = BlockDirichletBC([bcUin, bcUS, bcUD])

# ********  Define weak forms ********** #

af   = 2.0 * nu * inner(epsilon(uS), epsilon(vS)) * dx(stokes) + alpha*nu/sqrt(kappa)*dot(uS('+'),t('+'))*dot(vS('+'),t('+'))*dS(interf)

ap   = nu/kappa * dot(uD, vD) * dx(darcy)

bf1t = - pS * div(vS) * dx(stokes)
bf1  = - qS * div(uS) * dx(stokes)
bp1t = - pD * div(vD) * dx(darcy)
bp1  = - qD * div(uD) * dx(darcy)

bI1t = avg(la) * dot(vS('+'),n('+')) * dS(interf)
bI1  = avg(ze) * dot(uS('+'),n('+')) * dS(interf)
bI2t = avg(la) * dot(vD('-'),n('-')) * dS(interf)
bI2  = avg(ze) * dot(uD('-'),n('-')) * dS(interf)

FvS = dot(Constant((0,0)),vS)*dx(stokes)

FvD = 0
GqS = 0
GqD = 0

# ****** Assembly and solution of linear system ******** #

rhs = [FvS, FvD, GqS, GqD, 0]

# ordering of the unknowns 
#        uS   uD   pS     pD  la   beta0
lhs = [[ af,   0, bf1t,    0, bI1t ],  #vS
       [  0,  ap,    0, bp1t, bI2t ],  #vD
       [bf1,   0,    0,    0,    0 ],  #qS
       [  0, bp1,    0,    0,    0 ],  #qD
       [bI1, bI2,    0,    0,    0 ]]

AA = block_assemble(lhs)
FF = block_assemble(rhs)
bcs.apply(AA)
bcs.apply(FF)

sol = BlockFunction(Hh)
block_solve(AA, sol.block_vector(), FF, linear_solver="mumps")
uS_h, uD_h, pS_h, pD_h, la_h = block_split(sol)


epD = assemble((pD_h - p_exD)**2*dx(darcy))
elam = assemble((la_h - lam_ex)**2*dS(interf))

error:

   elam = assemble((la_h - lam_ex)**2*dS(interf))
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
    return Form(form,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 50, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 100, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/compute_form_data.py", line 321, in compute_form_data
    form = apply_restrictions(form)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/apply_restrictions.py", line 170, in apply_restrictions
    return map_integrand_dags(rules, expression,
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 46, in map_integrand_dags
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 27, in map_integrands
    mapped_integrals = [map_integrands(function, itg, only_integral_type)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 27, in <listcomp>
    mapped_integrals = [map_integrands(function, itg, only_integral_type)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 35, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 46, in <lambda>
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/usr/lib/python3/dist-packages/ufl_legacy/corealg/map_dag.py", line 36, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/usr/lib/python3/dist-packages/ufl_legacy/corealg/map_dag.py", line 97, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/apply_restrictions.py", line 93, in reference_value
    g = self(f)
  File "/usr/lib/python3/dist-packages/ufl_legacy/corealg/multifunction.py", line 89, in __call__
    return self._handlers[o._ufl_typecode_](o, *args)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/apply_restrictions.py", line 140, in coefficient
    return self._require_restriction(o)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/apply_restrictions.py", line 52, in _require_restriction
    error("Discontinuous type %s must be restricted." % o._ufl_class_.__name__)
  File "/usr/lib/python3/dist-packages/ufl_legacy/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl_legacy.log.UFLException: Discontinuous type Coefficient must be restricted.

Jesus_Vellojin · May 30, 2024, 6:37pm

Hi.

As stated by @francesco-ballarin, you are missing a restriction in the integration over the inner facets. The error is in the line

elam = assemble((la_h - lam_ex)**2*dS(interf))

because la_h is not restricted to + or -.

With respect to consider DG1 or DG0, you could use DGTk (Discontinuous Galerkin Traces) to define the element over the facets of the restriction (similar to the Stokes-Darcy tutorial of the legacy Multiphenics).

The following should work:

from __future__ import print_function
from dolfin import *
from multiphenics import *


def epsilon(u):
    return 0.5*(grad(u)+grad(u).T)


nu = Constant(1)
alpha = Constant(1.)
kappa = Constant(1.)


# ******* construct mesh and define normal, tangent ****** #
darcy = 10
stokes = 13
outlet = 14
inlet = 15
interf = 16
wallS = 17
wallD = 18
nps = 2
mesh = RectangleMesh(Point(0, -1.0), Point(1.0, 1.0), nps, 2*nps)
subdomains = MeshFunction("size_t", mesh, 2)
subdomains.set_all(0)
boundaries = MeshFunction("size_t", mesh, 1)
boundaries.set_all(0)


class Top(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[1], 1.0) and on_boundary)


class SRight(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 1.0) and between(x[1], (0.0, 1.0)) and on_boundary)


class SLeft(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 0.0) and between(x[1], (0.0, 1.0)) and on_boundary)


class DRight(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 1.0) and between(x[1], (-1.0, 0.0)) and on_boundary)


class DLeft(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 0.0) and between(x[1], (-1.0, 0.0)) and on_boundary)


class Bot(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[1], -1.0) and on_boundary)


class MStokes(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 0.


class MDarcy(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= 0.


class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)


MDarcy().mark(subdomains, darcy)
MStokes().mark(subdomains, stokes)
Interface().mark(boundaries, interf)

Top().mark(boundaries, inlet)
SRight().mark(boundaries, wallS)
SLeft().mark(boundaries, wallS)
DRight().mark(boundaries, wallD)
DLeft().mark(boundaries, wallD)
Bot().mark(boundaries, outlet)


n = FacetNormal(mesh)
t = as_vector((-n[1], n[0]))

# ******* set subdomains, boundaries, and interface ****** #

OmS = MeshRestriction(mesh, MStokes())
OmD = MeshRestriction(mesh, MDarcy())
Sig = MeshRestriction(mesh, Interface())

dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dS = Measure("dS", domain=mesh, subdomain_data=boundaries)

p_exD = Expression("x[0]*x[1]*(1-x[0])*(1-x[1])", degree=3, domain=mesh)

# extra terms due to manufactured solutions
lam_ex = -p_exD('-')


# ***** Global FE spaces and their restrictions ****** #

P2v = VectorFunctionSpace(mesh, "CG", 2)
P1 = FunctionSpace(mesh, "CG", 1)
RT0 = FunctionSpace(mesh, "RT", 1)
P0 = FunctionSpace(mesh, "DG", 0)
P0T = FunctionSpace(mesh, "DGT", 0) # it can be 1 and also converges
#                         uS,   uD, pS, pD, la, beta0
Hh = BlockFunctionSpace([P2v, RT0, P1, P0, P0T], restrict=[
                        OmS, OmD, OmS, OmD, Sig])

trial = BlockTrialFunction(Hh)
uS, uD, pS, pD, la = block_split(trial)
test = BlockTestFunction(Hh)
vS, vD, qS, qD, ze = block_split(test)

print("DoFs = ", Hh.dim())

# ******** Other parameters and BCs ************* #

bcUin = DirichletBC(Hh.sub(0), Constant((0, 0)), boundaries, inlet)
bcUS = DirichletBC(Hh.sub(0), Constant((0, 0)), boundaries, wallS)
bcUD = DirichletBC(Hh.sub(3), Constant(0), boundaries, wallD)
bcs = BlockDirichletBC([bcUin, bcUS, bcUD])

# ********  Define weak forms ********** #

af = 2.0 * nu * inner(epsilon(uS), epsilon(vS)) * dx(stokes) + alpha * \
    nu/sqrt(kappa)*dot(uS('+'), t('+'))*dot(vS('+'), t('+'))*dS(interf)

ap = nu/kappa * dot(uD, vD) * dx(darcy)

bf1t = - pS * div(vS) * dx(stokes)
bf1 = - qS * div(uS) * dx(stokes)
bp1t = - pD * div(vD) * dx(darcy)
bp1 = - qD * div(uD) * dx(darcy)

bI1t = avg(la) * dot(vS('+'), n('+')) * dS(interf)
bI1 = avg(ze) * dot(uS('+'), n('+')) * dS(interf)
bI2t = avg(la) * dot(vD('-'), n('-')) * dS(interf)
bI2 = avg(ze) * dot(uD('-'), n('-')) * dS(interf)

FvS = dot(Constant((0, 0)), vS)*dx(stokes)

FvD = 0
GqS = 0
GqD = 0

# ****** Assembly and solution of linear system ******** #

rhs = [FvS, FvD, GqS, GqD, 0]

# ordering of the unknowns
#        uS   uD   pS     pD  la   beta0
lhs = [[af,   0, bf1t,    0, bI1t],  # vS
       [0,  ap,    0, bp1t, bI2t],  # vD
       [bf1,   0,    0,    0,    0],  # qS
       [0, bp1,    0,    0,    0],  # qD
       [bI1, bI2,    0,    0,    0]]

AA = block_assemble(lhs)
FF = block_assemble(rhs)
bcs.apply(AA)
bcs.apply(FF)

sol = BlockFunction(Hh)
block_solve(AA, sol.block_vector(), FF, linear_solver="mumps")
uS_h, uD_h, pS_h, pD_h, la_h = block_split(sol)


epD = assemble((pD_h - p_exD)**2*dx(darcy))
elam = assemble((avg(la_h) - lam_ex)**2*dS(interf)) # la_h is restricted trough avg()

Hope this helps.

Topic		Replies	Views
Discontinuous type Argument must be restricted Errors	4	2081	July 10, 2021
[MixedFunctionSpace][FEniCS] Error in building the mapping Errors mesh , legacy-dolfin	4	62	December 11, 2024
1d dG advection/diffusion probem dolfinx boundary-conditions	3	91	November 12, 2024
Deal with Error: too many values to unpack while solving mixed element space Linear Algebra	7	1993	June 23, 2020
Error norm of the solution to the variational problem	3	1197	June 24, 2019

Discontinuous type Coefficient must be restricted type error

Related topics