You can subscribe to this list here.
2010 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}
(1) 
_{Nov}
(1) 
_{Dec}


2011 
_{Jan}

_{Feb}

_{Mar}
(1) 
_{Apr}

_{May}
(1) 
_{Jun}
(9) 
_{Jul}
(2) 
_{Aug}
(6) 
_{Sep}
(11) 
_{Oct}
(15) 
_{Nov}
(4) 
_{Dec}
(9) 
2012 
_{Jan}
(7) 
_{Feb}
(14) 
_{Mar}
(11) 
_{Apr}

_{May}

_{Jun}
(3) 
_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}
(3) 
_{Dec}

2013 
_{Jan}
(8) 
_{Feb}

_{Mar}
(1) 
_{Apr}
(3) 
_{May}

_{Jun}

_{Jul}
(3) 
_{Aug}
(2) 
_{Sep}

_{Oct}

_{Nov}

_{Dec}
(4) 
2014 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}
(2) 
_{Jun}
(2) 
_{Jul}

_{Aug}

_{Sep}
(5) 
_{Oct}
(3) 
_{Nov}
(1) 
_{Dec}
(4) 
2015 
_{Jan}
(5) 
_{Feb}
(10) 
_{Mar}
(12) 
_{Apr}
(14) 
_{May}
(26) 
_{Jun}
(7) 
_{Jul}

_{Aug}

_{Sep}
(4) 
_{Oct}

_{Nov}

_{Dec}

2016 
_{Jan}
(1) 
_{Feb}
(9) 
_{Mar}
(3) 
_{Apr}
(3) 
_{May}
(12) 
_{Jun}
(2) 
_{Jul}
(16) 
_{Aug}
(4) 
_{Sep}
(18) 
_{Oct}
(3) 
_{Nov}
(4) 
_{Dec}
(10) 
2017 
_{Jan}

_{Feb}

_{Mar}
(9) 
_{Apr}
(14) 
_{May}
(12) 
_{Jun}
(17) 
_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

From: Rafael Vazquez <vazquez@im...>  20170621 10:29:52

Dear Tony, you are working in the interface space sp_interface = space.sp_patch{patch}.boundary(side) The size of the vector for each interface is rhs_int = zeros(sp_interface.ndof,1); After you compute rhs_int, you have to pass from the local numbering of the interface to the global numbering of the patch (with boundary.dofs), and then from the numbering of the patch to the global multipatch numbering (with gnum) num = space.gnum{patch}(sp_interface.dofs) rhs(num) = rhs(num)  rhs_int; Best regards, Rafa On 21. 06. 17 05:48, tju_zmx wrote: > Dear Rafa > > Now I can get the equivalent load of the interface pressure correctly. > But the problem is that I don't know how to assemble rhs. > How to deal with following code? > > rhs_press = zeros (space.boundary.ndof, 1); > ...... > rhs_press(space.boundary.gnum{iref}) = > rhs_press(space.boundary.gnum{iref}) + op_pn_v (sp_side, msh_side, pval); > ...... > rhs(space.boundary.dofs) = rhs(space.boundary.dofs)  rhs_press; > > Yours Sincerely, > Tony > > 
From: tju_zmx <tju_zmx@12...>  20170621 03:49:14

Dear Rafa Now I can get the equivalent load of the interface pressure correctly. But the problem is that I don't know how to assemble rhs. How to deal with following code? rhs_press = zeros (space.boundary.ndof, 1); ...... rhs_press(space.boundary.gnum{iref}) = rhs_press(space.boundary.gnum{iref}) + op_pn_v (sp_side, msh_side, pval); ...... rhs(space.boundary.dofs) = rhs(space.boundary.dofs)  rhs_press; Yours Sincerely, Tony 
From: Rafael Vazquez <vazquez@im...>  20170620 14:00:59

