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}
(14) 
_{Oct}

_{Nov}

_{Dec}

From: Rafael Vazquez <vazquez@im...>  20160912 11:22:50

Hi Are, to be honest, I don't remember why it is commented. The same lines are uncommented in "msh_evaluate_col" and "msh_evaluate_element_list". I'll check it, in the meantime you can uncomment them, at your own risk. ;) Note that we are computing the element length (the standard "h" in FEM). To compute the area or volume, you should elevate that quantity to ndim. Or you can simply use Wilfred's suggestion, summing over the quadrature points: elem_area = sum (msh.quad_weights .* msh.jacdet, 1); Best, Rafa On 12/09/2016 11:28, Are Simonsen wrote: > Hi Rafael, > > Thanks for the answers you have provided. They have been valuable! > One more question: I am interested in the surface area vector of the elements. I have the normal, but now I need the element area as well. I see from the manual that if calling the "old" versions routines through "msh_precompute" I can get the element size (which I assume that for a 2D side is the area of the element in the physical domain...). However when calling msh_"precompute" I did not get this field and looking into the function (mfile) I see that the calculation of the "element_size" is commented out. Could you comment on this? Is there another way to calculate the element area already in the code? > > Best regards, > Are > > > ####################################### > #### Are J. Simonsen, PhD. > #### Research scientist > #### > #### SINTEF MATERIALS&CHEMISTRY, Process technology > #### SP Andersens v. 15B, room 305, 3'rd floor > #### NO7034 Trondheim, Norway > #### > #### Email: Are.Simonsen@...<mailto:Are.Simonsen@...> > #### Phone: +47 92424656 (M) > ####################################### > >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20160912 10:58:55

This may be a stupid suggestion, so Rafael, please correct me if I'm wrong. In the function sp_l2_err.m, the solution is integrated over the entire space: >>snip<< w = msh.quad_weights .* msh.jacdet; errl2_elem = sum (reshape (sum ((valu  valex).^2, 1), [msh.nqn, msh.nel]) .* w); >>end<< It would appear to me that the "w" are the surface area of each element (in the physical space). Wilfred On Monday, September 12, 2016 6:28 PM, Are Simonsen <Are.Simonsen@...> wrote: Hi Rafael, Thanks for the answers you have provided. They have been valuable! One more question: I am interested in the surface area vector of the elements. I have the normal, but now I need the element area as well. I see from the manual that if calling the "old" versions routines through "msh_precompute" I can get the element size (which I assume that for a 2D side is the area of the element in the physical domain...). However when calling msh_"precompute" I did not get this field and looking into the function (mfile) I see that the calculation of the "element_size" is commented out. Could you comment on this? Is there another way to calculate the element area already in the code? Best regards, Are ####################################### #### Are J. Simonsen, PhD. #### Research scientist #### #### SINTEF MATERIALS&CHEMISTRY, Process technology #### SP Andersens v. 15B, room 305, 3'rd floor #### NO7034 Trondheim, Norway #### #### Email: Are.Simonsen@...<mailto:Are.Simonsen@...> #### Phone: +47 92424656 (M) #######################################  _______________________________________________ Geopdesusers mailing list Geopdesusers@... https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Are Simonsen <Are.S<imonsen@si...>  20160912 09:28:21

Hi Rafael, Thanks for the answers you have provided. They have been valuable! One more question: I am interested in the surface area vector of the elements. I have the normal, but now I need the element area as well. I see from the manual that if calling the "old" versions routines through "msh_precompute" I can get the element size (which I assume that for a 2D side is the area of the element in the physical domain...). However when calling msh_"precompute" I did not get this field and looking into the function (mfile) I see that the calculation of the "element_size" is commented out. Could you comment on this? Is there another way to calculate the element area already in the code? Best regards, Are ####################################### #### Are J. Simonsen, PhD. #### Research scientist #### #### SINTEF MATERIALS&CHEMISTRY, Process technology #### SP Andersens v. 15B, room 305, 3'rd floor #### NO7034 Trondheim, Norway #### #### Email: Are.Simonsen@...<mailto:Are.Simonsen@...> #### Phone: +47 92424656 (M) ####################################### 
From: Rafael Vázquez <vazquez@im...>  20160911 17:37:55

