- firedrake - imperial.ac.uk

Re: [firedrake] Firedrake '17 Workshop
by David Ham 28 Feb '17

28 Feb '17

Hi Francis, I plan to get a schedule up for the workshop early next week (when I get back from SIAM CSE). I'm also in communication with a venue for the dinner on Monday. I think that informally getting together on the Sunday night is a great idea. We'll probably just designate a pub near Imperial and anyone who is in town already can come along. Regards, David On Sun, 26 Feb 2017 at 22:13 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello everyone, > > I wanted to ask a couple of questions about the workshop. > > Does anyone know when the schedule might be posted? > > Are there any plans to get together for dinner/drinks the Sunday before? > I will be arriving from across the pond that morning and not sure what > state I'll be in but would be happy to get together with people beforehand. > > Cheers, Francis > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 <(519)%20888-4567> > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Re: [firedrake] defining an expression piecewise
by David Ham 28 Feb '17

28 Feb '17

Hi Francis, There are two things here. The first is the proximate cause of your exception. If f is a FunctionSpace then it doesn't have an interpolate method. Instead, you would call the interpolate function: c = interpolate(MyExpression(), f). However, creating a subclass of Expression is not the best way to solve this problem. Instead you should interpolate a UFL expression: x = SpatialCoordinate(f.mesh()) c = interpolate(conditional(le(x[0], 0.5), 1.0, 0.0), f) Regards, David On Sun, 26 Feb 2017 at 22:06 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I am trying to define an expression piecewise and see that I need to do > that by defining a class. Below is something that I tried, building on > what I read, but it has an error, see below. > > class MyExpression(Expression): > def eval(self, value, x): > if x[0] <= 0.5: > value[0] = 1.0 > else: > value[0] = 0.0 > > def value_shape(self): > return (1,) > > f.interpolate(MyExpression()) > > It fails with the following error: > > AttributeError: 'FunctionSpace' object has no attribute 'interpolate' > > > Can someone maybe point me in the right direction on how to do tis > properly? > > Cheers, Francis > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 <(519)%20888-4567> > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Re: [firedrake] defining an expression piecewise
by Colin Cotter 27 Feb '17

27 Feb '17

Hi Francis, What is the function space you are using? You can build a piecewise constant function using the C-string Expression via the ?: ternary operator in C. all the best --cjc On 26 February 2017 at 22:05, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I am trying to define an expression piecewise and see that I need to do > that by defining a class. Below is something that I tried, building on > what I read, but it has an error, see below. > > class MyExpression(Expression): > def eval(self, value, x): > if x[0] <= 0.5: > value[0] = 1.0 > else: > value[0] = 0.0 > > def value_shape(self): > return (1,) > > f.interpolate(MyExpression()) > > It fails with the following error: > > AttributeError: 'FunctionSpace' object has no attribute 'interpolate' > > > Can someone maybe point me in the right direction on how to do tis > properly? > > Cheers, Francis > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 <(519)%20888-4567> > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Firedrake '17 Workshop
by Francis Poulin 27 Feb '17

27 Feb '17

Hello everyone, I wanted to ask a couple of questions about the workshop. Does anyone know when the schedule might be posted? Are there any plans to get together for dinner/drinks the Sunday before? I will be arriving from across the pond that morning and not sure what state I'll be in but would be happy to get together with people beforehand. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

1 0

defining an expression piecewise
by Francis Poulin 27 Feb '17

27 Feb '17

Hello, I am trying to define an expression piecewise and see that I need to do that by defining a class. Below is something that I tried, building on what I read, but it has an error, see below. class MyExpression(Expression): def eval(self, value, x): if x[0] <= 0.5: value[0] = 1.0 else: value[0] = 0.0 def value_shape(self): return (1,) f.interpolate(MyExpression()) It fails with the following error: AttributeError: 'FunctionSpace' object has no attribute 'interpolate' Can someone maybe point me in the right direction on how to do tis properly? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

1 0

Re: [firedrake] RectangleMesh
by David Ham 24 Feb '17

24 Feb '17