Dear Tony, there is ongoing research to impose high continuity across patches, either imposing it weakly (with Nitsche's method, for instance), or strongly (see for instance the works by Sangalli and Takacs, Jorg Peters, or Mario Kapl). None of them has been implemented in geopdes, and I am not planning to do it in the near future. About imposing the pressure on the interface, I don't know which kind of problem you want to solve. In any case, you can access the interfaces from msh.msh_patch{patch_number}.boundary(side_number) space.sp_patch{patch_number}.boundary(side_number) Use it at your own risk. ;) Best regards, Rafa On 20. 06. 17 14:15, tju_zmx wrote: > Dear Rafa > > To get the continuous stress, is it possible to control the regularity > of the solution between different patches? > > If I want to apply the pressure on the interface between patches, how > to assemble the righthand side vector correctly? > Do you think it necessary that I create the msh.interface (similar to > msh.boundary)? Is there another way to realise it ? > > Yours Sincerely, > Tony > > 
From: tju_zmx <tju_zmx@12...>  20170620 12:15:52

Dear Rafa To get the continuous stress, is it possible to control the regularity of the solution between different patches? If I want to apply the pressure on the interface between patches, how to assemble the righthand side vector correctly? Do you think it necessary that I create the msh.interface (similar to msh.boundary)? Is there another way to realise it ? Yours Sincerely, Tony 
From: Rafael Vazquez <vazquez@im...>  20170619 15:04:35

Dear Hadi, in this article there is an approximation of the airfoil with Bsplines. http://www.sciencedirect.com/science/article/pii/S0045782514004708?via%3Dihub You will not get an exact representation with Bsplines, but I don't know any work where it has been done with NURBS. Best regards, Rafael 
From: Hadi Minbashian <minbashian@gm...>  20170617 12:54:48

Dear Professor Vazquez, If you could please let me know how I can represent the geometry of an airfoil (e.g. NACA0012) in GeoPDEs either with Bsplines or NURBS. What I have at my hand is a set of points on the boundary of the airfoil and I need to produce a mesh similar to the one depicted in the attached figure in order to do CFD computations. Best regards  Hadi Minbashian (Mr.) PhD Student in Applied Mathematics, Dept. of Mathematics, Faculty of Math. & Comp. Science, Amirkabir University of Technology Hafez Ave., Tehran, Iran Office: Kharazmi Building, Room 338 Tel: (+98)2164542545 Fax: (+98)2166497930 Email: minbashian@... 
From: Rafael Vazquez <vazquez@im...>  20170614 17:04:04

Dear Tony, good to hear that it works. The code looks reasonable to me. About the discontinuous stress, that is normal. Between patches the displacement is only continuous, as in FEM. So the stress, which is computed from the first derivatives, is discontinuous. Best regards, Rafa 
From: tju_zmx <tju_zmx@12...>  20170614 08:01:44

Dear Rafa Firstly, I have successfully applied the pressure boundary condition on multipatch. I also made minor changes to the code as follow. I think it maybe right, is this code reasonable or not? % Apply pressure conditions if (exist ('press_sides', 'var')) for iref = press_sides ... pval = reshape (p (x{:}, iref), msh_side.nqn, msh_side.nel); rhs_press(space.boundary.gnum{iref}) = rhs_press(space.boundary.gnum{iref}) + op_pn_v (sp_side, msh_side, pval); ... end Second, I have another question to consult you. How to deal with discontinuous results of stress field？What factors might contribute to this result (Annex 1)? Yours Sincerely, Tony 
From: tju_zmx <tju_zmx@12...>  20170608 12:04:11

Dear Rafa Thanks for your patience and help. Now that you've got the correct solution, I think there must be something wrong with my operation. I will find out the problem and then solve it. Thank you again professor. Yours Sincerely, Tony 
From: Rafael Vazquez <vazquez@im...>  20170607 08:42:16

Dear Tony, in my case I get the correct solution. Are you sure to be removing the corresponding pressure side from the list of Dirichlet sides? Best, Rafa On 07. 06. 17 05:02, tju_zmx wrote: > Dear Rafa > > I'm sorry for the unclear expression in the previous email. This time > the detailed description of the problems I have encountered is added > in the attachment. > > Yours Sincerely, > Tony 
From: tju_zmx <tju_zmx@12...>  20170607 03:03:38

Dear Rafa I'm sorry for the unclear expression in the previous email. This time the detailed description of the problems I have encountered is added in the attachment. Yours Sincerely, Tony 
From: Rafael Vazquez <vazquez@im...>  20170606 16:35:51

Dear Tony, sorry, but I don't understand the problem. Can you try to explain this with an example, writing explicitly how you choose the cboundary conditions? Best, Rafa 
From: tju_zmx <tju_zmx@12...>  20170606 11:43:58

