Extend FESystem for more than a unique face #12647
Conversation
| if (this->n_unique_faces() > 1) | ||
| return; |
There was a problem hiding this comment.
Not ideal but we currently don't support hanging-nodes for meshes with wedges/pyramids anyway.
There was a problem hiding this comment.
Hm, yes, not ideal. How can we catch this error in a better way? Is there a better place where this assertion could be put?
|
|
||
| for (unsigned int i = 0; i < points.size(); ++i) | ||
| for (unsigned int i = 0; i < values.size(); ++i) |
There was a problem hiding this comment.
No, this isn't right. If the function is called with a different number of input and output arguments, then the place where the function is called needs to be fixed.
| /** | ||
| * GenericDoFsPerObject object this class was initialized with. | ||
| */ | ||
| const internal::GenericDoFsPerObject data; |
There was a problem hiding this comment.
First, can we give this variable a better name? Second, why is this necessary?
| fes[i]->n_dofs_per_quad(0) * multiplicities[i]; | ||
| multiplied_dofs_per_hex += | ||
| fes[i]->n_dofs_per_hex() * multiplicities[i]; | ||
| if (initialized == false) |
There was a problem hiding this comment.
If I read this right, what you are doing here is initialize the dpo object the first time you find an element with multiplicity>0. I would much rather you do this before the current for loop, as in
for (i=...)
if (multiplicity>0)
{
...do it here...
break;
}
I think that would be easier to read, and avoid the initialized variable altogether.
| ++j) | ||
| for (unsigned int k = 0; | ||
| k < fes[i]->data.dofs_per_object_exclusive[j].size(); | ||
| ++k) |
There was a problem hiding this comment.
It seems to me like you're repeating many of these loop statements. Could their bodies be merged?
| @@ -1140,6 +1140,8 @@ ReferenceCell::face_to_cell_lines(const unsigned int face, | |||
| AssertIndexRange(face, n_faces()); | |||
| AssertIndexRange(line, face_reference_cell(face).n_lines()); | |||
|
|
|||
| static unsigned int X = numbers::invalid_unsigned_int; | |||
There was a problem hiding this comment.
| static unsigned int X = numbers::invalid_unsigned_int; | |
| static constexpr unsigned int X = numbers::invalid_unsigned_int; |
| Assert(1 <= degree && degree <= 2, ExcNotImplemented()) | ||
|
|
||
| for (double z = 0.0; z <= 1.0; z += (degree == 1 ? 1.0 : 0.5)) | ||
| { |
| if (degree == 1) | ||
| Assert(1 <= degree && degree <= 2, ExcNotImplemented()) | ||
|
|
||
| for (double z = 0.0; z <= 1.0; z += (degree == 1 ? 1.0 : 0.5)) |
There was a problem hiding this comment.
| for (double z = 0.0; z <= 1.0; z += (degree == 1 ? 1.0 : 0.5)) | |
| for (double z = 0.0; z <= 1.0; z += 1./degree) |
This is going to run into trouble for values of degree for which 1./degree is not exactly representable in floating point arithmetic. The test z<=1 will then either work or not work, by accident. It's better to write this as
for (unsigned int layer=0; layer<=degree; ++layer)
{
const double delta_z = 1./degree;
const double z = layer * delta_z;
...
| this->unit_support_points.emplace_back(1.0, 0.0, z); | ||
| this->unit_support_points.emplace_back(0.0, 1.0, z); | ||
|
|
||
| if (degree == 2) |
There was a problem hiding this comment.
Maybe not for now, but if you wanted to leave a comment: These points are just the tensor product of the triangle points with z. We could simply create a FE_SimplexP<2>(degree) and ask for its support points; then you don't have to repeat the code for degree>2 and it has to be implemented in only one place.
|
|
||
| this->unit_face_support_points.resize(5); | ||
|
|
||
| for (unsigned int f = 0; f < 5; ++f) |
There was a problem hiding this comment.
same thing here with the bare 5
| template <int dim, int spacedim> | ||
| void | ||
| FE_WedgePoly<dim, spacedim>:: | ||
| convert_generalized_support_point_values_to_dof_values( |
There was a problem hiding this comment.
same here -- put it into a separate PR to make this one simpler.
...references #12646