Module regpy.vecsps.ngsolve

Finite element vector spaces using NGSolve

This module implements a VectorSpaceBase instance for NGSolve spaces. This gives the basic interface to use FES spaces defined in ngsolve to be used as VectorSpaceBasesin regpy.
Operators are using such spaces are implemented in the regpy.operators.ngsolve module. Hilbert spaces and Functionals defined on such spaces can be found in regpy.hilbert.ngsolve and regpy.functionals.ngsolve respectively.

Classes

class NgsBaseVector (vec: ngsolve.la.BaseVector, make_copy: bool | None = False)

Expand source code

@dataclass 
class NgsBaseVector:
    vec: ngs.la.BaseVector
    make_copy: Optional[bool] = field(default=False)

    __array_ufunc__ = None

    def copy(self):
        return copy(self)

    def __post_init__(self):
        if isinstance(self.vec,ngs.la.BaseVector):
            if self.make_copy:
                self.vec = deepcopy(self.vec)
            pass
        elif isinstance(self.vec,ngs.la.DynamicVectorExpression):
            if self.make_copy:
                self.vec = deepcopy(self.vec.Evaluate())
            else:
                self.vec = self.vec.Evaluate()
        else:
            raise TypeError("Could not treat {} type only ngs.la.BaseVector or ngs.la.DynamicVectorExpression".format(type(self.vec)))
        self.size = self.vec.size
        self.is_complex = self.vec.is_complex

    def conj(self):
        z = self.vec.CreateVector()
        for i in range(self.size): 
            z[i] = self.vec[i].real - 1j*self.vec[i].imag
        return NgsBaseVector(z)
    
    @property
    def real(self):
        z = self.vec.CreateVector()
        for i in range(self.size): 
            z[i] = self.vec[i].real
        return NgsBaseVector(z)
    
    @property
    def imag(self):
        z = self.vec.CreateVector()
        for i in range(self.size): 
            z[i] = self.vec[i].imag
        return NgsBaseVector(z)

    def __iadd__(self,other):
        assert isinstance(other,NgsBaseVector) and other.size == self.vec.size 
        self.vec.data += other.vec
        return self

    def __isub__(self,other):
        assert isinstance(other,NgsBaseVector) and other.size == self.vec.size 
        self.vec.data -= other.vec
        return self
    
    def __add__(self,other):
        assert isinstance(other,NgsBaseVector) and other.size == self.vec.size 
        v = self.vec.CreateVector()
        v.data = self.vec + other.vec
        return NgsBaseVector(v)
    
    def __radd__(self,other):
        return self + other
    
    def __sub__(self,other):
        return self + (-1*other)
    
    def __rsub__(self,other):
        return (-1*self) + other
    
    def __imul__(self,other):
        assert isinstance(other,float) or isinstance(other,int)
        self.vec.data *= other
        return self
    
    def __itruediv__(self,other):
        assert isinstance(other,float) or isinstance(other,int)
        self.vec.data /= other
        return self
    
    def __mul__(self,other):
        assert isinstance(other,float) or isinstance(other,int)
        v = self.vec.CreateVector()
        v.data = other * self.vec
        return NgsBaseVector(v)
        
    def __rmul__(self,other):
        return self * other

    def __truediv__(self,other):
        assert isinstance(other,float) or isinstance(other,int)
        v = self.vec.CreateVector()
        v.data = (1/other)*self.vec
        return NgsBaseVector(v)
    
    def __getitem__(self,i):
        return self.vec[i]
    
    def __setitem__(self,i,val):
        self.vec[i] = val

    def __iter__(self):
        return self.vec

    def iter_basis(self):
        v = NgsBaseVector(self.vec.CreateVector())
        for i in self.size:
            v[i] = 1 
            yield v
            if self.is_complex:
                v[i] = 1j
                yield v
            v[i] = 0                

    def __and__(self,x,y):
        assert isinstance(y,NgsBaseVector) and x.size == y.size
        return (x_i == y_i for x_i,y_i in zip(x,y))
    
    def __or__(self,x,y):
        assert isinstance(y,NgsBaseVector) and x.size == y.size
        return (x_i != y_i for x_i,y_i in zip(x,y))

    def __xor__(self,x,y):
        assert isinstance(y,NgsBaseVector) and x.size == y.size
        return (x_i != y_i for x_i,y_i in zip(x,y))
    
    def __copy__(self):
        return deepcopy(self)
    
    def __deepcopy__(self, memo):
        cls = self.__class__
        result = cls.__new__(cls)
        memo[id(self)] = result
        for k, v in self.__dict__.items():
            setattr(result, k, deepcopy(v, memo))
        return result

NgsBaseVector(vec: ngsolve.la.BaseVector, make_copy: Optional[bool] = False)

Instance variables

var vec : ngsolve.la.BaseVector

