evaluation function and its gradient after rotating the triangulation

[Mathieu]

Hello everyone,

I am working on a 2d triangulation which has a rectangular domain (see tria.png) perpendicular to the x- and y-axis (see tria.png).

I rotated this triangulation using GridTools::rotate(alpha, tria) and shifted it so as that the bottom left corner after these transformations has the original coordinates (see tria_rotated.png). 

The rest of what I do with the triangulation is identical.

In both screenshots, the function

f(x,y) over the respective triangulations is plotted. 

In my case, the solution vector associated with f(x,y) is obtained from

VectorTools::interpolate(___).

I first wanted to check that f(x,y) is identical at arbitrary points contained in both triangulations, however, there are small differences. 

As en example, see the point with the white cross -- f(x,y) at the rotated triangulation is 2.81005 and f(x,y) at the original triangulation is 2.8047. 

f(x,y) at the support points is correct and matches the value returned from the function object (Function::value) passed to the interpolate call. 

I am surprised about the differences since I expected that rotating the triangulation has no effect on f(x,y), if evaluated at the same point. 

My second question pertains to the gradient, which is a Tensor<1,2>.

I need to compute the derivatives of f(x,y) in x-and y-direction. 

Say, I call get_function_gradients(res,f(x,y)) based on the rotated triangulation, 

stores 'res' the derivatives wrt x and y, 

or do I have to make a transformation?

I expected to not make a transformation, since all what I did 

is to rotate the triangulation, but not the x-y coordinate system.

However, what get_function_gradients returns after the transformation 

is clearly wrong. 

df(x,y) / dy = 0 per design, which is correct in the original triangulation, 

but not in the rotated triangulation. 

Do I miss something here?

Thank you,

Math

[Wolfgang Bangerth]

> In both screenshots, the function
> f(x,y) over the respective triangulations is plotted.

No, it isn't. What you are plotting is the *interpolated* version of the
function. As a consequence...

> As en example, see the point with the white cross -- f(x,y) at the rotated
> triangulation is 2.81005 and f(x,y) at the original triangulation is 2.8047.
> f(x,y) at the support points is correct and matches the value returned from
> the function object (Function::value) passed to the interpolate call.
> I am surprised about the differences since I expected that rotating the
> triangulation has no effect on f(x,y), if evaluated at the same point.

...this should not be a surprise unless you can convince that the
*interpolated* function should be the same in both cases. This is a question
about the function f(x,y) you use: Is it in the function space you are using
or not?

> My second question pertains to the gradient, which is a Tensor<1,2>.
> I need to compute the derivatives of f(x,y) in x-and y-direction.
> Say, I call get_function_gradients(res,f(x,y)) based on the rotated
> triangulation,
> stores 'res' the derivatives wrt x and y,
> or do I have to make a transformation?

It returns the derivatives with regards to the global coordinate system within
which the triangulation lives, not the reference coordinate system of the
cells (or the coordinate system in which the triangulation lived before it was
rotated -- the triangulation does not retain a memory of its rotation).

> df(x,y) / dy = 0 per design, which is correct in the original triangulation,
> but not in the rotated triangulation.
> Do I miss something here?

Like before, what you need to think about is whether the *interpolated*
version of f(x,y) should have a zero derivative.

Best
W.

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

[Mathieu]

" No, it isn't. What you are plotting is the *interpolated* version of the

function. As a consequence...

...this should not be a surprise unless you can convince that the
*interpolated* function should be the same in both cases. This is a question
about the function f(x,y) you use: Is it in the function space you are using
or not? "

Makes totally sense!

The function I interpolate using an FE_Q<1> element is

f(x,y) = 0.12x + 0.002x^4 + sqrt(y), 3.0<{x,y}<4.0

Lets call the FE function f^h(x,y) which lives in the space H^1,

hence it lives in a different space than f(x,y).

Clearly, df^h(x,y) / dx and df^h(x,y) / dy are discontinuous and

I project them using global smoothing as implemented in VectorTools::project.

, hence df^h(x,y) / dx and df^h(x,y) / dy live again in H^1 and I can

differentiate them again to calculate the "second derivatives" of f^h(x,y).

The mixed second partial derivatives of f(x,y) are zero per design,

which is perfectly captured by f^h(x,y) in the unrotated grid with

only 9 degrees of freedom (refining once globally a GridGenerator::general_cell).

After rotating the grid as shown in my question, the mixed second partial derivatives

of f^h(x,y) are no longer zero,

although f^h(x,y) and the smoothed first derivatives show good agreement with the analytical counterparts.

I think this is because the (lines of the) cells are not aligned with the global coordinate system and

things get mixed up due to the interpolation.

If the cells are aligned with the global coordinate system,

the decoupling of x and y can be captured on a coarse grid

