#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on April 10th 07:43:08 2024
@author: nchevaug
"""
from __future__ import annotations
from warnings import warn
from collections.abc import Iterable
from typing import Callable, Mapping, Dict
import numpy as np
from pyxfem.core.config import FloatArray, IndexArray, IndexType
import pyxfem.core.function_spaces as fs
import pyxfem.core.mesh_simplex_one_level as sm
import pyxfem.core.quadrature as quad
[docs]
def get_space_entity(edge_groupe_names:Iterable[str],
mesh:sm.SMesh)-> Dict[str, IndexArray]:
''' get the space's edge id corresponding to the edge groupe name'''
edge_groupe_names2eid: Dict[str, IndexArray] = {}
for name in edge_groupe_names:
edge_groupe_names2eid[name] = mesh.get_edgegroupe(name)[0]
return edge_groupe_names2eid
[docs]
def get_dofs_and_interpolation_points(
mesh: sm.SMesh,
edge_groupe_names: Iterable[str],
space)->tuple[Dict[str, IndexArray], Dict[str, FloatArray]]:
''' get the space's degree of freedom id and the corresponding interpolation point
corresponding to the edge_groupe_names (alist of names). Usefull for Dirichlet BC for example'''
groupe_name2dofid:Dict[str, IndexArray]= {}
groupe_name2dofxy:Dict[str, FloatArray] = {}
func_name = get_dofs_and_interpolation_points.__name__
# if space is not None :
# interpat = space.interpolatoryAt
# vid2dofid = space.vertexid2dofid
# if vid2dofid is None and interpat=='vertex':
# raise ValueError()
for phys_name in edge_groupe_names:
eid, ___ = mesh.get_edgegroupe(phys_name)
if eid.size == 0:
warn('In '+func_name+' no edge found in groupe '+phys_name+'.')
if space.interpolatoryAt == 'vertex':
vertices = np.unique((mesh.e2v()[eid]).flatten())
coord = mesh.vcoord()[vertices]
dofsid = space.vertexid2dofid(vertices)
groupe_name2dofid[phys_name] = dofsid
groupe_name2dofxy[phys_name] = coord #.reshape((-1, 2))
elif space.interpolatoryAt == 'midEdge':
groupe_name2dofid[phys_name] = space.edgeid2dofid(eid)
groupe_name2dofxy[phys_name] = mesh.ecog()[eid]
else: raise ValueError(func_name+" not coded for space.interpolatoryAt.")
return groupe_name2dofid, groupe_name2dofxy
[docs]
def set_values(vxy: FloatArray,
evaluator: float|Callable[[FloatArray], FloatArray]
)->FloatArray:
""" A simple function that return values associated to the points in array vxy
according to the evaluator. The evaluator can be a scalar or a function
"""
func_name = set_values.__name__
if isinstance(evaluator, float):
values = np.ones(len(vxy))*evaluator
elif callable(evaluator):
values = evaluator(vxy)
else:
raise ValueError('evaluator of '+str(type(evaluator))+' unknown in '+func_name+'.')
return values
dirname2dirindex = {'x':0, 'y':1, 'z':2}
[docs]
def compute_dirichlet(space,
physname2dirichlet,
groupename2dofid:Mapping[str, IndexArray],
groupename2dofxy:Mapping[str, FloatArray]
)-> tuple[IndexArray, IndexArray, FloatArray]:
''' return the list of free dofs ids, fixed dofs ids and fixed dofs values '''
func_name = compute_dirichlet.__name__
fixed_values = np.array([], dtype=np.float64)
fixed_dofs = np.array([], dtype=np.int64)
for phys_name, dir2value in physname2dirichlet.items():
dofid = groupename2dofid[phys_name]
dofxy = groupename2dofxy[phys_name]
for dirname, value in dir2value.items():
direction = dirname2dirindex[dirname]
#print('ggg',dirname, direction)
if direction is None:
raise ValueError('directionName '+dirname+' unknown in '+func_name+'.')
fixed_dofs = np.append(fixed_dofs, dofid[:, direction]).astype(np.int64)
fixed_values = np.append(fixed_values, set_values(dofxy, value)).astype(np.float64)
#nbfixed = len(fixed_dofs)
fixed_dofs, index = np.unique(fixed_dofs, return_index=True)
# if nbfixed != len(fixed_dofs):
# warn(' Some dofs where fixed Twice in '+func_name+'.')
#nbfixed = len(fixed_dofs)
fixed_values = np.array(fixed_values[index], dtype = np.float64)
free_dofs = np.setdiff1d(np.arange(space.size()), fixed_dofs).astype(np.int64)
free_dofs = np.asarray(free_dofs, dtype = IndexType)
fixed_dofs = np.asarray(fixed_dofs, dtype = IndexType)
return free_dofs, fixed_dofs, fixed_values
[docs]
def compute_neumann(space:fs.FEMSpace, phys_name2neumann):
''' return the vector F such as F.T.dot(U) = int_neumann ( t.u_h) dGamma_h '''
fun_name = compute_neumann.__name__
F = np.zeros((space.scalar_size(), space.vecdim)) # pylint: disable=C0103
edgequad = quad.Edge_gauss(0)
for phys_name, dir2value in phys_name2neumann.items():
eids, ___ = space.m.get_edgegroupe(phys_name)
N = space.operator_dof2val_edge(eids, edgequad.s) # pylint: disable=C0103
xy = space.m.ecog(eids)
l = space.m.elength(eids)
t = np.zeros((eids.size, 2))
for direction_name, value in dir2value.items():
direction = dirname2dirindex[direction_name]
if direction is None:
raise ValueError('directionName '+direction_name+' unknown in '+fun_name)
t[:, direction] = l*set_values(xy, value)
F += (edgequad.w*N.T.dot(t.flatten())).reshape((-1, 2)) # pylint: disable=C0103
return F.flatten()
[docs]
def compute_dirichlet_fem_p1_vect2d(nvertex, nodex, valuex, nodey, valuey):
''' more direct way to compute dirichlet bc if node list is already known'''
fixed_dofs, index = np.unique(np.append(nodex*2, nodey*2+1), return_index=True)
fixed_dofs = fixed_dofs.astype(np.int64)
fixed_values = np.append(valuex, valuey)[index]
free_dofs = np.setdiff1d(np.arange(2*nvertex), fixed_dofs).astype(np.int64)
return free_dofs, fixed_dofs, fixed_values