Hi Are, the computation in *) is wrong. The vector of degrees of freedom should be associated to the space, which is the boundary space in this case. So, you need to collect the boundary degrees of freedom, like this: > ind = space.sp_patch{iptc}.boundary(iside).dofs; > indices = space.gnum{iptc}(ind); > [eu2, F2] = sp_eval_msh (u(indices), sp_side,msh_side); I have decided to add an error message when the length of the vector U and the number of degrees of freedom (sp_side.ndof) are different. With this change, the value of the function computed with *) and **) should be the same, although arranged differently. However, if you try to compute the gradient in a similar way, *) will give you the surface gradient, and **) the "volumetric" gradient, including the normal component, which is the one you want. To get the values arranged in the same way, instead of **) you can use the mesh given in "msh_boundary_side_from_interior", as in the example of sp_weak_drchlt_bc. Best, Rafa > Hi Rafael, > > Thank you for your prompt response late on a Friday evening...:) Amazing! > > I would definitely prefer to use the global numbering as a black box(I did > look into the files you mentioned and it would take some effort to fully > get the grasp of it...), but I was a little unsure on what to put into the > functions working only at a boundary. I have the lines > > *) > msh_side = msh_eval_boundary_side(msh_iptc, iside); > sp_side = sp_eval_boundary_side (space.sp_patch{iptc}, msh_side);% > [eu2, F2] = sp_eval_msh (u(space.gnum{iptc}), sp_side,msh_side); > *) > To calculate values on a side. I was unsure on the values going into > "sp_eval_msh" u(space.gnum{iptc}, where iptc is the patch number. If I > understand you right all the degrees of freedom. U values, should be given > as input and then sp_eval_msh" finds those values it needs by itself... > From the lines above I got the normals, F and eu2, but the eu2 values I > get are not the correct ones. I get the correct ones when I use > > **) > [eu, F] = sp_eval(u(space.gnum{iptc}), space.sp_patch{iptc}, > geometry(iptc), vtk_pts); > [grad_eu, F] = sp_eval (u(space.gnum{iptc}), space.sp_patch{iptc}, > geometry(iptc), vtk_pts,'gradient'); > **) > > Where vtk_pts are defined for the desired side vtk_pts={0 [xi_1] [xi_2]}. > The F I get from both *, and ** lines are the same (although they are > arranged differently). However the eu and eu2 values, that I expected to > be the same are not! This is more than likely due to a mistake by me, and > I thought that one source where that I fed the wrong u values into the > function... > > I see that you suggest to use the "sp_precompute" and " sp_eval_msh" so I > will look into that next. > > Thanks again! > > Best regards, > Are > > > > > Original Message > From: Rafael Vázquez [mailto:vazquez@...] > Sent: 9. september 2016 22:21 > To: Are Simonsen <Are.Simonsen@...> > Cc: Geopdesusers@... > Subject: Re: [Geopdesusers] solution gradient times geometry gradient at > a given side in a multipatch geometry > > By the way, since it may be useful for other things, let me explain how to > get the degrees of freedom of a boundary side for a multipatch domain. > > Assume that your patch is "iptc", and the side number in the patch is > "iside". First, you pass from the boundary to the interior of the patch > using boundary.dofs, and then from the patch numbering to the global > numbering using "gnum": > > ind = space.sp_patch{iptc}.boundary(iside).dofs; > indices = space.gnum{iptc}(ind); > > About the rule to define the global numbering, I first number the internal > degrees of freedom of every patch, and then the interfaces, recursively > searching for dofs shared by more than two patches. If you want to see the > details you can try to dig into the "mp_interface" function, but I would > use it as a black box. > > Best, > Rafa > > 
From: Are Simonsen <Are.S<imonsen@si...>  20160909 20:43:32

Hi Rafael, Thank you for your prompt response late on a Friday evening...:) Amazing! I would definitely prefer to use the global numbering as a black box(I did look into the files you mentioned and it would take some effort to fully get the grasp of it...), but I was a little unsure on what to put into the functions working only at a boundary. I have the lines *) msh_side = msh_eval_boundary_side(msh_iptc, iside); sp_side = sp_eval_boundary_side (space.sp_patch{iptc}, msh_side);% [eu2, F2] = sp_eval_msh (u(space.gnum{iptc}), sp_side,msh_side); *) To calculate values on a side. I was unsure on the values going into "sp_eval_msh" u(space.gnum{iptc}, where iptc is the patch number. If I understand you right all the degrees of freedom. U values, should be given as input and then sp_eval_msh" finds those values it needs by itself... >From the lines above I got the normals, F and eu2, but the eu2 values I get are not the correct ones. I get the correct ones when I use **) [eu, F] = sp_eval(u(space.gnum{iptc}), space.sp_patch{iptc}, geometry(iptc), vtk_pts); [grad_eu, F] = sp_eval (u(space.gnum{iptc}), space.sp_patch{iptc}, geometry(iptc), vtk_pts,'gradient'); **) Where vtk_pts are defined for the desired side vtk_pts={0 [xi_1] [xi_2]}. The F I get from both *, and ** lines are the same (although they are arranged differently). However the eu and eu2 values, that I expected to be the same are not! This is more than likely due to a mistake by me, and I thought that one source where that I fed the wrong u values into the function... I see that you suggest to use the "sp_precompute" and " sp_eval_msh" so I will look into that next. Thanks again! Best regards, Are Original Message From: Rafael Vázquez [mailto:vazquez@...] Sent: 9. september 2016 22:21 To: Are Simonsen <Are.Simonsen@...> Cc: Geopdesusers@... Subject: Re: [Geopdesusers] solution gradient times geometry gradient at a given side in a multipatch geometry By the way, since it may be useful for other things, let me explain how to get the degrees of freedom of a boundary side for a multipatch domain. Assume that your patch is "iptc", and the side number in the patch is "iside". First, you pass from the boundary to the interior of the patch using boundary.dofs, and then from the patch numbering to the global numbering using "gnum": ind = space.sp_patch{iptc}.boundary(iside).dofs; indices = space.gnum{iptc}(ind); About the rule to define the global numbering, I first number the internal degrees of freedom of every patch, and then the interfaces, recursively searching for dofs shared by more than two patches. If you want to see the details you can try to dig into the "mp_interface" function, but I would use it as a black box. Best, Rafa 
From: Rafael Vázquez <vazquez@im...>  20160909 20:21:07

By the way, since it may be useful for other things, let me explain how to get the degrees of freedom of a boundary side for a multipatch domain. Assume that your patch is "iptc", and the side number in the patch is "iside". First, you pass from the boundary to the interior of the patch using boundary.dofs, and then from the patch numbering to the global numbering using "gnum": ind = space.sp_patch{iptc}.boundary(iside).dofs; indices = space.gnum{iptc}(ind); About the rule to define the global numbering, I first number the internal degrees of freedom of every patch, and then the interfaces, recursively searching for dofs shared by more than two patches. If you want to see the details you can try to dig into the "mp_interface" function, but I would use it as a black box. Best, Rafa 
From: Rafael Vázquez <vazquez@im...>  20160909 20:11:05

Hi Are, this is trickier than it seems, even for a single patch. From boundary objects you can only compute quantities that you get from boundary information. For instance, when computing the gradient you will get the surface gradient (equal to n \times grad u \times n), and the normal component will be lost. I tried to explain this (not in detail) in Remark 3.5 of the new paper. So, to compute the gradient including the normal component you need to define an auxiliary "mesh" in the whole domain, setting one of the parametric coordinates of the points to 0 or 1, depending on the side. This is already done in "msh_boundary_side_from_interior". Then, you have to construct the space on this mesh, using the space constructor. You have an example of this procedure in @sp_vector/sp_weak_drchlt_bc, or in @sp_multipatch/sp_weak_drchlt_bc. The example there is for vector valued spaces, but it works exactly the same for scalars. The last step is to evaluate your field using sp_eval_msh. Since this function works with structures, and not with objects of the class, you need to call msh_precompute and sp_precompute before that. The two structures that you get as output can be passed to sp_eval_msh, and you will get the gradient evaluated at the boundary points, including the normal components. Note that you are using objects defined in the whole domain, so you need to pass the solution vector "u" with all the coefficients, without restricting it to the boundary. I would like to make this thing easier in a future version, but probably it will take time. Best regards, Rafa > I am doing some multiplication of the gradient of a scalar with the normal > vector of the surface in a multipath geometry. I see that I can use the > "msh_eval_boundary_side" to calculate the normals at the desired side > while "sp_eval_boundary_side" can be used to calculate the values and the > geometry points. As far as I understand I will then have to give the > correct weights for the solution (lets call these weights u) for that > specific side into "sp_eval_boundary_side". That is I have to extract the > correct local values for the side from the global solution weights, u. If > this is correct how do I find the correct u's for a given side in a > multipatch geometry. If I am not correct and there is a much easier way to > do this I would appreciate your advice! > > I know another strategy is to use "sp_eval" at the correct side. Then I > would have to calculate the normals of the geometry. This is done in > already in "msh_eval_boundary_side" which is why I thought this would be a > good way to do this... > > What is the rule for global numbering of solution weights (u's) in a mp > geometry? > > Best regards, > Are Simonsen > > > > ####################################### > #### Are J. Simonsen, PhD. > #### Research scientist > #### > #### SINTEF MATERIALS&CHEMISTRY, Process technology > #### SP Andersens v. 15B, room 305, 3'rd floor > #### NO7034 Trondheim, Norway > #### > #### Email: > Are.Simonsen@...<mailto:Are.Simonsen@...> > #### Phone: +47 92424656 (M) > ####################################### > >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers > 
From: Are Simonsen <Are.S<imonsen@si...>  20160909 16:35:44

