Hi Tomasz, There are two things here. The first one is cell subdomains. There is a pull request in at PETSc right now which will make this work to some extent. We'll let you know once the PETSc people accept that pull. The other issue is what you can pass to DirichletBC. That is always a boundary marker, not a cell subdomain marker - so the precise syntax you used there won't quite work. The functionality will hopefully be there soon though. On Tue, 24 May 2016 at 09:44 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Dear all,
Can I invoke in Firedrake a subdomain id in a mesh generated by gmsh?
It would be useful for me in constructing indicator functions.
The line:
...
bc = DirichletBC( V, 1, 101 )
...
causes an error:
LookupError: 101 is not a valid marker
although 101 is defined in gmsh as a Physical Volume. I had no problems with Physical Surface in the same context.
Regards,
Tomasz ------------------------------ *From:* firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> *Sent:* 27 April 2016 17:32:56
*To:* firedrake *Subject:* Re: [firedrake] Laplace equation on subdomain
That did the trick, I have a coupled solver working now, thank you!
Tomasz ------------------------------ *From:* firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of David Ham <David.Ham@imperial.ac.uk> *Sent:* 26 April 2016 22:29:44 *To:* firedrake *Subject:* Re: [firedrake] Laplace equation on subdomain Well the following reproduces the reported error for me:
from firedrake import * m = UnitSquareMesh(2, 2) fs = FunctionSpace(m, "CG", 1) v = TestFunction(fs) assemble(v*dS)
On Tue, 26 Apr 2016 at 21:27 Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
I didn't think this was true, unless something changed "recently".
I can't test right now, but I'm quite confident I used to use continuous functions in interior facet integrals on simplex meshes without needing a restriction.
However, I recall being forced to use restrictions on extruded meshes; I never got to the bottom of why. Perhaps Miklós can comment (or state that my recollection is nonsense)?
On 26 April 2016 at 21:12, David Ham <David.Ham@imperial.ac.uk> wrote:
Hi Tomasz,
This is a subtlety of UFL. *All* arguments and coefficients in an interior facet integral have to be "restricted". This even applies to ones which are continuous and for which the "+" and "-" restrictions are therefore the same. You can just arbitrarily decide to restrict the continuous things to the plus side, or take avg() of the lack of symmetry offends your sense of aesthetics.
Regards,
David
On Tue, 26 Apr 2016 at 15:46 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Hi,
Previous problem continued.
I have a working solver for separate subsystems in mixed function space now. I'd like to introduce the coupling terms at the interface, but they cause a UFL exception "Discontinuous type Argument must be restricted."
Full code is here: https://bitbucket.org/mmtjs/lin_coupled_fd_2d
<https://bitbucket.org/mmtjs/lin_coupled_fd_2d>
The relevant part from lib/coupled_system.py is:
...
# define the normal vector at the interface
self.n = fd.FacetNormal(self.mesh) self.n_int = (self.B.I("+")*self.W.I("-") - self.B.I("-")*self.W.I("+"))*self.n("+") # btw how can write this to file to verify?
...
# the second line causes exception (and later similar terms too)
a_U = fd.dot( self.trial_s * self.B.I, self.v_s * self.B.I) * fd.dx \ + fd.dot( self.v_s, self.n_int ) * self.trial_f * fd.dS
Regards,
Tomasz
------------------------------ *From:* firedrake-bounces@imperial.ac.uk < firedrake-bounces@imperial.ac.uk> on behalf of Tomasz Salwa [RPG] < mmtjs@leeds.ac.uk> *Sent:* 21 April 2016 11:04
*To:* firedrake *Subject:* Re: [firedrake] Laplace equation on subdomain
Cool, works now, dynamic case too. I'll do the same for the structural part.
Cheers,
Tomasz
------------------------------ *From:* firedrake-bounces@imperial.ac.uk < firedrake-bounces@imperial.ac.uk> on behalf of David Ham < David.Ham@imperial.ac.uk> *Sent:* 21 April 2016 09:38 *To:* firedrake *Subject:* Re: [firedrake] Laplace equation on subdomain Hi Tomasz,
See the new pull request. In essence, the projection is not safe in the indicator function case.
Regards,
David
On Wed, 20 Apr 2016 at 15:53 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Hi David,
Thank you, now it converges properly.
However, I noticed that the error depends on the mesh resolution beyond the fluid subdomain, which I think should not be the case (at least when I refine it in x direction). Perhaps the solution is ok, but the source is in the way I impose the analytical one, or compute the error. Can you diagnose this?
Bitbucket's repo is updated.
Regards,
Tomasz
------------------------------ *From:* firedrake-bounces@imperial.ac.uk < firedrake-bounces@imperial.ac.uk> on behalf of David Ham < David.Ham@imperial.ac.uk> *Sent:* 20 April 2016 10:05
*To:* firedrake *Subject:* Re: [firedrake] Laplace equation on subdomain Hi Tomasz,
I have just put in a pull request fixing your problem. In essence, you were only applying half of the indicator function fix. In order to correctly use indicator functions to solve problems over a subdomain, it is necessary to *both* use the custom boundary condition to force the solution in the rest of the domain to zero *and* multiply the integrands in your forms by the indicator function. You were only doing the former and not the latter. Making this fix causes the error to drop to 10^-5. I haven't done a convergence test to check that the error converges at the right rate, but you should really do so.
Regards,
David
On Tue, 19 Apr 2016 at 11:52 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Hi David,
Here's the code:
https://bitbucket.org/mmtjs/laplace_extended_domain
It's roughly the same as before, just with a gmsh-generated mesh.
Regards,
Tomasz
------------------------------ *From:* firedrake-bounces@imperial.ac.uk < firedrake-bounces@imperial.ac.uk> on behalf of David Ham < David.Ham@imperial.ac.uk> *Sent:* 19 April 2016 09:54 *To:* firedrake
*Subject:* Re: [firedrake] Laplace equation on subdomain Hi Tomasz,
Can you please post the code you are using for this (preferably by pointing us to a git repo we can clone). I can't diagnose what is going wrong just by looking at the visualisation, I'm afraid.
Regards,
David
On Mon, 18 Apr 2016 at 12:12 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
> Hi, > > I plotted the difference between analytical and numerical solutions > (see attached). It doesn't converge to the right solution with mesh > refinement. Is it possible that the Neumann BC is not applied at the > interface during the computation, although the extended Dirichlet BCs > enforce resulting zero solution beyond the fluid subdomain? > > Cheers, > Tomasz > > ________________________________________ > From: firedrake-bounces@imperial.ac.uk < > firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell < > lawrence.mitchell@imperial.ac.uk> > Sent: 15 April 2016 16:00 > To: firedrake@imperial.ac.uk > Subject: Re: [firedrake] Laplace equation on subdomain > > > On 15 Apr 2016, at 15:54, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> > wrote: > > > > Hi David, > > > > It works partially, i.e. the nodes beyond the fluid subdomain are > "turned off" now, but it looks like the interface is "smeared". I suspect > that the assembled matrix's entries involve also integration over elements > one layer beyond the subdomain boundary. Is it possible to exclude them? > > If you have neumann conditions on your subdomain boundary you've > only constrained the gradient of the solution. Therefore on the boundary > nodes you may have non-zero solution values. > > Because the domain you visualise in paraview extends outside the > subdomain you'll have (in 1D) something like: > > > ........ > \............ > ^ boundary node > > Note how the cell at the boundary has a non-zero coefficient value > are one end, but a zero value at the other (assuming piecewise liners). > Paraview will linearly interpolate this basis function when visualising, > which may be the artefact you're seeing. > > Can you construct an analytic solution and determine that you > converge at the correct rate as you refine the mesh? This is a safer test > than eyeballing the solution. > > Lawrence >