Skip to content

Commit 5b3d76e

Browse files
dhyan1272dhyan1272
committed
Strain Rate at MP calcualtion
1 parent f881586 commit 5b3d76e

File tree

9 files changed

+234
-52
lines changed

9 files changed

+234
-52
lines changed

‎src/pmpo_MPMesh.cpp‎

Lines changed: 71 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,10 +7,76 @@ namespace polyMPO{
77

88
void printVTP_mesh(MPMesh& mpMesh, int printVTPIndex=-1);
99

10+
void MPMesh::calculateStrain(){
11+
auto MPsPosition = p_MPs->getPositions();
12+
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
13+
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
14+
15+
//Mesh Fields
16+
auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<MeshF_TanLatVertexRotatedOverRadius>();
17+
auto gnomProjVtx = p_mesh->getMeshField<MeshF_VtxGnomProj>();
18+
auto gnomProjElmCenter = p_mesh->getMeshField<MeshF_ElmCenterGnomProj>();
19+
auto elm2VtxConn = p_mesh->getElm2VtxConn();
20+
auto velField = p_mesh->getMeshField<MeshF_Vel>();
21+
bool isRotated = p_mesh->getRotatedFlag();
22+
23+
auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
24+
if(mask){
25+
int numVtx = elm2VtxConn(elm,0);
26+
27+
Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
28+
if(isRotated){
29+
position3d[0] = -MPsPosition(mp, 2);
30+
position3d[2] = MPsPosition(mp, 0);
31+
}
32+
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
33+
double mpProjX, mpProjY;
34+
35+
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
36+
37+
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
38+
39+
double basisByArea[maxVtxsPerElm] = {0.0};
40+
initArray(basisByArea,maxVtxsPerElm,0.0);
41+
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
42+
initArray(gradBasisByArea,maxVtxsPerElm,0.0);
43+
44+
compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);
45+
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, gradBasisByArea);
46+
47+
double v11 = 0.0;
48+
double v12 = 0.0;
49+
double v21 = 0.0;
50+
double v22 = 0.0;
51+
double uTanOverR = 0.0;
52+
double vTanOverR = 0.0;
53+
for (int i = 0; i < numVtx; i++){
54+
int iVertex = elm2VtxConn(elm, i+1)-1;
55+
v11 = v11 + gradBasisByArea[i*2 + 0] * velField(iVertex, 0);
56+
v12 = v12 + gradBasisByArea[i*2 + 1] * velField(iVertex, 0);
57+
v21 = v21 + gradBasisByArea[i*2 + 0] * velField(iVertex, 1);
58+
v22 = v22 + gradBasisByArea[i*2 + 1] * velField(iVertex, 1);
59+
uTanOverR = uTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 0);
60+
vTanOverR = vTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 1);
61+
62+
if(MPsAppID(mp)==0){
63+
//printf("i %d Vtx %d Velocity %.15e %.15e \n", i, elm2VtxConn(elm,i+1)-1, velField(elm2VtxConn(elm,i+1)-1, 0), velField(elm2VtxConn(elm,i+1)-1, 1));
64+
//printf("Grads %.15e %.15e \n", gradBasisByArea[i*2 + 0], gradBasisByArea[i*2 + 1]);
65+
printf("Vs %.15e %.15e %.15e %.15e \n",v11, v12, v21, v22);
66+
printf("uvTan %.15e %.15e \n", uTanOverR, vTanOverR);
67+
}
68+
}
69+
}
70+
};
71+
p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate");
72+
}
73+
1074
void MPMesh::calcBasis(bool use3DArea) {
1175
assert(p_mesh->getGeomType() == geom_spherical_surf);
76+
1277
auto MPsPosition = p_MPs->getPositions();
13-
auto mp_basis_field = p_MPs->getData<MPF_Basis_Vals>(); // we can implement MPs->getBasisVals() like MPs->getPositions()
78+
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
79+
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
1480

1581
auto elm2VtxConn = p_mesh->getElm2VtxConn();
1682
auto vtxCoords = p_mesh->getMeshField<MeshF_VtxCoords>();
@@ -38,6 +104,9 @@ void MPMesh::calcBasis(bool use3DArea) {
38104

39105
double basisByArea[maxVtxsPerElm] = {0.0};
40106
initArray(basisByArea,maxVtxsPerElm,0.0);
107+
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
108+
initArray(gradBasisByArea,maxVtxsPerElm,0.0);
109+
41110

42111
if(!use3DArea){
43112
double mpProjX, mpProjY;
@@ -55,7 +124,7 @@ void MPMesh::calcBasis(bool use3DArea) {
55124
}
56125
// fill step
57126
for(int i=0; i<= numVtx; i++){
58-
mp_basis_field(mp,i) = basisByArea[i];
127+
MPsBasis(mp,i) = basisByArea[i];
59128
}
60129
}
61130
};

‎src/pmpo_MPMesh.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,8 @@ class MPMesh{
112112
void assembleField(int vtxPerElm, int nCells, int nVerticesSolve, int nVertices, double* array_sub, double* array_full);
113113

114114
void printVTP_mesh(int printVTPIndex);
115+
116+
void calculateStrain();
115117
};
116118

117119
}//namespace polyMPO end

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 12 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -335,29 +335,24 @@ void MPMesh::reconstruct_coeff_full(){
335335
if(p_mesh->getGeomType() == geom_spherical_surf)
336336
radius=p_mesh->getSphereRadius();
337337

338-
bool scaling=false;
339-
340338
//Assemble matrix for each vertex
341339
auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
342340
if(mask) { //if material point is 'active'/'enabled'
343341
int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
344342
for(int i=0; i<nVtxE; i++){
345343
int vID = elm2VtxConn(elm,i+1)-1; //vID = vertex id
346344
double w_vtx=weight(mp,i);
347-
double mScale=1;
348-
if(scaling)
349-
mScale=sqrt(dual_triangle_area(vID,0))/radius;
350-
351-
Kokkos::atomic_add(&vtxMatrices(vID,0), w_vtx*mScale*mScale);
352-
Kokkos::atomic_add(&vtxMatrices(vID,1), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))/radius);
353-
Kokkos::atomic_add(&vtxMatrices(vID,2), w_vtx*mScale*(-vtxCoords(vID,1)+mpPos(mp,1))/radius);
354-
Kokkos::atomic_add(&vtxMatrices(vID,3), w_vtx*mScale*(-vtxCoords(vID,2)+mpPos(mp,2))/radius);
355-
Kokkos::atomic_add(&vtxMatrices(vID,4), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,0)+mpPos(mp,0))/(radius*radius));
356-
Kokkos::atomic_add(&vtxMatrices(vID,5), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,1)+mpPos(mp,1))/(radius*radius));
357-
Kokkos::atomic_add(&vtxMatrices(vID,6), w_vtx*mScale*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
358-
Kokkos::atomic_add(&vtxMatrices(vID,7), w_vtx*mScale*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,1)+mpPos(mp,1))/(radius*radius));
359-
Kokkos::atomic_add(&vtxMatrices(vID,8), w_vtx*mScale*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
360-
Kokkos::atomic_add(&vtxMatrices(vID,9), w_vtx*mScale*(-vtxCoords(vID,2)+mpPos(mp,2))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
345+
346+
Kokkos::atomic_add(&vtxMatrices(vID,0), w_vtx);
347+
Kokkos::atomic_add(&vtxMatrices(vID,1), w_vtx*(-vtxCoords(vID,0)+mpPos(mp,0))/radius);
348+
Kokkos::atomic_add(&vtxMatrices(vID,2), w_vtx*(-vtxCoords(vID,1)+mpPos(mp,1))/radius);
349+
Kokkos::atomic_add(&vtxMatrices(vID,3), w_vtx*(-vtxCoords(vID,2)+mpPos(mp,2))/radius);
350+
Kokkos::atomic_add(&vtxMatrices(vID,4), w_vtx*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,0)+mpPos(mp,0))/(radius*radius));
351+
Kokkos::atomic_add(&vtxMatrices(vID,5), w_vtx*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,1)+mpPos(mp,1))/(radius*radius));
352+
Kokkos::atomic_add(&vtxMatrices(vID,6), w_vtx*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
353+
Kokkos::atomic_add(&vtxMatrices(vID,7), w_vtx*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,1)+mpPos(mp,1))/(radius*radius));
354+
Kokkos::atomic_add(&vtxMatrices(vID,8), w_vtx*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
355+
Kokkos::atomic_add(&vtxMatrices(vID,9), w_vtx*(-vtxCoords(vID,2)+mpPos(mp,2))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
361356
}
362357
}
363358
};
@@ -589,7 +584,7 @@ void MPMesh::reconstruct_full() {
589584
//Kokkos::fence();
590585
//assert(cudaDeviceSynchronize() == cudaSuccess);
591586
Kokkos::parallel_for("printSymmetricBlock", numVertices, KOKKOS_LAMBDA(const int vtx){
592-
if (vtx == 2630) {
587+
if (vtx >= 10 && vtx <= 10) {
593588
printf("Field in %d: ", vtx);
594589
for (int k=0; k<numEntries; k++)
595590
printf(" %.15e ", meshField(vtx, k));

‎src/pmpo_c.cpp‎

Lines changed: 19 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -628,30 +628,13 @@ void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs,
628628
}
629629

630630
//TODO: implement these
631-
void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn){
631+
void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh){
632632
checkMPMeshValid(p_mpmesh);
633-
auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
634-
PMT_ALWAYS_ASSERT(nComps == 6); //TODO: mp_sym_mat3d_t
635-
PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
636-
PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
637-
638-
auto mpStrainRate = p_MPs->getData<polyMPO::MPF_Vel>();
639-
auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
640-
kkViewHostU<const double**> mpStrainRateIn_h(mpStrainRateIn,nComps,numMPs);
641-
Kokkos::View<double**> mpStrainRateIn_d("mpStrainRateDevice",nComps,numMPs);
642-
Kokkos::deep_copy(mpStrainRateIn_d, mpStrainRateIn_h);
643-
auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
644-
if(mask){
645-
mpStrainRate(mp,0) = mpStrainRateIn_d(0, mpAppID(mp));
646-
mpStrainRate(mp,1) = mpStrainRateIn_d(1, mpAppID(mp));
647-
mpStrainRate(mp,2) = mpStrainRateIn_d(2, mpAppID(mp));
648-
mpStrainRate(mp,3) = mpStrainRateIn_d(3, mpAppID(mp));
649-
mpStrainRate(mp,4) = mpStrainRateIn_d(4, mpAppID(mp));
650-
mpStrainRate(mp,5) = mpStrainRateIn_d(5, mpAppID(mp));
651-
}
652-
};
653-
p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate");
633+
auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh);
634+
mpMesh->calculateStrain();
654635
}
636+
637+
655638
void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStrainRateHost){
656639
checkMPMeshValid(p_mpmesh);
657640
std::cerr << "Error: This routine is not implemented yet\n";
@@ -1184,6 +1167,20 @@ void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c
11841167
Kokkos::deep_copy(arrayHost, array_d);
11851168
}
11861169

