Second shape derivatives transform implementation #11253
Conversation
|
Is this intended to be overridden for different geometries? If not then I suggest to move it to the GeomUtils, to not further overload the Geometry |
|
I do not think so, it is analogous to the calculation of the first derivatives on integration points |
There was a problem hiding this comment.
I agree @jgonzalezusua observation. For the user it would be a bit counterintuitive to have the 1st derivatives in the Geometry and 2nd ones in a separate class. I think we should throw an error in those geometries that would require an especific treatment (i.e. non-square Jacobian ones).
In any case, I'd like @RiccardoRossi to have a special look to this one.
kratos/geometries/geometry.h
Outdated
| KRATOS_ERROR_IF_NOT(this->WorkingSpaceDimension() == this->LocalSpaceDimension()) | ||
| << "\'ShapeFunctionsIntegrationPointsSecondDerivatives\' is not defined for current geometry type as second derivatives are only defined in the local space." << std::endl; |
There was a problem hiding this comment.
OK, better than the check I was proposing.
|
@KratosMultiphysics/technical-committee would like to ask why the GlobalSpaceDerivatives won't work in this context. Another concern is the performance as we think that doing the loop over the Gauss points in the element would be more efficient rather than creating a vector of vectors and returning the results. |
As far as I know in |
|
@KratosMultiphysics/technical-committee when do you think we can merge this PR? Do you think I need to change the implementation? |
|
@pooyan-dadvand, did the committe reject the current implementation? |
|
we are pending @sunethwarna to speak with @tteschemacher or @oberbichler about how that funciton should be used |
|
Thank you @jgonzalezusua for your implementation. Generally, that looks like a good feature with good code quality. I do have a couple of comments. Please consider that I am not @KratosMultiphysics/technical-committee , therefore, they will have the last word. I think some of my points are at least important to be discussed. Sorry for my very late reply. I will try to stay in the loop to finish that up quickly.
The GlobalSpaceDerivatives would return the evaluated base vectors and their derivatives in global space. Basically, evaluated shape functions * node coordinates. Your function computes in local space, please add a comment to make it clear? It would be handy, if you could add your function to the GlobalSpaceDerivatives it. However, it is complementary to your function. If you consider doing this, please do this in a later PR. I do have a general question, is your function valid for any type of geometry, or is it restrained to linear mesh elements? Do you need any given polynomial degree (which?). If so, you could add a check to it. |
|
@tteschemacher first of all thank you for taking the time to review this. One of our issues was that by looking at the signature virtual void GlobalSpaceDerivatives(
std::vector<CoordinatesArrayType>& rGlobalSpaceDerivatives,
const CoordinatesArrayType& rLocalCoordinates,
const SizeType DerivativeOrder) constwe do not reallt understand how the output is organized. Can u explain what is the content of rGlobalSpaceDerivatives on output? imagine for example that we have a quadrilateral and we ask for the third derivatives. How would they be organized in the output vector? (i am aware they would be zero...but it is really not clear to us how they would be stored since that output vector has variable size (it is a std:.vector), but always contains entries of size 3 ... |
|
The function itself contains the comment about the order: Kratos/kratos/geometries/geometry.h Lines 2455 to 2472 in 88a8a45 It is: Unfortunately, there's nothing to ensure that this is done right. If that is still difficult to understand, I should update the comment. |
|
Hi @tteschemacher, this function is not restraint to linear elements. In fact, if you use triangle or tetrahedra, the output of this function should be 0. I have implemented a test to compute the second derivatives of a two order polynomial using a quadratic hexahedra and it works well. |
|
@tteschemacher the point is that i am not clear with what you mean by that imagine i have the shape functions N1, N2, N3, N4 referring to the nodes 1,2,3,4 respectively where do i find d^3 N1 / dxdydz ?? where is d^3 N3 / dxdydz ?? |
That vector does not contain shape functions per se. It contains its evaluation towards location and base vectors. Kratos/kratos/geometries/geometry.h Lines 2504 to 2507 in 88a8a45 All that should be in x, y, z coordinates in the global system (That is in type of CoordinatesArrayType). That is why I think the two functions are complemantary and it would be nice to complete them. |
|
sorry @tteschemacher but i really do not understand what you mean by telling "That vector does not contain shape functions per se. It contains its evaluation towards location and base vectors" could you be a bit more explicative about that... |
|
@oberbichler, I hope you could share your opinion about this implementation |
|
Sorry if I am a little tiresome but I have openned this PR last month, and it seems stucked. @KratosMultiphysics/technical-committee, do you have more comments about it? What do we need to do to approve this PR? |
Dear @jgonzalezusua, as I already told @sunethwarna, I am working on new projects and have no time to deal with Kratos. Moreover, I have never worked with Kratos - apart from some short experiments - and have never used Kratos for my research. Therefore, I can neither comment on your code nor explain any Kratos functions. Reading the comments, however, there seems to be a big misunderstanding. @tteschemacher probably has NURBS in mind, while @RiccardoRossi thinks of something else. Points of NURBS are linear combinations of the control points and therefore Best regards |
Ok thanks. If you are not involved in Kratos, I do not know why they request your presence here. Anyway, I hope we can conclude soon this PR. |
kratos/geometries/geometry.h
Outdated
| virtual void ShapeFunctionsIntegrationPointsSecondDerivatives( | ||
| DenseVector<DenseVector<Matrix>>& rResult, | ||
| IntegrationMethod ThisMethod ) const | ||
| { |
There was a problem hiding this comment.
In most of those functions we have actually intended to separate betweem the functions which use the integration method and the functions which are called with all the information obtained from it.
That would allow to use the same function with own integration rules. In most of my cases that had been apparent. Therefore, I would like to ask you if that makes sense here. I know that this deducts slightly the efficiency. But as I checked through the code, it appears that this effect can be kept minor.
If the function is only capable to express the shapes of integration points from the pre-evaluated points, then you should add a comment.
There was a problem hiding this comment.
To compute the second derivatives you need to know the integration method in order to get the number of integration points. This function has been implemented preserving the same format as ShapeFunctionsIntegrationPointsGradients, which does the same for the case of first derivatives.
kratos/geometries/geometry.h
Outdated
| * https://scicomp.stackexchange.com/questions/25196/implementing-higher-order-derivatives-for-finite-element | ||
| */ | ||
| virtual void ShapeFunctionsIntegrationPointsSecondDerivatives( | ||
| DenseVector<DenseVector<Matrix>>& rResult, |
There was a problem hiding this comment.
What is the difference to this data type:
https://github.com/KratosMultiphysics/Kratos/blob/88a8a45caac20686fbd0ddb8f66bf161344c1c2d/kratos/geometries/geometry.h#L195C5-L195C5
There was a problem hiding this comment.
This data type is DenseVector<DenseVector< Matrix > > and the one you cite is DenseVector< Matrix >
kratos/geometries/geometry.h
Outdated
| * The method used can be found here. | ||
| * https://scicomp.stackexchange.com/questions/25196/implementing-higher-order-derivatives-for-finite-element | ||
| */ | ||
| virtual void ShapeFunctionsIntegrationPointsSecondDerivatives( |
There was a problem hiding this comment.
probably I would rename this function to ShapeFunctionsSecondDerivatives to overload the existing function. Or let me know what are the intended differences?
There was a problem hiding this comment.
The function ShapeFunctionsSecondDerivatives is not the same as ShapeFunctionsIntegrationPointsSecondDerivatives.
ShapeFunctionsSecondDerivativescomputes the isoparametric second derivatives of the shape functions and depends on the geometry you use (linear, bilinear, quadratic, biquadratic, etc). Then, it is called from the element you use.ShapeFunctionsIntegrationPointsSecondDerivativescomputes the mapping between the isoparametric and natural coordinates.
So it cannot be renamed as you said.
There was a problem hiding this comment.
Ok, it really looks like this is a problematic one.
I now took at least 1h to try to understand this PR and honestly i am lost.
To assess the current situation, the geometry already has a function to compute second (and third) order derivatives IN THE LOCAL SPACE (xi eta zeta).
your goal is to have that in the global space (x y z)? It that is so, the name or at least the documentation of the function should reflect that very clearly.
could we have a meeting, either online or in person to discuss this?
kratos/geometries/geometry.h
Outdated
| * The method used can be found here. | ||
| * https://scicomp.stackexchange.com/questions/25196/implementing-higher-order-derivatives-for-finite-element | ||
| */ | ||
| virtual void ShapeFunctionsIntegrationPointsSecondDerivatives( |
There was a problem hiding this comment.
Ok, it really looks like this is a problematic one.
I now took at least 1h to try to understand this PR and honestly i am lost.
To assess the current situation, the geometry already has a function to compute second (and third) order derivatives IN THE LOCAL SPACE (xi eta zeta).
your goal is to have that in the global space (x y z)? It that is so, the name or at least the documentation of the function should reflect that very clearly.
could we have a meeting, either online or in person to discuss this?
Yes, we can meet in person or online. If you are in the university we can have the meeting right now. |
|
I am in my office, just pass by if you can |
|
As @RiccardoRossi and @philbucher suggested, I have changes the transform of the second derivatives of the shape functions to |
|
Something happens with the ContactMechanics but I have not modified nothing from it. Others PRs have the same problem... |
|
@RiccardoRossi I did the changes as you suggested |
|
@KratosMultiphysics/technical-committee agrees with the change as it is now an utility and not changing the geometry base class. @RiccardoRossi will review it |
| @@ -1029,6 +1034,86 @@ bool GeometryTesterUtility::VerifyShapeFunctionsSecondDerivativesValues( | |||
| return true; | |||
| } | |||
|
|
|||
| bool GeometryTesterUtility::VerifyShapeFunctionsSecondDerivativesInterpolation( | |||
| GeometryType& rGeometry, | |||
| * @return true If teh test fails, true otherwise¡ | ||
| */ | ||
| bool VerifyShapeFunctionsSecondDerivativesInterpolation( | ||
| GeometryType& rGeometry, |
| /***********************************************************************************/ | ||
|
|
||
| void GeometryUtils::ShapeFunctionsSecondDerivativesTransformOnIntegrationPoint( | ||
| Matrix DN_DX, |
There was a problem hiding this comment.
If so you can pass a const reference to avoid the copy.
| rGeometry.ShapeFunctionsSecondDerivatives(DDN_DDe, rLocalIntegrationPointCoordinates); | ||
|
|
||
| Matrix A, Ainv; | ||
| double DetA; |
There was a problem hiding this comment.
this function will fail if WorkingSpaceDimension is different from LocalSpaceDimension.
i think you should thow an error if that is the case
There was a problem hiding this comment.
This is already done in line 527
| for (IndexType p = 0; p < rGeometry.PointsNumber(); ++p) { | ||
| Vector rhs, result; | ||
| if (rGeometry.WorkingSpaceDimension() == 2){ | ||
| rhs.resize(3); |
There was a problem hiding this comment.
you are resizing this once per every point number. put the resize outside of the loop?
also use the false in the resize
| rhs[2] = DDN_DDe[p](0,1) - DN_DX(p,0) * H[0](0,1) - DN_DX(p,1) * H[1](0,1); | ||
| } | ||
| else if (rGeometry.WorkingSpaceDimension() == 3){ | ||
| rhs.resize(6); |
|
|
||
| aux[p].resize(rGeometry.WorkingSpaceDimension(), rGeometry.WorkingSpaceDimension(), false ); | ||
|
|
||
| result = prod(Ainv, rhs); |
| @@ -345,6 +345,10 @@ bool GeometryTesterUtility::StreamTestQuadrilateral2D9N( | |||
|
|
|||
| array_1d<double,3> point_in(3,1.0/3.0); | |||
| if( !VerifyShapeFunctionsSecondDerivativesValues(geometry,point_in,rErrorMessage) ) successful = false; | |||
| Quadrilateral2D9<NodeType> exact_geometry( rModelPart.pGetNode(1), rModelPart.pGetNode(3), rModelPart.pGetNode(9), rModelPart.pGetNode(7), | |||
There was a problem hiding this comment.
as a general comment, not for this PR, why the geometry_tester is in the utilities? if it is only for testing it should be in the testing folder.
| aux[p](1,2) = result[4]; | ||
| } | ||
| } | ||
| rResult = aux; |
|
@RiccardoRossi, I have made the changes that you suggested. |
| @@ -518,7 +518,7 @@ void GeometryUtils::ShapeFunctionsSecondDerivativesTransformOnAllIntegrationPoin | |||
| /***********************************************************************************/ | |||
|
|
|||
| void GeometryUtils::ShapeFunctionsSecondDerivativesTransformOnIntegrationPoint( | |||
| Matrix DN_DX, | |||
| const Matrix DN_DX, | |||
There was a problem hiding this comment.
you are passing this by copy. is it intentional?
There was a problem hiding this comment.
No, I think we can pass it as a reference
| @@ -519,7 +519,7 @@ class KRATOS_API(KRATOS_CORE) GeometryTesterUtility | |||
| GeometryType& rGeometry, | |||
| GeometryType::CoordinatesArrayType& rGlobalCoordinates, | |||
There was a problem hiding this comment.
rGlobalCoordinates should also be const no?
There was a problem hiding this comment.
In this case I follow the same definition as VerifyShapeFunctionValues method defined just above. Why should we pass rGlobalCoordinates as a const in VerifyShapeFunctionsSecondDerivativesValues but not in VerifyShapeFunctionValues?
There was a problem hiding this comment.
well, mostly because if you do not change it then it is nice to know it is const ... if it is "incorrectly" not const in the other function i would not replicate the same mistake here ...
| @@ -691,7 +691,7 @@ class KRATOS_API(KRATOS_CORE) GeometryUtils | |||
| * @param rResult the transform of the second derivative of the shape function on the integration point | |||
| */ | |||
| static void ShapeFunctionsSecondDerivativesTransformOnIntegrationPoint( | |||
| Matrix DN_DX, | |||
There was a problem hiding this comment.
@RiccardoRossi, I have already changed it
| @@ -519,7 +519,7 @@ class KRATOS_API(KRATOS_CORE) GeometryTesterUtility | |||
| GeometryType& rGeometry, | |||
| GeometryType::CoordinatesArrayType& rGlobalCoordinates, | |||
There was a problem hiding this comment.
well, mostly because if you do not change it then it is nice to know it is const ... if it is "incorrectly" not const in the other function i would not replicate the same mistake here ...
jgonzalezusua
deleted the
second_shape_derivatives_transform_implementation
branch
📝 Description
This PR is aim to implement the transformation of the second derivatives of the shape functions from isoparametric coordinates. The implementation follows this method: https://scicomp.stackexchange.com/questions/25196/implementing-higher-order-derivatives-for-finite-element
Please mark the PR with appropriate tags:
🆕 Changelog
ShapeFunctionsIntegrationPointsSecondDerivativesingeometry.hVerifyShapeFunctionsSecondDerivativesInterpolationingeometry_tester.cpp