Fix normal pressure Neumann conditions on 2D solid lines - #2210
Conversation
Honor orthopressure and pseudo_orthopressure on plane solid lines. Scale pressure by thickness and add follower-load linearization.
d34ab8a to
7c64a2b
Compare
dragos-ana
left a comment
dragos-ana
left a comment
There was a problem hiding this comment.
Thanks for implementing this! I added some comments / suggestions
| void add_normal_pressure_load(const Core::Elements::ShapeFunctionsAndDerivatives<celltype>& shape, | ||
| const Core::LinAlg::Matrix<Core::FE::num_nodes(celltype), 2>& current_coordinates, | ||
| const double pressure_times_weight, const double reference_thickness, | ||
| const int num_dof_per_node, Core::LinAlg::SerialDenseVector& force, | ||
| Core::LinAlg::SerialDenseMatrix* load_linearization) | ||
| { | ||
| double radial_derivative = 0.0; | ||
| double axial_derivative = 0.0; | ||
| for (int node = 0; node < Core::FE::num_nodes(celltype); ++node) | ||
| { | ||
| radial_derivative += shape.derivatives(0, node) * current_coordinates(node, 0); | ||
| axial_derivative += shape.derivatives(0, node) * current_coordinates(node, 1); | ||
| } | ||
|
|
||
| const std::array<double, 2> normal_measure{ | ||
| reference_thickness * axial_derivative, -reference_thickness * radial_derivative}; | ||
| for (int node = 0; node < Core::FE::num_nodes(celltype); ++node) |
There was a problem hiding this comment.
The function implementation could be moved to the .cpp. Also, I think that it would benefit a lot from formulae added as comments for each step of the computation (along with doxygen for the parameters and so on ;) ).
Also: does it need to be exposed in the hpp, or could it be simply implemented as a helper function in an anonymous namespace in the corresponding cpp? As far as I see, it's only used within evaluate_normal_pressure, and within the unit test you added.
UPDATE: After looking over the unit test it's fine for me to keep it here, as it makes testing the functionality a bit more accessible as described in the respective comment.
| "A two-dimensional pressure condition requires at least two entries in ONOFF, VAL, and " | ||
| "FUNCT."); | ||
| FOUR_C_ASSERT_ALWAYS( | ||
| onoff[0] == 1 && onoff[1] == 0, "Normal pressure must be activated on the first dof only."); |
There was a problem hiding this comment.
I am not really acquainted with the orthopressure conditions, but this is the current convention right? This would mean that the function id and values should also be posed on the first dof, would it make sense to add these consistency conditions here as well?
| for (int gp = 0; gp < integration.num_points(); ++gp) | ||
| { | ||
| const auto xi = Core::Elements::evaluate_parameter_coordinate<celltype>(integration, gp); | ||
| Core::Elements::ElementNodes<celltype, 2> nodes{.coordinates = current_coordinates}; |
There was a problem hiding this comment.
This could be moved out of the GP loop.
| Core::Elements::ElementNodes<celltype, 2> nodes{.coordinates = current_coordinates}; | ||
| const auto shape = Core::Elements::evaluate_shape_functions_and_derivs<celltype>(xi, nodes); | ||
|
|
||
| Core::LinAlg::Matrix<2, 1> current_coordinate; |
There was a problem hiding this comment.
Same here, you could declare it out of the loop, and then update it via multiply_tn:
| Core::LinAlg::Matrix<2, 1> current_coordinate; | |
| current_coordinate.multiply_tn(1.0,current_coordinates, shape.values,0.0); |
| for (int gp = 0; gp < integration.num_points(); ++gp) | ||
| { | ||
| const auto xi = Core::Elements::evaluate_parameter_coordinate<celltype>(integration, gp); | ||
| const Core::Elements::ElementNodes<celltype, 2> nodes{.coordinates = evaluated_coordinates}; |
There was a problem hiding this comment.
Same here, could be moved out of the loop
| for (int gp = 0; gp < integration.num_points(); ++gp) | ||
| { | ||
| const auto xi = Core::Elements::evaluate_parameter_coordinate<celltype>(integration, gp); | ||
| const Core::Elements::ElementNodes<celltype, 2> nodes{.coordinates = evaluated_coordinates}; |
There was a problem hiding this comment.
This could be moved out of the Gauss point loop
| const double factor = shape.values(force_node) * pressure_times_weight * | ||
| reference_thickness * shape.derivatives(0, coordinate_node); | ||
| (*load_linearization)( | ||
| force_node* num_dof_per_node, coordinate_node* num_dof_per_node + 1) -= factor; | ||
| (*load_linearization)( | ||
| force_node* num_dof_per_node + 1, coordinate_node * num_dof_per_node) += factor; |
There was a problem hiding this comment.
Maybe I'm missing something, but doesn't the current computation assume that the applied pressure is spatially homogeneous? Otherwise we would have a contribution from
| { | ||
| using namespace FourC; | ||
|
|
||
| TEST(PlaneNormalPressure, ReferenceThicknessAndConsistentLinearization) |
There was a problem hiding this comment.
As far as I see, the unit test addresses specifically the "helper" function add_normal_pressure_load.
It may be more practical to test the function that really evaluates the Neumann condition (Discret::Elements::evaluate_normal_pressure_by_element), but I understand that it is more tedious because one would have to construct an element, a discretization, and so on for such a test.
I don't have a strong opinion on this and it's fine for me as it is - I just wanted to note that this current approach isolates a function that is ultimately an internal step of the Neumann condition evaluation.
2 participants
Description and Context
Currently
DESIGN LINE NEUMANN CONDITIONSdoes not consider theTYPEparameter, particularlyTYPE=[pseudo_]orthopressure.This issue is fixed here, i.e., orthopressure and pseudo_orthopressure are honored on plane solid lines; a test is added.
The pressure is scaled by thickness and add follower-load linearization.
Disclosure of AI assistance
GitHub Copilot CLI assisted with implementing the 2D solid-line
orthopressurehandling, follower-load linearization, thickness scaling, and associated unit and regression tests.The resulting code was compiled and validated using the solid-element unit suite and the QUAD9 quarter-cylinder regression test.
Reviewed and tested by me (to my best knowledge)