August 2015 - firedrake - imperial.ac.uk

Re: [firedrake] Specify boundary conditions at a corner
by David Ham 24 Aug '15

24 Aug '15

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 >

2 1

Re: [firedrake] Specify boundary conditions at a corner
by Andrew McRae 23 Aug '15

23 Aug '15

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 >

2 1

Re: [firedrake] Taylor again
by Colin Cotter 21 Aug '15

21 Aug '15

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.cotter www.cambridge.org/9781107663916

1 0

Taylor again
by Colin Cotter 21 Aug '15

21 Aug '15

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

1 0

Specify boundary conditions at a corner
by Justin Chang 21 Aug '15

21 Aug '15

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

1 0

Re: [firedrake] Solve a Variational problem in a part of the domain
by Colin Cotter 17 Aug '15

17 Aug '15

Hi Anna, Not as written, no. I think that we would need to design a bespoke solver for that. all the best --cjc On 17 August 2015 at 12:18, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Hi all, > > I have to solve a problem with the bilinear form similar to the following: > > a = ( v*u + dt*inner(grad(u),grad(v)) )*dx + ( v*dt*assemble(u*dx) )*dx > > The last term is essentially the product of the integral of u and the > integral of test function v. > Is it even possible to solve this in Firedrake? > > Best, > > Anna. > > > On 14/08/15 15:43, Lawrence Mitchell wrote: > > On 14 Aug 2015, at 11:55, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> <a.kalogirou(a)leeds.ac.uk> wrote: > > Hi, > > I was planning to define the heavyside function as a DG0 function. > > Any ideas about solving equation (1b), which contains a time-update for both a function and a scalar? > > I haven't been fully following along, but it looks like the formulation in (1f) introduces, effectively, a global coupling. Is that right? Assuming you didn't have any implementation constraints, how would you go about solving this? > > Lawrence > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > -- > > Dr Anna Kalogirou > Research Fellow > School of Mathematics > University of Leeds > > http://www1.maths.leeds.ac.uk/~matak/ > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] Solve a Variational problem in a part of the domain
by Colin Cotter 17 Aug '15

17 Aug '15