I am doing some multiplication of the gradient of a scalar with the normal vector of the surface in a multipath geometry. I see that I can use the "msh_eval_boundary_side" to calculate the normals at the desired side while "sp_eval_boundary_side" can be used to calculate the values and the geometry points. As far as I understand I will then have to give the correct weights for the solution (lets call these weights u) for that specific side into "sp_eval_boundary_side". That is I have to extract the correct local values for the side from the global solution weights, u. If this is correct how do I find the correct u's for a given side in a multipatch geometry. If I am not correct and there is a much easier way to do this I would appreciate your advice! I know another strategy is to use "sp_eval" at the correct side. Then I would have to calculate the normals of the geometry. This is done in already in "msh_eval_boundary_side" which is why I thought this would be a good way to do this... What is the rule for global numbering of solution weights (u's) in a mp geometry? Best regards, Are Simonsen ####################################### #### Are J. Simonsen, PhD. #### Research scientist #### #### SINTEF MATERIALS&CHEMISTRY, Process technology #### SP Andersens v. 15B, room 305, 3'rd floor #### NO7034 Trondheim, Norway #### #### Email: Are.Simonsen@...<mailto:Are.Simonsen@...> #### Phone: +47 92424656 (M) ####################################### 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20160908 14:20:00

Hey, I see my name being used in a message to this list. I feel honored ;)) To be honest, my case was a littlebit different: in my case, the source term in each step depends on the previous solution, and in the present case, the diffusion coefficient depends on the previous solution (I think  correct me if I'm wrong). The solution is not too difficult. Normally, the function op_gradu_gradv_tp has an (optional) argument "coeff" which is a function to calculate the diffusion coefficient. What you should do, is modify this routine as follows: When you know u_k, you can use the routine [eu, F] = sp_eval_msh() to determine the solution "eu" ("evaluated u") at each quadrature point. Then you must write a routine to calculate q(u_k) at each quadratute point. Finally, you should modify op_gradu_gradv_tp so that the argument "coeff" is not a function, but your calculated matrix q(u_k). Anyway  if you read the source code, it is not too difficult to understand (I think). Wilfred On Thursday, September 8, 2016 10:48 PM, Martina Bekrová <to.zapomenu@...> wrote: I'm sorry again... This happens to me all the time :/ Thank you for summarising and suggestions. I'll take a closer look at it later. I'm not a programmer, so it seems difficult for me, but I can try, if it's not too hard. I'll be happy to do something new. Thank you Martina P.S.: TeX file included for everyone 20160908 15:33 GMT+02:00 Rafael Vazquez <vazquez@...>: > Hi Martina, > thanks for the TeX file, but you only sent it to me. Let me summarize > your problem for other users. > > Given $u_k$, you want to compute $u_{k+1} = u_k + \delta u$, where > $\delta_u$ is solution of: > Find $\delta u$ in $H^1_0$ such that for all $v \in H^1_0$ holds > $a(\delta u, v) = L(v)$, where > > a(\delta u, v) = \int_\Omega q(u_k)\nabla \delta u \cdot \nabla v > (q(u_k))^3(\nabla u_k \cdot \nabla \delta u)( \nabla u_k \cdot \nabla v) > L(v) = \int_\Omega q(u_k) \nabla u_k \cdot \nabla v > > and q(u) = \frac{1}{\sqrt{1 + \nabla u^2}}. > > Let's analyze the problem term by term: > > 1) \int_\Omega q(u_k)\nabla \delta u \cdot \nabla v > > Here you would have to modify the function op_gradu_gradv_tp, to > consider a coefficient (q(u_k)) that depends on a previous solution. I > think Wilfred has already implemented this (see our discussion in the > mailing list, it was around May). Maybe he wants to share the code with > you. > > 2) \int_\Omega (q(u_k))^3 (\nabla u_k \cdot \nabla \delta u)( \nabla u_k > \cdot \nabla v) > > We don't have this operator in GeoPDEs. You will have to implement it, > but it shouldn't be difficult. First, you can probably reuse parts of > other operators (in particular op_vel_dot_gradu_v) to make it work for a > constant u_k. Then, you can make it work for a nonconstant coefficient, > that can be computed as in 1). I would also allow to pass two different > coefficients, one for q(u_k)^3, and the other for u_k, so the same code > can be used for other problems. > > 3) L(v) = \int_\Omega q(u_k) \nabla u_k \cdot \nabla v > > This operator is also missing, but you can avoid its implementation. > After you compute the matrix in step 1), you can multiply by the vector > u_k to get L(v). Of course, you will have to do it at the right moment, > to compute the integrals with q(u_k), and not q(u_{k1}). > > I hope this helps. Let me know if you need more details to start the > implementation. > > Best, > Rafa > > On 08/09/2016 12:22, Martina Bekrová wrote: > > Hello > > > > I would like to use GeoPDEs to solve Dirichlet problem for minimal > surface > > equation using Newton method, just nonlinear poisson equation, but I'm > not > > able to find out how to assemble the matrices for this variational > problem. > > Description of the problem is attached below. > > > > I already have this program in Fenics, but I would like to apply IGA for > > this problem too. Is here a way how I can do this? > > > > Thank you > > Martina > > > > > >  >  > > > > > > _______________________________________________ > > Geopdesusers mailing list > > Geopdesusers@... > > https://lists.sourceforge.net/lists/listinfo/geopdesusers > >  >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers >  _______________________________________________ Geopdesusers mailing list Geopdesusers@... https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Martina Bekrová <to.zapomenu@gm...>  20160908 13:48:31

I'm sorry again... This happens to me all the time :/ Thank you for summarising and suggestions. I'll take a closer look at it later. I'm not a programmer, so it seems difficult for me, but I can try, if it's not too hard. I'll be happy to do something new. Thank you Martina P.S.: TeX file included for everyone 20160908 15:33 GMT+02:00 Rafael Vazquez <vazquez@...>: > Hi Martina, > thanks for the TeX file, but you only sent it to me. Let me summarize > your problem for other users. > > Given $u_k$, you want to compute $u_{k+1} = u_k + \delta u$, where > $\delta_u$ is solution of: > Find $\delta u$ in $H^1_0$ such that for all $v \in H^1_0$ holds > $a(\delta u, v) = L(v)$, where > > a(\delta u, v) = \int_\Omega q(u_k)\nabla \delta u \cdot \nabla v > (q(u_k))^3(\nabla u_k \cdot \nabla \delta u)( \nabla u_k \cdot \nabla v) > L(v) = \int_\Omega q(u_k) \nabla u_k \cdot \nabla v > > and q(u) = \frac{1}{\sqrt{1 + \nabla u^2}}. > > Let's analyze the problem term by term: > > 1) \int_\Omega q(u_k)\nabla \delta u \cdot \nabla v > > Here you would have to modify the function op_gradu_gradv_tp, to > consider a coefficient (q(u_k)) that depends on a previous solution. I > think Wilfred has already implemented this (see our discussion in the > mailing list, it was around May). Maybe he wants to share the code with > you. > > 2) \int_\Omega (q(u_k))^3 (\nabla u_k \cdot \nabla \delta u)( \nabla u_k > \cdot \nabla v) > > We don't have this operator in GeoPDEs. You will have to implement it, > but it shouldn't be difficult. First, you can probably reuse parts of > other operators (in particular op_vel_dot_gradu_v) to make it work for a > constant u_k. Then, you can make it work for a nonconstant coefficient, > that can be computed as in 1). I would also allow to pass two different > coefficients, one for q(u_k)^3, and the other for u_k, so the same code > can be used for other problems. > > 3) L(v) = \int_\Omega q(u_k) \nabla u_k \cdot \nabla v > > This operator is also missing, but you can avoid its implementation. > After you compute the matrix in step 1), you can multiply by the vector > u_k to get L(v). Of course, you will have to do it at the right moment, > to compute the integrals with q(u_k), and not q(u_{k1}). > > I hope this helps. Let me know if you need more details to start the > implementation. > > Best, > Rafa > > On 08/09/2016 12:22, Martina Bekrová wrote: > > Hello > > > > I would like to use GeoPDEs to solve Dirichlet problem for minimal > surface > > equation using Newton method, just nonlinear poisson equation, but I'm > not > > able to find out how to assemble the matrices for this variational > problem. > > Description of the problem is attached below. > > > > I already have this program in Fenics, but I would like to apply IGA for > > this problem too. Is here a way how I can do this? > > > > Thank you > > Martina > > > > > >  >  > > > > > > _______________________________________________ > > Geopdesusers mailing list > > Geopdesusers@... > > https://lists.sourceforge.net/lists/listinfo/geopdesusers > >  >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers > 
From: Rafael Vazquez <vazquez@im...>  20160908 13:33:32

