Coverage for pyxfem/core/function_spaces.py: 76%
572 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3#pylint: disable=C0302
4"""
5Created on Mon Dec 4 23:30:08 2023
7@author: nchevaug
8"""
9#from warnings import warn
10#from __future__ import annotations
11from abc import ABC
12import unittest
13from typing import Union
14import numpy as np
15import numpy.typing as npt
16import scipy as sp
17import pyxfem.core.level_set_mesh_cut as lscut
18import pyxfem.core.mapping as mapping
19from pyxfem.core.mesh_simplex_one_level import SMesh
20from pyxfem.core.mesh_modifications import MeshHierarchicModification as MeshMod
21import pyxfem.core.quadrature as quad
23class DofManagerBase(ABC):
24 ''' Base class to help manage degree of freedoms '''
25 def scalar_size(self):
26 """ size (nbdofs) for a scalar space """
27 raise NotImplementedError()
29 def size(self):
30 """ size (nbdofs) taking into account the vector dim of the field """
31 raise NotImplementedError()
33 def vertexid2dofid(self, vids):
34 """ return the dofs number associated to the vertices id """
35 raise NotImplementedError()
37 def edge2edgeid(self, e2v, getxy=False):
38 """ from edge defined by there vertices id, recover the edgeid as seen by the space """
39 raise NotImplementedError()
41 def edgeid2dofid(self, edges):
42 """ from edge id recover the associated dofid
43 For FEMSpaceP1, return an empty array
44 """
45 raise NotImplementedError()
47class VertexSpaceDofManager():
48 """ Class to compute the degree of freedom associated to the mesh and the dirichlet BC """
49 # pylint: disable=C0103
50 def __init__(self, m:SMesh, vecdim = 3):
51 self.m = m
52 self.vid = self.m.used_v() # np.unique(tris.flatten())
53 self.vecdim = vecdim
54 self.interpolatoryAt = 'vertex'
55 if len(self.vid) != len(self.m.xy):
56 self.allnodes = False
57 self.vid2dofid = dict(zip(self.vid, np.arange(len(self.vid))))
58 self.a2dofs = np.array([[self.vid2dofid[vid] for vid in t]
59 for t in self.m.t2v()])
60 else:
61 self.allnodes = True
62 self.vid2dofid = np.arange(len(self.vid))
63 self.a2dofs = self.vid2dofid[self.m.t2v()]
64 self.nelem = len(self.a2dofs)
65 self.nscaldof = len(self.vid)
67 def vertexid2dofid(self, vids):
68 ''' return an array A of shape (nv, vecdim) containing such as A[i, j] is the de
69 index of degree of freedom associated to node i in vector direction j
70 '''
72 vid_space = np.array([self.vid2dofid[vid] for vid in vids])
73 return ((self.vecdim*vid_space)[:, np.newaxis] + np.arange(self.vecdim)).squeeze()
75 def size(self):
76 """ Return the total size of the space (number of (scalar) dof in the space )"""
77 return len(self.vid)*self.vecdim
80def op2sparse(NScalar, nscalar_dofs, a2scalardofs, vecdim, #pylint: disable=C0103
81 nshapeperelem, nptperelem, ndataperpoint, only_on=None):
82 """ Given a N array (shape function values at integration point for each element
83 and numbering of degree of freedom dat, construct the sparse N Operator (Assembly)
84 """
85 # nd dataperpoint : not a very good name.
86 # my retionnal is this : vecDim represent the number of scalar dof per node.
87 # For a scalar field it would be 1, For a 2d vector field, it would be 2,
88 # for a 3 d vector Field it would be 2. IT is also possible to have a scalar field, say T
89 # and a 3d vector field say, u stored in the same space. the vec dim would then be 4.
90 # (The function is meant to assemble a multid field, but with the same space for each dir)
91 # The nb data perpoint, is the number of value per vec dim per interpolation point.
92 # For the evaluation of the field, it would be 1, for it's gradient dNdu,
93 # it would be 2 on aface, 1 on an edge and 3 on a 3d element.
94 nelem = len(a2scalardofs)
95 neval = len(NScalar)
96 npts = nelem*nptperelem
97 nrow = npts*ndataperpoint*vecdim
98 ncol = nscalar_dofs*vecdim
99 shape = (nrow, ncol)
100 if only_on is None:
101 col = np.repeat(a2scalardofs, nptperelem*ndataperpoint, axis=0)
102 row_ptr = np.arange(0, (nrow+1)*nshapeperelem, nshapeperelem)
103 else:
104 col = np.repeat(a2scalardofs[only_on], nptperelem*ndataperpoint, axis=0)
105 row = np.repeat(np.arange(nrow).reshape((nelem, -1))[only_on].flatten(), nshapeperelem)
106 row_ptr = np.cumsum(np.bincount(row+1, minlength=nrow+1))
107 if vecdim > 1:
108 col = col.reshape((-1, nshapeperelem*ndataperpoint))
109 col = np.column_stack([vecdim*col + i for i in range(vecdim)]).flatten()
110 rshape = (neval, nptperelem, ndataperpoint, nshapeperelem)
111 N = np.tile(NScalar.reshape(rshape), [1,1,vecdim,1]).flatten() #pylint: disable=C0103
112 else:
113 col = col.flatten()
114 N = NScalar.flatten() #pylint: disable=C0103
115 N = sp.sparse.csr_array((N, col, row_ptr), shape) #pylint: disable=C0103
116 N.sort_indices()
117 return N
119def op2sparse_rmap(NScalar, nscalar_dofs, a2scalardofs, vecdim, #pylint: disable=C0103
120 nshapeperelem, nptperelem,
121 ndataperpoint, rmap, frmap=None, sort=True):
122 """ Given a N array (shape function values at integration point for each element
123 (or sub element)) and numbering of degree of freedom dat, construct
124 the sparse N Operator (Assembly)
125 Specialized version for the case (typically xfem) where element are cut in sub cell
126 """
127 #nelem = len(a2ScalarDofs)
128 #neval = len(NScalar)
129 #npts = nelem*nptperelem
130 #shape = (npts*ndataperpoint*vecDim, nScalarDofs*vecDim)
131 npts = rmap.nbevals*nptperelem
132 nscalrow = npts*ndataperpoint # nb row for a scalar field
133 nrow = nscalrow*vecdim
134 ncol = nscalar_dofs*vecdim
135 shape = (nrow, ncol)
136 if frmap is None:
137 col = np.repeat(a2scalardofs, rmap.appro2nbevals*nptperelem*ndataperpoint, axis=0)
138 row_ptr = np.arange(0, (nrow+1)*nshapeperelem, nshapeperelem)
139 if vecdim > 1:
140 col = col.reshape((-1, nshapeperelem*ndataperpoint))
141 col = np.column_stack([vecdim*col + i for i in range(vecdim)]).flatten()
142 rshape = (rmap.nbevals, nptperelem, ndataperpoint, nshapeperelem)
143 N = np.tile(NScalar.reshape(rshape), [1,1,vecdim,1]).flatten() #pylint: disable=C0103
144 else:
145 col = col.flatten()
146 N = NScalar.flatten() #pylint: disable=C0103
147 else:
148 fscalcol = a2scalardofs[frmap.only_on_parents]
149 col = np.repeat(fscalcol, frmap.appro2nbevals*nptperelem*ndataperpoint, axis=0)
150 row = np.arange(rmap.nbevals*nptperelem*ndataperpoint).reshape((rmap.nbevals, -1))
151 row = row[np.isin(rmap.eval2appro, frmap.only_on_parents)].flatten()
152 if vecdim > 1:
153 row = np.column_stack(
154 tuple(vecdim*row + i for i in range(vecdim))).flatten()
155 col = col.reshape((-1, nshapeperelem*ndataperpoint))
156 col = np.column_stack([vecdim*col + i for i in range(vecdim)]).flatten()
157 rshape = (frmap.nbevals, nptperelem, ndataperpoint, nshapeperelem)
158 N = np.tile(NScalar.reshape(rshape), [1,1,vecdim,1]).flatten() #pylint: disable=C0103
159 #warn("op2sparse_rmap Not fully tested for case vecDim > 1 and frmap != None")
160 else:
161 row = row.flatten()
162 col = col.flatten()
163 N = NScalar.flatten() #pylint: disable=C0103
164 row_ptr = np.cumsum(nshapeperelem*np.bincount(row+1, minlength=nrow+1))
165 N = sp.sparse.csr_array((N, col, row_ptr), shape) #pylint: disable=C0103
166 if sort:
167 N.sort_indices()
168 return N
170def meshcutdata2remap_helper(mesh_cut, hmin=0., alpha=0.1, p=1):
171 """ build a RemapHelper Object from mesh_cut data """
172 appro_mapping = mapping.MappingT3(
173 mesh_cut.xy[mesh_cut.tris], hmin=hmin, alpha=alpha, p=p)
174 appro2evals_index, appro2evals = mesh_cut.parents2sons()
175 eval_tris_xy = mesh_cut.xy_new_vertices[mesh_cut.new_tris]
176 eval_mapping = mapping.MappingT3(eval_tris_xy, hmin=hmin, alpha=alpha, p=p)
177 return mapping.RemapHelper(appro2evals_index, appro2evals, appro_mapping, eval_mapping)
179def operator_trival2edge(nedge, t2e, t2eflip, ndataperedge=None, datasize=1, dtype=np.float64):
180 """ return linear operator that push values at triangles to each edge. useful to compute
181 interelemnt jump for example """
182 ntri = len(t2e)
183 if ndataperedge is None:
184 push2edge_data = np.ones(ntri*datasize*3, dtype=dtype)
185 push2edge_row = np.column_stack( [(2*datasize*t2e.flatten()+i \
186 + (datasize)*t2eflip.flatten()) for i in range(datasize)]).flatten()
187 push2edge_col = np.repeat(np.arange(ntri*datasize).reshape((ntri,-1)), 3, axis = 0) \
188 .flatten()
189 shape = (2*nedge*datasize, ntri*datasize)
190 else:
191 push2edge_data = np.ones(ntri*3*ndataperedge*datasize, dtype=dtype)
192 push2edge_col = np.arange(ntri*3*ndataperedge*datasize)
193 noflip = 2*np.arange(ndataperedge *
194 datasize).reshape((ndataperedge, datasize))
195 doflip = noflip[::-1]+1
196 push2edge_row = (ndataperedge*datasize*2*t2e[:, :, np.newaxis, np.newaxis] +
197 np.where((t2eflip == 0)[:, :, np.newaxis, np.newaxis],
198 noflip[np.newaxis, np.newaxis, :, :],
199 doflip[np.newaxis, np.newaxis, :, :])
200 ).flatten()
201 shape = (nedge*ndataperedge*datasize*2, ntri*ndataperedge*3*datasize)
202 return sp.sparse.csr_array((push2edge_data, (push2edge_row, push2edge_col)), shape=shape)
204def operator_dof2value(space, vecdim, scalar_space_size, tri2dofs, evaltrisuv,
205 rmap=None, frmap=None, raw_data=False, only_on = None):
206 """ return a linear operator N that map dof value to value at interpolation point of each
207 element """
208 npts = len(evaltrisuv)
209 nshape = space.lenN
210 Neval = space.N(evaltrisuv) #pylint: disable=C0103
211 N = None #pylint: disable=C0103
212 if rmap is None and frmap is None:
213 if only_on is None:
214 N = np.broadcast_to(Neval, (len(tri2dofs), npts, nshape)) #pylint: disable=C0103
215 else:
216 N = np.broadcast_to(Neval, (len(only_on), npts, nshape)) #pylint: disable=C0103
217 if not raw_data:
218 N = op2sparse(N, scalar_space_size, tri2dofs, vecdim, #pylint: disable=C0103
219 nshape, npts, 1, only_on=only_on)
220 if rmap is not None:
221 rmap1 = rmap if frmap is None else frmap
222 sons_point_on_parents_uv = rmap1.remap_on_parent(evaltrisuv)
223 N = np.zeros((rmap1.nbevals, npts, nshape)) #pylint: disable=C0103
224 N[rmap1.noson2eval_index] = Neval
225 N[rmap1.son2eval_index] = space.N(sons_point_on_parents_uv)
226 if not raw_data:
227 N = op2sparse_rmap( #pylint: disable=C0103
228 N, scalar_space_size, tri2dofs, vecdim, nshape, npts, 1, rmap, frmap)
229 if N is None:
230 raise ValueError("N was not constructed due to wrong input")
231 return N
233def gradNTri(dNdu, uv, invF): #pylint: disable=C0103
234 ''' Compute the gradient with regard to space coordinates xy from gradient with regard to
235 element coordinate uv'''
236 dNdx = np.einsum('...ij,...ik->...jk', invF, dNdu) #pylint: disable=C0103
237 # map is constant and DNdU is constant.
238 # repeat the value as much as the number of points per elem
239 return np.repeat(dNdx[:, np.newaxis, :, :], uv.shape[-2], axis=1)
242def operator_dof2grad_tri(space, vecdim, geomdim, scalar_space_size, tri2dofs, trimap, evaltrisuv,
243 rmap=None, frmap=None, only_on=None, raw_data=False, limited=False):
244 """ helper function : return operator N that map the dof value to the gradient of field values
245 at each point uv for each triangle """
246 npts = len(evaltrisuv)
247 nshape = space.lenN
248 dNdx = None #pylint: disable=C0103
249 invF = trimap.inv_F if not limited else trimap.invFLimited #pylint: disable=C0103
250 points_uv = np.zeros((0, npts, 0))
251 if rmap is None and frmap is None:
252 invF = invF if only_on is None else invF[only_on] #pylint: disable=C0103
253 dNdx = gradNTri(space.dNdu, points_uv, invF) #pylint: disable=C0103
254 if not raw_data:
255 dNdx = op2sparse(dNdx, scalar_space_size, tri2dofs, #pylint: disable=C0103
256 vecdim, nshape, npts, geomdim, only_on=only_on)
257 elif vecdim!=1:
258 neval = len(dNdx)
259 nptperelem = npts
260 ndataperpoint = geomdim
261 nshapeperelem = nshape
262 rshape = (neval, nptperelem, ndataperpoint, nshapeperelem)
263 dNdx.reshape(rshape)
264 rshapev = (neval, nptperelem, ndataperpoint*vecdim, nshapeperelem*vecdim)
265 dNdxv = np.zeros(rshapev) #pylint: disable=C0103
266 for r in range(vecdim):
267 rs = ndataperpoint*r
268 re = ndataperpoint*(r+1)
269 rows = np.arange(rs,re,1)
270 cs = r
271 ce = nshapeperelem*vecdim
272 cols = np.arange(cs, ce, vecdim)
273 dNdxv[..., rows[:, np.newaxis], cols] = dNdx
274 dNdx = dNdxv #pylint: disable=C0103
275 if rmap is not None:
276 # Warning : all coordinate are put to zeros because here we don't care :
277 # dNdu is constant on the element
278 son_on_parent_uv = np.zeros((0, npts, 0))
279 parent_uv = np.zeros((0, npts, 0))
280 rmap1 = rmap if frmap is None else frmap
281 dNdu = space.dNdu #pylint: disable=C0103
282 dNdx = np.zeros((rmap1.nbevals, npts, geomdim, nshape)) #pylint: disable=C0103
283 dNdx[rmap1.noson2eval_index] = gradNTri(dNdu, parent_uv, invF[rmap1.approEQeval])
284 dNdx[rmap1.son2eval_index] = gradNTri(dNdu, son_on_parent_uv,
285 invF[rmap1.son2parent_elem_ids])
286 if not raw_data:
287 dNdx = op2sparse_rmap( #pylint: disable=C0103
288 dNdx, scalar_space_size, tri2dofs, vecdim, nshape, npts, geomdim, rmap, frmap)
289 if dNdx is None:
290 raise ValueError("dNdx was not constructed due to wrong input")
291 return dNdx
293class TestSparseOp(unittest.TestCase):
294 """ Class for unit test of operator on dofs """
295 def test_operator_trival2edge(self):
296 """ Test of the operator tha push triangles data to edges """
297 xy = np.array([[0.,0.], [1.,0.], [0.,1.], [1.,1.]])
298 tris = np.array([[0,1,2], [1,3,2]])
299 mesh = SMesh(xy, tris)
300 t2e, t2ef = mesh.t2e_flip()
301 nedges = mesh.nedges()
302 ntris = mesh.ntris()
303 # Testing operator_trival2edges, for datasize = 1
304 opt2e = operator_trival2edge(nedges, t2e, t2ef)
305 valtri = np.arange(1, ntris+1)
306 valedge = (opt2e@valtri).reshape((-1,2))
307 valedgeref = np.array([1., 0.,0.,1.,1.,2.,2.,0., 0.,2.]).reshape((-1,2))
308 err = float(np.linalg.norm(valedge - valedgeref))
309 self.assertAlmostEqual(err, 0.,\
310 msg = "test_operator_trival2edge do not return the expected result for case datasize=1")
311 # Testing operator_trival2edges, for datasize = 3
312 datasize = (3,)
313 opt2e = operator_trival2edge(nedges, t2e, t2ef, datasize= datasize[0])
314 valtri = np.arange(1, ntris*np.prod(datasize)+1).reshape((-1,) + datasize)
315 valedge = (opt2e@valtri.flatten()).reshape((-1,2) + datasize)
316 valedgeref = np.array([1.,2.,3.,0.,0.,0., 0.,0.,0.,1.,2.,3., 1.,2.,3.,4.,5.,6., \
317 4.,5.,6., 0.,0.,0., 0.,0.,0., 4.,5.,6.]).reshape((-1,2) + datasize)
318 err = float(np.linalg.norm(valedge - valedgeref))
319 self.assertAlmostEqual(err, 0.,\
320 msg = "test_operator_trival2edge do not return the expected result for case datasize=3")
323class FEMSpaceP1():
324 """ Classical Linear Finite Element space """
325 interpolatoryAt = 'vertex'
326 dNdu = np.array([[-1., 1., 0.], [-1., 0., 1.]])
327 lenN = 3
329 @staticmethod
330 def N(uv): #pylint: disable=C0103
331 """ Shape function for the triangle """
332 u = uv[..., 0]
333 v = uv[..., 1]
334 return np.stack((1.-u-v, u, v), axis=u.ndim)
336 def __init__(self, name:str, mesh:SMesh,
337 vecdim:int=1, tris_map:None|mapping.MappingT3 = None) -> None:
338 self.name = name
339 self.m = mesh
340 self.vecdim = vecdim
341 self._tris_map = tris_map
342 self.vid = self.m.used_v() # np.unique(tris.flatten())
343 print(len(self.vid) , len(self.m.vcoord()))
344 if len(self.vid) != len(self.m.vcoord()):
345 self.allnodes = False
346 self.vid2dofid = dict(zip(self.vid, np.arange(len(self.vid))))
347 self.a2dofs = np.array([[self.vid2dofid[vid] for vid in t]
348 for t in self.m.t2v()])
349 else:
350 self.allnodes = True
351 self.vid2dofid = np.arange(len(self.vid))
352 self.a2dofs = self.vid2dofid[self.m.t2v()]
353 def get_tris_mapping(self):
354 """ return the (linear mapping (see mapping.MappingT3)) of the mesh entities """
355 if self._tris_map is None:
356 self._tris_map = mapping.MappingT3(self.m.t2v_coord())
357 return self._tris_map
359 def scalar_size(self):
360 """ size (nbdofs) for a scalar space """
361 return len(self.vid)
363 def size(self):
364 """ size (nbdofs) taking into account the vector dim of the field """
365 return self.scalar_size()*self.vecdim
367 def vertexid2dofid(self, vids):
368 """ return the dofs number associated to the vertices id """
369 vid_space = np.array([self.vid2dofid[vid] for vid in vids])
370 return ((self.vecdim*vid_space)[:, np.newaxis] + np.arange(self.vecdim)).squeeze()
372 def edge2edgeid(self, e2v, getxy=False):#pylint: disable=W0613
373 """ from edge defined by there vertices id, recover the edgeid as seen by the space """
374 return np.array([], dtype=np.int64).reshape((0, self.vecdim)).squeeze()
376 def edgeid2dofid(self, edges):#pylint: disable=W0613
377 """ from edge id recover the associated dofid
378 For FEMSpaceP1, return an empty array
379 """
380 return np.array([], dtype=np.int64).reshape((0, self.vecdim)).squeeze()
382 def interpolate(self, f):
383 """ interpolate the function f on to mesh nodes """
384 return f(self.m.vcoord()[self.vid]).flatten()
385 # if len(fv.flatten()) != len(self.xy)*self.vecDim : raise
386 # return fv[self.vid].flatten()
388 def operator_dof2val_tri(self, evaltris_uv,
389 rmap=None, frmap=None, only_on=None, raw_data=False):
390 """ Return operator N that map the dof value to the field values at each point uv
391 for each triangle """
392 return operator_dof2value(FEMSpaceP1, self.vecdim, self.scalar_size(), self.a2dofs,
393 evaltris_uv, rmap=rmap, frmap=frmap,
394 only_on=only_on, raw_data=raw_data)
396 def operator_dof2val_edge(self, eids, evaledge_s=quad.Edge_gauss(0).s):
397 """ Return operator N that map the dof value to the field values at each point s of each
398 edge """
399 if not self.allnodes:
400 raise ValueError("member only coded when all nodes are connected to at \
401 least one triangle")
402 N = np.tile( #pylint: disable=C0103
403 np.stack([1.-evaledge_s, evaledge_s], axis=evaledge_s.ndim).flatten(),
404 eids.size).reshape((eids.size, evaledge_s.shape[0], 2))
405 return op2sparse(N, self.scalar_size(), self.m.e2v(eids),
406 self.vecdim, 2, evaledge_s.shape[0], 1)
408 def operator_dof2grad_tri(self, evaltrisuv,
409 rmap=None, frmap=None, only_on=None, raw_data=False,
410 limited=False): #pylint: disable=W0613
411 """ Return operator dNdx that map the dof value to the gradient of field values at
412 each point uv for each triangle """
413 return operator_dof2grad_tri(FEMSpaceP1, self.vecdim, self.m.geo_dim(),
414 self.scalar_size(), self.a2dofs, self.get_tris_mapping(),
415 evaltrisuv,
416 rmap=rmap, frmap=frmap, only_on=only_on, raw_data=raw_data)
418class FEMSpaceP1NC():
419 ''' Implementation of cruzeix-raviart Element '''
420 interpolatoryAt = 'midEdge'
421 dNdu = np.array([[0., 2., -2.], [-2., 2., 0.]])
422 lenN = 3
424 @staticmethod
425 def N(uv): #pylint: disable=C0103
426 """ Shape function for the triangle """
427 u = uv[..., 0]
428 v = uv[..., 1]
429 return np.stack((1. - 2.*v, 2.*(u + v) - 1, 1 - 2.*u), axis=u.ndim)
431 def __init__(self, name, mesh:SMesh, vecdim=1, tris_map=None) -> None:
432 self.name = name
433 self.m = mesh
434 self.vecdim = vecdim
435 self.geomdim = mesh.geo_dim()
436 self._tris_map = tris_map
437 self.a2dofs = mesh.t2e()
439 def get_tris_mapping(self):
440 """ return the (linear mapping (see mapping.MappingT3)) of the mesh entities """
441 if self._tris_map is None:
442 self._tris_map = mapping.MappingT3(self.m.t2v_coord())
443 return self._tris_map
445 def scalar_size(self):
446 """ size (nbdofs) for a scalar space """
447 return self.m.nedges()
449 def size(self):
450 """ size (nbdofs) taking into account the vector dim of the field """
451 return self.scalar_size()*self.vecdim
453 def vertexid2dofid(self, vids): #pylint: disable=W0613
454 """ return the dofs number associated to the vertices id """
455 return np.array([], dtype=np.int64).reshape((0, self.vecdim)).squeeze()
457 def edge2edgeid(self, e2v, getxy=False):
458 """ from edge defined by there vertices id, recover the edgeid as seen by the space """
459 eid, __ = self.m.find_edgeids(e2v)
460 if not getxy:
461 return eid
462 return eid, self.m.e2v_coord()[eid]
464 def edgeid2dofid(self, eidspace):
465 """ from edge id recover the associated dofid """
466 return ((self.vecdim*eidspace)[:, np.newaxis] + np.arange(self.vecdim))
468 def edge2dofid(self, e2v):
469 """ from edge defined by there vertices id, return the associated dofs """
470 return self.edgeid2dofid(self.edge2edgeid(e2v))
472 def interpolate(self, f)->npt.NDArray[np.float64]:
473 """ Compute values of u dofs so that Ni*ui is equal to f at mid edge """
474 fe = f(np.sum(self.m.e2v_coord(), axis=1)/2.)
475 if len(fe.flatten()) != self.m.nedges()*self.vecdim:
476 raise ValueError("size mismatch between f and number of edge dofs")
477 return fe.flatten()
479 def operator_dof2val_tri(self, evaltrisuv, rmap=None, frmap=None, only_on=None, raw_data=False):
480 """ Return operator N that map the dof value to the field values at each point uv
481 for each triangle """
482 return operator_dof2value(FEMSpaceP1NC, self.vecdim, self.scalar_size(), self.a2dofs,
483 evaltrisuv, rmap=rmap, frmap=frmap, only_on=only_on, raw_data=raw_data)
485 def operator_dof2grad_tri(self, evaltrisuv, rmap=None, frmap=None, only_on=None, raw_data=False,
486 limited=False): #pylint: disable=W0613
487 """ Return operator dNdx that map the dof value to the gradient of field values at
488 each point uv for each triangle """
489 return operator_dof2grad_tri(FEMSpaceP1NC, self.vecdim, self.geomdim, self.scalar_size(),
490 self.a2dofs, self.get_tris_mapping(), evaltrisuv,
491 rmap=rmap, frmap=frmap, only_on=only_on, raw_data=raw_data)
493 def operator_trival2edge(self, ndata_per_edge=None, data_size=1):
494 """ return linear operator that push values at triangles to each edge. useful to compute
495 interelemnt jump for example """
496 return operator_trival2edge(self.m.nedges(), *self.m.t2e_flip(),
497 ndataperedge=ndata_per_edge, datasize=data_size)
499 def operator_edgetrace(self, quad_edges):
500 """ return a linear sparse operator N that compute the trace on each edge at points
501 quad_edges from the dof values.
502 t = N@U is an np.array for size nedge*npt*vecdim*2 (it that order)
503 the entries i,j k of t.reshape((nedge, npt, vecdim, 2))
504 contain the 2 value at edge i, point j, direction k, which correspond to
505 the values at left and right of the edge (we assume here a 2d manifold)
506 """
507 npt_per_edge = len(quad_edges.s)
508 push2edge = self.operator_trival2edge(ndata_per_edge=npt_per_edge, data_size=self.vecdim)
509 evaltrisuv = quad.T3_edges(quad_edges).uv
510 N = self.operator_dof2val_tri(evaltrisuv) #pylint: disable=C0103
511 return push2edge.dot(N)
513 def operator_dof2val_edge(self, eids, evaledges=quad.Edge_gauss(0).s):
514 """ Return operator N that map the dof value to the field values at each point s of each
515 edge """
516 if evaledges != np.array([.5]):
517 raise ValueError("function only coded to eval at mid edge")
518 N = np.ones((eids.size, 1, 1)) #pylint: disable=C0103
519 return op2sparse(N, self.scalar_size(), eids.reshape((-1, 1)), self.vecdim, 1,
520 evaledges.shape[0], 1)
522FEMSpace = Union[FEMSpaceP1, FEMSpaceP1NC]
524class Heaviside():
525 """ Heaviside shape function defined on a cut mesh """
526 def __init__(self, name, mesh_cut_data, ls_tris, ls_tris_cut) -> None:
527 self.name = name
528 self.side_appro = (1.-lscut.compute_face_side(ls_tris))*0.5
529 self.side_eval = (1.-lscut.compute_face_side(ls_tris_cut))*0.5
530 if isinstance(mesh_cut_data,mapping.RemapHelper):
531 self.rmap = mesh_cut_data
532 elif isinstance(mesh_cut_data, MeshMod):
533 self.rmap = meshcutdata2remap_helper(mesh_cut_data)
534 else:
535 raise ValueError("Can't interprete mesh_cut_data")
537 def eval(self, evaltrisuv, only_on_parents=None, frmap=None):
538 """ evaluate the function at each point uv of each element or sub element"""
539 if (frmap is not None) and (only_on_parents is not None):
540 raise ValueError("can't have both frmap and only_parents be None")
541 if only_on_parents is not None:
542 frmap = mapping.RemapHelperFiltered(
543 self.rmap.index_parent2sons, self.rmap.sons, self.rmap.parent_mapping,
544 self.rmap.son_mapping, only_on_parents)
545 npts = len(evaltrisuv)
546 rmap = self.rmap
547 if frmap is None:
548 H = np.zeros((rmap.nbevals, npts)) #pylint: disable=C0103
549 H[rmap.noson2eval_index] = np.repeat(
550 self.side_appro[rmap.approEQeval, np.newaxis], npts, axis=1)
551 H[rmap.son2eval_index] = np.repeat(
552 self.side_eval[rmap.sons, np.newaxis], npts, axis=1)
553 else: # case where we filter the parent element with frmap.only_on_parent
554 Htmp = np.zeros((frmap.nbevals, npts)) #pylint: disable=C0103
555 Htmp[frmap.noson2eval_index] = np.repeat(
556 self.side_appro[frmap.approEQeval, np.newaxis], npts, axis=1)
557 Htmp[frmap.son2eval_index] = np.repeat(
558 self.side_eval[frmap.sons, np.newaxis], npts, axis=1)
559 H = np.zeros((rmap.nbevals, npts)) #pylint: disable=C0103
560 i = np.arange(rmap.nbevals)[np.isin(
561 rmap.eval2appro, frmap.only_on_parents)]
562 H[i] = Htmp
563 return H
565class Ridge():
566 """ Ridge shape function defined on a cut mesh (usefull for discontinuous gradient) """
567 def __init__(self, name, mesh_cut_data, ls_tris, ls_tris_cut) -> None:
568 self.name = name
569 self.side_appro = (1.-lscut.compute_face_side(ls_tris))*0.5
570 self.side_eval = (1.-lscut.compute_face_side(ls_tris_cut))*0.5
571 self.abs_phi_appro_vertices = np.abs(ls_tris)
572 self.abs_phi_eval_vertices = np.abs(ls_tris_cut)
573 if isinstance(mesh_cut_data, mapping.RemapHelper):
574 self.rmap = mesh_cut_data
575 elif isinstance(mesh_cut_data, MeshMod):
576 self.rmap = meshcutdata2remap_helper(mesh_cut_data)
577 else:
578 raise ValueError("Can't interprete mesh_cut_data")
580 def eval(self, evaltrisuv, only_on_parents=None, frmap=None, raw_data=False):
581 """ evaluate the function at each point uv of each element or sub element"""
582 #pylint: disable=C0103
583 if only_on_parents is not None:
584 frmap = mapping.RemapHelperFiltered(
585 self.rmap.index_parent2sons, self.rmap.sons, self.rmap.parent_mapping,
586 self.rmap.son_mapping, only_on_parents)
587 npts = len(evaltrisuv)
588 rmap = self.rmap
589 N = FEMSpaceP1.N
590 if frmap is None:
591 R = np.zeros((rmap.nbevals, npts))
592 sons_point_on_parents_uv = rmap.remap_on_parent(evaltrisuv)
593 Nabsphi_appro = np.einsum('lik, lk -> li', N(sons_point_on_parents_uv),
594 self.abs_phi_appro_vertices[rmap.son2parent_elem_ids])
595 Nabsphi_eval = np.einsum('...ik, ...jk -> ...ji', N(evaltrisuv),
596 self.abs_phi_eval_vertices[rmap.sons])
597 R[rmap.son2eval_index] = Nabsphi_appro - Nabsphi_eval
599 else:
600 Rtmp = np.zeros((frmap.nbevals, npts))
601 sons_point_on_parents_uv = frmap.remap_on_parent(evaltrisuv)
602 Nabsphi_appro = np.einsum('lik, lk -> li', N(sons_point_on_parents_uv),
603 self.abs_phi_appro_vertices[frmap.son2parent_elem_ids])
604 Nabsphi_eval = np.einsum('...ik, ...jk -> ...ji', N(evaltrisuv),
605 self.abs_phi_eval_vertices[frmap.sons])
606 Rtmp[frmap.son2eval_index] = Nabsphi_appro - Nabsphi_eval
607 if raw_data:
608 return Rtmp
609 R = np.zeros((rmap.nbevals, npts))
610 index = np.arange(rmap.nbevals)[np.isin(
611 rmap.eval2appro, frmap.only_on_parents)]
612 R[index] = Rtmp
613 return R
615 def evalGrad(self, evaltrisuv, only_on_parents=None, frmap=None):
616 """ evaluate the gradiant of thefunction at each point uv of each element or sub element"""
617 #pylint: disable=C0103
618 if only_on_parents is not None:
619 frmap = mapping.RemapHelperFiltered(
620 self.rmap.index_parent2sons, self.rmap.sons, self.rmap.parent_mapping,
621 self.rmap.son_mapping, only_on_parents)
622 npts = len(evaltrisuv)
623 rmap = self.rmap
624 inv_F_parent = rmap.parent_mapping.inv_F
625 inv_F_son = rmap.son_mapping.inv_F
626 dNdu = FEMSpaceP1.dNdu
627 dim = 2
628 if frmap is None:
629 dRdx = np.zeros((rmap.nbevals, npts, dim))
630 dNdxparent = gradNTri(dNdu, evaltrisuv, inv_F_parent[rmap.son2parent_elem_ids])
631 dNdxson = gradNTri(dNdu, evaltrisuv, inv_F_son[rmap.sons])
632 dNdx_absphi_appro = np.einsum('lijk, lk -> lij', dNdxparent,
633 self.abs_phi_appro_vertices[rmap.son2parent_elem_ids])
634 dNdx_absphi_eval = np.einsum('lijk, lk -> lij', dNdxson,
635 self.abs_phi_eval_vertices[rmap.sons])
636 dRdx[rmap.son2eval_index] = dNdx_absphi_appro - dNdx_absphi_eval
637 else:
638 dRdxtmp = np.zeros((frmap.nbevals, npts, dim))
639 dNdxparent = gradNTri(dNdu, evaltrisuv, inv_F_parent[frmap.son2parent_elem_ids])
640 dNdxson = gradNTri(dNdu, evaltrisuv, inv_F_son[frmap.sons])
641 dNdx_absphi_appro = np.einsum('lijk, lk -> lij', dNdxparent,
642 self.abs_phi_appro_vertices[frmap.son2parent_elem_ids])
643 dNdx_absphi_eval = np.einsum('lijk, lk -> lij', dNdxson,
644 self.abs_phi_eval_vertices[frmap.sons])
645 dRdxtmp[frmap.son2eval_index] = dNdx_absphi_appro - dNdx_absphi_eval
646 dRdx = np.zeros((rmap.nbevals, npts, dim))
647 i = np.arange(rmap.nbevals)[np.isin(
648 rmap.eval2appro, frmap.only_on_parents)]
649 dRdx[i] = dRdxtmp
650 return dRdx
653class XFEMSpace():
654 """ XFEM Function space : input an enrichment function and a base space
655 and a level set on which the enrichment is defined
656 """
657 def __init__(self, name, basespace, enrichment, mesh_cut, ls, ls_cut, rmap=None) -> None:
658 self.name = name
659 self.rmap = meshcutdata2remap_helper(mesh_cut) if rmap is None else rmap
660 self.enrich = enrichment(
661 name+'_enrich', self.rmap, ls[mesh_cut.tris], ls_cut[mesh_cut.new_tris])
662 self.base = basespace(name+'_base', SMesh(mesh_cut.xy, mesh_cut.tris),
663 tris_map=self.rmap.parent_mapping)
664 cut_elem = self.rmap.approNEQeval
665 if isinstance(self.enrich, Heaviside):
666 venr1 = mesh_cut.tris[cut_elem].flatten()
667 venr0 = np.array(np.where(ls == 0.)[0], dtype=venr1.dtype)
668 self.enriched_vertices = np.unique(np.hstack((venr0, venr1)))
669 self.enriched_support = np.where(
670 np.any(np.isin(mesh_cut.tris, self.enriched_vertices), axis=1))[0]
671 if isinstance(self.enrich, Ridge):
672 self.enriched_support = cut_elem
673 self.enriched_vertices = np.unique(mesh_cut.tris[cut_elem])
674 self.frmap = mapping.RemapHelperFiltered(self.rmap.index_parent2sons, self.rmap.sons,
675 self.rmap.parent_mapping, self.rmap.son_mapping, self.enriched_support)
677 def scalar_size(self):
678 """ size (nbdofs) for a scalar space """
679 return len(self.enriched_vertices)
681 def operator_dof2val_tri(self, evaltrisuv, rmap = None): #pylint: disable=W0613
682 """ Return operator N that map the dof value to the field values at each point uv
683 for each triangle """
684 # if (rmap != self.rmap):
685 # raise ValueError("Trying to eval XRidgeSpaceP1 with an rmap not \
686 # conforming to the constructor's rmap")
687 base_val = self.base.operator_dof2val_tri(evaltrisuv, self.rmap, self.frmap)
688 enr_val = (self.enrich.eval(evaltrisuv, frmap=self.frmap).flatten()).reshape((-1, 1))
689 return (enr_val*base_val[:, self.enriched_vertices]).tocsr()
691 def operator_dof2grad_tri(self, evaltrisuv, rmap = None): #pylint: disable=W0613
692 """ Return operator dNdx that map the dof value to the gradient of field values at
693 each point uv for each triangle """
694 # if (rmap != self.rmap):
695 # raise ValueError("Trying to eval XFEMSpace with an rmap not \
696 # conforming to the constructor's rmap")
697 dim = 2
698 base_grad = self.base.operator_dof2grad_tri(evaltrisuv, self.rmap, self.frmap)
699 enr_val = np.repeat(self.enrich.eval(
700 evaltrisuv, frmap=self.frmap).flatten(), dim).reshape((-1, 1))
701 dNdx = (enr_val*base_grad[:, self.enriched_vertices]).tocsr() # pylint: disable=C0103
702 # grad of heaviside is zero every where....
703 if not isinstance(self.enrich, Heaviside):
704 grad_enr = self.enrich.evalGrad(evaltrisuv).flatten()
705 grad_enr = sp.sparse.csr_array((grad_enr, (np.arange(len(grad_enr)), np.repeat(
706 np.arange(len(grad_enr)//2), 2))), shape=(len(grad_enr), len(grad_enr)//2))
707 base_val = self.base.operator_dof2val_tri(evaltrisuv, self.rmap, self.frmap)[
708 :, self.enriched_vertices]
709 dNdx = dNdx + grad_enr.dot(base_val) # pylint: disable=C0103
710 return dNdx
712class XRidgeSpaceP1():
713 """ A special enrichment space for Ridge enrichment """
714 def __init__(self, name, mesh_cut, ls, ls_cut, rmap=None, dosvd=False) -> None:
715 self.dosvd = dosvd
716 self.mesh_cut = mesh_cut
717 self.name = name
718 self.rmap = meshcutdata2remap_helper(mesh_cut) if rmap is None else rmap
719 cut_elem = self.rmap.approNEQeval
720 self.enrich = Ridge(name+'_enrich', self.rmap,
721 ls[mesh_cut.tris], ls_cut[mesh_cut.new_tris])
722 self.base = FEMSpaceP1(name+'_base', SMesh(mesh_cut.xy, mesh_cut.tris),
723 vecdim=1, tris_map=self.rmap.parent_mapping)
724 self.frmap = mapping.RemapHelperFiltered(self.rmap.index_parent2sons, self.rmap.sons,
725 self.rmap.parent_mapping, self.rmap.son_mapping, cut_elem)
727 def scalar_size(self):
728 """ size (nbdofs) for a scalar space """
729 return len(self.mesh_cut.cut_edge_vertices)
731 def operator_dof2val_tri(self, evaltrisuv, rmap=None, raw=False): #pylint: disable=W0613
732 """ Return operator N that map the dof value to the field values at each point uv
733 for each triangle """
734 # if (rmap != self.rmap):
735 # raise ValueError("Trying to eval XRidgeSpaceP1 with an rmap not \
736 # conforming to the constructor's rmap")
737 nodes_uv = quad.T3_nodes().uv
738 enrval_atnodes:np.ndarray = self.enrich.eval(nodes_uv, frmap=self.frmap,
739 raw_data=True) # .reshape((-1,1))
740 nodes_uv = quad.T3_nodes().uv
741 npts = len(nodes_uv)
742 cut_vertices = self.mesh_cut.cut_vertices
743 cut_edges = self.mesh_cut.cut_edge_vertices
744 cut_tris1 = self.mesh_cut.new_tris
745 edge_funtmp = np.zeros((self.frmap.nbevals, npts, len(cut_edges)))
746 for i, v in enumerate(cut_vertices):
747 edge_funtmp[self.frmap.son2eval_index, :, i] = np.where(cut_tris1 == v, 1., 0.)[
748 self.frmap.sons]
749 Nat_nodes_nzelm = enrval_atnodes[:, :, np.newaxis]*edge_funtmp # pylint: disable=C0103
750 nenr = Nat_nodes_nzelm.shape[-1]
751 if self.dosvd:
752 N = Nat_nodes_nzelm.reshape((-1, len(cut_edges))) # pylint: disable=C0103
753 svd = np.linalg.svd(N)
754 print('shape : ', N.shape, 'singularvarlues ridge N:', svd[1])
755 if np.all(nodes_uv != evaltrisuv):
756 Nat_pt_nzelm = np.einsum('ps, ise -> ipe', # pylint: disable=C0103
757 FEMSpaceP1.N(evaltrisuv), Nat_nodes_nzelm)
758 else:
759 Nat_pt_nzelm = Nat_nodes_nzelm # pylint: disable=C0103
760 if raw is True:
761 N = Nat_pt_nzelm # pylint: disable=C0103
762 else:
763 eval_elements = np.arange(self.rmap.nbevals)[np.isin(
764 self.rmap.eval2appro, self.frmap.only_on_parents)]
765 npts = len(evaltrisuv)
766 row = np.repeat(np.repeat(npts*eval_elements, npts) +
767 np.tile(np.arange(npts), len(eval_elements)), nenr)
768 col = np.tile(np.arange(nenr), npts*len(eval_elements))
769 N = sp.sparse.csr_array( # pylint: disable=C0103
770 (Nat_pt_nzelm.flatten(), (row, col)), (self.rmap.nbevals*npts, nenr))
771 return N
773 def operator_dof2grad_tri(self, evaltrisuv, rmap = None): #pylint: disable=W0613
774 """ Return operator dNdx that map the dof value to the gradient of field values at
775 each point uv for each triangle """
776 # if (rmap != self.rmap):
777 # raise ValueError("Trying to eval XRidgeSpaceP1 with an rmap not \
778 # conforming to the constructor's rmap")
779 dim = 2
780 npts = len(evaltrisuv)
781 nodes_uv = quad.T3_nodes().uv
782 nshape = 3
783 NRidge = self.operator_dof2val_tri(nodes_uv, raw=True) # pylint: disable=C0103
784 nenr = NRidge.shape[-1]
785 NRidge = NRidge.reshape((-1, 3, nenr)) # pylint: disable=C0103
786 dNdx = np.zeros((self.frmap.nbevals, npts, dim, nshape)) # pylint: disable=C0103
787 dNdx[self.frmap.son2eval_index] = gradNTri(
788 FEMSpaceP1.dNdu, evaltrisuv, self.frmap.son_mapping.inv_F[self.rmap.sons])
789 dNRidgedx = np.einsum('ijkl, ilm ->ijkm', dNdx, NRidge) # pylint: disable=C0103
790 eval_elements = np.arange(self.rmap.nbevals)[np.isin(
791 self.rmap.eval2appro, self.frmap.only_on_parents)]
792 npts = len(evaltrisuv)
793 row = np.repeat(np.repeat(npts*dim*eval_elements, npts*dim) +
794 np.tile(np.arange(npts*dim), len(eval_elements)), nenr)
795 col = np.tile(np.arange(nenr), dim*npts*len(eval_elements))
796 dN = sp.sparse.csr_array( # pylint: disable=C0103
797 (dNRidgedx.flatten(), (row, col)), (self.rmap.nbevals*npts*dim, nenr))
798 return dN
801class SumSpace():
802 """ make a space as a sum of different space """
803 def __init__(self, list_of_space) -> None:
804 self.list_of_space = list_of_space
805 def operator_dof2val_tri(self, evaltrisuv, rmap=None):
806 """ Return operator N that map the dof value to the field values at each point uv
807 for each triangle """
808 N = [space.operator_dof2val_tri(evaltrisuv, rmap=rmap) # pylint: disable=C0103
809 for space in self.list_of_space]
810 return sp.sparse.hstack(N)
811 def operator_dof2grad_tri(self, evaltrisuv, rmap=None):
812 """ Return operator dNdx that map the dof value to the gradient of field values at
813 each point uv for each triangle """
814 dNs = [space.operator_dof2grad_tri(evaltrisuv, rmap=rmap) # pylint: disable=C0103
815 for space in self.list_of_space]
816 return sp.sparse.hstack(dNs)
818class SomeScalarField():
819 """ an example of scalar field for test purpose """
820 def __init__(self, f=1.) -> None:
821 self.f = f
822 def fun(self, xy):
823 """ eval the function at points xy """
824 return xy[:, 1]*np.cos(self.f*np.pi*xy[:, 0])
826 def grad(self, xy):
827 """ eval the gradient of function at points xy """
828 gx = -xy[:, 1]*self.f*np.pi*np.sin(self.f*np.pi*xy[:, 0])
829 gy = np.cos(self.f*np.pi*xy[:, 0])
830 return np.stack((gx, gy), 1)
833class LinearScalarField():
834 " A linear field in 2d : u(x,y) = a*x + b*y"
835 def __init__(self, a=3, b=4) -> None:
836 self.a, self.b = a, b
838 def fun(self, xy):
839 """ eval the function at points xy """
840 return self.a*xy[:, 0]+self.b*xy[:, 1]
842 def grad(self, xy):
843 """ eval the gradient of function at points xy """
844 gx = self.a*np.ones(len(xy))
845 gy = self.b*np.ones(len(xy))
846 return np.stack((gx, gy), 1)
849class QuadraticScalarField():
850 " A Quadratic field for example purpose "
851 def __init__(self, hess = np.eye(2) ) -> None:
852 self.hess = 0.5*(hess + hess.T)
854 def fun(self, xy):
855 """ eval the function at points xy """
856 x = xy[:, 0]
857 y = xy[:, 1]
858 h = self.hess
859 return .5*h[0,0]*x**2 + (h[0,1]+h[1,0])*x*y+ .5*h[1,1]*y*y
861 def grad(self, xy):
862 """ eval the gradient of function at points xy """
863 x = xy[:, 0]
864 y = xy[:, 1]
865 h = self.hess
866 gx = h[0,0]*x + (h[0,1]+h[1,0])*x*y
867 gy = (h[0,1]+h[1,0])*x + h[1,1]*y
868 return np.stack((gx, gy), 1)
871class SomeVectorField():
872 """ an example of vector field for test purpose """
873 def __init__(self, f) -> None:
874 self.f = f
876 def fun(self, xy):
877 """ eval the function at points xy """
878 f = self.f
879 x = xy[:, 0]
880 y = xy[:,1]
881 return np.column_stack((y*np.cos(f*np.pi*x), *np.cos(f*np.pi*y)))
883 def grad(self, xy):
884 """ eval the gradient of function at points xy """
885 f = self.f
886 gxx = -xy[:, 1]*f*np.pi*np.sin(f*np.pi*xy[:, 0])
887 gxy = np.cos(f*np.pi*xy[:, 0])
888 gyx = np.cos(f*np.pi*xy[:, 1])
889 gyy = -xy[:, 0]*f*np.pi*np.sin(f*np.pi*xy[:, 1])
890 return np.column_stack((gxx, gxy, gyx, gyy)).reshape((-1, 2, 2))
893class QuadraticVectorField():
894 """ an example of Vector field for test purpose """
895 def __init__(self) -> None:
896 pass
898 def fun(self, xy):
899 """ eval the function at points xy """
900 return 0.1*np.column_stack((xy[:, 0]**2, xy[:, 1]**2))
902 def grad(self, xy):
903 """ eval the gradient of function at points xy """
904 gxx = 0.1*2*xy[:, 0]
905 gxy = np.zeros(len(xy))
906 gyx = np.zeros(len(xy))
907 gyy = 0.1*2*xy[:, 1]
908 return np.column_stack((gxx, gxy, gyx, gyy)).reshape((-1, 2, 2))
911class LinearVectorField():
912 """ an example of Vector field for test purpose """
913 def __init__(self) -> None:
914 pass
916 def fun(self, xy):
917 """ eval the function at points xy """
918 return np.column_stack((xy[:, 0], 2.*xy[:, 1]))
920 def grad(self, xy):
921 """ eval the gradient of function at points xy """
922 gxx = np.ones(len(xy))
923 gxy = np.zeros(len(xy))
924 gyx = np.zeros(len(xy))
925 gyy = 2.*np.ones(len(xy))
926 return np.column_stack((gxx, gxy, gyx, gyy)).reshape((-1, 2, 2))