* Authors: Luca Heltai, Bruno Blais, 2019
*/
+// @sect3{Include files}
+// Most of these have been introduced elsewhere, we'll comment only on the new
+// ones.
+
#include <deal.II/base/function.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/linear_operator_tools.h>
-#include <deal.II/particles/data_out.h>
-
#define FORCE_USE_OF_TRILINOS
namespace LA
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/vector_tools.h>
+// These are the only new include files w.r.t. step-60. In this tutorial,
+// the non-matching coupling between the solid and the fluid is computed using
+// an intermediate data structure that keeps track of how the quadrature points
+// of the solid evolve w.r.t. the fluid mesh. This data structure needs to keep
+// track of the position of the quadrature points on each cell describing the
+// solid domain, of the quadrature weights, and possibly of the normal vector
+// to each point, if the solid domain is of co-dimension one.
+//
+// Deal.II offers these facilities on the Particles namespace, through the
+// ParticleHandler class. ParticleHandler is a class that allows you to manage
+// a collection of particles (objects of type Particles::Particle), representing
+// a collection of points with some attached properties floating on a
+// parallel::distributed::Triangulation. The methods and classes on the
+// namespace Particles allows one to easily implement Particle In Cell methods
+// and particle tracing on distributed triangulations.
+//
+// We "abuse" this data structure to store information about the location of
+// solid quadrature points w.r.t. to the surrounding fluid grid, including
+// integration weights, and possibly surface normals. The reason why we use this
+// additional data structure is related to the fact that the solid and the fluid
+// grids are non-overlapping, and distributed independently among processes.
+//
+// In order to couple the two problems, we rely on the ParticleHandler class,
+// storing in each particle the position of a solid quadrature point (which is
+// in general not aligned to any of the fluid quadrature points), its weight,
+// and any other information that may be required to couple the two problems.
+//
+// Ownership of the solid quadrature points is inherited by the MPI partitioning
+// on the solid mesh itslef. The Particles so generated are later distributed to
+// the fluid mesh using the methods of the ParticleHandler class. This allows
+// transparent exchange of information between MPI processes about the
+// overlapping pattern between fluid cells and solid quadrature points.
+#include <deal.II/particles/data_out.h>
#include <deal.II/particles/generators.h>
#include <deal.II/particles/particle_handler.h>
{
using namespace dealii;
+ // REMOVE THIS FUNCTION ONCE #9891 is merged.
template <int dim,
int spacedim,
typename InputVectorType,
interpolated_field.compress(VectorOperation::add);
}
+ // Similiarly to what we have done in step-60, we set up a class that holds
+ // all the parameters of our problem and derive it from the ParameterAcceptor
+ // class to simplify the management and creation of parameter files.
+ //
+ // The ParameterAcceptor paradigm requires all parameters to be writeable by
+ // the ParameterAcceptor methods. In order to avoid bugs that would be very
+ // difficult to trace down (such as witing things like `time = 0` instead of
+ // `time == 0`), we declare all the parameters in an external class, which is
+ // initialized before the actual StokesImmersedProblem class, and pass it to
+ // the main class as a const reference.
template <int dim, int spacedim = dim>
class StokesImmersedProblemParameters : public ParameterAcceptor
{
public:
- StokesImmersedProblemParameters()
- : ParameterAcceptor("Stokes Immersed Problem/")
- , rhs("Right hand side", spacedim + 1)
- , angular_velocity("Angular velocity", spacedim == 3 ? spacedim : 1)
- {
- add_parameter("Velocity degree",
- velocity_degree,
- "",
- this->prm,
- Patterns::Integer(1));
-
- add_parameter("Number of time steps", number_of_time_steps);
- add_parameter("Output frequency", mod_output);
-
- add_parameter("Final time", final_time);
-
- add_parameter("Viscosity", viscosity);
-
- add_parameter("Nitsche penalty term", penalty_term);
-
- add_parameter("Initial fluid refinement",
- initial_fluid_refinement,
- "Initial mesh refinement used for the fluid domain Omega");
-
- add_parameter("Initial solid refinement",
- initial_solid_refinement,
- "Initial mesh refinement used for the solid domain Gamma");
-
- add_parameter(
- "Particle insertion refinement",
- particle_insertion_refinement,
- "Refinement of the volumetric mesh used to insert the particles");
-
- add_parameter(
- "Homogeneous Dirichlet boundary ids",
- homogeneous_dirichlet_ids,
- "Boundary Ids over which homogeneous Dirichlet boundary conditions are applied");
-
- enter_my_subsection(this->prm);
- this->prm.enter_subsection("Grid generation");
- this->prm.add_parameter("Grid one generator", name_of_grid1);
- this->prm.add_parameter("Grid one generator arguments",
- arguments_for_grid1);
-
- this->prm.add_parameter("Grid two generator", name_of_grid2);
- this->prm.add_parameter("Grid two generator arguments",
- arguments_for_grid2);
-
- this->prm.add_parameter("Particle grid generator", name_of_particle_grid);
- this->prm.add_parameter("Particle grid generator arguments",
- arguments_for_particle_grid);
- this->prm.leave_subsection();
-
- leave_my_subsection(this->prm);
-
-
-
- enter_my_subsection(this->prm);
- this->prm.enter_subsection("Refinement and remeshing");
- this->prm.add_parameter("Refinement step frequency", mod_refinement);
- this->prm.add_parameter("Refinement maximal level", max_level_refinement);
- this->prm.add_parameter("Refinement strategy",
- refinement_strategy,
- "",
- Patterns::Selection(
- "fixed_fraction|fixed_number"));
- this->prm.add_parameter("Refinement coarsening fraction",
- coarsening_fraction);
- this->prm.add_parameter("Refinement fraction", refinement_fraction);
- this->prm.add_parameter("Maximum number of cells", max_cells);
-
- this->prm.leave_subsection();
- leave_my_subsection(this->prm);
-
- // correct the default dimension for the functions
- rhs.declare_parameters_call_back.connect([&]() {
- Functions::ParsedFunction<spacedim>::declare_parameters(this->prm,
- spacedim + 1);
- });
- angular_velocity.declare_parameters_call_back.connect([&]() {
- Functions::ParsedFunction<spacedim>::declare_parameters(
- this->prm, spacedim == 3 ? spacedim : 1);
- });
- }
-
+ // The constructor is responsible for the connection between the members of
+ // this class and the corresponding entries in the ParameterHandler. Thanks
+ // to the use of the ParameterHandler::add_parameter() method, this
+ // connection is trivial, but requires all members of this class to be
+ // writeable
+ StokesImmersedProblemParameters();
+
+ // however, since this class will be passed as a const reference to the
+ // StokesImmersedProblem class, we have to make sure we can still set the
+ // time correctly in the objects derived by the Function class defined
+ // here. In order to do so, we declare both the
+ // StokesImmersedProblemParameters::rhs and
+ // StokesImmersedProblemParameters::angular_velocity members to be mutable,
+ // and define this little helper method that sets their time to the correct
+ // value.
void set_time(const double &time) const
{
rhs.set_time(time);
angular_velocity.set_time(time);
}
- unsigned int velocity_degree = 2;
- unsigned int number_of_time_steps = 1;
- double viscosity = 1.0;
- double final_time = 1.0;
- unsigned int initial_fluid_refinement = 3;
- unsigned int initial_solid_refinement = 3;
- unsigned int particle_insertion_refinement = 1;
- double penalty_term = 1e3;
- std::list<types::boundary_id> homogeneous_dirichlet_ids{0, 1, 2, 3};
- std::string name_of_grid1 = "hyper_cube";
- std::string arguments_for_grid1 = "-1: 1: false";
- std::string name_of_grid2 = "hyper_rectangle";
- std::string arguments_for_grid2 =
+ // We will use a Taylor-Hood function space of arbitrary order. This
+ // parameter is used to initialize the FiniteElement space with the corret
+ // FESystem object
+ unsigned int velocity_degree = 2;
+
+ // Instead of defining a time step increment, in this tutorial we prefer to
+ // let the user choose a final simulation time, and the number of steps in
+ // which we want to reach the final time
+ unsigned int number_of_time_steps = 1;
+ double final_time = 1.0;
+
+ // Instead of producing an output at every time step, we allow the user to
+ // select the frequency at which output is produced:
+ unsigned int output_frequency = 1;
+
+ // We allow every grid to be refined independently. In this tutorial, no
+ // physics is resolved on the solid grid, and its velocity is given as a
+ // datum. However it relatively straight forward to incorporate some
+ // elasticity model in this tutorial, and transform it in a fully fledged
+ // FSI solver.
+ unsigned int initial_fluid_refinement = 3;
+ unsigned int initial_solid_refinement = 3;
+ unsigned int particle_insertion_refinement = 1;
+
+ // The only two parameters used in the equations are the viscosity of the
+ // fluid, and the penalty term used in the Nitsche formulation:
+ double viscosity = 1.0;
+ double penalty_term = 1e3;
+
+ // By default, we create a hyper_cube without colorisation, and we use
+ // homogenous Dirichlet boundary conditions. In this set we store the
+ // boundary ids to use when setting the boundary conditions:
+ std::list<types::boundary_id> homogeneous_dirichlet_ids{0};
+
+ // We illustrate here another way to create a Triangulation from a parameter
+ // file, using the method GridGenerator::generate_from_name_and_arguments(),
+ // that takes the name of a function in the GridGenerator namespace, and its
+ // arguments as a single string representing the arguments as a tuple.
+ //
+ // The mechanism with which the arguments are parsed from and to a string is
+ // explained in detail in the Patterns::Tools::Convert class, which is
+ // used to translate from strings to most of the basic STL types (vectors,
+ // maps, tuples) and basic dealii types (Point, Tensor, BoundingBox, etc.).
+ //
+ // In general objects that can be represented by rank 1 uniform elements
+ // (i.e., std::vector<double>, Point<dim>, std::set<int>, etc.) are comma
+ // separated. Additional ranks take a semicolon, allowing you to parse
+ // strings into objects of type `std::vector<std::vector<double>>`, or,
+ // for example, `std::vector<Point<dim>>`, as `0.0, 0.1; 0.1, 0.2`. This
+ // string could be interpreted as a vector of two Point objects, or a vector
+ // of vector of doubles.
+ //
+ // When the entries are not uniform, as in the tuple case, we use a colon
+ // to separate the various entries. For example, a string like `5: 0.1, 0.2`
+ // could be used to parse an object of type `std::pair<int, Point<2>>` or a
+ // `std::tuple<int, std::vector<double>>`.
+ //
+ // In our case most of the arguments are Point objects (representing
+ // centers, corners, subdivision elements, etc.), integer values (number of
+ // subdivisions), double values (radius, lengths, etc.), or boolean options
+ // (such as the `colorize` option that many GridGenerator functions take).
+ //
+ // In the example below, we set reasonable default values, but these can be
+ // changed at run time by selecting any other supported function of the
+ // GridGenerator namespace.
+ //
+ // We do this for each of the generated grids, to be as generic as possible:
+ std::string name_of_grid1 = "hyper_cube";
+ std::string arguments_for_grid1 = "-1: 1: false";
+ std::string name_of_grid2 = "hyper_rectangle";
+ std::string arguments_for_grid2 =
dim == 2 ? "-.5, -.1: .5, .1: false" : "-.5, -.1, -.1: .5, .1, .1: false";
std::string name_of_particle_grid = "hyper_ball";
std::string arguments_for_particle_grid =
dim == 2 ? "0.3, 0.3: 0.1: false" : "0.3, 0.3, 0.3 : 0.1: false";
- // Refinement parameters
+ // Similarly, we allow for different local refinement strategies. In
+ // particular, we limit the maximum number of refinement levels, in order
+ // to control the minimum size of the fluid grid, and guarantee that it is
+ // compatible with the solid grid, and we perform local refinement based
+ // on standard error estimators on the fluid velocity field.
+ //
+ // We permit the user to choose between the
+ // two most common refinement strategies, namely `fixed_number` or
+ // `fixed_fraction`, that refer to the methods
+ // GridRefinement::refine_and_coarsen_fixed_fraction() and
+ // GridRefinement::refine_and_coarsen_fixed_number().
+ //
+ // Refinement may be done every few time steps, instead of continuosly, and
+ // we control this value by the `refinement_frequency` parameter:
int max_level_refinement = 5;
std::string refinement_strategy = "fixed_fraction";
double coarsening_fraction = 0.3;
double refinement_fraction = 0.3;
unsigned int max_cells = 1000;
- int mod_refinement = 5;
- int mod_output = 1;
-
+ int refinement_frequency = 5;
+
+ // These two functions are used to control the source term of Stokes flow
+ // and the angular velocity at which we move solid. In a more realistic
+ // simulation, the solid velocity or its deformation would come from the
+ // solution of an auxiliary problem on the solid domain. In this example
+ // step we leave this part aside, and simply impose a fixed rotational
+ // velocity field on the immersed solid, governed by function that can be
+ // specified in the parameter file:
mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>> rhs;
mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
angular_velocity;
}; // namespace Step70
+ // Once the angular velocity is provided as a Function object, we reconstruct
+ // the pointwise solid velocity thrugh the following class.
template <int spacedim>
class SolidVelocity : public Function<spacedim>
{
public:
SolidVelocity(const Functions::ParsedFunction<spacedim> &angular_velocity)
: angular_velocity(angular_velocity)
- {}
+ {
+ static_assert(spacedim > 1,
+ "Cannot instatiate SolidVelocity for spacedim == 1");
+ }
virtual double value(const Point<spacedim> &p,
unsigned int component = 0) const
velocity = cross_product_3d(p, omega);
}
-
- if (spacedim == 2)
+ else if (spacedim == 2)
{
double omega = angular_velocity.value(p, 0);
const Functions::ParsedFunction<spacedim> &angular_velocity;
};
-
+ // Similarly, we assume that the incremental solid displacement can be
+ // computed simply by a one step time integration process (here using a
+ // trivial forward Euler method), so that at each time step, the solid simply
+ // displaces by `v*dt`.
template <int spacedim>
class SolidDisplacement : public Function<spacedim>
{
: Function<spacedim>(spacedim)
, angular_velocity(angular_velocity)
, time_step(time_step)
- {}
+ {
+ static_assert(spacedim > 1,
+ "Cannot instatiate SolidDisplacement for spacedim == 1");
+ }
virtual double value(const Point<spacedim> &p,
unsigned int component = 0) const
double time_step;
};
+ // We are now ready to introduce the main class of our tutorial program.
template <int dim, int spacedim = dim>
class StokesImmersedProblem
{
StokesImmersedProblem(
const StokesImmersedProblemParameters<dim, spacedim> &par);
+ // As usual, we leave a single public entry point to the user: the run
+ // method. Everything else is left private, and accessed through the run
+ // method itself.
void run();
private:
void make_grid();
+
+ // These two methods are new w.r.t. previous examples, and initiliaze the
+ // ParticleHandler objects used in this class. We have two such objects: one
+ // is a passive tracer, used to plot the trajectories of fluid particles,
+ // while the the other is composed of the actual solid quadrature points,
+ // and represent material particles of the solid.
void setup_tracer_particles();
void setup_solid_particles();
+
+ // The setup is split in two parts: create all objects that are needed once
+ // per simulation,
void initial_setup();
+ // followed by all objects that need to be reinitialized at every refinement
+ // step.
void setup_dofs();
+
+ // The assembly rutine is identical to other Stokes assembly rutines,
void assemble_stokes_system();
+ // with the exception of the Nistche restriction part, which exploits one of
+ // the particle handlers to integrate on a non-matching part of the fluid
+ // domain, corresponding to the position of the solid.
void assemble_nitche_restriction();
+
void solve();
- void refine_grid();
+
+ // The refine_and_transfer() method is called only every
+ // `refinement_frequency` steps, and makes sure that all the fields
+ // that were computed on the time step before refinement are transfered
+ // correctly to the new grid. This includes vector fields, as well as
+ // particle information.
+ void refine_and_transfer();
void output_results(const unsigned int cycle, const double time) const;
void
const unsigned int iter,
const double time) const;
+ // As noted before, make sure we cannot modify this object from within this
+ // class, by making it a const reference.
const StokesImmersedProblemParameters<dim, spacedim> ∥
MPI_Comm mpi_communicator;
- std::unique_ptr<FiniteElement<spacedim>> fe1;
- std::unique_ptr<FiniteElement<dim, spacedim>> fe2;
-
- parallel::distributed::Triangulation<spacedim> tria1;
- parallel::distributed::Triangulation<dim, spacedim> tria2;
-
- DoFHandler<spacedim> dh1;
- DoFHandler<dim, spacedim> dh2;
-
- std::unique_ptr<MappingFEField<dim, spacedim>> mapping2;
-
- std::vector<IndexSet> owned1;
- std::vector<IndexSet> owned2;
-
- std::vector<IndexSet> relevant1;
- std::vector<IndexSet> relevant2;
-
- IndexSet owned_tracer_particles;
- IndexSet relevant_tracer_particles;
+ // For the current implemenation, only `fluid_fe` would be really necessary.
+ // For completeness, and to allow easy extension, we keep also the
+ // `solid_fe` around, which is however initialized to a FE_Nothing finite
+ // element space, i.e., one that has no degrees of freedom.
+ //
+ // We declare both finite element spaces as unique pointers, to allow their
+ // generation after StokesImmersedProblemParameters has been initialized. In
+ // particular, we assume that they are filled only after initial_setup() has
+ // been called.
+ std::unique_ptr<FiniteElement<spacedim>> fluid_fe;
+ std::unique_ptr<FiniteElement<dim, spacedim>> solid_fe;
+
+ // This is one of the main novelty w.r.t. the tutorial step-60. Here we
+ // assume that both the solid and the fluid are fully distributed
+ // triangulations. This allows the problem to scale to a very large number
+ // of degrees of freedom, at the cost of communicating all the overlapping
+ // regions between non matching triangulations. This is especially tricky,
+ // since we make no assumptions on the relative position or distribution of
+ // the various subdomains. In particular, we assume that ever process owns
+ // only a part of the solid_tria, and only a part of the fluid_tria, not
+ // necessarily in the same physical region, and not necessarily overlapping.
+ //
+ // In order to couple the overlapping regions, we exploit the facilities
+ // implemented in the ParticleHandler class.
+ parallel::distributed::Triangulation<spacedim> fluid_tria;
+ parallel::distributed::Triangulation<dim, spacedim> solid_tria;
+
+ DoFHandler<spacedim> fluid_dh;
+ DoFHandler<dim, spacedim> solid_dh;
+
+ std::unique_ptr<MappingFEField<dim, spacedim>> solid_mapping;
+
+ // Similarly to how things are done in step-32, we use a block system to
+ // treat the Stokes part of the problem, and follow very closely what was
+ // done there.
+ std::vector<IndexSet> fluid_owned_dofs;
+ std::vector<IndexSet> solid_owned_dofs;
+
+ std::vector<IndexSet> fluid_relevant_dofs;
+ std::vector<IndexSet> solid_relevant_dofs;
AffineConstraints<double> constraints;
LA::MPI::BlockVector locally_relevant_solution;
LA::MPI::BlockVector system_rhs;
+ // For every tracer particle, we need to compute the velocity field in its
+ // current position, and update its position using a discrete time stepping
+ // scheme. We do this using distributed linear algebra objects, where the
+ // owner of a particle is set to be equal to the process that generated that
+ // particle at time t=0. This information is stored for every process in the
+ // `owned_tracer_particles` IndexSet, that indicates which particles the
+ // current process owns.
+ //
+ // Once the particles have been distributed around to match the process that
+ // owns the region where the particle lives, we will need read access from
+ // that process on the corresponding velocity field. We achieve this by
+ // filling a read only velocity vector field, that contains the relevant
+ // information in ghost entries. This is achieved using the
+ // `relevant_tracer_particles` IndexSet, that keeps track of how things
+ // change during the simulation, i.e., it keeps track of where particles
+ // that I own have ended up being, and who owns the particles that ended up
+ // in my subdomain.
+ //
+ // While this is not the most efficient strategy, we keep it this way to
+ // illustrate how things would work in a real FSI problem. If a particle
+ // is linked to a specific solid degree of freedom, we are not free to
+ // choose who owns it, and we have to communicate this information around.
+ // We illustrate this here, and show that the communication pattern is
+ // point-to-point, and negligible in terms of total cost of the algorithm.
+ IndexSet owned_tracer_particles;
+ IndexSet relevant_tracer_particles;
+
+ // These vectors are used to store the particles velocities (read-only, with
+ // ghost entries) and their displacement (read/write, no ghost entries).
LA::MPI::Vector tracer_particle_velocities;
LA::MPI::Vector relevant_tracer_particle_displacements;
+ // We fix once the quadrature formula that is used to integrate the solid
+ // domain.
std::unique_ptr<Quadrature<dim>> quadrature_formula;
Particles::ParticleHandler<dim, spacedim> tracer_particle_handler;
const StokesImmersedProblemParameters<dim, spacedim> &par)
: par(par)
, mpi_communicator(MPI_COMM_WORLD)
- , tria1(mpi_communicator,
- typename Triangulation<spacedim>::MeshSmoothing(
- Triangulation<spacedim>::smoothing_on_refinement |
- Triangulation<spacedim>::smoothing_on_coarsening))
- , tria2(mpi_communicator,
- typename Triangulation<dim, spacedim>::MeshSmoothing(
- Triangulation<dim, spacedim>::smoothing_on_refinement |
- Triangulation<dim, spacedim>::smoothing_on_coarsening))
- , dh1(tria1)
- , dh2(tria2)
+ , fluid_tria(mpi_communicator,
+ typename Triangulation<spacedim>::MeshSmoothing(
+ Triangulation<spacedim>::smoothing_on_refinement |
+ Triangulation<spacedim>::smoothing_on_coarsening))
+ , solid_tria(mpi_communicator,
+ typename Triangulation<dim, spacedim>::MeshSmoothing(
+ Triangulation<dim, spacedim>::smoothing_on_refinement |
+ Triangulation<dim, spacedim>::smoothing_on_coarsening))
+ , fluid_dh(fluid_tria)
+ , solid_dh(solid_tria)
, pcout(std::cout,
(Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
, computing_timer(mpi_communicator,
template <int dim, int spacedim>
void StokesImmersedProblem<dim, spacedim>::make_grid()
{
- GridGenerator::generate_from_name_and_arguments(tria1,
+ GridGenerator::generate_from_name_and_arguments(fluid_tria,
par.name_of_grid1,
par.arguments_for_grid1);
- tria1.refine_global(par.initial_fluid_refinement);
+ fluid_tria.refine_global(par.initial_fluid_refinement);
- GridGenerator::generate_from_name_and_arguments(tria2,
+ GridGenerator::generate_from_name_and_arguments(solid_tria,
par.name_of_grid2,
par.arguments_for_grid2);
- tria2.refine_global(par.initial_solid_refinement);
+ solid_tria.refine_global(par.initial_solid_refinement);
}
template <int dim, int spacedim>
particles_dof_handler.distribute_dofs(particles_fe);
// Create the particle handler associated with the fluid triangulation
- tracer_particle_handler.initialize(tria1,
+ tracer_particle_handler.initialize(fluid_tria,
StaticMappingQ1<spacedim>::mapping);
// The generation of the global bounding boxes requires an all-to-all
// communication
auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
- tria1, IteratorFilters::LocallyOwnedCell());
+ fluid_tria, IteratorFilters::LocallyOwnedCell());
auto global_bounding_boxes =
Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
relevant_tracer_particles = owned_tracer_particles;
// Now make sure that upon refinement, particles are correctly transferred
- tria1.signals.pre_distributed_refinement.connect(std::bind(
+ fluid_tria.signals.pre_distributed_refinement.connect(std::bind(
&Particles::ParticleHandler<spacedim>::register_store_callback_function,
&tracer_particle_handler));
- tria1.signals.post_distributed_refinement.connect(std::bind(
+ fluid_tria.signals.post_distributed_refinement.connect(std::bind(
&Particles::ParticleHandler<dim,
spacedim>::register_load_callback_function,
&tracer_particle_handler,
template <int dim, int spacedim>
void StokesImmersedProblem<dim, spacedim>::setup_solid_particles()
{
- QGauss<dim> quadrature(fe1->degree + 1);
+ QGauss<dim> quadrature(fluid_fe->degree + 1);
// In codimension one case, we store also the normal, else only the
// quadrature weight.
const unsigned int n_properties = (dim == spacedim) ? 1 : spacedim + 1;
- solid_particle_handler.initialize(tria1,
+ solid_particle_handler.initialize(fluid_tria,
StaticMappingQ1<dim>::mapping,
n_properties);
std::vector<Point<spacedim>> quadrature_points_vec(
- quadrature.size() * tria2.n_locally_owned_active_cells());
+ quadrature.size() * solid_tria.n_locally_owned_active_cells());
std::vector<std::vector<double>> properties(
- quadrature.size() * tria2.n_locally_owned_active_cells(),
+ quadrature.size() * solid_tria.n_locally_owned_active_cells(),
std::vector<double>(n_properties));
UpdateFlags flags = update_JxW_values | update_quadrature_points;
if (spacedim > dim)
flags |= update_normal_vectors;
- FEValues<dim, spacedim> fe_v(*fe2, quadrature, flags);
+ FEValues<dim, spacedim> fe_v(*solid_fe, quadrature, flags);
unsigned int cell_index = 0;
- for (const auto &cell : dh2.active_cell_iterators())
+ for (const auto &cell : solid_dh.active_cell_iterators())
if (cell->is_locally_owned())
{
fe_v.reinit(cell);
// Distribute the local points to the processor that owns
// them on the triangulation
auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
- tria1, IteratorFilters::LocallyOwnedCell());
+ fluid_tria, IteratorFilters::LocallyOwnedCell());
auto global_bounding_boxes =
Utilities::MPI::all_gather(mpi_communicator, my_bounding_box);
// Now make sure that upon refinement, particles are correctly transferred
- tria1.signals.pre_distributed_refinement.connect(std::bind(
+ fluid_tria.signals.pre_distributed_refinement.connect(std::bind(
&Particles::ParticleHandler<spacedim>::register_store_callback_function,
&solid_particle_handler));
- tria1.signals.post_distributed_refinement.connect(std::bind(
+ fluid_tria.signals.post_distributed_refinement.connect(std::bind(
&Particles::ParticleHandler<dim,
spacedim>::register_load_callback_function,
&solid_particle_handler,
{
TimerOutput::Scope t(computing_timer, "initial setup");
- fe1 =
+ fluid_fe =
std::make_unique<FESystem<spacedim>>(FE_Q<spacedim>(par.velocity_degree),
spacedim,
FE_Q<spacedim>(par.velocity_degree -
1);
- fe2 = std::make_unique<FE_Nothing<dim, spacedim>>();
- dh2.distribute_dofs(*fe2);
+ solid_fe = std::make_unique<FE_Nothing<dim, spacedim>>();
+ solid_dh.distribute_dofs(*solid_fe);
quadrature_formula = std::make_unique<QGauss<dim>>(par.velocity_degree + 1);
}
{
TimerOutput::Scope t(computing_timer, "setup dofs");
- dh1.distribute_dofs(*fe1);
+ fluid_dh.distribute_dofs(*fluid_fe);
std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
stokes_sub_blocks[dim] = 1;
- DoFRenumbering::component_wise(dh1, stokes_sub_blocks);
+ DoFRenumbering::component_wise(fluid_dh, stokes_sub_blocks);
auto dofs_per_block =
- DoFTools::count_dofs_per_fe_block(dh1, stokes_sub_blocks);
+ DoFTools::count_dofs_per_fe_block(fluid_dh, stokes_sub_blocks);
const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
- pcout << " Number of degrees of freedom: " << dh1.n_dofs() << " (" << n_u
- << '+' << n_p << " -- " << solid_particle_handler.n_global_particles()
- << '+' << tracer_particle_handler.n_global_particles() << ')'
- << std::endl;
+ pcout << " Number of degrees of freedom: " << fluid_dh.n_dofs() << " ("
+ << n_u << '+' << n_p << " -- "
+ << solid_particle_handler.n_global_particles() << '+'
+ << tracer_particle_handler.n_global_particles() << ')' << std::endl;
- owned1.resize(2);
- owned1[0] = dh1.locally_owned_dofs().get_view(0, n_u);
- owned1[1] = dh1.locally_owned_dofs().get_view(n_u, n_u + n_p);
+ fluid_owned_dofs.resize(2);
+ fluid_owned_dofs[0] = fluid_dh.locally_owned_dofs().get_view(0, n_u);
+ fluid_owned_dofs[1] =
+ fluid_dh.locally_owned_dofs().get_view(n_u, n_u + n_p);
IndexSet locally_relevant_dofs;
- DoFTools::extract_locally_relevant_dofs(dh1, locally_relevant_dofs);
- relevant1.resize(2);
- relevant1[0] = locally_relevant_dofs.get_view(0, n_u);
- relevant1[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p);
+ DoFTools::extract_locally_relevant_dofs(fluid_dh, locally_relevant_dofs);
+ fluid_relevant_dofs.resize(2);
+ fluid_relevant_dofs[0] = locally_relevant_dofs.get_view(0, n_u);
+ fluid_relevant_dofs[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p);
{
constraints.reinit(locally_relevant_dofs);
FEValuesExtractors::Vector velocities(0);
- DoFTools::make_hanging_node_constraints(dh1, constraints);
- VectorTools::interpolate_boundary_values(dh1,
- 0,
- ZeroFunction<spacedim>(spacedim +
- 1),
- constraints,
- fe1->component_mask(velocities));
+ DoFTools::make_hanging_node_constraints(fluid_dh, constraints);
+ VectorTools::interpolate_boundary_values(
+ fluid_dh,
+ 0,
+ ZeroFunction<spacedim>(spacedim + 1),
+ constraints,
+ fluid_fe->component_mask(velocities));
constraints.close();
}
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
- DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false);
+ DoFTools::make_sparsity_pattern(
+ fluid_dh, coupling, dsp, constraints, false);
SparsityTools::distribute_sparsity_pattern(
dsp,
- dh1.compute_locally_owned_dofs_per_processor(),
+ fluid_dh.compute_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
- system_matrix.reinit(owned1, dsp, mpi_communicator);
+ system_matrix.reinit(fluid_owned_dofs, dsp, mpi_communicator);
}
{
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
- DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false);
+ DoFTools::make_sparsity_pattern(
+ fluid_dh, coupling, dsp, constraints, false);
SparsityTools::distribute_sparsity_pattern(
dsp,
- dh1.compute_locally_owned_dofs_per_processor(),
+ fluid_dh.compute_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
- preconditioner_matrix.reinit(owned1, dsp, mpi_communicator);
+ preconditioner_matrix.reinit(fluid_owned_dofs, dsp, mpi_communicator);
}
- locally_relevant_solution.reinit(owned1, relevant1, mpi_communicator);
- system_rhs.reinit(owned1, mpi_communicator);
- solution.reinit(owned1, mpi_communicator);
+ locally_relevant_solution.reinit(fluid_owned_dofs,
+ fluid_relevant_dofs,
+ mpi_communicator);
+ system_rhs.reinit(fluid_owned_dofs, mpi_communicator);
+ solution.reinit(fluid_owned_dofs, mpi_communicator);
}
TimerOutput::Scope t(computing_timer, "Stokes_assembly");
- FEValues<spacedim> fe_values(*fe1,
+ FEValues<spacedim> fe_values(*fluid_fe,
*quadrature_formula,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values);
- const unsigned int dofs_per_cell = fe1->dofs_per_cell;
+ const unsigned int dofs_per_cell = fluid_fe->dofs_per_cell;
const unsigned int n_q_points = quadrature_formula->size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(spacedim);
- for (const auto &cell : dh1.active_cell_iterators())
+ for (const auto &cell : fluid_dh.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_matrix = 0;
}
const unsigned int component_i =
- fe1->system_to_component_index(i).first;
+ fluid_fe->system_to_component_index(i).first;
cell_rhs(i) += fe_values.shape_value(i, q) *
rhs_values[q](component_i) * fe_values.JxW(q);
}
SolidVelocity<spacedim> solid_velocity(par.angular_velocity);
- std::vector<types::global_dof_index> dof_indices1(fe1->dofs_per_cell);
+ std::vector<types::global_dof_index> dof_indices1(fluid_fe->dofs_per_cell);
- FullMatrix<double> local_matrix(fe1->dofs_per_cell, fe1->dofs_per_cell);
- dealii::Vector<double> local_rhs(fe1->dofs_per_cell);
+ FullMatrix<double> local_matrix(fluid_fe->dofs_per_cell,
+ fluid_fe->dofs_per_cell);
+ dealii::Vector<double> local_rhs(fluid_fe->dofs_per_cell);
auto particle = solid_particle_handler.begin();
while (particle != solid_particle_handler.end())
{
local_matrix = 0;
local_rhs = 0;
- const auto &cell = particle->get_surrounding_cell(tria1);
+ const auto &cell = particle->get_surrounding_cell(fluid_tria);
const auto &dh_cell =
- typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &dh1);
+ typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &fluid_dh);
dh_cell->get_dof_indices(dof_indices1);
const auto pic = solid_particle_handler.particles_in_cell(cell);
const auto real_q = p.get_location();
const auto properties = p.get_properties();
const auto &JxW = properties[0];
- for (unsigned int i = 0; i < fe1->dofs_per_cell; ++i)
+ for (unsigned int i = 0; i < fluid_fe->dofs_per_cell; ++i)
{
- const auto comp_i = fe1->system_to_component_index(i).first;
+ const auto comp_i =
+ fluid_fe->system_to_component_index(i).first;
if (comp_i < spacedim)
{
- for (unsigned int j = 0; j < fe1->dofs_per_cell; ++j)
+ for (unsigned int j = 0; j < fluid_fe->dofs_per_cell; ++j)
{
const auto comp_j =
- fe1->system_to_component_index(j).first;
+ fluid_fe->system_to_component_index(j).first;
if (comp_i == comp_j)
local_matrix(i, j) +=
- par.penalty_term * fe1->shape_value(i, ref_q) *
- fe1->shape_value(j, ref_q) * JxW;
+ par.penalty_term * fluid_fe->shape_value(i, ref_q) *
+ fluid_fe->shape_value(j, ref_q) * JxW;
}
local_rhs(i) += par.penalty_term *
solid_velocity.value(real_q, comp_i) *
- fe1->shape_value(i, ref_q) * JxW;
+ fluid_fe->shape_value(i, ref_q) * JxW;
}
}
}
locally_relevant_solution = solution;
const double mean_pressure =
- VectorTools::compute_mean_value(dh1,
+ VectorTools::compute_mean_value(fluid_dh,
QGauss<spacedim>(par.velocity_degree + 2),
locally_relevant_solution,
spacedim);
template <int dim, int spacedim>
- void StokesImmersedProblem<dim, spacedim>::refine_grid()
+ void StokesImmersedProblem<dim, spacedim>::refine_and_transfer()
{
TimerOutput::Scope t(computing_timer, "refine");
const FEValuesExtractors::Vector velocity(0);
- Vector<float> error_per_cell(tria1.n_active_cells());
- KellyErrorEstimator<dim>::estimate(dh1,
+ Vector<float> error_per_cell(fluid_tria.n_active_cells());
+ KellyErrorEstimator<dim>::estimate(fluid_dh,
QGauss<dim - 1>(par.velocity_degree + 1),
{},
locally_relevant_solution,
error_per_cell,
- fe1->component_mask(velocity));
+ fluid_fe->component_mask(velocity));
if (par.refinement_strategy == "fixed_fraction")
{
parallel::distributed::GridRefinement::
- refine_and_coarsen_fixed_fraction(tria1,
+ refine_and_coarsen_fixed_fraction(fluid_tria,
error_per_cell,
par.refinement_fraction,
par.coarsening_fraction);
else if (par.refinement_strategy == "fixed_number")
{
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
- tria1,
+ fluid_tria,
error_per_cell,
par.refinement_fraction,
par.coarsening_fraction,
par.max_cells);
}
- for (const auto &cell : tria1.active_cell_iterators())
+ for (const auto &cell : fluid_tria.active_cell_iterators())
if (cell->refine_flag_set() && cell->level() == par.max_level_refinement)
cell->clear_refine_flag();
parallel::distributed::SolutionTransfer<dim, LA::MPI::BlockVector> transfer(
- dh1);
- tria1.prepare_coarsening_and_refinement();
+ fluid_dh);
+ fluid_tria.prepare_coarsening_and_refinement();
transfer.prepare_for_coarsening_and_refinement(locally_relevant_solution);
- tria1.execute_coarsening_and_refinement();
+ fluid_tria.execute_coarsening_and_refinement();
setup_dofs();
transfer.interpolate(solution);
constraints.distribute(solution);
DataComponentInterpretation::component_is_scalar);
DataOut<spacedim> data_out;
- data_out.attach_dof_handler(dh1);
+ data_out.attach_dof_handler(fluid_dh);
data_out.add_data_vector(locally_relevant_solution,
solution_names,
DataOut<spacedim>::type_dof_data,
data_component_interpretation);
LA::MPI::BlockVector interpolated;
- interpolated.reinit(owned1, MPI_COMM_WORLD);
- VectorTools::interpolate(dh1,
+ interpolated.reinit(fluid_owned_dofs, MPI_COMM_WORLD);
+ VectorTools::interpolate(fluid_dh,
ConstantFunction<spacedim>(1.0, spacedim + 1),
interpolated);
- LA::MPI::BlockVector interpolated_relevant(owned1,
- relevant1,
+ LA::MPI::BlockVector interpolated_relevant(fluid_owned_dofs,
+ fluid_relevant_dofs,
MPI_COMM_WORLD);
interpolated_relevant = interpolated;
{
}
- Vector<float> subdomain(tria1.n_active_cells());
+ Vector<float> subdomain(fluid_tria.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
- subdomain(i) = tria1.locally_owned_subdomain();
+ subdomain(i) = fluid_tria.locally_owned_subdomain();
data_out.add_data_vector(subdomain, "subdomain");
data_out.build_patches();
}
{
TimerOutput::Scope t(computing_timer, "Set tracer particle motion");
- interpolate_field_on_particles(dh1,
+ interpolate_field_on_particles(fluid_dh,
tracer_particle_handler,
locally_relevant_solution,
tracer_particle_velocities,
assemble_nitche_restriction();
solve();
- if (cycle % par.mod_output == 0)
+ if (cycle % par.output_frequency == 0)
{
static unsigned int output_cycle = 0;
output_results(output_cycle, time);
}
++output_cycle;
}
- if (cycle % par.mod_refinement == 0 &&
+ if (cycle % par.refinement_frequency == 0 &&
cycle != par.number_of_time_steps - 1)
- refine_grid();
+ refine_and_transfer();
}
}
+
+ template <int dim, int spacedim>
+ StokesImmersedProblemParameters<dim,
+ spacedim>::StokesImmersedProblemParameters()
+ : ParameterAcceptor("Stokes Immersed Problem/")
+ , rhs("Right hand side", spacedim + 1)
+ , angular_velocity("Angular velocity", spacedim == 3 ? spacedim : 1)
+ {
+ // We split the parameters in various cathegories, by putting them in
+ // different sections of the ParameterHandler class. We begin by declaring
+ // all the global parameters used by StokesImmersedProblem in the global
+ // scope:
+ add_parameter(
+ "Velocity degree", velocity_degree, "", this->prm, Patterns::Integer(1));
+
+ add_parameter("Number of time steps", number_of_time_steps);
+ add_parameter("Output frequency", output_frequency);
+
+ add_parameter("Final time", final_time);
+
+ add_parameter("Viscosity", viscosity);
+
+ add_parameter("Nitsche penalty term", penalty_term);
+
+ add_parameter("Initial fluid refinement",
+ initial_fluid_refinement,
+ "Initial mesh refinement used for the fluid domain Omega");
+
+ add_parameter("Initial solid refinement",
+ initial_solid_refinement,
+ "Initial mesh refinement used for the solid domain Gamma");
+
+ add_parameter(
+ "Particle insertion refinement",
+ particle_insertion_refinement,
+ "Refinement of the volumetric mesh used to insert the particles");
+
+ add_parameter(
+ "Homogeneous Dirichlet boundary ids",
+ homogeneous_dirichlet_ids,
+ "Boundary Ids over which homogeneous Dirichlet boundary conditions are applied");
+
+ // Next section is dedicated to the parameters used to create the various
+ // grids. We will need three different triangulations: `Grid one` is used
+ // to define the fluid domain, `Grid two` defines the solid domain, and
+ // `Particle grid` is used to distribute some tracer particles, that are
+ // advected with the velocity and only used as passive tracers.
+ enter_my_subsection(this->prm);
+ this->prm.enter_subsection("Grid generation");
+ this->prm.add_parameter("Grid one generator", name_of_grid1);
+ this->prm.add_parameter("Grid one generator arguments",
+ arguments_for_grid1);
+
+ this->prm.add_parameter("Grid two generator", name_of_grid2);
+ this->prm.add_parameter("Grid two generator arguments",
+ arguments_for_grid2);
+
+ this->prm.add_parameter("Particle grid generator", name_of_particle_grid);
+ this->prm.add_parameter("Particle grid generator arguments",
+ arguments_for_particle_grid);
+ this->prm.leave_subsection();
+
+ leave_my_subsection(this->prm);
+
+
+
+ enter_my_subsection(this->prm);
+ this->prm.enter_subsection("Refinement and remeshing");
+ this->prm.add_parameter("Refinement step frequency", refinement_frequency);
+ this->prm.add_parameter("Refinement maximal level", max_level_refinement);
+ this->prm.add_parameter("Refinement strategy",
+ refinement_strategy,
+ "",
+ Patterns::Selection("fixed_fraction|fixed_number"));
+ this->prm.add_parameter("Refinement coarsening fraction",
+ coarsening_fraction);
+ this->prm.add_parameter("Refinement fraction", refinement_fraction);
+ this->prm.add_parameter("Maximum number of cells", max_cells);
+
+ this->prm.leave_subsection();
+ leave_my_subsection(this->prm);
+
+ // correct the default dimension for the functions
+ rhs.declare_parameters_call_back.connect([&]() {
+ Functions::ParsedFunction<spacedim>::declare_parameters(this->prm,
+ spacedim + 1);
+ });
+ angular_velocity.declare_parameters_call_back.connect([&]() {
+ Functions::ParsedFunction<spacedim>::declare_parameters(
+ this->prm, spacedim == 3 ? spacedim : 1);
+ });
+ }
+
} // namespace Step70
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
StokesImmersedProblemParameters<2> par;
- par.declare_all_parameters();
- std::ofstream out("default.prm");
- par.prm.print_parameters(out, ParameterHandler::ShortText);
ParameterAcceptor::initialize("parameters.prm", "used_parameters.prm");
StokesImmersedProblem<2> problem(par);
projects / dealii.git / commitdiff