[firedrake] Specify boundary conditions at a corner

[Justin Chang]

David,

Okay yes this issue is resolve now, thank you very much.

Lawrence,

So if I use a GMSH file, the third argument in DirichletBC(...) would
correspond to whatever the Physical IDS are? Would it correspond to
"Face set" values if I were to use an exodusii file?

Thanks,
Justin

On Wed, Aug 26, 2015 at 9:17 AM, David Ham <David.Ham at imperial.ac.uk> wrote:
> Hi Justin,
>
> Thanks for that, you uncovered a bug in some of how our internal imports
> work. I have just merged a fix into Firedrake master so hopefully this
> should now work.
>
> Cheers,
>
> David
>
> On Wed, 26 Aug 2015 at 10:29 Justin Chang <jychang48 at gmail.com> wrote:
>>
>> David,
>>
>> I am still getting the same error. Attached the code, can you have a
>> look at it and see what's going on? Note that if you comment out the
>> point dirichlet stuff and uncomment the nullspace related operations,
>> it outputs the correct solution.
>>
>> Thanks,
>> Justin
>>
>> On Tue, Aug 25, 2015 at 3:03 AM, David Ham <David.Ham at imperial.ac.uk>
>> wrote:
>> > Hi Justin,
>> >
>> > The imports in bcs.py are relative imports (which is slightly bad style
>> > but...) In order to access correct module from outside you need to use:
>> >
>> > from firedrake import utils
>> >
>> > Regards,
>> >
>> > David
>> >
>> > On Tue, 25 Aug 2015 at 00:48 Justin Chang <jychang48 at gmail.com> wrote:
>> >>
>> >> David,
>> >>
>> >> I get the following error:
>> >>
>> >> Traceback (most recent call last):
>> >>   File "Test_LSAug.py", line 41, in <module>
>> >>     class PointDirichletBC(DirichletBC):
>> >>   File "Test_LSAug.py", line 42, in PointDirichletBC
>> >>     @utils.cached_property
>> >> AttributeError: 'module' object has no attribute 'cached_property'
>> >>
>> >> When I comment out that line, i get a long error. When I try to import
>> >> utils and all the other modules that bcs.py imports, i get errors
>> >> saying those modules don't exist. Know what's going on?
>> >>
>> >> Thanks,
>> >> Justin
>> >>
>> >> On Mon, Aug 24, 2015 at 3:14 AM, David Ham <David.Ham at imperial.ac.uk>
>> >> wrote:
>> >> > Hi Justin,
>> >> >
>> >> > Having subclassed DirichletBC, what you now have is a new type of
>> >> > Dirichlet
>> >> > "boundary" object, so you just use it as a Dirichlet condition:
>> >> >
>> >> > bc3 = PointDirichletBC(W.sub(1), 0, 0)
>> >> >
>> >> > bc_all = [bc1, bc2, bc3]
>> >> >
>> >> > However Andrew is right, removing the nullspace is almost certainly a
>> >> > better
>> >> > solution in this case.
>> >> >
>> >> > Cheers,
>> >> >
>> >> > David
>> >> >
>> >> > On Sun, 23 Aug 2015 at 05:23 Justin Chang <jychang48 at gmail.com>
>> >> > wrote:
>> >> >>
>> >> >> Andrew,
>> >> >>
>> >> >> That actually did the trick, thank you very much.
>> >> >>
>> >> >> But I still would like to know the answer to my original question.
>> >> >> Or
>> >> >> perhaps, instead of a single corner, I would like to set (x=0,y<0.1)
>> >> >> and (x<0.1,y=0) to a specific value while everywhere else is set to
>> >> >> a
>> >> >> different value.
>> >> >>
>> >> >> Thanks,
>> >> >> Justin
>> >> >>
>> >> >> On Sat, Aug 22, 2015 at 11:58 AM, Andrew McRae
>> >> >> <a.mcrae12 at imperial.ac.uk>
>> >> >> wrote:
>> >> >> > This doesn't answer your main question, but setting an appropriate
>> >> >> > nullspace
>> >> >> > might be more appropriate than pinning a single value; see
>> >> >> > http://www.firedrakeproject.org/solving-interface.html
>> >> >> >
>> >> >> > [If I'm wrong, I'll let the usual suspects correct my
>> >> >> > misinformation]
>> >> >> >
>> >> >> > On 22 August 2015 at 18:47, Justin Chang <jychang48 at gmail.com>
>> >> >> > wrote:
>> >> >> >>
>> >> >> >> David,
>> >> >> >>
>> >> >> >> How exactly do I use that class in my code? Say I have the
>> >> >> >> following
>> >> >> >> function spaces/discretization based on the least-squares finite
>> >> >> >> element method:
>> >> >> >>
>> >> >> >> # Mesh
>> >> >> >> mesh = UnitSquareMesh(seed,seed)
>> >> >> >> V = VectorFunctionSpace(mesh,"CG",1)
>> >> >> >> Q = FunctionSpace(mesh,"CG",1)
>> >> >> >> W = V*Q
>> >> >> >> v,p = TrialFunctions(W)
>> >> >> >> w,q = TestFunctions(W)
>> >> >> >>
>> >> >> >> # Weak form
>> >> >> >> g = Function(V)
>> >> >> >>
>> >> >> >>
>> >> >> >>
>> >> >> >>
>> >> >> >> g.interpolate(Expression(("cos(pi*x[0])*sin(pi*x[1])+2*pi*cos(2*pi*x[0])*sin(2*pi*x[1])","-sin(pi*x[0])*cos(pi*x[1])+2*pi*sin(2*pi*x[0])*cos(2*pi*x[1])")))
>> >> >> >> a = dot(v+grad(p),w+grad(q))*dx + div(v)*div(w)*dx
>> >> >> >> L = dot(w+grad(q),g)*dx
>> >> >> >>
>> >> >> >> # Boundary conditions
>> >> >> >> bc1 = DirichletBC(W.sub(0).sub(0),
>> >> >> >> Expression("cos(pi*x[0])*sin(pi*x[1])"), (1,2))
>> >> >> >> bc2 = DirichletBC(W.sub(0).sub(1),
>> >> >> >> Expression("-sin(pi*x[0])*cos(pi*x[1])"), (3,4))
>> >> >> >> bc_all = [bc1,bc2]
>> >> >> >>
>> >> >> >> # Solve
>> >> >> >> cg_parameters = {
>> >> >> >>     'ksp_type': 'cg',
>> >> >> >>     'pc_type': 'bjacobi'
>> >> >> >>     }
>> >> >> >> solution = Function(W)
>> >> >> >> A = assemble(a,bcs=bc_all)
>> >> >> >> b = assemble(L,bcs=bc_all)
>> >> >> >> solver =
>> >> >> >>
>> >> >> >> LinearSolver(A,solver_parameters=cg_parameters,options_prefix="cg_")
>> >> >> >> solver.solve(solution,b)
>> >> >> >>
>> >> >> >> If I run the above code I get an error saying 'LinearSolver
>> >> >> >> failed
>> >> >> >> to
>> >> >> >> converge after %d iterations with reason: %s', 196,
>> >> >> >> 'DIVERGED_INDEFINITE_MAT'. Which I am guessing has to do with the
>> >> >> >> lack
>> >> >> >> of a boundary condition for the Q space, thus I want to ensure a
>> >> >> >> unique solution by prescribing the bottom left constraint to a
>> >> >> >> zero
>> >> >> >> value.
>> >> >> >>
>> >> >> >> Thanks,
>> >> >> >> Justin
>> >> >> >>
>> >> >> >>
>> >> >> >> On Fri, Aug 21, 2015 at 4:52 AM, David Ham
>> >> >> >> <David.Ham at imperial.ac.uk>
>> >> >> >> wrote:
>> >> >> >> > Hi Justin,
>> >> >> >> >
>> >> >> >> > The nice way of doing this would require better subdomain
>> >> >> >> > support
>> >> >> >> > than
>> >> >> >> > we
>> >> >> >> > have right now. However there is a slightly hacky way of doing
>> >> >> >> > it
>> >> >> >> > which
>> >> >> >> > I
>> >> >> >> > think will cover your case nicely.
>> >> >> >> >
>> >> >> >> > If you take a look at the DirichletBC class (in bcs.py), you'll
>> >> >> >> > notice
>> >> >> >> > that
>> >> >> >> > the set of nodes at which the BC should be applied is
>> >> >> >> > calculated
>> >> >> >> > in
>> >> >> >> > DirichletBC.nodes . So you could simply subclass DirichletBC
>> >> >> >> > and
>> >> >> >> > replace
>> >> >> >> > nodes with a function which returns the index of the zero node.
>> >> >> >> > For
>> >> >> >> > example
>> >> >> >> > (and I confess this is a sketch code which I haven't tried to
>> >> >> >> > run):
>> >> >> >> >
>> >> >> >> > class PointDirichletBC(DirichletBC):
>> >> >> >> >     @utils.cached_property
>> >> >> >> >     def nodes(self):
>> >> >> >> >         # Find the array of coordinate values.
>> >> >> >> >         x =
>> >> >> >> > self.function_space().mesh().coordinates.dat.data_ro
>> >> >> >> >         # Find the location of the zero rows in that
>> >> >> >> >         return np.where(~x.any(axis=1))[0]
>> >> >> >> >
>> >> >> >> > Does that work for you?
>> >> >> >> >
>> >> >> >> > Cheers,
>> >> >> >> >
>> >> >> >> > David
>> >> >> >> >
>> >> >> >> > On Fri, 21 Aug 2015 at 03:32 Justin Chang <jychang48 at gmail.com>
>> >> >> >> > wrote:
>> >> >> >> >>
>> >> >> >> >> Hi all,
>> >> >> >> >>
>> >> >> >> >> If I create a mesh using UnitSquareMesh or UnitCubeMesh, is
>> >> >> >> >> there
>> >> >> >> >> a
>> >> >> >> >> way to subject a single point (as opposed to an entire
>> >> >> >> >> edge/face)
>> >> >> >> >> to
>> >> >> >> >> a
>> >> >> >> >> DirichletBC? I want to subject the the location x=0,y=0 to
>> >> >> >> >> some
>> >> >> >> >> value.
>> >> >> >> >>
>> >> >> >> >> Thanks,
>> >> >> >> >> Justin
>> >> >> >> >>
>> >> >> >> >> _______________________________________________
>> >> >> >> >> firedrake mailing list
>> >> >> >> >> firedrake at imperial.ac.uk
>> >> >> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >> >> >> >
>> >> >> >> >
>> >> >> >> > _______________________________________________
>> >> >> >> > firedrake mailing list
>> >> >> >> > firedrake at imperial.ac.uk
>> >> >> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >> >> >> >
>> >> >> >>
>> >> >> >> _______________________________________________
>> >> >> >> firedrake mailing list
>> >> >> >> firedrake at imperial.ac.uk
>> >> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >> >> >
>> >> >> >
>> >> >> >
>> >> >> > _______________________________________________
>> >> >> > firedrake mailing list
>> >> >> > firedrake at imperial.ac.uk
>> >> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >> >> >
>> >> >>
>> >> >> _______________________________________________
>> >> >> firedrake mailing list
>> >> >> firedrake at imperial.ac.uk
>> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >> >
>> >> >
>> >> > _______________________________________________
>> >> > firedrake mailing list
>> >> > firedrake at imperial.ac.uk
>> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >> >
>> >>
>> >> _______________________________________________
>> >> firedrake mailing list
>> >> firedrake at imperial.ac.uk
>> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >
>> >
>> > _______________________________________________
>> > firedrake mailing list
>> > firedrake at imperial.ac.uk
>> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>> >
>
>
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>