Dear Ed,
Here’s an attempt to lay out how boundary conditions work, which might enable you to do what you’re after.
Let’s make some basic data structures:
from firedrake
import *
m = UnitSquareMesh(2,2)
fs = FunctionSpace(m, "CG",
2)
bc = DirichletBC(fs, 1.0,
1)
Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function.
The (process-local) indices of the nodes in bc are given by:
print(bc.nodes)
array([ 0, 3, 4, 12, 14], dtype=int32)
So let’s suppose we want that indicator function:
f = Function(fs)
bc.apply(f)
If we actually look at the node values, we now find:
print(f.dat.data_ro)
array([1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0.])
Note that the fact that the nodal values are 1 is a function of the node definition. It will be true for CG and DG spaces, because the node definition for these cases is point evaluation. For more esoteric finite elements such as those
in H(div) and H(curl), the relationship between nodal values and Function value are more complex.
If you want the physical locations of these nodes then you’ll need to interpolate the mesh coordinates into the VectorFunctionSpace corresponding to the FunctionSpace you’re working in. This process is described here:
https://www.firedrakeproject.org/interpolation.html#interpolation-from-external-data
In order to only perform the operation over the boundary nodes, you’d simply replace the last line there with:
f.dat.data[bc.nodes]
= mydata(X.dat.data_ro[bc.nodes])
Regards,
David
From:
<firedrake-bounces@imperial.ac.uk> on behalf of Ed Bueler <elbueler@alaska.edu>
Reply-To: firedrake <firedrake@imperial.ac.uk>
Date: Saturday, 12 May 2018 at 01:36
To: firedrake <firedrake@imperial.ac.uk>
Subject: Re: [firedrake] mesh coordinates and boundaries questions
> 1) How does one get
the coordinates of the points on a boundary? That is, suppose I know the id
> number of that part, for example "id" as in "ds(id)" used in forming boundary
terms in the weak form.
> What is the corresponding idiom for getting the spatial coordinates?
Something like "mesh.coordinates(id)"
> (though that does not work)? In parallel I guess this function will only
be defined on a subset of the processors?
Generating a Function or PETSc Vec defined only on the desired part of the
boundary, and not defined elsewhere, would generally imply that some processes had no values, which can be awkward. In any case I don't need this any longer.
My much-reduced question has become: How to generate the indicator function
of part of the boundary? This can be locally-generated and is distributed the same way as a Function on the whole mesh. The FEniCS book shows an example of generating
MeshFunctions which indicate boundary parts, but perhaps this is just Function in Firedrake? What is the idiom for using the subdomain id to generate the indicator function?
> 2) Given a function along a boundary, how does one extend that to the whole mesh,
> by holding constant in one variable? In the case I care about the domain is of the special form
> D = {(x,y) | f(x) < y < g(x)}
> and one has values w(x) along the y=g(x) boundary. How does one build
a Function on the
> whole mesh which is
> W(x,y) = w(x)
> (I understand that my domain is special, namely logically rectangular.
And I can imagine an
> answer specific to an extruded mesh ... but I would like to not go there
now.) Again I am interested
> in a solution which works in parallel. Maybe I have to commit to solving
a first-order PDE just to do this?
I was missing something deeper here. In parallel each process only owns (knows about) its part of the mesh. Other than solve(), Firedrake
operations are almost all local in the sense that each process works on its part of the mesh, or across a halo, but in any case without global knowledge. By contrast, solve() is a well-defined interface for spreading certain effects from some locales (e.g.
the boundary) to all processes. One needs precise equations to solve ... but such equations give precise meaning to "spreading certain effects". Thus solve() is the key to global influence on the entries of a Function and/or its underlying mpi-type Vec.
My question could therefore answered this way: Set up a weak form to solve a differential equation problem describing the field you want
with the boundary values you want. Then call solve(). This approach can propagate information to all processes in a well-supported way.
In fact it turns out that what I wanted could be done by setting up and solving a scalar, linear, 2nd-order PDE problem. It works fine.
It is computationally-cheap compared to the "main" PDE I was solving, and scales better. In the longer term I almost certainly will use an extruded mesh, which may be even easier.
Ed
On Wed, May 9, 2018 at 10:54 PM, Ed Bueler <elbueler@alaska.edu>
wrote: