Hi Mike,
Thanks for the reply.
On Fri, Aug 25, 2017 at 6:17 PM, Mike Marchywka <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>
> 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> 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>
> > wrote:
> >
> >>
> >>
> >> On Fri, Aug 25, 2017 at 10:34 AM, Renato Poli <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
> 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
