Hi Daniel

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 
> <https://www.dealii.org/8.5.0/doxygen/deal.II/namespaceVectorTools.html#a05db6c8cebf924b417dd92f525efe3db>.

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 
> <https://www.dealii.org/8.5.0/doxygen/deal.II/classFiniteElement.html#a27220a135402b96c7e6eecbb04acda56>
> .
> 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
Quadrature<dim> dummy_q(unit_points.size()); 
MappingQGeneric<dim> mapping(1);
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, 
FE_Q<dim>(degree), 1) 

*       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!

Thank you

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 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Reply via email to