Good morning everyone,
I am attempting to write a program which projects a solution vector from
mesh A to mesh B. Both meshes are parallel distributed and have different
geometries (the first is adaptively refined, the latter has uniform grid
size). In addition, the solution uses Q1 elements on mesh A and is to be
piecewise constant on mesh B.
I have attempted to use VectorTools::project() for this purpose.
- - I create a child of the Function class which takes in a
solution vector and its corresponding mesh (both parallel distributed).
This function allows the solution vector to be evaluated at any coordinate
within the domain (and should also consider ghost elements).
- - I insert this as a parameter into the project() function.
This is projected onto a non-ghosted vector, which is then assigned to a
ghosted vector and is subsequently compressed.
The resulting vector correctly projects to some elements, however other
*elements
appear to be set to 0* (incorrectly). In addition, the number of elements
that are zero and their *locations vary with the number of processors used*.
I suspect that it may be something to do with project() not being able to
handle ghost elements but I’m not sure. It is very possible that the
Function used to evaluate the input vector at any coordinate is incorrect,
i have included it in a text file below. I have tried using both
FEFieldFunction() and later find_all_active_cells_around_point() but to no
avail (the latter appears to return empty for these zeroed elements, even
after I perturb the coordinate by a significant amount when no cell is
found or use a high tolerance).
I then attempted the following crude workaround:
- serialize the vector and triangulation
- project them to a new serial vector using project()
- scatter them back onto the parallel mesh
Unfortunately, this also failed. I am having difficulty copying the
parallel mesh onto a serial one, which appears far more involved than I
originally thought.
On a sidenote I have attempted to use interpolate() as well and was still
unsuccessful (I encountered the same issue with elements set to 0 at the
same locations for a given number of cores used).
Any help would be very, very much appreciated
Best wishes,
Daniel Rivlin
--
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.
To view this discussion visit
https://groups.google.com/d/msgid/dealii/1fface4c-15f1-4afd-b667-2153ccfcb8c1n%40googlegroups.com.
template <int dim>
class LinearSolutionFunction : public Function<dim>
{
public:
/**
* Constructor.
*
* @param dof_handler_in The DoFHandler corresponding to the linear finite
element space.
* @param solution_in The solution vector.
*/
LinearSolutionFunction(const DoFHandler<dim> &dof_handler_in,
const LA::MPI::Vector &solution_in)
: Function<dim>(1),
dof_handler_linear(dof_handler_in),
solution_linear(solution_in),
mapping_linear(1)
{}
/**
* Evaluate the solution at point p by:
* - Locating all active cells around p (within a given tolerance).
* - For each cell that is locally available (owned or ghost),
* compute the reference coordinate and use FEValues to perform linear
interpolation.
* - Return the average value from all valid cells.
*
* Returns 0.0 if no valid cell is found.
*/
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override
{
// Choose a tolerance that might need tuning.
double tolerance = 1e-6; // Try increasing tolerance for adaptive meshes
// Find all active cells around the point p within the specified
tolerance.
auto cell_and_refs =
GridTools::find_all_active_cells_around_point(mapping_linear,
dof_handler_linear,
p,
tolerance);
// If no cells found, optionally try a slight perturbation.
if (cell_and_refs.empty())
{
Point<dim> p_perturbed = p;
p_perturbed[0] += 1e-5;
p_perturbed[1] += 1e-5;
cell_and_refs =
GridTools::find_all_active_cells_around_point(mapping_linear,
dof_handler_linear,
p_perturbed,
tolerance);
if (cell_and_refs.empty())
return 0.0; // No cell found even after perturbation.
}
double result = 0.0;
unsigned int valid_cell_count = 0;
for (const auto &cell_and_ref : cell_and_refs)
{
const auto &cell = cell_and_ref.first;
// Only use cells that are locally available (owned or ghost).
if (!(cell->is_locally_owned() || cell->is_ghost()))
continue;
const auto &xi = cell_and_ref.second;
const FiniteElement<dim> &fe = cell->get_fe();
const unsigned int dofs_per_cell = fe.dofs_per_cell;
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
// Create a quadrature rule with the computed reference coordinate xi.
std::vector<Point<dim>> quadrature_points(1, xi);
Quadrature<dim> quadrature(quadrature_points);
FEValues<dim> fe_values(mapping_linear, fe, quadrature,
update_values);
fe_values.reinit(cell);
double local_value = 0.0;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
local_value += solution_linear[local_dof_indices[i]] *
fe_values.shape_value(i, 0);
}
result += local_value;
++valid_cell_count;
}
return (valid_cell_count > 0 ? result / valid_cell_count : 0.0);
}
private:
const DoFHandler<dim> &dof_handler_linear;
const LA::MPI::Vector &solution_linear;
MappingQ<dim> mapping_linear;
};