Hi Martina, thanks for the TeX file, but you only sent it to me. Let me summarize your problem for other users. Given $u_k$, you want to compute $u_{k+1} = u_k + \delta u$, where $\delta_u$ is solution of: Find $\delta u$ in $H^1_0$ such that for all $v \in H^1_0$ holds $a(\delta u, v) = L(v)$, where a(\delta u, v) = \int_\Omega q(u_k)\nabla \delta u \cdot \nabla v (q(u_k))^3(\nabla u_k \cdot \nabla \delta u)( \nabla u_k \cdot \nabla v) L(v) = \int_\Omega q(u_k) \nabla u_k \cdot \nabla v and q(u) = \frac{1}{\sqrt{1 + \nabla u^2}}. Let's analyze the problem term by term: 1) \int_\Omega q(u_k)\nabla \delta u \cdot \nabla v Here you would have to modify the function op_gradu_gradv_tp, to consider a coefficient (q(u_k)) that depends on a previous solution. I think Wilfred has already implemented this (see our discussion in the mailing list, it was around May). Maybe he wants to share the code with you. 2) \int_\Omega (q(u_k))^3 (\nabla u_k \cdot \nabla \delta u)( \nabla u_k \cdot \nabla v) We don't have this operator in GeoPDEs. You will have to implement it, but it shouldn't be difficult. First, you can probably reuse parts of other operators (in particular op_vel_dot_gradu_v) to make it work for a constant u_k. Then, you can make it work for a nonconstant coefficient, that can be computed as in 1). I would also allow to pass two different coefficients, one for q(u_k)^3, and the other for u_k, so the same code can be used for other problems. 3) L(v) = \int_\Omega q(u_k) \nabla u_k \cdot \nabla v This operator is also missing, but you can avoid its implementation. After you compute the matrix in step 1), you can multiply by the vector u_k to get L(v). Of course, you will have to do it at the right moment, to compute the integrals with q(u_k), and not q(u_{k1}). I hope this helps. Let me know if you need more details to start the implementation. Best, Rafa On 08/09/2016 12:22, Martina Bekrová wrote: > Hello > > I would like to use GeoPDEs to solve Dirichlet problem for minimal surface > equation using Newton method, just nonlinear poisson equation, but I'm not > able to find out how to assemble the matrices for this variational problem. > Description of the problem is attached below. > > I already have this program in Fenics, but I would like to apply IGA for > this problem too. Is here a way how I can do this? > > Thank you > Martina > > >  > > > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Rafael Vazquez <vazquez@im...>  20160908 11:02:18

