Generalize identity determination in FE_SimplexP #14607
Conversation
| // dofs are located along lines, so two dofs are identical if they are | ||
| // located at identical positions. | ||
| // Therefore, read the points in unit_support_points for the | ||
| // first coordinate direction. For FE_SimplexP, they are currently | ||
| // hard-coded and we iterate over points on the first line which begin | ||
| // after the 3 vertex points in the complete list of unit support points | ||
| std::vector<std::pair<unsigned int, unsigned int>> identities; | ||
|
|
||
| for (unsigned int i = 0; i < this->degree - 1; ++i) | ||
| for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j) | ||
| if (std::fabs(this->unit_support_points[i + this->reference_cell() | ||
| .n_vertices()][0] - | ||
| fe_p_other->unit_support_points | ||
| [i + this->reference_cell().n_vertices()][0]) < 1e-14) | ||
| identities.emplace_back(i, j); |
There was a problem hiding this comment.
This looks reasonable, but I think that the comment actually was helpful in understanding the code. Why did you remove it?
Separately, do we fill face_unit_support_points for FE_SimplexP? If so, it might be easier to use because then the indexing becomes simpler.
| else if (const auto *fe_p_other = | ||
| dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other)) | ||
| { | ||
| std::vector<std::pair<unsigned int, unsigned int>> identities; | ||
|
|
||
| for (unsigned int i = 0; i < this->degree - 1; ++i) | ||
| for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j) | ||
| if (std::fabs(this->unit_support_points[i + 3][0] - | ||
| fe_p_other->unit_support_points[i + 3][0]) < 1e-14) | ||
| if (std::fabs(this->unit_support_points[i + this->reference_cell() | ||
| .n_vertices()][0] - | ||
| fe_p_other->get_unit_support_points() | ||
| [i + this->reference_cell().n_vertices()][1]) < 1e-14) |
There was a problem hiding this comment.
This block could use a comment that says that if this function is called with both an FE_SimplexP and a FE_PyramidP, then we are clearly at a triangular face of the pyramid and can treat it the same way we did above for two FE_SimplexP meeting.
| else | ||
| { | ||
| // If nodes are not located in the same place, we have to | ||
| // interpolate. This is then not handled through the | ||
| // current function, but via interpolation matrices that | ||
| // result in constraints, rather than identities. Since | ||
| // that happens in a different function, there is nothing | ||
| // for us to do here. | ||
| } |
There was a problem hiding this comment.
I added this in #14444 because I was wondering how this is working if the if above is not satisfied. I think the comment is still relevant and if you don't mind too much, please leave it here and duplicate it in the FE_SimplexP case above.
| else if (const auto *fe_p_other = | ||
| dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other)) |
There was a problem hiding this comment.
I believe that you can avoid the code duplication in the same way as you did for the vertex case above: Just have a sequence of dynamic_cast<...> || dynamic_cast<...> questions in the if, and then access the information via fe_other rather than via fe_p_other. I believe that nothing you do in the if block actually requires access to information of the derived object.
tamiko
changed the title
|
@peterrum If you address the comments, I will gladly look at it again! |
@sebproell Maybe you could use this as starting point.