Setting smooth coefficients in space

[Konstantin Ladutenko]

Hi all,

What is the correct way to set smooth coefficients in space for high p-order of elements simulation in deal.ii?

E.g. in Step-6 nonconstant coefficient is a step function  alpha=20 for R<0.5 and alpha=1 otherwise (0<R<1).  How to set it to be a smooth function, e.g. some polynomial function of order 'n' like alpha=20*(1-R^n)+1? So it should be enough to use few high-order elements for the whole model to converge to the smooth solution.

I`ve checked a number of tutorials (e.g. http://www.dealii.org/developer/doxygen/deal.II/step_4.html#Righthandsideandboundaryvalues , equation data section in step-27, and few others - all of them use a point-wise approximation of an analytic function  (so it is tuned for h-refinements), which is practically the same as a step function for a given mesh. 

This is actually very important when solving any open-boundary problem using perfectly-matching layer (PML), as soon as the latter needs to have a parameter grading (conductivity). Any ideas\examples how it can be done using deal.ii?

As an external reference, there is a PhD thesis http://arizona.openrepository.com/arizona/handle/10150/195940 with a description of a single cell PML. In a short, it is preferable to use one high-p  element to set PML from the point of view of performance (wave reflection) to DoF ratio.

Best regards,

Konstantin Ladutenko 

[Jean-Paul Pelteret]

Hi Konstantin,

One option would be to employ a combination of a transfer function and level set functions. With this approach you effectively provide a regularised expression for the material parameters. You can see how this is done equation 18 in this paper and equation 32 in this paper. If you can't access either of the articles then let me know and I can provide some further details.

I hope that this helps!

Regards,

J-P

[Wolfgang Bangerth]

Konstantin,

> What is the correct way to set smooth coefficients in space for high p-order
> of elements simulation in deal.ii?
>
> E.g. in Step-6 nonconstant coefficient is a step function alpha=20 for R<0.5
> and alpha=1 otherwise (0<R<1). How to set it to be a smooth function, e.g.
> some polynomial function of order 'n' like alpha=20*(1-R^n)+1? So it should be
> enough to use few high-order elements for the whole model to converge to the
> smooth solution.
>
> I`ve checked a number of tutorials
> (e.g. http://www.dealii.org/developer/doxygen/deal.II/step_4.html#Righthandsideandboundaryvalues
> , equation data section in step-27, and few others - all of them use a
> point-wise approximation of an analytic function (so it is tuned for
> h-refinements), which is practically the same as a step function for a given
> mesh.

No, that's not true. When you say "pointwise", what I think you really mean is
that we only evaluate the coefficient at individual (quadrature) points. But
that's good enough, because when you do quadrature -- say for an integrand on
a cell K,

A_{ij}^K = \sum_q a(x_q) \nabla phi_i(x_q) \nabla phi_j(x_q) JxW(q)

then that is equivalent to computing the integral of a polynomial that is
equal to

a(x_q) \nabla phi_i(x_q) \nabla phi_j(x_q)

i.e., that *interpolates* these point values. This polynomial is of course
smooth. In particular, this corresponds to a *smooth* interpolation of the
coefficient a(x_q), even though you only evaluate it at quadrature points.

What kind of polynomial you use to interpolate the integrand depends on how
many quadrature points you use. If you need a very accurate interpolation of
the coefficient, just use a higher order Gauss formula.


So in some sense, you've got it exactly the wrong way around: If, in step-6,
we use a discontinuous coefficient, then the finite element method using
quadrature really sees only a polynomial approximation of that discontinuous
coefficient on every mesh. (Though that polynomial approximation of course
tends to a discontinuous function under mesh refinement.)

Best
W.


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

[Konstantin Ladutenko]

Dear J-P and Wolfgang,

Thank you for your answers! As I see the actual problem is to define a step function for material parameters. In FEM simulations mesh is often aligned to the problem`s inner interfaces (like in a "Better mesh" extension of step-6). Is it possible to define for such a case with a discontinuous coefficient using continuous element?  E.g. is it valid to take coefficient evaluation out of the q_index cycle and use   

current_coefficient = coefficient<dim> (cell->center())   

in step-6 code for the case of mesh aligned to the interface of coefficient drop?

Best regards,

Kostya

[Jean-Paul Pelteret]

Hi Konstantin,

Thank you for your answers!

You're welcome!  

 

As I see the actual problem is to define a step function for material parameters. In FEM simulations mesh is often aligned to the problem`s inner interfaces (like in a "Better mesh" extension of step-6). Is it possible to define for such a case with a discontinuous coefficient using continuous element?  E.g. is it valid to take coefficient evaluation out of the q_index cycle and use   

current_coefficient = coefficient<dim> (cell->center())   

in step-6 code for the case of mesh aligned to the interface of coefficient drop?

Sure, it obviously depends on what you're trying to model but in many cases I would say that this is legitimate. The underlying material description is, for most problems (e.g. those with no discontinuities internal to a cell), completely independent of the chosen basis functions. If your mesh exactly discretises your geometry into material patches, then could divide up your triangulation (note: nothing to do with DoFs) into these different material zones using cell->set_material_id(). You could then choose your material parameters based on a cell->material_id(), which is exactly the same as what you had proposed. You are, in essence, assuming that the material parameters are piecewise constant over each cell. Whether or not this is a good approximation of how your material parameters vary in space is entirely up to you to decide, but I use this assumption (and the approach that I laid out) almost exclusively in my own work. A counter-example for which you'd have to be more careful are anisotropic (e.g. fibrous) media, like is shown here.

 Best regards,

J-P

[Konstantin Ladutenko]

Dear J-P and Wolfgang,

I`ve tried to incorporate your answers to the documentation, see https://github.com/dealii/dealii/pull/4339

Best regards,

Kostya

вторник, 2 мая 2017 г., 20:05:28 UTC+3 пользователь Jean-Paul Pelteret написал:

[Jean-Paul Pelteret]

Hi Konstantin,

Thanks a lot for the contribution! I'll be taking a look at it later in the day.

Best,

J-P