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.405 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

[gertjanvanzwieten]

Hi Michael, apologies for the late reply!
Your issue is the effect of the ndivisions argument of the trim method, which at its default value of 8 causes the edge intersections to be rounded to a uniform partition of 2^8=256 intervals. The reason for this rounding operation is twofold: 1. it helps to match edge intersections in neighbouring elements, and 2. it is a measure to discard nearly vanishing elements as they round down to zero. In general the rounding error is negligible against that of the linear interpolation of the levelset function, but this does not apply in your case where the interpolation is exact. What should help though is to increase ndivisions, for example to 32, which decreases the rounding error by a factor 6e8. Indeed, when I run your code with that modification I obtain the much improved values 0.4049999999084103 and 0.40500000009978976.
Hope this helps!

[MichaelLoibl]

Hi Gertjan,
thank you for the helpful explanation. I am also able to reproduce your answer. That definitely helped me to understand better what is actually going on.
Best,
Michael