Calculate mesh face areas

src/Ume/CMakeLists.txt

@@ -28,6 +28,7 @@ set(UME_INCLUDE_FILES
   SOA_Idx_Zones.hh
   Timer.hh
   VecN.hh
+  face_area.hh
   gradient.hh
   soa_idx_helpers.hh
   utils.hh
@@ -48,6 +49,7 @@ add_library(Ume
   SOA_Idx_Points.cc
   SOA_Idx_Sides.cc
   SOA_Idx_Zones.cc
+  face_area.cc
   gradient.cc
   utils.cc
   )

src/Ume/VecN.hh

@@ -165,5 +165,10 @@ inline void normalize(Vec3 &a) {
   a /= std::sqrt(mag);
 }
 
+constexpr double vectormag(Vec3 const &a) {
+  double const mag = std::sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
+  return mag;
+}
+
 } // namespace Ume
 #endif

src/Ume/face_area.cc

@@ -0,0 +1,56 @@
+/*
+  Copyright (c) 2023, Triad National Security, LLC. All rights reserved.
+
+  This is open source software; you can redistribute it and/or modify it under
+  the terms of the BSD-3 License. If software is modified to produce derivative
+  works, such modified software should be clearly marked, so as not to confuse
+  it with the version available from LANL. Full text of the BSD-3 License can be
+  found in the LICENSE.md file, and the full assertion of copyright in the
+  NOTICE.md file.
+*/
+
+/*!
+  \file Ume/face_area.cc
+*/
+#include "Ume/face_area.hh"
+
+namespace Ume {
+
+using Mesh = SOA_Idx::Mesh;
+using DBLV_T = DS_Types::DBLV_T;
+using INTV_T = DS_Types::INTV_T;
+using VEC3V_T = DS_Types::VEC3V_T;
+
+void calc_face_area(Mesh &mesh, DBLV_T &face_area) {
+  auto const &side_type = mesh.sides.mask;
+  auto const &face_comm_type = mesh.faces.comm_type;
+  auto const &s_to_f_map = mesh.ds->caccess_intv("m:s>f");
+  auto const &s_to_s2_map = mesh.ds->caccess_intv("m:s>s2");
+  auto const &surz = mesh.ds->caccess_vec3v("side_surz");
+
+  int const sll = mesh.sides.size();
+  int const sl = mesh.sides.local_size();
+
+  std::fill(face_area.begin(), face_area.end(), 0.0);
+  INTV_T side_tag(sll, 0);
+
+  for (int s = 0; s < sl; ++s) {
+    if (side_type[s] < 1)
+      continue; // We want internal sides only
+    if (side_tag[s] == 1)
+      continue; // Already added this side via s2
+
+    int const f = s_to_f_map[s];
+    if (face_comm_type[f] < 3) { // Internal or master face
+      double const side_area = vectormag(surz[s]); // Flat area
+      face_area[f] += side_area;
+
+      int const s2 = s_to_s2_map[s];
+      side_tag[s2] = 1;
+    }
+  }
+
+  mesh.faces.scatter(face_area);
+}
+
+} // namespace Ume

src/Ume/face_area.hh

@@ -0,0 +1,31 @@
+/*
+  Copyright (c) 2023, Triad National Security, LLC. All rights reserved.
+
+  This is open source software; you can redistribute it and/or modify it under
+  the terms of the BSD-3 License. If software is modified to produce derivative
+  works, such modified software should be clearly marked, so as not to confuse
+  it with the version available from LANL. Full text of the BSD-3 License can be
+  found in the LICENSE.md file, and the full assertion of copyright in the
+  NOTICE.md file.
+*/
+
+/*!
+   \file Ume/face_area.hh
+*/
+#ifndef UME_FACE_AREA_HH
+#define UME_FACE_AREA_HH 1
+
+#include "Ume/SOA_Idx_Mesh.hh"
+
+namespace Ume {
+
+//! Calculate the flat face area
+/* This method computes the face areas for each face in the mesh by
+ * accumulating side surface area vector magnitudes of internal sides
+ * corresponding to internal or master faces in the MPI comm stencil,
+ * making sure not to double count sides. */
+void calc_face_area(SOA_Idx::Mesh &mesh, DS_Types::DBLV_T &face_area);
+
+} // namespace Ume
+
+#endif

src/ume_mpi.cc

@@ -25,6 +25,7 @@
 #include "Ume/Comm_MPI.hh"
 #include "Ume/SOA_Idx_Mesh.hh"
 #include "Ume/Timer.hh"
+#include "Ume/face_area.hh"
 #include "Ume/gradient.hh"
 #include "Ume/utils.hh"
 #include <cassert>
@@ -107,7 +108,7 @@ int main(int argc, char *argv[]) {
   if (comm.pe() == 0) {
     std::cout << "Original algorithm took: " << orig_time.seconds() << "s\n";
     std::cout << "Inverted algorithm took: " << invert_time.seconds() << "s\n";
-    std::cout << "Checking result..." << std::endl;
+    std::cout << "Checking gradient result..." << std::endl;
   }
 
   if (zgrad != zgrad_invert) {
@@ -159,6 +160,20 @@ int main(int argc, char *argv[]) {
               << grad_points.size() << " expected " << z2p.size(czi) << '\n';
   }
 
+  // Do a face area calculation
+  if (comm.pe() == 0)
+    std::cout << "Calculating face areas..." << std::endl;
+
+  // Create a result vector and initialize to impossible value
+  DBLV_T face_area(mesh.faces.size(), -100000.0);
+
+  orig_time.start();
+  Ume::calc_face_area(mesh, face_area);
+  orig_time.stop();
+
+  if (comm.pe() == 0)
+    std::cout << "Face area calculation took: " << orig_time.seconds() << "s\n";
+
   if (comm.pe() == 0)
     std::cout << "Done." << std::endl;
   comm.stop();