Accesing cell from cell index

[Victoria Ponce]

Hi,

I’m working with a Deal code and I have a doubt in the implementation which hopefully one of  you can help me.

I have certain active cells that have property A. I want to loop over these cells and then  loop over its neighbours so I can assign them property A. 

The way that I wanted to implement this is by:

1:Create a vector "Cells_with_A" with the cell indices of the cells that have property A 

2: create a loop with respect to this vector and get the cells from the index and once I have the cell call its neighbours.

However, I can’t find which function access the cell from the cell index..

I’m looking for something like 

cell_from_index(cell->index())

but I haven’t find it in the documentation. Is there a function like that? 

Thanks!

[Bruno Turcksin]

Victoria,

On Monday, May 9, 2016 at 8:01:00 AM UTC-4, Victoria Ponce wrote:

I have certain active cells that have property A. I want to loop over these cells and then  loop over its neighbours so I can assign them property A. 

The way that I wanted to implement this is by:

1:Create a vector "Cells_with_A" with the cell indices of the cells that have property A 

It would be simpler to keep the iterator instead of the cell indices. Also the index of the cell is useful only if the mesh is uniform or if you know the level of the cell.

2: create a loop with respect to this vector and get the cells from the index and once I have the cell call its neighbours.

However, I can’t find which function access the cell from the cell index..

I’m looking for something like 

cell_from_index(cell->index())

You need to know both the level and the index of the cell. Then, you can use this constructor https://dealii.org/developer/doxygen/deal.II/classCellAccessor.html#a41aa4117adc0c7beeb75c9169318886d

Best

Bruno

[Wolfgang Bangerth]

Victoria,

> I have certain active cells that have property A. I want to loop over these
> cells and then loop over its neighbours so I can assign them property A.
>
>
> The way that I wanted to implement this is by:
>
> 1:Create a vector "Cells_with_A" with the cell indices of the cells that have
> property A
>
> 2: create a loop with respect to this vector and get the cells from the index
> and once I have the cell call its neighbours.
>
>
> However, I can’t find which function access the cell from the cell index..
>
>
> I’m looking for something like
>
>
> cell_from_index(cell->index())
>
>
> but I haven’t find it in the documentation. Is there a function like that?

There is no such function because you can't reconstruct a cell just from its
index -- cells are only uniquely described by the index and the level.

But why not do something of the following kind (pseudo-code):
std::vector<active_cell_iterator> cells_with_prop_A;
for (cell=...)
if (cell has property A)
cells_with_prop_A.push_back (cell);

for (std::vector<a_c_i>::iterator p=cells_with_prop_A.begin();
p!=cells_with_prop_A.end(); ++p)
for (face=...)
if ((*p)->at_boundary(face) == false)
set property A on (*p)->neighbor(face)

Best
W.

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

[Martin Kronbichler]

Dear Victoria,

> However, I can’t find which function access the cell from the cell
> index..
>
> I’m looking for something like
>
> cell_from_index(cell->index())
>

If you have meshes that are adaptively refined, cell->index() does not
uniquely qualify a cell. You need to have the pair of cell->level() and
cell->index().

Assume you have the cell_index given (on a uniform mesh), or the pair of
the two values. Then you construct a cell_iterator by the following
simple code:
typename DoFHandler<dim>::cell_iterator cell(&triangulation,
cell_level, cell_index, &dof_handler);

The constructor you are interested in gets a triangulation, the two
integers for level and index within the level, and the DoFHandler, see here:
https://www.dealii.org/developer/doxygen/deal.II/classDoFAccessor.html

Best,
Martin

[Victoria Ponce]

Thanks everyone for the replies!

I decided to change the pseudo code for my particular case, however since I want a fast** way to see if a cell has value A I decided to have a vector with length the num of active cells that acts as an indicator function "has_A" and  just have  has_A[cell->index()]=1; for the ones that have it 

I don't have problem when I call has_A[(*p)->neighbor(face_no)->index()]==0 to check if the value A is assigned 

however when I want to see if the neighbor of neighbor value is assigned has_A[(*p)->neighbor(face_no)->neighbor(face_no2)->index()]; then it tells me I'm dereferencing a cell iterator... but I think I'm doing exactly the same as before..

I hope anyone has an idea of what I'm doing wrong.

**I think it is faster to just check if has_A[cell] is 0 or 1 rather than checking if cell is in cells_with_prop_A. The downside is that the vector has_A is now big and sparse.....

[Victoria Ponce]

I'm using a uniform mesh so the cells are well defined by their index so there is no problem there..

[Bruno Turcksin]

