From: alexis papagiannopoulos <alexis.papagiannopoulos@gm...>  20111012 11:49:14

Dear geopdesusers, I tried to test some convergenece rates by applying several refinement strategies such krefinement and h refinement.Although h refinement has relatively a smooth convergence , krefinement presents a lot of oscillations .After creating a krefinement m file as such: function geo = krefinement( geo,degr,knts) nurbs=geo.nurbs; degelev=max(degr(nurbs.order1),0); nurbs=nrbdegelev(nurbs,degelev); [rknots,zeta,nknots]=kntrefine(nurbs.knots,knts,nurbs.order1,nurbs.order2); nurbs=nrbkntins(nurbs,nknots); geo=geo_load(nurbs); end I then apply homo_kerror.m to the 'geo_ring.txt' domain function errors = homo_kerror(geometry,i ) errors=[]; for j=1:1:i [dofs(1) errors(1)] = homo_poisson( geometry ); refined_geo(j)=krefinement(geometry,[j j],[1 1]); [dofs(j+1) errors(j+1)]=homo_poisson(refined_geo(j)); [dofs3(1) errors3(1)] = homo_poisson( geometry ); refined_geo2(j)=hrefinement(geometry,[j j]); [dofs2(j+1) errors2(j+1)]=homo_poisson(refined_geo2(j)); end loglog(dofs,errors,'ko',dofs2,errors2,'go') legend('krefinement','hrefinement') xlabel('dofs'); ylabel('error'); grid on end where homo_poisson.m is the problem presented in the geopdes report (article 4.1) function [dofs ,err] = homo_poisson( geometry ) nurbs=geometry.nurbs; geometry=geo_load(nurbs); knots = geometry.nurbs.knots; [qn, qw] = msh_set_quad_nodes (knots, msh_gauss_nodes (geometry.nurbs.order)); msh = msh_2d (knots, qn, qw, geometry); space = sp_nurbs_2d (geometry.nurbs, msh); dofs=space.ndof; mat = op_gradu_gradv_tp (space, space, msh, @(x, y) ones (size (x))); rhs = op_f_v_tp (space, msh, @(x, y) (89*sqrt(x.^2+y.^2)).*sin(2*atan(y./x))./(x.^2+y.^2)); drchlt_dofs = []; for iside = 1:4 drchlt_dofs = union (drchlt_dofs, space.boundary(iside).dofs); end int_dofs = setdiff (1:space.ndof, drchlt_dofs); u = zeros (space.ndof, 1); u(int_dofs) = mat(int_dofs, int_dofs) \ rhs(int_dofs); sp_to_vtk (u, space, geometry, [20 20], 'laplace_solution.vts', 'u') err = sp_l2_error (space, msh, u, @(x,y)(x.^2+y.^23*sqrt(x.^2+y.^2)+2).*sin(2.*atan(y./x))) end I tried several domains and problems and I get similar results. In krefinement the order elevates by 1 and then a new knot is inserted(in both directions) at each step. any suggestions why this occurs? Thank you, Alexis Papagiannopoulos 
From: Rafael Vazquez <vazquez@im...>  20111013 08:28:30

Hi Alexis, I've tried to do the same modifying the file ex_laplace_iso_ring (see below). The oscillations appear around degree 17, that is much higher than what you usually do in IGA. The first thing I must say is that this shouldn't be considered krefinement: you have only one element, and the continuity is not playing any role. What you are doing is some kind of spectral method with NURBS. You can try to read some works on spectral methods (I'm not an expert) and check if they need some special technique to obtain the right convergence, maybe a particular quadrature rule. But I'm also thinking that the error may come from the NURBS toolbox. I know that some functions in the toolbox will add rounding errors to the computation of the shape functions, and I guess that these errors become much higher for high degrees. Since I've always worked with degree lower than 8, I never had any problem with this. You can check if the problem is coming from there. Regards, Rafa The problem_data is taken from ex_laplace_iso_ring. The postprocessing part is removed. The method data and the call to the solver are like this method_data.nsub = [1 1]; % In this example, this is only one element! for ii = 2:20 method_data.degree = [ii ii]; method_data.regularity = [ii1 ii1]; method_data.nquad = [ii+1 ii+1]; [geometry, msh, space, u] = solve_laplace_2d_iso (problem_data, method_data); [error_h1(ii), error_l2(ii)] = ... sp_h1_error (space, msh, u, problem_data.uex, problem_data.graduex); ndofs(ii) = space.ndof; end loglog(ndofs, error_l2); 
From: Rafael Vazquez <vazquez@im...>  20111028 10:39:24

Dear GeoPDEs users, for some days I kept thinking about the krefinement. I think this is an interesting topic, so I send you a summary of what I did. First of all, Alexis was right about the krefinement. I misunderstood his code, and thought that he only had one element, but in fact he was using a 2x2 mesh. So the functions are not globally smooth, as I thought, and the continuity is important. Well, I asked a colleague to repeat the same test in his own code (thanks, Massimiliano). The results are the same, except for rounding errors (that grow with the degree), and the oscillations appear from a very low degree. At least it is clear that the problem is not a bug in the code. We have also tried increasing the number of elements, and the oscillations appear much later, or don't appear at all. The explanation comes from the theory. When doing h or prefinement, the error decreases at each refinement step, because the spaces are nested. But when doing krefinement, since you increase the degree and the continuity at the same time, the spaces are not nested, and nothing assures you that the error will be reduced. Numerically, krefinement seems to work better when the mesh is fine enough, and in this particular example, for coarse meshes, odd degree discretizations seem to work better than the even degree ones. Regards, Rafa 