Coverage for pyxfem/solidmechanics/boundary_conditions.py: 37%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Created on April 10th 07:43:08 2024 

5 

6@author: nchevaug 

7""" 

8from __future__ import annotations 

9from warnings import warn 

10from collections.abc import Iterable 

11from typing import Callable, Mapping, Dict 

12import numpy as np 

13from pyxfem.core.config import FloatArray, IndexArray, IndexType 

14import pyxfem.core.function_spaces as fs 

15import pyxfem.core.mesh_simplex_one_level as sm 

16import pyxfem.core.quadrature as quad 

17 

18def get_space_entity(edge_groupe_names:Iterable[str], 

19 mesh:sm.SMesh)-> Dict[str, IndexArray]: 

20 ''' get the space's edge id corresponding to the edge groupe name''' 

21 edge_groupe_names2eid: Dict[str, IndexArray] = {} 

22 for name in edge_groupe_names: 

23 edge_groupe_names2eid[name] = mesh.get_edgegroupe(name)[0] 

24 return edge_groupe_names2eid 

25 

26def get_dofs_and_interpolation_points( 

27 mesh: sm.SMesh, 

28 edge_groupe_names: Iterable[str], 

29 space)->tuple[Dict[str, IndexArray], Dict[str, FloatArray]]: 

30 ''' get the space's degree of freedom id and the corresponding interpolation point 

31 corresponding to the edge_groupe_names (alist of names). Usefull for Dirichlet BC for example''' 

32 groupe_name2dofid:Dict[str, IndexArray]= {} 

33 groupe_name2dofxy:Dict[str, FloatArray] = {} 

34 func_name = get_dofs_and_interpolation_points.__name__ 

35 # if space is not None : 

36 # interpat = space.interpolatoryAt 

37 # vid2dofid = space.vertexid2dofid 

38 # if vid2dofid is None and interpat=='vertex': 

39 # raise ValueError() 

40 for phys_name in edge_groupe_names: 

41 eid, ___ = mesh.get_edgegroupe(phys_name) 

42 if eid.size == 0: 

43 warn('In '+func_name+' no edge found in groupe '+phys_name+'.') 

44 if space.interpolatoryAt == 'vertex': 

45 vertices = np.unique((mesh.e2v()[eid]).flatten()) 

46 coord = mesh.vcoord()[vertices] 

47 dofsid = space.vertexid2dofid(vertices) 

48 groupe_name2dofid[phys_name] = dofsid 

49 groupe_name2dofxy[phys_name] = coord #.reshape((-1, 2)) 

50 elif space.interpolatoryAt == 'midEdge': 

51 groupe_name2dofid[phys_name] = space.edgeid2dofid(eid) 

52 groupe_name2dofxy[phys_name] = mesh.ecog()[eid] 

53 else: raise ValueError(func_name+" not coded for space.interpolatoryAt.") 

54 return groupe_name2dofid, groupe_name2dofxy 

55 

56def set_values(vxy: FloatArray, 

57 evaluator: float|Callable[[FloatArray], FloatArray] 

58 )->FloatArray: 

59 """ A simple function that return values associated to the points in array vxy 

60 

61 according to the evaluator. The evaluator can be a scalar or a function 

62 """ 

63 func_name = set_values.__name__ 

64 if isinstance(evaluator, float): 

65 values = np.ones(len(vxy))*evaluator 

66 elif callable(evaluator): 

67 values = evaluator(vxy) 

68 else: 

69 raise ValueError('evaluator of '+str(type(evaluator))+' unknown in '+func_name+'.') 

70 return values 

71 

72dirname2dirindex = {'x':0, 'y':1, 'z':2} 

73def compute_dirichlet(space, 

74 physname2dirichlet, 

75 groupename2dofid:Mapping[str, IndexArray], 

76 groupename2dofxy:Mapping[str, FloatArray] 

77 )-> tuple[IndexArray, IndexArray, FloatArray]: 

78 ''' return the list of free dofs ids, fixed dofs ids and fixed dofs values ''' 

79 func_name = compute_dirichlet.__name__ 

80 fixed_values = np.array([], dtype=np.float64) 

81 fixed_dofs = np.array([], dtype=np.int64) 

82 for phys_name, dir2value in physname2dirichlet.items(): 

83 dofid = groupename2dofid[phys_name] 

84 dofxy = groupename2dofxy[phys_name] 

85 for dirname, value in dir2value.items(): 

86 direction = dirname2dirindex[dirname] 

87 #print('ggg',dirname, direction) 

88 if direction is None: 

89 raise ValueError('directionName '+dirname+' unknown in '+func_name+'.') 

90 fixed_dofs = np.append(fixed_dofs, dofid[:, direction]).astype(np.int64) 

91 fixed_values = np.append(fixed_values, set_values(dofxy, value)).astype(np.float64) 

92 #nbfixed = len(fixed_dofs) 

93 fixed_dofs, index = np.unique(fixed_dofs, return_index=True) 

94 # if nbfixed != len(fixed_dofs): 

95 # warn(' Some dofs where fixed Twice in '+func_name+'.') 

96 #nbfixed = len(fixed_dofs) 

97 fixed_values = np.array(fixed_values[index], dtype = np.float64) 

98 free_dofs = np.setdiff1d(np.arange(space.size()), fixed_dofs).astype(np.int64) 

99 free_dofs = np.asarray(free_dofs, dtype = IndexType) 

100 fixed_dofs = np.asarray(fixed_dofs, dtype = IndexType) 

101 return free_dofs, fixed_dofs, fixed_values 

102 

103def compute_neumann(space:fs.FEMSpace, phys_name2neumann): 

104 ''' return the vector F such as F.T.dot(U) = int_neumann ( t.u_h) dGamma_h ''' 

105 fun_name = compute_neumann.__name__ 

106 F = np.zeros((space.scalar_size(), space.vecdim)) # pylint: disable=C0103 

107 edgequad = quad.Edge_gauss(0) 

108 for phys_name, dir2value in phys_name2neumann.items(): 

109 eids, ___ = space.m.get_edgegroupe(phys_name) 

110 N = space.operator_dof2val_edge(eids, edgequad.s) # pylint: disable=C0103 

111 xy = space.m.ecog(eids) 

112 l = space.m.elength(eids) 

113 t = np.zeros((eids.size, 2)) 

114 for direction_name, value in dir2value.items(): 

115 direction = dirname2dirindex[direction_name] 

116 if direction is None: 

117 raise ValueError('directionName '+direction_name+' unknown in '+fun_name) 

118 t[:, direction] = l*set_values(xy, value) 

119 F += (edgequad.w*N.T.dot(t.flatten())).reshape((-1, 2)) # pylint: disable=C0103 

120 return F.flatten() 

121 

122def compute_dirichlet_fem_p1_vect2d(nvertex, nodex, valuex, nodey, valuey): 

123 ''' more direct way to compute dirichlet bc if node list is already known''' 

124 fixed_dofs, index = np.unique(np.append(nodex*2, nodey*2+1), return_index=True) 

125 fixed_dofs = fixed_dofs.astype(np.int64) 

126 fixed_values = np.append(valuex, valuey)[index] 

127 free_dofs = np.setdiff1d(np.arange(2*nvertex), fixed_dofs).astype(np.int64) 

128 return free_dofs, fixed_dofs, fixed_values