[MPM] Compute Cauchy stress vector on grid nodes from mp#14437
[MPM] Compute Cauchy stress vector on grid nodes from mp#14437lauramoremar wants to merge 23 commits into
Conversation
…ng but not working in VTK)
| def _InitializeNodalCauchyStressVector(self): | ||
| if self._GetDomainSize() == 3: | ||
| stress_size = 6 | ||
| elif self.settings["axis_symmetric_flag"].GetBool(): | ||
| stress_size = 4 | ||
| else: | ||
| stress_size = 3 | ||
|
|
||
| zero_stress_vector = KratosMultiphysics.Vector(stress_size) | ||
| for i in range(stress_size): | ||
| zero_stress_vector[i] = 0.0 | ||
|
|
||
| for node in self.grid_model_part.Nodes: | ||
| node.SetSolutionStepValue(KratosMultiphysics.CAUCHY_STRESS_VECTOR, 0, zero_stress_vector) | ||
|
|
There was a problem hiding this comment.
Is it possible to replace this code with a call to the cpp function MPMNodalCauchyStressUtility::ResetNodalCauchyStress?
There was a problem hiding this comment.
I considered reusing MPMNodalCauchyStressUtility::ResetNodalCauchyStress, but I kept the current implementation since ResetNodalCauchyStress is currently a private method and is only internally called by CalculateNodalCauchyStress, which performs additional operations beyond the reset itself. In other words, the utility would need a refactoring to reuse it directly from Python.
Do you think it would be a good idea to make that change?
| for (std::size_t i = 0; i < number_of_nodes; ++i) { | ||
|
|
||
| Vector nodal_cauchy_stress_vector = cauchy_stress_values[0]; | ||
| const double stress_weight = r_N(0,i) * mass_values[0]; |
There was a problem hiding this comment.
Even though you are calling this utility in the scheme FinalizeSolutionStep (so after material points positions have been updated), the values of the the shape functions refer to the position of the material points at the beginning of the time step, no? Is it the intended behavior?
There was a problem hiding this comment.
Yes, you are right. The shape functions correspond to the latest P2G/search, which is performed at the beginning of the step. The current implementation intentionally projects the updated MP stress using that same grid association, consistently with the available NODAL_MASS. A projection using the final MP positions would require an additional search/update of the MP quadrature geometry together with a recomputation of the nodal weights.
The utility is called in FinalizeSolutionStep because we want to project the updated MP stress of the current time step, rather than the stress from the previous one. However, this does not imply that the shape functions correspond to the final MP positions, but to the latest search/P2G mapping.
I don't know if it is preferable to compute at the beginning of the time step, but in that case the MP stress is the corresponding to the previous one, if I am not wrong.
| void FinalizeSolutionStep( | ||
| ModelPart& rModelPart, | ||
| TSystemMatrixType& rA, | ||
| TSystemVectorType& rDx, | ||
| TSystemVectorType& rb) override | ||
| { | ||
| BossakBaseType::FinalizeSolutionStep(rModelPart, rA, rDx, rb); | ||
| const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo(); | ||
| const bool compute_nodal_cauchy_stress = (r_current_process_info.Has(COMPUTE_NODAL_CAUCHY_STRESS)) | ||
| ? r_current_process_info.GetValue(COMPUTE_NODAL_CAUCHY_STRESS) | ||
| : false; | ||
|
|
||
| if(mFrictionIsActive) { | ||
| block_for_each(mGridModelPart.Nodes(), [&](Node& rNode) | ||
| { | ||
| const Node& rConstNode = rNode; // const Node reference to avoid issues with previously unset GetValue() | ||
| if( mRotationTool.IsConformingSlip(rConstNode) && rConstNode.GetValue(FRICTION_COEFFICIENT) > 0 ) | ||
| rNode.FastGetSolutionStepValue(REACTION).clear(); | ||
| }); | ||
|
|
||
| mRotationTool.ComputeFrictionAndResetFlags(rModelPart); | ||
| } | ||
|
|
||
| block_for_each(mGridModelPart.Nodes(), [&](Node& rNode) { | ||
| const Node& rConstNode = rNode; // const Node reference to avoid issues with previously unset GetValue() | ||
|
|
||
| // rotate forces stored in REACTION to global coordinates on conforming boundaries | ||
| if (mRotationTool.IsConformingSlip(rConstNode) ) { | ||
| mRotationTool.RotateVector(rNode.FastGetSolutionStepValue(REACTION), rConstNode, true); | ||
| } | ||
| }); | ||
|
|
||
| if (compute_nodal_cauchy_stress) { | ||
| MPMNodalCauchyStressUtility::CalculateNodalCauchyStress(rModelPart, mGridModelPart); | ||
| } | ||
|
|
||
| } |
There was a problem hiding this comment.
I understand that the computation performed by this utility is post processing but P2G (particle to grid) mapping is usually performed at the beginning of the time step and not at the end
What do you think @andimkatili?
There was a problem hiding this comment.
I think this comment is related to the previous one (regarding why the projection is computed at the end of the time step).
I am open to changing the position of this calculation, keeping in mind that in that case the postprocessed CAUCHY_STRESS_VECTOR would correspond to the previous time step instead of the current one.
3 participants
📝 Description
This PR adds optional projection of material point Cauchy stresses to the background grid nodes in the MPMApplication for post-processing purposes.
When
compute_nodal_cauchy_stressis enabled, the material pointMP_CAUCHY_STRESS_VECTORis projected to the grid nodes and stored as the historical core variableCAUCHY_STRESS_VECTOR. This allows writing the projected stress on the background grid with the standard VTK output process.Key changes
compute_nodal_cauchy_stresssolver setting and correspondingCOMPUTE_NODAL_CAUCHY_STRESSprocess-info flag.MPMNodalCauchyStressUtilityto projectMP_CAUCHY_STRESS_VECTORfrom material points to grid nodes.CAUCHY_STRESS_VECTORnodal variable for the projected stress instead of introducing a new nodal stress variable.CAUCHY_STRESS_VECTORis written with a valid vector size.Validation
Two new tests have been added:
1 - applications/MPMApplication/tests/test_nodal_cauchy_stress.py
This test checks the irreducible and UP formulations using constant material point Cauchy stress vector values.
2 - A C++ test has also been added:
test_mpm_nodal_cauchy_stress_utility.cpp
In this test, one triangular element with two material points is considered. Each material point has different MP_CAUCHY_STRESS_VECTOR values, which are interpolated to the nodes.