[firedrake] Specify boundary conditions at a corner

[Justin Chang]

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
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Test_LSAug.py
Type: text/x-python-script
Size: 2193 bytes
Desc: not available
URL: <http://mailman.ic.ac.uk/pipermail/firedrake/attachments/20150826/f2fbc792/attachment.bin>