Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/block_assemble.py: 70%

1from __future__ import division 

2from __future__ import absolute_import 

3from builtins import map 

4from builtins import range 

5from . import * 

6from .block_util import block_tensor, isscalar, wrap_in_list, create_vec_from 

7 

8def block_assemble(lhs, rhs=None, bcs=None, 

9 symmetric=False, signs=None, symmetric_mod=None): 

10 """ 

11 Assembles block matrices, block vectors or block systems. 

12 Input can be arrays of variational forms or block matrices/vectors. 

13 

14 Arguments: 

15 

16 symmetric : Boundary conditions are applied so that symmetry of the system 

17 is preserved. If only the left hand side of the system is given, 

18 then a matrix represententing the rhs corrections is returned 

19 along with a symmetric matrix. NOTE: Symmetric assembly not 

20 supported in parallel (MPI). 

21 

22 symmetric_mod : Matrix describing symmetric corrections for assembly of the 

23 of the rhs of a variational system. 

24 

25 signs : An array to specify the signs of diagonal blocks. The sign 

26 of the blocks are computed if the argument is not provided. 

27 """ 

28 error_msg = {'incompatibility' : 'A and b do not have compatible dimensions.', 

29 'symm_mod error' : 'symmetric_mod argument only accepted when assembling a vector', 

30 'not square' : 'A must be square for symmetric assembling', 

31 'invalid bcs' : 'Expecting a list or list of lists of DirichletBC.', 

32 'invalid signs' : 'signs should be a list of length n containing only 1 or -1'} 

33 # Check arguments 

34 from numpy import ndarray 

35 has_rhs = True if isinstance(rhs, ndarray) else rhs != None 

36 has_lhs = True if isinstance(rhs, ndarray) else rhs != None 

37 

38 if has_lhs and has_rhs: 

39 A, b = list(map(block_tensor,[lhs,rhs])) 

40 n, m = A.blocks.shape 

41 if not ( isinstance(b,block_vec) and len(b.blocks) is m): 41 ↛ 42line 41 didn't jump to line 42, because the condition on line 41 was never true

42 raise TypeError(error_msg['incompatibility']) 

43 else: 

44 A, b = block_tensor(lhs), None 

45 if isinstance(A,block_vec): 

46 A, b = None, A 

47 n, m = 0, len(b.blocks) 

48 else: 

49 n,m = A.blocks.shape 

50 if A and symmetric and (m is not n): 50 ↛ 51line 50 didn't jump to line 51, because the condition on line 50 was never true

51 raise RuntimeError(error_msg['not square']) 

52 if symmetric_mod and ( A or not b ): 52 ↛ 53line 52 didn't jump to line 53, because the condition on line 52 was never true

53 raise RuntimeError(error_msg['symmetric_mod error']) 

54 # First assemble everything needing assembling. 

55 from dolfin import assemble 

56 assemble_if_form = lambda x: assemble(x, keep_diagonal=True) if _is_form(x) else x 

57 if A: 

58 A.blocks.flat[:] = list(map(assemble_if_form,A.blocks.flat)) 

59 if b: 

60 #b.blocks.flat[:] = map(assemble_if_form, b.blocks.flat) 

61 b = block_vec(list(map(assemble_if_form, b.blocks.flat))) 

62 # If there are no boundary conditions then we are done. 

63 if bcs is None: 

64 return [A,b] if (A and b) else A or b 

65 

66 # Apply supplied RHS BCs if we are only assembling the right hand side. If we are 

67 # assembling A anyway, we use that for symmetry preservation instead. 

68 if A is None and symmetric_mod: 68 ↛ 69line 68 didn't jump to line 69, because the condition on line 68 was never true

69 symmetric_mod.apply(b) 

70 return b 

71 

72 # check if arguments are forms, in which case bcs have to be split 

73 from ufl import Form 

74 lhs_bcs = (block_bc.from_mixed if isinstance(lhs, Form) else block_bc)(bcs, symmetric=symmetric, signs=signs) 

75 

76 result = [] 

77 if A: 

78 rhs_bcs = lhs_bcs.apply(A) 

79 result.append(A) 

80 else: 

81 rhs_bcs = lhs_bcs.rhs(None) 

82 if symmetric and A: 82 ↛ 83line 82 didn't jump to line 83, because the condition on line 82 was never true

83 result.append(rhs_bcs) 

84 if b: 

85 rhs_bcs.apply(b) 

86 result.append(b) 

87 return result[0] if len(result)==1 else result 

88 

89 

90def block_symmetric_assemble(forms, bcs): 

91 return block_assemble(forms,bcs=bcs,symmetric=True) 

92 

93def _is_form(form): 

94 from dolfin import Form as cpp_Form 

95 from ufl.form import Form as ufl_Form 

96 return isinstance(form, (cpp_Form, ufl_Form)) 

97 

98def _new_square_matrix(bc, val): 

99 from dolfin import TrialFunction, TestFunction, FunctionSpace 

100 from dolfin import assemble, Constant, inner, dx 

101 import numpy 

102 

103 V = FunctionSpace(bc.function_space()) 

104 

105 u,v = TrialFunction(V),TestFunction(V) 

106 Z = assemble(Constant(0)*inner(u,v)*dx) 

107 if val != 0.0: 

108 lrange = list(range(*Z.local_range(0))) 

109 idx = numpy.ndarray(len(lrange), dtype=numpy.intc) 

110 idx[:] = lrange 

111 Z.ident(idx) 

112 if val != 1.0: 

113 Z *= val 

114 return Z