Skip to content

Commit 89adb70

Browse files
dhyan1272dhyan1272
committed
Starting gnomonic projection in polyMPO
1 parent 9e7a093 commit 89adb70

File tree

4 files changed

+99
-19
lines changed

4 files changed

+99
-19
lines changed

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -684,6 +684,12 @@ void MPMesh::startCommunication(){
684684

685685
pumipic::RecordTime("Start Communication" + std::to_string(self), timer.seconds());
686686

687+
//Temporary
688+
double radius = p_mesh->getSphereRadius();
689+
bool isRotated=false;
690+
//p_mesh->setGnomonicProjection(isRotated, radius);
691+
692+
687693
if (p_MPs->getOpMode() != polyMPO::MP_DEBUG)
688694
return;
689695
printf("Rank %d Owners %d Halos %d Total %d \n", self, numOwnersTot, numHalosTot, numEntities);

‎src/pmpo_c.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,7 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a
7575
void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array);
7676
void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
7777
void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array);
78-
78+
void polympo_setGnomonicProjection(MPMesh_ptr p_mpmesh);
7979

8080
//Mesh fields
8181
int polympo_getMeshFVtxType_f();

‎src/pmpo_mesh.cpp‎

Lines changed: 59 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,6 @@ namespace polyMPO{
4444

4545
auto dualTriangleAreaEntry = meshFields2TypeAndString.at(MeshF_DualTriangleArea);
4646
PMT_ALWAYS_ASSERT(dualTriangleAreaEntry.first == MeshFType_VtxBased);
47-
dualTriangleArea_ = MeshFView<MeshF_DualTriangleArea>(dualTriangleAreaEntry.second,numVtxs_);
4847
}
4948