Hi Floriane, you can do this using a par_loop with a custom kernel. I don't have time to work out the precise syntax for this today so I've raised an issue on GitHub so we don't lose this question. See: https://github.com/firedrakeproject/firedrake/issues/1022 I hope to get back to this next week, or maybe someone else will get to it faster. Regards, David On Thu, 23 Feb 2017 at 15:03 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> wrote: > > Dear David, > > > Thanks for your fast answer. > > > I change the coordinates for saving purpose only (once the WF are solved). > > So my objective is to solve the WF in a fixed 2D RectangleMesh, and then > save the functions in a 3D mesh in which I apply the coordinate transforms > to get moving boundaries. > > I managed to do this when Ny>1, and I changed my code so that it > works for Ny=1, but I would like a general command that works whatever the > value of Ny. > > > So my question would be : how can I make the link between the > x-coordinates in a RectangleMesh, and the x-coordinates in the 3D mesh > obtained by extrusion of the RectangleMesh ? > > > For instance, for a function h_2D defined in the RectangleMesh, if I want > to change the z-coordinate of the extruded mesh as z_new = z*h_2D/H0 so > that the top boundary of the 3D domain for a given (xi,yi) coordinates > moves as h_2D(xi,yi), how can I find the indices "..." to substitute into > the following expression : > > > mesh_3D.coordinates.dat.data[ ..., 2]*= h_2D.dat.data[ ... ]/H0 ? > > > Thank you very much, > > Floriane > > > > > > > > > > > ------------------------------ > *De :* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > de la part de David Ham <David.Ham(a)imperial.ac.uk> > *Envoyé :* jeudi 23 février 2017 14:16:06 > *À :* firedrake > *Objet :* Re: [firedrake] RectangleMesh > > Dear Floriane, > > The thing to understand is that there is no guaranteed relationship > between vertex number and the location of that vertex. Vertex numbers are > (as far as the user is concerned), completely arbitrary. This is fine > because none of the correct ways to access vertices depend on what the > order is. > > When you say that you are applying a change of coordinates, what do you > mean? Do you mean that you are changing the coordinate values before you > compute? If so, please tell us mathematically what you are attempting to do > and we can tell you the right way to do it. > > Regards, > > David > > On Thu, 23 Feb 2017 at 14:08 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> > wrote: > > Dear all, > > > I have a question concerning the mesh definition with RectangleMesh: > > > I define my mesh as a Rectangle mesh, with quadrilaterals as: > > hor_mesh = RectangleMesh(Nx, Ny, 10.0, 1.0, quadrilateral=True). > > > As my functions are only (x,t) dependent, I'd like to get only one element > in y, i.e. Ny=1. > > However I noticed that if I set Ny=1, then the x-coordinates are swapped > from 10.0 to 0.0 . Basically, I obtain > > hor_mesh.coordinates.dat.data[:,0] = [10.0 10.0 9.9 9.9 ...... 0.0] > > instead of > > hor_mesh.coordinates.dat.data[:,0] = [0.0 0.0 0.1 0.1 ........ 10.0] (that > I obtain for Ny>1) > > > Can someone explain me what happens there ? Is there a way to avoid that ? > > And if I use the integral on ds(1), will it be at x=0 of x=10 in that case? > > > The reason why I want to understand this is because I apply a change of > coordinate in some part of the domain (let's say for x<L<Lx) and the > coordinate transform does not give the right result for Ny=1 (while it does > for Ny>1). > > > Thanks in advance for your help, > > Floriane > > > > -- > Dr David Ham > Department of Mathematics > Imperial College London > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Re: [firedrake] function at a point
by David Ham 23 Feb '17

23 Feb '17

Aha. I see the issue. I need to give this some thought (although right now I need to go and teach). Regards, David On Thu, 23 Feb 2017 at 15:28 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Attached, see weak form in equation (20). > > xp used in equation (20) is define in equation (19). > > Hope that settles the question? > > You should be familiar with some of this at it also emerged in your OMAE > paper! > > Regarding a time discretisation: any, but just do first order at the > moment! > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of David Ham <David.Ham(a)imperial.ac.uk> > *Sent:* Thursday, February 23, 2017 2:17 PM > *To:* firedrake > *Subject:* Re: [firedrake] function at a point > > Hi Onno, > > I'm afraid I don't think we understand what you are asking for. Can you > show us the maths please? > > Regards, > > David > > On Thu, 23 Feb 2017 at 08:07 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > > Dear Firedrake(rs), > > > I am considering a "standard" FEM problem with the following variables > h(chi,tau) and phi(chi,tau) as function of transformed coordinate chi and > time tau in chi=[0,L]. > > > In the weak formulation, the point value h(Lp,tau), so the function value at > 0<chi=Lp<L > > as well as, of course, the function h(chi,tau) appear in the space > integral of the weak form. > > In essence, due to a transformed moving boundary this point evaluation > arises. > > > I can formulate the detailed (nonlinear) matrix-vector FEM and its time > discretisation > > because after the expansion of h(chi,tau) the coefficients h_j(tau) appear > > and h(Lp,tau) = h_Np(tau) for j =Np, say, such that the distinction > between function > > and point values disappears. > > > I am struggling to see how I can do this on the weak form level in > firedrake. > > It would work in integral form by using a delta function after > introducing h(Lp,tau) > > as an auxiliary scalar variable but I assume that is not available? > > > What is the FD-tactic here which one can employ? > > > Thank you and best wishes, > > > Onno > > > > > -- > Dr David Ham > Department of Mathematics > Imperial College London > -- Dr David Ham Department of Mathematics Imperial College London

1 0

Re: [firedrake] function at a point
by David Ham 23 Feb '17

23 Feb '17

Hi Onno, I'm afraid I don't think we understand what you are asking for. Can you show us the maths please? Regards, David On Thu, 23 Feb 2017 at 08:07 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Dear Firedrake(rs), > > > I am considering a "standard" FEM problem with the following variables > h(chi,tau) and phi(chi,tau) as function of transformed coordinate chi and > time tau in chi=[0,L]. > > > In the weak formulation, the point value h(Lp,tau), so the function value at > 0<chi=Lp<L > > as well as, of course, the function h(chi,tau) appear in the space > integral of the weak form. > > In essence, due to a transformed moving boundary this point evaluation > arises. > > > I can formulate the detailed (nonlinear) matrix-vector FEM and its time > discretisation > > because after the expansion of h(chi,tau) the coefficients h_j(tau) appear > > and h(Lp,tau) = h_Np(tau) for j =Np, say, such that the distinction > between function > > and point values disappears. > > > I am struggling to see how I can do this on the weak form level in > firedrake. > > It would work in integral form by using a delta function after > introducing h(Lp,tau) > > as an auxiliary scalar variable but I assume that is not available? > > > What is the FD-tactic here which one can employ? > > > Thank you and best wishes, > > > Onno > > > > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Re: [firedrake] RectangleMesh
by David Ham 23 Feb '17

23 Feb '17

Dear Floriane, The thing to understand is that there is no guaranteed relationship between vertex number and the location of that vertex. Vertex numbers are (as far as the user is concerned), completely arbitrary. This is fine because none of the correct ways to access vertices depend on what the order is. When you say that you are applying a change of coordinates, what do you mean? Do you mean that you are changing the coordinate values before you compute? If so, please tell us mathematically what you are attempting to do and we can tell you the right way to do it. Regards, David On Thu, 23 Feb 2017 at 14:08 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> wrote: > Dear all, > > > I have a question concerning the mesh definition with RectangleMesh: > > > I define my mesh as a Rectangle mesh, with quadrilaterals as: > > hor_mesh = RectangleMesh(Nx, Ny, 10.0, 1.0, quadrilateral=True). > > > As my functions are only (x,t) dependent, I'd like to get only one element > in y, i.e. Ny=1. > > However I noticed that if I set Ny=1, then the x-coordinates are swapped > from 10.0 to 0.0 . Basically, I obtain > > hor_mesh.coordinates.dat.data[:,0] = [10.0 10.0 9.9 9.9 ...... 0.0] > > instead of > > hor_mesh.coordinates.dat.data[:,0] = [0.0 0.0 0.1 0.1 ........ 10.0] (that > I obtain for Ny>1) > > > Can someone explain me what happens there ? Is there a way to avoid that ? > > And if I use the integral on ds(1), will it be at x=0 of x=10 in that case? > > > The reason why I want to understand this is because I apply a change of > coordinate in some part of the domain (let's say for x<L<Lx) and the > coordinate transform does not give the right result for Ny=1 (while it does > for Ny>1). > > > Thanks in advance for your help, > > Floriane > > > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Firedrake .PVD visualisation
by Robert Long [RPG] 23 Feb '17

23 Feb '17

Hello Firedrake peoples, I have ran a simulation which outputs a time dependent series of results files, with a .pvd extension, and a load of .vtu files. When i try to visualise these in firedrake i just see the x-domain (its a 1d problem) with a function of intensities on it. DO you know how i can take the values of the function and plot them against x for an x-y plot? If not, are there alternative output options? Best, Rob

2 1