Hi,
I find a better way to realize this by using the MPI distributed vector
like the finite element solution vector, where I save the values in each
cell dof by dof using cell->set_dof_values() and then read the dof values
using cell->get_dof_values(). The latter function can exactly read the dof
values at the inner partition boundary. Only one problem that at a few
dofs, the values read are different from the values those were set(I’m not
sure these dofs are on the partition boundary or other positions), but I
can check these dof values and then adjust them to the right ones.
在 2016年11月9日星期三 UTC+8上午10:06:07，Jack写道：
>
> Hi, everyone
>
> I’m not familiar to ASPECT, could someone tell me how to get the data of
> ghost cells?
> In the ASPECT, I find that it is implemented through the following
> functions to obtain ghost particles
> in ~/aspect/source/particle/world.cc
>
> template <int dim>
>
> std::map<types::subdomain_id,unsigned int>
>
> World<dim>::get_subdomain_id_to_neighbor_map() const
>
> {
>
> }
>
> template <int dim>
>
> void
>
> World<dim>::exchange_ghost_particles()
> {
>
> }
>
> template <int dim>
>
> void
>
> World<dim>::send_recv_particles
> (const std::vector<std::vector<std::pair<types::LevelInd,Particle <dim> >
> > > &send_particles,
>
> std::vector<std::pair<types::LevelInd,
> Particle<dim> > > &received_particles)
>
> {
>
> }
>
> *Presently, I just have a simple data(for example, an object of
> vector<valuetype> on each cell), which is not so complex like the particle,
> how can I realize this this easily?*
>
>
>
> Thanks very much!
>
>
> May modify directly in the following codes
>
>
> template <int dim>
> std::map<types::subdomain_id, unsigned int>
> World<dim>::get_subdomain_id_to_neighbor_map() const
> {
> std::map<types::subdomain_id, unsigned int>
> subdomain_id_to_neighbor_map;
> const std::set<types::subdomain_id> ghost_owners =
> this->get_triangulation().ghost_owners();
> std::set<types::subdomain_id>::const_iterator ghost_owner =
> ghost_owners.begin();
>
> for (unsigned int neighbor_id=0; neighbor_id<ghost_owners.size();
> ++neighbor_id,++ghost_owner)
> {
>
> subdomain_id_to_neighbor_map.insert(std::make_pair(*ghost_owner,neighbor_id));
> }
> return subdomain_id_to_neighbor_map;
> }
>
>
>
>
>
> template <int dim>
> void
> World<dim>::exchange_ghost_particles()
> {
> TimerOutput::Scope timer_section(this->get_computing_timer(),
> "Particles: Exchange ghosts");
>
> // First clear the current ghost_particle information
> ghost_particles.clear();
>
> const std::map<types::subdomain_id, unsigned int>
> subdomain_to_neighbor_map(get_subdomain_id_to_neighbor_map());
>
> std::vector<std::vector<std::pair<types::LevelInd, Particle<dim> > >
> > ghost_particles_by_domain(subdomain_to_neighbor_map.size());
> std::vector<std::set<unsigned int> >
> vertex_to_neighbor_subdomain(this->get_triangulation().n_vertices());
>
> typename DoFHandler<dim>::active_cell_iterator
> cell = this->get_dof_handler().begin_active(),
> endc = this->get_dof_handler().end();
> for (; cell != endc; ++cell)
> {
> if (cell->is_ghost())
> for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell;
> ++v)
>
> vertex_to_neighbor_subdomain[cell->vertex_index(v)].insert(cell->subdomain_id());
> }
>
> cell = this->get_triangulation().begin_active();
> for (; cell != endc; ++cell)
> {
> if (!cell->is_ghost())
> {
> std::set<unsigned int> cell_to_neighbor_subdomain;
> for (unsigned int v=0;
> v<GeometryInfo<dim>::vertices_per_cell; ++v)
> {
>
> cell_to_neighbor_subdomain.insert(vertex_to_neighbor_subdomain[cell->vertex_index(v)].begin(),
>
> vertex_to_neighbor_subdomain[cell->vertex_index(v)].end());
> }
>
> if (cell_to_neighbor_subdomain.size() > 0)
> {
> const std::pair< const typename
> std::multimap<types::LevelInd,Particle <dim> >::iterator,
> const typename
> std::multimap<types::LevelInd,Particle <dim> >::iterator>
> particle_range_in_cell =
> particles.equal_range(std::make_pair(cell->level(),cell->index()));
> for (std::set<types::subdomain_id>::const_iterator
> domain=cell_to_neighbor_subdomain.begin();
> domain != cell_to_neighbor_subdomain.end();
> ++domain)
> {
> const unsigned int neighbor_id =
> subdomain_to_neighbor_map.find(*domain)->second;
> ghost_particles_by_domain[neighbor_id].insert(
> ghost_particles_by_domain[neighbor_id].end(),
> particle_range_in_cell.first,
> particle_range_in_cell.second);
> }
> }
> }
> }
>
> std::vector<std::pair<types::LevelInd, Particle<dim> > >
> received_ghost_particles;
> send_recv_particles(ghost_particles_by_domain,
> received_ghost_particles);
>
> ghost_particles.insert(received_ghost_particles.begin(),
> received_ghost_particles.end());
> }
> *I'm not clear about the particles.equal_range().*
> *The following send_receive particles is mots difficult for me to
> understand*
>
> template <int dim>
> void
> World<dim>::send_recv_particles(const
> std::vector<std::vector<std::pair<types::LevelInd,Particle <dim> > > >
> &send_particles,
> std::vector<std::pair<types::LevelInd,
> Particle<dim> > > &received_particles)
> {
> // Determine the communication pattern
> const std::vector<types::subdomain_id> neighbors
> (this->get_triangulation().ghost_owners().begin(),
>
> this->get_triangulation().ghost_owners().end());
> const unsigned int n_neighbors = neighbors.size();
>
> Assert(n_neighbors == send_particles.size(),
> ExcMessage("The particles to send to other processes should
> be sorted into a vector "
> "containing as many vectors of particles as there
> are neighbor processes. This "
> "is not the case for an unknown reason. Contact
> the developers if you encounter "
> "this error."));
>
> unsigned int n_send_particles = 0;
> for (unsigned int i=0; i<n_neighbors; ++i)
> n_send_particles += send_particles[i].size();
>
> #if DEAL_II_VERSION_GTE(8,5,0)
> const unsigned int cellid_size = sizeof(CellId::binary_type);
> const unsigned int particle_size =
> property_manager->get_particle_size() + integrator->get_data_size() +
> cellid_size;
> #else
> const unsigned int particle_size =
> property_manager->get_particle_size() + integrator->get_data_size();
> #endif
>
> // Determine the amount of data we will send to other processors
> std::vector<unsigned int> n_send_data(n_neighbors);
> std::vector<unsigned int> n_recv_data(n_neighbors);
>
> std::vector<unsigned int> send_offsets(n_neighbors);
> std::vector<unsigned int> recv_offsets(n_neighbors);
>
> // Allocate space for sending and receiving particle data
> std::vector<char> send_data(n_send_particles * particle_size);
> void *data = static_cast<void *> (&send_data.front());
>
> for (types::subdomain_id neighbor_id = 0; neighbor_id < n_neighbors;
> ++neighbor_id)
> {
> send_offsets[neighbor_id] = reinterpret_cast<std::size_t> (data)
> - reinterpret_cast<std::size_t> (&send_data.front());
>
> // Copy the particle data into the send array
> typename std::vector<std::pair<types::LevelInd,Particle <dim> >
> >::const_iterator
> cell_particle = send_particles[neighbor_id].begin(),
> end_particle = send_particles[neighbor_id].end();
> for (; cell_particle != end_particle; ++cell_particle)
> {
> #if DEAL_II_VERSION_GTE(8,5,0)
> const typename
> parallel::distributed::Triangulation<dim>::cell_iterator cell
> (&this->get_triangulation(),
>
> cell_particle->first.first,
>
> cell_particle->first.second);
> const CellId::binary_type cellid = cell->id().template
> to_binary<dim>();
> memcpy(data, &cellid, cellid_size);
> data = static_cast<char *>(data) + cellid_size;
> #endif
> cell_particle->second.write_data(data);
> data = integrator->write_data(data,
> cell_particle->second.get_id());
> }
> n_send_data[neighbor_id] = reinterpret_cast<std::size_t> (data)
> - send_offsets[neighbor_id] - reinterpret_cast<std::size_t>
> (&send_data.front());
> }
>
> // Notify other processors how many particles we will send
> std::vector<MPI_Request> n_requests(2*n_neighbors);
> for (unsigned int i=0; i<n_neighbors; ++i)
> MPI_Irecv(&(n_recv_data[i]), 1, MPI_INT, neighbors[i], 0,
> this->get_mpi_communicator(), &(n_requests[2*i]));
> for (unsigned int i=0; i<n_neighbors; ++i)
> MPI_Isend(&(n_send_data[i]), 1, MPI_INT, neighbors[i], 0,
> this->get_mpi_communicator(), &(n_requests[2*i+1]));
> MPI_Waitall(2*n_neighbors,&n_requests[0],MPI_STATUSES_IGNORE);
>
> // Determine how many particles and data we will receive
> unsigned int total_recv_data = 0;
> for (unsigned int neighbor_id=0; neighbor_id<n_neighbors;
> ++neighbor_id)
> {
> recv_offsets[neighbor_id] = total_recv_data;
> total_recv_data += n_recv_data[neighbor_id];
> }
>
> // Set up the space for the received particle data
> std::vector<char> recv_data(total_recv_data);
>
> // Exchange the particle data between domains
> std::vector<MPI_Request> requests(2*n_neighbors);
> unsigned int send_ops = 0;
> unsigned int recv_ops = 0;
>
> for (unsigned int i=0; i<n_neighbors; ++i)
> if (n_recv_data[i] > 0)
> {
> MPI_Irecv(&(recv_data[recv_offsets[i]]), n_recv_data[i],
> MPI_CHAR, neighbors[i], 1,
> this->get_mpi_communicator(),&(requests[send_ops]));
> send_ops++;
> }
>
> for (unsigned int i=0; i<n_neighbors; ++i)
> if (n_send_data[i] > 0)
> {
> MPI_Isend(&(send_data[send_offsets[i]]), n_send_data[i],
> MPI_CHAR, neighbors[i], 1,
> this->get_mpi_communicator(),&(requests[send_ops+recv_ops]));
> recv_ops++;
> }
> MPI_Waitall(send_ops+recv_ops,&requests[0],MPI_STATUSES_IGNORE);
>
> // Put the received particles into the domain if they are in the
> triangulation
> const void *recv_data_it = static_cast<const void *>
> (&recv_data.front());
>
> #if !DEAL_II_VERSION_GTE(8,5,0)
> const std::vector<std::set<typename
> Triangulation<dim>::active_cell_iterator> >
>
> vertex_to_cells(GridTools::vertex_to_cell_map(this->get_triangulation()));
> const std::vector<std::vector<Tensor<1,dim> > >
> vertex_to_cell_centers(vertex_to_cell_centers_directions(vertex_to_cells));
>
> std::vector<unsigned int> neighbor_permutation;
> #endif
>
> while (reinterpret_cast<std::size_t> (recv_data_it) -
> reinterpret_cast<std::size_t> (&recv_data.front()) < total_recv_data)
> {
> #if DEAL_II_VERSION_GTE(8,5,0)
> CellId::binary_type binary_cellid;
> memcpy(&binary_cellid, recv_data_it, cellid_size);
> const CellId id(binary_cellid);
> recv_data_it = static_cast<const char *> (recv_data_it) +
> cellid_size;
>
> const typename
> parallel::distributed::Triangulation<dim>::active_cell_iterator cell =
> id.to_cell(this->get_triangulation());
> const Particle<dim>
> recv_particle(recv_data_it,property_manager->get_particle_size());
> recv_data_it = integrator->read_data(recv_data_it,
> recv_particle.get_id());
> #else
> Particle<dim>
> recv_particle(recv_data_it,property_manager->get_particle_size());
> recv_data_it = integrator->read_data(recv_data_it,
> recv_particle.get_id());
> const std::pair<const typename
> parallel::distributed::Triangulation<dim>::active_cell_iterator,
> Point<dim> > current_cell_and_position =
> GridTools::find_active_cell_around_point<>
> (this->get_mapping(),
>
> this->get_triangulation(),
>
> recv_particle.get_location());
> typename
> parallel::distributed::Triangulation<dim>::active_cell_iterator cell =
> current_cell_and_position.first;
>
> recv_particle.set_reference_location(current_cell_and_position.second);
>
> // GridTools::find_active_cell_around_point can find a different
> cell than
> // particle_is_in_cell if the particle is very close to the
> boundary
> // therefore, we might get a cell here that does not belong to
> us.
> // But then at least one of its neighbors belongs to us, and the
> particle
> // is extremely close to the boundary of these two cells. Look
> in the
> // neighbor cells for the particle.
> if (!cell->is_locally_owned())
> {
> const unsigned int closest_vertex =
> get_closest_vertex_of_cell(cell,recv_particle.get_location());
> const unsigned int closest_vertex_index =
> cell->vertex_index(closest_vertex);
> Tensor<1,dim> vertex_to_particle =
> recv_particle.get_location() - cell->vertex(closest_vertex);
> vertex_to_particle /= vertex_to_particle.norm();
>
> const unsigned int n_neighbor_cells =
> vertex_to_cells[closest_vertex_index].size();
>
> neighbor_permutation.resize(n_neighbor_cells);
> for (unsigned int i=0; i<n_neighbor_cells; ++i)
> neighbor_permutation[i] = i;
>
> std::sort(neighbor_permutation.begin(),
> neighbor_permutation.end(),
> std_cxx11::bind(&compare_particle_association<dim>,
> std_cxx11::_1,
> std_cxx11::_2,
>
> std_cxx11::cref(vertex_to_particle),
>
> std_cxx11::cref(vertex_to_cell_centers[closest_vertex_index])));
>
> // Search all of the cells adjacent to the closest vertex of
> the previous cell
> // Most likely we will find the particle in them.
> for (unsigned int i=0; i<n_neighbor_cells; ++i)
> {
> try
> {
> typename std::set<typename
> Triangulation<dim>::active_cell_iterator>::const_iterator neighbor_cell =
> vertex_to_cells[closest_vertex_index].begin();
> std::advance(neighbor_cell,neighbor_permutation[i]);
> const Point<dim> p_unit =
> this->get_mapping().transform_real_to_unit_cell(*neighbor_cell,
>
> recv_particle.get_location());
> if (GeometryInfo<dim>::is_inside_unit_cell(p_unit))
> {
> cell = *neighbor_cell;
> recv_particle.set_reference_location(p_unit);
> break;
> }
> }
> catch (typename Mapping<dim>::ExcTransformationFailed &)
> {}
> }
> }
> #endif
>
> const types::LevelInd found_cell =
> std::make_pair(cell->level(),cell->index());
>
> received_particles.push_back(std::make_pair(found_cell,
> recv_particle));
> }
>
> AssertThrow(recv_data_it == &recv_data.back()+1,
> ExcMessage("The amount of data that was read into new
> particles "
> "does not match the amount of data sent
> around."));
> }
>
> Thanks so much!
>
> Best
>
> Jack
>
>
>
>
>>
>> Just for reference, Rene Gassmoeller implemented something similar for
>> particles recently: There, we store particles on every locally owned
>> cell, but
>> there are times when one also wants to access particles on ghost cells.
>> This
>> is implemented here:
>> https://github.com/geodynamics/aspect/pull/1216
>> https://github.com/geodynamics/aspect/pull/1253
>>
>> It's not terribly complicated because one just has to gather the
>> information
>> on cells that are ghosts on some other processor, and then send this
>> information to all relevant processors.
>>
>> Best
>> W.
>>
>> --
>> ------------------------------------------------------------------------
>> Wolfgang Bangerth email: [email protected]
>> www: http://www.math.colostate.edu/~bangerth/
>>
>>
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