5049
void Mesh::setMeshElmBasedFieldSize(){
@@ -53,7 +52,12 @@ namespace polyMPO{
5352
auto elmMassMapEntry = meshFields2TypeAndString.at(MeshF_ElmMass);
5453
PMT_ALWAYS_ASSERT(elmMassMapEntry.first == MeshFType_ElmBased);
5554
elmMass_ = MeshFView<MeshF_ElmMass>(elmMassMapEntry.second,numElms_);
56-
elmCenterXYZ_ = MeshFView<MeshF_ElmCenterXYZ>(meshFields2TypeAndString.at(MeshF_ElmCenterXYZ).second,numElms_);
55+
56+
elmCenterXYZ_ = MeshFView<MeshF_ElmCenterXYZ>(meshFields2TypeAndString.at(MeshF_ElmCenterXYZ).second, numElms_);
57+
58+
elmCenterGnomProj_= MeshFView<MeshF_ElmCenterGnomProj>(meshFields2TypeAndString.at(MeshF_ElmCenterGnomProj).second, numElms_);
59+
60+
vtxGnomProj_ = MeshFView<MeshF_VtxGnomProj>(meshFields2TypeAndString.at(MeshF_VtxGnomProj).second, numElms_);
5761
}
5862

5963
void Mesh::computeRotLatLonIncr(){
@@ -74,4 +78,57 @@ namespace polyMPO{
7478
pumipic::RecordTime("PolyMPO_computeRotLatLonIncr", timer.seconds());
7579
}
7680

81+
KOKKOS_INLINE_FUNCTION
82+
void computeGnomonicProjectionAtPoint( const Vec3d& vtxCoord, const double gnomProjElmCenter[4], double& outX, double& outY) {
83+
const double iDen = 1.0 / ( gnomProjElmCenter[1] * gnomProjElmCenter[3] * vtxCoord[0] +
84+
gnomProjElmCenter[0] * gnomProjElmCenter[3] * vtxCoord[1] +
85+
gnomProjElmCenter[2] * vtxCoord[2]);
86+
outX = iDen * (vtxCoord[1] * gnomProjElmCenter[1] -
87+
vtxCoord[1] * gnomProjElmCenter[0]);
88+
outY = iDen * (vtxCoord[2] * gnomProjElmCenter[3] - vtxCoord[1] * gnomProjElmCenter[2] * gnomProjElmCenter[0] -
89+
vtxCoord[0] * gnomProjElmCenter[1] * gnomProjElmCenter[2]);
90+
}
91+
92+
void Mesh::setGnomonicProjection(bool isRotated, double radius){
93+
std::cout<<__FUNCTION__<<std::endl;
94+
auto gnomProjVtx = getMeshField<MeshF_VtxGnomProj>();
95+
auto gnomProjElmCenter = getMeshField<MeshF_ElmCenterGnomProj>();
96+
97+
auto vtxCoords = getMeshField<MeshF_VtxCoords>();
98+
auto elmCenters = getMeshField<MeshF_ElmCenterXYZ>();
99+
auto elm2VtxConn = getElm2VtxConn();
100+
101+
102+
Kokkos::parallel_for("setGnomprojCenter", numElms_, KOKKOS_LAMBDA(const int iElm){
103+
104+
Vec3d elmCenter(elmCenters(iElm, 0), elmCenters(iElm, 1), elmCenters(iElm, 2));
105+
if(isRotated){
106+
107+
}
108+
auto cos2LatR = elmCenter[0]*elmCenter[0] + elmCenter[1]*elmCenter[1];
109+
auto invR = 1.0/ sqrt(cos2LatR + elmCenter[2]*elmCenter[2]);
110+
auto cosLatR = sqrt(cos2LatR);
111+
112+
gnomProjElmCenter(iElm, 0) = elmCenter[1]/cosLatR;
113+
gnomProjElmCenter(iElm, 1) = elmCenter[0]/cosLatR;
114+
gnomProjElmCenter(iElm, 2) = invR*elmCenter[2];
115+
gnomProjElmCenter(iElm, 3) = invR*cosLatR;
116+
117+
int nVtxE = elm2VtxConn(iElm,0);
118+
for(int i=0; i<nVtxE; i++){
119+
int vID = elm2VtxConn(iElm, i+1) - 1;
120+
Vec3d vtxCord(vtxCoords(vID, 0), vtxCoords(vID, 1), vtxCoords(vID, 2));
121+
if (isRotated){
122+
123+
}
124+
125+
double outX, outY;
126+
computeGnomonicProjectionAtPoint(vtxCord, &gnomProjElmCenter(iElm, 0), outX, outY);
127+
128+
gnomProjVtx(iElm, i, 0) = outX;
129+
gnomProjVtx(iElm, i, 1) = outY;
130+
}
131+
});
132+
}
133+
77134
} // namespace polyMPO

‎src/pmpo_mesh.hpp‎

Lines changed: 33 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,9 @@ enum MeshFieldIndex{
2626
MeshF_ElmMass,
2727
MeshF_OnSurfVeloIncr,
2828
MeshF_OnSurfDispIncr,
29-
MeshF_RotLatLonIncr
29+
MeshF_RotLatLonIncr,
30+
MeshF_VtxGnomProj,
31+
MeshF_ElmCenterGnomProj
3032
};
3133
enum MeshFieldType{
3234
MeshFType_Invalid = -2,
@@ -46,24 +48,28 @@ template <> struct meshFieldToType < MeshF_ElmMass > { using type = Ko
4648
template <> struct meshFieldToType < MeshF_OnSurfVeloIncr > { using type = Kokkos::View<vec2d_t*>; };
4749
template <> struct meshFieldToType < MeshF_OnSurfDispIncr > { using type = Kokkos::View<vec2d_t*>; };
4850
template <> struct meshFieldToType < MeshF_RotLatLonIncr > { using type = Kokkos::View<vec2d_t*>; };
51+
template <> struct meshFieldToType < MeshF_VtxGnomProj > { using type = Kokkos::View<double*[maxVtxsPerElm][2]>; };
52+
template <> struct meshFieldToType < MeshF_ElmCenterGnomProj > { using type = Kokkos::View<double*[4]>; };
4953

5054
template <MeshFieldIndex index>
5155
using MeshFView = typename meshFieldToType<index>::type;
5256

53-
const std::map<MeshFieldIndex, std::pair<MeshFieldType,
54-
std::string>> meshFields2TypeAndString =
55-
{{MeshF_Invalid, {MeshFType_Invalid,"MeshField_InValid!"}},
56-
{MeshF_Unsupported, {MeshFType_Unsupported,"MeshField_Unsupported"}},
57-
{MeshF_VtxCoords, {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
58-
{MeshF_VtxRotLat, {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
59-
{MeshF_ElmCenterXYZ, {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}},
60-
{MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}},
61-
{MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}},
62-
{MeshF_VtxMass, {MeshFType_VtxBased,"MeshField_VerticesMass"}},
63-
{MeshF_ElmMass, {MeshFType_ElmBased,"MeshField_ElementsMass"}},
64-
{MeshF_OnSurfVeloIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}},
65-
{MeshF_OnSurfDispIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}},
66-
{MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}}};
57+
const std::map<MeshFieldIndex, std::pair<MeshFieldType, std::string>> meshFields2TypeAndString = {
58+
{MeshF_Invalid, {MeshFType_Invalid,"MeshField_InValid!"}},
59+
{MeshF_Unsupported, {MeshFType_Unsupported,"MeshField_Unsupported"}},
60+
{MeshF_VtxCoords, {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
61+
{MeshF_VtxRotLat, {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
62+
{MeshF_ElmCenterXYZ, {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}},
63+
{MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}},
64+
{MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}},
65+
{MeshF_VtxMass, {MeshFType_VtxBased,"MeshField_VerticesMass"}},
66+
{MeshF_ElmMass, {MeshFType_ElmBased,"MeshField_ElementsMass"}},
67+
{MeshF_OnSurfVeloIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}},
68+
{MeshF_OnSurfDispIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}},
69+
{MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}},
70+
{MeshF_VtxGnomProj, {MeshFType_ElmBased,"MeshField_VertexGnomonicProjection"}},
71+
{MeshF_ElmCenterGnomProj,{MeshFType_ElmBased,"MeshField_ElementCenterGnomonicprojection"}}
72+
};
6773

6874
enum mesh_type {mesh_unrecognized_lower = -1,
6975
mesh_general_polygonal, //other meshes
@@ -102,6 +108,9 @@ class Mesh {
102108
MeshFView<MeshF_OnSurfVeloIncr> vtxOnSurfVeloIncr_;
103109
MeshFView<MeshF_OnSurfDispIncr> vtxOnSurfDispIncr_;
104110
MeshFView<MeshF_RotLatLonIncr> vtxRotLatLonIncr_;
111+
//GnomonicProjection
112+
MeshFView<MeshF_VtxGnomProj> vtxGnomProj_;
113+
MeshFView<MeshF_ElmCenterGnomProj> elmCenterGnomProj_;
105114
//DoubleMat2DView vtxStress_;
106115

107116
public:
@@ -126,7 +135,7 @@ class Mesh {
126135
setMeshElmBasedFieldSize();
127136
meshEdit_ = false;
128137
vtxCoords_ = vtxCoords;
129-
}
138+
}
130139

131140
bool meshEditable(){ return meshEdit_; }
132141
bool checkMeshType(int meshType);
@@ -177,6 +186,8 @@ class Mesh {
177186
IntView getElmGlobal() {return globalElm_;}
178187
IntView getVtxGlobal() {return globalVtx_;}
179188

189+
void setGnomonicProjection(bool isRotated, double radius);
190+
180191
void computeRotLatLonIncr();
181192
};
182193

@@ -220,6 +231,12 @@ auto Mesh::getMeshField(){
220231
else if constexpr (index==MeshF_RotLatLonIncr){
221232
return vtxRotLatLonIncr_;
222233
}
234+
else if constexpr (index==MeshF_VtxGnomProj){
235+
return vtxGnomProj_;
236+
}
237+
else if constexpr (index==MeshF_ElmCenterGnomProj){
238+
return elmCenterGnomProj_;
239+
}
223240
fprintf(stderr,"Mesh Field Index error!\n");
224241
exit(1);
225242
}

0 commit comments

Comments
 (0)