Hello Everyone. Hope you all are good.
The weak formulation of my problem is to seek (q_{h},u_{h},\lambda_{h}) \in V_{h}\times W_{h} \times M_{h} such that
\int_{\Omega} (1+u_{h}^{2})q_{h}v_{h} dx - \sum_{T}\int_{T} u_{h}div(v_{h})dx + \sum_{e\in \mathcal{E}_{i,h}} \int_{e} \lambda_{h}[[v_{h}]]ds = \int_{\partial{\Omega}} g[[v_{h}]]ds \forall v_{h}\in V_{h}
\sum_{T}\int_{T}w_{h}div{q_{h}}dx = \int_{\Omega} fw_{h}dx \forall w_{h}\in W_{h}
\sum_{e \in \mathcal{E}_{i,h}} \int_{e} m_{h}[[q_{h}]] ds = 0 \forall m_{h}\in M_{h}
where
V_{h}=\left\{v \in L^{2}(\Omega)^{2} \mbox{such that} v|_{T} \in P^{k}(T)\times P^{k}(T)+xP^{k}(T) \forall T\in \mathcal{T}_{h}\right\}.
Here V_{h} is not a subspace of H(div,\Omega).
W_{h} = \left\{ w \in L^{2}(\Omega) \mbox{such that} w|_{T} \in P^{k}(T) \forall T \in \mathcal{T}_{h}\right\},
M_{h}=\left\{w\in L^{2}(\mathcal{E}_{i,h}) \mbox{such that} w|_{e}\in P^{k}(e) \forall e\in \mathcal{E}_{i,h}\right\} ,
\mathcal{E}_{i,h} be the collection of edges of \mathcal{T}_{h} that are in the interior of the domain \Omega.
from dolfin import *
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(8,8)
V = FunctionSpace ( mesh, "BDM", 1)
W = FunctionSpace ( mesh, "DG", 1)
M = FunctionSpace ( mesh, "DG", 1 , restriction = 'facet')
Ve = FiniteElement ( "BDM", triangle, 1 )
We = FiniteElement ( "DG", triangle, 1 )
Me = FiniteElement ( "DG", triangle, 1)
VWMe = MixedElement ( Ve, We, Me )
VWMs = FunctionSpace ( mesh, VWMe )
(q, u, e) = TrialFunctions(VWMs)
(qt, ut, et) = TestFunctions(VWMs)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
n = FacetNormal(mesh)
F = inner((1 + u**2)*q, qt)*dx - u*div(qt)*dx + e*jump(qt,n)*dS + ut*div(q)*dx - f*ut*dx + et*jump(q,n)*dS
u_D = Constant(0.0)
bc = DirichletBC ( VWMs.sub(1), u_D, DomainBoundary())
file_u1 = File('hybridmixed/q.pvd')
file_u2 = File('hybridmixed/u.pvd')
file_u3 = File('hybridmixed/e.pvd')
que = Function(VWMs)
snes_solver_parameters = {"nonlinear_solver": "snes",
"snes_solver": {"linear_solver": "lu",
"maximum_iterations": 20,
"report": True,
"error_on_nonconvergence": False}}
problem = NonlinearVariationalProblem(F , que, bc)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
info(solver.parameters, True)
(q, u, e) = que.split()
file_u1 << q
file_u2 << u
file_u3 << e
This is the error what I got.
---------------------------------------------------------------------------
UFLException Traceback (most recent call last)
<ipython-input-92-d658d2e193cb> in <module>
7 "error_on_nonconvergence": False}}
8 #solve(F == 0, u , bc)
----> 9 problem = NonlinearVariationalProblem(F , que, bc)
10 solver = NonlinearVariationalSolver(problem)
11 solver.parameters.update(snes_solver_parameters)
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py in __init__(self, F, u, bcs, J, form_compiler_parameters)
155
156 # Wrap forms
--> 157 F = Form(F, form_compiler_parameters=form_compiler_parameters)
158 if J is not None:
159 J = Form(J, form_compiler_parameters=form_compiler_parameters)
/usr/lib/petsc/lib/python3/dist-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
/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in mpi_jit(*args, **kwargs)
48 # Just call JIT compiler when running in serial
49 if MPI.size(mpi_comm) == 1:
---> 50 return local_jit(*args, **kwargs)
51
52 # Default status (0 == ok, 1 == fail)
/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in ffc_jit(ufl_form, form_compiler_parameters)
98 p.update(dict(parameters["form_compiler"]))
99 p.update(form_compiler_parameters or {})
--> 100 return ffc.jit(ufl_form, parameters=p)
101
102
/usr/lib/python3/dist-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
/usr/lib/python3/dist-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
/usr/lib/python3/dist-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)
/usr/lib/python3/dist-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,
/usr/lib/python3/dist-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
/usr/lib/python3/dist-packages/ffc/compiler.py in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
183 # Stage 1: analysis
184 cpu_time = time()
--> 185 analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
186 _print_timing(1, time() - cpu_time)
187
/usr/lib/python3/dist-packages/ffc/analysis.py in analyze_ufl_objects(ufl_objects, kind, parameters)
88 # Analyze forms
89 form_datas = tuple(_analyze_form(form, parameters)
---> 90 for form in forms)
91
92 # Extract unique elements accross all forms
/usr/lib/python3/dist-packages/ffc/analysis.py in <genexpr>(.0)
88 # Analyze forms
89 form_datas = tuple(_analyze_form(form, parameters)
---> 90 for form in forms)
91
92 # Extract unique elements accross all forms
/usr/lib/python3/dist-packages/ffc/analysis.py in _analyze_form(form, parameters)
172 do_apply_geometry_lowering=True,
173 preserve_geometry_types=(Jacobian,),
--> 174 do_apply_restrictions=True)
175 elif r == "tsfc":
176 try:
/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
319 # Propagate restrictions to terminals
320 if do_apply_restrictions:
--> 321 form = apply_restrictions(form)
322
323 # If in real mode, remove any complex nodes introduced during form processing.
/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py in apply_restrictions(expression)
163 rules = RestrictionPropagator()
164 return map_integrand_dags(rules, expression,
--> 165 only_integral_type=integral_types)
166
167
/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in map_integrand_dags(function, form, only_integral_type, compress)
45 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
46 return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
---> 47 form, only_integral_type)
/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in map_integrands(function, form, only_integral_type)
26 if isinstance(form, Form):
27 mapped_integrals = [map_integrands(function, itg, only_integral_type)
---> 28 for itg in form.integrals()]
29 nonzero_integrals = [itg for itg in mapped_integrals
30 if not isinstance(itg.integrand(), Zero)]
/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in <listcomp>(.0)
26 if isinstance(form, Form):
27 mapped_integrals = [map_integrands(function, itg, only_integral_type)
---> 28 for itg in form.integrals()]
29 nonzero_integrals = [itg for itg in mapped_integrals
30 if not isinstance(itg.integrand(), Zero)]
/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in map_integrands(function, form, only_integral_type)
33 itg = form
34 if (only_integral_type is None) or (itg.integral_type() in only_integral_type):
---> 35 return itg.reconstruct(function(itg.integrand()))
36 else:
37 return itg
/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in <lambda>(expr)
44
45 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
---> 46 return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
47 form, only_integral_type)
/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dag(function, expression, compress)
24 Return the result of the final function call.
25 """
---> 26 result, = map_expr_dags(function, [expression], compress=compress)
27 return result
28
/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dags(function, expressions, compress)
71 # Cache miss: Get transformed operands, then apply transformation
72 if cutoff_types[v._ufl_typecode_]:
---> 73 r = handlers[v._ufl_typecode_](v)
74 else:
75 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py in reference_value(self, o)
86 f, = o.ufl_operands
87 assert f._ufl_is_terminal_
---> 88 g = self(f)
89 if isinstance(g, Restricted):
90 side = g.side()
/usr/lib/python3/dist-packages/ufl/corealg/multifunction.py in __call__(self, o, *args)
87 def __call__(self, o, *args):
88 "Delegate to handler function based on typecode of first argument."
---> 89 return self._handlers[o._ufl_typecode_](o, *args)
90
91 def undefined(self, o, *args):
/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py in _require_restriction(self, o)
45 "Restrict a discontinuous quantity to current side, require a side to be set."
46 if self.current_restriction is None:
---> 47 error("Discontinuous type %s must be restricted." % o._ufl_class_.__name__)
48 return o(self.current_restriction)
49
/usr/lib/python3/dist-packages/ufl/log.py in error(self, *message)
156 "Write error message and raise an exception."
157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))
159
160 def begin(self, *message):
UFLException: Discontinuous type Argument must be restricted.
Kindly help me to resolve it.