Hi Martina, I will try to fix the problem with PDF files. I can see them as the list administrator, but other users can't. Can you please send again the variational formulation of your problem, written as TeX? Best, Rafa On 08/09/2016 12:22, Martina Bekrová wrote: > Hello > > I would like to use GeoPDEs to solve Dirichlet problem for minimal surface > equation using Newton method, just nonlinear poisson equation, but I'm not > able to find out how to assemble the matrices for this variational problem. > Description of the problem is attached below. > > I already have this program in Fenics, but I would like to apply IGA for > this problem too. Is here a way how I can do this? > > Thank you > Martina > > >  > > > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20160908 10:49:17

Please send the description of the problem ;)) Wilfred On Thursday, September 8, 2016 7:30 PM, Martina Bekrová <to.zapomenu@...> wrote: Hello I would like to use GeoPDEs to solve Dirichlet problem for minimal surface equation using Newton method, just nonlinear poisson equation, but I'm not able to find out how to assemble the matrices for this variational problem. Description of the problem is attached below. I already have this program in Fenics, but I would like to apply IGA for this problem too. Is here a way how I can do this? Thank you Martina  _______________________________________________ Geopdesusers mailing list Geopdesusers@... https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Martina Bekrová <to.zapomenu@gm...>  20160908 10:22:45

Hello I would like to use GeoPDEs to solve Dirichlet problem for minimal surface equation using Newton method, just nonlinear poisson equation, but I'm not able to find out how to assemble the matrices for this variational problem. Description of the problem is attached below. I already have this program in Fenics, but I would like to apply IGA for this problem too. Is here a way how I can do this? Thank you Martina 
From: Windsor, Ben <B.W<indsor@wa...>  20160825 11:40:28

