<commit_before># -*- coding: utf8 -*-
#
# Module ELEMENT
#
# Part of Nutils: open source numerical utilities for Python. Jointly developed
# by HvZ Computational Engineering, TU/e Multiscale Engineering Fluid Dynamics,
# and others. More info at http://nutils.org <<EMAIL>>. (c) 2014

"""
The element module defines reference elements such as the :class:`QuadElement`
and :class:`TriangularElement`, but also more exotic objects like the
:class:`TrimmedElement`. A set of (interconnected) elements together form a
:mod:`nutils.topology`. Elements have edges and children (for refinement), which
are in turn elements and map onto self by an affine transformation. They also
have a well defined reference coordinate system, and provide pointsets for
purposes of integration and sampling.
"""

from . import log, util, numpy, core, numeric, function, cache, transform, rational, _
import re, warnings


## ELEMENT

class Element( object ):

  __slots__ = 'transform', 'reference', 'opposite'

  def __init__( self, reference, transform, opposite=None ):
    assert transform.fromdims == reference.ndims
    self.reference = reference
    self.transform = transform
    self.opposite = opposite or transform

  def __eq__( self, other ):
    return self is other or isinstance(other,Element) \
      and self.reference == other.reference \
      and self.transform == other.transform \
      and self.opposite == other.opposite

  @property
  def vertices( self ):
    return self.transform.apply( self.reference.vertices )

  @property
  def ndims( self ):
    return self.reference.ndims

  @property
  def nverts( self ):
    return self.reference.nverts

  @property
  def edges( self ):
    return [ self.edge(i) for i in range(self.nverts) ]
    
  def edge( self, iedge ):
    trans, edge = self.reference.edges[iedge]
    return Element( edge, trans >> self.transform, trans >> self.opposite )

  @property
  def children( self ):
    return [ Element( child, trans >> self.transform, trans >> self.opposite )
      for trans, child in self.reference.children ]

  def trim( self, levelset, maxrefine, numer ):
    'trim element along levelset'

    assert rational.isrational( numer )
    pos = []
    neg = []
    if maxrefine <= 0:
      repeat = True
      levels = levelset.eval( self, 'bezier2' )
      while repeat: # set almost-zero points to zero if cutoff within eps
        repeat = False
        assert levels.shape == (self.reference.nverts,)
        if numpy.greater_equal( levels, 0 ).all():
          return self, None
        if numpy.less_equal( levels, 0 ).all():
          return None, self
        isects = []
        for ribbon in self.reference.ribbon2vertices:
          a, b = levels[ribbon]
          if a * b < 0: # strict sign change
            x = int( numer * a / float(a-b) + .5 ) # round to [0,1,..,numer]
            if 0 < x < numer:
              isects.append(( x, ribbon ))
            else: # near intersection of vertex
              v = ribbon[ (0,numer).index(x) ]
              log.debug( 'rounding vertex #%d from %f to 0' % ( v, levels[v] ) )
              levels[v] = 0
              repeat = True
      coords = self.reference.vertices
      if isects:
        coords = numpy.vstack([
          self.reference.vertices * int(numer),
          [ numpy.dot( (int(numer)-x,x), self.reference.vertices[ribbon] ) for x, ribbon in isects ]
        ])
      assert coords.dtype == int
      simplex = SimplexReference( self.ndims )
      for tri in util.delaunay( coords ):
        ispos = isneg = False
        for ivert in tri:
          if ivert < self.reference.nverts:
            sign = levels[ivert]
            ispos = ispos or sign > 0
            isneg = isneg or sign < 0
        assert ispos is not isneg, 'domains do not separate in two convex parts'
        offset = coords[tri[0]]
        matrix = ( coords[tri[1:]] - offset ).T
        if numpy.linalg.det( matrix.astype(float) ) < 0:
          tri[-2:] = tri[-1], tri[-2]
          matrix = ( coords[tri[1:]] - offset ).T
        trans = transform.linear(matrix,numer) >> transform.shift(offset,numer)
        ( pos if ispos else neg ).append(( trans, simplex ))
    else:
      for trans, child in self.reference.children:
        child = Element( child, trans >> self.transform, trans >> self.opposite )
        poschild, negchild = child.trim( levelset, maxrefine-1, numer )
        if poschild:
          pos.append( (trans,poschild.reference) )
        if negchild:
          neg.append( (trans,negchild.reference) )
    if not neg:
      return self, None
    if not pos:
      return None, self
    return Element( MosaicReference( self.ndims, tuple(pos) ), self.transform, self.opposite ), \
           Element( MosaicReference( self.ndims, tuple(neg) ), self.transform, self.opposite )

  @property
  def simplices( self ):
    return [ Element( reference, trans >> self.transform, trans >> self.opposite )
      for trans, reference in self.reference.simplices ]

  def __str__( self ):
    return 'Element(%s)' % self.vertices


## REFERENCE ELEMENTS

