Thank you for your help!
The code below suggests that you want to constraint all degrees of freedom.
> Is this really what you intend? If so, it might be better to just
> interpolate the values using VectorTools::interpolate
Yes, I want to constrain the velocity of the fluid in a certain region. The
problem arises in my fluid-structure interaction code using immersed finite
element method where I want to interpolate the solid velocity to the fluid
grid which is overlapped with the solid as Dirichlet BCs. I find the
to-be-constrained region on the fly as the solid moves. Since the dofs to
be constrained are not on the actual boundary and are constantly changing,
I can't use VectorTools::interpolate.
> Don't rely on any particular order and use the whole finite element
> instead of just the first base element in your code! To find out if a
> degree of freedom corresponds to the velocity FiniteElement you can then
> call system_to_component_index
> The problem with the approach above is that there (probably?) is no
> DoFHandler that is responsible for the velocity degrees of freedoms alone.
> Hence, a cell iterator (based on a DoFHandler) can only give you the global
> degrees of freedoms
> for all the components at the same time. By using
> system_to_component_index you can then decide which to use.
Following the idea you described, I came up with the following code, but
something went wrong...
const vector<Point<dim>>& unit_points = fe.get_unit_support_points(); //
the size should be equal to dofs_per_cell
FEValues<dim> dummy_fe_values(mapping, fe, dummy_q,
std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell); // need
this to find the global dof index
for (auto cell = dof_handler.begin_active(); cell != dof_handler.end();
*cell->get_dof_indices(dof_indices); // querying all the global indicies
associated with this cell*
for (unsigned int i = 0; i < unit_points.size(); ++i)
// fe is defined as FESystem(FE_Q<dim>(degree+1), dim,
* auto index = fe.system_to_component_index(i); // find the component
number and base number of current dof i if (index.first < dim) // Is
this the correct method to tell if the dof i is a velocity dof rather than
a pressure dof??*
Point<dim> mapped_point = mapping.transform_unit_to_real_cell(cell,
unit_points[i]); // get the real coordinates
Vector<double> value = f(mapped_point); // Calculate (interpolate)
the constrained value which is a *vector*
// set the inhomogeneity constraints on the velocity
*auto line = dof_indices[i]; // I think this is the global dof number
corresponding to dof i*
*constraints.set_inhomogeneity(**line, value[index.first]); // the
component number should be the correct velocity component to constrain *
I think the logic should be good, but I might have misunderstood the return
value of system_to_component.. Please let me know if you can see the error.
> No, as far as I know we don't support writing the results as 2nd order
> elements. The approach we take instead is to split a cell into several
> (most of the time degree^dim) using DataOut::build_patches(degree) to
> obtain a linear interpolation of the correct order.
I can't believe I did not know DataOut::build_patches can do this trick!
This is exactly what I need!
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
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
For more options, visit https://groups.google.com/d/optout.