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

1import dolfin 

2from .block_mat import block_mat 

3from .block_vec import block_vec 

4from .block_util import wrap_in_list, create_diagonal_matrix 

5from .helpers import supports_mpi 

6from .splitting import split_bcs 

7import itertools 

8import numpy 

9 

10class block_bc: 

11 """ 

12 This class applies Dirichlet BCs to a block matrix. It is not a block 

13 operator itself. 

14 

15 Creating bcs: 

16 

17 >>> bcs = block_bc([...], symmetric=...) 

18 

19 Basic inplace usage, with static BCs: 

20 

21 >>> rhs_bcs = bcs.apply(A) 

22 >>> rhs_bcs.apply(b) 

23 

24 with the shortcut 

25 

26 >>> bcs.apply(A,b) 

27 

28 If boundary conditions are applied multiple times, it may be useful 

29 to keep the original assembled mat/vec unchanged (for example, to 

30 avoid having to re-assemble the RHS if the BC is time dependent but 

31 the form itself isn't): 

32 

33 >>> Ainv = conjgrad(bcs(A)) 

34 >>> while True: 

35 >>> source.t = t 

36 >>> x = Ainv * bcs(b) 

37 """ 

38 def __init__(self, bcs, symmetric=False, signs=None, subspace_bcs=None): 

39 # Clean up self, and check arguments 

40 self.symmetric = symmetric 

41 bcs = [wrap_in_list(bc, dolfin.DirichletBC) for bc in bcs] 

42 if subspace_bcs is not None: 

43 subspace_bcs = split_bcs(subspace_bcs, None) 

44 combined_bcs = [] 

45 for ss,ns in itertools.zip_longest(subspace_bcs, bcs): 

46 combined_bcs.append([]) 

47 if ns is not None: combined_bcs[-1] += ns 

48 if ss is not None: combined_bcs[-1] += ss 

49 self.bcs = combined_bcs 

50 else: 

51 self.bcs = bcs 

52 self.signs = signs or [1]*len(self.bcs) 

53 if not all(s in [-1,1] for s in self.signs): 53 ↛ 54line 53 didn't jump to line 54, because the condition on line 53 was never true

54 raise ValueError('signs should be a list of length n containing only 1 or -1') 

55 supports_mpi(not symmetric, 'Symmetric BCs not supported in parallel') 

56 

57 @classmethod 

58 def from_mixed(cls, bcs, *args, **kwargs): 

59 return cls([], *args, **kwargs, subspace_bcs=bcs) 

60 

61 def __call__(self, other, A=None): 

62 if isinstance(other, block_mat): 

63 self._A = other # opportunistic, expecting a block_vec later 

64 A = other.copy() 

65 self.apply(A) 

66 return A 

67 elif isinstance(other, block_vec): 

68 if A is None: 

69 A = getattr(self, '_A', None) # opportunistic 

70 if self.symmetric and A is None: 

71 raise ValueError('A not available. Call with A first, or pass A=A.') 

72 rhs = self.rhs(A) 

73 return rhs(other) 

74 else: 

75 raise TypeError('BCs can only act on block_mat (A) or block_vec (b)') 

76 

77 def apply(self, A, b=None): 

78 """ 

79 Apply BCs to a block_mat LHS, and return an object which modifies the 

80 corresponding block_vec RHS. Typical use: 

81 

82 b_bcs = A_bcs.apply(A); b_bcs.apply(b) 

83 """ 

84 if not isinstance(A, block_mat): 84 ↛ 85line 84 didn't jump to line 85, because the condition on line 84 was never true

85 raise RuntimeError('A is not a block matrix') 

86 

87 A_orig = A.copy() if self.symmetric else A 

88 self._apply(A) 

89 

90 # Create rhs_bc with a copy of A before any symmetric modifications 

91 rhs = self.rhs(A_orig) 

92 if b is not None: 92 ↛ 93line 92 didn't jump to line 93, because the condition on line 92 was never true

93 rhs.apply(b) 

94 return rhs 

95 

96 def rhs(self, A): 

97 return block_rhs_bc(self.bcs, A, symmetric=self.symmetric, signs=self.signs) 

98 

99 def _apply(self, A): 

100 if self.symmetric: 

101 # dummy vec, required by dolfin api -- we don't use this 

102 # corrections directly since it may be time dependent 

103 b = A.create_vec(dim=0) 

104 for i,bcs in enumerate(self.bcs): 

105 if bcs: 

106 for j in range(len(A)): 

107 if i==j: 

