From c0b02fba776eb7b4b1c3ef91d57db47b966fdf3b Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Thu, 20 Oct 2016 11:39:06 +0200
Subject: [PATCH] fixed: open edge constraints in ASMu2D

the assumption of the corners being the first / last node
in the list was wrong.
---
 src/ASM/LR/ASMu2D.C          | 36 +++++++++++++++------
 src/ASM/LR/Test/TestASMu2D.C | 61 ++++++++++++++++++++++++++++++++++++
 2 files changed, 87 insertions(+), 10 deletions(-)

diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C
index 71418157..cd34c882 100644
--- a/src/ASM/LR/ASMu2D.C
+++ b/src/ASM/LR/ASMu2D.C
@@ -616,14 +616,33 @@ std::vector<int> ASMu2D::getEdgeNodes (int edge, int basis) const
 void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
 {
   LR::parameterEdge edge;
+  int c1, c2;
   switch (dir) {
-  case -2: edge = LR::SOUTH; break;
-  case -1: edge = LR::WEST; break;
-  case  1: edge = LR::EAST; break;
-  case  2: edge = LR::NORTH; break;
+  case -2:
+    edge = LR::SOUTH;
+    c1 = this->getCorner(-1, -1, basis);
+    c2 = this->getCorner( 1, -1, basis);
+    break;
+  case -1:
+    edge = LR::WEST;
+    c1 = this->getCorner(-1, -1, basis);
+    c2 = this->getCorner(-1,  1, basis);
+    break;
+  case  1:
+    edge = LR::EAST;
+    c1 = this->getCorner( 1, -1, basis);
+    c2 = this->getCorner( 1,  1, basis);
+    break;
+  case  2:
+    edge = LR::NORTH;
+    c1 = this->getCorner(-1,  1, basis);
+    c2 = this->getCorner( 1,  1, basis);
+    break;
   default: return;
   }
   std::vector<int> nodes = this->getEdgeNodes(edge,basis);
+  nodes.erase(std::find(nodes.begin(), nodes.end(), c1));
+  nodes.erase(std::find(nodes.begin(), nodes.end(), c2));
 
   int bcode = code;
   if (code > 0) {
@@ -634,18 +653,15 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
     bcode = -code;
 
   // Skip the first and last function if we are requesting an open boundary.
-  // I here assume the edgeFunctions are ordered such that the physical
-  // end points are represented by the first and last edgeFunction.
   if (!open)
-    this->prescribe(nodes.front(),dof,bcode);
+    this->prescribe(c1,dof,bcode);
 
-
-  for (size_t i = 1; i < nodes.size()-1; i++)
+  for (size_t i = 0; i < nodes.size(); i++)
     if (this->prescribe(nodes[i],dof,-code) == 0 && code > 0)
       dirich.back().nodes.push_back(std::make_pair(i,nodes[i]));
 
   if (!open)
-    this->prescribe(nodes.back(),dof,bcode);
+    this->prescribe(c2,dof,bcode);
 
   if (code > 0)
     if (dirich.back().nodes.empty())
diff --git a/src/ASM/LR/Test/TestASMu2D.C b/src/ASM/LR/Test/TestASMu2D.C
index 1ef262ec..3b862e5b 100644
--- a/src/ASM/LR/Test/TestASMu2D.C
+++ b/src/ASM/LR/Test/TestASMu2D.C
@@ -16,6 +16,21 @@
 
 #include "gtest/gtest.h"
 
+
+struct EdgeTest {
+  int edge;
+  int edgeIdx;
+  std::array<int,2> c1;
+  std::array<int,2> c2;
+};
+
+
+class TestASMu2D : public testing::Test,
+                   public testing::WithParamInterface<EdgeTest>
+{
+};
+
+
 TEST(TestASMu2D, BoundaryNodesE1)
 {
   SIM2D sim(1);
@@ -78,3 +93,49 @@ TEST(TestASMu2D, BoundaryNodesE4)
   for (int i = 0; i < 4; ++i)
     ASSERT_EQ(vec[i], 5+i);
 }
+
+
+TEST_P(TestASMu2D, ConstrainEdge)
+{
+  SIM2D sim(1);
+  sim.opt.discretization = ASM::LRSpline;
+  ASSERT_TRUE(sim.read("src/ASM/LR/Test/refdata/boundary_nodes.xinp"));
+  ASMu2D* pch = static_cast<ASMu2D*>(sim.getPatch(1));
+  pch->constrainEdge(GetParam().edgeIdx, false, 1, 1, 1);
+  std::vector<int> glbNodes;
+  pch->getBoundaryNodes(GetParam().edge, glbNodes, 1);
+  for (auto& it : glbNodes)
+    ASSERT_TRUE(pch->findMPC(it, 1) != nullptr);
+}
+
+
+TEST_P(TestASMu2D, ConstrainEdgeOpen)
+{
+  SIM2D sim(1);
+  sim.opt.discretization = ASM::LRSpline;
+  ASSERT_TRUE(sim.read("src/ASM/LR/Test/refdata/boundary_nodes.xinp"));
+  ASMu2D* pch = static_cast<ASMu2D*>(sim.getPatch(1));
+  pch->constrainEdge(GetParam().edgeIdx, true, 1, 1, 1);
+  std::vector<int> glbNodes;
+  pch->getBoundaryNodes(GetParam().edge, glbNodes, 1);
+  int crn = pch->getCorner(GetParam().c1[0], GetParam().c1[1], 1);
+  ASSERT_TRUE(pch->findMPC(crn, 1) == nullptr);
+  glbNodes.erase(std::find(glbNodes.begin(), glbNodes.end(), crn));
+  crn = pch->getCorner(GetParam().c2[0], GetParam().c2[1], 1);
+  ASSERT_TRUE(pch->findMPC(crn, 1) == nullptr);
+  glbNodes.erase(std::find(glbNodes.begin(), glbNodes.end(), crn));
+  for (auto& it : glbNodes)
+    ASSERT_TRUE(pch->findMPC(it, 1) != nullptr);
+}
+
+
+static const std::vector<EdgeTest> edgeTestData =
+        {{1, -1, {-1, -1}, {-1 , 1}},
+         {2,  1, { 1, -1}, { 1 , 1}},
+         {3, -2, {-1, -1}, { 1, -1}},
+         {4,  2, {-1,  1}, { 1,  1}}};
+
+
+
+INSTANTIATE_TEST_CASE_P(TestASMu2D, TestASMu2D,
+                        testing::ValuesIn(edgeTestData));