Dear Rafael I just give the pressure condition like this: problem_data.p = @(x, y, ind) 9.81 * ones(size(x)); What I want to do is not apply the pressure condition on all the sides. For example, if press_sides=1(or other side), the correct result will be calculated. But when I change press_sides into other sides , where other kind of boundary condition are not applied. So I think the code maybe have some details to be modified. Yours Sincerely, Tony At 20170606 18:15:47, "Rafael Vazquez" <vazquez@...> wrote: Dear Tony, the code should work, provided that you give in press_sides all the boundaries where you want to apply the pressure condition. If you apply the pressure condition on all the sides, then your solution is not unique. Best, Rafa On 06. 06. 17 04:21, tju_zmx wrote: Dear Rafael Thanks for your detailed solution and code. I have understood the logical code and have implemented it. Unexpectedly, this code cannot apply boundary condition to all boundaries of multipatch. I tried to modify it, but I found nothing wrong with it. So I have to ask you for help :) Yours Sincerely, Tony 
From: Rafael Vazquez <vazquez@im...>  20170606 10:15:58

Dear Tony, the code should work, provided that you give in press_sides all the boundaries where you want to apply the pressure condition. If you apply the pressure condition on all the sides, then your solution is not unique. Best, Rafa On 06. 06. 17 04:21, tju_zmx wrote: > Dear Rafael > > Thanks for your detailed solution and code. I have understood the > logical code and have implemented it. > Unexpectedly, this code cannot apply boundary condition to > allboundaries of multipatch. I tried to modify it, but I found nothing > wrong with it. > So I have to ask you for help :) > Yours Sincerely, > Tony > > 
From: tju_zmx <tju_zmx@12...>  20170606 02:23:43

Dear Rafael Thanks for your detailed solution and code. I have understood the logical code and have implemented it. Unexpectedly, this code cannot apply boundary condition to all boundaries of multipatch. I tried to modify it, but I found nothing wrong with it. So I have to ask you for help :) Yours Sincerely, Tony 
From: Rafael Vazquez <vazquez@im...>  20170602 13:38:22

Dear Tony, sorry, I now realized that for the pressure conditions you also need the normal vector, which is not computed when you use msh.boundary. You need to do that using msh_eval_boundary_side. This makes the implementation for Neumann and pressure conditions a little bit different. I realized that it was faster to implement it myself than to explain how to do it. So, here is the code. I have also pushed it to the github repository. if (exist ('press_sides', 'var')) for iref = press_sides rhs_press = zeros (space.boundary.ndof, 1); for iside = 1:numel(boundaries(iref).nsides) patch = boundaries(iref).patches(iside); side = boundaries(iref).faces(iside); msh_side = msh_eval_boundary_side (msh.msh_patch{patch}, side); sp_side = sp_eval_boundary_side (space.sp_patch{patch}, msh_side); x = cell (msh_side.rdim, 1); for idim = 1:msh_side.rdim x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel); end pval = reshape (p (x{:}, iside), msh_side.nqn, msh_side.nel); rhs_press(space.boundary.gnum{patch}) = rhs_press(space.boundary.gnum{patch}) + op_pn_v (sp_side, msh_side, pval); end rhs(space.boundary.dofs) = rhs(space.boundary.dofs)  rhs_press; end end Best regards, Rafa 
From: tju_zmx <tju_zmx@12...>  20170601 12:08:41

Dear Rafael To add the pressure conditions of multipatch, I modify the code along the following lines of thinking. But it didn't achieve the desired result. So please help me to check if there is any mistake in the train of thought. Firstly, I add following codes into the function mp_solve_linear_elasticity. Nbnd = cumsum ([0, boundaries.nsides]); for iref = press_sides iref_patch_list = Nbnd(iref)+1:Nbnd(iref+1); pref = @(varargin) p(varargin{:},iref); rhs_press = op_pn_v_mp (space.boundary, msh.boundary, pref, iref_patch_list); rhs(space.boundary.dofs) = rhs(space.boundary.dofs)  rhs_press; end Second, for convenience of calculations, functions op_f_v_mp and op_f_v_tp were replaced by op_pn_v_mp and op_pn_v_tp. op_pn_v_mp is very similar to op_f_v_mp, op_pn_v_tp is a little different to op_f_v_tp. I modify op_f_v_tp use following code: x = cell (msh.rdim, 1); for idim = 1:msh.rdim x{idim} = squeeze (msh_col.geo_map(idim,:,:)); end pval = reshape (press (x{:}), msh_col.nqn, msh_col.nel); rhs = rhs + op_pn_v (sp_col, msh_col, pval); Is this idea feasible? If not, please point out the mistakes and I will make a further revision. Yours Sincerely, Tony 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20170531 11:48:01