1170+
void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){
1171+
//chech validity
1172+
checkMPMeshValid(p_mpmesh);
1173+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1174+
1175+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1176+
//copy the host array to the device
1177+
auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<polyMPO::MeshF_TanLatVertexRotatedOverRadius>();
1178+
auto h_tanLatVertexRotatedOverRadius = Kokkos::create_mirror_view(tanLatVertexRotatedOverRadius);
1179+
for(int i=0; i<nVertices; i++)
1180+
h_tanLatVertexRotatedOverRadius(i, 0) = array[i];
1181+
Kokkos::deep_copy(tanLatVertexRotatedOverRadius, h_tanLatVertexRotatedOverRadius);
1182+
}
1183+
11871184
bool polympo_push1P_f(MPMesh_ptr p_mpmesh){
11881185
checkMPMeshValid(p_mpmesh);
11891186
bool is_migrating=((polyMPO::MPMesh*)p_mpmesh)->push1P();

‎src/pmpo_c.h‎

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -22,13 +22,9 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm);
2222
void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int numElms, const int numMPs, int* mpsPerElm, const int* mp2Elm, const int* isMPActive);
2323
void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, const int numMPs, const int* allTgtMpElmIn, const int* addedMPMask);
2424
void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, const int size1, const int* arg1, const int size2, const int size3, int* arg2, int* arg3);
25-
2625
int polympo_getMPCount_f(MPMesh_ptr p_mpmesh);
27-
2826
void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh);
29-
3027
void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs);
31-
3228
void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs);
3329
void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs);
3430
void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag);
@@ -45,12 +41,12 @@ void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int n
4541
void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost);
4642
void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn);
4743
void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost);
48-
44+
//MP fields
4945
void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn);
5046
void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpMassHost);
5147
void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn);
5248
void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost);
53-
void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn);
49+
void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh);
5450
void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStrainRateHost);
5551
void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn);
5652
void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost);
@@ -98,7 +94,9 @@ void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* x
9894
//Area Triangle
9995
void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* areaTriangle);
10096
void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle);
101-
97+
//TanLatRotVertexOverR
98+
void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
99+
102100
// calculations
103101
void polympo_push_f(MPMesh_ptr p_mpmesh);
104102
void polympo_push_ahead_f(MPMesh_ptr p_mpmesh);
@@ -112,7 +110,7 @@ void polympo_setReconstructionOfVel_f(MPMesh_ptr p_mpmesh, const int order, cons
112110
void polympo_setReconstructionOfStrainRate_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType);
113111
void polympo_setReconstructionOfStress_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType);
114112
void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh);
115-
113+
//Simpler/cleaner way of calling reconstruction
116114
void polympo_reconstruct_coeff_with_MPI_f(MPMesh_ptr p_mpmesh);
117115
void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh);
118116
void polympo_reconstruct_velocity_with_MPI_f(MPMesh_ptr p_mpmesh);