The type of the None singleton.

var make_copy : bool | None

The type of the None singleton.

prop real

Expand source code

@property
def real(self):
    z = self.vec.CreateVector()
    for i in range(self.size): 
        z[i] = self.vec[i].real
    return NgsBaseVector(z)

prop imag

Expand source code

@property
def imag(self):
    z = self.vec.CreateVector()
    for i in range(self.size): 
        z[i] = self.vec[i].imag
    return NgsBaseVector(z)

Methods

def copy(self)

Expand source code

def copy(self):
    return copy(self)

def conj(self)

Expand source code

def conj(self):
    z = self.vec.CreateVector()
    for i in range(self.size): 
        z[i] = self.vec[i].real - 1j*self.vec[i].imag
    return NgsBaseVector(z)

def iter_basis(self)

Expand source code

def iter_basis(self):
    v = NgsBaseVector(self.vec.CreateVector())
    for i in self.size:
        v[i] = 1 
        yield v
        if self.is_complex:
            v[i] = 1j
            yield v
        v[i] = 0                

class NgsVectorSpace (fes, bdr=None)

Expand source code

class NgsVectorSpace(VectorSpaceBase):
    """A vector space wrapping an `ngsolve.FESpace`.

    Parameters
    ----------
    fes : ngsolve.FESpace
       The wrapped NGSolve vector space.
    bdr : 
        Boundary of the NGSolve vector space.
    """

    log = classlogger

    def __init__(self, fes, bdr=None):
        assert isinstance(fes, ngs.FESpace)
        self.fes = fes
        self.bdr = bdr
        super().__init__(vec_type=NgsBaseVector, shape=(fes.ndof,), complex=fes.is_complex)
        # Checks if FES is Vector valued and stores the dimension in self.codim
        from netgen.libngpy._meshing import NgException
        try:
            self.codim = len(fes.components)
            assert self.codim == fes.mesh.dim
            self._fes_util = ngs.VectorL2(self.fes.mesh, order=0, complex = self.is_complex)
        except NgException:
            self.codim = 1
            self._fes_util = ngs.L2(self.fes.mesh, order=0, complex = self.is_complex)
        self._gfu_util = ngs.GridFunction(self._fes_util)
        self._gfu_fes = ngs.GridFunction(self.fes)
        self._help_x = NgsBaseVector(self._gfu_fes.vec)

    def zeros(self):
        h = self._gfu_fes.vec.CreateVector()
        h *= 0
        return NgsBaseVector(h,make_copy=True)
    
    def ones(self):
        self._gfu_fes.Set(1)
        return NgsBaseVector(ngs.Projector(self.fes.FreeDofs(), range=True).Project(self._gfu_fes.vec),make_copy=True)
    
    def empty(self):
        h = self._gfu_fes.vec.CreateVector()
        h *= 0
        return NgsBaseVector(h,make_copy=True)
    
    def rand(self,random_generator = None):
        random_generator = random_generator or np.random.random_sample 
        r = random_generator(self._fes_util.ndof)
        if self.is_complex and not is_complex_dtype(r.dtype):
            c = np.empty(self._fes_util.ndof, dtype=complex)
            c.real = r
            c.imag = random_generator(self._fes_util.ndof)
            self._gfu_util.vec.FV().NumPy()[:] = c            
        else:
            self._gfu_util.vec.FV().NumPy()[:] = r
        self._gfu_fes.Set(self._gfu_util)
        return NgsBaseVector(ngs.Projector(self.fes.FreeDofs(), range=True).Project(self._gfu_fes.vec),make_copy=True)
    
    def poisson(self,x):
        assert not self.is_complex
        self._gfu_fes.vec.FV().NumPy[:] =  np.random.poisson(x.vec.FV().NumPy())
        return NgsBaseVector(ngs.Projector(self.fes.FreeDofs(), range=True).Project(self._gfu_fes.vec),make_copy=True)

    def __contains__(self,x):
        if not isinstance(x,NgsBaseVector):
            return False
        elif x.size != self.fes.ndof:
            return False
        elif x.is_complex:
            return self.is_complex
        else:
            return True
        
    def vdot(self, x, y):
        return ngs.InnerProduct(x.vec,y.vec)

    def complex_space(self):
        if self.is_complex:
            return copy(self)
        return NgsVectorSpace(type(fes)(fes.mesh,order=fes.globalorder,bdr=self.bdr,complex=True),bdr=self.bdr)

    def real_space(self):
        if not self.is_complex:
            return copy(self)
        return NgsVectorSpace(type(fes)(fes.mesh,order=fes.globalorder,bdr=self.bdr,complex=False),bdr=self.bdr)
    
    def flatten(self, x:NgsBaseVector) -> np.ndarray:
        if self.is_complex:
            return np.concatenate([x.vec.FV().NumPy().real,x.vec.FV().NumPy().imag])
        else:
            return x.vec.FV().NumPy().copy()

    def fromflat(self, vec:np.ndarray) -> NgsBaseVector:
        if vec.ndim == 1:
            if self.is_complex and vec.size == self.shape[0] * 2:
                x = self.zeros()
                x.vec.FV().NumPy()[:] = vec[:self.shape[0]] + 1j*vec[self.shape[0]]
            elif vec.size == self.shape[0]:
                x = self.zeros()
                x.vec.FV().NumPy()[:] = vec
            else:
                raise ValueError("provided vector has non fitting shape.")
        else:
            raise ValueError("Provided vector must be one dimensional")
        return x
    
    def __eq__(self, other: object) -> bool:
        if not isinstance(other,type(self)):
            return False
        return self.fes == other.fes
    
    def iter_basis(self):
        r"""Generator iterating over the standard basis of the vector space. For efficiency,
        the same array is returned in each step, and subsequently modified in-place. If you need
        the array longer than that, perform a copy. In case of complex a vector space after each
        each array modified in its place with a real one it returns the same vector with \(1i\)
        in its place.   
        """
        elm = self.zeros()
        for idx in range(self.shape[0]):
            elm[idx] = 1
            yield elm
            if self.is_complex:
                elm[idx] = 1j
                yield elm
            elm[idx] = 0

    def logical_and(self,x,y):
        return x & y 
    
    def logical_or(self,x,y):
        return x | y 
    
    def logical_not(self,x):
        return not x
    
    def logical_xor(self,x,y):
        return x ^ y

    def is_on_boundary(self,x):
        if self.bdr is None:
            return False
        t = x.vec.CreateVector()
        t.data = x.vec
        ngs.Projector(self.fes.FreeDofs(), range=True).Project(t)
        return np.all(t.FV().NumPy() == 0)
    
    def to_gf(self, x):
        gf = ngs.GridFunction(self.fes)
        gf.vec.data = x.vec
        return gf
    
    def from_ngs(self, ngs_elem, definedon : ngs.comp.Region|None = None):
        if isinstance(ngs_elem,ngs.comp.GridFunction):
            return NgsBaseVector(ngs_elem.vec,make_copy=True)
        else:
            self._gfu_fes.Set(ngs_elem,definedon=definedon)
            return self._help_x.copy()

