New subject: Specify boundary conditions at a corner

25 Aug 2015

      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@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@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@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@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@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@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@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@imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

Re: [firedrake] Specify boundary conditions at a corner

David Ham

Justin Chang

tags

participants (2)