Hi everyone
I am trying to implement a template Dirac delta function, which takes a
source vector and its associated DoFHandler, then interpolates the source
vector onto an arbitrary point. Different from VectorTools::point_value,
the target point does not have to be in the source grid. Instead, every
node on the source grid has an "influential domain", as long as the target
point is in this domain, the source node will contribute to the target
value. The Dirac delta function is defined as follows:
phi(r) = (|r| <= 2 ? 0.25 * (1 + cos(M_PI * r / 2) : 0);
delta(x) = phi(x1/h) * phi(x2/h) * phi(x3/h) / h^3
where x is the vector formed by the target point and source point, and h is
a given mesh size.
Because the same interpolation will be done for multiple vectors, I wrote a
class to find and cache all the "influential nodes" and corresponding
weights of a given target point, then does the interpolation using any
source vector. Below is my declaration:
template <int dim, typename VectorType>
class DiracDeltaInterpolator
{
public:
DiracDeltaInterpolator(const DoFHandler<dim> &, const Point<dim> &,
double);
void interpolate(const VectorType &,
Vector<typename VectorType::value_type> &);
private:
/// Source DoFHandler
const DoFHandler<dim> &dof_handler;
/// Target point to interpolate to
const Point<dim> ⌖
/// Background mesh size
double h;
/**
* Source nodes that contribute to the target point,
* denoted as pairs of cell iterator and supporting point id,
* as there is no convenient way to express "global node id".
*/
std::vector<std::pair<DoFHandler<dim>::active_cell_iterator,
unsigned int>> sources;
};
And the following is my implementation:
template <int dim, typename VectorType>
DiracDeltaInterpolator<dim, VectorType>::DiracDeltaInterpolator(
const DoFHandler<dim> &dof_handler, const Point<dim> &point, double h)
: dof_handler(dof_handler), target(point), h(h)
{
const FiniteElement<dim> &fe = dof_handler.get_fe();
const std::vector<Point<dim>> &unit_points =
fe.get_unit_support_points();
Quadrature<dim> dummy_q(unit_points.size());
MappingQGeneric<dim> mapping(1);
FEValues<dim> dummy_fe_values(
mapping, fe, dummy_q, update_quadrature_points);
std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
for (auto cell : dof_handler.active_cell_iterators())
{
dummy_fe_values.reinit(cell);
cell->get_dof_indices(dof_indices);
// Real coordinates of the current cell's support points
auto support_points = dummy_fe_values.get_quadrature_points();
for (unsigned int v = 0; v < unit_points.size(); ++v)
{
// Compute the weight of each support point
double weight = 1;
for (unsigned int d = 0; d < dim; ++d)
{
double r = unit_points[v][d] - target[d];
weight *= (std::abs(r) <= 2 ? 0.25 / h * (1 +
std::cos(M_PI * r / 2)) : 0 );
}
// If weight is nonzero then it is an influential point
if (weight > 0)
{
sources.push_back({cell, v});
}
}
}
}
*Question 1*: if a node is shared by multiple cells then there will be
duplicates in the above code. Is there a way to avoid this other than uing
a std::map with Point<dim> as keys?
*Question 2*: now I have pre-computed all the influential support points,
denoted as pairs of cell pointer and id, how to interpolate?
template <int dim, typename VectorType>
void DiracDeltaInterpolator<dim, VectorType>::interpolate(
const VectorType &source_vector,
Vector<typename VectorType::value_type> &target_vector)
{
}
This problem originates from Immersed Boundary Method, which is a bit ugly
but should still be doable in deal.II framework. Thanks for any suggestions!
Jie
--
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.
For more options, visit https://groups.google.com/d/optout.
