Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[Core] Modify coordinate_transformation_utilities.h to allow transformation back to global coordinates #12327

Merged
merged 10 commits into from
Jul 22, 2024
Merged

[Core] Modify coordinate_transformation_utilities.h to allow transformation back to global coordinates #12327

Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
a389d45
Refactor coordinate_transformation_utilities.h to allow transformatio…
gzjing Apr 27, 2024
6ecadbf
Update kratos/utilities/coordinate_transformation_utilities.h
VeronikaSinger Jul 11, 2024
bec795b
RotateRHSAux private, Indentation, rename TToGlobalCoord
VeronikaSinger Jul 11, 2024
be4d308
Merge branch 'master' into mpm/refactor_rotation_utility
gzjing Jul 13, 2024
73ff95f
Minor edit
gzjing Jul 13, 2024
7fdb38d
Update kratos/utilities/coordinate_transformation_utilities.h
gzjing Jul 20, 2024
2b84a26
Update kratos/utilities/coordinate_transformation_utilities.h
gzjing Jul 20, 2024
8626a5d
fix bug from variable renaming
gzjing Jul 20, 2024
dd03f83
Add documentations
gzjing Jul 20, 2024
48515dd
Update kratos/utilities/coordinate_transformation_utilities.h
gzjing Jul 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
285 changes: 190 additions & 95 deletions kratos/utilities/coordinate_transformation_utilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -461,99 +461,31 @@ class CoordinateTransformationUtils {

}