Does this argumentation make sense?

I would have to refine the rotated grid roughly seven times to approximate this decoupling accurately.

Unfortunately, I can refine the grid globally at most twice.

Do you have an idea how to capture the above decoupling also at the rotated grid better --

maybe by generating the grid differently?

The input to my grid are four points which represent the minimal bounding box of a set of points in 2d

and the grid should roughly represent this box.

Thank you very much for your input,

Math

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/34518363-1cb9-2912-5a5e-1e6a5fcdb00c%40colostate.edu.

[Wolfgang Bangerth]

On 6/9/23 09:14, Mathieu wrote:
>
> I would have to refine the rotated grid roughly seven times to approximate
> this decoupling accurately.
> Unfortunately, I can refine the grid globally at most twice.
> Do you have an idea how to capture the above decoupling also at the rotated
> grid better --
> maybe by generating the grid differently?
> The input to my grid are four points which represent the minimal bounding box
> of a set of points in 2d
> and the grid should roughly represent this box.

I don't actually understand why you would want to do that. If I understand
correctly, then you have a function
f(x,y) = g(x) + h(x)
which you interpolate or project onto a grid without very specific properties
to get
f^h(x,y)
and you observe that f^h is no longer a sum of functions of x and y separately.

What isn't clear to me why that bothers you. If you don't like it, just
project g and h to the mesh separately, or even better project them onto 1d
meshes. If you have a situation where you are dealing with functions of one
argument, treat them as functions of one argument.

[Mathieu]

"I don't actually understand why you would want to do that. If I understand
correctly, then you have a function
   f(x,y) = g(x) + h(x)
which you interpolate or project onto a grid without very specific properties
to get
   f^h(x,y)
and you observe that f^h is no longer a sum of functions of x and y separately.

What isn't clear to me why that bothers you. If you don't like it, just
project g and h to the mesh separately, or even better project them onto 1d
meshes. If you have a situation where you are dealing with functions of one
argument, treat them as functions of one argument."

f(x,y) = g(x)+h(y) 

is just for this example that I showed you. 

The function f is in general dependent on x and y, and in most versions of f, they are not separated.  

The fact that f^h is no longer separated changes the solution of my pde drastically:

If I run my program with the analytical f or with f^h on the axis parallel grid, the output is quite similar. This is not so if I run it with f^h defined on the rotated grid. 

And since the axis parallel grid is able to capture the separation, 

I hoped the rotated grid can capture it too. 

But if I understood you correctly, there is no obvious workaround for this issue , right? 

Best,

Math

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---

You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/QopdaFOjy04/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a29f5219-fea7-9cbf-efd8-07aa48ffed49%40colostate.edu.

[Wolfgang Bangerth]

On 6/9/23 12:52, Mathieu wrote:
>
> The fact that f^h is no longer separated changes the solution of my pde
> drastically:
> If I run my program with the analytical f or with f^h on the axis parallel
> grid, the output is quite similar. This is not so if I run it with f^h defined
> on the rotated grid.
> And since the axis parallel grid is able to capture the separation,
> I hoped the rotated grid can capture it too.
> But if I understood you correctly, there is no obvious workaround for this
> issue , right?

When you discretize a function via interpolation of projection, you generally
lose most structural properties of the function -- for example, it will have a
different mean value, it may lose symmetry if your mesh is not symmetric, it
will only be in H^1 even though the original function may have been in
C^infty, and vector functions may no longer be divergence free after
projection. Some of these issues can be addressed by using specialized
operations (e.g., divergence free projections), but you generally have to
expect that whatever special property your original function has, the
interpolated or projected function will not any longer possess.

That happens to be one of the many things you just have to live with when
working with discrete functions.

[Mathieu]

" When you discretize a function via interpolation of projection, you generally

lose most structural properties of the function -- for example, it will have a
different mean value, it may lose symmetry if your mesh is not symmetric, it
will only be in H^1 even though the original function may have been in
C^infty, and vector functions may no longer be divergence free after
projection. Some of these issues can be addressed by using specialized
operations (e.g., divergence free projections), but you generally have to
expect that whatever special property your original function has, the
interpolated or projected function will not any longer possess."

So in my case there is no special operation in your mind to make sure f^h at the rotated grid is separated in x and y (e.g., ensuring a certain aspect ratio of the cells,...)? 

-Math

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---

You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/f85c0f37-b442-45e7-39de-dab5b86ff121%40colostate.edu.

[Wolfgang Bangerth]

On 6/9/23 15:26, Mathieu wrote:
>
> So in my case there is no special operation in your mind to make sure f^h at
> the rotated grid is separated in x and y (e.g., ensuring a certain aspect
> ratio of the cells,...)?

I don't have any suggestions, sorry.