Thanks a lot for your reply.
It seems I can not upload files so I just paste the code here. I hope this is not too messy to look at. And I have described what I really want to do in the code.
the plugin cc file:
#ifndef _aspect_particle_property_my_particle_h
#define _aspect_particle_property_my_particle_h
#include <aspect/simulator_access.h>
#include <aspect/particle/property/interface.h>
#include <aspect/postprocess/visualization/error_indicator.h>
namespace aspect
{
namespace Particle
{
namespace Property
{
template <int dim>
class MyParticle : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
void
initialize_one_particle_property (const Point<dim> &position,
std::vector<double> &particle_properties) const override;
virtual
void
update_particle_property (const unsigned int data_position,
const Vector<double> &solution,
const std::vector<Tensor<1,dim> > &gradients,
typename ParticleHandler<dim>::particle_iterator &particle) const;
UpdateTimeFlags
need_update () const override;
UpdateFlags
get_needed_update_flags () const override;
std::vector<std::pair<std::string, unsigned int> >
get_property_information() const override;
};
}
}
}
#endif
namespace aspect
{
namespace Particle
{
namespace Property
{
template <int dim>
void
MyParticle<dim>::initialize_one_particle_property(const Point<dim> &position,
std::vector<double> &data) const
{
data.push_back(0.0);
}
template <int dim>
void
MyParticle<dim>::update_particle_property (const unsigned int data_position,
const Vector<double> &solution,
const std::vector<Tensor<1,dim> > &gradients,
typename ParticleHandler<dim>::particle_iterator &particle) const
{
const Postprocess::VisualizationPostprocessors::ErrorIndicator<dim> &error_indicator =
this->get_postprocess_manager().template get_matching_postprocessor<Postprocess::VisualizationPostprocessors::ErrorIndicator<dim> >();
// Do anything with error_indicator, I don't really care. Or just don't do anything.
// For example, do the following code to intepolate error indicator for every particle:
/*
const std::pair<std::string, Vector<float> *> error_indicator_cell_data = error_indicator->execute();
typename parallel::distributed::Triangulation<dim>::cell_iterator found_cell;
found_cell = particle->get_surrounding_cell(this->get_triangulation());
data[manager.get_data_info().get_position_by_field_name("my particle")] = *error_indicator_cell_data.second(found_cell.active_cell_index());
*/
// What I really want to do is a little bit complex than this:
// 1st: I have already defined a properly working new way of writing visualization postprocessors in './source/postprocess/visualization.cc'
// which out put a field data. It is similar to DataPostprocessor in the way that it can be evaluate on every point.
// But it is different from DataPostprocessor in that it is not a point-wise function but needs all other cells' info to caculate it.
// 2nd: Then I want to intepolate the value for every particle.
//
// But I think I'll be fine as long as I can run this simple test case.
}
template <int dim>
UpdateTimeFlags
MyParticle<dim>::need_update() const
{
return update_time_step;
}
template <int dim>
UpdateFlags
MyParticle<dim>::get_needed_update_flags() const
{
return update_values;
}
template <int dim>
std::vector<std::pair<std::string, unsigned int>>
MyParticle<dim>::get_property_information() const
{
std::vector<std::pair<std::string, unsigned int>> property_information(1, std::make_pair("my particle", 1));
return property_information;
}
} // namespace Property
} // namespace Particle
} // namespace aspect
// explicit instantiations
namespace aspect
{
namespace Particle
{
namespace Property
{
ASPECT_REGISTER_PARTICLE_PROPERTY(MyParticle,
"my particle",
"A test case.")
}
} // namespace Particle
} // namespace aspect
the prm input file:
set Additional shared libraries = usr_plugin/libmyParticle.so
set Dimension = 2
set Output directory = testCase
set Resume computation = false
set End time = 500e6
set Use years in output instead of seconds = true
set CFL number = 1
subsection Geometry model
set Model name = box
subsection Box
set X extent = 990e3
set Y extent = 660e3
set X repetitions = 12
set Y repetitions = 8
end
end
subsection Material model
set Model name = simple
end
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 9.8
end
end
subsection Boundary velocity model
set Tangential velocity boundary indicators = bottom, left, top, right
end
subsection Initial temperature model
set List of model names = function
subsection Function
set Function constants = u0=0.020, kappa=31.6, T0=273, T1=1623, Lmax=9999e3, H=660e3, L=1320e3, t0=100e6, x0=0e3, w=10e3
set Function expression = T1 - (T1-T0) * erfc((H-y)/2.0/sqrt(kappa*max(3e6, x/u0)))
end
end
subsection Boundary temperature model
set Fixed temperature boundary indicators = top, bottom
set List of model names = initial temperature
end
##################### Mesh refinement #########################
subsection Mesh refinement
set Initial global refinement = 3
set Initial adaptive refinement = 0
set Time steps between mesh refinement = 0
end
##################### Postprocessing ########################
subsection Postprocess
set List of postprocessors = visualization, particles
subsection Visualization
set List of output variables = error indicator
set Time between graphical output = 0
set Output format = vtu
end
subsection Particles
set Data output format = vtu
set Number of particles = 5e4
set Load balancing strategy = add particles # remove and add particles
set Minimum particles per cell = 20
set Time between data output = 0
set List of particle properties = my particle
end
end
subsection Termination criteria
set Checkpoint on termination = true
set Termination criteria = user request, end step
set End step = 0
subsection User request
set File name = terminate
end
end