I have not looked into the problem in detail, but my "sources" in the transport community assure me that the calculation proceeds "elementbyelement". I think that the problem is not quite so complicated as it sounds:  suppose there are two neighbor patches, P1 and P2, and the boundary, B12. Suppose 2D geometry.  On P1, the quadrature nodes are given as uq1, vq1, and the weights, uw1 and vw1; in physical space, the quadrature nodes are indicated as p1(x, y)  On P2, the quadrature nodes are given as uq2, vq2, and the weights, uw2 and vw2; in physical space, the quadrature nodes are indicated as p2(x, y) Suppose that B12 corresponds to the boundary (u1 = 1, v1 = v) and (u2 = 0, v2 = v)  the solution proceeds, and suppose that the solution for P1 has been determined; thus the solution can be "reconstructed" on any point (x, y), including the boundary B12  the only thing you need to know is the physical location of the quadrature nodes for P2, as applied to the "incoming" boundary B12. I think that the "incoming" flux can be created with an interpolation of the "outgoing" flux  but again, I have not yet written down the equations. Later,Wilfred On Wednesday, May 31, 2017 12:33 AM, Rafael Vazquez <vazquez@...> wrote: Hello Wilfred, yes, I know a bit about DG. But as far as I know you don't compute things element by element. What you do is to add in the variational formulation some terms containing the jumps and the averages at the edges, but at the end you need to solve a global system. In any case, I know that the group of Ulrich Langer has been applying DG techniques between patches in IGA. In GeoPDEs there is something for the divconforming splines in Stokes problem, in the files "mp_solve_stokes_div_conforming" and "mp_dg_penalty". All the computations there assume that the mesh is the same in both sides. I don't know if that can be of any help. Best, Rafa On 30. 05. 17 14:20, Wilfred van Rooijen wrote: Hello Rafael, Thank you for your reply, even if it came a bit later than expected :)) Yes, the method is well known, and it is commonly called "Discontinuous Galerkin Finite Element Method" (DGFEM)  at least in the universe of the people working with neutron transport theory. I can send you some references later, I don't have my notes at hand. In the mean time, I programmed most of the DGFEM in FORTRAN, and I have GeoPDEs extensively, to create the NURBS domains, and also to check my matrix calculations. Without GeoPDEs it would have been impossible :)) So, onward I go, into the Great Unknown. Wilfred On Tuesday, May 30, 2017 7:05 PM, Rafael Vazquez <vazquez@...> wrote: Hi Wilfred, sorry for the late response, your message finished in the spam folder, and I only saw it very recently. :( To be honest, I would not call the method you describe a "finite element method". In FEM, you solve your problem in the global domain at once, and not element by element. Do you have a reference for that method? To me it resembles more the "finite volume method" (although I'm not an expert on that). About your question, as you said, if the two patches match conformingly (same knot vector, same control points), then the basis functions from each side coincide exactly, and the coefficients for both sides are the same. Instead, if the two meshes do not match, you have to use some kind of projector. In the IGA literature probably the most common way to do that is to use mortar methods or Lagrange multipliers. Both of them will require to compute your fields on an intermediate mesh. This is not difficult in 2D, but may become more complicated in 3D. There is also the possibility to pass your solution from one side to the next one using quasiinterpolants, but in any case, you need to compute the functions from one side in the next one, and they use different meshes. The SP_DRCHLT_L2_PROJ is a particular implementation, where we impose boundary conditions using the L^2 projection. There are two important things to note. First, the function works for an input function "h" that evaluates the function you want to project. In your case, you would need to change this function. Second, the L^2 projection is not a good choice for the transport equation, because it is not local. I can't recommend that. I hope that helps. Best, Rafa 
From: Rafael Vazquez <vazquez@im...>  20170530 15:33:59

