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.