@@ -43,67 +43,6 @@ void sphericalInterpolation(MPMesh& mpMesh){
4343 pumipic::RecordTime (" PolyMPO_sphericalInterpolation" , timer.seconds ());
4444}
4545
46- KOKKOS_INLINE_FUNCTION
47- void compute2DplanarTriangleArea (int numVtx,
48- const Kokkos::View<double [maxVtxsPerElm][2 ], Kokkos::LayoutStride,
49- Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnom_vtx_subview,
50- double mpProjX, double mpProjY, double * basis){
51-
52- double vertCoords[2 ][maxVtxsPerElm + 1 ];
53- for (int i = 0 ; i < numVtx; ++i) {
54- vertCoords[0 ][i] = gnom_vtx_subview (i, 0 );
55- vertCoords[1 ][i] = gnom_vtx_subview (i, 1 );
56- }
57- vertCoords[0 ][numVtx] = vertCoords[0 ][0 ];
58- vertCoords[1 ][numVtx] = vertCoords[1 ][0 ];
59-
60- // Helper lambda for 2D triangle area
61- auto triArea = [&](const double p1[2 ], const double p2[2 ], const double p3[2 ]) -> double {
62- return 0.5 * (p1[0 ] * (p2[1 ] - p3[1 ]) - p2[0 ] * (p1[1 ] - p3[1 ]) + p3[0 ] * (p1[1 ] - p2[1 ]));
63- };
64-
65- // Compute areaV and areaXV
66- double areaV[maxVtxsPerElm];
67- double areaXV[maxVtxsPerElm];
68- double xy[2 ] = {mpProjX, mpProjY};
69-
70- // Special case
71- double p1[2 ] = { vertCoords[0 ][numVtx - 1 ], vertCoords[1 ][numVtx - 1 ] };
72- double p2[2 ] = { vertCoords[0 ][0 ], vertCoords[1 ][0 ] };
73- double p3[2 ] = { vertCoords[0 ][1 ], vertCoords[1 ][1 ] };
74- areaV[0 ] = triArea (p1, p2, p3);
75- double q1[2 ] = { vertCoords[0 ][0 ], vertCoords[1 ][0 ] };
76- double q3[2 ] = { vertCoords[0 ][1 ], vertCoords[1 ][1 ] };
77- areaXV[0 ] = triArea (q1, xy, q3);
78-
79- for (int i = 1 ; i < numVtx; ++i) {
80- double p1[2 ] = { vertCoords[0 ][i - 1 ], vertCoords[1 ][i - 1 ] };
81- double p2[2 ] = { vertCoords[0 ][i], vertCoords[1 ][i] };
82- double p3[2 ] = { vertCoords[0 ][i + 1 ], vertCoords[1 ][i + 1 ] };
83- areaV[i] = triArea (p1, p2, p3);
84- double q1[2 ] = { vertCoords[0 ][i], vertCoords[1 ][i] };
85- double q3[2 ] = { vertCoords[0 ][i + 1 ], vertCoords[1 ][i + 1 ] };
86- areaXV[i] = triArea (q1, xy, q3);
87- }
88-
89- // Wachspress weights
90- double denominator = 0.0 ;
91- for (int i = 0 ; i < numVtx; ++i){
92- double product = areaV[i];
93- for (int j = 0 ; j < numVtx - 2 ; ++j) {
94- int ind1 = (i + j + 1 ) % numVtx;
95- product *= areaXV[ind1];
96- }
97- basis[i] = product;
98- denominator += product;
99- }
100- // Normalize
101- for (int i = 0 ; i < numVtx; ++i){
102- basis[i] /= denominator;
103- // printf("i %d basis %.15e \n", i, basis[i]);
104- }
105- }
106-
10746KOKKOS_INLINE_FUNCTION
10847void wachpress_weights_grads_2D (int numVtx, const Kokkos::View<double [maxVtxsPerElm][2 ],
10948 Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnom_vtx_subview,
0 commit comments