Coverage for src/methodsnm/formint.py: 89%

1from abc import ABC, abstractmethod 

2import numpy as np 

3from numpy import array, einsum 

4from methodsnm.intrule import * 

5from methodsnm.fe import * 

6from numpy.linalg import det, inv 

7from methodsnm.meshfct import ConstantFunction 

8 

9class FormIntegral(ABC): 

10 

11 @abstractmethod 

12 def __init__(self): 

13 raise NotImplementedError("Not implemented") 

14 

15class LinearFormIntegral(FormIntegral): 

16 

17 @abstractmethod 

18 def __init__(self): 

19 raise NotImplementedError("Not implemented") 

20 

21 def compute_element_vector(self, fe, trafo): 

22 raise NotImplementedError("Not implemented") 

23 

24class SourceIntegral(LinearFormIntegral): 

25 

26 def __init__(self, coeff=ConstantFunction(1)): 

27 self.coeff = coeff 

28 

29 def compute_element_vector(self, fe, trafo, intrule = None): 

30 if intrule is None: 

31 intrule = select_integration_rule(2*fe.order, fe.eltype) 

32 shapes = fe.evaluate(intrule.nodes) 

33 coeffs = self.coeff.evaluate(intrule.nodes, trafo) 

34 weights = [w*c*abs(det(trafo.jacobian(ip))) for ip,w,c in zip(intrule.nodes,intrule.weights,coeffs)] 

35 return np.dot(shapes.T, weights) 

36 

37class BilinearFormIntegral(FormIntegral): 

38 

39 @abstractmethod 

40 def __init__(self): 

41 raise NotImplementedError("Not implemented") 

42 

43 def compute_element_matrix(self, fe, trafo): 

44 raise NotImplementedError("Not implemented") 

45 

46class MassIntegral(BilinearFormIntegral): 

47 

48 def __init__(self, coeff=ConstantFunction(1)): 

49 self.coeff = coeff 

50 

51 def compute_element_matrix(self, fe_test, fe_trial, trafo, intrule = None): 

52 if intrule is None: 

53 if fe_test.eltype != fe_trial.eltype: 

54 raise Exception("Finite elements must have the same el. type") 

55 intrule = select_integration_rule(fe_test.order + fe_trial.order, fe_test.eltype) 

56 shapes_test = fe_test.evaluate(intrule.nodes) 

57 shapes_trial = fe_test.evaluate(intrule.nodes) 

58 coeffs = self.coeff.evaluate(intrule.nodes, trafo) 

59 F = trafo.jacobian(intrule.nodes) 

60 adetF = array([abs(det(F[i,:,:])) for i in range(F.shape[0])]) 

61 ret = einsum("ij,ik,i,i,i->jk", shapes_test, shapes_trial, adetF, coeffs, intrule.weights) 

62 return ret 

63 

64class LaplaceIntegral(BilinearFormIntegral): 

65 

66 def __init__(self, coeff=ConstantFunction(1)): 

67 self.coeff = coeff 

68 

69 def compute_element_matrix(self, fe_test, fe_trial, trafo, intrule = None): 

70 if intrule is None: 

71 if fe_test.eltype != fe_trial.eltype: 

72 raise Exception("Finite elements must have the same el. type") 

73 intrule = select_integration_rule(fe_test.order + fe_trial.order, fe_test.eltype) 

74 dshapes_ref_test = fe_test.evaluate(intrule.nodes, deriv=True) 

75 dshapes_ref_trial = fe_test.evaluate(intrule.nodes, deriv=True) 

76 F = trafo.jacobian(intrule.nodes) 

77 invF = array([inv(F[i,:,:]) for i in range(F.shape[0])]) 

78 adetF = array([abs(det(F[i,:,:])) for i in range(F.shape[0])]) 

79 coeffs = self.coeff.evaluate(intrule.nodes, trafo) 

80 ret = einsum("ijk,imn,ijl,iml,i,i,i->kn", dshapes_ref_test, dshapes_ref_trial, invF, invF, adetF, coeffs, intrule.weights) 

81 return ret