class Reference( object ):
  'reference element'

  __metaclass__ = cache.Meta

  def __init__( self, vertices ):
    self.vertices = numpy.asarray( vertices )
    assert self.vertices.dtype == int
    self.nverts, self.ndims = self.vertices.shape

  @property
  def simplices( self ):
    return [ (transform.identity(self.ndims),self) ]

  def getischeme( self, ischeme ):
    if self.ndims == 0:
      return numpy.zeros([1,0]), numpy.array([1.])
    match = re.match( '([a-zA-Z]+)(.+)', ischeme )
    assert match, 'cannot parse integration scheme %r' % ischeme
    ptype, args = match.groups()
    get = getattr( self, 'getischeme_'+ptype )
    return get( eval(args) ) if args else get()

  def __mul__( self, other ):
    assert isinstance( other, Reference )
    return other if self.ndims == 0 \
      else self if other.ndims == 0 \
      else TensorReference( self, other )

  def __pow__( self, n ):
    assert isinstance( n, int ) and n >= 0
    return SimplexReference(0) if n == 0 \
      else self if n == 1 \
      else self * self**(n-1)

class SimplexReference( Reference ):

  def __init__( self, ndims ):
    assert ndims >= 0
    vertices = numpy.concatenate( [ numpy.zeros(ndims,dtype=int)[_,:],
                                    numpy.eye(ndims,dtype=int) ], axis=0 )
    Reference.__init__( self, vertices )
    self._bernsteincache = [] # TEMPORARY

  def stdfunc( self, degree ):
    if self.ndims == 1:
      if len(self._bernsteincache) <= degree or self._bernsteincache[degree] is None:
        self._bernsteincache += [None] * (degree-len(self._bernsteincache))
        self._bernsteincache.append( PolyLine( PolyLine.bernstein_poly(degree) ) )
      return self._bernsteincache[degree]
    if self.ndims == 2:
      return PolyTriangle(degree)
    raise NotImplementedError

  @cache.property
  def ribbon2vertices( self ):
    return numpy.array([ (i,j) for i in range( self.ndims+1 ) for j in range( i+1, self.ndims+1 ) ])

  def getischeme_contour( self, n ):
    assert self.ndims == 2
    p = numpy.arange( n+1, dtype=float ) / (n+1)
    z = numpy.zeros_like( p )
    return numpy.hstack(( [1-p,p], [z,1-p], [p,z] )).T, None<commit_msg>move trim algorithm into reference<commit_after># -*- coding: utf8 -*-
#
# Module ELEMENT
#
# Part of Nutils: open source numerical utilities for Python. Jointly developed
# by HvZ Computational Engineering, TU/e Multiscale Engineering Fluid Dynamics,
# and others. More info at http://nutils.org <<EMAIL>>. (c) 2014

"""
The element module defines reference elements such as the :class:`QuadElement`
and :class:`TriangularElement`, but also more exotic objects like the
:class:`TrimmedElement`. A set of (interconnected) elements together form a
:mod:`nutils.topology`. Elements have edges and children (for refinement), which
are in turn elements and map onto self by an affine transformation. They also
have a well defined reference coordinate system, and provide pointsets for
purposes of integration and sampling.
"""

from . import log, util, numpy, core, numeric, function, cache, transform, rational, _
import re, warnings


## ELEMENT

class Element( object ):

  __slots__ = 'transform', 'reference', 'opposite'

  def __init__( self, reference, transform, opposite=None ):
    assert transform.fromdims == reference.ndims
    self.reference = reference
    self.transform = transform
    self.opposite = opposite or transform

  def __eq__( self, other ):
    return self is other or isinstance(other,Element) \
      and self.reference == other.reference \
      and self.transform == other.transform \
      and self.opposite == other.opposite

  @property
  def vertices( self ):
    return self.transform.apply( self.reference.vertices )

  @property
  def ndims( self ):
    return self.reference.ndims

  @property
  def nverts( self ):
    return self.reference.nverts

  @property
  def edges( self ):
    return [ self.edge(i) for i in range(self.nverts) ]
    
  def edge( self, iedge ):
    trans, edge = self.reference.edges[iedge]
    return Element( edge, trans >> self.transform, trans >> self.opposite )

  @property
  def children( self ):
    return [ Element( child, trans >> self.transform, trans >> self.opposite )
      for trans, child in self.reference.children ]

  def trim( self, levelset, maxrefine, numer ):
    'trim element along levelset'

    pos, neg = self.reference.trim( self.transform, levelset, maxrefine, numer )
    if not neg:
      return self, None
    if not pos:
      return None, self
    return Element( pos, self.transform, self.opposite ), \
           Element( neg, self.transform, self.opposite )

  @property
  def simplices( self ):
    return [ Element( reference, trans >> self.transform, trans >> self.opposite )
      for trans, reference in self.reference.simplices ]

  def __str__( self ):
    return 'Element(%s)' % self.vertices


## REFERENCE ELEMENTS

