Hi Jorge, I have encountered the same hurdle. In my problem, the "periodic" dimension is the time, not some space dimension. Quite obviously, deal.II has not been designed for such setups, and I have found at least three theoretical reasons why a problem with time dependent variables can be solved faster when converted into an initial-value problem which deal.II supports directly. However, I like hacking my own way through. As a consequence, I have some experience to share here.
Things would be much easier if deal.ii's Triangulation class supported the notion of (topologically) identical vertices with different coordinates. Otherwise, combining periodic boundary conditions with mesh refinement seems to require significant amounts of code on the application side. I will call the periodic dimension "time" here; please feel free to generalize. My approach is inductive: Start with equally discretized time boundaries, then maintain equality of time-boundary discretizations across refinements. I use prepare_coarsening_and_refinement for grid smoothing and then fix any discrepancy in refine/coarsen requests between the time boundaries. This is repeated until the refine/coarsen flags do not change any more. The fix needs a conflict resolution strategy. My approach is to favour local refinements over local non-refinements. Each time-boundary cell is compared with its corresponding cell from the other time boundary, and when the refine/coarsen requests differ between the two cells, the requests are set to the finer version. Caveat 1: You have to avoid the possibility of endless loops where the mesh smoother tries to un-refine some cell that you have fixed for refinement. Therefore, I have disabled the option eliminate_refined_boundary_islands when initializing the triangulation object. I have attached code snippets, essentially two member function bodies, along with auxiliary stuff. The code presumes that you have "colorized" the boundaries such that time boundaries have recognizable boundary_indicator() values different from those for space boundaries. You will have to adapt the code, of course. (Or just glance at it, and then write your own. I did not clean it up much.) Caveat 2: You may want (or, depending on PDE type, you must avoid) to refine time and space independently. Therefore, the code snippets have an (untested) option to work with anisotropic refinement. If it is enabled, the fixes maintaining equally discretized time boundaries will refine boundary cells only along the respective time boundary, i.e. subdivision will be applied to the concerned space dimensions, but never to the time dimension. Otherwise, time-boundary cells may be refined even along their time dimension, in order to preserve refinement isotropy. Regards, -- Christian Cornelssen
// C++ code snippets for dealing with periodic boundary conditions and
// local grid coarsening or refinement.
// Written 2008, 2009 by Christian Cornelssen.
// Published under BSD-like license terms.
// deal.II includes
#include <base/function.h>
#include <base/geometry_info.h>
#include <base/smartpointer.h>
#include <lac/block_sparse_matrix.h>
#include <lac/vector.h>
#include <grid/tria.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/grid_refinement.h>
#include <fe/fe_system.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_constraints.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_renumbering.h>
#include <dofs/dof_tools.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <numerics/matrices.h>
#include <numerics/solution_transfer.h>
#include <numerics/vectors.h>
// Standard C++ includes
#include <cassert> // assert
#include <cmath> // std::ldexp
#include <iostream> // std::cout, std::cerr, std::ostream::operator<<
#include <map> // std::map, std::multimap
#include <vector> // std::vector
// We need a lexicographic comparison operator for points
// in order to build multimaps with points as keys
namespace dealii
{
template< int dim >
bool operator< (const Point<dim> &a, const Point<dim> &b)
{
for (unsigned int d = 0; d < dim; ++d)
{
if (a[d] < b[d]) return true;
if (b[d] < a[d]) return false;
}
return false;
}
}
/*
* Lots of class definitions for modelling the physical problem (not the
* approach for solving) have been removed. The remaining code occasionally
* references those removed definitions, but should be comprehensible
* nevertheless. Just note that
* n_coords == space_dim + 1 (for the time dimension);
* typedef Tensor<1,space_dim> Space; typedef Point<space_dim> SpacePoint;
* typedef Tensor<1,n_coords> SpaceTime; typedef Point<n_coords> SpaceTimePoint;
* and some auxiliary functions convert between those types.
*/
// Model of the approach for the solution
template< int n_coords >
class RectFlowT : public Problem<n_coords>
{
public:
// Members of base class templates must be referenced as such explicitly.
// Caution: When redeclaring, current access specifications apply.
// Using "typedef typename" because "using typename" does not help with g++.
using SpaceTimeInfo<n_coords>::space_dim;
typedef typename SpaceTimeInfo<n_coords>::Coord Coord;
using SpaceTimeInfo<n_coords>::space_coord_first;
using SpaceTimeInfo<n_coords>::space_coord_last;
using SpaceTimeInfo<n_coords>::time_coord;
typedef typename SpaceTimeInfo<n_coords>::Space Space;
typedef typename SpaceTimeInfo<n_coords>::SpacePoint SpacePoint;
typedef typename SpaceTimeInfo<n_coords>::SpaceTime SpaceTime;
typedef typename SpaceTimeInfo<n_coords>::SpaceTimePoint SpaceTimePoint;
using SpaceTimeInfo<n_coords>::space;
using SpaceTimeInfo<n_coords>::time;
using SpaceTimeInfo<n_coords>::separate_space_time;
using SpaceTimeInfo<n_coords>::spacetime;
using Given<n_coords>::domain;
using Given<n_coords>::fluid;
// Constructor
RectFlowT (const Problem<n_coords> &problem);
~RectFlowT ();
void run ();
private:
// ... Off-topic stuff removed ...
class HashSupportInfo; // hash function (Point<n_coords>, component)
// ... Off-topic member functions removed ...
void allocate_structures_for_dof_couplings ();
// ... Off-topic member functions removed ...
void refine_grid ();
// ... Off-topic member functions removed ...
dealii::Triangulation<n_coords> grid;
dealii::FESystem<n_coords> fe;
dealii::DoFHandler<n_coords> dof_handler;
std::vector<unsigned int> dofs_per_block;
dealii::ConstraintMatrix constraints;
dealii::BlockSparsityPattern system_matrix_structure;
dealii::BlockSparseMatrix<double> system_matrix;
dealii::BlockVector<double> system_rhs;
// Solution vector and its previous iterate
dealii::BlockVector<double> solution, old_solution;
};
template< int n_coords >
RectFlowT<n_coords>::RectFlowT (const Problem<n_coords> &problem) :
Given<n_coords> (problem),
Problem<n_coords> (problem),
// eliminate_refined_boundary_islands could lead to endless loops in
// refine_grid() because domain wrapping might add such islands again.
grid
( static_cast<typename dealii::Triangulation<n_coords>::MeshSmoothing>
( dealii::Triangulation<n_coords>::maximum_smoothing &~
dealii::Triangulation<n_coords>::eliminate_refined_boundary_islands )
),
fe( /* ... irrelevant FE setup args removed ... */ ),
dof_handler (grid),
dofs_per_block (fe.n_blocks())
{}
// Function object that hashes its arguments into one double value
template< int n_coords >
class RectFlowT<n_coords>::HashSupportInfo : public dealii::Function<n_coords>
{
protected:
// Roughly the number of bits taken from each significant part.
// Negative value makes the component index have least significance.
static const int precision = -16;
const Problem<n_coords> *problem;
// Protected default constructor, unusable for outsiders
HashSupportInfo () {};
// The function called for every (support point, component index)
double value_for_single_support_info
(const dealii::Point<n_coords> &p, const unsigned int component) const
{
double hash = double (component) / double (n_components - 1);
for (int s = 0; s < space_dim; ++s)
{
const unsigned int st = space_coord_first + s;
hash = std::ldexp (hash, precision) +
p[st] / problem->domain.spacetime_box_length()[st];
}
// time components are not hashed.
return hash;
}
public:
HashSupportInfo (const Problem<n_coords> *problem) :
dealii::Function<n_coords> (RectFlowT<n_coords>::n_components),
problem (problem)
{};
/*
* dealii::Function interface
*/
virtual double value
(const dealii::Point<n_coords> &p, const unsigned int component = 0) const
{
return value_for_single_support_info (p, component);
}
virtual void vector_value
(const dealii::Point<n_coords> &p, dealii::Vector<double> &values) const
{
for (unsigned int c = 0; c < n_components; ++c)
{
values(c) = value_for_single_support_info (p, c);
}
}
virtual void value_list
( const std::vector< dealii::Point<n_coords> > &points,
std::vector<double> &values, const unsigned int component = 0 ) const
{
for (unsigned int q = 0; q < points.size(); ++q)
{
values[q] = value_for_single_support_info (points[q], component);
}
}
virtual void vector_value_list
( const std::vector< dealii::Point<n_coords> > &points,
std::vector< dealii::Vector<double> > &values ) const
{
for (unsigned int q = 0; q < points.size(); ++q)
{
for (unsigned int c = 0; c < n_components; ++c)
{
values[q](c) = value_for_single_support_info (points[q], c);
}
}
}
};
// Member function to be run after remeshing and before assemble_system()
// Requires: New mesh's dof_handler.distribute_dofs(fe) has been run
// Provides: Settings of data members: constraints, system_matrix_structure
template< int n_coords >
void RectFlowT<n_coords>::allocate_structures_for_dof_couplings ()
{
constraints.clear ();
// Setup constraints due to hanging nodes
dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints);
// Setup constraints due to periodic boundary conditions
{
// Find time-boundary nodes and the associated DoF indices.
// Trick: Use a library function that associates boundary DoFs with values.
// Thus, map each (location,component) tuple to a hash value.
typedef std::map<unsigned int,double> DofToSupport;
typedef std::multimap<double,unsigned int> SupportToDoF;
DofToSupport dof_to_support;
{
HashSupportInfo hash_support_info (this);
typename dealii::FunctionMap<n_coords>::type
hash_support_info_at_time_boundaries;
hash_support_info_at_time_boundaries[2*time_coord+1] =
hash_support_info_at_time_boundaries[2*time_coord] =
&hash_support_info;
dealii::VectorTools::interpolate_boundary_values( domain.mapping(),
dof_handler, hash_support_info_at_time_boundaries, dof_to_support );
}
// Reverse the map. This should yield two DoFs per support item.
SupportToDoF support_to_dofs;
for ( DofToSupport::iterator
i = dof_to_support.begin(), e = dof_to_support.end();
i != e; dof_to_support.erase (i++) )
{
support_to_dofs.insert (SupportToDoF::value_type (i->second, i->first));
}
// Equate the DoFs with equal support data (space coords, component).
// We use the fact that multimaps are traversed in non-descending key order.
for ( SupportToDoF::const_iterator
i = support_to_dofs.begin(), e = support_to_dofs.end(), j;
i != e; ++(i = j) )
{
// Check that there are exactly two distinct DoFs per hash value.
// Otherwise the hash function would need improvement.
assert (i == support_to_dofs.begin() || j->first != i->first);
++(j = i);
assert (j != e);
assert
( std::abs (j->first - i->first) <= std::ldexp
(std::max (std::abs (i->first), std::abs (j->first)), -40) );
unsigned int dof_i = i->second;
unsigned int dof_j = j->second;
assert (dof_i != dof_j);
std::cerr << ":"; // for debugging
// Equal discretization of time boundaries implies equal constraints.
assert
( constraints.is_constrained (dof_i) ==
constraints.is_constrained (dof_j) );
// At hanging nodes, the DoFs may be constrained already. Skipping.
if (constraints.is_constrained (dof_j)) continue;
// Equate the DoFs, thus enforcing periodic boundary conditions.
constraints.add_line (dof_j);
constraints.add_entry (dof_j, dof_i, 1.0);
}
std::cerr << std::endl;
}
constraints.close ();
dealii::BlockCompressedSparsityPattern
temp_system_matrix_structure
(dofs_per_block, dofs_per_block);
dealii::DoFTools::make_sparsity_pattern
(dof_handler, temp_system_matrix_structure, constraints);
system_matrix.clear();
system_matrix_structure.copy_from (temp_system_matrix_structure);
}
// Refine grid, update dof_handler, solution, old_solution
template< int n_coords >
void RectFlowT<n_coords>::refine_grid ()
{
// Refinement parameters
const bool anisotropic_boundary_refinement = false;
const unsigned int max_children_per_cell =
dealii::GeometryInfo<n_coords>::
max_children_per_cell;
const double desired_cell_count_factor = 2.0;
const double coarsen_fraction = 0.01 * max_children_per_cell;
const double refine_fraction = std::min
( 1.0 - coarsen_fraction,
coarsen_fraction / max_children_per_cell +
(desired_cell_count_factor - 1.0) /
(max_children_per_cell - 1) );
// For accurate implementation of periodic time-boundary conditions,
// the time boundaries must be equally discretized in space.
// To this end, first find out which cells are at time boundaries
// and group them in pairs with equal space locations (and distinct time
// coords). A suitable container is a multimap from space to cell iterators.
typedef typename dealii::DoFHandler<n_coords>::active_cell_iterator
CellIterator;
typedef typename dealii::DoFHandler<n_coords>::active_face_iterator
FaceIterator;
typedef std::multimap<SpacePoint,CellIterator>
LocationToCells;
LocationToCells location_to_cells;
for ( CellIterator
cell = dof_handler.begin_active(), endc = dof_handler.end();
cell != endc; ++cell )
{
if (!cell->at_boundary()) continue;
for ( unsigned int face_no = 0;
face_no < dealii::GeometryInfo<n_coords>::faces_per_cell;
++face_no )
{
FaceIterator face = cell->face (face_no);
unsigned char b = face->boundary_indicator();
// Examine time boundary face(s) only.
// Note: If the time axis is not subdivided,
// there are two (opposite) time boundary faces per cell.
if (!(b == 2*time_coord || b == 2*time_coord+1)) continue;
location_to_cells.insert
( typename LocationToCells::value_type
(space (face->center()), cell) );
}
}
// Check that each space location maps to exactly two different cells
// (one per time boundary).
// Otherwise there are numeric errors, or the space discretizations
// of the time boundaries are different even before refinement.
for ( typename LocationToCells::const_iterator
i = location_to_cells.begin(), e = location_to_cells.end(), j;
i != e; ++(i = j) )
{
assert (i == location_to_cells.begin() || j->first != i->first);
++(j = i);
assert (j != e);
assert
( (j->first - i->first).square() <= std::ldexp
(std::max ((i->first).square(), (j->first).square()), -40) );
// Actually, the corresponding cell(s) may be identical as long as
// they have two separate time boundary faces.
// For simpler code checking, that sophistication is neglected here.
// If you deal with a non-subdivided time axis, remove the following
// check and rely on the buildup procedure for location_to_cells
// to enter a cell only once per time-boundary face.
assert (j->second != i->second);
}
// Estimate the approximation error for each cell in the domain
dealii::Vector<float> estimated_error_per_cell
(grid.n_active_cells());
// ... Off-topic code removed ...
// Refine and coarsen the grid according to the estimated error distribution
dealii::GridRefinement::refine_and_coarsen_fixed_number
(grid, estimated_error_per_cell, refine_fraction, coarsen_fraction);
// Postprocess the grid to preserve mesh smoothness
// and enable maintenance of periodic boundary conditions
bool grid_smoothed = false;
bool domain_wrapped = false;
while (!(grid_smoothed && domain_wrapped))
{
// Ensure that the refine/coarsen flag settings respect grid smoothness
if (!grid_smoothed)
{
if (grid.prepare_coarsening_and_refinement ()) // actually smoothed
{
domain_wrapped = false;
}
grid_smoothed = true;
}
// If there are refined or coarsened cells at a time boundary, make sure
// that the corresponding cells at the other time boundary have the same
// refine/coarsen status, except perhaps along the time axis.
// Note: exclusion of the time axis may introduce anisotropic refinement.
if (!domain_wrapped)
{
for ( typename LocationToCells::const_iterator
i = location_to_cells.begin(), e = location_to_cells.end(), j;
i != e; ++(i = j) )
{
++(j = i);
CellIterator cell_i = i->second,
cell_j = j->second;
// In case of discrepancies in refine/coarsen states, decide
// in favor of refinement. Otherwise, grid smoothing and domain
// wrapping could undo each other's work and thus cause endless
// postprocessing.
if (cell_i->coarsen_flag_set() != cell_j->coarsen_flag_set())
{
cell_i->clear_coarsen_flag ();
cell_j->clear_coarsen_flag ();
grid_smoothed = false;
}
const dealii::RefinementCase<n_coords>
refine_i = cell_i->refine_flag_set(),
refine_j = cell_j->refine_flag_set(),
refine_mask =
anisotropic_boundary_refinement ?
~ dealii::RefinementCase<n_coords>::
cut_axis (time_coord) :
dealii::RefinementPossibilities<n_coords>::
isotropic_refinement;
if ((refine_i & refine_mask) != (refine_j & refine_mask))
{
const dealii::RefinementCase<n_coords>
refine = (refine_i | refine_j) & refine_mask;
cell_i->set_refine_flag (refine_i | refine);
cell_j->set_refine_flag (refine_j | refine);
grid_smoothed = false;
}
}
domain_wrapped = true;
}
}
// Carry over the solution from the old grid to the new grid
dealii::SolutionTransfer< n_coords, dealii::BlockVector<double> >
solution_field (dof_handler);
solution_field.prepare_for_coarsening_and_refinement (solution);
grid.execute_coarsening_and_refinement ();
dof_handler.distribute_dofs (fe);
dealii::DoFRenumbering::component_wise (dof_handler);
// dealii::DoFRenumbering::Cuthill_McKee (dof_handler);
dealii::DoFTools::count_dofs_per_block (dof_handler, dofs_per_block);
old_solution.reinit (dofs_per_block);
solution_field.interpolate (solution, old_solution);
solution = old_solution;
}
signature.asc
Description: PGP signature
_______________________________________________