Hi Anna, It *was* interesting! How do you intend to compute the integrals against the Heaviside function? The Heaviside function is non-polynomial so can't be computed exactly by quadrature. You are right that our standard setup can't cope with this situation currently. What we really need are function spaces that are defined to be constant in an entire column. We should try to sketch out with David and Lawrence what might be required in the infrastructure to support this. I can see it being useful for many other free surface type problems as well. all the best --Colin On 13 August 2015 at 15:36, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Hi Colin, > > I attached a two-page document, where the system of equations I was > talking about in the previous email is (1a)-(1e). I am puzzled with > equation (1b) really. I can also eliminate (1b) & (1d) and solve for the > remaining (1a), (1c) & (1f), but in this case I still don't know how to > deal with (1f) since there is an integral into the integral, which contains > the unknown function lambda^{n+1/2}. > > Thanks, Anna. > > > On 13/08/15 13:47, Colin Cotter wrote: > > Hi Anna, > Sounds interesting. Please could you provide a bit more detail? > > all the best > --Colin > > On 13 August 2015 at 12:09, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> > wrote: > >> Dear all, >> >> I have a system of equations to solve, which involves three space >> dependent functions phi, eta, lambda and two constants/scalars Z, W. >> These need to be solved simultaneously because all the equations involve >> at least 2 unknowns. >> >> How do I solve that, considering that one of the scalar equations >> includes a spacial integral of one of the (still unknown) functions? Is >> it best to define the scalars as Constants? >> >> This problem goes away when I write down the system in a standard FEM >> formulation, introducing the mass matrix etc. In this case, it is clear >> that I could solve for that function and then consecutively solve each >> of the remaining 4 equations. It is not that obvious how I can do that >> using Firedrake, that is why I thought I would have to solve >> simultaneously, but then I have the problem described above. >> >> Regards, >> >> Anna. >> >> >> On 06/08/15 11:04, Lawrence Mitchell wrote: >> > -----BEGIN PGP SIGNED MESSAGE----- >> > Hash: SHA1 >> > >> > Hi Anna, >> > >> > On 06/08/15 10:37, Anna Kalogirou wrote: >> >> Dear all, >> >> >> >> I have a rather simple question but I would like to get some >> >> feedback from someone in the Firedrake team. >> >> >> >> I am now working on the problem which includes solving a >> >> variational problem in a part of the domain only (a water wave >> >> problem which includes a floating body). >> >> >> >> I can think of a couple of possible solutions on how to solve this: >> >> 1. Define two domains and solve the problem separately in each >> >> domain. However, I will have to deal with nonzero boundary >> >> conditions on the common boundary. >> >> >> >> 2. I prefer solving the problem in the whole domain, since most of >> >> the equations/functions are valid everywhere. Then I can define a >> >> Heavyside step function which will be 0 in one part and 1 in the >> >> part of the domain I am interested in (under the floating body). I >> >> will essentially write down a variational problem valid everywhere, >> >> but will actually be zero in a part of the domain. >> >> >> >> Is the 2nd step a good approach? The question essentially is how >> >> to split a mass matrix M_kl which is defined everywhere, and solve >> >> and integral form in a part of the domain only. >> > I think step two is a fine approach. However, note the following >> > issues at present. The way you would have to do this currently is as >> > follows: >> > >> > Define a DG0 field to hold your indicator function >> > >> > indicator = Function(DG0) >> > >> > # Set it to 1 in the appropriate part of the domain >> > indicator.interpolate(...) >> > >> > # Now use this extra field everywhere when defining your variational >> > # problem. >> > >> > However, all your integrals are still over the whole domain, you just >> > pick up lots of zeros. >> > >> > Steadily climbing our todo list (and at an increasing pace, I feel), >> > is the ability to define proper sub domains in a mesh, and then be >> > able to perform integrals over them. Until that time, I think >> > approach 2 is somewhat easier to do than approach 1. >> > >> > Cheers, >> > >> > LAwrence >> > -----BEGIN PGP SIGNATURE----- >> > Version: GnuPG v1 >> > >> > iQEcBAEBAgAGBQJVwzElAAoJECOc1kQ8PEYvGswIAN+Mr6hD+U0NmsnZAGminaq+ >> > j/SqmI+o7GlbXRSgPqEfrpyWwFiiFVN6DnPmTboRooZr3h93ZhGDhECQ/O12XKA/ >> > 00yNuCiDn4YBaXJamWBM6u8ILdZYGMyVNTG1JOaZZIMBoiIkrTyYpXatGApcE4sD >> > 4T+DoAYSUEY7u6MV4eJmRc1lzZEd8Rw43YP+Viy0itCs2jKV+HsDpvpvu1DU/FSQ >> > zMeE/cRR8It6aImk1L3OON8Zo2ZB0HNirK6YuD3sMkIoqLEYANnNl8qaeU3LL227 >> > 3TL4ddkALXQLGaYGo5ETgtMKI6afNDIsfjwGUcVRhTTxHM3GKOMQgWV6miNtgv4= >> > =wE25 >> > -----END PGP SIGNATURE----- >> > >> > _______________________________________________ >> > firedrake mailing list >> > firedrake(a)imperial.ac.uk >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> -- >> >> Dr Anna Kalogirou >> Research Fellow >> School of Mathematics >> University of Leeds >> >> http://www1.maths.leeds.ac.uk/~matak/ >> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > -- > > Dr Anna Kalogirou > Research Fellow > School of Mathematics > University of Leeds > > http://www1.maths.leeds.ac.uk/~matak/ > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 3

Re: [firedrake] Solve a Variational problem in a part of the domain
by Colin Cotter 14 Aug '15

