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