Skip to content

Eg/add mesh vtx strain rate #62

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

Open
EamonnGlynn wants to merge 8 commits into
base: cws/pumipicDps
Choose a base branch
from
+111 −14
Open

Eg/add mesh vtx strain rate #62

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
5de5e5e
Added StrainRate to c.cpp, c.h, mesh.cpp and mesh.hpp
EamonnGlynn May 20, 2024
a756fe0
added a coma in a list
EamonnGlynn May 20, 2024
cf6783e
fixed naming error
EamonnGlynn May 20, 2024
3b81cd4
sam issue as previous commit but in a different file
EamonnGlynn May 20, 2024
3f4c3b4
Added test cases in fortan
EamonnGlynn May 27, 2024
efb2d77
Added test cases for VtxStrainRate functions
EamonnGlynn Jun 1, 2024
4b2f050
fixed line length in fortran file
EamonnGlynn Jun 1, 2024
3dfd842
fixed indentation, change strain rate functions and tests to take 1 2…
EamonnGlynn Jun 9, 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
42 changes: 42 additions & 0 deletions src/pmpo_c.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -699,6 +699,48 @@ void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons
Kokkos::deep_copy(arrayHost, array_d);
}

void polympo_setMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, const double *xNormal, const double *yNormal, const double *zNormal, const double *xyShear, const double *xzShear, const double *yzShear){
Copy link

@onkarsahni onkarsahni Jun 4, 2024

Choose a reason for hiding this comment

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

Need to discuss if each component (of the 6 components of symmetric tensor) will be given as a separate array.

Copy link
Author

@EamonnGlynn EamonnGlynn Jun 9, 2024

Choose a reason for hiding this comment

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

It has been converted to take one 2d array instead

//check mpMesh is valid
checkMPMeshValid(p_mpmesh);
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;

//check the size
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);

auto strainRate = p_mesh->getMeshField<polyMPO::MeshF_VtxStrainRate>();
auto h_strainRate = Kokkos::create_mirror_view(strainRate);
for(int i = 0; i < nVertices; i++){
h_strainRate(i,0) = xNormal[i];
h_strainRate(i,1) = yNormal[i];
h_strainRate(i,2) = zNormal[i];
h_strainRate(i,3) = xyShear[i];
h_strainRate(i,4) = xzShear[i];
h_strainRate(i,5) = yzShear[i];
}
Copy link

@onkarsahni onkarsahni Jun 4, 2024

Choose a reason for hiding this comment

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

Need to discuss this numbering, i.e., for 6 components of symmetric tensor per vertex (i.e., 0 through 5).

Copy link
Author

@EamonnGlynn EamonnGlynn Jun 9, 2024

Choose a reason for hiding this comment

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

function now uses one 2d array instead of six 1d arrays.

Kokkos::deep_copy(strainRate, h_strainRate);
}

void polympo_getMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, double *xNormal, double *yNormal, double *zNormal, double *xyShear, double *xzShear, double *yzShear){
//check mpMesh is valid
checkMPMeshValid(p_mpmesh);
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;

//check the size
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);

auto strainRate = p_mesh->getMeshField<polyMPO::MeshF_VtxStrainRate>();
auto h_strainRate = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), strainRate);
for(int i = 0; i < nVertices; i++){
xNormal[i] = h_strainRate(i,0);
yNormal[i] = h_strainRate(i,1);
zNormal[i] = h_strainRate(i,2);
xyShear[i] = h_strainRate(i,3);
xzShear[i] = h_strainRate(i,4);
yzShear[i] = h_strainRate(i,5);
}
}


void polympo_push_f(MPMesh_ptr p_mpmesh){
checkMPMeshValid(p_mpmesh);
((polyMPO::MPMesh*)p_mpmesh) ->push();
Expand Down
2 changes: 2 additions & 0 deletions src/pmpo_c.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@ void polympo_setMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons
void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
void polympo_setMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* xNormal, const double* yNormal, const double* zNormal, const double* xyShear, const double* xzShear, const double* yzShear);
void polympo_getMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, double* xNormal, double* yNormal, double* zNormal, double* xyShear, double* xzShear, double* yzShear);

