Adaptive solve does not work with both cell domains and facet domains

I'm not quite sure how to best resolve this. (BTW: it doesn't really
have anything to do with both cell domains and facet domains. Try
removing the facet domains)

For the error control, a bunch of forms are created, in particular, one
that looks more or less like

     # u = Function(...)
     # v = L's TestFunction
     # R_T = Function(...)

     L_dT = L - action(a, u) - inner(R_T, v)*dx

Here, the last term (let's call it 'L_R') is supposed to be an integral
over the entire domain. However, UFL complains (with good reason) since
the domain data in the cell integral in 'a' (cell_markers) and in 'L_R'
(None) doesn't match.

Possible solutions:

(1) Redesign error control generation so that all terms in automatically
created forms are stored and assembled separately and added on the
linear algebra level. For instance L_1 = L, L_2 = - action(a, u), L_3 =
-inner(R_T, v)*dx in the example above. (Will become a mess and involves
a significant amount of work, but is very robust)

(2) Hack the generation of 'L_R', so that -inner(R_T, v)*dx(i) is added
for i in range-of-values-used-in-a-and-or-L-to-indicate-subdomains

(3) Keep rethinking a more flexible way of specifying subdomains and
assembly over such.

(4) ... ?

--
Marie

On 10/11/11 06:00, Johan Hake wrote:
> Public bug reported:
>
> When I define both cell and facet domains in an equation I ship to
> adaptive solve I get a weird error:
>
> ###################################
> from dolfin import *
>
> mesh = UnitSquare(2,2)
> face_markers = FacetFunction("uint", mesh, 1)
> cell_markers = CellFunction("uint", mesh, 1)
> V = FunctionSpace(mesh, "CG", 1)
> v = TestFunction(V)
> u = TrialFunction(V)
>
> dxx = dx[cell_markers]
> dss = ds[face_markers]
>
> a = u*v*dxx(1) + u*v*dss(1)
> L = v*dss(1)
>
> u = Function(V)
>
> solve(a==L, u, tol=0.001, M=u*dxx(1))
>
> ###############################
>
> UFLException: Found two domain data objects for same domain type.
>
> ** Affects: dolfin
> Importance: Undecided
> Status: New
>

On 11 October 2011 16:24, Marie Rognes <email address hidden> wrote:
> I'm not quite sure how to best resolve this. (BTW: it doesn't really
> have anything to do with both cell domains and facet domains. Try
> removing the facet domains)
>
> For the error control, a bunch of forms are created, in particular, one
> that looks more or less like
>
>     # u = Function(...)
>     # v = L's TestFunction
>     # R_T = Function(...)
>
>     L_dT = L - action(a, u) - inner(R_T, v)*dx
>
> Here, the last term (let's call it 'L_R') is supposed to be an integral
> over the entire domain.

Which it is not, even ignoring the bug reported by Johan here... It is
an integral over \Omega_0 not \Omega, as we've talked about before.
Not very nice, but that's what it is now anyway.

> However, UFL complains (with good reason) since
> the domain data in the cell integral in 'a' (cell_markers) and in 'L_R'
> (None) doesn't match.
>
> Possible solutions:
>
> (1) Redesign error control generation so that all terms in automatically
> created forms are stored and assembled separately and added on the
> linear algebra level. For instance L_1 = L, L_2 = - action(a, u), L_3 =
> -inner(R_T, v)*dx in the example above. (Will become a mess and involves
> a significant amount of work, but is very robust)
>

I don't think you can make this very robust without some work on the
entire subdomain modeling, starting with UFL and propagating through
other components.

> (2) Hack the generation of 'L_R', so that -inner(R_T, v)*dx(i) is added
> for i in range-of-values-used-in-a-and-or-L-to-indicate-subdomains
>

Mhm. Actually, more like
for i in range-of-values-used-in-mesh-to-indicate-subdomains,
which is information only available in the pydolfin context.

> (3) Keep rethinking a more flexible way of specifying subdomains and
> assembly over such.
>

Yes, basically what I wrote above.

> (4) ... ?

(5) Profit?

:)

Martin

