how to get the real coordinates of interpolation nodes when use higher-order FEM?

[FanBai]

Dear all,

When we use linear FEM, in each cell the interpolation nodes are at the vertices, so we could use cell->vertex(cell->get_fe().system_to_component_index(i).second) to find out the real coordinate of the four nodes.

however, if we use higher-order FEM, the interpolation nodes are not only at vertices, is there a way to know the real coordinates of all of them?

Best,

Fan

[Matthias Maier]

Hi,

Am 28. Apr 2013, 12:18 schrieb FanBai <baifa...@gmail.com>:

> Dear all,
>
> When we use linear FEM, in each cell the interpolation nodes are at the
> vertices, so we could use cell->vertex(cell->get_fe().system_to_component_index
> (i).second) to find out the real coordinate of the four nodes.

You're looking for FiniteElement::get_unit_support_points(). This will
give you the support points on the unit cell.

If you want to have the support points on a real cell, one possibility
is to generate a custom quadrature rule with the support points and use
FEValues::get_quadrature_points() to get the transformed support points.

(Alternatively, use Mapping::transform directly)

Best,
Matthias

[Fan Bai]

Hi Matthias,

Thanks for the advice. What I did is to create a Vector with 9 points (for second order element functions), and define every entry as: for the first 4, I define V(0)=cell->vertex(0)...V(3)=cell->vertex(3); and for the 5th-8th, I define V()=0.5*(cell->vertex(.)+cell->vertex(.)); for the last one: V(8)=0.25*(cell->vertex()+cell->vertex()+cell->vertex()+cell->vertex()).

For the 3rd order and 4th order I did the same. But apparently this is not a good way to code..

I think my way is kind of similar to the "custom quadrature rule" method.

I will try by using  FiniteElement::get_unit_support_points(), and then map these points to the real cell.

Best,

Fan

--
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.

[Wolfgang Bangerth]

> however, if we use higher-order FEM, the interpolation nodes are not only at
> vertices, is there a way to know the real coordinates of all of them?

Use
DoFTools::map_dofs_to_support_points
to get their locations all at once.

Best
W.


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

[Fan Bai]

Hi Wolfgang,

I am still confused on how to use DoFTools::map_dofs_to_support_points, I learnt from the step 34 that this function

can obtain real coordinates of all global support points.

However, my problem is I just need the coordinates of support points in certain cells. That means I do not need the global

dof_handler, I need local dofs. Can I still use this function to get coordinates of these support points? Thanks for the help!

Best,
Fan

--
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]

> I am still confused on how to use DoFTools::map_dofs_to_support_points, I
> learnt from the step 34 that this function
> can obtain real coordinates of all global support points.
>
> However, my problem is I just need the coordinates of support points in
> certain cells. That means I do not need the global
> dof_handler, I need local dofs. Can I still use this function to get
> coordinates of these support points? Thanks for the help!

Well, the function gives you a list of coordinates for *all* degrees of
freedom. If you want only some of them for a single cell, you can always get
their coordinates out of the list of all coordinates.

[Fan Bai]

OK, I see..Thanks!

[Denis Davydov]

Hi Fan,

I think you can get the map cell-wise following these steps:

typename dealii::DoFHandler<dim>::active_cell_iterator

      cell = dof_handler.begin_active(),

      endc = dof_handler.end();

const std::vector< dealii::Point< dim > > sp = fe.get_unit_support_points ();

std::vector< unsigned int > dof_indices(fe.dofs_per_cell);

for (; cell!=endc; ++cell) {

      //skip cells you are not interested...

      //get DoFs...

      cell->get_dof_indices(dof_indices);

      //get support points...

      std::vector< dealii::Point<dim>> rp;

      dealii::Quadrature<dim> q(sp);

      dealii::FEValues<dim> fe_values (fe, q, dealii::update_q_points);

      fe_values.reinit (cell);

      rp = fe_values.get_quadrature_points();

//sizes should match:

      assert (rp.size() == dof_indices.size());

}

From here you can either build dof2point or points2dof maps;

If i did everything correct, that should work for arbitrary interpolation order and probably for other finite elements 

with support points associated to DoFs.

Regards,

Denis.

To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.

[Wolfgang Bangerth]

> const std::vector< dealii::Point< dim > > sp = fe.get_unit_support_points();
> std::vector< unsigned int > dof_indices(fe.dofs_per_cell);
>
> for (; cell!=endc; ++cell) {
> //skip cells you are not interested...
> //get DoFs...
> cell->get_dof_indices(dof_indices);
>
> //get support points...
> std::vector< dealii::Point<dim>> rp;
> dealii::Quadrature<dim> q(sp);
> dealii::FEValues<dim> fe_values (fe, q, dealii::update_q_points);

Since 'sp' (and consequently 'q') is the same on every cell, you should pull
the declaration and initialization of 'q' and 'fe_values' out of the loop.
Creating an FEValues object is much more expensive than calling
fe_values.reinit().

[Denis Davydov]

Thanks for advices, Wolfgang.
That was just a quickly done debug version. I will definitely move those out of the cell loop.

31.05.2013, в 23:02, Wolfgang Bangerth <bang...@math.tamu.edu> написал(а):