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

1from abc import ABC, abstractmethod 

2import numpy as np 

3from numpy import array 

4from methodsnm.trafo import * 

5from methodsnm.mesh import * 

6from scipy import sparse 

7 

8class Form(ABC): 

9 integrals = None 

10 fes = None 

11 def __init__(self): 

12 raise NotImplementedError("Not implemented") 

13 

14 def add_integrator(self, integrator): 

15 self.integrals.append(integrator) 

16 

17 def __iadd__(self, integrator): 

18 self.add_integrator(integrator) 

19 return self 

20 

21class LinearForm(Form): 

22 

23 vector = None 

24 def __init__(self, fes=None): 

25 self.fes = fes 

26 self.integrals = [] 

27 

28 def assemble(self): 

29 self.vector = np.zeros(self.fes.ndof) 

30 mesh = self.fes.mesh 

31 for elnr, verts in enumerate(mesh.elements()): 

32 trafo = mesh.trafo(elnr) 

33 fe = self.fes.finite_element(elnr) 

34 dofs = self.fes.element_dofs(elnr) 

35 for integral in self.integrals: 

36 self.vector[dofs] += integral.compute_element_vector(fe, trafo) 

37 

38 

39class BilinearForm(Form): 

40 

41 matrix = None 

42 def __init__(self, fes=None, fes_test=None, fes_trial=None): 

43 if fes is None and (fes_test is None or fes_trial is None) or all([f is not None for f in [fes,fes_test,fes_trial]]): 

44 raise Exception("Invalid arguments, specify either `fes` or `fes_test` and `fes_trial`") 

45 if fes_test is None: 

46 fes_test = fes 

47 fes_trial = fes 

48 

49 self.fes_trial = fes_trial 

50 self.fes_test = fes_test 

51 self.fes = fes_test 

52 self.integrals = [] 

53 

54 def assemble(self): 

55 self.matrix = sparse.lil_matrix((self.fes_test.ndof, self.fes_trial.ndof)) 

56 mesh = self.fes.mesh 

57 for elnr, verts in enumerate(mesh.elements()): 

58 trafo = mesh.trafo(elnr) 

59 fe_test = self.fes_test.finite_element(elnr) 

60 dofs_test = self.fes_test.element_dofs(elnr) 

61 fe_trial = self.fes_trial.finite_element(elnr) 

62 dofs_trial = self.fes_trial.element_dofs(elnr) 

63 elmat = np.zeros((len(dofs_test), len(dofs_trial))) 

64 for integral in self.integrals: 

65 elmat += integral.compute_element_matrix(fe_test, fe_trial, trafo) 

66 for i, dofi in enumerate(dofs_test): 

67 for j, dofj in enumerate(dofs_trial): 

68 self.matrix[dofi,dofj] += elmat[i,j] 

69 self.matrix = self.matrix.tocsr() 

70 

71 

72from methodsnm.intrule import select_integration_rule 

73from numpy.linalg import det 

74def compute_difference_L2(uh, uex, mesh, intorder=5): 

75 sumint = 0 

76 for elnr in range(len(mesh.elements())): 

77 trafo = mesh.trafo(elnr) 

78 intrule = select_integration_rule(intorder, trafo.eltype) 

79 uhvals = uh.evaluate(intrule.nodes, trafo) 

80 uexvals = uex.evaluate(intrule.nodes, trafo) 

81 diff = np.zeros(len(intrule.nodes)) 

82 diff = (uhvals - uexvals)**2 

83 F = trafo.jacobian(intrule.nodes) 

84 w = array([abs(det(F[i,:,:])) * intrule.weights[i] for i in range(F.shape[0])]) 

85 sumint += np.dot(w, diff) 

86 return np.sqrt(sumint)