14 Aug '15

Hi Anna, Ah, I see. So you "snap" the Heaviside function to the nearest gridpoint. all the best --cjc On 14 August 2015 at 11:55, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Hi, > > I was planning to define the heavyside function as a DG0 function. > > Any ideas about solving equation (1b), which contains a time-update for > both a function and a scalar? > > Thanks, > Anna. > > > > On 13/08/15 18:35, Colin Cotter wrote: > > Hi Anna, > It *was* interesting! > > How do you intend to compute the integrals against the Heaviside function? > The Heaviside function is non-polynomial so can't be computed exactly by > quadrature. > > You are right that our standard setup can't cope with this situation > currently. What we really need are function spaces that are defined to be > constant in an entire column. We should try to sketch out with David and > Lawrence what might be required in the infrastructure to support this. I > can see it being useful for many other free surface type problems as well. > > all the best > --Colin > > On 13 August 2015 at 15:36, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> > wrote: > >> Hi Colin, >> >> I attached a two-page document, where the system of equations I was >> talking about in the previous email is (1a)-(1e). I am puzzled with >> equation (1b) really. I can also eliminate (1b) & (1d) and solve for the >> remaining (1a), (1c) & (1f), but in this case I still don't know how to >> deal with (1f) since there is an integral into the integral, which contains >> the unknown function lambda^{n+1/2}. >> >> Thanks, Anna. >> >> >> On 13/08/15 13:47, Colin Cotter wrote: >> >> Hi Anna, >> Sounds interesting. Please could you provide a bit more detail? >> >> all the best >> --Colin >> >> On 13 August 2015 at 12:09, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> >> wrote: >> >>> Dear all, >>> >>> I have a system of equations to solve, which involves three space >>> dependent functions phi, eta, lambda and two constants/scalars Z, W. >>> These need to be solved simultaneously because all the equations involve >>> at least 2 unknowns. >>> >>> How do I solve that, considering that one of the scalar equations >>> includes a spacial integral of one of the (still unknown) functions? Is >>> it best to define the scalars as Constants? >>> >>> This problem goes away when I write down the system in a standard FEM >>> formulation, introducing the mass matrix etc. In this case, it is clear >>> that I could solve for that function and then consecutively solve each >>> of the remaining 4 equations. It is not that obvious how I can do that >>> using Firedrake, that is why I thought I would have to solve >>> simultaneously, but then I have the problem described above. >>> >>> Regards, >>> >>> Anna. >>> >>> >>> On 06/08/15 11:04, Lawrence Mitchell wrote: >>> > -----BEGIN PGP SIGNED MESSAGE----- >>> > Hash: SHA1 >>> > >>> > Hi Anna, >>> > >>> > On 06/08/15 10:37, Anna Kalogirou wrote: >>> >> Dear all, >>> >> >>> >> I have a rather simple question but I would like to get some >>> >> feedback from someone in the Firedrake team. >>> >> >>> >> I am now working on the problem which includes solving a >>> >> variational problem in a part of the domain only (a water wave >>> >> problem which includes a floating body). >>> >> >>> >> I can think of a couple of possible solutions on how to solve this: >>> >> 1. Define two domains and solve the problem separately in each >>> >> domain. However, I will have to deal with nonzero boundary >>> >> conditions on the common boundary. >>> >> >>> >> 2. I prefer solving the problem in the whole domain, since most of >>> >> the equations/functions are valid everywhere. Then I can define a >>> >> Heavyside step function which will be 0 in one part and 1 in the >>> >> part of the domain I am interested in (under the floating body). I >>> >> will essentially write down a variational problem valid everywhere, >>> >> but will actually be zero in a part of the domain. >>> >> >>> >> Is the 2nd step a good approach? The question essentially is how >>> >> to split a mass matrix M_kl which is defined everywhere, and solve >>> >> and integral form in a part of the domain only. >>> > I think step two is a fine approach. However, note the following >>> > issues at present. The way you would have to do this currently is as >>> > follows: >>> > >>> > Define a DG0 field to hold your indicator function >>> > >>> > indicator = Function(DG0) >>> > >>> > # Set it to 1 in the appropriate part of the domain >>> > indicator.interpolate(...) >>> > >>> > # Now use this extra field everywhere when defining your variational >>> > # problem. >>> > >>> > However, all your integrals are still over the whole domain, you just >>> > pick up lots of zeros. >>> > >>> > Steadily climbing our todo list (and at an increasing pace, I feel), >>> > is the ability to define proper sub domains in a mesh, and then be >>> > able to perform integrals over them. Until that time, I think >>> > approach 2 is somewhat easier to do than approach 1. >>> > >>> > Cheers, >>> > >>> > LAwrence >>> > -----BEGIN PGP SIGNATURE----- >>> > Version: GnuPG v1 >>> > >>> > iQEcBAEBAgAGBQJVwzElAAoJECOc1kQ8PEYvGswIAN+Mr6hD+U0NmsnZAGminaq+ >>> > j/SqmI+o7GlbXRSgPqEfrpyWwFiiFVN6DnPmTboRooZr3h93ZhGDhECQ/O12XKA/ >>> > 00yNuCiDn4YBaXJamWBM6u8ILdZYGMyVNTG1JOaZZIMBoiIkrTyYpXatGApcE4sD >>> > 4T+DoAYSUEY7u6MV4eJmRc1lzZEd8Rw43YP+Viy0itCs2jKV+HsDpvpvu1DU/FSQ >>> > zMeE/cRR8It6aImk1L3OON8Zo2ZB0HNirK6YuD3sMkIoqLEYANnNl8qaeU3LL227 >>> > 3TL4ddkALXQLGaYGo5ETgtMKI6afNDIsfjwGUcVRhTTxHM3GKOMQgWV6miNtgv4= >>> > =wE25 >>> > -----END PGP SIGNATURE----- >>> > >>> > _______________________________________________ >>> > firedrake mailing list >>> > firedrake(a)imperial.ac.uk >>> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> -- >>> >>> Dr Anna Kalogirou >>> Research Fellow >>> School of Mathematics >>> University of Leeds >>> >>> http://www1.maths.leeds.ac.uk/~matak/ >>> >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >> >> >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> -- >> >> Dr Anna Kalogirou >> Research Fellow >> School of Mathematics >> University of Leeds >> >> http://www1.maths.leeds.ac.uk/~matak/ >> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] Solve a Variational problem in a part of the domain
by Colin Cotter 13 Aug '15

