diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index 01eccadce6fcedb1cfaea71d87e959781190cc25..6a525a6aa731baff2b901fe6086a8309cc1ad3e4 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -105,7 +105,25 @@ class MeshData : public IMeshData if constexpr (Dimension == 1) { static_assert(Dimension != 1, "Nlr does not make sense in 1d"); } else if constexpr (Dimension == 2) { - throw NotImplementedError("Nlr are not yet computed in 2d"); + const NodeValue<const Rd>& xr = m_mesh.xr(); + + NodeValuePerFace<Rd> Nlr(m_mesh.connectivity()); + const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix(); + + parallel_for( + m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId l) { + const auto& face_nodes = face_to_node_matrix[l]; + + const Rd xr0 = xr[face_nodes[0]]; + const Rd xr1 = xr[face_nodes[1]]; + const Rd dx = xr1 - xr0; + + const Rd Nr = 0.5 * Rd{dx[1], -dx[0]}; + + Nlr(l, 0) = Nr; + Nlr(l, 1) = Nr; + }); + m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr); } else { const NodeValue<const Rd>& xr = m_mesh.xr();