Dear Rafa Thank you so much for the detailed reply, that will definitely get me started in the C^0 and C^1 cases. The project will remain on curves it seems so that should keep things a bit simpler. Many thanks for the help, all the best Ben ________________________________ From: Rafael Vazquez <vazquez@...> Sent: 25 August 2016 12:05:53 To: geopdesusers@... Subject: Re: [Geopdesusers] Periodic NURBS space Dear Ben, periodic boundary conditions are not implemented yet, and I don't think I will implement them shortly. But all what you need should be already in GeoPDEs. Let me start with the simple C^0 periodic case and an open knot vector. In this case you only need to identify the first and last functions, changing the ID array. In Matlab I suggest you not to change the connectivity array, but to change the array for matrix construction. First, notice that the operators (op_gradu_gradv_tp, and so on) can give you as the output the rows, columns and values to construct a sparse matrix. You can impose periodicity in this way (only valid for curves): [rows, cols, vals] = op_gradu_gradv_tp (space, space, msh); rows(rows == space.ndof) = 1; cols(cols == space.ndof) = 1; mat = sparse (rows, cols, vals, space.ndof1, space.ndof1); You must also modify the righthand side: rhs = op_f_v_tp (space, msh, coeff); rhs(1) = rhs(1) + rhs(end); rhs(space.ndof) = []; And after you solve the system, to be able to use all the GeoPDEs functionality, you must also modify the computed solution: u = mat \ rhs; % A vector of length space.ndof1 u(space.ndof) = u(1); % A vector of length space.ndof, that you can now use for postprocessing That was the easy case. For the more difficult C^1 (or higher case), you need to do two things: 1. Unclamp the open knot vector, to get basis functions as in Fig. 1(b) of the article, with the correct weights. 2. Impose the periodicity in the matrix and righthand side, considering more than one index. For the first part, you can call the function nrbunclamp of the NURBS toolbox. This will give you the same parametrization with a nonopen knot vector, and with the correct weights. You can now use the new knot vector and weights to construct the space with sp_nurbs. Unfortunately, some of the functions of the toolbox (like nrbderiv) do not work with nonopen knot vectors, so you will have to call nrbunclamp *after geo_load*. In other words, the geometry is created with the open knot vector, but the space is defined with the nonopen knot vector. The second part is very similar to the C^0 case, and you can do it with a loop from 1 to k, where C^{k1} is the aimed continuity, identifying now the indices 1:k with space.ndofk+1:space.ndof. Some comments: 1. I have not studied the paper you cited, so I'm not completely sure that the approach is completely equivalent. 2. With this approach, you must have periodicity also in the parametrization, like in Fig. 17(a) of the paper. If you move one of the control points, without moving the corresponding control point on the other side, you would lose C^1 continuity in the physical domain. I'm not sure if this is the same in their approach. 3. If you are trying to solve a high order PDE, like CahnHilliard, you will need to avoid C^0 functions. Using nrbcirc, C^0 functions appear if you use an angle greater than pi/2. That's probably the reason why they used one sixth of the ring (Fig. 9). 4. I explained how to impose the periodic conditions for curves. For surfaces and volumes one can use the indices in boundary.dofs (C^0 case) and boundary.adjacent_dofs (C^1 case). 5. I think that periodic conditions may be already implemented in PetIGA (written in C). Best, Rafa On 24/08/2016 18:09, Windsor, Ben wrote: > Apologies, I actually meant the paper: > > 'Isogeometric analysis of the advective Cahn–Hilliard equation: Spinodal decomposition under shear flow' > > Sorry about that, > Ben > > > > ________________________________ > From: Windsor, Ben <B.Windsor@...> > Sent: 24 August 2016 16:38:12 > To: geopdesusers@... > Subject: [Geopdesusers] Periodic NURBS space > > Dear geoPDEs users, > > > I am a second year mathematics student at Warwick university on a summer research placement. For this project we have been using the geoPDEs library with MATLAB to solve a PDE on a closed curve (just first using the nrbcirc function to start with). However to do this it would be most useful to have a periodic function space to search instead of the normal nurbs space. > > > I was just wondering if this periodic basis space was something implemented in the libary? > > > I have had a look through the sp_nurbs and the @sp_scalar code and most of the library but cant see an option to add conditions onto the functions that make up the basis space or to somehow alter the underlying basis functions. I have seen such functions discussed in papers like 'Isogeometric Finite Element Data Structures based on Bezier Extraction of NURBS' by Michael J. Borden, Michael A. Scott, John A. Evans, and Thomas J.R. Hughes but have yet to see any code implementation of these periodic functions that is available, though apologies if I have overlooked some functionality in the library! > > > Many thanks, > > Ben >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers  _______________________________________________ Geopdesusers mailing list Geopdesusers@... https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Rafael Vazquez <vazquez@im...>  20160825 11:06:26

