|
Hi All, I think there is a bug when trying to read a MeshFunction from a file that was generated on a SubMesh. I've attached some python based on the SubMesh demo, that demonstrates the problem. I think the issue is that there is an extra MeshFunction embedded in a SubMesh to describe the "global vertex indices" that gets read instead of the intended mesh-function. On a related note, it's easy enough to extract a SubMesh from a larger Mesh using a MeshFunction, but what is the best way to extract just the piece of a discrete function on the SubMesh, when the Function is originally defined for a Function Space on the larger mesh. (does this make sense?) These issues arise from a multiphysics problem that I'm working with overlapping domains where I basically want to solve an embedded porous media problem on a subdomain of a larger stokes flow (this is slightly different from the fluid-structure demo where the two domains only share a boundary.) All help greatly appreciated as always cheers marc p.s. all current dev repositories p.p.s. another official release would be great...a lot has changed since 2009-04 |
"""demonstration of error reading a MeshFunction defined on a submesh from a file demo is based on dolfin demo/mesh/submesh/python/demo.py """
__author__ = "Marc Spiegelman (msp...@ldeo.columbia.edu)"
__date__ = "26 Jun 2009 14:38:56"
__copyright__ = "Copyright (C) 2009 Marc Spiegelman"
__license__ = "GNU LGPL Version 2.1"
from dolfin import *
from numpy import *
import sys
# Structure sub domain
class Structure(SubDomain):
def inside(self, x, on_boundary):
b= array([1.4, 1.6, 0., 0.6 ])
return x[0] > b[0] - DOLFIN_EPS and x[0] < b[1] + DOLFIN_EPS and x[1] < b[3] + DOLFIN_EPS
# Boundary of structure mesh
class StructureBoundary(SubDomain):
def inside(self,x,on_boundary):
return on_boundary
def main(args):
# Create mesh
mesh = Rectangle(0.0, 0.0, 3.0, 1.0, 60, 20)
# Create sub domain markers and mark everaything as 0
sub_domains = MeshFunction("uint", mesh, mesh.topology().dim())
sub_domains.set_all(0)
# Mark structure domain as 1
structure = Structure()
structure.mark(sub_domains, 1)
# Extract sub mesh
structure_mesh = SubMesh(mesh, sub_domains, 1)
# mark the boundary facets on the structure mesh
sub_domain_boundary = MeshFunction("uint",structure_mesh,structure_mesh.topology().dim()-1)
sub_domain_boundary.set_all(0)
structure_boundary = StructureBoundary()
structure_boundary.mark(sub_domain_boundary,1)
# so far so good...the sub_domain_boundary mesh function is correct
print 'sub_domain_boundary.values()'
print sub_domain_boundary.values()
# save sub_domains to a file and read back in
f = File("subdomain_boundaries.xml")
f << sub_domain_boundary
sub_domain_boundary_from_file = MeshFunction("uint",structure_mesh,"subdomain_boundaries.xml")
# this is broken...the mesh function seems to actually be the global vertex indices of the submesh
print '\nsub_domain_boundary_from_file.values()'
print sub_domain_boundary_from_file.values()
if __name__ == "__main__":
main(sys.argv[1:])
---------------------------------------------------- Marc Spiegelman Lamont-Doherty Earth Observatory Dept. of Applied Physics/Applied Math Columbia University tel: 845 704 2323 (SkypeIn) ---------------------------------------------------- |
_______________________________________________ DOLFIN-dev mailing list DOLFIN-dev@fenics.org http://www.fenics.org/mailman/listinfo/dolfin-dev