// calculations
void polympo_push_f(MPMesh_ptr p_mpmesh);
Expand Down
26 changes: 26 additions & 0 deletions src/pmpo_fortran.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -537,6 +537,32 @@ subroutine polympo_getMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
type(c_ptr), value :: array
end subroutine
!---------------------------------------------------------------------------
!> @brief set the verte strain rate from a host array
!> @param mpmesh(in/out) MPMesh object
!> @param nVertices(in) numVertices
!> @param x/y/zNormal(in) input mesh normal strain rate 1D array (numVtx)
!> @param xy/xz/yzShear(in) input mesh shear strain rate 1D array (numVtx)
subroutine polympo_setMeshVtxStrainRate(mpMesh, nVertices, xNormal, yNormal, zNormal, xyShear, xzShear, yzShear) &
bind(C, NAME='polympo_setMeshVtxStrainRate_f')
use :: iso_c_binding
type(c_ptr), value :: mpMesh
integer(c_int), value :: nVertices
type(c_ptr), intent(in), value :: xNormal, yNormal, zNormal, xyShear, xzShear, yzShear
end subroutine
!---------------------------------------------------------------------------
!> @brief get the verte strain rate from polyMPO
!> @param mpmesh(in/out) MPMesh object
!> @param nVertices(in) numVertices
!> @param x/y/zNormal(in/out) output mesh normal strain rate 1D array (numVtx), allocated by user
!> @param xy/xz/yzShear(in/out) output mesh shear strain rate 1D array (numVtx), allocated by user
subroutine polympo_getMeshVtxStrainRate(mpMesh, nVertices, xNormal, yNormal, zNormal, xyShear, xzShear, yzShear) &
bind(C, NAME='polympo_getMeshVtxStrainRate_f')
use :: iso_c_binding
type(c_ptr), value :: mpMesh
integer(c_int), value :: nVertices
type(c_ptr), value :: xNormal, yNormal, zNormal, xyShear, xzShear, yzShear
end subroutine
!---------------------------------------------------------------------------
!> @brief calculate the MPs from given mesh vertices rotational latitude
!> longitude, update the MP slices
!> MPs MUST have rotated flag set to True(>0)
Expand Down
4 changes: 4 additions & 0 deletions src/pmpo_mesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,10 @@ namespace polyMPO{
auto vtxRotLatLonIncrMapEntry = meshFields2TypeAndString.at(MeshF_RotLatLonIncr);
PMT_ALWAYS_ASSERT(vtxRotLatLonIncrMapEntry.first == MeshFType_VtxBased);
vtxRotLatLonIncr_ = DoubleVec2dView(vtxRotLatLonIncrMapEntry.second,numVtxs_);

auto vtxStrainRateMapEntry = meshFields2TypeAndString.at(MeshF_VtxStrainRate);
PMT_ALWAYS_ASSERT(vtxStrainRateMapEntry.first == MeshFType_VtxBased);
vtxStrainRate_ = DoubleSymMat3dView(vtxStrainRateMapEntry.second,numVtxs_);
}

void Mesh::computeRotLatLonIncr(){
Expand Down
12 changes: 9 additions & 3 deletions src/pmpo_mesh.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,8 @@ enum MeshFieldIndex{
MeshF_Vel,
MeshF_OnSurfVeloIncr,
MeshF_OnSurfDispIncr,
MeshF_RotLatLonIncr
MeshF_RotLatLonIncr,
MeshF_VtxStrainRate
};
enum MeshFieldType{
MeshFType_Invalid = -2,
Expand All @@ -39,10 +40,11 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType,
{MeshF_Unsupported, {MeshFType_Unsupported,"MeshField_Unsupported"}},
{MeshF_VtxCoords, {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
{MeshF_VtxRotLat, {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
{MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}},
{MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}},
{MeshF_OnSurfVeloIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}},
{MeshF_OnSurfDispIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}},
{MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}}};
{MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}},
{MeshF_VtxStrainRate, {MeshFType_VtxBased,"MeshField_StrainRate"}}};

enum mesh_type {mesh_unrecognized_lower = -1,
mesh_general_polygonal, //other meshes
Expand Down Expand Up @@ -74,6 +76,7 @@ class Mesh {
DoubleVec2dView vtxOnSurfVeloIncr_;
DoubleVec2dView vtxOnSurfDispIncr_;
DoubleVec2dView vtxRotLatLonIncr_;
DoubleSymMat3dView vtxStrainRate_;
//DoubleMat2DView vtxStress_;

public:
Expand Down Expand Up @@ -161,6 +164,9 @@ auto Mesh::getMeshField(){
else if constexpr (index==MeshF_RotLatLonIncr){
return vtxRotLatLonIncr_;
}
else if constexpr (index==MeshF_VtxStrainRate){
return vtxStrainRate_;
}
fprintf(stderr,"Mesh Field Index error!\n");
exit(1);
}
Expand Down
22 changes: 22 additions & 0 deletions test/testFortran.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,28 @@ program main
call assert((xArray(i) .eq. i+value1), "Assert MeshVel u-component Velocity Fail")
call assert((yArray(i) .eq. value2-i), "Assert MeshVel v-component Velocity Fail")
end do

!VtxStrainRate needs 6 arrays, using x,y,zArrays for all 6
Copy link

@onkarsahni onkarsahni Jun 4, 2024

Choose a reason for hiding this comment

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

Only in this file, change 2d array Mesharray name to Mesharray_vel and use a new 2d array Mesharray_strainrate(numCompsStrainRate,nverts) for "VtxStrainRate" testing with 6 values per vertex

Copy link
Author

@EamonnGlynn EamonnGlynn Jun 9, 2024

Choose a reason for hiding this comment

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

Changes have been made.

do i = 1, nverts
xArray(i) = i + value1
yArray(i) = value2 - i
zArray(i) = i + (value1 * value2)
end do
call polympo_setMeshVtxStrainRate(mpMesh, nverts, c_loc(xArray),c_loc(yArray), &
c_loc(zArray),c_loc(xArray),c_loc(yArray),c_loc(zArray))
xArray = -1
yArray = -1
zArray = -1
call polympo_getMeshVtxStrainRate(mpMesh, nverts, c_loc(xArray),c_loc(yArray), &
c_loc(zArray),c_loc(xArray),c_loc(yArray),c_loc(zArray))
do i = 1, nverts
call assert((xArray(i) .eq. i+value1), &
"Assert MeshVtxStrainRate xx-component Velocity Fail")
call assert((yArray(i) .eq. value2-i), &
"Assert MeshVtxStrainRate yy-component Velocity Fail")
call assert((zArray(i) .eq. i+(value1*value2)), &
"Assert MeshVtxStrainRate zz-component Velocity Fail")
end do

deallocate(MParray)
deallocate(Mesharray)
Expand Down