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, update_quadrature_points); 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) { dummy_fe.reinit(cell); *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.add_line(line); *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 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.