2016-05-09 12:45 GMT-04:00 Victoria Ponce <anavicto...@gmail.com>:
> however when I want to see if the neighbor of neighbor value is assigned
> has_A[(*p)->neighbor(face_no)->neighbor(face_no2)->index()]; then it tells
> me I'm dereferencing a cell iterator... but I think I'm doing exactly the
> same as before..
>
> I hope anyone has an idea of what I'm doing wrong.

Please send the exact error message and the corresponding lines of code.

Best,

Bruno

[Victoria Ponce]

Property A is having a velocity assigned, so has_A=has_vel cells_with_A=cells_with_vel and I also have cell_vel_value

if Q is a cell that has assigned a velocity strength I'm extending the velocity strength to a neighbor cell of it N as the mean value of velocity from its neighbor cells (which include Q and sometimes others)

The code is the following:

for ( typename   std::vector<typename hp::DoFHandler<dim>::active_cell_iterator>::iterator p=cells_with_vel.begin();  p!=cells_with_vel.end(); ++p)

{

   for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)

   {

   if ((*p)->at_boundary(face_no) == false)  //if the face is not at the boundary

   {   if ( has_vel[(*p)->neighbor(face_no)->index()]==0)//if the neighbor at this face does not have a velocity assigned

       {

       double  vel_value=0;

       int n_neighbors_non_zero_vel= 0;

       for (unsigned int face_no2=0; face_no2<GeometryInfo<dim>::faces_per_cell; ++face_no2)

       {// loop to assign the mean value of the neighbor cells that have velocity assigned

          vel_value+=cells_vel_value[(*p)->neighbor(face_no)->neighbor(face_no2)->index()];

          n_neighbors_non_zero_vel+=has_vel[(*p)->neighbor(face_no)->neighbor(face_no2)->index()];

       }

       vel_value/=vel_value/n_neighbors_non_zero_vel;

       }

   }

  }

}

The error is this:

An error occurred in line <934> of file </Users/victoria/software/dealii-dev/include/deal.II/grid/tria_iterator.h> in function

    Accessor &dealii::TriaRawIterator<dealii::DoFCellAccessor<dealii::hp::DoFHandler<2, 2>, false> >::operator*() [Accessor = dealii::DoFCellAccessor<dealii::hp::DoFHandler<2, 2>, false>]

The violated condition was: 

    Accessor::structure_dimension!=Accessor::dimension || state() == IteratorState::valid

The name and call sequence of the exception was:

    ExcDereferenceInvalidCell(accessor)

Additional Information: 

You tried to dereference a cell iterator for which this is not possible. More information on this iterator: level=-1, index=-1, state=past_the_end

[Bruno Turcksin]

2016-05-09 15:21 GMT-04:00 Victoria Ponce <anavicto...@gmail.com>:>
> The name and call sequence of the exception was:
>
> ExcDereferenceInvalidCell(accessor)
>
> Additional Information:
>
> You tried to dereference a cell iterator for which this is not possible.
> More information on this iterator: level=-1, index=-1, state=past_the_end

You went too far, the cell does not exist. You check that the current
cell is not at the boundary but you also need to check that the
neighbor cell is not on the boundary if you want to ask for its own
neighbor.

Best,

Bruno

[Victoria Ponce]

You are right!! Just had to check with respect to the boundary and it works! Thanks!

Thanks all and sorry for the trouble!

> --
> 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/azGWeZrIgR0/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[hank...@gmail.com]

Dr. Bangerth, Ponce, Turcksin...

I have a similar problem.

Could you help me...?

In my case, property is not velocity but material ids. (In my problem, there are two material, one is plasma the other is vacuum)

The definition of boundary for these two material is contour of psi_l   

where let say I know value of psi_l. Actually, psi_l comes from |grad(solution)|=0 point.

Anyway, If the condition is satisfied, I will set material ids.