13 Aug '15

Hi Anna, Sounds interesting. Please could you provide a bit more detail? all the best --Colin On 13 August 2015 at 12:09, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Dear all, > > I have a system of equations to solve, which involves three space > dependent functions phi, eta, lambda and two constants/scalars Z, W. > These need to be solved simultaneously because all the equations involve > at least 2 unknowns. > > How do I solve that, considering that one of the scalar equations > includes a spacial integral of one of the (still unknown) functions? Is > it best to define the scalars as Constants? > > This problem goes away when I write down the system in a standard FEM > formulation, introducing the mass matrix etc. In this case, it is clear > that I could solve for that function and then consecutively solve each > of the remaining 4 equations. It is not that obvious how I can do that > using Firedrake, that is why I thought I would have to solve > simultaneously, but then I have the problem described above. > > Regards, > > Anna. > > > On 06/08/15 11:04, Lawrence Mitchell wrote: > > -----BEGIN PGP SIGNED MESSAGE----- > > Hash: SHA1 > > > > Hi Anna, > > > > On 06/08/15 10:37, Anna Kalogirou wrote: > >> Dear all, > >> > >> I have a rather simple question but I would like to get some > >> feedback from someone in the Firedrake team. > >> > >> I am now working on the problem which includes solving a > >> variational problem in a part of the domain only (a water wave > >> problem which includes a floating body). > >> > >> I can think of a couple of possible solutions on how to solve this: > >> 1. Define two domains and solve the problem separately in each > >> domain. However, I will have to deal with nonzero boundary > >> conditions on the common boundary. > >> > >> 2. I prefer solving the problem in the whole domain, since most of > >> the equations/functions are valid everywhere. Then I can define a > >> Heavyside step function which will be 0 in one part and 1 in the > >> part of the domain I am interested in (under the floating body). I > >> will essentially write down a variational problem valid everywhere, > >> but will actually be zero in a part of the domain. > >> > >> Is the 2nd step a good approach? The question essentially is how > >> to split a mass matrix M_kl which is defined everywhere, and solve > >> and integral form in a part of the domain only. > > I think step two is a fine approach. However, note the following > > issues at present. The way you would have to do this currently is as > > follows: > > > > Define a DG0 field to hold your indicator function > > > > indicator = Function(DG0) > > > > # Set it to 1 in the appropriate part of the domain > > indicator.interpolate(...) > > > > # Now use this extra field everywhere when defining your variational > > # problem. > > > > However, all your integrals are still over the whole domain, you just > > pick up lots of zeros. > > > > Steadily climbing our todo list (and at an increasing pace, I feel), > > is the ability to define proper sub domains in a mesh, and then be > > able to perform integrals over them. Until that time, I think > > approach 2 is somewhat easier to do than approach 1. > > > > Cheers, > > > > LAwrence > > -----BEGIN PGP SIGNATURE----- > > Version: GnuPG v1 > > > > iQEcBAEBAgAGBQJVwzElAAoJECOc1kQ8PEYvGswIAN+Mr6hD+U0NmsnZAGminaq+ > > j/SqmI+o7GlbXRSgPqEfrpyWwFiiFVN6DnPmTboRooZr3h93ZhGDhECQ/O12XKA/ > > 00yNuCiDn4YBaXJamWBM6u8ILdZYGMyVNTG1JOaZZIMBoiIkrTyYpXatGApcE4sD > > 4T+DoAYSUEY7u6MV4eJmRc1lzZEd8Rw43YP+Viy0itCs2jKV+HsDpvpvu1DU/FSQ > > zMeE/cRR8It6aImk1L3OON8Zo2ZB0HNirK6YuD3sMkIoqLEYANnNl8qaeU3LL227 > > 3TL4ddkALXQLGaYGo5ETgtMKI6afNDIsfjwGUcVRhTTxHM3GKOMQgWV6miNtgv4= > > =wE25 > > -----END PGP SIGNATURE----- > > > > _______________________________________________ > > firedrake mailing list > > firedrake(a)imperial.ac.uk > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > -- > > Dr Anna Kalogirou > Research Fellow > School of Mathematics > University of Leeds > > http://www1.maths.leeds.ac.uk/~matak/ > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Solve a Variational problem in a part of the domain
by Anna Kalogirou 13 Aug '15

13 Aug '15

Dear all, I have a rather simple question but I would like to get some feedback from someone in the Firedrake team. I am now working on the problem which includes solving a variational problem in a part of the domain only (a water wave problem which includes a floating body). I can think of a couple of possible solutions on how to solve this: 1. Define two domains and solve the problem separately in each domain. However, I will have to deal with nonzero boundary conditions on the common boundary. 2. I prefer solving the problem in the whole domain, since most of the equations/functions are valid everywhere. Then I can define a Heavyside step function which will be 0 in one part and 1 in the part of the domain I am interested in (under the floating body). I will essentially write down a variational problem valid everywhere, but will actually be zero in a part of the domain. Is the 2nd step a good approach? The question essentially is how to split a mass matrix M_kl which is defined everywhere, and solve and integral form in a part of the domain only. Thank you. Regards, Anna.

3 2