1) Does this region line up with cells of your mesh? If so, my normal way
is to interpolate that expression into a P0 function:
P0 = Function(mesh, "DG", 0)
f = Function(P0)
f.interpolate(Expression("(x[0] >= 0.25 && x[0] <= 0.30 && x[1] >= 0.25 &&
x[1] <= 0.30) ? 1.0 : 0.0"))
and then use f in calculations.
2) Pass, someone else can answer this
3) "Yes", but can you be more specific? SUPG is just modifying the test
function, so my code looked something like
...
u = ...
alpha = ...
phi = TestFunction(V)
phiSUPG = phi + alpha*dot(u, grad(phi))
a = phiSUPG*....*dx + ....
On 31 August 2015 at 22:31, Justin Chang <jychang48(a)gmail.com> wrote:
> Hi all,
>
> 1) How do I define a volumetric source/sink for a specific region? For
> example, if I am solving a diffusion equation with homogeneous BCs and have
> a volumetric source f(x,y) such that
>
> if (x[0] >= 0.25 && x[0] <= 0.30 && x[1] >= 0.25 && x[1] <= 0.30):
> f == 1
> else:
> f == 0
>
> 2) I want to play around with PETSc's Variational Inequality
> <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/SNES/SNESVINEWTONSS…>
> solver. In my own PETSc code all I had to do was provide a Jacobian matrix,
> a residual vector, pass in the options -snes_type vinewtonssls (or it might
> be vinewotnrsls), and pass in SNESVINEWTONSSLS
> <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/SNES/SNESVINEWTONSS…> .
> From my last discussion with some PETSc folks it seems this VI feature is
> relatively now, so I was wondering if that last item can be invoked through
> the version of petsc4py you guys are using?
>
> 3) Is it possible to employ the SUPG formulation for advection-diffusion
> problems in firedrake?
>
> Thanks,
> Justin
>
Hi All,
I have a branch going through testing now which will handle
PETSC_CONFIGURE_OPTIONS in a smarter way. Basically we'll just make sure
that the things which are required are there, and we'll also honour
anything the user has already set. I've also updated docs to say that
that's what happens.
I am also open to the suggestion that we could add more configuration
options to the default set. Would anyone like to suggest what they should
be?
Cheers,
David
On Fri, 28 Aug 2015 at 10:25 Mitchell, Lawrence <
lawrence.mitchell(a)imperial.ac.uk> wrote:
>
> > On 27 Aug 2015, at 23:57, Justin Chang <jychang48(a)gmail.com> wrote:
> >
> > David,
> >
> > A few comments:
> >
> > 1) If I run that script (w/ --developer) and have the latest installed
> > petsc-dev and established $PETSC_DIR, the installation of petsc4py
> > will produce an error similar to the one reported in this thread:
> > http://lists.mcs.anl.gov/pipermail/petsc-users/2015-July/026404.html
> > to which I had to either: a) use the patch Satish provided or b) clone
> > from https://bitbucket.org/petsc/petsc4py and not mapdes/petsc4py. If
> > I unset $PETSC_DIR then the installation works fine. Not a major
> > problem or anything but it seems to me that the installation script is
> > using a somewhat outdated petsc4py so I was wondering if this was
> > intentional or not.
>
>
> The installation uses a matching set of PETSc/petsc4py, that are
> "up-to-date" and known good. But unfortunately it appears that if
> PETSC_DIR is set to a more recent PETSc, then the downloaded petsc4py might
> not match. The script provides a warning in this case, but it was maybe
> missed in all the other output. Perhaps if PETSC_DIR is set the script
> should exit with an error (overridable if you definitely want to use the
> external PETSc).
>
> > 2) Some of my problems require some PETSc external packages like
> > hypre, but I see that the default configure options only include
> > chaco, ctegen, and triangle. I had to manually add the
> > --download-hypre to the firedrake-install to get what I wanted, but do
> > you think y'all can also include other important packages like mumps,
> > superlu_dist, hdf5, metis, parmetis, etc?
>
> You can add options to the petsc configure call by setting the environment
> variable PETSC_CONFIGURE_OPTIONS to appropriate values. Make sure you
> include the necessary packages (chaco, ctetgen and triangle). Perhaps
> instead one should be able specify additional PETSc configure options with
> something like:
>
> firedrake-install ... --petsc-options="..."
>
>
> Cheers,
>
> Lawrence
>
Dear Firedrakers,
I am happy to be able to announce the release of our new install and
upgrade system for Firedrake. The new system makes installing and upgrading
Firedrake much easier than the previous instruction set. Particular
features:
* single command install on MacOS X, Ubuntu (and sometimes other Linux).
* support for per user or global installations.
* support for installation into a Python virtualenv to ensure no
interference between Firedrake and other software (in particular this makes
installing Firedrake alongside Dolfin safer and easier).
* Support for "developer mode" installs where core Firedrake components run
from the source directories.
The core Firedrake developers will be using this install, and it's the
install system used by our automated build testing, so you can be assured
that it will always be well tested and emerging bugs will be fixed fast.
Installation instructions are available now from
http://www.firedrakeproject.org/download.html
I look forward to your feedback on the new system.
Cheers,
David
Hi Justin,
Thanks for that, you uncovered a bug in some of how our internal imports
work. I have just merged a fix into Firedrake master so hopefully this
should now work.
Cheers,
David
On Wed, 26 Aug 2015 at 10:29 Justin Chang <jychang48(a)gmail.com> wrote:
> 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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)imperial.ac.uk
> >> >> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >> >> >
> >> >> >> >
> >> >> >> > _______________________________________________
> >> >> >> > firedrake mailing list
> >> >> >> > firedrake(a)imperial.ac.uk
> >> >> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >> >> >
> >> >> >>
> >> >> >> _______________________________________________
> >> >> >> firedrake mailing list
> >> >> >> firedrake(a)imperial.ac.uk
> >> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >> >
> >> >> >
> >> >> >
> >> >> > _______________________________________________
> >> >> > firedrake mailing list
> >> >> > firedrake(a)imperial.ac.uk
> >> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >> >
> >> >>
> >> >> _______________________________________________
> >> >> firedrake mailing list
> >> >> firedrake(a)imperial.ac.uk
> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >
> >> >
> >> > _______________________________________________
> >> > firedrake mailing list
> >> > firedrake(a)imperial.ac.uk
> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >
> >>
> >> _______________________________________________
> >> firedrake mailing list
> >> firedrake(a)imperial.ac.uk
> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
> >
> > _______________________________________________
> > firedrake mailing list
> > firedrake(a)imperial.ac.uk
> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
>
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(a)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(a)imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)imperial.ac.uk
> >> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >> >
> >> >> >
> >> >> > _______________________________________________
> >> >> > firedrake mailing list
> >> >> > firedrake(a)imperial.ac.uk
> >> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >> >
> >> >>
> >> >> _______________________________________________
> >> >> firedrake mailing list
> >> >> firedrake(a)imperial.ac.uk
> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >
> >> >
> >> >
> >> > _______________________________________________
> >> > firedrake mailing list
> >> > firedrake(a)imperial.ac.uk
> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >
> >>
> >> _______________________________________________
> >> firedrake mailing list
> >> firedrake(a)imperial.ac.uk
> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
> >
> > _______________________________________________
> > firedrake mailing list
> > firedrake(a)imperial.ac.uk
> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
>
> _______________________________________________
> firedrake mailing list
> firedrake(a)imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
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(a)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(a)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(a)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(a)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(a)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(a)imperial.ac.uk
> >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >
> >> >
> >> > _______________________________________________
> >> > firedrake mailing list
> >> > firedrake(a)imperial.ac.uk
> >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >> >
> >>
> >> _______________________________________________
> >> firedrake mailing list
> >> firedrake(a)imperial.ac.uk
> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
> >
> >
> > _______________________________________________
> > firedrake mailing list
> > firedrake(a)imperial.ac.uk
> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
>
> _______________________________________________
> firedrake mailing list
> firedrake(a)imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
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(a)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(a)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(a)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(a)imperial.ac.uk
> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
> >
> > _______________________________________________
> > firedrake mailing list
> > firedrake(a)imperial.ac.uk
> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
> >
>
> _______________________________________________
> firedrake mailing list
> firedrake(a)imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
Don't worry, found it. It was a bug in FIAT, that doesn't behave right if
the DOF is the integral of a scalar function.
all the best
--cjc
On 21 August 2015 at 17:52, Cotter, Colin J <colin.cotter(a)imperial.ac.uk>
wrote:
> Dear all,
> I realised that the Taylor basis used in FEM actually doesn't use the
> function evaluation for the lowest order part, it uses the element mean of
> the function. I
> started a new branch of FIAT to try to fix this, taylor-dg, but something
> I did is wrong.
>
> When running the test (just execute FIAT/discontinuous_taylor.py), I get
>
> Traceback (most recent call last):
> File "discontinuous_taylor.py", line 70, in <module>
> element = DiscontinuousTaylor(T, 1)
> File "discontinuous_taylor.py", line 65, in DiscontinuousTaylor
> return HigherOrderDiscontinuousTaylor( ref_el, degree )
> File "discontinuous_taylor.py", line 59, in __init__
> finite_element.FiniteElement.__init__( self, poly_set, dual, degree,
> formdegree )
> File "/home/cjc1/firedrake/fiat/FIAT/finite_element.py", line 47, in
> __init__
> dualmat = dual.to_riesz( poly_set )
> File "/home/cjc1/firedrake/fiat/FIAT/dual_set.py", line 64, in to_riesz
> self.mat[i][:] = self.nodes[i].to_riesz( poly_set )
> File "/home/cjc1/firedrake/fiat/FIAT/functional.py", line 330, in
> to_riesz
> result[self.comp, :] = numpy.dot(bfs, wts)
> IndexError: too many indices for array
>
> Can anyone see what is wrong?
>
> all the best
> --cjc
>
> --
> http://www.imperial.ac.uk/people/colin.cotter
>
> www.cambridge.org/9781107663916
>
>
>
--
http://www.imperial.ac.uk/people/colin.cotterwww.cambridge.org/9781107663916
Dear all,
I realised that the Taylor basis used in FEM actually doesn't use the
function evaluation for the lowest order part, it uses the element mean of
the function. I
started a new branch of FIAT to try to fix this, taylor-dg, but something I
did is wrong.
When running the test (just execute FIAT/discontinuous_taylor.py), I get
Traceback (most recent call last):
File "discontinuous_taylor.py", line 70, in <module>
element = DiscontinuousTaylor(T, 1)
File "discontinuous_taylor.py", line 65, in DiscontinuousTaylor
return HigherOrderDiscontinuousTaylor( ref_el, degree )
File "discontinuous_taylor.py", line 59, in __init__
finite_element.FiniteElement.__init__( self, poly_set, dual, degree,
formdegree )
File "/home/cjc1/firedrake/fiat/FIAT/finite_element.py", line 47, in
__init__
dualmat = dual.to_riesz( poly_set )
File "/home/cjc1/firedrake/fiat/FIAT/dual_set.py", line 64, in to_riesz
self.mat[i][:] = self.nodes[i].to_riesz( poly_set )
File "/home/cjc1/firedrake/fiat/FIAT/functional.py", line 330, in to_riesz
result[self.comp, :] = numpy.dot(bfs, wts)
IndexError: too many indices for array
Can anyone see what is wrong?
all the best
--cjc
--
http://www.imperial.ac.uk/people/colin.cotterwww.cambridge.org/9781107663916