It should not be too hard to just write the value u = s1*u1+s2*u2 ... and differentiate s_i to get the grad at that q point and see if the code matches the paper and look for any obvious missing factors or terms. I'm just doing QUAD4 - linear shape functions and rectangles with a lot of my own code- but any polynomial not much harder. It looks like you have a sum of values that nearly cancel out. I know I did some interpolation and the linear worked better than quadratic due to nonphysical excess curvature ( often leading to nonphysical negative values between nodes, In some cases I thought interpolating the logs would work better but never got back to it) . But with an analytical results again comparisons should be informative. Also I guess with an integral over a closed surface you probably need to check the normal sign but that should be obvious just plotting the flux at each node if you expect it to be uniform rather than oscillating :)
________________________________________ From: Renato Poli <[email protected]> Sent: Friday, August 25, 2017 7:46 PM To: Mike Marchywka Cc: John Peterson; libmesh-users Subject: Re: [Libmesh-users] Trouble with boundary integration On Fri, Aug 25, 2017 at 8:08 PM, Mike Marchywka <[email protected]<mailto:[email protected]>> wrote: >Any hint on how to easily reproduce the shape function calculations by hand >(that is, manually producing the dphi_face values)? Very carefully :) I ended up writing some polynomial classes and doing the integrals and derivatives for QUAD4 that way. With my rational number class it was easy to get exact results and make generated code for the simple things like Laplacian of rectangle. That seems a good idea. I am using TRI6. The results I am getting are so off the expected that I won't matter getting approximations for now as reference. This is the reason why I think I got a conceptual issue. It is really far from target. I guess you could probably do the same with MatLab or Mathematica or maybe R. I saw the point by point quadrature thing in the examples but AFAICT in many cases you can do it exactly and just skip a lot of that. The Laplacian is just an integral of a second derivative which would be zero for linear shape functions except for integration by parts which then gives you an easily computable integral. I'm using all analytical integrals now- just putting numbers into the system matrix and rhs with QUAD4 elements- but have a possible model problem that probably needs point by point integration :) But do the derivatives look right? The normals? The solution? Does it work without AMR ? p? h? etc... The normals look right. I have simplified the mesh to have right angles and thus simple normals (1,0,0) or (-1,0,0) etc. I am not using AMR at all. Mesh is dynamically generated with a couple of parameters (maximum triangle size and internal angles, for example). What does work, if anything? Everything works in a wrong way. :-) Look this example: +++ Triangle coordinates (each dof) [0] (x,y,z)=( 10, 5, 0) [1] (x,y,z)=( 10, 10, 0) [2] (x,y,z)=( 6.73058, 6.91831, 0) [3] (x,y,z)=( 10, 7.5, 0) [4] (x,y,z)=( 8.36529, 8.45915, 0) [5] (x,y,z)=( 8.36529, 5.95915, 0) See the highlighted system solutions for dof[1] and dof[2]. It shows a delta of 0.1206E7 in (dx=10-6.73=3.27). Something like grad_x=0.0369E7. I dont get why the shape function and the summation sort of kills this gradient... Still struggling here... +++ intermediate calculations qface_point = (x,y,z)=( 10, 7.5, 0) == i=0 === dof_indices[0]: 12 dphi_face[0][0]: (x,y,z)=(0.188516, -0.2, -0) _sys_solution[0->12]: 3e+07 _sys_solution[0->12] * dphi_face[0]: (x,y,z)=(5.65548e+06, -6e+06, -0) == i=1 === dof_indices[1]: 13 dphi_face[1][0]: (x,y,z)=(0.117348, 0.2, 0) _sys_solution[1->13]: 3e+07 /// DOF[1] => 3E7 @ (10,5) /// <=== _sys_solution[1->13] * dphi_face[1]: (x,y,z)=(3.52045e+06, 6e+06, 0) == i=2 === dof_indices[2]: 14 dphi_face[2][0]: (x,y,z)=(0.305864, -0, -0) _sys_solution[2->14]: 3.12065e+07 /// DOF[2] => 3.12E7 @ (6.73, 6.91) /// <=== _sys_solution[2->14] * dphi_face[2]: (x,y,z)=(9.54495e+06, -0, -0) == i=3 === dof_indices[3]: 15 dphi_face[3][0]: (x,y,z)=(0.611729, 0, 0) _sys_solution[3->15]: 3e+07 _sys_solution[3->15] * dphi_face[3]: (x,y,z)=(1.83519e+07, 0, 0) == i=4 === dof_indices[4]: 16 dphi_face[4][0]: (x,y,z)=(-0.611729, 0, 0) _sys_solution[4->16]: 3.02508e+07 _sys_solution[4->16] * dphi_face[4]: (x,y,z)=(-1.85053e+07, 0, 0) == i=5 === dof_indices[5]: 17 dphi_face[5][0]: (x,y,z)=(-0.611729, 0, 0) _sys_solution[5->17]: 3.03524e+07 _sys_solution[5->17] * dphi_face[5]: (x,y,z)=(-1.85674e+07, 0, 0) ======= JxW_face[0]: 5 grad_at_qn= (x,y,z)=(-6.29947e-06, -7.7514e-06, 0) normals[0]= (x,y,z)=( 1, 0, -0) apply_normal= -3.14973e-05 ________________________________________ From: Renato Poli <[email protected]<mailto:[email protected]>> Sent: Friday, August 25, 2017 6:12 PM To: Mike Marchywka Cc: John Peterson; libmesh-users Subject: Re: [Libmesh-users] Trouble with boundary integration Hi Mike, Thanks for the reply. Any hint on how to easily reproduce the shape function calculations by hand (that is, manually producing the dphi_face values)? On Fri, Aug 25, 2017 at 6:17 PM, Mike Marchywka <[email protected]<mailto:[email protected]><mailto:[email protected]<mailto:[email protected]>>> wrote: I've been trying to learn this myself and encounter a variety of issues like this. If you have an exact solution, why not just print all the pieces- values and normals etc- and see how they differ point by point instead of trying to do an inverse problem of guessing what is wrong with an integral ? Well - that is what I am trying to do. I isolated a single element and, as the problem is symmetrical, it should produce a non-zero result. There is not a lot of freedom in the calculation. It should be a simple multiplication of a bunch of numbers... Varying the mesh size seems to be useful sometimes. How do the element sizes compare to the radius of curvature ? Can you test the code on a simpler geometry such as between two planes? Agreed. I spent two days on that, varying all sort of stuff. I tried from very small elements, to huge ones. I think my problem might be related to the order of the shape function or the gauss quadrature. Really not sure. Any hint on how to easily reproduce the shape function calculations by hand (that is, manually producing the dphi_face values)? AFAICT in the general case doing integrals of normal components over surfaces is not always easy with p and h refinements etc. Although maybe someone can clarify a bit. I am using CGAL to generate and refine the mesh. It is quite flexible. I am really considering that I have a conceptual issue. It seems there is something really close to me that is fooling me ... I will keep trying. Please let me know if anyone has any suggestion ... Thanks, Renato Thanks. ________________________________________ From: Renato Poli <[email protected]<mailto:[email protected]><mailto:[email protected]<mailto:[email protected]>>> Sent: Friday, August 25, 2017 3:02 PM To: John Peterson Cc: libmesh-users Subject: Re: [Libmesh-users] Trouble with boundary integration Hi John, I am still struggling and running out of ideas here. I wrote down some debugging - perhaps you can find something weird. The calculations seem correct, at least to what I could understand so far. But things are zero'ing out. If, in the mesh creation, I declare the variable as "FIRST" order, things do not zero out. The integration of the gradient only zeros out if I use: sys.add_variable("u", SECOND); For face integration, I am using first order gauss quadrature: QGauss qface(dim-1, FIRST); Any hint? Any debugging suggestion? Rgds, Renato # Debug[1]: Elem: 4 # Debug[1]: (x,y,z)=( 5, 10, 0) # Debug[1]: (x,y,z)=( 0, 10, 0) # Debug[1]: (x,y,z)=( 3.08169, 6.73058, 0) # Debug[1]: (x,y,z)=( 2.5, 10, 0) # Debug[1]: (x,y,z)=( 1.54085, 8.36529, 0) # Debug[1]: (x,y,z)=( 4.04085, 8.36529, 0) # Debug[1]: num_dofs=6 # Debug[1]: Side: 0 # Debug[1]: Boundary: 1 # Debug[1]: qface_point = (x,y,z)=( 2.5, 10, 0) # Debug[1]: == i=0 === # Debug[1]: dof_indices[0]: 23 # Debug[1]: dphi_face[0][0]: (x,y,z)=( 0.2, 0.188516, -0) # Debug[1]: _sys_solution[0->23]: 3e+07 # Debug[1]: _sys_solution[0->23] * dphi_face[0]: (x,y,z)=( 6e+06, 5.65548e+06, -0) # Debug[1]: == i=1 === # Debug[1]: dof_indices[1]: 24 # Debug[1]: dphi_face[1][0]: (x,y,z)=( -0.2, 0.117348, 0) # Debug[1]: _sys_solution[1->24]: 3e+07 # Debug[1]: _sys_solution[1->24] * dphi_face[1]: (x,y,z)=( -6e+06, 3.52045e+06, 0) # Debug[1]: == i=2 === # Debug[1]: dof_indices[2]: 25 # Debug[1]: dphi_face[2][0]: (x,y,z)=( 0, 0.305864, -0) # Debug[1]: _sys_solution[2->25]: 3.12065e+07 # Debug[1]: _sys_solution[2->25] * dphi_face[2]: (x,y,z)=( 0, 9.54495e+06, -0) # Debug[1]: == i=3 === # Debug[1]: dof_indices[3]: 26 # Debug[1]: dphi_face[3][0]: (x,y,z)=( -0, 0.611729, 0) # Debug[1]: _sys_solution[3->26]: 3e+07 # Debug[1]: _sys_solution[3->26] * dphi_face[3]: (x,y,z)=( -0, 1.83519e+07, 0) # Debug[1]: == i=4 === # Debug[1]: dof_indices[4]: 27 # Debug[1]: dphi_face[4][0]: (x,y,z)=( 0, -0.611729, 0) # Debug[1]: _sys_solution[4->27]: 3.02508e+07 # Debug[1]: _sys_solution[4->27] * dphi_face[4]: (x,y,z)=( 0, -1.85053e+07, 0) # Debug[1]: == i=5 === # Debug[1]: dof_indices[5]: 28 # Debug[1]: dphi_face[5][0]: (x,y,z)=( 0, -0.611729, 0) # Debug[1]: _sys_solution[5->28]: 3.03524e+07 # Debug[1]: _sys_solution[5->28] * dphi_face[5]: (x,y,z)=( 0, -1.85674e+07, 0) # Debug[1]: ======= # Debug[1]: JxW_face[0]: 5 # Debug[1]: grad_at_qn= (x,y,z)=(8.95187e-06, -6.02752e-06, 0) # Debug[1]: normals[0]= (x,y,z)=( 0, 1, -0) # Debug[1]: JxW_face[0] * grad_at_qn * normals[0] = -3.01376e-05 On Fri, Aug 25, 2017 at 2:25 PM, Renato Poli <[email protected]<mailto:[email protected]><mailto:[email protected]<mailto:[email protected]>>> wrote: > Hi John, > > _sys_solution is the system->solution vector. > > It is fed by: > _system.solution->localize_to_one( _sys_solution ); > > > On Fri, Aug 25, 2017 at 1:48 PM, John Peterson > <[email protected]<mailto:[email protected]><mailto:[email protected]<mailto:[email protected]>>> > wrote: > >> >> >> On Fri, Aug 25, 2017 at 10:34 AM, Renato Poli >> <[email protected]<mailto:[email protected]><mailto:[email protected]<mailto:[email protected]>>> >> wrote: >> >>> Hi, >>> >>> I need some help here ... >>> >>> I am having trouble to understand what is going on with a boundary >>> integral. >>> I have solved a Poisson problem for pressure in a 2D mesh, as a model >>> for a >>> incompressible water flow in porous media. >>> The mesh (Tri6 elements) has a near-circular polygon in the middle >>> representing an injection well. >>> >>> RHS of the equations are zero except in boundaries dofs. Dirichlet >>> boundary >>> conditions are imposed using a penalty both on the well and on the >>> external >>> boundaries to represent fixed pressure. (the solution looks ok - seems to >>> work) >>> >>> To validate the solution, I am trying to estimate the water flow both in >>> the well and in the external boundaries (they must be identical and match >>> the analytical solution). I therefore calculate the normal pressure >>> gradient and integrate in the whole boundary path. >>> >>> When I use first order variable, the results are "almost correct". But >>> when >>> I use second order , it is absolutely off. >>> >>> I am using the code below to do the integration. I tried many >>> alternatives, >>> no success. >>> >>> Any idea of what I am doing wrong? >>> >> >> The code seems OK to me, although you should not need to call >> fe_face->reinit() inside the loop over boundary ids. >> >> What is _sys_solution? >> >> -- >> John >> > > ------------------------------------------------------------------------------ Check out the vibrant tech community on one of the world's most engaging tech sites, Slashdot.org! http://sdm.link/slashdot _______________________________________________ Libmesh-users mailing list [email protected]<mailto:[email protected]><mailto:[email protected]<mailto:[email protected]>> https://lists.sourceforge.net/lists/listinfo/libmesh-users ------------------------------------------------------------------------------ Check out the vibrant tech community on one of the world's most engaging tech sites, Slashdot.org! http://sdm.link/slashdot _______________________________________________ Libmesh-users mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-users