I attached initial mesh for my calculation (I'm going to use picard iteration, so boundary between two materials would change as iterated)

In the picture, yellow zone represents vacuum and dark green zone represents plasma.

In dark green zone(plasma), solution must be always larger than psi_l.

But, the problem is that in yellow zone(vacuum) solution could be either lager or smaller than psi_l.

So, I can't simply use the condition that solution >= psi_l to define plasma zone.

The way I work out to solve this problem is following...

1.find the position of center of plasma.

2.starting from center of plasma, neighbor cells would be accessed by using cell->neighbor()

3.check whether these neighbor cells satisfy condition for solution>=psi_l on all quadrature points 

4.If it is satisfied, set new material id (that would be iteration number)  as index for plasma, else set material id as zero(zero represents vacuum) 

5.remove previous material id (that is iteration-1)

the code that I thought is as follows...

 FEValues<dim> fe_values (fe, quadrature_formula,

                           update_values    |  update_gradients |

                           update_quadrature_points  |  update_JxW_values);

std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> next_cell;

for(unsigned int face_no=0;face_no<GeometryInfo<dim>::faces_per_cell;++face_no)

     next_cell.push_back(center_cell->neighbor(face_no));                                                           //center_cell is kind of starting point for the bellow loop (it points to center of plasma)

while(condition)                                                                                                             //acutually I don't know what condition could stop this loop appropriately

{

    std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> neigh_cell;             //renew the neigh_cell for next calculation

    neigh_cell.resize(next_cell);

    std::copy(next_cell.begin(), next_cell.end(), neigh_cell.begin());

    std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> next_cell;            //after next_cell is copied to neigh_cell, renew the next_cell for next calculation

    for(std::vector<typename hp::DoFHandler<dim>::active_cell_iterator>::iterator cell=neigh_cell.begin(); cell!=neigh_cell.end();++cell)    //loop over the all neighbor cell

    {

           fe_values.reinit(cell);

           fe_value.get_function_values(solution, solution_tmp);

           int determine = 0;                                                           //this int wolud be used to determine whether this cell is plasma or not

           for(unsigned int q_index=0; q_index<quadrature_formula.size(); q_index++)

           {

                if(solution_tmp[q_index]>=psi_l)

                        determine++;

                else

                {

                       cell->set_material_ids(0);                            //this cell belong to vacuum

                       break;

                }

           }

           if(determine==quadrature_formula.size())                  //if all quadrature points satisfy the condition for solution >= psi_l

           {

                cell->set_material_id(iteration);                           //I'm going to use iteration number as new index for plasma

                for(unsigned int face_no=0;face_no<GeometryInfo<dim>::faces_per_cell;++face_no)

                {

                    if(cell->neighbor(face_no)->material_id() != iteration)       //exclude already checked cell and let's refer to next neighbor cells

                       next_cell.push_back(cell->neighbor(face_no));       

                }

           }

    }

                                                                                  // I'm going to remove previous material id that is iteration-1

  typename DoFHandler<dim>::active_cell_iterator

  cell = dof_handler.begin_active(),

  endc = dof_handler.end(); 

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

  {

       if(cell->material_id != 0 && cell->material_id !=iteration)

                cell->set_material_id(0);                                      //this cell now belong to vacuum

  }

}

So, my question is....

1.Do you think it is right way to set new material index as it is iterated?

2.As indicated in the code, Could you give me some comments for appropriate condition for me to stop this loop when all cells in the plasma are accessed?

3.Do you think the my code could also be used when the locally refined mesh is used...?   

Thanks in advance for your answer.

Kyusik.

[Jean-Paul Pelteret]

Dear Kyusik,

Two points:

1. Please create a new thread for a new question. Even though yours might be related to this one, this question has been asked and answered and its definitely tangential to the discussion here.

2. You have asked a question wanting general commentary on code that you have written. That is not the purpose of this forum - I ask to to please read the page on posting guidelines, in particular the section on What sort of questions are appropriate for the forum?. 

Your questions (1) and (2) pertain directly to my point (2). Your point (3) suggests that you have not taken the time to test your code thoroughly nor look at the examples in the documentation to assess for yourself whether this might be plausible.

The developers have given you many tools and much information in the way of the documentation and tutorials in order to piece together a solution to your problem. Yes, there is not an example of absolutely every permutation of problem that can be solved with deal.II. It is your responsibility as a researcher to perform your due diligence and assess and tackle your own problem as best you can. When you have particular problems in doing so then here's the place to ask. Its extremely difficult and, often more problematically, time consuming for any one of us to analyse the problem that you've put forward, consider what you've done and what the implications of each step are especially when there's no running code to test. I recall that in the past you've said that you need "quick results" for your supervisor, but its not a realistic nor sustainable expectation that someone else is going to sort out all of these minor details for you. Sometimes you just have to work through it yourself until you've run out of ideas.

So as far as this problem goes, I'm going to throw it back to you by asking the standard questions that any reasonable bystander would ask:

- Have you produced a simple test case to test this part of code? Does it work as you expect it to? 

- If not, how does it not work? Is the program crashing or is there a problem with your algorithm? 

- If the former, what is the error message? If the latter, then you obviously need to think a bit more about how you've structured your code to tackle your problem. Simplify it more if possible. 

When you've solved your simplified problem then rejoice for a moment and implement the next complexity. The cycle starts over again.

J-P

[hank...@gmail.com]

Dear Pelteret,

I'm sorry to have bothered you...

I should have been more considerate before asking something.

Thank you.

Kyusik.