On Fri, Aug 25, 2017 at 8:08 PM, Mike Marchywka <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>
> 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
> -us...@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
