@@ -11,19 +11,19 @@ void MPMesh::calculateStrain(){
1111 auto MPsPosition = p_MPs->getPositions ();
1212 auto MPsBasis = p_MPs->getData <MPF_Basis_Vals>();
1313 auto MPsAppID = p_MPs->getData <MPF_MP_APP_ID>();
14-
14+
1515 // Mesh Fields
1616 auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField <MeshF_TanLatVertexRotatedOverRadius>();
1717 auto gnomProjVtx = p_mesh->getMeshField <MeshF_VtxGnomProj>();
1818 auto gnomProjElmCenter = p_mesh->getMeshField <MeshF_ElmCenterGnomProj>();
1919 auto elm2VtxConn = p_mesh->getElm2VtxConn ();
2020 auto velField = p_mesh->getMeshField <MeshF_Vel>();
2121 bool isRotated = p_mesh->getRotatedFlag ();
22-
22+
2323 auto setMPStrainRate = PS_LAMBDA (const int & elm, const int & mp, const int & mask){
2424 if (mask){
2525 int numVtx = elm2VtxConn (elm,0 );
26-
26+
2727 Vec3d position3d (MPsPosition (mp,0 ),MPsPosition (mp,1 ),MPsPosition (mp,2 ));
2828 if (isRotated){
2929 position3d[0 ] = -MPsPosition (mp, 2 );
@@ -33,14 +33,14 @@ void MPMesh::calculateStrain(){
3333 double mpProjX, mpProjY;
3434 computeGnomonicProjectionAtPoint (position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
3535 auto gnom_vtx_subview = Kokkos::subview (gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
36-
36+
3737 double basisByArea[maxVtxsPerElm] = {0.0 };
3838 initArray (basisByArea,maxVtxsPerElm,0.0 );
3939 double gradBasisByArea[2 *maxVtxsPerElm] = {0.0 };
4040 initArray (gradBasisByArea,maxVtxsPerElm,0.0 );
41-
41+
4242 wachpress_weights_grads_2D (numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea, gradBasisByArea);
43-
43+
4444 double v11 = 0.0 ;
4545 double v12 = 0.0 ;
4646 double v21 = 0.0 ;
@@ -55,13 +55,6 @@ void MPMesh::calculateStrain(){
5555 v22 = v22 + gradBasisByArea[i*2 + 1 ] * velField (iVertex, 1 );
5656 uTanOverR = uTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius (iVertex, 0 ) * velField (iVertex, 0 );
5757 vTanOverR = vTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius (iVertex, 0 ) * velField (iVertex, 1 );
58-
59- if (MPsAppID (mp)==0 ){
60- // 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));
61- // printf("Grads %.15e %.15e \n", gradBasisByArea[i*2 + 0], gradBasisByArea[i*2 + 1]);
62- printf (" Vs %.15e %.15e %.15e %.15e \n " ,v11, v12, v21, v22);
63- printf (" uvTan %.15e %.15e \n " , uTanOverR, vTanOverR);
64- }
6558 }
6659 }
6760 };
@@ -70,11 +63,11 @@ void MPMesh::calculateStrain(){
7063
7164void MPMesh::calcBasis () {
7265 assert (p_mesh->getGeomType () == geom_spherical_surf);
73-
66+
7467 auto MPsPosition = p_MPs->getPositions ();
7568 auto MPsBasis = p_MPs->getData <MPF_Basis_Vals>();
7669 auto MPsAppID = p_MPs->getData <MPF_MP_APP_ID>();
77-
70+
7871 auto elm2VtxConn = p_mesh->getElm2VtxConn ();
7972 auto vtxCoords = p_mesh->getMeshField <MeshF_VtxCoords>();
8073 double radius = p_mesh->getSphereRadius ();
@@ -134,10 +127,10 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
134127 auto elm2VtxConn = p_mesh->getElm2VtxConn ();
135128 auto elm2ElmConn = p_mesh->getElm2ElmConn ();
136129 auto MPs2Elm = p_MPs->getData <MPF_Tgt_Elm_ID>();
137- const auto vtxCoords = p_mesh->getMeshField <polyMPO::MeshF_VtxCoords>();
130+ const auto vtxCoords = p_mesh->getMeshField <polyMPO::MeshF_VtxCoords>();
138131 auto mpPositions = p_MPs->getData <MPF_Cur_Pos_XYZ>();
139132 Kokkos::View<Vec2d*[maxVtxsPerElm]> edgeCenters (" EdgeCenters" ,numElms);
140-
133+
141134 Kokkos::parallel_for (" calcEdgeCenter" , numElms, KOKKOS_LAMBDA (const int elm){
142135 int numVtx = elm2VtxConn (elm,0 );
143136 int v[maxVtxsPerElm];
@@ -150,7 +143,7 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
150143 edgeCenters (elm,i) = (v_ip1 + v_i)*0.5 ;
151144 }
152145 });
153-
146+
154147 auto CVTEdgeTracking = PS_LAMBDA (const int & elm, const int & mp, const int & mask){
155148 Vec2d MP (mpPositions (mp,0 ),mpPositions (mp,1 ));// XXX:the input is XYZ, but we only support 2d vector
156149 if (mask){
@@ -173,7 +166,7 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
173166 }
174167 if (currentDistSq < minDistSq){
175168 edgeIndex = i+1 ;
176- minDistSq = currentDistSq;
169+ minDistSq = currentDistSq;
177170 }
178171 }
179172 if (edgeIndex <0 ){
@@ -232,7 +225,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
232225
233226 Vec3d center (elmCenter (iElm, 0 ), elmCenter (iElm, 1 ), elmCenter (iElm, 2 ));
234227 Vec3d delta = MPnew - center;
235-
228+
236229 double minDistSq = delta[0 ]*delta[0 ] + delta[1 ]*delta[1 ] + delta[2 ]*delta[2 ];
237230 int closestElm = -1 ;
238231 // go through all the connected elm, calc distance
0 commit comments