Why fe.shape_value( const unsigned int i, const Point< dim > & p ) return a negative number?

[Zhao Yidong]

Hi everyone!

I use fe, which is FESystem class, but I find fe.shape_value( const unsigned int i, const Point< dim > & p ) (the document is here: https://www.dealii.org/8.5.0/doxygen/deal.II/classFE__Poly.html#a8c769a449d25d54e756fc67cca5dd30b) sometimes returns a negative number which seems not right.

The Point I send into the function locates in the reference space ( [-1, 1]x[-1, 1] ), I also try the Point in the real space, but also get negative numbers sometimes.

Thanks any response!

[Denis Davydov]

Hi,

it's perfectly fine for a general FE shape functions to be negative at some points, that's clearly the case for quadratic

http://hplgit.github.io/INF5620/doc/pub/sphinx-fem/._main_fem003.html#fem-approx-fe-fig-p2 

Denis.

[Zhao Yidong]

Thanks a lot, it's really clear to understand by the help of the website you share!

Indeed, my problem is this:

1. I have fe which is a FiniteElement class, and I have fe_values which is a FEValues class. I also have QGauss<dim> quadrature to create fe_values.

2. I iterate on every quadrature point in a iteration on every cell, store the quadrature points.

3. Use fe_values.value(i, q) to calculate shape functions, and fe.shape_value(i, point) to calculate another shape function, but find they are different.

So I want to know what's the difference between this two methods? Because I want to separate the fe and quadrature points, so can I use fe.shape_value to calculate the shape functions?

Thanks again!

在 2018年5月6日星期日 UTC+8上午6:08:55，Denis Davydov写道：

[Wolfgang Bangerth]

On 05/05/2018 09:11 AM, Zhao Yidong wrote:
>
> The Point I send into the function locates in the reference space ( [-1,
> 1]x[-1, 1] ),

By the way, in deal.II the reference cell is [0,1]x[0,1]. So even bilinear
shape functions will be negative if you evaluate them at (-1,-1), for example.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Wolfgang Bangerth]

On 05/06/2018 09:05 AM, Zhao Yidong wrote:
>
> Indeed, my problem is this:
> 1. I have fe which is a FiniteElement class, and I have fe_values which is a
> FEValues class. I also have QGauss<dim> quadrature to create fe_values.
> 2. I iterate on every quadrature point in a iteration on every cell, store the
> quadrature points.
> 3. Use fe_values.value(i, q) to calculate shape functions, and
> fe.shape_value(i, point) to calculate another shape function, but find they
> are different.
>
> So I want to know what's the difference between this two methods? Because I
> want to separate the fe and quadrature points, so can I use fe.shape_value to
> calculate the shape functions?

FiniteElement defines shape functions on the reference element. FEValues
describes these shape functions on the real cell. For many elements, you get
the same result whether you
* evaluate the shape function on the reference cell at a quadrature point
also defined on the reference cell (via FiniteElement::shape_value and
Quadrature::point)
* evaluate the shape function on the real cell at a quadrature point
also defined on the real cell (via FEValues::shape_value, with quadrature
point locations also available from FEValues).

But there are also elements for which this is not true. It's not clear to me
what exactly you do in your comparison -- can you show a piece of code?

(You will have to use FEValues anyway at one point if you also want the
*gradients* of the shape functions -- these are never the same on reference
and real cell.)

[Zhao Yidong]

Thank you very much!

My fe is FE_Q(1) type:

FE_Q<dim>(1)

I make comparison like this:

1.Create fe_values like this:

FEValues<dim> fe_values(fe, quadrature,

                               update_values |

                               update_quadrature_points |

                               update_gradients |
                               update_JxW_values);

2.Iterate on every quadrature point in every cell, store them. (These points are in real space right?)

for (unsigned q=0; q<quadrature.size(); ++q) {

               Point<dim> point = fe_values.quadrature_point(q);

3. Use fe.shape_value() to get shape value at every degree of a cell, store_point is the point I store in the 2nd step.

double shapeValue = fe.shape_value(i, store_point); // shape value

4. Compare with shape value calculated by fe_values.

double N = fe_values.value(i, q);

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Because above gives me different answers, so I think maybe the points I store are in real space, so I transfer them into reference space( here [-1, 1]x[-1, 1] ) as this:

1. With the help of TriaAccessor, I can get cell's every node.

// there are 8 nodes every cell
for (unsigned numberPoint=0; numberPoint<8; ++numberPoint) {

            vertexPoint = triaAccessorP->vertex(numberPoint);

Use these nodes, I transfer the real point into a reference point.

2. I use these points in reference space to calculate shape function:

double shapeValue = fe.shape_value(i, reference_point); // shape value

This also gives me different answer with fe_values

I will check can I get the same answer if I transfer the real point into reference point(in [0,1]x[0,1]).

(My understanding of gradient is differential, so can I calculate it by finite difference?

I can get forward_point = current_point + Point<dim>(epsilon, 0, 0)

and backward_point = current_point - Point<dim>(-epsilon, 0, 0)

to calculate gradient[0] = (shape_function(forward_point) - shape_function(backward_point)) / 2*epsilon,

and the same way for gradient[1], gradient[2] in 3d space ?)

Thanks again!

Best

Yidong

在 2018年5月7日星期一 UTC+8下午1:20:34，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 05/07/2018 02:20 PM, Zhao Yidong wrote:
>
> 2.Iterate on every quadrature point in every cell, store them. (These points
> are in real space right?)
> |
>

> for(unsignedq=0;q<quadrature.size();++q){

>
> Point<dim>point =fe_values.quadrature_point(q);
>
> |
>
> 3. Use fe.shape_value() to get shape value at every degree of a cell,
> store_point is the point I store in the 2nd step.
> |
>
> doubleshapeValue =fe.shape_value(i,store_point);// shape value
>
> |

This is your problem. 'store_point' lives in real space, but 'fe.shape_value'
wants a point on the reference cell. You need to call
fe.shape_value (quadrature.point(q));

[Zhao Yidong]

Thanks a lot! I get the same answer by using quadrature.point(q), and I get the similar answer after I transfer the real space point into reference space( [0, 1]x[0, 1] ), and then use fe.shape_value(reference_point).

Best

Yidong

在 2018年5月7日星期一 UTC+8下午2:37:22，Wolfgang Bangerth写道：