Discontinuous Taylor basis functions

[Praveen C]

Hello

Does deal.II have discontinuous Taylor basis functions. In 1-D they are defined like this

xc = centroid of cell

h = width of cell

The basis functions are

1

(x-xc)/h

(x-xc)^2/h^2 - 1/12

etc.

The basis functions are defined in terms of physical coordinates. Except the first one, all the others have zero integral on the cell.

If they are not available in deal.II is it possible to implement such an FE space ?

Thanks

praveen

[kevin...@gmail.com]

I believe that's what the FE_DGP classes are for.

http://www.dealii.org/8.0.0/doxygen/deal.II/classFE__DGP.html

[Praveen C]

FE_DGP is quite close to the Taylor basis. But FE_DGP basis functions are defined in terms of the coordinates in the reference element. I want them to be defined directly in terms of the physical coordinates. 

praveen

--
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.
For more options, visit https://groups.google.com/groups/opt_out.

[Tobias Leicht]

Hello Praveen,

Conceptually, the closest thing to a Taylor basis is
FE_DGPNonparametric. Here, the functions are defined in physical
coordinates, but the resulting basis is not orthogonal, in general. In
multiple dimensions and under the influence of mappings the
orthogonality is difficult to achieve. It typically requires a numerical
orthogonalization on each element, which is either compute or storage
intensive and destroys some of the abstraction in deal.II, but it can be
done.
However, FE_DGPNonparametric is probably not what you want to use,
either, because, citing from the documentation:
"Warning: this class does not work properly, yet. Don't use it! [..] The
purpose of this class is experimental, therefore the implementation will
remain incomplete."

Best,
Tobias

[Praveen C]

Hello Tobias

Orthogonality would be nice to have but it is not absolutely essential. I need all the basis functions to have zero integral except for one. More importantly I want to be able to evaluate basis functions defined on one cell in a neighbouring cell. This is possible if the basis functions are defined directly in physical coordinates. This is related to another question I asked

https://groups.google.com/forum/#!topic/dealii/gaongNBSXfY

If I do not demand orthogonality, is it possible to implement such an fespace. What all classes would I need to implement for this to work ?

Thanks

praveen

[Praveen C]

On Mon, Aug 19, 2013 at 5:33 PM, Praveen C <cpra...@gmail.com> wrote:

Hello Tobias

Orthogonality would be nice to have but it is not absolutely essential. I need all the basis functions to have zero integral except for one. More importantly I want to be able to evaluate basis functions defined on one cell in a neighbouring cell. This is possible if the basis functions are defined directly in physical coordinates. This is related to another question I asked

https://groups.google.com/forum/#!topic/dealii/gaongNBSXfY

If I do not demand orthogonality, is it possible to implement such an fespace. What all classes would I need to implement for this to work ?

Thanks

An example of Taylor basis in two dimensions which is not orthogonal is given here

http://www.sciencedirect.com/science/article/pii/S0021999108003604

See pages 8878-8879. 

Also on page 1180 here

http://onlinelibrary.wiley.com/doi/10.1002/fld.3707/abstract

Thanks

praveen

[Tobias Leicht]

Praveen,

> Orthogonality would be nice to have but it is not absolutely essential.
> I need all the basis functions to have zero integral except for one.

But even this amounts to an orthogonalization, even if you only require
orthogonality w.r.t. the first basis function, which is the constant
function in this case.

> More importantly I want to be able to evaluate basis functions defined
> on one cell in a neighbouring cell.

That is not a problem.

> If I do not demand orthogonality, is it possible to implement such an
> fespace.

It depends on what you want to hear:
a) Sure, I've done so, even with ortho-normalization. For DG with the
BR2 discretization of viscous terms this is a really good choice.
b) Not really, because it violates the conceptual separation of mesh,
mapping, quadrature and finite element space. The basis you are looking
for requires information about the mapped mesh element which is
different from the information that is required for reference element
basis functions. This is problematic, for example, if you evaluate basis
functions on faces: here you only know a face quadrature rule (via
FEFaceValues), but you need a cell quadrature to evaluate integrals over
the cell (which is how you define the cell centroids). deal.II is really
built around the idea of reference elements.

> What all classes would I need to implement for this to work ?

This is reaching into the past, so the following might be incorrect, but
I think that I only had to change FESystem (pass the mapped quadrature
points correctly to the base FEs) and implement the new finite element
space class itself. A good point to start is making the
FE_DGPNonparametric work.

Best,
Tobias

PS: My implementation is in a version of deal.II that is quite far from
the trunk so I cannot simply provide a patch.

[Praveen C]

Thanks. I tried the DGPNonparametric for a scalar convection problem using RKDG and it works. This comes closest to what I want.

To clarify my understanding, the basis function in this fe space are (in 2d)

1

x

y

etc. Is that so ? I would like to do some shifting and scaling of the basis functions wrt the cell center (xc,yc) and cell size (h = cell->diameter())

1

(x-xc)/h

(y-yc)/h

