Differing results for same trimmed structure on different background meshes

[MichaelLoibl]

Dear all,

I am using the trimming functionality of Nutils. When playing around with different background meshes, I obtained different results even though I expected from the underlying methodology that they should be the same. I add a simple example below which computes the area of a triangle in different background meshes. The nodes of the triangle are: [0,0.9;0.9,0;0.9,0.9]. Am I doing a mistake in the definition of the background mesh or the evaluation of the integral?

from nutils import mesh, cli
from nutils.expression_v2 import Namespace
import numpy as np

def main(nelems: int, degree: int, maxrefine: int):
    '''
    Horizontally loaded linear elastic plate with FCM hole.

    .. arguments::

       nelems [2]
         Number of elements along edge.
       degree [1]
         Polynomial degree.
       maxrefine [1]
         Number or refinement levels used for the finite cell method.
    '''

    # reference value for area
    # triangle (0,0.9;0.9,0;0.9,0.9)
    measure_ref = 0.5*0.9**2 # 0.045

    # geometry 1: background domain [-0.1,1.6;-0.1,1.6]
    domain0, geom = mesh.rectilinear([np.linspace(-0.1, 1.6, nelems+1), np.linspace(-0.1, 1.6, nelems+1)])
    ns = Namespace()
    ns.x = geom
    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
    x,y = geom
    domain = domain0.trim(-(- (2**(1/2)*x)/2 - (2**(1/2)*(y - 9/10))/2), maxrefine=maxrefine, name='a1')
    domain = domain.trim(-(x - 9/10), maxrefine=maxrefine, name='a2')
    domain = domain.trim(-(y - 9/10), maxrefine=maxrefine, name='a3')
    gauss = domain.sample('gauss',degree+1)
    measure1 = gauss.integral('dV' @ ns)
    abserror1 = measure1-measure_ref
    print(measure1.eval())
    print(abserror1.eval())

    # geometry 2: background domain [0.0,1.6;0.0,1.6]
    domain0, geom = mesh.rectilinear([np.linspace(0.0, 1.6, nelems+1), np.linspace(0.0, 1.6, nelems+1)])
    ns = Namespace()
    ns.x = geom
    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
    x,y = geom
    domain = domain0.trim(-(- (2**(1/2)*x)/2 - (2**(1/2)*(y - 9/10))/2), maxrefine=maxrefine, name='a1')
    domain = domain.trim(-(x - 9/10), maxrefine=maxrefine, name='a2')
    domain = domain.trim(-(y - 9/10), maxrefine=maxrefine, name='a3')
    gauss = domain.sample('gauss',degree+1)
    measure2 = gauss.integral('dV' @ ns)
    abserror2 = measure2-measure_ref
    print(measure2.eval())
    print(abserror2.eval())
    
    # geometry 3: background domain [-0.5,1.6;-0.5,1.6]
    domain0, geom = mesh.rectilinear([np.linspace(-0.5, 1.6, nelems+1), np.linspace(-0.5, 1.6, nelems+1)])
    ns = Namespace()
    ns.x = geom
    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
    x,y = geom
    domain = domain0.trim(-(- (2**(1/2)*x)/2 - (2**(1/2)*(y - 9/10))/2), maxrefine=maxrefine, name='a1')
    domain = domain.trim(-(x - 9/10), maxrefine=maxrefine, name='a2')
    domain = domain.trim(-(y - 9/10), maxrefine=maxrefine, name='a3')
    gauss = domain.sample('gauss',degree+1)
    measure3 = gauss.integral('dV' @ ns)
    abserror3 = measure3-measure_ref
    print(measure3.eval())
    print(abserror3.eval())

The reference area is 0.045 which I obtain with the second background mesh definition. For the other two cases, I obtain 0.4034640609100461 and 0.406677201144048, respectively.

The geometry and the quadrature points for the first background mesh look as follows. Which seems to be fine, leaving it unclear to me whether there is a problem with the weights respectively the mapping.

Any help is appreciated.
Best, Michael Loibl