>
> --
> Marie
>
> On 10/11/11 06:00, Johan Hake wrote:
>> Public bug reported:
>>
>> When I define both cell and facet domains in an equation I ship to
>> adaptive solve I get a weird error:
>>
>> ###################################
>> from dolfin import *
>>
>> mesh = UnitSquare(2,2)
>> face_markers = FacetFunction("uint", mesh, 1)
>> cell_markers = CellFunction("uint", mesh, 1)
>> V = FunctionSpace(mesh, "CG", 1)
>> v = TestFunction(V)
>> u = TrialFunction(V)
>>
>> dxx = dx[cell_markers]
>> dss = ds[face_markers]
>>
>> a = u*v*dxx(1) + u*v*dss(1)
>> L = v*dss(1)
>>
>> u = Function(V)
>>
>> solve(a==L, u, tol=0.001, M=u*dxx(1))
>>
>> ###############################
>>
>> UFLException: Found two domain data objects for same domain type.
>>
>> ** Affects: dolfin
>>       Importance: Undecided
>>           Status: New
>>
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/872105
>
> Title:
>  Adaptive solve does not work with both cell domains and facet domains
>
> Status in DOLFIN:
>  Confirmed
>
> Bug description:
>  When I define both cell and facet domains in an equation I ship to
>  adaptive solve I get a weird error:
>
>  ############################...

On Tuesday October 11 2011 07:24:11 Marie Rognes wrote:
> I'm not quite sure how to best resolve this. (BTW: it doesn't really
> have anything to do with both cell domains and facet domains. Try
> removing the facet domains)

Ok, so it is all about cell domains then. Feel free to change the name of the
bug then.

> For the error control, a bunch of forms are created, in particular, one
> that looks more or less like
>
> # u = Function(...)
> # v = L's TestFunction
> # R_T = Function(...)
>
> L_dT = L - action(a, u) - inner(R_T, v)*dx
>
> Here, the last term (let's call it 'L_R') is supposed to be an integral
> over the entire domain. However, UFL complains (with good reason) since
> the domain data in the cell integral in 'a' (cell_markers) and in 'L_R'
> (None) doesn't match.
>
> Possible solutions:
>
> (1) Redesign error control generation so that all terms in automatically
> created forms are stored and assembled separately and added on the
> linear algebra level. For instance L_1 = L, L_2 = - action(a, u), L_3 =
> -inner(R_T, v)*dx in the example above. (Will become a mess and involves
> a significant amount of work, but is very robust)
>
> (2) Hack the generation of 'L_R', so that -inner(R_T, v)*dx(i) is added
> for i in range-of-values-used-in-a-and-or-L-to-indicate-subdomains
>
> (3) Keep rethinking a more flexible way of specifying subdomains and
> assembly over such.

Is L_R always over the entire domain? Then I think 3 is the best option, where
we create a default domain, which is interpreted to beover the whole domain.
Let say dx(-1). We have talked about this before. It is a handy feature which
can be used outside this special case.

Do not think 2 will do it if L_R should be integrated over the whole domain.
Then you need the union of all domain_i to fill the whole domain, which might
not be the case.

Johan

> (4) ... ?
>
>
> --
> Marie
>
> On 10/11/11 06:00, Johan Hake wrote:
> > Public bug reported:
> >
> > When I define both cell and facet domains in an equation I ship to
> > adaptive solve I get a weird error:
> >
> > ###################################
> > from dolfin import *
> >
> > mesh = UnitSquare(2,2)
> > face_markers = FacetFunction("uint", mesh, 1)
> > cell_markers = CellFunction("uint", mesh, 1)
> > V = FunctionSpace(mesh, "CG", 1)
> > v = TestFunction(V)
> > u = TrialFunction(V)
> >
> > dxx = dx[cell_markers]
> > dss = ds[face_markers]
> >
> > a = u*v*dxx(1) + u*v*dss(1)
> > L = v*dss(1)
> >
> > u = Function(V)
> >
> > solve(a==L, u, tol=0.001, M=u*dxx(1))
> >
> > ###############################
> >
> > UFLException: Found two domain data objects for same domain type.
> >
> > ** Affects: dolfin
> >
> > Importance: Undecided
> >
> > Status: New

Moving this to 1.1, will need infrastructure.

I just bumped into the exact same error, so I guess a fix hasn't quite made it into 1.1.0.
Has anyone been able to come up with a workaround of some sort?

Looks like this has been fixed.

Well, at least the error is no longer there (thanks to Martin's work on domains in UFL). I have no idea about the veracity of the results ...