dolfin team mailing list archive

[mspieg <mspieg@xxxxxxxxxxxxxxxxx>]

"""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 (mspieg@xxxxxxxxxxxxxxxxx)"
__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:])