#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#pylint: disable=C0302
"""
Created on Mon Dec 4 23:30:08 2023
@author: nchevaug
"""
#from warnings import warn
#from __future__ import annotations
from abc import ABC
import unittest
from typing import Union
import numpy as np
import numpy.typing as npt
import scipy as sp
import pyxfem.core.level_set_mesh_cut as lscut
import pyxfem.core.mapping as mapping
from pyxfem.core.mesh_simplex_one_level import SMesh
from pyxfem.core.mesh_modifications import MeshHierarchicModification as MeshMod
import pyxfem.core.quadrature as quad
[docs]
class DofManagerBase(ABC):
''' Base class to help manage degree of freedoms '''
[docs]
def scalar_size(self):
""" size (nbdofs) for a scalar space """
raise NotImplementedError()
[docs]
def size(self):
""" size (nbdofs) taking into account the vector dim of the field """
raise NotImplementedError()
[docs]
def vertexid2dofid(self, vids):
""" return the dofs number associated to the vertices id """
raise NotImplementedError()
[docs]
def edge2edgeid(self, e2v, getxy=False):
""" from edge defined by there vertices id, recover the edgeid as seen by the space """
raise NotImplementedError()
[docs]
def edgeid2dofid(self, edges):
""" from edge id recover the associated dofid
For FEMSpaceP1, return an empty array
"""
raise NotImplementedError()
[docs]
class VertexSpaceDofManager():
""" Class to compute the degree of freedom associated to the mesh and the dirichlet BC """
# pylint: disable=C0103
def __init__(self, m:SMesh, vecdim = 3):
self.m = m
self.vid = self.m.used_v() # np.unique(tris.flatten())
self.vecdim = vecdim
self.interpolatoryAt = 'vertex'
if len(self.vid) != len(self.m.xy):
self.allnodes = False
self.vid2dofid = dict(zip(self.vid, np.arange(len(self.vid))))
self.a2dofs = np.array([[self.vid2dofid[vid] for vid in t]
for t in self.m.t2v()])
else:
self.allnodes = True
self.vid2dofid = np.arange(len(self.vid))
self.a2dofs = self.vid2dofid[self.m.t2v()]
self.nelem = len(self.a2dofs)
self.nscaldof = len(self.vid)
[docs]
def vertexid2dofid(self, vids):
''' return an array A of shape (nv, vecdim) containing such as A[i, j] is the de
index of degree of freedom associated to node i in vector direction j
'''
vid_space = np.array([self.vid2dofid[vid] for vid in vids])
return ((self.vecdim*vid_space)[:, np.newaxis] + np.arange(self.vecdim)).squeeze()
[docs]
def size(self):
""" Return the total size of the space (number of (scalar) dof in the space )"""
return len(self.vid)*self.vecdim
[docs]
def op2sparse(NScalar, nscalar_dofs, a2scalardofs, vecdim, #pylint: disable=C0103
nshapeperelem, nptperelem, ndataperpoint, only_on=None):
""" Given a N array (shape function values at integration point for each element
and numbering of degree of freedom dat, construct the sparse N Operator (Assembly)
"""
# nd dataperpoint : not a very good name.
# my retionnal is this : vecDim represent the number of scalar dof per node.
# For a scalar field it would be 1, For a 2d vector field, it would be 2,
# for a 3 d vector Field it would be 2. IT is also possible to have a scalar field, say T
# and a 3d vector field say, u stored in the same space. the vec dim would then be 4.
# (The function is meant to assemble a multid field, but with the same space for each dir)
# The nb data perpoint, is the number of value per vec dim per interpolation point.
# For the evaluation of the field, it would be 1, for it's gradient dNdu,
# it would be 2 on aface, 1 on an edge and 3 on a 3d element.
nelem = len(a2scalardofs)
neval = len(NScalar)
npts = nelem*nptperelem
nrow = npts*ndataperpoint*vecdim
ncol = nscalar_dofs*vecdim
shape = (nrow, ncol)
if only_on is None:
col = np.repeat(a2scalardofs, nptperelem*ndataperpoint, axis=0)
row_ptr = np.arange(0, (nrow+1)*nshapeperelem, nshapeperelem)
else:
col = np.repeat(a2scalardofs[only_on], nptperelem*ndataperpoint, axis=0)
row = np.repeat(np.arange(nrow).reshape((nelem, -1))[only_on].flatten(), nshapeperelem)
row_ptr = np.cumsum(np.bincount(row+1, minlength=nrow+1))
if vecdim > 1:
col = col.reshape((-1, nshapeperelem*ndataperpoint))
col = np.column_stack([vecdim*col + i for i in range(vecdim)]).flatten()
rshape = (neval, nptperelem, ndataperpoint, nshapeperelem)
N = np.tile(NScalar.reshape(rshape), [1,1,vecdim,1]).flatten() #pylint: disable=C0103
else:
col = col.flatten()
N = NScalar.flatten() #pylint: disable=C0103
N = sp.sparse.csr_array((N, col, row_ptr), shape) #pylint: disable=C0103
N.sort_indices()
return N
[docs]
def op2sparse_rmap(NScalar, nscalar_dofs, a2scalardofs, vecdim, #pylint: disable=C0103
nshapeperelem, nptperelem,
ndataperpoint, rmap, frmap=None, sort=True):
""" Given a N array (shape function values at integration point for each element
(or sub element)) and numbering of degree of freedom dat, construct
the sparse N Operator (Assembly)
Specialized version for the case (typically xfem) where element are cut in sub cell
"""
#nelem = len(a2ScalarDofs)
#neval = len(NScalar)
#npts = nelem*nptperelem
#shape = (npts*ndataperpoint*vecDim, nScalarDofs*vecDim)
npts = rmap.nbevals*nptperelem
nscalrow = npts*ndataperpoint # nb row for a scalar field
nrow = nscalrow*vecdim
ncol = nscalar_dofs*vecdim
shape = (nrow, ncol)
if frmap is None:
col = np.repeat(a2scalardofs, rmap.appro2nbevals*nptperelem*ndataperpoint, axis=0)
row_ptr = np.arange(0, (nrow+1)*nshapeperelem, nshapeperelem)
if vecdim > 1:
col = col.reshape((-1, nshapeperelem*ndataperpoint))
col = np.column_stack([vecdim*col + i for i in range(vecdim)]).flatten()
rshape = (rmap.nbevals, nptperelem, ndataperpoint, nshapeperelem)
N = np.tile(NScalar.reshape(rshape), [1,1,vecdim,1]).flatten() #pylint: disable=C0103
else:
col = col.flatten()
N = NScalar.flatten() #pylint: disable=C0103
else:
fscalcol = a2scalardofs[frmap.only_on_parents]
col = np.repeat(fscalcol, frmap.appro2nbevals*nptperelem*ndataperpoint, axis=0)
row = np.arange(rmap.nbevals*nptperelem*ndataperpoint).reshape((rmap.nbevals, -1))
row = row[np.isin(rmap.eval2appro, frmap.only_on_parents)].flatten()
if vecdim > 1:
row = np.column_stack(
tuple(vecdim*row + i for i in range(vecdim))).flatten()
col = col.reshape((-1, nshapeperelem*ndataperpoint))
col = np.column_stack([vecdim*col + i for i in range(vecdim)]).flatten()
rshape = (frmap.nbevals, nptperelem, ndataperpoint, nshapeperelem)
N = np.tile(NScalar.reshape(rshape), [1,1,vecdim,1]).flatten() #pylint: disable=C0103
#warn("op2sparse_rmap Not fully tested for case vecDim > 1 and frmap != None")
else:
row = row.flatten()
col = col.flatten()
N = NScalar.flatten() #pylint: disable=C0103
row_ptr = np.cumsum(nshapeperelem*np.bincount(row+1, minlength=nrow+1))
N = sp.sparse.csr_array((N, col, row_ptr), shape) #pylint: disable=C0103
if sort:
N.sort_indices()
return N
[docs]
def meshcutdata2remap_helper(mesh_cut, hmin=0., alpha=0.1, p=1):
""" build a RemapHelper Object from mesh_cut data """
appro_mapping = mapping.MappingT3(
mesh_cut.xy[mesh_cut.tris], hmin=hmin, alpha=alpha, p=p)
appro2evals_index, appro2evals = mesh_cut.parents2sons()
eval_tris_xy = mesh_cut.xy_new_vertices[mesh_cut.new_tris]
eval_mapping = mapping.MappingT3(eval_tris_xy, hmin=hmin, alpha=alpha, p=p)
return mapping.RemapHelper(appro2evals_index, appro2evals, appro_mapping, eval_mapping)
[docs]
def operator_trival2edge(nedge, t2e, t2eflip, ndataperedge=None, datasize=1, dtype=np.float64):
""" return linear operator that push values at triangles to each edge. useful to compute
interelemnt jump for example """
ntri = len(t2e)
if ndataperedge is None:
push2edge_data = np.ones(ntri*datasize*3, dtype=dtype)
push2edge_row = np.column_stack( [(2*datasize*t2e.flatten()+i \
+ (datasize)*t2eflip.flatten()) for i in range(datasize)]).flatten()
push2edge_col = np.repeat(np.arange(ntri*datasize).reshape((ntri,-1)), 3, axis = 0) \
.flatten()
shape = (2*nedge*datasize, ntri*datasize)
else:
push2edge_data = np.ones(ntri*3*ndataperedge*datasize, dtype=dtype)
push2edge_col = np.arange(ntri*3*ndataperedge*datasize)
noflip = 2*np.arange(ndataperedge *
datasize).reshape((ndataperedge, datasize))
doflip = noflip[::-1]+1
push2edge_row = (ndataperedge*datasize*2*t2e[:, :, np.newaxis, np.newaxis] +
np.where((t2eflip == 0)[:, :, np.newaxis, np.newaxis],
noflip[np.newaxis, np.newaxis, :, :],
doflip[np.newaxis, np.newaxis, :, :])
).flatten()
shape = (nedge*ndataperedge*datasize*2, ntri*ndataperedge*3*datasize)
return sp.sparse.csr_array((push2edge_data, (push2edge_row, push2edge_col)), shape=shape)
[docs]
def operator_dof2value(space, vecdim, scalar_space_size, tri2dofs, evaltrisuv,
rmap=None, frmap=None, raw_data=False, only_on = None):
""" return a linear operator N that map dof value to value at interpolation point of each
element """
npts = len(evaltrisuv)
nshape = space.lenN
Neval = space.N(evaltrisuv) #pylint: disable=C0103
N = None #pylint: disable=C0103
if rmap is None and frmap is None:
if only_on is None:
N = np.broadcast_to(Neval, (len(tri2dofs), npts, nshape)) #pylint: disable=C0103
else:
N = np.broadcast_to(Neval, (len(only_on), npts, nshape)) #pylint: disable=C0103
if not raw_data:
N = op2sparse(N, scalar_space_size, tri2dofs, vecdim, #pylint: disable=C0103
nshape, npts, 1, only_on=only_on)
if rmap is not None:
rmap1 = rmap if frmap is None else frmap
sons_point_on_parents_uv = rmap1.remap_on_parent(evaltrisuv)
N = np.zeros((rmap1.nbevals, npts, nshape)) #pylint: disable=C0103
N[rmap1.noson2eval_index] = Neval
N[rmap1.son2eval_index] = space.N(sons_point_on_parents_uv)
if not raw_data:
N = op2sparse_rmap( #pylint: disable=C0103
N, scalar_space_size, tri2dofs, vecdim, nshape, npts, 1, rmap, frmap)
if N is None:
raise ValueError("N was not constructed due to wrong input")
return N
[docs]
def gradNTri(dNdu, uv, invF): #pylint: disable=C0103
''' Compute the gradient with regard to space coordinates xy from gradient with regard to
element coordinate uv'''
dNdx = np.einsum('...ij,...ik->...jk', invF, dNdu) #pylint: disable=C0103
# map is constant and DNdU is constant.
# repeat the value as much as the number of points per elem
return np.repeat(dNdx[:, np.newaxis, :, :], uv.shape[-2], axis=1)
[docs]
def operator_dof2grad_tri(space, vecdim, geomdim, scalar_space_size, tri2dofs, trimap, evaltrisuv,
rmap=None, frmap=None, only_on=None, raw_data=False, limited=False):
""" helper function : return operator N that map the dof value to the gradient of field values
at each point uv for each triangle """
npts = len(evaltrisuv)
nshape = space.lenN
dNdx = None #pylint: disable=C0103
invF = trimap.inv_F if not limited else trimap.invFLimited #pylint: disable=C0103
points_uv = np.zeros((0, npts, 0))
if rmap is None and frmap is None:
invF = invF if only_on is None else invF[only_on] #pylint: disable=C0103
dNdx = gradNTri(space.dNdu, points_uv, invF) #pylint: disable=C0103
if not raw_data:
dNdx = op2sparse(dNdx, scalar_space_size, tri2dofs, #pylint: disable=C0103
vecdim, nshape, npts, geomdim, only_on=only_on)
elif vecdim!=1:
neval = len(dNdx)
nptperelem = npts
ndataperpoint = geomdim
nshapeperelem = nshape
rshape = (neval, nptperelem, ndataperpoint, nshapeperelem)
dNdx.reshape(rshape)
rshapev = (neval, nptperelem, ndataperpoint*vecdim, nshapeperelem*vecdim)
dNdxv = np.zeros(rshapev) #pylint: disable=C0103
for r in range(vecdim):
rs = ndataperpoint*r
re = ndataperpoint*(r+1)
rows = np.arange(rs,re,1)
cs = r
ce = nshapeperelem*vecdim
cols = np.arange(cs, ce, vecdim)
dNdxv[..., rows[:, np.newaxis], cols] = dNdx
dNdx = dNdxv #pylint: disable=C0103
if rmap is not None:
# Warning : all coordinate are put to zeros because here we don't care :
# dNdu is constant on the element
son_on_parent_uv = np.zeros((0, npts, 0))
parent_uv = np.zeros((0, npts, 0))
rmap1 = rmap if frmap is None else frmap
dNdu = space.dNdu #pylint: disable=C0103
dNdx = np.zeros((rmap1.nbevals, npts, geomdim, nshape)) #pylint: disable=C0103
dNdx[rmap1.noson2eval_index] = gradNTri(dNdu, parent_uv, invF[rmap1.approEQeval])
dNdx[rmap1.son2eval_index] = gradNTri(dNdu, son_on_parent_uv,
invF[rmap1.son2parent_elem_ids])
if not raw_data:
dNdx = op2sparse_rmap( #pylint: disable=C0103
dNdx, scalar_space_size, tri2dofs, vecdim, nshape, npts, geomdim, rmap, frmap)
if dNdx is None:
raise ValueError("dNdx was not constructed due to wrong input")
return dNdx
[docs]
class TestSparseOp(unittest.TestCase):
""" Class for unit test of operator on dofs """
[docs]
def test_operator_trival2edge(self):
""" Test of the operator tha push triangles data to edges """
xy = np.array([[0.,0.], [1.,0.], [0.,1.], [1.,1.]])
tris = np.array([[0,1,2], [1,3,2]])
mesh = SMesh(xy, tris)
t2e, t2ef = mesh.t2e_flip()
nedges = mesh.nedges()
ntris = mesh.ntris()
# Testing operator_trival2edges, for datasize = 1
opt2e = operator_trival2edge(nedges, t2e, t2ef)
valtri = np.arange(1, ntris+1)
valedge = (opt2e@valtri).reshape((-1,2))
valedgeref = np.array([1., 0.,0.,1.,1.,2.,2.,0., 0.,2.]).reshape((-1,2))
err = float(np.linalg.norm(valedge - valedgeref))
self.assertAlmostEqual(err, 0.,\
msg = "test_operator_trival2edge do not return the expected result for case datasize=1")
# Testing operator_trival2edges, for datasize = 3
datasize = (3,)
opt2e = operator_trival2edge(nedges, t2e, t2ef, datasize= datasize[0])
valtri = np.arange(1, ntris*np.prod(datasize)+1).reshape((-1,) + datasize)
valedge = (opt2e@valtri.flatten()).reshape((-1,2) + datasize)
valedgeref = np.array([1.,2.,3.,0.,0.,0., 0.,0.,0.,1.,2.,3., 1.,2.,3.,4.,5.,6., \
4.,5.,6., 0.,0.,0., 0.,0.,0., 4.,5.,6.]).reshape((-1,2) + datasize)
err = float(np.linalg.norm(valedge - valedgeref))
self.assertAlmostEqual(err, 0.,\
msg = "test_operator_trival2edge do not return the expected result for case datasize=3")
[docs]
class FEMSpaceP1():
""" Classical Linear Finite Element space """
interpolatoryAt = 'vertex'
dNdu = np.array([[-1., 1., 0.], [-1., 0., 1.]])
lenN = 3
[docs]
@staticmethod
def N(uv): #pylint: disable=C0103
""" Shape function for the triangle """
u = uv[..., 0]
v = uv[..., 1]
return np.stack((1.-u-v, u, v), axis=u.ndim)
def __init__(self, name:str, mesh:SMesh,
vecdim:int=1, tris_map:None|mapping.MappingT3 = None) -> None:
self.name = name
self.m = mesh
self.vecdim = vecdim
self._tris_map = tris_map
self.vid = self.m.used_v() # np.unique(tris.flatten())
print(len(self.vid) , len(self.m.vcoord()))
if len(self.vid) != len(self.m.vcoord()):
self.allnodes = False
self.vid2dofid = dict(zip(self.vid, np.arange(len(self.vid))))
self.a2dofs = np.array([[self.vid2dofid[vid] for vid in t]
for t in self.m.t2v()])
else:
self.allnodes = True
self.vid2dofid = np.arange(len(self.vid))
self.a2dofs = self.vid2dofid[self.m.t2v()]
[docs]
def get_tris_mapping(self):
""" return the (linear mapping (see mapping.MappingT3)) of the mesh entities """
if self._tris_map is None:
self._tris_map = mapping.MappingT3(self.m.t2v_coord())
return self._tris_map
[docs]
def scalar_size(self):
""" size (nbdofs) for a scalar space """
return len(self.vid)
[docs]
def size(self):
""" size (nbdofs) taking into account the vector dim of the field """
return self.scalar_size()*self.vecdim
[docs]
def vertexid2dofid(self, vids):
""" return the dofs number associated to the vertices id """
vid_space = np.array([self.vid2dofid[vid] for vid in vids])
return ((self.vecdim*vid_space)[:, np.newaxis] + np.arange(self.vecdim)).squeeze()
[docs]
def edge2edgeid(self, e2v, getxy=False):#pylint: disable=W0613
""" from edge defined by there vertices id, recover the edgeid as seen by the space """
return np.array([], dtype=np.int64).reshape((0, self.vecdim)).squeeze()
[docs]
def edgeid2dofid(self, edges):#pylint: disable=W0613
""" from edge id recover the associated dofid
For FEMSpaceP1, return an empty array
"""
return np.array([], dtype=np.int64).reshape((0, self.vecdim)).squeeze()
[docs]
def interpolate(self, f):
""" interpolate the function f on to mesh nodes """
return f(self.m.vcoord()[self.vid]).flatten()
# if len(fv.flatten()) != len(self.xy)*self.vecDim : raise
# return fv[self.vid].flatten()
[docs]
def operator_dof2val_tri(self, evaltris_uv,
rmap=None, frmap=None, only_on=None, raw_data=False):
""" Return operator N that map the dof value to the field values at each point uv
for each triangle """
return operator_dof2value(FEMSpaceP1, self.vecdim, self.scalar_size(), self.a2dofs,
evaltris_uv, rmap=rmap, frmap=frmap,
only_on=only_on, raw_data=raw_data)
[docs]
def operator_dof2val_edge(self, eids, evaledge_s=quad.Edge_gauss(0).s):
""" Return operator N that map the dof value to the field values at each point s of each
edge """
if not self.allnodes:
raise ValueError("member only coded when all nodes are connected to at \
least one triangle")
N = np.tile( #pylint: disable=C0103
np.stack([1.-evaledge_s, evaledge_s], axis=evaledge_s.ndim).flatten(),
eids.size).reshape((eids.size, evaledge_s.shape[0], 2))
return op2sparse(N, self.scalar_size(), self.m.e2v(eids),
self.vecdim, 2, evaledge_s.shape[0], 1)
[docs]
def operator_dof2grad_tri(self, evaltrisuv,
rmap=None, frmap=None, only_on=None, raw_data=False,
limited=False): #pylint: disable=W0613
""" Return operator dNdx that map the dof value to the gradient of field values at
each point uv for each triangle """
return operator_dof2grad_tri(FEMSpaceP1, self.vecdim, self.m.geo_dim(),
self.scalar_size(), self.a2dofs, self.get_tris_mapping(),
evaltrisuv,
rmap=rmap, frmap=frmap, only_on=only_on, raw_data=raw_data)
[docs]
class FEMSpaceP1NC():
''' Implementation of cruzeix-raviart Element '''
interpolatoryAt = 'midEdge'
dNdu = np.array([[0., 2., -2.], [-2., 2., 0.]])
lenN = 3
[docs]
@staticmethod
def N(uv): #pylint: disable=C0103
""" Shape function for the triangle """
u = uv[..., 0]
v = uv[..., 1]
return np.stack((1. - 2.*v, 2.*(u + v) - 1, 1 - 2.*u), axis=u.ndim)
def __init__(self, name, mesh:SMesh, vecdim=1, tris_map=None) -> None:
self.name = name
self.m = mesh
self.vecdim = vecdim
self.geomdim = mesh.geo_dim()
self._tris_map = tris_map
self.a2dofs = mesh.t2e()
[docs]
def get_tris_mapping(self):
""" return the (linear mapping (see mapping.MappingT3)) of the mesh entities """
if self._tris_map is None:
self._tris_map = mapping.MappingT3(self.m.t2v_coord())
return self._tris_map
[docs]
def scalar_size(self):
""" size (nbdofs) for a scalar space """
return self.m.nedges()
[docs]
def size(self):
""" size (nbdofs) taking into account the vector dim of the field """
return self.scalar_size()*self.vecdim
[docs]
def vertexid2dofid(self, vids): #pylint: disable=W0613
""" return the dofs number associated to the vertices id """
return np.array([], dtype=np.int64).reshape((0, self.vecdim)).squeeze()
[docs]
def edge2edgeid(self, e2v, getxy=False):
""" from edge defined by there vertices id, recover the edgeid as seen by the space """
eid, __ = self.m.find_edgeids(e2v)
if not getxy:
return eid
return eid, self.m.e2v_coord()[eid]
[docs]
def edgeid2dofid(self, eidspace):
""" from edge id recover the associated dofid """
return ((self.vecdim*eidspace)[:, np.newaxis] + np.arange(self.vecdim))
[docs]
def edge2dofid(self, e2v):
""" from edge defined by there vertices id, return the associated dofs """
return self.edgeid2dofid(self.edge2edgeid(e2v))
[docs]
def interpolate(self, f)->npt.NDArray[np.float64]:
""" Compute values of u dofs so that Ni*ui is equal to f at mid edge """
fe = f(np.sum(self.m.e2v_coord(), axis=1)/2.)
if len(fe.flatten()) != self.m.nedges()*self.vecdim:
raise ValueError("size mismatch between f and number of edge dofs")
return fe.flatten()
[docs]
def operator_dof2val_tri(self, evaltrisuv, rmap=None, frmap=None, only_on=None, raw_data=False):
""" Return operator N that map the dof value to the field values at each point uv
for each triangle """
return operator_dof2value(FEMSpaceP1NC, self.vecdim, self.scalar_size(), self.a2dofs,
evaltrisuv, rmap=rmap, frmap=frmap, only_on=only_on, raw_data=raw_data)
[docs]
def operator_dof2grad_tri(self, evaltrisuv, rmap=None, frmap=None, only_on=None, raw_data=False,
limited=False): #pylint: disable=W0613
""" Return operator dNdx that map the dof value to the gradient of field values at
each point uv for each triangle """
return operator_dof2grad_tri(FEMSpaceP1NC, self.vecdim, self.geomdim, self.scalar_size(),
self.a2dofs, self.get_tris_mapping(), evaltrisuv,
rmap=rmap, frmap=frmap, only_on=only_on, raw_data=raw_data)
[docs]
def operator_trival2edge(self, ndata_per_edge=None, data_size=1):
""" return linear operator that push values at triangles to each edge. useful to compute
interelemnt jump for example """
return operator_trival2edge(self.m.nedges(), *self.m.t2e_flip(),
ndataperedge=ndata_per_edge, datasize=data_size)
[docs]
def operator_edgetrace(self, quad_edges):
""" return a linear sparse operator N that compute the trace on each edge at points
quad_edges from the dof values.
t = N@U is an np.array for size nedge*npt*vecdim*2 (it that order)
the entries i,j k of t.reshape((nedge, npt, vecdim, 2))
contain the 2 value at edge i, point j, direction k, which correspond to
the values at left and right of the edge (we assume here a 2d manifold)
"""
npt_per_edge = len(quad_edges.s)
push2edge = self.operator_trival2edge(ndata_per_edge=npt_per_edge, data_size=self.vecdim)
evaltrisuv = quad.T3_edges(quad_edges).uv
N = self.operator_dof2val_tri(evaltrisuv) #pylint: disable=C0103
return push2edge.dot(N)
[docs]
def operator_dof2val_edge(self, eids, evaledges=quad.Edge_gauss(0).s):
""" Return operator N that map the dof value to the field values at each point s of each
edge """
if evaledges != np.array([.5]):
raise ValueError("function only coded to eval at mid edge")
N = np.ones((eids.size, 1, 1)) #pylint: disable=C0103
return op2sparse(N, self.scalar_size(), eids.reshape((-1, 1)), self.vecdim, 1,
evaledges.shape[0], 1)
FEMSpace = Union[FEMSpaceP1, FEMSpaceP1NC]
[docs]
class Heaviside():
""" Heaviside shape function defined on a cut mesh """
def __init__(self, name, mesh_cut_data, ls_tris, ls_tris_cut) -> None:
self.name = name
self.side_appro = (1.-lscut.compute_face_side(ls_tris))*0.5
self.side_eval = (1.-lscut.compute_face_side(ls_tris_cut))*0.5
if isinstance(mesh_cut_data,mapping.RemapHelper):
self.rmap = mesh_cut_data
elif isinstance(mesh_cut_data, MeshMod):
self.rmap = meshcutdata2remap_helper(mesh_cut_data)
else:
raise ValueError("Can't interprete mesh_cut_data")
[docs]
def eval(self, evaltrisuv, only_on_parents=None, frmap=None):
""" evaluate the function at each point uv of each element or sub element"""
if (frmap is not None) and (only_on_parents is not None):
raise ValueError("can't have both frmap and only_parents be None")
if only_on_parents is not None:
frmap = mapping.RemapHelperFiltered(
self.rmap.index_parent2sons, self.rmap.sons, self.rmap.parent_mapping,
self.rmap.son_mapping, only_on_parents)
npts = len(evaltrisuv)
rmap = self.rmap
if frmap is None:
H = np.zeros((rmap.nbevals, npts)) #pylint: disable=C0103
H[rmap.noson2eval_index] = np.repeat(
self.side_appro[rmap.approEQeval, np.newaxis], npts, axis=1)
H[rmap.son2eval_index] = np.repeat(
self.side_eval[rmap.sons, np.newaxis], npts, axis=1)
else: # case where we filter the parent element with frmap.only_on_parent
Htmp = np.zeros((frmap.nbevals, npts)) #pylint: disable=C0103
Htmp[frmap.noson2eval_index] = np.repeat(
self.side_appro[frmap.approEQeval, np.newaxis], npts, axis=1)
Htmp[frmap.son2eval_index] = np.repeat(
self.side_eval[frmap.sons, np.newaxis], npts, axis=1)
H = np.zeros((rmap.nbevals, npts)) #pylint: disable=C0103
i = np.arange(rmap.nbevals)[np.isin(
rmap.eval2appro, frmap.only_on_parents)]
H[i] = Htmp
return H
[docs]
class Ridge():
""" Ridge shape function defined on a cut mesh (usefull for discontinuous gradient) """
def __init__(self, name, mesh_cut_data, ls_tris, ls_tris_cut) -> None:
self.name = name
self.side_appro = (1.-lscut.compute_face_side(ls_tris))*0.5
self.side_eval = (1.-lscut.compute_face_side(ls_tris_cut))*0.5
self.abs_phi_appro_vertices = np.abs(ls_tris)
self.abs_phi_eval_vertices = np.abs(ls_tris_cut)
if isinstance(mesh_cut_data, mapping.RemapHelper):
self.rmap = mesh_cut_data
elif isinstance(mesh_cut_data, MeshMod):
self.rmap = meshcutdata2remap_helper(mesh_cut_data)
else:
raise ValueError("Can't interprete mesh_cut_data")
[docs]
def eval(self, evaltrisuv, only_on_parents=None, frmap=None, raw_data=False):
""" evaluate the function at each point uv of each element or sub element"""
#pylint: disable=C0103
if only_on_parents is not None:
frmap = mapping.RemapHelperFiltered(
self.rmap.index_parent2sons, self.rmap.sons, self.rmap.parent_mapping,
self.rmap.son_mapping, only_on_parents)
npts = len(evaltrisuv)
rmap = self.rmap
N = FEMSpaceP1.N
if frmap is None:
R = np.zeros((rmap.nbevals, npts))
sons_point_on_parents_uv = rmap.remap_on_parent(evaltrisuv)
Nabsphi_appro = np.einsum('lik, lk -> li', N(sons_point_on_parents_uv),
self.abs_phi_appro_vertices[rmap.son2parent_elem_ids])
Nabsphi_eval = np.einsum('...ik, ...jk -> ...ji', N(evaltrisuv),
self.abs_phi_eval_vertices[rmap.sons])
R[rmap.son2eval_index] = Nabsphi_appro - Nabsphi_eval
else:
Rtmp = np.zeros((frmap.nbevals, npts))
sons_point_on_parents_uv = frmap.remap_on_parent(evaltrisuv)
Nabsphi_appro = np.einsum('lik, lk -> li', N(sons_point_on_parents_uv),
self.abs_phi_appro_vertices[frmap.son2parent_elem_ids])
Nabsphi_eval = np.einsum('...ik, ...jk -> ...ji', N(evaltrisuv),
self.abs_phi_eval_vertices[frmap.sons])
Rtmp[frmap.son2eval_index] = Nabsphi_appro - Nabsphi_eval
if raw_data:
return Rtmp
R = np.zeros((rmap.nbevals, npts))
index = np.arange(rmap.nbevals)[np.isin(
rmap.eval2appro, frmap.only_on_parents)]
R[index] = Rtmp
return R
[docs]
def evalGrad(self, evaltrisuv, only_on_parents=None, frmap=None):
""" evaluate the gradiant of thefunction at each point uv of each element or sub element"""
#pylint: disable=C0103
if only_on_parents is not None:
frmap = mapping.RemapHelperFiltered(
self.rmap.index_parent2sons, self.rmap.sons, self.rmap.parent_mapping,
self.rmap.son_mapping, only_on_parents)
npts = len(evaltrisuv)
rmap = self.rmap
inv_F_parent = rmap.parent_mapping.inv_F
inv_F_son = rmap.son_mapping.inv_F
dNdu = FEMSpaceP1.dNdu
dim = 2
if frmap is None:
dRdx = np.zeros((rmap.nbevals, npts, dim))
dNdxparent = gradNTri(dNdu, evaltrisuv, inv_F_parent[rmap.son2parent_elem_ids])
dNdxson = gradNTri(dNdu, evaltrisuv, inv_F_son[rmap.sons])
dNdx_absphi_appro = np.einsum('lijk, lk -> lij', dNdxparent,
self.abs_phi_appro_vertices[rmap.son2parent_elem_ids])
dNdx_absphi_eval = np.einsum('lijk, lk -> lij', dNdxson,
self.abs_phi_eval_vertices[rmap.sons])
dRdx[rmap.son2eval_index] = dNdx_absphi_appro - dNdx_absphi_eval
else:
dRdxtmp = np.zeros((frmap.nbevals, npts, dim))
dNdxparent = gradNTri(dNdu, evaltrisuv, inv_F_parent[frmap.son2parent_elem_ids])
dNdxson = gradNTri(dNdu, evaltrisuv, inv_F_son[frmap.sons])
dNdx_absphi_appro = np.einsum('lijk, lk -> lij', dNdxparent,
self.abs_phi_appro_vertices[frmap.son2parent_elem_ids])
dNdx_absphi_eval = np.einsum('lijk, lk -> lij', dNdxson,
self.abs_phi_eval_vertices[frmap.sons])
dRdxtmp[frmap.son2eval_index] = dNdx_absphi_appro - dNdx_absphi_eval
dRdx = np.zeros((rmap.nbevals, npts, dim))
i = np.arange(rmap.nbevals)[np.isin(
rmap.eval2appro, frmap.only_on_parents)]
dRdx[i] = dRdxtmp
return dRdx
[docs]
class XFEMSpace():
""" XFEM Function space : input an enrichment function and a base space
and a level set on which the enrichment is defined
"""
def __init__(self, name, basespace, enrichment, mesh_cut, ls, ls_cut, rmap=None) -> None:
self.name = name
self.rmap = meshcutdata2remap_helper(mesh_cut) if rmap is None else rmap
self.enrich = enrichment(
name+'_enrich', self.rmap, ls[mesh_cut.tris], ls_cut[mesh_cut.new_tris])
self.base = basespace(name+'_base', SMesh(mesh_cut.xy, mesh_cut.tris),
tris_map=self.rmap.parent_mapping)
cut_elem = self.rmap.approNEQeval
if isinstance(self.enrich, Heaviside):
venr1 = mesh_cut.tris[cut_elem].flatten()
venr0 = np.array(np.where(ls == 0.)[0], dtype=venr1.dtype)
self.enriched_vertices = np.unique(np.hstack((venr0, venr1)))
self.enriched_support = np.where(
np.any(np.isin(mesh_cut.tris, self.enriched_vertices), axis=1))[0]
if isinstance(self.enrich, Ridge):
self.enriched_support = cut_elem
self.enriched_vertices = np.unique(mesh_cut.tris[cut_elem])
self.frmap = mapping.RemapHelperFiltered(self.rmap.index_parent2sons, self.rmap.sons,
self.rmap.parent_mapping, self.rmap.son_mapping, self.enriched_support)
[docs]
def scalar_size(self):
""" size (nbdofs) for a scalar space """
return len(self.enriched_vertices)
[docs]
def operator_dof2val_tri(self, evaltrisuv, rmap = None): #pylint: disable=W0613
""" Return operator N that map the dof value to the field values at each point uv
for each triangle """
# if (rmap != self.rmap):
# raise ValueError("Trying to eval XRidgeSpaceP1 with an rmap not \
# conforming to the constructor's rmap")
base_val = self.base.operator_dof2val_tri(evaltrisuv, self.rmap, self.frmap)
enr_val = (self.enrich.eval(evaltrisuv, frmap=self.frmap).flatten()).reshape((-1, 1))
return (enr_val*base_val[:, self.enriched_vertices]).tocsr()
[docs]
def operator_dof2grad_tri(self, evaltrisuv, rmap = None): #pylint: disable=W0613
""" Return operator dNdx that map the dof value to the gradient of field values at
each point uv for each triangle """
# if (rmap != self.rmap):
# raise ValueError("Trying to eval XFEMSpace with an rmap not \
# conforming to the constructor's rmap")
dim = 2
base_grad = self.base.operator_dof2grad_tri(evaltrisuv, self.rmap, self.frmap)
enr_val = np.repeat(self.enrich.eval(
evaltrisuv, frmap=self.frmap).flatten(), dim).reshape((-1, 1))
dNdx = (enr_val*base_grad[:, self.enriched_vertices]).tocsr() # pylint: disable=C0103
# grad of heaviside is zero every where....
if not isinstance(self.enrich, Heaviside):
grad_enr = self.enrich.evalGrad(evaltrisuv).flatten()
grad_enr = sp.sparse.csr_array((grad_enr, (np.arange(len(grad_enr)), np.repeat(
np.arange(len(grad_enr)//2), 2))), shape=(len(grad_enr), len(grad_enr)//2))
base_val = self.base.operator_dof2val_tri(evaltrisuv, self.rmap, self.frmap)[
:, self.enriched_vertices]
dNdx = dNdx + grad_enr.dot(base_val) # pylint: disable=C0103
return dNdx
[docs]
class XRidgeSpaceP1():
""" A special enrichment space for Ridge enrichment """
def __init__(self, name, mesh_cut, ls, ls_cut, rmap=None, dosvd=False) -> None:
self.dosvd = dosvd
self.mesh_cut = mesh_cut
self.name = name
self.rmap = meshcutdata2remap_helper(mesh_cut) if rmap is None else rmap
cut_elem = self.rmap.approNEQeval
self.enrich = Ridge(name+'_enrich', self.rmap,
ls[mesh_cut.tris], ls_cut[mesh_cut.new_tris])
self.base = FEMSpaceP1(name+'_base', SMesh(mesh_cut.xy, mesh_cut.tris),
vecdim=1, tris_map=self.rmap.parent_mapping)
self.frmap = mapping.RemapHelperFiltered(self.rmap.index_parent2sons, self.rmap.sons,
self.rmap.parent_mapping, self.rmap.son_mapping, cut_elem)
[docs]
def scalar_size(self):
""" size (nbdofs) for a scalar space """
return len(self.mesh_cut.cut_edge_vertices)
[docs]
def operator_dof2val_tri(self, evaltrisuv, rmap=None, raw=False): #pylint: disable=W0613
""" Return operator N that map the dof value to the field values at each point uv
for each triangle """
# if (rmap != self.rmap):
# raise ValueError("Trying to eval XRidgeSpaceP1 with an rmap not \
# conforming to the constructor's rmap")
nodes_uv = quad.T3_nodes().uv
enrval_atnodes:np.ndarray = self.enrich.eval(nodes_uv, frmap=self.frmap,
raw_data=True) # .reshape((-1,1))
nodes_uv = quad.T3_nodes().uv
npts = len(nodes_uv)
cut_vertices = self.mesh_cut.cut_vertices
cut_edges = self.mesh_cut.cut_edge_vertices
cut_tris1 = self.mesh_cut.new_tris
edge_funtmp = np.zeros((self.frmap.nbevals, npts, len(cut_edges)))
for i, v in enumerate(cut_vertices):
edge_funtmp[self.frmap.son2eval_index, :, i] = np.where(cut_tris1 == v, 1., 0.)[
self.frmap.sons]
Nat_nodes_nzelm = enrval_atnodes[:, :, np.newaxis]*edge_funtmp # pylint: disable=C0103
nenr = Nat_nodes_nzelm.shape[-1]
if self.dosvd:
N = Nat_nodes_nzelm.reshape((-1, len(cut_edges))) # pylint: disable=C0103
svd = np.linalg.svd(N)
print('shape : ', N.shape, 'singularvarlues ridge N:', svd[1])
if np.all(nodes_uv != evaltrisuv):
Nat_pt_nzelm = np.einsum('ps, ise -> ipe', # pylint: disable=C0103
FEMSpaceP1.N(evaltrisuv), Nat_nodes_nzelm)
else:
Nat_pt_nzelm = Nat_nodes_nzelm # pylint: disable=C0103
if raw is True:
N = Nat_pt_nzelm # pylint: disable=C0103
else:
eval_elements = np.arange(self.rmap.nbevals)[np.isin(
self.rmap.eval2appro, self.frmap.only_on_parents)]
npts = len(evaltrisuv)
row = np.repeat(np.repeat(npts*eval_elements, npts) +
np.tile(np.arange(npts), len(eval_elements)), nenr)
col = np.tile(np.arange(nenr), npts*len(eval_elements))
N = sp.sparse.csr_array( # pylint: disable=C0103
(Nat_pt_nzelm.flatten(), (row, col)), (self.rmap.nbevals*npts, nenr))
return N
[docs]
def operator_dof2grad_tri(self, evaltrisuv, rmap = None): #pylint: disable=W0613
""" Return operator dNdx that map the dof value to the gradient of field values at
each point uv for each triangle """
# if (rmap != self.rmap):
# raise ValueError("Trying to eval XRidgeSpaceP1 with an rmap not \
# conforming to the constructor's rmap")
dim = 2
npts = len(evaltrisuv)
nodes_uv = quad.T3_nodes().uv
nshape = 3
NRidge = self.operator_dof2val_tri(nodes_uv, raw=True) # pylint: disable=C0103
nenr = NRidge.shape[-1]
NRidge = NRidge.reshape((-1, 3, nenr)) # pylint: disable=C0103
dNdx = np.zeros((self.frmap.nbevals, npts, dim, nshape)) # pylint: disable=C0103
dNdx[self.frmap.son2eval_index] = gradNTri(
FEMSpaceP1.dNdu, evaltrisuv, self.frmap.son_mapping.inv_F[self.rmap.sons])
dNRidgedx = np.einsum('ijkl, ilm ->ijkm', dNdx, NRidge) # pylint: disable=C0103
eval_elements = np.arange(self.rmap.nbevals)[np.isin(
self.rmap.eval2appro, self.frmap.only_on_parents)]
npts = len(evaltrisuv)
row = np.repeat(np.repeat(npts*dim*eval_elements, npts*dim) +
np.tile(np.arange(npts*dim), len(eval_elements)), nenr)
col = np.tile(np.arange(nenr), dim*npts*len(eval_elements))
dN = sp.sparse.csr_array( # pylint: disable=C0103
(dNRidgedx.flatten(), (row, col)), (self.rmap.nbevals*npts*dim, nenr))
return dN
[docs]
class SumSpace():
""" make a space as a sum of different space """
def __init__(self, list_of_space) -> None:
self.list_of_space = list_of_space
[docs]
def operator_dof2val_tri(self, evaltrisuv, rmap=None):
""" Return operator N that map the dof value to the field values at each point uv
for each triangle """
N = [space.operator_dof2val_tri(evaltrisuv, rmap=rmap) # pylint: disable=C0103
for space in self.list_of_space]
return sp.sparse.hstack(N)
[docs]
def operator_dof2grad_tri(self, evaltrisuv, rmap=None):
""" Return operator dNdx that map the dof value to the gradient of field values at
each point uv for each triangle """
dNs = [space.operator_dof2grad_tri(evaltrisuv, rmap=rmap) # pylint: disable=C0103
for space in self.list_of_space]
return sp.sparse.hstack(dNs)
[docs]
class SomeScalarField():
""" an example of scalar field for test purpose """
def __init__(self, f=1.) -> None:
self.f = f
[docs]
def fun(self, xy):
""" eval the function at points xy """
return xy[:, 1]*np.cos(self.f*np.pi*xy[:, 0])
[docs]
def grad(self, xy):
""" eval the gradient of function at points xy """
gx = -xy[:, 1]*self.f*np.pi*np.sin(self.f*np.pi*xy[:, 0])
gy = np.cos(self.f*np.pi*xy[:, 0])
return np.stack((gx, gy), 1)
[docs]
class LinearScalarField():
" A linear field in 2d : u(x,y) = a*x + b*y"
def __init__(self, a=3, b=4) -> None:
self.a, self.b = a, b
[docs]
def fun(self, xy):
""" eval the function at points xy """
return self.a*xy[:, 0]+self.b*xy[:, 1]
[docs]
def grad(self, xy):
""" eval the gradient of function at points xy """
gx = self.a*np.ones(len(xy))
gy = self.b*np.ones(len(xy))
return np.stack((gx, gy), 1)
[docs]
class QuadraticScalarField():
" A Quadratic field for example purpose "
def __init__(self, hess = np.eye(2) ) -> None:
self.hess = 0.5*(hess + hess.T)
[docs]
def fun(self, xy):
""" eval the function at points xy """
x = xy[:, 0]
y = xy[:, 1]
h = self.hess
return .5*h[0,0]*x**2 + (h[0,1]+h[1,0])*x*y+ .5*h[1,1]*y*y
[docs]
def grad(self, xy):
""" eval the gradient of function at points xy """
x = xy[:, 0]
y = xy[:, 1]
h = self.hess
gx = h[0,0]*x + (h[0,1]+h[1,0])*x*y
gy = (h[0,1]+h[1,0])*x + h[1,1]*y
return np.stack((gx, gy), 1)
[docs]
class SomeVectorField():
""" an example of vector field for test purpose """
def __init__(self, f) -> None:
self.f = f
[docs]
def fun(self, xy):
""" eval the function at points xy """
f = self.f
x = xy[:, 0]
y = xy[:,1]
return np.column_stack((y*np.cos(f*np.pi*x), *np.cos(f*np.pi*y)))
[docs]
def grad(self, xy):
""" eval the gradient of function at points xy """
f = self.f
gxx = -xy[:, 1]*f*np.pi*np.sin(f*np.pi*xy[:, 0])
gxy = np.cos(f*np.pi*xy[:, 0])
gyx = np.cos(f*np.pi*xy[:, 1])
gyy = -xy[:, 0]*f*np.pi*np.sin(f*np.pi*xy[:, 1])
return np.column_stack((gxx, gxy, gyx, gyy)).reshape((-1, 2, 2))
[docs]
class QuadraticVectorField():
""" an example of Vector field for test purpose """
def __init__(self) -> None:
pass
[docs]
def fun(self, xy):
""" eval the function at points xy """
return 0.1*np.column_stack((xy[:, 0]**2, xy[:, 1]**2))
[docs]
def grad(self, xy):
""" eval the gradient of function at points xy """
gxx = 0.1*2*xy[:, 0]
gxy = np.zeros(len(xy))
gyx = np.zeros(len(xy))
gyy = 0.1*2*xy[:, 1]
return np.column_stack((gxx, gxy, gyx, gyy)).reshape((-1, 2, 2))
[docs]
class LinearVectorField():
""" an example of Vector field for test purpose """
def __init__(self) -> None:
pass
[docs]
def fun(self, xy):
""" eval the function at points xy """
return np.column_stack((xy[:, 0], 2.*xy[:, 1]))
[docs]
def grad(self, xy):
""" eval the gradient of function at points xy """
gxx = np.ones(len(xy))
gxy = np.zeros(len(xy))
gyx = np.zeros(len(xy))
gyy = 2.*np.ones(len(xy))
return np.column_stack((gxx, gxy, gyx, gyy)).reshape((-1, 2, 2))