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

1from __future__ import division 

2 

3from builtins import str 

4from builtins import range 

5from builtins import object 

6"""Utility functions for plotting, boundaries, etc.""" 

7 

8import os 

9from dolfin import * 

10import time as timer 

11 

12class BoxBoundary(object): 

13 def __init__(self, mesh): 

14 c = mesh.coordinates() 

15 self.c_min, self.c_max = c.min(0), c.max(0) 

16 dim = len(self.c_min) 

17 

18 sd = self._compile(west = self._boundary(0, self.c_min) if dim>1 else '0', 

19 east = self._boundary(0, self.c_max) if dim>1 else '0', 

20 south = self._boundary(1, self.c_min) if dim>2 else '0', 

21 north = self._boundary(1, self.c_max) if dim>2 else '0', 

22 bottom= self._boundary(dim-1, self.c_min), 

23 top = self._boundary(dim-1, self.c_max), 

24 ew = self._boundary(0) if dim>1 else '0', 

25 ns = self._boundary(1) if dim>1 else '0', 

26 tb = self._boundary(dim-1), 

27 all = 'on_boundary') 

28 for name,subdomain in sd: 

29 setattr(self, name, subdomain) 

30 

31 def _boundary(self, idx, coords=None): 

32 if coords is not None: 

33 return 'on_boundary && near(x[{idx}], {coord})' \ 

34 .format(idx=idx, coord=coords[idx]) 

35 else: 

36 return 'on_boundary && (near(x[{idx}], {min}) || near(x[{idx}], {max}))' \ 

37 .format(idx=idx, min=self.c_min[idx], max=self.c_max[idx]) 

38 

39 def _compile(self, **kwargs): 

40 # Make sure all expressions sent to compile_subdomains are different 

41 expr_to_code = {} 

42 for expr in list(kwargs.values()): 

43 expr_to_code[expr] = None 

44 

45 #print expr_to_code.keys() 

46 #compiled = CompiledSubDomain(expr_to_code.keys()) 

47 for i, expr in enumerate(expr_to_code.keys()): 

48 expr_to_code[expr] = CompiledSubDomain(expr) 

49 

50 return [(name, expr_to_code[expr]) for name, expr in list(kwargs.items())] 

51 

52class update(object): 

53 """Plot and save given functional(s). Example: 

54 u = problem.solve() 

55 update.set_args(displacement={'mode': 'displacement'}) 

56 update(displacement=u, volumetric=tr(sigma(u))) 

57 """ 

58 files = {} 

59 plots = {} 

60 kwargs = {} 

61 projectors = {} 

62 functions = {} 

63 

64 def _extract_function_space(self, expression, mesh=None): 

65 """Try to extract a suitable function space for projection of 

66 given expression. Copied from dolfin/fem/projection.py""" 

67 import ufl 

68 

69 # Extract functions 

70 functions = ufl.algorithms.extract_coefficients(expression) 

71 

72 # Extract mesh from functions 

73 if mesh is None: 

74 for f in functions: 

75 if isinstance(f, Function): 

76 mesh = f.function_space().mesh() 

77 if mesh is not None: 

78 break 

79 if mesh is None: 

80 raise RuntimeError 

81 

82 # Create function space 

83 shape = expression.shape() 

84 if shape == (): 

85 V = FunctionSpace(mesh, "CG", 1) 

86 elif len(shape) == 1: 

87 V = VectorFunctionSpace(mesh, "CG", 1, dim=shape[0]) 

88 elif len(shape) == 2: 

89 V = TensorFunctionSpace(mesh, "CG", 1, shape=shape) 

90 else: 

91 raise RuntimeError 

92 

93 return V 

94 

95 def project(self, f, name, V, mesh=None): 

96 if V is None: 

97 # If trying to project an Expression 

98 if isinstance(f, Expression): 

99 if isinstance(mesh, cpp.Mesh): 

100 V = FunctionSpaceBase(mesh, v.ufl_element()) 

101 else: 

102 raise TypeError 

103 else: 

104 V = self._extract_function_space(f, mesh) 

105 key = str(V) 

106 v = TestFunction(V) 