Thanks

praveen

--
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+unsubscribe@googlegroups.com.

[Wolfgang Bangerth]

On 08/30/2013 11:08 AM, Praveen C wrote:
> Thanks. I tried the DGPNonparametric for a scalar convection problem using
> RKDG and it works. This comes closest to what I want.
>
> To clarify my understanding, the basis function in this fe space are (in 2d)
>
> 1
> x
> y
>
> etc. Is that so ?

I believe that is correct.

> I would like to do some shifting and scaling of the basis
> functions wrt the cell center (xc,yc) and cell size (h = cell->diameter())
>
> 1
> (x-xc)/h
> (y-yc)/h

I'm sure that's possible. The challenge will be to get information about the
cell's diameter and center into the place where the finite element class fills
the data structures of FEValues. This would be in the fill_fe_values()
function and friends of the FEDGPNonparametric class. I don't know the element
well enough to tell you in more detail what to do there, but that's the place
where you'll have to look.

Of course, given the fact that your shape functions are so simple, another
option would be to use the finite element when telling the DoFHandler and
matrices about how many DoFs there are and where; then, whenever you need to
assemble something, just bypass FEValues and the implement and replace by hand
fe_values.shape_value/grad by the definition of the function above. When using
higher order elements, instead of coding these functions by hand, the various
polynomial classes can probably do some of that for you. I think that's the
way I'd go -- just bypass FEValues if your shape functions aren't mapped from
a reference cell.

Best
W.

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

[Praveen C]

Hello Wolfgang

I actually started writing my own FEValues class. But then I could not use other deal.ii functions like in VectorTools or DataOut.

I have modified DGP_Nonparametric to shift and scale the basis functions. They look like this

template <int dim, int spacedim>

double

FE_DGPNonparametric<dim,spacedim>::shape_value (const typename                  Triangulation<dim,spacedim>::cell_iterator & cell, const unsigned int i,

                                                const Point<dim> &p) const

{

  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));

  const Point<dim> pp = (p - cell->center()) / cell->diameter();

  return polynomial_space.compute_value(i, pp);

}

template <int dim, int spacedim>

void

FE_DGPNonparametric<dim,spacedim>::fill_fe_values (

  const Mapping<dim,spacedim> &,

  const typename Triangulation<dim,spacedim>::cell_iterator & cell,

  const Quadrature<dim> &,

  typename Mapping<dim,spacedim>::InternalDataBase &,

  typename Mapping<dim,spacedim>::InternalDataBase &fedata,

  FEValuesData<dim,spacedim> &data,

  CellSimilarity::Similarity &/*cell_similarity*/) const

{

  // convert data object to internal

  // data for this class. fails with

  // an exception if that is not

  // possible

  Assert (dynamic_cast<InternalData *> (&fedata) != 0,

          ExcInternalError());

  InternalData &fe_data = static_cast<InternalData &> (fedata);

  const UpdateFlags flags(fe_data.current_update_flags());

  Assert (flags & update_quadrature_points, ExcInternalError());

  const unsigned int n_q_points = data.quadrature_points.size();

  double h = cell->diameter();

  if (flags & (update_values | update_gradients))

    for (unsigned int i=0; i<n_q_points; ++i)

      {

        const Point<dim> p = (data.quadrature_points[i] - cell->center())/h;

        polynomial_space.compute(p, /* data.quadrature_points[i] */

                                 fe_data.values, fe_data.grads, fe_data.        grad_grads);

        for (unsigned int k=0; k<this->dofs_per_cell; ++k)

          {

            if (flags & update_values)

              data.shape_values[k][i] = fe_data.values[k];

            if (flags & update_gradients)

              data.shape_gradients[k][i] = fe_data.grads[k]/h;

            if (flags & update_hessians)

              data.shape_hessians[k][i] = fe_data.grad_grads[k]/h/h;

          }

      }

}

This seems to work for my example. I have similarly modified the other fill_fe_face_values and fill_fe_subface_values functions. Are these changes correct and sufficient, or do I need to change anything else ?

Thanks
praveen

[Wolfgang Bangerth]

> I have modified DGP_Nonparametric to shift and scale the basis functions. They
> look like this
>
> template <int dim, int spacedim>
> double
> FE_DGPNonparametric<dim,spacedim>::shape_value (const typename
> Triangulation<dim,spacedim>::cell_iterator & cell, const unsigned int i,
> const Point<dim> &p) const
> {
> Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
> const Point<dim> pp = (p - cell->center()) / cell->diameter();
> return polynomial_space.compute_value(i, pp);
> }
>
> template <int dim, int spacedim>
> void
> FE_DGPNonparametric<dim,spacedim>::fill_fe_values (
> const Mapping<dim,spacedim> &,
> const typename Triangulation<dim,spacedim>::cell_iterator & cell,

Ah, yes, I had forgotten that fill_fe_values actually takes a reference to the
cell you are working on as argument.

I haven't looked at your implementing in every detail, but it seems reasonable.