‎src/pmpo_fortran.f90‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -341,6 +341,13 @@ subroutine polympo_getMPVel(mpMesh, nComps, numMPs, array) &
341341
integer(c_int), value :: nComps, numMPs
342342
type(c_ptr), value :: array
343343
end subroutine
344+
345+
subroutine polympo_setMPStrainRate(mpMesh) &
346+
bind(C, NAME='polympo_setMPStrainRate_f')
347+
use :: iso_c_binding
348+
type(c_ptr), value :: mpMesh
349+
end subroutine
350+
344351
!---------------------------------------------------------------------------
345352
!> @brief Enable the setting of mesh topology (number of entities and entity adjacencies).
346353
!> By default, modifying the mesh topology without calling this function first will result
@@ -785,6 +792,15 @@ subroutine polympo_getMeshVtxOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
785792
integer(c_int), value :: nComps, nVertices
786793
type(c_ptr), value :: array
787794
end subroutine
795+
796+
subroutine polyMPO_setTanLatVertexRotatedOverRadius(mpMesh, nVertices, array) &
797+
bind(C, NAME=' polyMPO_setTanLatVertexRotatedOverRadius_f')
798+
use :: iso_c_binding
799+
type(c_ptr), value :: mpMesh
800+
integer(c_int), value :: nVertices
801+
type(c_ptr), value :: array
802+
end subroutine
803+
788804
!---------------------------------------------------------------------------
789805
!> @brief set the owning process array
790806
!> @param mpmesh(in/out) MPMesh object