Dear Ben, periodic boundary conditions are not implemented yet, and I don't think I will implement them shortly. But all what you need should be already in GeoPDEs. Let me start with the simple C^0 periodic case and an open knot vector. In this case you only need to identify the first and last functions, changing the ID array. In Matlab I suggest you not to change the connectivity array, but to change the array for matrix construction. First, notice that the operators (op_gradu_gradv_tp, and so on) can give you as the output the rows, columns and values to construct a sparse matrix. You can impose periodicity in this way (only valid for curves): [rows, cols, vals] = op_gradu_gradv_tp (space, space, msh); rows(rows == space.ndof) = 1; cols(cols == space.ndof) = 1; mat = sparse (rows, cols, vals, space.ndof1, space.ndof1); You must also modify the righthand side: rhs = op_f_v_tp (space, msh, coeff); rhs(1) = rhs(1) + rhs(end); rhs(space.ndof) = []; And after you solve the system, to be able to use all the GeoPDEs functionality, you must also modify the computed solution: u = mat \ rhs; % A vector of length space.ndof1 u(space.ndof) = u(1); % A vector of length space.ndof, that you can now use for postprocessing That was the easy case. For the more difficult C^1 (or higher case), you need to do two things: 1. Unclamp the open knot vector, to get basis functions as in Fig. 1(b) of the article, with the correct weights. 2. Impose the periodicity in the matrix and righthand side, considering more than one index. For the first part, you can call the function nrbunclamp of the NURBS toolbox. This will give you the same parametrization with a nonopen knot vector, and with the correct weights. You can now use the new knot vector and weights to construct the space with sp_nurbs. Unfortunately, some of the functions of the toolbox (like nrbderiv) do not work with nonopen knot vectors, so you will have to call nrbunclamp *after geo_load*. In other words, the geometry is created with the open knot vector, but the space is defined with the nonopen knot vector. The second part is very similar to the C^0 case, and you can do it with a loop from 1 to k, where C^{k1} is the aimed continuity, identifying now the indices 1:k with space.ndofk+1:space.ndof. Some comments: 1. I have not studied the paper you cited, so I'm not completely sure that the approach is completely equivalent. 2. With this approach, you must have periodicity also in the parametrization, like in Fig. 17(a) of the paper. If you move one of the control points, without moving the corresponding control point on the other side, you would lose C^1 continuity in the physical domain. I'm not sure if this is the same in their approach. 3. If you are trying to solve a high order PDE, like CahnHilliard, you will need to avoid C^0 functions. Using nrbcirc, C^0 functions appear if you use an angle greater than pi/2. That's probably the reason why they used one sixth of the ring (Fig. 9). 4. I explained how to impose the periodic conditions for curves. For surfaces and volumes one can use the indices in boundary.dofs (C^0 case) and boundary.adjacent_dofs (C^1 case). 5. I think that periodic conditions may be already implemented in PetIGA (written in C). Best, Rafa On 24/08/2016 18:09, Windsor, Ben wrote: > Apologies, I actually meant the paper: > > 'Isogeometric analysis of the advective Cahn–Hilliard equation: Spinodal decomposition under shear flow' > > Sorry about that, > Ben > > > > ________________________________ > From: Windsor, Ben <B.Windsor@...> > Sent: 24 August 2016 16:38:12 > To: geopdesusers@... > Subject: [Geopdesusers] Periodic NURBS space > > Dear geoPDEs users, > > > I am a second year mathematics student at Warwick university on a summer research placement. For this project we have been using the geoPDEs library with MATLAB to solve a PDE on a closed curve (just first using the nrbcirc function to start with). However to do this it would be most useful to have a periodic function space to search instead of the normal nurbs space. > > > I was just wondering if this periodic basis space was something implemented in the libary? > > > I have had a look through the sp_nurbs and the @sp_scalar code and most of the library but cant see an option to add conditions onto the functions that make up the basis space or to somehow alter the underlying basis functions. I have seen such functions discussed in papers like 'Isogeometric Finite Element Data Structures based on Bezier Extraction of NURBS' by Michael J. Borden, Michael A. Scott, John A. Evans, and Thomas J.R. Hughes but have yet to see any code implementation of these periodic functions that is available, though apologies if I have overlooked some functionality in the library! > > > Many thanks, > > Ben >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers >  > _______________________________________________ > Geopdesusers mailing list > Geopdesusers@... > https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Windsor, Ben <B.W<indsor@wa...>  20160824 16:24:49

Apologies, I actually meant the paper: 'Isogeometric analysis of the advective Cahn–Hilliard equation: Spinodal decomposition under shear flow' Sorry about that, Ben ________________________________ From: Windsor, Ben <B.Windsor@...> Sent: 24 August 2016 16:38:12 To: geopdesusers@... Subject: [Geopdesusers] Periodic NURBS space Dear geoPDEs users, I am a second year mathematics student at Warwick university on a summer research placement. For this project we have been using the geoPDEs library with MATLAB to solve a PDE on a closed curve (just first using the nrbcirc function to start with). However to do this it would be most useful to have a periodic function space to search instead of the normal nurbs space. I was just wondering if this periodic basis space was something implemented in the libary? I have had a look through the sp_nurbs and the @sp_scalar code and most of the library but cant see an option to add conditions onto the functions that make up the basis space or to somehow alter the underlying basis functions. I have seen such functions discussed in papers like 'Isogeometric Finite Element Data Structures based on Bezier Extraction of NURBS' by Michael J. Borden, Michael A. Scott, John A. Evans, and Thomas J.R. Hughes but have yet to see any code implementation of these periodic functions that is available, though apologies if I have overlooked some functionality in the library! Many thanks, Ben  _______________________________________________ Geopdesusers mailing list Geopdesusers@... https://lists.sourceforge.net/lists/listinfo/geopdesusers 
From: Windsor, Ben <B.W<indsor@wa...>  20160824 15:53:56

Dear geoPDEs users, I am a second year mathematics student at Warwick university on a summer research placement. For this project we have been using the geoPDEs library with MATLAB to solve a PDE on a closed curve (just first using the nrbcirc function to start with). However to do this it would be most useful to have a periodic function space to search instead of the normal nurbs space. I was just wondering if this periodic basis space was something implemented in the libary? I have had a look through the sp_nurbs and the @sp_scalar code and most of the library but cant see an option to add conditions onto the functions that make up the basis space or to somehow alter the underlying basis functions. I have seen such functions discussed in papers like 'Isogeometric Finite Element Data Structures based on Bezier Extraction of NURBS' by Michael J. Borden, Michael A. Scott, John A. Evans, and Thomas J.R. Hughes but have yet to see any code implementation of these periodic functions that is available, though apologies if I have overlooked some functionality in the library! Many thanks, Ben 
From: Rafael Vázquez <vazquez@im...>  20160728 15:09:57