/// RHS only version of Rotate
/**
* @brief RHS only version of Rotate
* @details Rotates the local RHS vector (rLocalVector) to a tangential-normal coordinate system defined by the
* normals of the associated rGeometry
* @param rLocalVector Local RHS vector associated with rGeometry, to be transformed
* @param rGeometry Geometry associated with rLocalVector
*/
virtual void Rotate(TLocalVectorType& rLocalVector,
GeometryType& rGeometry) const
{
//const unsigned int LocalSize = rLocalVector.size(); // We expect this to work both with elements (4 nodes) and conditions (3 nodes)

//unsigned int Index = 0;

if (rLocalVector.size() > 0)
{
if(mBlockSize != mDomainSize) //Monolithic case
{
for(unsigned int j = 0; j < rGeometry.PointsNumber(); ++j)
{
if( this->IsSlip(rGeometry[j]) )
{
if(mDomainSize == 3)
{
array_1d<double,4> aux,aux1;
BoundedMatrix<double,4,4> rRot;
LocalRotationOperator3D<4>(rRot,rGeometry[j]);

for(unsigned int k=0; k<4; k++)
aux[k] = rLocalVector[j*mBlockSize+k];

noalias(aux1) = prod(rRot,aux);

for(unsigned int k=0; k<4; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
else
{
array_1d<double,3> aux,aux1;
BoundedMatrix<double,3,3> rRot;
LocalRotationOperator2D<3>(rRot,rGeometry[j]);

for(unsigned int k=0; k<3; k++)
{
aux[k] = rLocalVector[j*mBlockSize+k];
}

noalias(aux1) = prod(rRot,aux);

for(unsigned int k=0; k<3; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
}
//Index += mBlockSize;
}

}
else //fractional step case
{
for(unsigned int j = 0; j < rGeometry.PointsNumber(); ++j)
{
if( this->IsSlip(rGeometry[j]) )
{
if(mDomainSize == 3)
{
array_1d<double,3> aux,aux1;
BoundedMatrix<double,3,3> rRot;
LocalRotationOperatorPure(rRot,rGeometry[j]);

for(unsigned int k=0; k<3; k++)
aux[k] = rLocalVector[j*mBlockSize+k];

noalias(aux1) = prod(rRot,aux);

for(unsigned int k=0; k<3; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
else
{
array_1d<double,2> aux,aux1;
BoundedMatrix<double,2,2> rRot;
LocalRotationOperatorPure(rRot,rGeometry[j]);

for(unsigned int k=0; k<2; k++)
aux[k] = rLocalVector[j*mBlockSize+k];

noalias(aux1) = prod(rRot,aux);

for(unsigned int k=0; k<2; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
}
//Index += mBlockSize;
}

}

}
RotateRHSAux(rLocalVector, rGeometry, false);
}

/**
* @brief RHS only version of Rotate which reverts the rotation to the original frame
* @details Rotates the local RHS vector (rLocalVector) in the tangential-normal coordinate system back to the
* original coordinate system
* @param rLocalVector Local RHS vector associated with rGeometry, to be transformed
* @param rGeometry Geometry associated with rLocalVector
*/
virtual void RevertRotate(
TLocalVectorType& rLocalVector,
GeometryType& rGeometry) const
{
RotateRHSAux(rLocalVector, rGeometry, true);
}

/// Apply slip boundary conditions to the rotated local contributions.
Expand Down Expand Up @@ -753,7 +685,21 @@ class CoordinateTransformationUtils {
///@name Protected Operations
///@{

template<unsigned int TDim, unsigned int TBlockSize, unsigned int TSkip = 0>
/**
* @brief Auxilliary function to rotate local system contributions between the original and tangential-normal
* frames of reference, to be used when pressure etc. is present
* @details Transforms the local system based on the normal associated with the rGeometry into tangential-
* normal coordinates if TRevertRotation is false. Otherwise, the rotation matrix used is transposed and the
* inverse transformation is applied -- this transforms the system back to the original coordinate system.
* @tparam TDim Number of spatial dimensions
* @tparam TBlockSize Total number of DoFs associated with each node
* @tparam TSkip Auxiliary value to shift the entries of the rotated array
* @tparam TRevertRotation Boolean indicating if the transformation is reverted (true -> revert transformation)
* @param rLocalMatrix Local LHS matrix contribution associated with rGeometry
* @param rLocalVector Local RHS vector contribution associated with rGeometry
* @param rGeometry Geometry associated with local system
*/
template<unsigned int TDim, unsigned int TBlockSize, unsigned int TSkip = 0, bool TRevertRotation = false>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is no longer needed right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unfortunately, the rotate operation on the LHS in MPMApplication is implemented differently from the core. We thus require the two versions of RotateAux with the ability to revert the rotation applied. Since this is already an auxillary function, making the changes here directly makes more sense in my opinion (as compared to extracting the contents to yet another aux function). I can add some documentation here if you think the changes might be confusing.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense for me. I just wanted to confirm that this was not a leftover.

void RotateAux(TLocalMatrixType& rLocalMatrix,
TLocalVectorType& rLocalVector,
GeometryType& rGeometry) const
Expand All @@ -766,6 +712,8 @@ class CoordinateTransformationUtils {
DenseVector<bool> NeedRotation( NumBlocks, false);

std::vector< BoundedMatrix<double,TBlockSize,TBlockSize> > rRot(NumBlocks);
BoundedMatrix<double,TBlockSize,TBlockSize> tmp;

for(unsigned int j = 0; j < NumBlocks; ++j)
{
if( this->IsSlip(rGeometry[j]) )
Expand All @@ -775,14 +723,20 @@ class CoordinateTransformationUtils {

if constexpr (TDim == 2) LocalRotationOperator2D<TBlockSize,TSkip>(rRot[j],rGeometry[j]);
else LocalRotationOperator3D<TBlockSize,TSkip>(rRot[j],rGeometry[j]);

if constexpr (TRevertRotation)
{
noalias(tmp) = trans(rRot[j]);
rRot[j] = tmp;
}
}

//Index += TBlockSize;
}

if(rotations_needed > 0)
{
BoundedMatrix<double,TBlockSize,TBlockSize> mat_block, tmp;
BoundedMatrix<double,TBlockSize,TBlockSize> mat_block;
array_1d<double,TBlockSize> aux, aux1;

for(unsigned int i=0; i<NumBlocks; i++)
Expand Down Expand Up @@ -832,8 +786,19 @@ class CoordinateTransformationUtils {
}
}

//to be used when there is only velocity (no additional pressure or other var block)
template<unsigned int TDim>
/**
* @brief Auxilliary function to rotate local system contributions between the original and tangential-normal
* frames of reference, to be used when there is only velocity (no additional pressure or other var block)
* @details Transforms the local system based on the normal associated with the rGeometry into tangential-
* normal coordinates if TRevertRotation is false. Otherwise, the rotation matrix used is transposed and the
* inverse transformation is applied -- this transforms the system back to the original coordinate system.
* @tparam TDim Number of spatial dimensions
* @tparam TRevertRotation Boolean indicating if the transformation is reverted (true -> revert transformation)
* @param rLocalMatrix Local LHS matrix contribution associated with rGeometry
* @param rLocalVector Local RHS vector contribution associated with rGeometry
* @param rGeometry Geometry associated with local system
*/
template<unsigned int TDim, bool TRevertRotation = false>
rubenzorrilla marked this conversation as resolved.
Show resolved Hide resolved
void RotateAuxPure(TLocalMatrixType& rLocalMatrix,
TLocalVectorType& rLocalVector,
GeometryType& rGeometry) const
Expand All @@ -846,6 +811,8 @@ class CoordinateTransformationUtils {
DenseVector<bool> NeedRotation( NumBlocks, false);

std::vector< BoundedMatrix<double,TDim,TDim> > rRot(NumBlocks);
BoundedMatrix<double,TDim,TDim> tmp;

for(unsigned int j = 0; j < NumBlocks; ++j)
{
if( this->IsSlip(rGeometry[j]) )
Expand All @@ -854,14 +821,20 @@ class CoordinateTransformationUtils {
rotations_needed++;

LocalRotationOperatorPure(rRot[j],rGeometry[j]);
}

if constexpr (TRevertRotation)
{
noalias(tmp) = trans(rRot[j]);
rRot[j] = tmp;
}
}

//Index += mBlockSize;
}

if(rotations_needed > 0)
{
BoundedMatrix<double,TDim,TDim> mat_block, tmp;
BoundedMatrix<double,TDim,TDim> mat_block;
array_1d<double,TDim> aux, aux1;

for(unsigned int i=0; i<NumBlocks; i++)
Expand Down Expand Up @@ -1052,6 +1025,128 @@ class CoordinateTransformationUtils {
///@name Private Operations
///@{

/**
* @brief Auxilliary function to implement a reversible rotation operation on the local RHS vector
* @details Transforms the local RHS vector based on the normal associated with the rGeometry into tangential-
* normal coordinates if RevertRotation == false. Otherwise, the rotation matrix used is transposed and the
* inverse transformation is applied -- this transforms the RHS vector back to the original coordinate system.
* @param rLocalVector Local RHS vector associated with rGeometry
* @param rGeometry Geometry associated with rLocalVector
* @param RevertRotation Boolean indicating if the transformation is reverted (true -> revert transformation)
*/
void RotateRHSAux(
TLocalVectorType& rLocalVector,
GeometryType& rGeometry,
const bool RevertRotation) const
{
if (rLocalVector.size() > 0)
{
if(mBlockSize != mDomainSize) //Monolithic case
{
for(unsigned int j = 0; j < rGeometry.PointsNumber(); ++j)
{
if( this->IsSlip(rGeometry[j]) )
{
if(mDomainSize == 3)
{
array_1d<double,4> aux,aux1;
BoundedMatrix<double,4,4> rRot;
LocalRotationOperator3D<4>(rRot,rGeometry[j]);

for(unsigned int k=0; k<4; k++)
aux[k] = rLocalVector[j*mBlockSize+k];

if(RevertRotation)
{
noalias(aux1) = prod(trans(rRot),aux);
} else{
noalias(aux1) = prod(rRot,aux);
}

for(unsigned int k=0; k<4; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
else
{
array_1d<double,3> aux,aux1;
BoundedMatrix<double,3,3> rRot;
LocalRotationOperator2D<3>(rRot,rGeometry[j]);

for(unsigned int k=0; k<3; k++)
{
aux[k] = rLocalVector[j*mBlockSize+k];
}

if(RevertRotation)
{
noalias(aux1) = prod(trans(rRot),aux);
} else{
noalias(aux1) = prod(rRot,aux);
}

for(unsigned int k=0; k<3; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
}
//Index += mBlockSize;
}

}
else //fractional step case
{
for(unsigned int j = 0; j < rGeometry.PointsNumber(); ++j)
{
if( this->IsSlip(rGeometry[j]) )
{
if(mDomainSize == 3)
{
array_1d<double,3> aux,aux1;
BoundedMatrix<double,3,3> rRot;
LocalRotationOperatorPure(rRot,rGeometry[j]);

for(unsigned int k=0; k<3; k++)
aux[k] = rLocalVector[j*mBlockSize+k];

if(RevertRotation)
{
noalias(aux1) = prod(trans(rRot),aux);
} else{
noalias(aux1) = prod(rRot,aux);
}

for(unsigned int k=0; k<3; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
else
{
array_1d<double,2> aux,aux1;
BoundedMatrix<double,2,2> rRot;
LocalRotationOperatorPure(rRot,rGeometry[j]);

for(unsigned int k=0; k<2; k++)
aux[k] = rLocalVector[j*mBlockSize+k];

if(RevertRotation)
{
noalias(aux1) = prod(trans(rRot),aux);
} else{
noalias(aux1) = prod(rRot,aux);
}

for(unsigned int k=0; k<2; k++)
rLocalVector[j*mBlockSize+k] = aux1[k];
}
}
//Index += mBlockSize;
}

}

}

}


// /// Compute a rotation matrix to transform values from the cartesian base to one oriented with the node's normal
// /**
// * The normal is read from solution step data NORMAL. Use NormalCalculationUtils::CalculateOnSimplex to
Expand Down
Loading