‎src/pmpo_mesh.cpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,11 @@ namespace polyMPO{
4545
auto dualTriangleAreaEntry = meshFields2TypeAndString.at(MeshF_DualTriangleArea);
4646
PMT_ALWAYS_ASSERT(dualTriangleAreaEntry.first == MeshFType_VtxBased);
4747
dualTriangleArea_ = MeshFView<MeshF_DualTriangleArea>(dualTriangleAreaEntry.second,numVtxs_);
48+
49+
auto tanLatVertexRotOverRadiusEntry = meshFields2TypeAndString.at(MeshF_TanLatVertexRotatedOverRadius);
50+
PMT_ALWAYS_ASSERT(tanLatVertexRotOverRadiusEntry.first == MeshFType_VtxBased);
51+
tanLatVertexRotatedOverRadius_ = MeshFView<MeshF_TanLatVertexRotatedOverRadius>(tanLatVertexRotOverRadiusEntry.second, numVtxs_);
52+
4853
}
4954

5055
void Mesh::setMeshElmBasedFieldSize(){

‎src/pmpo_mesh.hpp‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,8 @@ enum MeshFieldIndex{
2828
MeshF_OnSurfDispIncr,
2929
MeshF_RotLatLonIncr,
3030
MeshF_VtxGnomProj,
31-
MeshF_ElmCenterGnomProj
31+
MeshF_ElmCenterGnomProj,
32+
MeshF_TanLatVertexRotatedOverRadius
3233
};
3334
enum MeshFieldType{
3435
MeshFType_Invalid = -2,
@@ -50,6 +51,7 @@ template <> struct meshFieldToType < MeshF_OnSurfDispIncr > { using type = Ko
5051
template <> struct meshFieldToType < MeshF_RotLatLonIncr > { using type = Kokkos::View<vec2d_t*>; };
5152
template <> struct meshFieldToType < MeshF_VtxGnomProj > { using type = Kokkos::View<double*[maxVtxsPerElm][2]>; };
5253
template <> struct meshFieldToType < MeshF_ElmCenterGnomProj > { using type = Kokkos::View<double*[4]>; };
54+
template <> struct meshFieldToType < MeshF_TanLatVertexRotatedOverRadius > { using type = Kokkos::View<doubleSclr_t*>; };
5355

5456
template <MeshFieldIndex index>
5557
using MeshFView = typename meshFieldToType<index>::type;
@@ -68,7 +70,8 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType, std::string>> meshFields
6870
{MeshF_OnSurfDispIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}},
6971
{MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}},
7072
{MeshF_VtxGnomProj, {MeshFType_ElmBased,"MeshField_VertexGnomonicProjection"}},
71-
{MeshF_ElmCenterGnomProj,{MeshFType_ElmBased,"MeshField_ElementCenterGnomonicprojection"}}
73+
{MeshF_ElmCenterGnomProj,{MeshFType_ElmBased,"MeshField_ElementCenterGnomonicprojection"}},
74+
{MeshF_TanLatVertexRotatedOverRadius, {MeshFType_VtxBased,"MeshField_TanLatVertexRotatedOverRadius"}},
7275
};
7376