Hello Wilfred, for your information, mp_geo_load can also give as an output a list of subdomains/materials. You have to give this list in the text file of the geometry (see geo_specs_mp_v21.txt, in the doc folder). It can be a starting point. Best, Rafa 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20160727 12:30:04

I tested the code and it seems to work correctly. The next steps require a bit of software architecture. Indeed, the next question is how to set coefficients for each patch. One option could be to add a "material index" (integer) to each "geometry" structure. Then you loop over the patches, and geometry(iptc).mat_index would give the index where to find the coefficients for the patch. It is probably easiest when the user sets the geometry material index. In our case, in the end, we need different coefficients for each patch, and we need multigroup, i.e. we solve the equation for many values of the neutron energy. It will be necessary to make special routines. I will start to make some software architecture next week or so. The problem is to keep everything simple enough for the user (i.e. my student) but at the same time general enough to allow relevant calculations. OK, for now, I am happy with the multipatch thingie. I know how it works now, so the next step is to create the different "materials" for each patch. When that works, we will go to multigroup. Thanks, as always,Wilfred On Wednesday, July 27, 2016 7:59 PM, Rafael Vázquez <vazquez@...> wrote: Hello Wilfred, the code seems ok to me, but you (and your student) should be the ones to test it. :) A more general way to compute the coefficients, to allow the patchdependency and the computation from a previous solution, is already in my TODO list. But I would like to do that without adding a new function for each case. Best, Rafa 
From: Rafael Vázquez <vazquez@im...>  20160727 10:59:34

Hello Wilfred, the code seems ok to me, but you (and your student) should be the ones to test it. :) A more general way to compute the coefficients, to allow the patchdependency and the computation from a previous solution, is already in my TODO list. But I would like to do that without adding a new function for each case. Best, Rafa 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20160726 13:20:42

Hello Rafael, Thanks again. Would you be willing to check my code? In this code, the solution u of the previous iteration is used to calculate a new source term for the diffusion equation. Basically, there is a loop over all the patches; for each patch, there is a loop over the columns; for each column, the solution is calculated and used to determine the new source term: function out = calc_fission_source_mp(space, msh, u) out = zeros(space.ndof, 1); % % Loop over the patches in a multipatch domain % for iptc = 1:msh.npatch % % Loop over the columns of the mesh (slices in 3D); for each column, set % up the mesh (msh_col) and the basis functions (sp_col). % [eu, F] = sp_eval_msh(u, sp_col, msh_col), where eu is the evaluated % solution at the quadrature nodes and F are the locations of the % quadrature nodes in the actual domain (same as msh_col.geo_map but % not reshaped) % for iel = 1:msh.msh_patch{iptc}.nel_dir(1) msh_col = msh_evaluate_col(msh.msh_patch{iptc}, iel); sp_col = sp_evaluate_col(space.sp_patch{iptc}, msh_col); [eu, F] = sp_eval_msh(u(space.gnum{iptc}), sp_col, msh_col); for idim = 1:msh.rdim x{idim} = reshape (msh_col.geo_map(idim,:,:), msh_col.nqn, msh_col.nel); end nsf = nu_sigma_f(x{:}); out(space.gnum{iptc}) = out(space.gnum{iptc}) + ... op_f_v(sp_col, msh_col, nsf .* eu); end end end Obviously, in the future we will need to make an extra argument for the function nu_sigma_f(x, iptc) so that we can set different values for each patch. Thanks,Wilfred On Tuesday, July 26, 2016 6:27 PM, Rafael Vazquez <vazquez@...> wrote: Hello Wilfred, no, there is no equivalent of sp_eval for multipatch problems, but as you said one can use sp_eval for each patch. To get the correct degrees of freedom, the local space and the geometry, you can respectively use: u(space.gnum{patch}) space.sp_patch{patch} geometry(patch) You can find a similar example in the function sp_to_vtk, that already exists in the multipatch case. Best, Rafa 
From: Rafael Vazquez <vazquez@im...>  20160726 09:27:40

Hello Wilfred, no, there is no equivalent of sp_eval for multipatch problems, but as you said one can use sp_eval for each patch. To get the correct degrees of freedom, the local space and the geometry, you can respectively use: u(space.gnum{patch}) space.sp_patch{patch} geometry(patch) You can find a similar example in the function sp_to_vtk, that already exists in the multipatch case. Best, Rafa 
From: Wilfred van Rooijen <wvanrooijen@ya...>  20160726 09:13:37

Hello Rafael, Is there a multipatch equivalent of sp_eval()? For our work, we perform an iteration. In each iteration, the source is determined by the solution from the previous iteration. In earlier work we had only one patch, so that I used sp_eval() to get the solution on the quadrature points. For multipatch, I guess I could make a loop over the patches and then use sp_eval() on each patch but I am a bit hesitant: how to get the correct part of the solution u, and the correct info for space and geometry on each patch from the "overall" info? Thanks,Wilfred 
From: Rafael Vazquez <vazquez@im...>  20160725 14:54:56

Hi Nitel, they model the shells as 3D solids with very few elements in the throughthickness direction, and solve a linear elasticity problem. If you want to reproduce exactly their results, everything should be ready in GeoPDEs except the imposition of the force in a single point (as in the pinched hemisphere and the pinched cylinder), but that should not be difficult. Of course, for the geometry you should give the 3D solid, not just the surface, but these examples should be easy to generate with some functions in the NURBS toolbox (nrbruled, nrbcirc, nrbrevolve, nrbextrude). Instead, if you want a surface shell model, it has not been implemented yet. Best, Rafa 