I've found what was wrong in this piece of code. It turned out that
index.first < dim does not guarantee this dof is a velocity dof. Changing
it to* if (i < fe.base_element(0).dofs_per_cell && index.first < dim)*
solves the problem.
But there are two things that I can not understand:
1. using i < fe.base_element(0).dofs_per_cell to make sure the current dof
is a velocity dof which I want to constrain, but this is obviously not good
practice as it assumes the velocity dofs are always in front of pressure
dofs. What is the correct way?
2. index.first should return the component number: if the dof is a
velocity dof then it would be in range [0, dim), if it is a pressure dof
then it should be 0. However, based on my test index.first can be greater
than dim, why? Then how to get the velocity component number (0 for x, 1
for y)?
* 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 *
> }
>
Thanks
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.
