Thanks for guidance and resources, after going through it, I wrote couple of lines -
x, y = SpatialCoordinate(mesh)
condition = lt(y, 0.5)
below = 1.*T
above = 1./T
varsigma = conditional(condition, below, above)
But running into the error -
True value of Condtional should only be one function: {((1, 'nzc3[k]', 3, 6),): Product([Symbol('FE0_C0[ip][k]', BASIS)])}
---------------------------------------------------------------------------
Exception Traceback (most recent call last)
/tmp/ipykernel_3265/318594624.py in <module>
29
30 while t < tMax :
---> 31 solve(Form == 0,unkn,bc,J=Gain ,solver_parameters = solver_parameters,form_compiler_parameters = form_compiler_parameters)
32 unkn0.assign(unkn)
33 t += Dt
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
231 # tolerance)
232 elif isinstance(args[0], ufl.classes.Equation):
--> 233 _solve_varproblem(*args, **kwargs)
234
235 # Default case, just call the wrapped C++ solve function
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
306
307 # Create problem
--> 308 problem = NonlinearVariationalProblem(eq.lhs, u, bcs, J,
309 form_compiler_parameters=form_compiler_parameters)
310
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py in __init__(self, F, u, bcs, J, form_compiler_parameters)
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)
160
161 # Initialize C++ base class
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py in __init__(self, form, **kwargs)
41 form_compiler_parameters = {"external_include_dirs": dolfin_pc["include_dirs"]}
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])
/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)
128
129 # Carry out jit compilation, calling jit_generate only if needed
--> 130 module, signature = dijitso.jit(jitable=ufl_object,
131 name=module_name,
132 params=params,
/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)
63 compile_object = compile_coordinate_mapping
64
---> 65 code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
66 prefix=module_name, parameters=parameters, jit=True)
67
/usr/lib/python3/dist-packages/ffc/compiler.py in compile_form(forms, object_names, prefix, parameters, jit)
140 prefix="Form", parameters=None, jit=False):
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
/usr/lib/python3/dist-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
/usr/lib/python3/dist-packages/ffc/representation.py in compute_ir(analysis, prefix, parameters, jit)
180 # Compute and flatten representation of integrals
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))
/usr/lib/python3/dist-packages/ffc/representation.py in <listcomp>(.0)
180 # Compute and flatten representation of integrals
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))
/usr/lib/python3/dist-packages/ffc/representation.py in _compute_integral_ir(form_data, form_id, prefix, element_numbers, classnames, parameters, jit)
453
454 # Compute representation
--> 455 ir = r.compute_integral_ir(itg_data,
456 form_data,
457 form_id, # FIXME: Can we remove this?
/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py in compute_integral_ir(itg_data, form_data, form_id, element_numbers, classnames, parameters)
102 # Transform integrals.
103 cell = itg_data.domain.ufl_cell()
--> 104 ir["trans_integrals"] = _transform_integrals_by_type(ir, transformer,
105 integrals_dict,
106 itg_data.integral_type,
/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py in _transform_integrals_by_type(ir, transformer, integrals_dict, integral_type, cell)
180 info("Transforming cell integral")
181 transformer.update_cell()
--> 182 terms = _transform_integrals(transformer, integrals_dict, integral_type)
183
184 elif integral_type == "exterior_facet":
/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py in _transform_integrals(transformer, integrals, integral_type)
227 for point, integral in sorted(integrals.items()):
228 transformer.update_points(point)
--> 229 terms = transformer.generate_terms(integral.integrand(), integral_type)
230 transformed_integrals.append((point, terms, transformer.function_data,
231 {}, transformer.coordinate,
/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py in generate_terms(self, integrand, integral_type)
805
806 # Get terms
--> 807 terms = self.visit(integrand)
808
809 # Get formatting
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
101 # No, this is a handler that handles its own children
102 # (arguments self and o, where self is already bound)
--> 103 r = h(o)
104
105 # Update stack and return
/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py in index_sum(self, o)
591 for i in range(o.dimension()):
592 self._index2value.push(index, i)
--> 593 ops.append(self.visit(summand))
594 self._index2value.pop()
595
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
101 # No, this is a handler that handles its own children
102 # (arguments self and o, where self is already bound)
--> 103 r = h(o)
104
105 # Update stack and return
/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py in indexed(self, o)
556
557 # Visit expression subtrees and generate code.
--> 558 code = self.visit(indexed_expr)
559
560 # Remove component again
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
101 # No, this is a handler that handles its own children
102 # (arguments self and o, where self is already bound)
--> 103 r = h(o)
104
105 # Update stack and return
/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py in component_tensor(self, o)
760
761 # Visit expression subtrees and generate code.
--> 762 code = self.visit(component_expr)
763
764 # Remove the index map from the StackDict
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in <listcomp>(.0)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ufl/algorithms/transformer.py in visit(self, o)
97 if visit_children_first:
98 # Yes, visit all children first and then call h.
---> 99 r = h(o, *[self.visit(op) for op in o.ufl_operands])
100 else:
101 # No, this is a handler that handles its own children
/usr/lib/python3/dist-packages/ffc/quadrature/optimisedquadraturetransformer.py in conditional(self, o, *operands)
253 ffc_assert(len(cond) == 1 and firstkey(cond) == (),
254 "Condtion should only be one function: " + repr(cond))
--> 255 ffc_assert(len(true) == 1 and firstkey(true) == (),
256 "True value of Condtional should only be one function: " + repr(true))
257 ffc_assert(len(false) == 1 and firstkey(false) == (),
/usr/lib/python3/dist-packages/ffc/log.py in ffc_assert(condition, *message)
44 def ffc_assert(condition, *message):
45 "Assert that condition is true and otherwise issue an error with given message."
---> 46 condition or error(*message)
47
48
/usr/lib/python3/dist-packages/ffc/log.py in <lambda>(*message)
/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):
Exception: True value of Condtional should only be one function: {((1, 'nzc3[k]', 3, 6),): Product([Symbol('FE0_C0[ip][k]', BASIS)])}
Here is the MWE to reproduce the error -
from dolfin import *
import numpy as np
from mshr import *
from ufl import indices
parameters["form_compiler"]["representation"] = 'quadrature'
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)
domain = Rectangle(Point(0, 0),Point(1, 1))
mesh = generate_mesh(domain, 10)
N = FacetNormal(mesh)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
Vector = VectorFunctionSpace(mesh , 'P' , 1)
Tensor = TensorFunctionSpace(mesh , 'P' , 1)
Vs = [Scalar , Scalar]
Space = FunctionSpace(mesh, MixedElement(Q, Q))
cells = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
da = Measure( 'ds' , domain=mesh , subdomain_data=facets)
dv = Measure( 'dx' , domain=mesh , subdomain_data=cells)
left = CompiledSubDomain( "near(x[0] , l ) && on_boundary" , l=0.)
right = CompiledSubDomain( "near(x[0] , l ) && on_boundary" , l=1.)
facets.set_all(0)
bc = [ DirichletBC(Space.sub(0) , 0.0 , left) , DirichletBC(Space.sub(0) , 0.1 ,right)]
dunkn = TrialFunction(Space)
test = TestFunction(Space)
del_phi , del_T = split(test)
Tref = 300.
Tamb = Tref
t_alpha = 3.9*(10**-3)
varsigma0 = 5.8*(10**+7)
rho = 8960.
kappa = 400.0
capacity = 390.
h = 10.
tMax = 2500.0
Dt = 50.0
t = 0.0
unkn = Function(Space)
unkn0 = Function(Space)
unkn_init = Expression(( '0.0' , 'T_ini' ) , degree = 0, T_ini=Tref )
unkn0 = interpolate(unkn_init , Space)
unkn = interpolate(unkn_init , Space)
phi , T = split(unkn)
phi0 , T0 = split(unkn0)
i,j,k,l = indices(4)
delta = Identity(3)
x, y = SpatialCoordinate(mesh)
condition = lt(y, 0.5)
below = 1.*T
above = 1./T
varsigma = conditional(condition, below, above)
J = -varsigma*grad(phi)
q = -kappa*grad(T)
r = Constant(0.0)
Form = (-J[i]*del_phi.dx(i) + rho*capacity*(T - T0)/Dt*(del_T/T)\
- q[i]*(del_T/T).dx(i) - rho*r*(del_T/T) + J[i]*phi.dx(i)*(del_T/T))*dv\
+ h*(T - Tamb)*(del_T/T)*da
Gain = derivative(Form , unkn , dunkn)
solver_parameters = {
"newton_solver": {
"linear_solver": "umfpack",
"relative_tolerance": 1e-5
}
}
form_compiler_parameters = {
"cpp_optimize": True,
"representation": "quadrature",
"quadrature_degree": 2
}
while t < tMax :
solve(Form == 0,unkn,bc,J=Gain ,solver_parameters = solver_parameters,form_compiler_parameters = form_compiler_parameters)
unkn0.assign(unkn)
t += Dt
Any further suggestions would be very much helpful.