Coverage for src/methodsnm/trafo.py: 81%

1from abc import ABC, abstractmethod 

2import numpy as np 

3from numpy import array 

4 

5from methodsnm.mesh import Mesh 

6 

7class Transformation(ABC): 

8 dim_domain = None 

9 dim_range = None 

10 eltype = None 

11 

12 def __init__(self, interval): 

13 pass 

14 

15 @abstractmethod 

16 def _map(self, ip): 

17 raise NotImplementedError("Not implemented") 

18 

19 def _map_array(self, ips): 

20 ret = np.empty(ips.shape) 

21 for i in range(ips.shape[0]): 

22 ret[i,:] = self._map(ips[i]) 

23 return ret 

24 

25 @abstractmethod 

26 def _jacobian(self, ip): 

27 raise NotImplementedError("Not implemented") 

28 

29 def _jacobian_array(self, ips): 

30 ret = np.empty((ips.shape[0], self.dim_range, self.dim_domain)) 

31 for i in range(ips.shape[0]): 

32 ret[i,:,:] = self._jacobian(ips[i]) 

33 return ret 

34 

35 def map(self, ip): 

36 if ip.ndim == 1: 

37 return self._map(ip) 

38 else: 

39 return self._map_array(ip) 

40 

41 def jacobian(self, ip): 

42 if ip.ndim == 1: 

43 return self._jacobian(ip) 

44 else: 

45 return self._jacobian_array(ip) 

46 

47 def __call__(self, ip): 

48 return self.map(ip) 

49 

50class ElementTransformation(Transformation): 

51 mesh = None 

52 elnr = None 

53 def __init__(self, mesh, elnr): 

54 self.elnr = elnr 

55 self.mesh = mesh 

56 

57class IntervalTransformation(ElementTransformation): 

58 # This class represents a transformation from the interval [0,1] to the given interval. 

59 def __init__(self, mesh, elnr): 

60 super().__init__(mesh, elnr) 

61 self.interval = tuple(mesh.points[mesh.elements()[elnr]]) 

62 self.dim_range = 1 

63 self.dim_domain = 1 

64 self.eltype = "segment" 

65 

66 def _map(self, ip): 

67 a,b = self.interval 

68 return a+(b-a)*ip[0] 

69 

70 def _jacobian(self, ip): 

71 a,b = self.interval 

72 return np.array([[b-a]]) 

73 

74class TriangleTransformation(ElementTransformation): 

75 points = None 

76 def calculate_jacobian(self): 

77 a,b,c = self.points 

78 return array([b-a,c-a]).T 

79 

80 def __init__(self, mesh_or_points, elnr=None): 

81 if isinstance(mesh_or_points, Mesh): 

82 mesh = mesh_or_points 

83 if elnr is None: 

84 raise ValueError("TriangleTransformation needs an element number") 

85 super().__init__(mesh, elnr) 

86 self.points = tuple(mesh.points[mesh.elements()[elnr]]) 

87 else: 

88 super().__init__(mesh = None, elnr = -1) 

89 self.points = mesh_or_points 

90 self.jac = self.calculate_jacobian() 

91 self.dim_range = 2 

92 self.dim_domain = 2 

93 self.eltype = "triangle" 

94 

95 def _map(self, ip): 

96 a,b,c = self.points 

97 return a+(b-a)*ip[0]+(c-a)*ip[1] 

98 

99 def _jacobian(self, ip): 

100 return self.jac 

101