Dear all,


Thank you, the code seems to work fine, which is also good news for Firedrake since my new team here (at INRIA Bordeaux, France) is now interested to use Firedrake.

I'd like to verify my results with an analytical solution of the form u(r, theta) = r*exp(i theta), but I'm having some issues that you might be able to understand. As my mesh is cartesian I define

r = sqrt(x[0]*x[0] + x[1]*x[1]), and theta  = atan(x[1]/x[0]) (ensuring that x[0]!=0),

but have an issue when computing

u.interpolate(r*exp(1j*theta))


Traceback (most recent call last):

  File "sol_analytique.py", line 147, in <module>

    u.interpolate(r*exp(1j*theta))

TypeError: unsupported operand type(s) for *: 'complex' and 'Function'


I also tried

u.interpolate(r*exp(complex(0,theta)))

But I get a UFL error:

Traceback (most recent call last):

  File "sol_analytique.py", line 147, in <module>

    u.interpolate(r*exp(complex(0,theta)))

  File "/Users/mmfg/firedrake/src/ufl/ufl/operators.py", line 584, in exp

    return _mathfunction(f, Exp)

  File "/Users/mmfg/firedrake/src/ufl/ufl/operators.py", line 570, in _mathfunction

    f = as_ufl(f)

  File "/Users/mmfg/firedrake/src/ufl/ufl/constantvalue.py", line 417, in as_ufl

    " to any UFL type." % str(expression))

ufl.log.UFLValueError: Invalid type conversion: -1j can not be converted to any UFL type.


It might be a very basic question but how can I deal with complex functions ?
Thank you,
Floriane

De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk>
Envoyé : mercredi 10 octobre 2018 15:41
À : firedrake
Objet : Re: [firedrake] local discontinuity
 

Yes. Those are equivalent. You should be able to work out which of your results is correct by just integrating the unit function along the line and observing which solution approximates pi*d

 

From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk>
Reply-To: firedrake <firedrake@imperial.ac.uk>
Date: Wednesday, 10 October 2018 at 15:39
To: firedrake <firedrake@imperial.ac.uk>
Subject: Re: [firedrake] local discontinuity

 

Dear all,

 

Concerning the integral on the interface between the disk and the plane, is it not equivalent to write

 

*4*avg(I_d)*avg(I_c)*dS

and

*dS(interface_id)

 

where interface_id is the id of the interface ?

I don't get the same results when using one or the other.

 

Thank you,

Floriane

 

 

De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk>
Envoyé : lundi 8 octobre 2018 14:48
À : firedrake
Objet : Re: [firedrake] local discontinuity

 

Note also that using avg in a dx integral doesn’t make much sense (thought it doesn’t actually change the answer)

 

From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk>
Reply-To: firedrake <firedrake@imperial.ac.uk>
Date: Monday, 8 October 2018 at 14:43
To: firedrake <firedrake@imperial.ac.uk>
Subject: Re: [firedrake] local discontinuity

 

Dear Thomas,

 

That should indeed be dS and not dx, thank you!

 

Floriane

 

De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Gibson, Thomas <t.gibson15@imperial.ac.uk>
Envoyé : lundi 8 octobre 2018 14:34
À : firedrake
Objet : Re: [firedrake] local discontinuity

 

Dear Floriane,

 

The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker.

 

Best regards,

- Thomas

From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk>
Sent: 08 October 2018 14:18:09
To: firedrake
Subject: Re: [firedrake] local discontinuity

 

Dear David,

 

I tried what you suggested and I have the following error when solving the Laplace equation:

ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal.


I defined:

V_d = FunctionSpace(mesh, "CG", 1)  

V_p = FunctionSpace(mesh, "CG", 1)  

W = V_d*V_p

 

trial_d, trial_p = TrialFunction(W)

v_d, v_p = TestFunction(W)

 

result_mixed = Function(W)

u_d, u_p = split(result_mixed)

 

I define the normal pointing outward the disk interface as:

n_vec = FacetNormal(mesh)

n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-")

 

where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows:

 

A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d

A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p

 

L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface 

L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) -   sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx  # right BC and interface 

 

A_u = A_u_d + A_u_p

L_u = L_u_d + L_u_p

 

cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed)

cell_solver = LinearVariationalSolver(cell_pb)

 

cell_solver.solve()

 

The full error is pasted below. Do you know where it comes from ?

Thank you,

Floriane

 

UFL:ERROR Integral of type cell cannot contain a ReferenceNormal.

Traceback (most recent call last):

  File "lin_cell.py", line 196, in <module>

    cell_solver = LinearVariationalSolver(cell_pb)

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__

    super(LinearVariationalSolver, self).__init__(*args, **kwargs)

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__

    pre_function_callback=pre_f_callback)

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__

    form_compiler_parameters=fcp)

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable

    collect_loops=True)

  File "<decorator-gen-280>", line 2, in _assemble

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper

    return f(*args, **kwargs)

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble

    kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse)

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form

    number_map).kernels

  File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__

    obj = make_obj()

  File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj

    obj.__init__(*args, **kwargs)

  File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__

    tree = tsfc_compile_form(form, prefix=name, parameters=parameters)

  File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form

    fd = ufl_utils.compute_form_data(form)

  File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data

    do_estimate_degrees=do_estimate_degrees,

  File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data

    _check_facet_geometry(self.integral_data)

  File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry

    error("Integral of type %s cannot contain a %s." % (it, cls.__name__))

  File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error

    raise self._exception_type(self._format_raw(*message))

ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal.

 

 

De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk>
Envoyé : vendredi 5 octobre 2018 15:56
À : firedrake
Objet : Re: [firedrake] local discontinuity

 

Dear Floriane,

 

I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction.

 

We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk.

 

We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d.

 

You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain.

 

Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the  *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d).

 

If

W = V_d * V_p

and

u_d, u_p = TrialFunctions(W)

 

then you can write the jump as:

 

avg(u_d) – avg(u_p)

 

The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)).

 

Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html


www.firedrakeproject.org

Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake.

 

 

I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications.

 

Regards,

 

David

 

 

From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk>
Reply-To: firedrake <firedrake@imperial.ac.uk>
Date: Friday, 5 October 2018 at 15:33
To: firedrake <firedrake@imperial.ac.uk>
Subject: [firedrake] local discontinuity

 

Dear Firedrakers,

 

I would like to solve equations in a domain made of two subdomains:

- a disk, inside which test functions are continuous (Omega_d)

- a plane around the disk, in which test functions are continuous (Omega_p)

However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane.

This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is

Delta u = 0 in Omega_d

Delta u = 0 in Omega_p

partial_t [u] = ...   on Gamma, where [u]= u(+)-u(-) is the jump

 

Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ?

 

Thank you for your help,

 

Floriane