L-shape mesh from gmsh and boundary conditions
Dear All, I have the following mesh from gmsh (see attached). I can read it just fine: mesh = Mesh("lshape.msh") but when I try to solve equation on this mesh it fails to apply boundary conditions (because boundaries isn't marked in the msh file, I guess): --------------------------------------------------------------------------- LookupError Traceback (most recent call last) <ipython-input-17-ae4ad3326f99> in <module>() 6 fun0 = None 7 cff = 0.002 ----> 8 sol = get_sol(frm, rhs) 9 fun0 = get_functional(sol, f) <ipython-input-15-e7cda2e05bae> in get_sol(frm, rhs) 4 sol = Function(V) 5 #solve(frm, sol, rhs, solver_parameters={'ksp_type': 'cg'}) ----> 6 solve(frm, sol, rhs, solver_parameters={"ksp_type": "preonly", "pc_type": "lu"}) 7 8 return sol /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/solving.pyc in solve(*args, **kwargs) 121 else: 122 # Solve pre-assembled system --> 123 return _la_solve(*args, **kwargs) 124 125 /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/solving.pyc in _la_solve(A, x, b, **kwargs) 209 solver = ls.LinearSolver(A, solver_parameters=solver_parameters, 210 nullspace=nullspace, --> 211 options_prefix=options_prefix) 212 213 solver.solve(x, b) /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/linear_solver.pyc in __init__(self, A, P, solver_parameters, nullspace, options_prefix) 87 # applied 88 # Force evaluation here ---> 89 self.ksp.setOperators(A=self.A.M.handle, P=self.P.M.handle) 90 91 @cached_property /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/matrix.pyc in M(self) 143 just need a handle on the :class:`pyop2.Mat` object it's 144 wrapping, use :attr:`_M` instead.""" --> 145 self.assemble() 146 # User wants to see it, so force the evaluation. 147 self._M._force_evaluation() <decorator-gen-414> in assemble(self) /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/utils.pyc in wrapper(f, *args, **kwargs) 60 opts["type_check"] = safe 61 try: ---> 62 return f(*args, **kwargs) 63 finally: 64 opts["type_check"] = check /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/matrix.pyc in assemble(self) 72 return 73 self._bcs_at_point_of_assembly = copy.copy(self.bcs) ---> 74 self._assembly_callback(self.bcs) 75 self._assembled = True 76 /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assembly_cache.pyc in inner(bcs) 356 return r 357 --> 358 r = thunk(bcs) 359 if isinstance(r, float): 360 # 0-form case /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in thunk(bcs) 328 tensor_arg = mat(lambda s: get_map(s, tsbc, decoration), 329 lambda s: get_map(s, trbc, decoration), --> 330 i, j) 331 elif is_vec: 332 tensor_arg = vec(lambda s: get_map(s), i) /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in mat(testmap, trialmap, i, j) 189 def mat(testmap, trialmap, i, j): 190 return tensor[i, j](op2.INC, --> 191 (testmap(test.function_space()[i])[op2.i[0]], 192 trialmap(trial.function_space()[j])[op2.i[1]]), 193 flatten=True) /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in <lambda>(s) 326 # Output argument 327 if is_mat: --> 328 tensor_arg = mat(lambda s: get_map(s, tsbc, decoration), 329 lambda s: get_map(s, trbc, decoration), 330 i, j) /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in get_map(x, bcs, decoration) 280 281 def get_map(x, bcs=None, decoration=None): --> 282 return x.cell_node_map(bcs) 283 284 elif integral_type in ("exterior_facet", "exterior_facet_vert"): /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/functionspaceimpl.pyc in cell_node_map(self, bcs) 411 "cell_node", 412 self.offset, --> 413 parent) 414 415 def interior_facet_node_map(self, bcs=None): /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/functionspacedata.pyc in get_map(self, V, entity_set, map_arity, bcs, name, offset, parent) 485 # Cache miss. 486 # Any top and bottom bcs (for the extruded case) are handled elsewhere. --> 487 nodes = [bc.nodes for bc in lbcs if bc.sub_domain not in ['top', 'bottom']] 488 decorate = False 489 if nodes: /home/george/polygon/firedrake/local/lib/python2.7/site-packages/pyop2/utils.pyc in __get__(self, obj, cls) 62 if obj is None: 63 return self ---> 64 obj.__dict__[self.__name__] = result = self.fget(obj) 65 return result 66 /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/bcs.pyc in nodes(self) 152 fs.exterior_facet_boundary_node_map( 153 self.method).values_with_halo.take( --> 154 fs._mesh.exterior_facets.subset(self.sub_domain).indices, 155 axis=0)) 156 /home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/mesh.pyc in subset(self, markers) 125 raise LookupError( 126 '{0} is not a valid marker'. --> 127 format(marker)) 128 129 # build a list of indices corresponding to the subsets selected by LookupError: 4 is not a valid marker So, I need to somehow mark boundaries, but I cannot find it out, how. I've tried https://openfoamwiki.net/index.php/2D_Mesh_Tutorial_using_GMSH but step 4 does not work for me -- I simply cannot select line with gui, and it is not clear how to select boundaries editing mesh file. Sincerely, George
Hi George, You probably need to mark the boundaries with integers in your GMSH geo file like this: Physical Line(1) = {2}; Physical Line(2) = {4}; Physical Line(3) = {1,3}; You can then use those integers {1,2,3} to identify boundaries in firedrake. See the attached file for a simple example. - Tuomas On 03/31/2016 01:00 PM, George Ovchinnikov wrote:
Dear All,
I have the following mesh from gmsh (see attached).
I can read it just fine: mesh = Mesh("lshape.msh")
but when I try to solve equation on this mesh it fails to apply boundary conditions (because boundaries isn't marked in the msh file, I guess):
--------------------------------------------------------------------------- LookupError Traceback (most recent call last) <ipython-input-17-ae4ad3326f99> in <module>() 6 fun0 = None 7 cff = 0.002 ----> 8 sol = get_sol(frm, rhs) 9 fun0 = get_functional(sol, f)
<ipython-input-15-e7cda2e05bae> in get_sol(frm, rhs) 4 sol = Function(V) 5 #solve(frm, sol, rhs, solver_parameters={'ksp_type': 'cg'}) ----> 6 solve(frm, sol, rhs, solver_parameters={"ksp_type": "preonly", "pc_type": "lu"}) 7 8 return sol
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/solving.pyc in solve(*args, **kwargs) 121 else: 122 # Solve pre-assembled system --> 123 return _la_solve(*args, **kwargs) 124 125
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/solving.pyc in _la_solve(A, x, b, **kwargs) 209 solver = ls.LinearSolver(A, solver_parameters=solver_parameters, 210 nullspace=nullspace, --> 211 options_prefix=options_prefix) 212 213 solver.solve(x, b)
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/linear_solver.pyc in __init__(self, A, P, solver_parameters, nullspace, options_prefix) 87 # applied 88 # Force evaluation here ---> 89 self.ksp.setOperators(A=self.A.M.handle, P=self.P.M.handle) 90 91 @cached_property
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/matrix.pyc in M(self) 143 just need a handle on the :class:`pyop2.Mat` object it's 144 wrapping, use :attr:`_M` instead.""" --> 145 self.assemble() 146 # User wants to see it, so force the evaluation. 147 self._M._force_evaluation()
<decorator-gen-414> in assemble(self)
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/utils.pyc in wrapper(f, *args, **kwargs) 60 opts["type_check"] = safe 61 try: ---> 62 return f(*args, **kwargs) 63 finally: 64 opts["type_check"] = check
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/matrix.pyc in assemble(self) 72 return 73 self._bcs_at_point_of_assembly = copy.copy(self.bcs) ---> 74 self._assembly_callback(self.bcs) 75 self._assembled = True 76
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assembly_cache.pyc in inner(bcs) 356 return r 357 --> 358 r = thunk(bcs) 359 if isinstance(r, float): 360 # 0-form case
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in thunk(bcs) 328 tensor_arg = mat(lambda s: get_map(s, tsbc, decoration), 329 lambda s: get_map(s, trbc, decoration), --> 330 i, j) 331 elif is_vec: 332 tensor_arg = vec(lambda s: get_map(s), i)
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in mat(testmap, trialmap, i, j) 189 def mat(testmap, trialmap, i, j): 190 return tensor[i, j](op2.INC, --> 191 (testmap(test.function_space()[i])[op2.i[0]], 192 trialmap(trial.function_space()[j])[op2.i[1]]), 193 flatten=True)
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in <lambda>(s) 326 # Output argument 327 if is_mat: --> 328 tensor_arg = mat(lambda s: get_map(s, tsbc, decoration), 329 lambda s: get_map(s, trbc, decoration), 330 i, j)
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/assemble.pyc in get_map(x, bcs, decoration) 280 281 def get_map(x, bcs=None, decoration=None): --> 282 return x.cell_node_map(bcs) 283 284 elif integral_type in ("exterior_facet", "exterior_facet_vert"):
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/functionspaceimpl.pyc in cell_node_map(self, bcs) 411 "cell_node", 412 self.offset, --> 413 parent) 414 415 def interior_facet_node_map(self, bcs=None):
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/functionspacedata.pyc in get_map(self, V, entity_set, map_arity, bcs, name, offset, parent) 485 # Cache miss. 486 # Any top and bottom bcs (for the extruded case) are handled elsewhere. --> 487 nodes = [bc.nodes for bc in lbcs if bc.sub_domain not in ['top', 'bottom']] 488 decorate = False 489 if nodes:
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/pyop2/utils.pyc in __get__(self, obj, cls) 62 if obj is None: 63 return self ---> 64 obj.__dict__[self.__name__] = result = self.fget(obj) 65 return result 66
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/bcs.pyc in nodes(self) 152 fs.exterior_facet_boundary_node_map( 153 self.method).values_with_halo.take( --> 154 fs._mesh.exterior_facets.subset(self.sub_domain).indices, 155 axis=0)) 156
/home/george/polygon/firedrake/local/lib/python2.7/site-packages/firedrake/mesh.pyc in subset(self, markers) 125 raise LookupError( 126 '{0} is not a valid marker'. --> 127 format(marker)) 128 129 # build a list of indices corresponding to the subsets selected by
LookupError: 4 is not a valid marker
So, I need to somehow mark boundaries, but I cannot find it out, how. I've tried https://openfoamwiki.net/index.php/2D_Mesh_Tutorial_using_GMSH but step 4 does not work for me -- I simply cannot select line with gui, and it is not clear how to select boundaries editing mesh file.
Sincerely, George
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
- 
                
                George Ovchinnikov
- 
                
                Tuomas Karna