Thanks for you reply. Here is a minimal example that does 3 three things: 
(i) interpolate a solution (some simple function in this case, but would 
actually be the result of solving a PDE) onto the DOFs of a mesh, (ii) 
interpolate the solution onto randomly placed particles, and (iii) try to 
associated the interpolated solution with particle properties.  For me, 
this code compiles but does not run. Do you know why?

I have two other related questions.

(i) In the code below I could replace the counter part with iter->get_id(). 
Is this the intended use of the ID or would we expect this to fail in more 
complicated examples where particles may be added/removed or distributed 
across processors?

(ii) Is it necessary to copy the particle properties into an std::vector or 
can we just give each particle an iterator the the 
interolatedParticleQuantities vector? The copy step seems unnecessary and 
potentially slow/burdensome if there are lots of particles.

#include <deal.II/grid/grid_generator.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>

#include <deal.II/numerics/vector_tools.h>

#include <deal.II/particles/generators.h>
#include <deal.II/particles/particle_handler.h>
#include <deal.II/particles/utilities.h>

using namespace dealii;

/// The "solution" to some PDE
template<unsigned int dim>
class SolutionFunction : public Function<dim> {
public:
  inline SolutionFunction() : Function<dim>(dim+1) {}

  inline virtual void vector_value(Point<dim> const& p, Vector<double>& 
values) const override {
    for( unsigned int i=0; i<dim; ++i ) {
      values[i] = p[i];
    }
    values[dim] = 1.0;
  }

private:
};

int main() {
  // spatial dimension
  constexpr unsigned int dim = 2;

  // number of "unknowns"
  const unsigned int ncomponents = dim+1;

  // number of particles to randomly place
  const unsigned int nparticles = 10;

  // create the mesh
  const MappingQ1<dim> mapping;
  Triangulation<dim> triangulation;
  DoFHandler<dim> dofHandler(triangulation);
  GridGenerator::hyper_cube(triangulation, 0.0, 1.0);
  triangulation.refine_global(5);
  const dealii::FESystem<dim> fe(dealii::FE_Q<dim>(1), ncomponents);
  dofHandler.distribute_dofs(fe);

  // "solve" the pde
  SolutionFunction<dim> function;
  Vector<double> solution(dofHandler.n_dofs());
  VectorTools::interpolate(dofHandler, function, solution);

  // create the particle handler
  Particles::ParticleHandler<dim> particleHandler(triangulation, mapping, 
ncomponents);

  // randomly place particles in the domain
  ConstantFunction<dim> pdf(1.0);
  Particles::Generators::probabilistic_locations(triangulation, pdf, true, 
nparticles, particleHandler, mapping);

  // map the conserved properties onto the particles
  dealii::Vector<double> 
interpolatedParticleQuantities(ncomponents*nparticles);
  Particles::Utilities::interpolate_field_on_particles(dofHandler, 
particleHandler, solution, interpolatedParticleQuantities);

  // check to make sure we get what we expect
  unsigned int part = 0;
  for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); 
++iter, ++part ) {
    std::cout << "particle " << part << " quantities: ";
    for( unsigned int i=0; i<ncomponents; ++i ) {
      std::cout << interpolatedParticleQuantities[part*ncomponents + i] << 
" ";
    }
    std::cout << std::endl;
    std::cout << "expected quantities: ";
    for( unsigned int i=0; i<dim; ++i ) {
      std::cout << iter->get_location()[i] << " ";
    }
    std::cout << " 1" << std::endl;
    std::cout << std::endl;
  }

  // now try to set the interpolated quantities as particle properties
  part = 0;
  for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); 
++iter, ++part ) {
    std::vector<double> quantities(ncomponents);
    for( unsigned int i=0; i<ncomponents; ++i ) {
      quantities[i] = interpolatedParticleQuantities[part*ncomponents + i];
    }

    iter->set_properties(quantities);
  }
}

On Tuesday, June 16, 2020 at 10:25:01 PM UTC-4, Wolfgang Bangerth wrote:
>
> On 6/16/20 1:46 PM, Andrew Davis wrote: 
> > 
> > and I have gotten what I expect.  I have also tried attaching the 
> particle 
> > properties using: 
> > 
> > unsigned int npart = 0; 
> > for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); 
> ++iter, 
> > ++part ) { 
> >    dealii::ArrayView<double> data(interpolatedParticleQuantities.begin() 
> + 
> > part*ncomponents, ncomponents); 
> >    iter->set_properties(data); 
> > } 
> > 
> > but I get the compiler error: 
> > 
> > *error: *cannot convert ‘*dealii::ArrayView<double>*’ to ‘*const 
> > std::vector<double, std::allocator<double> >&*’ 
> > 
> > 353 | iter->set_properties(*data*); 
> > 
> > 
> > Does anyone know what I am doing wrong? Or is there is a better way to 
> do 
> > this? I'm sure this is a dumb question, but I couldn't find anything in 
> the 
> > tutorials/examples. 
>
> You've already solved this one problem by converting the ArrayView to a 
> std::vector by hand in your next mail. 
>
>
>  > The violated condition was: 
>  > 
>  > property_pool != nullptr 
>  > 
>  > Additional information: 
>  > 
>  > This exception -- which is used in many places in the library -- 
> usually 
>
> That's more complicated. The exception happens in the line 
>    iter->set_properties(...) 
> and tells you that the particle 'iter' points to isn't associated with a 
> PropertyPool, i.e., the object that stores particle properties. I'm not 
> sure 
> why that is so. I believe that from your code, you set up the 
> ParticleHandler 
> to hold properties, but I can't see where the ParticleHandler's 
> PropertyPool 
> would be propagated to the newly added particles. 
>
> Can you create a minimal program that shows the issue? It doesn't have to 
> do 
> anything useful, just illustrate the error. 
>
> Best 
>   W. 
>
>
>
>
> -- 
> ------------------------------------------------------------------------ 
> Wolfgang Bangerth          email:                 bang...@colostate.edu 
> <javascript:> 
>                             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 dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a6e20774-6525-438d-b1c9-12ca256f6193o%40googlegroups.com.

Reply via email to