107 if not key in self.projectors: 

108 # Create mass matrix 

109 u = TrialFunction(V) 

110 a = inner(v,u) * dx 

111 solver = LinearSolver("direct") 

112 solver.set_operator(assemble(a)) 

113 #solver.parameters['preconditioner']['reuse'] = True 

114 self.projectors[key] = solver 

115 # Use separate function objects for separate quantities, since this 

116 # object is used as key by viper 

117 if not name in self.functions: 

118 self.functions[name] = Function(V) 

119 

120 solver, Pf = self.projectors[key], self.functions[name] 

121 b = assemble(inner(v,f) * dx) 

122 

123 # Solve linear system for projection 

124 solver.solve(Pf.vector(), b) 

125 

126 return Pf 

127 

128 

129 def set_args(self, **kwargs): 

130 """Set additional kwargs to pass to plot for a given name. 

131 

132 In addition to the kwargs for plot, these are accepted: 

133 'plot' (bool) -- plot to screen [True] 

134 'save' (bool) -- save to file [True] 

135 'functionspace' (FunctionSpace) -- space to project to [CG(1)]""" 

136 self.kwargs.update(kwargs) 

137 

138 def save_to_file(self, name, data, time): 

139 if not os.path.exists('data'): 

140 os.mkdir('data') 

141 if not name in self.files: 

142 self.files[name] = File('data/%s.pvd'%name) 

143 if time is not None: 

144 self.files[name] << (data, time) 

145 else: 

146 self.files[name] << data 

147 

148 def plot(self, name, title, data, time): 

149 kwargs = self.kwargs.get(name, {}) 

150 if not name in self.plots: 

151 self.plots[name] = plot(data, title=title, size=(400,400), 

152 axes=True, warpscalar=False, 

153 **kwargs) 

154 else: 

155 self.plots[name].update(data, title=title, **kwargs) 

156 

157 def __call__(self, time=None, postfix="", **functionals): 

158 for name,func in sorted(functionals.items()): 

159 args = self.kwargs.get(name, {}) 

160 if 'functionspace' in args or not isinstance(func, Function): 

161 func = self.project(func, name, args.get('functionspace')) 

162 if hasattr(func, 'rename'): 

163 func.rename(name+postfix, name+postfix) 

164 if args.get('plot', True): 

165 self.plot(name, name+postfix, func, time) 

166 if args.get('save', True): 

167 self.save_to_file(name+postfix, func, time) 

168 

169update = update() # singleton 

170 

171def orthogonalize(v, basis): 

172 """basis vectors are assumed to be normalized -- w.inner(w)==1.0""" 

173 for w in basis: 

174 v -= w.inner(v)*w 

175 

176def rigid_body_modes(V, show_plot=False): 

177 """Compute orthogonal rigid body modes of a function space.""" 

178 T = timer.time() 

179 

180 mesh = V.mesh() 

181 if mesh.geometry().dim() == 2: 181 ↛ 184line 181 didn't jump to line 184, because the condition on line 181 was never false

182 modes = [Constant((1, 0)), Constant((0, 1)), Expression(('x[1]', '-x[0]'), degree=1)] 

183 else: 

184 translations = [Constant((1, 0, 0)), Constant((0, 1, 0)), Constant((0, 0, 1))] 

185 rotations = [Expression(('0', 'x[2]', '-x[1]'), degree=1), 

186 Expression(('-x[2]', '0', 'x[0]'), degree=1), 

187 Expression(('x[1]', '-x[0]', '0'), degree=1)] 

188 

189 modes = translations + rotations 

190 

191 modes = [interpolate(m, V).vector() for m in modes] 

192 basis = VectorSpaceBasis(modes) 

193 basis.orthonormalize() 

194 

195 if show_plot: 195 ↛ 196line 195 didn't jump to line 196, because the condition on line 195 was never true

196 import matplotlib.pyplot as plt 

197 for mode in modes: 

198 plt.figure() 

199 plot(Function(V,mode)) 

200 plt.show() 

201 

202 info("computed rigid body modes in %.2f s"%(timer.time()-T)) 

203 return modes