Re: [firedrake] Specify boundary conditions at a corner
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
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
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 21/08/15 11:52, David Ham wrote:
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]
I believe David is fixing the error in using this kind of idea. However, I note in passing that this exact form may not do exactly what you want. It assumes that the function space that you're wanting to apply the boundary condition to is the same space that is holding the coordinate field (namely piecewise linears, P1 on simplices, Q1 on quads). If the space you want to pin is not P1, you'll have to work harder (this is the subdomain problem David alluded to). If you want to have a boundary marker that applies to a small subset of the boundary of a mesh, I recommend, rather than using the built in utility meshes, to build the mesh using Gmsh. There you can provide physical IDs to whatever subset of the boundary you like and then use these ID numbers in your boundary condition objects. Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJV3bdwAAoJECOc1kQ8PEYvJbEH/j3pgwIdWig5IHQqVGSI+OUf xsu1XZMGhf4MCRzXYlqlQQZ3HGNTO7tCOXHN/DATeKW4ZYB21bSrBQEsPUr3P+iL EteoaXYv5nDbMBJzQcxnXI/VQyPRhVbSd2mORK8pf/6xyVUceVrBgBCxaG5hbI8e Axubo+ugrg/32NkX1CXISaMCn7fdF2qkZKJA2YZhoZqmiW0v1n3TiYuzDUlhRoC+ vrdAVUOCfyEvU/5Q/GMSlq76YVvGIBgwcpnzvUYK3gPPnNrj1PCMeNty1LyHqThJ zPZX1ijUfXtjsz37LkDJxu4bWUkD5tXJfKUJXex3YfoqIYpyRiDs8IPBRcXsEl0= =1qnA -----END PGP SIGNATURE-----
participants (3)
-
David Ham
-
Justin Chang
-
Lawrence Mitchell