From eeb7a99cb24a16884c30ebc65361a35f20e2ee11 Mon Sep 17 00:00:00 2001
From: Joost van Zwieten <joostvanzwieten@evalf.com>
Date: Thu, 14 May 2020 08:42:02 +0200
Subject: [PATCH] WIP: fix flipping normal of manifold refs

---
 nutils/element.py        | 12 ++++++++++++
 nutils/topology.py       |  5 +++--
 tests/test_finitecell.py | 14 +++++++++++++-
 3 files changed, 28 insertions(+), 3 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 5372f73f7e..10f9aa3a36 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -313,6 +313,10 @@ def withmanifoldedges(self):
     else:
       return self
 
+  @property
+  def flipped(self):
+    return self
+
 strictreference = types.strict[Reference]
 
 class EmptyLike(Reference):
@@ -933,6 +937,10 @@ def get_poly_coeffs(self, basis, **kwargs):
   def get_edge_dofs(self, degree, iedge):
     return self.baseref.get_edge_dofs(degree, iedge)
 
+  @property
+  def flipped(self):
+    return OwnChildReference(self.baseref.flipped)
+
 class WithChildrenReference(Reference):
   'base reference with explicit children'
 
@@ -1362,6 +1370,10 @@ def inside(self, point, eps=0):
   def slice(self, levelfunc, ndivisions):
     return ManifoldReference(self.ref.slice(lambda vertices: levelfunc(self.trans.apply(vertices)), ndivisions), self.trans)
 
+  @property
+  def flipped(self):
+    return ManifoldReference(self.ref, self.trans.flipped)
+
 ## UTILITY FUNCTIONS
 
 def parse_legacy_ischeme(ischeme):
diff --git a/nutils/topology.py b/nutils/topology.py
index 004436e4ea..9bcf3ab281 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -820,7 +820,8 @@ class OppositeTopology(Topology):
 
   def __init__(self, basetopo):
     self.basetopo = basetopo
-    super().__init__(basetopo.roots, basetopo.references, basetopo.opposites, basetopo.transforms)
+    refs = elementseq.asreferences((ref.flipped for ref in basetopo.references), self.basetopo.ndims)
+    super().__init__(basetopo.roots, refs, basetopo.opposites, basetopo.transforms)
 
   def getitem(self, item):
     return ~(self.basetopo.getitem(item))
@@ -2592,7 +2593,7 @@ def addnewedge(ielem, etrans):
               iface = aname, bname
             else:
               iface = bname, aname
-              aetrans, betrans = betrans, aetrans
+              aetrans, betrans, aeref, beref = betrans, aetrans, beref, aeref
             newreferences[iface].append(aeref)
             newtransforms[iface].append(addnewedge(ibase, aetrans.separate(todims)))
             newopposites[iface].append(addnewedge(ibase, betrans.separate(todims)))
diff --git a/tests/test_finitecell.py b/tests/test_finitecell.py
index 30c49dc08b..940ff2c1a5 100644
--- a/tests/test_finitecell.py
+++ b/tests/test_finitecell.py
@@ -371,7 +371,8 @@ class partialtrim(TestCase):
   # +-----+-----+
 
   def setUp(self):
-    self.topo, geom = mesh.rectilinear([2,2])
+    self.topo, self.geom = mesh.rectilinear([2,2])
+    geom = self.geom
     if self.method == 'trim':
       self.topoA = self.topo.trim(geom[0]-1+geom[1]*(geom[1]-.5), maxrefine=1)
       self.topoB = self.topo - self.topoA
@@ -415,5 +416,16 @@ def test_baseboundaries(self):
       alttopo = topology.ConnectedTopology(topo.roots, topo.references, topo.transforms, topo.opposites, topo.connectivity)
       self.assertEqual(dict(zip(alttopo.boundary.transforms, alttopo.boundary.references)), dict(zip(topo.boundary.transforms, topo.boundary.references)))
 
+  def test_volumes(self):
+    geom = self.geom
+    f = ((0.5 - geom)**2).sum(axis=0)
+    lhs = self.topoA.integrate(f.grad(geom)*function.J(geom), ischeme='gauss2')
+    rhs = self.topoA.boundary.integrate(f*function.normal(geom)*function.J(geom), ischeme='gauss2')
+    numpy.testing.assert_array_almost_equal(lhs, rhs)
+    lhs = self.topoB.integrate(f.grad(geom)*function.J(geom), ischeme='gauss2')
+    rhs = self.topoB.boundary.integrate(f*function.normal(geom)*function.J(geom), ischeme='gauss2')
+    numpy.testing.assert_array_almost_equal(lhs, rhs)
+
+
 partialtrim(method='trim')
 partialtrim(method='partition')