Hi everybody.
I am writing a stabilized Navier Stokes Solver, and as stabilizing parameter tau I chose this element dependent function:
tau = ( v* (G* v) + c * ni**2 * ( G:G) ) ** -0.5
where:
v is the velocity
c,ni constants (not element-dependent)
G is the metric tensor, defined as:
G = J^T * J
being J the Jacobi Matrix of the element.
I implemented this part
# Stabilizing parameter
c = 9.0;
ELEM_cell = mesh.ufl_cell()
Elem_Jac = Elem_cell.J
Elem_Metric = dot( transpose(Elem_Jac) , Elem_Jac)
tau = (ni**2 * c * inner(Elem_Metric,Elem_Metric) )**(-0.5)
When defining the form, tau is multiplying some terms
RESIDUAL = tau*dot(grad(p), grad(q))*dx + [OMITTED]
After that, define a system Jacobian matrix (not to be confused with the Element's Jacobi), and the nonlinear variational problem.
JAC = derivative(RESIDUAL,vp , dvp)
NLproblem = NonlinearVariationalProblem(RESIDUAL,vp, bcs, JAC)
NLsolver = NonlinearVariationalSolver(NLproblem)
However i get:
Visiting GeometricQuantity: GeometryJacobi(Cell('triangle', 2))
This type of GeometricQuantity is not supported (yet).
Traceback (most recent call last):
File "Leaky_Cavity_NS_Steady.py", line 221, in <module>
NLproblem = NonlinearVariationalProblem(RES,vp, bcs, JAC)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 116, in __init__
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/form.py", line 56, in __init__
common_cell)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 60, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 122, in jit
return jit_compile(form, parameters=p, common_cell=common_cell)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 78, in jit
return jit_form(ufl_object, parameters, common_cell)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 208, in jit_form
parameters=parameters)
File "/usr/lib/python2.7/dist-packages/ffc/compiler.py", line 156, in compile_form
ir = compute_ir(analysis, parameters)
File "/usr/lib/python2.7/dist-packages/ffc/representation.py", line 96, in compute_ir
for (i, fd) in enumerate(form_datas)]
File "/usr/lib/python2.7/dist-packages/ffc/representation.py", line 216, in _compute_integral_ir
parameters)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturerepresentation.py", line 87, in compute_integral_ir
itg_data.domain_type, form_data.cell)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturerepresentation.py", line 145, in _transform_integrals_by_type
terms = _transform_integrals(transformer, integrals_dict, domain_type)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturerepresentation.py", line 442, in _transform_integrals
terms = transformer.generate_terms(integral.integrand(), domain_type)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 774, in generate_terms
terms = self.visit(integrand)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformer.py", line 251, in power
base_code = self.visit(base)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 552, in index_sum
ops.append(self.visit(summand))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 552, in index_sum
ops.append(self.visit(summand))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 515, in indexed
code = self.visit(indexed_expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 728, in component_tensor
code = self.visit(component_expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 552, in index_sum
ops.append(self.visit(summand))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 515, in indexed
code = self.visit(indexed_expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 728, in component_tensor
code = self.visit(component_expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 515, in indexed
code = self.visit(indexed_expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 105, in visit
r = h(o)
File "/usr/lib/python2.7/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 223, in geometric_quantity
error("This type of GeometricQuantity is not supported (yet).")
File "<string>", line 1, in <lambda>
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 154, in error
raise self._exception_type(self._format_raw(*message))
Exception: This type of GeometricQuantity is not supported (yet).
It seems then that this feature is not yet implemented.
Is there any other way to be able to "extract" the element Jacobi? For example, I saw a "get_vertex_coordinates()" in the class "dolfin.cpp.mesh.Cell", which is an ingredient (as long as the derivatives of the basis functions in respect to the isoparametric coordinates) to build the element's Jacobi.
I ask here because I am not an expert at all.
Thank you in advance,
R.
