finite difference method to calculate the gradient of density in cells by reference file of cook_membrane.cc

24 views
Skip to first unread message

Lance Zhang

unread,
Jul 2, 2023, 8:36:50 AM7/2/23
to deal.II User Group
Hello ,

I would like to use finite difference method to calculate the gradient of density in cells.

But I got issue when the program was processed at the part of compute_gradient();

Here is the code below  about density and gradient.


//>>>>>>>>>>>>>>>Density>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
template <int dim,typename NumberType>

void Solid<dim,NumberType>::compute_density()
{
  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler_ref.begin_active(),
                                                 endc = dof_handler_ref.end();

  for (; cell != endc; ++cell)
  {
    // Compute the grid density for each cell,
    density[cell->active_cell_index()] = 1.0 / cell->measure();
  }
}



//>>>>>>>>>>>>>>>>>>>Gradient>>>>>>>>>>>>>>>>>>>>>>>>
template <int dim,typename NumberType>


void Solid<dim,NumberType>::compute_gradient()
{
  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler_ref.begin_active(),
                                                 endc = dof_handler_ref.end();

  for (; cell != endc; ++cell)
  {
    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
    {
      if (cell->face(face)->at_boundary())
      {
        // Compute the gradient at each boundary face by differencing the neighboring cell densities
        const typename DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face);
        gradient[cell->face(face)->index()] = density[neighbor->active_cell_index()] - density[cell->active_cell_index()];
      }
    }
  }
}
//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

The error comes out at the part of compute_gradient().
----------------------------------------------------------------------------------------------------------------
Error information:
Thread 1 "cook_membrane" received signal SIGSEGV, Segmentation fault.
0x00007ffff4bfcdaa in dealii::CellAccessor<3, 3>::active_cell_index() const () from /lib/x86_64-linux-gnu/libdeal.ii.so.9.3.2
----------------------------------------------------------------------------------------------------------------
I tried to reduce the numbers of grids and dofs ,but I sill got the same issue.

in the part of void Solid<dim,NumberType>::run(),
I add  the code below  

 template <int dim,typename NumberType>
  void Solid<dim,NumberType>::run(){
 .//other codes
 .//other codes
 .//other codes
    density.reinit(dof_handler_ref.n_dofs());
    gradient.reinit(dof_handler_ref.n_dofs());
    //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    print_vertical_tip_displacement();

    std::cout<<"Density claculation"<<std::endl;
    compute_density();//sucessfully
    std::cout<<"Gradient calculation"<<std::endl;
    compute_gradient();//error comes out
}

Could anyone provide any hint or suggestions?

Thanks in advance!
Best regards
Lance

Lance Zhang

unread,
Jul 2, 2023, 9:30:45 AM7/2/23
to deal.II User Group

Hello,
 
later , I tried to use the number of active cells,
density.reinit(triangulation.n_active_cells());
gradient.reinit(triangulation.n_active_cells());

But I still got the same issue.

In last code pasted in the first email I think I have to replace dof_handler_ref.n_dofs() with triangulation.n_active_cells(),because dof_handler_ref.n_dofs()  is about the number of dofs instead of activae cells.


Could anyone provide any suggestions?

Best regards
Lance

Wolfgang Bangerth

unread,
Jul 2, 2023, 11:10:27 AM7/2/23
to dea...@googlegroups.com

Lance:
let me start by saying that you need to learn to debug problems yourself. This
forum can help you with questions that are specific to deal.II, but we do not
have the resources to help you with general programming questions. Finding the
root cause of a segmentation fault is one of these issue: You need to learn to
figure this out yourself.


> But I got issue when the program was processed at the part of compute_gradient();
>
> Here is the code below  about density and gradient.
> [...]
>         gradient[cell->face(face)->index()] =
> density[neighbor->active_cell_index()] - density[cell->active_cell_index()];

One hint is here: If the index to density is the number of a cell, then the
size of that vector clearly needs to be the number of cells. If the index to
the gradient vector is the number of a face, then the size of that vector
needs to be the number of faces.

I don't actually know how you would get the number of faces of a
triangulation. More importantly, perhaps, is that I don't know what you would
do with such a vector. But that's a question for you to answer.

Best
W.

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


Reply all
Reply to author
Forward
0 new messages