108 if numpy.isscalar(A[i,i]): 

109 # Convert to a diagonal matrix, so that the individual rows can be modified 

110 A[i,i] = create_diagonal_matrix(dolfin.FunctionSpace(bc.function_space()), A[i,i]) 

111 if self.symmetric: 

112 for bc in bcs: 

113 bc.zero_columns(A[i,i], b[i], self.signs[i]) 

114 else: 

115 if self.signs[i] != 1: 115 ↛ 116line 115 didn't jump to line 116, because the condition on line 115 was never true

116 A[i,i] *= 1/self.signs[i] 

117 for bc in bcs: 

118 bc.apply(A[i,i]) 

119 if self.signs[i] != 1: 119 ↛ 120line 119 didn't jump to line 120, because the condition on line 119 was never true

120 A[i,i] *= self.signs[i] 

121 else: 

122 if numpy.isscalar(A[i,j]): 

123 if A[i,j] != 0.0: 123 ↛ 124line 123 didn't jump to line 124, because the condition on line 123 was never true

124 dolfin.error("can't modify nonzero scalar off-diagonal block (%d,%d)" % (i,j)) 

125 else: 

126 for bc in bcs: 

127 bc.zero(A[i,j]) 

128 if self.symmetric: 

129 if numpy.isscalar(A[j,i]): 

130 if A[j,i] != 0.0: 130 ↛ 131line 130 didn't jump to line 131, because the condition on line 130 was never true

131 dolfin.error("can't modify nonzero scalar off-diagonal block (%d,%d)" % (i,j)) 

132 else: 

133 for bc in bcs: 

134 bc.zero_columns(A[j,i], b[j]) 

135 

136class block_rhs_bc: 

137 """ 

138 This class applies Dirichlet BCs to a block block vector. It can be used 

139 as a block operator; but since it is nonlinear, it can only operate 

140 directly on a block vector, and not combine with other operators. 

141 """ 

142 def __init__(self, bcs, A, symmetric, signs): 

143 if symmetric: 

144 assert A is not None 

145 self.bcs = bcs 

146 self.A = A 

147 self.symmetric = symmetric 

148 self.signs = signs 

149 

150 def apply(self, b): 

151 """Apply Dirichlet boundary conditions statically. In the 

152 inhomogeneous symmetric case b is changed globally; if BCs are mutable 

153 (time dependent for example), it is then more convenient to use the 

154 callable form -- rhs_bc(b) -- to preserve the original contents of b 

155 for repeated application without reassembly. 

156 """ 

157 if not isinstance(b, block_vec): 157 ↛ 158line 157 didn't jump to line 158, because the condition on line 157 was never true

158 raise TypeError('not a block vector') 

159 

160 try: 

161 b.allocate(self.A, dim=0, alternative_templates=[self.bcs]) 

162 except Exception: 

163 raise ValueError('Failed to allocate block vector, call b.allocate(something) first') 

164 

165 # Apply lifting to correct for the matrix columns zeroed by symmetricization 

166 if self.symmetric: 

167 # First, collect a vector containing all non-zero BCs. 

168 g = self.A.create_vec(dim=0) 

169 g.zero() 

170 for i,bcs in enumerate(self.bcs): 

171 for bc in bcs: 

172 bc.apply(g[i]) 

173 # The non-zeroes of g are now exactly the inhomogeneous BC 

174 # values. We can thus create the necessary modifications to b by 

175 # just multiplying with the un-symmetricised original matrix. The 

176 # bc values are overwritten below, hence only the non-BC rows of A 

177 # matter. 

178 b -= self.A * g 

179 

180 # Apply the actual BC dofs to b. (This must be done after the symmetric 

181 # correction above, since A*g is not necessarily zero at BC dofs.) 

182 # If the sign is negative, we negate twice to effectively apply -bc. 

183 for i,bcs in enumerate(self.bcs): 

184 if self.signs[i] != 1 and bcs: 184 ↛ 185line 184 didn't jump to line 185, because the condition on line 184 was never true

185 b[i] *= 1/self.signs[i] 

186 for bc in bcs: 

187 bc.apply(b[i]) 

188 if self.signs[i] != 1 and bcs: 188 ↛ 189line 188 didn't jump to line 189, because the condition on line 188 was never true

189 b[i] *= self.signs[i] 

190 

191 return self 

192 

193 def __call__(self, other): 

194 if not isinstance(other, block_vec): 

195 raise TypeError() 

196 b = other.copy() 

197 self.apply(b) 

198 return b