class Reference( object ):
  'reference element'

  __metaclass__ = cache.Meta

  def __init__( self, vertices ):
    self.vertices = numpy.asarray( vertices )
    assert self.vertices.dtype == int
    self.nverts, self.ndims = self.vertices.shape

  @property
  def simplices( self ):
    return [ (transform.identity(self.ndims),self) ]

  def getischeme( self, ischeme ):
    if self.ndims == 0:
      return numpy.zeros([1,0]), numpy.array([1.])
    match = re.match( '([a-zA-Z]+)(.+)', ischeme )
    assert match, 'cannot parse integration scheme %r' % ischeme
    ptype, args = match.groups()
    get = getattr( self, 'getischeme_'+ptype )
    return get( eval(args) ) if args else get()

  def __mul__( self, other ):
    assert isinstance( other, Reference )
    return other if self.ndims == 0 \
      else self if other.ndims == 0 \
      else TensorReference( self, other )

  def __pow__( self, n ):
    assert isinstance( n, int ) and n >= 0
    return SimplexReference(0) if n == 0 \
      else self if n == 1 \
      else self * self**(n-1)

  def trim( self, trans, levelset, maxrefine, numer ):
    'trim element along levelset'

    assert rational.isrational( numer )
    pos = []
    neg = []
    if maxrefine <= 0:
      repeat = True
      levels = levelset.eval( Element(self,trans), 'bezier2' )
      while repeat: # set almost-zero points to zero if cutoff within eps
        repeat = False
        assert levels.shape == (self.nverts,)
        if numpy.greater_equal( levels, 0 ).all():
          return self, None
        if numpy.less_equal( levels, 0 ).all():
          return None, self
        isects = []
        for ribbon in self.ribbon2vertices:
          a, b = levels[ribbon]
          if a * b < 0: # strict sign change
            x = int( numer * a / float(a-b) + .5 ) # round to [0,1,..,numer]
            if 0 < x < numer:
              isects.append(( x, ribbon ))
            else: # near intersection of vertex
              v = ribbon[ (0,numer).index(x) ]
              log.debug( 'rounding vertex #%d from %f to 0' % ( v, levels[v] ) )
              levels[v] = 0
              repeat = True
      coords = self.vertices
      if isects:
        coords = numpy.vstack([
          self.vertices * int(numer),
          [ numpy.dot( (int(numer)-x,x), self.vertices[ribbon] ) for x, ribbon in isects ]
        ])
      assert coords.dtype == int
      simplex = SimplexReference( self.ndims )
      for tri in util.delaunay( coords ):
        ispos = isneg = False
        for ivert in tri:
          if ivert < self.nverts:
            sign = levels[ivert]
            ispos = ispos or sign > 0
            isneg = isneg or sign < 0
        assert ispos is not isneg, 'domains do not separate in two convex parts'
        offset = coords[tri[0]]
        matrix = ( coords[tri[1:]] - offset ).T
        if numpy.linalg.det( matrix.astype(float) ) < 0:
          tri[-2:] = tri[-1], tri[-2]
          matrix = ( coords[tri[1:]] - offset ).T
        trans = transform.linear(matrix,numer) >> transform.shift(offset,numer)
        ( pos if ispos else neg ).append(( trans, simplex ))
    else:
      for ctrans, child in self.children:
        poschild, negchild = child.trim( ctrans >> trans, levelset, maxrefine-1, numer )
        if poschild:
          pos.append( (ctrans,poschild) )
        if negchild:
          neg.append( (ctrans,negchild) )
    if not neg:
      return self, None
    if not pos:
      return None, self
    return MosaicReference( self.ndims, tuple(pos) ), \
           MosaicReference( self.ndims, tuple(neg) )

class SimplexReference( Reference ):

  def __init__( self, ndims ):
    assert ndims >= 0
    vertices = numpy.concatenate( [ numpy.zeros(ndims,dtype=int)[_,:],
                                    numpy.eye(ndims,dtype=int) ], axis=0 )
    Reference.__init__( self, vertices )
    self._bernsteincache = [] # TEMPORARY

  def stdfunc( self, degree ):
    if self.ndims == 1:
      if len(self._bernsteincache) <= degree or self._bernsteincache[degree] is None:
        self._bernsteincache += [None] * (degree-len(self._bernsteincache))
        self._bernsteincache.append( PolyLine( PolyLine.bernstein_poly(degree) ) )
      return self._bernsteincache[degree]
    if self.ndims == 2:
      return PolyTriangle(degree)
    raise NotImplementedError

  @cache.property
  def ribbon2vertices( self ):
    return numpy.array([ (i,j) for i in range( self.ndims+1 ) for j in range( i+1, self.ndims+1 ) ])

  def getischeme_contour( self, n ):
    assert self.ndims == 2
    p = numpy.arange( n+1, dtype=float ) / (n+1)
    z = numpy.zeros_like( p )
    return numpy.hstack(( [1-p,p], [z,1-p], [p,z] )).T, None
{'fraction_non_alphanumeric': 0.11733824231267471, 'fraction_numerical': 0.007090747937546874, 'mean_word_length': 2.650572424091588, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4778123', 'n_tokens_mistral': 4608, 'n_tokens_neox': 4477, 'n_words': 1413}