Hi Mike, Thanks once again. I will do the calculations.
The problem is not in the normals, at least not only. I am pretty sure the gradient at the quadrature point is wrong. It must be around 1E+5. Must not cancel out. The values in the dofs are reasonable. Shape function issue. Must be. I dont believe in any non-physical stuff in this case. It is such a simple case! Will let you know about my findings! Thanks Renato On Sat, Aug 26, 2017 at 10:43 AM, Mike Marchywka <marchy...@hotmail.com> wrote: > > 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 <rebp...@gmail.com> > 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 <marchy...@hotmail.com > <mailto:marchy...@hotmail.com>> 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 <rebp...@gmail.com<mailto:rebp...@gmail.com>> > 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 <marchy...@hotmail.com > <mailto:marchy...@hotmail.com><mailto:marchy...@hotmail.com<mailto:m > archy...@hotmail.com>>> 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 <rebp...@gmail.com<mailto:rebp...@gmail.com><mailto: > rebp...@gmail.com<mailto:rebp...@gmail.com>>> > 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 <rebp...@gmail.com<mailto: > rebp...@gmail.com><mailto:rebp...@gmail.com<mailto:rebp...@gmail.com>>> > 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 <jwpeter...@gmail.com > <mailto:jwpeter...@gmail.com><mailto:jwpeter...@gmail.com<mailto:jwpe > ter...@gmail.com>>> > > wrote: > > > >> > >> > >> On Fri, Aug 25, 2017 at 10:34 AM, Renato Poli <rebp...@gmail.com > <mailto:rebp...@gmail.com><mailto:rebp...@gmail.com<mailto:rebpoli@ > gmail.com>>> 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 > Libmesh-users@lists.sourceforge.net<mailto:Libmesh > -us...@lists.sourceforge.net><mailto:Libmesh-users@lists.sourceforge.net > <mailto:Libmesh-users@lists.sourceforge.net>> > 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 Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