A vector space wrapping an ngsolve.FESpace.

Parameters

fes : ngsolve.FESpace

 

The wrapped NGSolve vector space.

bdr

Boundary of the NGSolve vector space.

Ancestors

• VectorSpaceBase

Methods

def iter_basis(self)

Expand source code

def iter_basis(self):
    r"""Generator iterating over the standard basis of the vector space. For efficiency,
    the same array is returned in each step, and subsequently modified in-place. If you need
    the array longer than that, perform a copy. In case of complex a vector space after each
    each array modified in its place with a real one it returns the same vector with \(1i\)
    in its place.   
    """
    elm = self.zeros()
    for idx in range(self.shape[0]):
        elm[idx] = 1
        yield elm
        if self.is_complex:
            elm[idx] = 1j
            yield elm
        elm[idx] = 0

Generator iterating over the standard basis of the vector space. For efficiency, the same array is returned in each step, and subsequently modified in-place. If you need the array longer than that, perform a copy. In case of complex a vector space after each each array modified in its place with a real one it returns the same vector with $1i$ in its place.

def is_on_boundary(self, x)

Expand source code

def is_on_boundary(self,x):
    if self.bdr is None:
        return False
    t = x.vec.CreateVector()
    t.data = x.vec
    ngs.Projector(self.fes.FreeDofs(), range=True).Project(t)
    return np.all(t.FV().NumPy() == 0)

def to_gf(self, x)

Expand source code

def to_gf(self, x):
    gf = ngs.GridFunction(self.fes)
    gf.vec.data = x.vec
    return gf

def from_ngs(self, ngs_elem, definedon: ngsolve.comp.Region | None = None)

Expand source code

def from_ngs(self, ngs_elem, definedon : ngs.comp.Region|None = None):
    if isinstance(ngs_elem,ngs.comp.GridFunction):
        return NgsBaseVector(ngs_elem.vec,make_copy=True)
    else:
        self._gfu_fes.Set(ngs_elem,definedon=definedon)
        return self._help_x.copy()

Inherited members

• VectorSpaceBase:
  • complex_space
  • empty
  • flatten
  • fromflat
  • identity
  • is_complex
  • log
  • logical_and
  • logical_not
  • logical_or
  • logical_xor
  • masked_space
  • ndim
  • ones
  • poisson
  • rand
  • randn
  • real_space
  • realsize
  • shape
  • size
  • vdot
  • vec_type
  • zeros