7477
enum mesh_type {mesh_unrecognized_lower = -1,
@@ -112,6 +115,7 @@ class Mesh {
112115
MeshFView<MeshF_VtxGnomProj> vtxGnomProj_;
113116
MeshFView<MeshF_ElmCenterGnomProj> elmCenterGnomProj_;
114117
//DoubleMat2DView vtxStress_;
118+
MeshFView<MeshF_TanLatVertexRotatedOverRadius> tanLatVertexRotatedOverRadius_;
115119
bool isRotatedFlag = false;
116120

117121
public:
@@ -245,6 +249,9 @@ auto Mesh::getMeshField(){
245249
else if constexpr (index==MeshF_ElmCenterGnomProj){
246250
return elmCenterGnomProj_;
247251
}
252+
else if constexpr (index==MeshF_TanLatVertexRotatedOverRadius){
253+
return tanLatVertexRotatedOverRadius_;
254+
}
248255
fprintf(stderr,"Mesh Field Index error!\n");
249256
exit(1);
250257
}
@@ -270,7 +277,6 @@ void computeGnomonicProjectionAtPoint(const Vec3d& Coord,
270277
outY = iDen * (Coord[2] * gnomProjElmCenter_sub(3) - Coord[1] * gnomProjElmCenter_sub(2) * gnomProjElmCenter_sub(0) -
271278
Coord[0] * gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(2));
272279
}
273-
274280
}
275281

276282
#endif

0 commit comments

Comments
 (0)