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 ?
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?
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.
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