Hello Wilfred, yes, I know a bit about DG. But as far as I know you don't compute things element by element. What you do is to add in the variational formulation some terms containing the jumps and the averages at the edges, but at the end you need to solve a global system. In any case, I know that the group of Ulrich Langer has been applying DG techniques between patches in IGA. In GeoPDEs there is something for the divconforming splines in Stokes problem, in the files "mp_solve_stokes_div_conforming" and "mp_dg_penalty". All the computations there assume that the mesh is the same in both sides. I don't know if that can be of any help. Best, Rafa On 30. 05. 17 14:20, Wilfred van Rooijen wrote: > Hello Rafael, > > Thank you for your reply, even if it came a bit later than expected :)) > > Yes, the method is well known, and it is commonly called > "Discontinuous Galerkin Finite Element Method" (DGFEM)  at least in > the universe of the people working with neutron transport theory. I > can send you some references later, I don't have my notes at hand. > > In the mean time, I programmed most of the DGFEM in FORTRAN, and I > have GeoPDEs extensively, to create the NURBS domains, and also to > check my matrix calculations. Without GeoPDEs it would have been > impossible :)) > > So, onward I go, into the Great Unknown. > > Wilfred > > > On Tuesday, May 30, 2017 7:05 PM, Rafael Vazquez > <vazquez@...> wrote: > > > Hi Wilfred, > sorry for the late response, your message finished in the spam folder, > and I only saw it very recently. :( > > To be honest, I would not call the method you describe a "finite element > method". In FEM, you solve your problem in the global domain at once, > and not element by element. Do you have a reference for that method? To > me it resembles more the "finite volume method" (although I'm not an > expert on that). > > About your question, as you said, if the two patches match conformingly > (same knot vector, same control points), then the basis functions from > each side coincide exactly, and the coefficients for both sides are the > > same. > > > Instead, if the two meshes do not match, you have to use some kind of > projector. In the IGA literature probably the most common way to do that > is to use mortar methods or Lagrange multipliers. Both of them will > require to compute your fields on an intermediate mesh. This is not > difficult in 2D, but may become more complicated in 3D. There is also > the possibility to pass your solution from one side to the next one > using quasiinterpolants, but in any case, you need to compute the > functions from one side in the next one, and they use different meshes. > > The SP_DRCHLT_L2_PROJ is a particular implementation, where we impose > boundary conditions using the L^2 projection. There are two important > things to note. First, the function works for an input function "h" that > evaluates the function you want to project. In your case, you would need > to change this function. Second, the L^2 projection is not a good choice > for the transport equation, because it is not local. I can't recommend > that. > > I hope that helps. > Best, > Rafa > > > 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20170530 12:24:44

Hello Rafael, Thank you for your reply, even if it came a bit later than expected :)) Yes, the method is well known, and it is commonly called "Discontinuous Galerkin Finite Element Method" (DGFEM)  at least in the universe of the people working with neutron transport theory. I can send you some references later, I don't have my notes at hand. In the mean time, I programmed most of the DGFEM in FORTRAN, and I have GeoPDEs extensively, to create the NURBS domains, and also to check my matrix calculations. Without GeoPDEs it would have been impossible :)) So, onward I go, into the Great Unknown. Wilfred On Tuesday, May 30, 2017 7:05 PM, Rafael Vazquez <vazquez@...> wrote: Hi Wilfred, sorry for the late response, your message finished in the spam folder, and I only saw it very recently. :( To be honest, I would not call the method you describe a "finite element method". In FEM, you solve your problem in the global domain at once, and not element by element. Do you have a reference for that method? To me it resembles more the "finite volume method" (although I'm not an expert on that). About your question, as you said, if the two patches match conformingly (same knot vector, same control points), then the basis functions from each side coincide exactly, and the coefficients for both sides are the same. Instead, if the two meshes do not match, you have to use some kind of projector. In the IGA literature probably the most common way to do that is to use mortar methods or Lagrange multipliers. Both of them will require to compute your fields on an intermediate mesh. This is not difficult in 2D, but may become more complicated in 3D. There is also the possibility to pass your solution from one side to the next one using quasiinterpolants, but in any case, you need to compute the functions from one side in the next one, and they use different meshes. The SP_DRCHLT_L2_PROJ is a particular implementation, where we impose boundary conditions using the L^2 projection. There are two important things to note. First, the function works for an input function "h" that evaluates the function you want to project. In your case, you would need to change this function. Second, the L^2 projection is not a good choice for the transport equation, because it is not local. I can't recommend that. I hope that helps. Best, Rafa 
From: Rafael Vazquez <vazquez@im...>  20170530 10:05:57

