How to compute derivatives of the composite function with the reference cell mapping

[Anton]

Hello,

I would like to use the NonMatching::QuadratureGenerator<dim> class directly, as I need to generate multiple (as in many!, unfortunately) nonmatching quadratures for the same cell.  Therefore I would like to avoid the route of interpolating the level set function onto the whole mesh etc, as described in CutFEM tutorial.

Therefore I am thinking about calling the "generate" function, which only needs a level set function and a bounding box as inputs.  According to the documentation the level set function should provide the gradient and the Hessian.  I would like to do this in the reference cell coordinates, so that I can simply use this quadrature later via FEView as one normally does.

 

 * The bounding box I presume is the hypercube [0,1]^{dim} for the meshes I am interested in. (Can one request them for the reference cell and not the triangulation cell?)

 * The level set function should also be reasonably straightforward.  I know the function in the physical coordinates (say, f(x)) and all its derivatives.  I can also get the map from reference to physical coordinates via FEView, let us call this map x=M(y).  The issue is how to extract the derivatives of this map so that the chain rule can be applied?  I am mostly working with Q1 maps, so of course I could figure out the derivatives "by hand".  I am just wondering if there is a more intelligent (code-wise) way of doing this?

[Simon Sticko]

Hi.

I have done this in the past. I dug up the wrapper function I had for it, in case we decide that we should add it to the library:

https://github.com/dealii/dealii/pull/15838

The problem with this approach is that evaluating the reference space level set function gets very expensive. QuadratureGenerator uses many function calls and we have to use the Jacobian and its derivative in the transformation. This makes this approach significantly slower than using an interpolated level set function.

If you are working with a perfectly Cartesian background I would expect that it would actually be faster to generate the quadrature in real space and then transform the quadrature to reference space before passing it to FEValues. If you want to do this the functions cell->bounding_box() and BoundingBox::real_to_unit(..) are useful.

Best,
Simon

> --
> The deal.II project is located at http://www.dealii.org/ <http://www.dealii.org/>
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en <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 <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/4ac38563-e467-4f1d-8cb4-ee35379515a6n%40googlegroups.com <https://groups.google.com/d/msgid/dealii/4ac38563-e467-4f1d-8cb4-ee35379515a6n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Anton Evgrafov]

Dear Simon,

Thank you very much for your reply. The wrapper code looks
embarrassingly simple!

> If you are working with a perfectly Cartesian background I would expect that it would actually be faster to generate the quadrature in real space and then transform the quadrature to reference space before passing it to FEValues. If you want to do this the functions cell->bounding_box() and BoundingBox::real_to_unit(..) are useful.

That I did not think about for some reason :)
If my element is not aligned with the axes, won't I risk getting
quadrature points outside of the element though? Or is this what you
meant by the "perfect Cartesian" situation? Also I am not 100% clear
at the moment about how to change the quadrature weights after the
physical->unit transformation, as I would need access to Jacobian
determinants.

Another option I was contemplating is to apply the
NonMatching::DiscreteQuadratureGenerator<dim> on a fake triangulation
with just one element. Then I simply would need to interpolate
(evaluate) my real-space function at this one element and let the
class' implementation do the magic for me. This is probably cheaper
and about equally accurate?

/Anton

> 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/FONV6StZZus/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/f1642630-6453-2028-1345-3db743d814f2%40gmail.com.



--
Anton Evgrafov

[Simon Sticko]

Hi.

> If my element is not aligned with the axes, [...]
> is this what meant by the "perfect Cartesian" situation?

Yes, exactly. If this is not what you have, I'm afraid the idea of generating the quadratures in real space is not an option.

> Another option I was contemplating is to apply the
> NonMatching::DiscreteQuadratureGenerator<dim> on a fake triangulation
> with just one element. Then I simply would need to interpolate
> (evaluate) my real-space function at this one element and let the
> class' implementation do the magic for me. This is probably cheaper
> and about equally accurate?

That sounds doable and cheaper. If you are going for this, FEFieldFunction might be useful:
https://www.dealii.org/developer/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html

Regarding accuracy, you typically get the same order of accuracy. That is, say you that you use Lagrange elements of order p in the finite element space for you pde solution, and compare using the exact level set function with one that has been interpolated (also onto a space of elements of order p), then you typically see that

Exact level set function => error ~ C_1 h^q
Interpolated level set function => error ~ C_2 h^q

(error in some global norm)

However, I know one expert who have claimed that the difference between the constants C_1 and C_2 can be surprisingly large for some problems. But I am unaware of any publication that has actually compared this.

Just out of curiosity (and if it is not a research secret), what is your use case? Why do you need to generate many quadratures for the same cell?

Best,
Simon

[Anton Evgrafov]

Simon,

> That sounds doable and cheaper. If you are going for this, FEFieldFunction might be useful:
> https://www.dealii.org/developer/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html

I have now compared the implementation using UnitToReal function
wrapper + QuadratureGenerator vs using an auxiliary triangulation with
one reference cell onto which I interpolate my "level set function",
so to speak, + DiscreteQuadratureGenerator. The latter is waaay
faster.

> Just out of curiosity (and if it is not a research secret), what is your use case? Why do you need to generate many quadratures for the same cell?

Well, I am certainly misusing deal.II. I have a problem involving a
field living in R^{dim^2} and a field in R^{dim}. The weak form
involves products between the two, and integrals are over
dim^2-dimensional domains, which may be decomposed into two iterated
dim-dimensional integrals. The inner integral then goes over a ball
which cuts elements. This is where non-matching quadrature comes
along, but I have to compute them for each outer quadrature point, so
to speak. Also the integrand is not polynomial all, so some
quadratures are of higher degree.

The proper way of doing this would be to utilize a FEM library that
supports mesh/elements/quadratures in dim^2. Can one do this with
deal.ii? Most of the code seems to be dimension-independent, but then
there are plenty of places where e.g. Point<dim> constructors are
specialized for dimensions 1..3.

Right now I am doing a "quick and dirty" test implementation, where I
am basically manually working on a tensor product mesh with tensor
product FEM etc, which is highly inelegant. If you are
doing/interested in research, I would be happy to discuss the problem
in detail/collaborate.


/Anton

[Wolfgang Bangerth]

On 8/5/23 04:30, Anton Evgrafov wrote:
> I have a problem involving a
> field living in R^{dim^2} and a field in R^{dim}.

This may or may not be relevant to you, but it might be worth taking a look at
hyper.deal:
https://github.com/hyperdeal/hyperdeal

Best
W.

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

[Anton Evgrafov]

Thank you Wolfgang, this does look relevant! /Anton

> --
> 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/FONV6StZZus/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/0f072f81-1297-fa96-4744-d4f22908b19b%40colostate.edu.



--
Anton Evgrafov