>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. 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... What does work, if anything? ________________________________________ From: Renato Poli <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>> 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>> 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>> 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>> > wrote: > >> >> >> On Fri, Aug 25, 2017 at 10:34 AM, Renato Poli >> <rebp...@gmail.com<mailto:rebp...@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-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
