Dear Wolfgang
Thanks for your attention. The physical background is that I want to
interpolate fluid velocity (v^f) onto an arbitrary point (which is in fact
a solid support point). Instead of doing a standard finite element
interpolation, I want to do a "delta function interpolation":
v^s = \int_\Omega v^f \delta(x - x^s) d\Omega
where \delta is the "delta function" defined in the previous post.
I've tried to read your question a couple of times, but I can't seem to
> figure
> out what exactly you want to do. I understand the definition of the delta
> function located around a given point 'p', but what is (in mathematical
> formulas) what you want to do in the interpolation process? You have an
> input
> and an output vector, but I don't think I understand how they are related.
>
The logic of the Interpolator class is: for a specific target point,
iterate over all the nodes on the source grid, compute the distance between
the node and the target point, then the weight. Next, given a source
vector, the Interpolator does the interpolation using the
previously-computed weights.
1. However, there is no convenient way to iterate over all the "nodes" on a
grid and query their corresponding dofs. Therefore I have to iterate over
all cells, then iterate over all the support points and find which support
points are contributing to the target point.
2. In addition, I may want to interpolate more than one vector onto a
specific point, so I want to cache the influential nodes, compute their
weights, and reuse them.
3. The tuple is used to mark the "node" that is influential to the target
point. Intuitively one would use a "global node id" for that, but deal.II
doesn't work this way. So I use {cell_pointer, support_point_id, weight},
which essentially identifies a "node". While doing interpolation, I must be
careful that I don't use the same node for multiple times.
Yes, system_to_component_index() returns a std::pair, that was a bug.
What I did is as following:
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 tuples of cell iterator, supporting point id,
* and weight
* as there is no convenient way to express "global node id".
*/
std::vector<std::tuple<typename DoFHandler<dim>::active_cell_iterator,
unsigned int,
double>>
sources;
};
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 xbar = (unit_points[v][d] - target[d]) / h;
weight *= (std::abs(xbar) <= 2
? xbar * 0.25 * (1 + std::cos(M_PI * xbar / 2))
: 0);
}
// If weight is nonzero then it is an influential point
if (weight > 0)
{
sources.push_back(std::make_tuple(cell, v, weight));
}
}
}
}
template <int dim, typename VectorType>
void DiracDeltaInterpolator<dim, VectorType>::interpolate(
const VectorType &source_vector,
Vector<typename VectorType::value_type> &target_vector)
{
std::vector<bool> bs(dof_handler.n_dofs(), false);
target_vector = 0;
const FiniteElement<dim> &fe = dof_handler.get_fe();
std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
for (auto p : sources)
{
std::get<0>(p)->get_dof_indices(dof_indices);
// Vector component of this support point
auto d = fe.system_to_component_index(std::get<1>(p)).first;
AssertThrow(d < dim, ExcMessage("Vector component not less than
dim!"));
// Global dof index of this support point
auto index = dof_indices[std::get<1>(p)];
AssertThrow(index < dof_handler.n_dofs(),
ExcMessage("Wrong index of global dof!"));
if (!bs[index])
{
// Interpolate
target_vector[d] += std::get<2>(p) * source_vector[index];
bs[index] = true;
}
}
}
Note that I used a bool vector of size n_dofs() to make sure every
influential node is used only once. Besides, I used get_dof_indices and
support point id to query the dof value. Are these methods correct?
Thank you
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.
