[firedrake] Specify boundary conditions at a corner

[Justin Chang]

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
>