Hi Wilfred, sorry for the late response, your message finished in the spam folder, and I only saw it very recently. :( To be honest, I would not call the method you describe a "finite element method". In FEM, you solve your problem in the global domain at once, and not element by element. Do you have a reference for that method? To me it resembles more the "finite volume method" (although I'm not an expert on that). About your question, as you said, if the two patches match conformingly (same knot vector, same control points), then the basis functions from each side coincide exactly, and the coefficients for both sides are the same. Instead, if the two meshes do not match, you have to use some kind of projector. In the IGA literature probably the most common way to do that is to use mortar methods or Lagrange multipliers. Both of them will require to compute your fields on an intermediate mesh. This is not difficult in 2D, but may become more complicated in 3D. There is also the possibility to pass your solution from one side to the next one using quasiinterpolants, but in any case, you need to compute the functions from one side in the next one, and they use different meshes. The SP_DRCHLT_L2_PROJ is a particular implementation, where we impose boundary conditions using the L^2 projection. There are two important things to note. First, the function works for an input function "h" that evaluates the function you want to project. In your case, you would need to change this function. Second, the L^2 projection is not a good choice for the transport equation, because it is not local. I can't recommend that. I hope that helps. Best, Rafa 
From: Rafael Vazquez <vazquez@im...>  20170530 09:30:35

Dear Tony, if I remember correctly, it should not be difficult to modify the function mp_solve_linear_elasticity to add the pressure conditions. The changes are very similar to those already done to apply Neumann conditions. Please, be more specific if you encounter a problem to implement them. Best, Rafa 
From: tju_zmx <tju_zmx@12...>  20170527 08:55:13

Dear geopdes users： For the 2D linear elasticity problem, how can I apply pressure conditions on the multipatch model? How to modify the function mp_solve_linear_elasticity (in folder multipatch) to implement it? Is there any thoughts to solve this kind of problem? Any suggestions are gratefully accepted. Yours Sincerely, Tony. 
From: tju_zmx <tju_zmx@12...>  20170526 11:13:37

Dear professor, I appreciate your detail explanation and patiently answer. I will try to solve the problems encountered. Best regards, Tony At 20170526 16:18:15, "Rafael Vazquez" <vazquez@...> wrote: Dear Tony, for the first, I don't know exactly where you are using it, but the best thing is to do as you did for the assembly operator: inside the loop for the patches, generate the function handles lambda_patch and mu_patch, and then call sp_eval_msh using these two handles, instead of the original ones. For the second, the variable "ind" indicates the number of the boundary side. For different constant pressures, you can store the pressures in an array, and then use problem_data.p = @(x, y, ind) P(ind) * ones (size (x)); If you want to apply nonconstant pressures, you can use a different file for the function, like in the Neumann condition (g) of the example "ex_laplace_ring_mixed_bc". Or you can use logic operators to check the boundary side, for instance: problem_data.p = @(x,y,ind) 10 * ones(size(x)) * (ind == 3) + y.^2 * (ind == 1); Best, Rafa 
From: Rafael Vazquez <vazquez@im...>  20170526 08:18:28

Dear Tony, for the first, I don't know exactly where you are using it, but the best thing is to do as you did for the assembly operator: inside the loop for the patches, generate the function handles lambda_patch and mu_patch, and then call sp_eval_msh using these two handles, instead of the original ones. For the second, the variable "ind" indicates the number of the boundary side. For different constant pressures, you can store the pressures in an array, and then use problem_data.p = @(x, y, ind) P(ind) * ones (size (x)); If you want to apply nonconstant pressures, you can use a different file for the function, like in the Neumann condition (g) of the example "ex_laplace_ring_mixed_bc". Or you can use logic operators to check the boundary side, for instance: problem_data.p = @(x,y,ind) 10 * ones(size(x)) * (ind == 3) + y.^2 * (ind == 1); Best, Rafa 