Skip to content

Commit 533c8ab

Browse files
dhyan1272dhyan1272
committed
CleanUp_v3
1 parent f58ddb6 commit 533c8ab

File tree

4 files changed

+242
-261
lines changed

4 files changed

+242
-261
lines changed

‎src/pmpo_MPMesh.cpp‎

Lines changed: 57 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -31,18 +31,15 @@ void MPMesh::calculateStrain(){
3131
}
3232
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
3333
double mpProjX, mpProjY;
34-
3534
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
36-
3735
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
3836

3937
double basisByArea[maxVtxsPerElm] = {0.0};
4038
initArray(basisByArea,maxVtxsPerElm,0.0);
4139
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
4240
initArray(gradBasisByArea,maxVtxsPerElm,0.0);
4341

44-
compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);
45-
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, gradBasisByArea);
42+
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea, gradBasisByArea);
4643

4744
double v11 = 0.0;
4845
double v12 = 0.0;
@@ -71,64 +68,64 @@ void MPMesh::calculateStrain(){
7168
p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate");
7269
}
7370

74-
void MPMesh::calcBasis(bool use3DArea) {
75-
assert(p_mesh->getGeomType() == geom_spherical_surf);
71+
void MPMesh::calcBasis() {
72+
assert(p_mesh->getGeomType() == geom_spherical_surf);
7673

77-
auto MPsPosition = p_MPs->getPositions();
78-
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
79-
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
74+
auto MPsPosition = p_MPs->getPositions();
75+
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
76+
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
8077

81-
auto elm2VtxConn = p_mesh->getElm2VtxConn();
82-
auto vtxCoords = p_mesh->getMeshField<MeshF_VtxCoords>();
83-
double radius = p_mesh->getSphereRadius();
84-
//For Gnomonic Projection
85-
auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
86-
auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
87-
88-
bool isRotated = p_mesh->getRotatedFlag();
89-
90-
auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
91-
if(mask) { //if material point is 'active'/'enabled'
92-
Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
93-
// formating
94-
Vec3d v3d[maxVtxsPerElm+1];
95-
int numVtx = elm2VtxConn(elm,0);
96-
for(int i = 1; i<=numVtx; i++){
97-
v3d[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0);
98-
v3d[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1);
99-
v3d[i-1][2] = vtxCoords(elm2VtxConn(elm,i)-1,2);
100-
}
101-
v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0);
102-
v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1);
103-
v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2);
104-
105-
double basisByArea[maxVtxsPerElm] = {0.0};
106-
initArray(basisByArea,maxVtxsPerElm,0.0);
107-
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
108-
initArray(gradBasisByArea,maxVtxsPerElm,0.0);
78+
auto elm2VtxConn = p_mesh->getElm2VtxConn();
79+
auto vtxCoords = p_mesh->getMeshField<MeshF_VtxCoords>();
80+
double radius = p_mesh->getSphereRadius();
81+
//For Gnomonic Projection
82+
auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
83+
auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
84+
85+
bool isRotated = p_mesh->getRotatedFlag();
86+
87+
auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
88+
if(mask) { //if material point is 'active'/'enabled'
89+
int numVtx = elm2VtxConn(elm,0);
90+
Vec3d position3d(MPsPosition(mp, 0),MPsPosition(mp, 1),MPsPosition(mp, 2));
91+
if(isRotated){
92+
position3d[0] = -MPsPosition(mp, 2);
93+
position3d[2] = MPsPosition(mp, 0);
94+
}
95+
96+
double mpProjX, mpProjY;
97+
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
98+
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
99+
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
100+
101+
double basisByArea[maxVtxsPerElm] = {0.0};
102+
initArray(basisByArea,maxVtxsPerElm, 0.0);
103+
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
104+
initArray(gradBasisByArea,maxVtxsPerElm, 0.0);
109105

110-
111-
if(!use3DArea){
112-
double mpProjX, mpProjY;
113-
auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
114-
if(isRotated){
115-
position3d[0] = -MPsPosition(mp, 2);
116-
position3d[2] = MPsPosition(mp, 0);
117-
}
118-
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
119-
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
120-
compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);
121-
}
122-
else{
123-
getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea);
124-
}
125-
// fill step
126-
for(int i=0; i<= numVtx; i++){
127-
MPsBasis(mp,i) = basisByArea[i];
128-
}
129-
}
130-
};
131-
p_MPs->parallel_for(calcbasis, "calcbasis");
106+
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea, gradBasisByArea);
107+
108+
for(int i=0; i<= numVtx; i++){
109+
MPsBasis(mp,i) = basisByArea[i];
110+
}
111+
112+
//Old method where basis functions calculated using 3D Area
113+
/*
114+
Vec3d v3d[maxVtxsPerElm+1];
115+
int numVtx = elm2VtxConn(elm,0);
116+
for(int i = 1; i<=numVtx; i++){
117+
v3d[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0);
118+
v3d[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1);
119+
v3d[i-1][2] = vtxCoords(elm2VtxConn(elm,i)-1,2);
120+
}
121+
v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0);
122+
v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1);
123+
v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2);
124+
getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea);
125+
*/
126+
}
127+
};
128+
p_MPs->parallel_for(calcbasis, "calcbasis");
132129
}
133130

134131
void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
@@ -415,11 +412,9 @@ void MPMesh::push_ahead(){
415412
p_mesh->computeRotLatLonIncr();
416413

417414
//Interpolates latitude longitude, mesh velocity increments to MPs
418-
bool use3DArea=false;
419-
calcBasis(use3DArea);
415+
calcBasis();
420416
sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
421417
sphericalInterpolation<MeshF_OnSurfVeloIncr>(*this);
422-
//sphericalInterpolationDispVelIncr(*this);
423418

424419
//Push the MPs
425420
p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius(), p_mesh->getRotatedFlag());

‎src/pmpo_MPMesh.hpp‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -18,9 +18,6 @@ template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_OnSurfVeloI
1818
#define maxMPsPerElm 8
1919

2020
class MPMesh{
21-
private:
22-
23-
bool isPreComputed;
2421

2522
public:
2623

@@ -62,7 +59,7 @@ class MPMesh{
6259

6360
//Used before advection to interpolate fields from mesh to MPs
6461
//And also before reconstruction
65-
void calcBasis(bool use3DArea);
62+
void calcBasis();
6663

6764
//Reconstruction
6865
DoubleView assemblyV0();

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -110,8 +110,7 @@ void MPMesh::reconstruct_coeff_full(){
110110
auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
111111

112112
//Material Points
113-
bool use3DArea=false;
114-
calcBasis(use3DArea);
113+
calcBasis();
115114

116115
auto weight = p_MPs->getData<MPF_Basis_Vals>();
117116
auto mpPos = p_MPs->getData<MPF_Cur_Pos_XYZ>();

0 commit comments

Comments
 (0)