subdomain marking problem

Bug #1206559 reported by Peter Kekenes-Huskey
6
This bug affects 1 person
Affects Status Importance Assigned to Milestone
DOLFIN
New
Undecided
Unassigned

Bug Description

Hi,
I tried to apply the protocol for marking cells using the based on the poisson problem with two materials example (http://fenicsproject.org/pub/book/tutorial/stationary/poisson/mat2_p2D.py) for a 2D mesh I created with gmsh (attached). It appears that the classes for marking the subdomains do not work for my mesh (several elements are missed); unfortunately it is not obvious to me why this is failing. Is there something incorrect about the mesh I am using (it runs fine with other dolfin PDE solves) or the code I have below?

Thanks for your help
pete

from dolfin import *
import sys, math, numpy
import numpy as np

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] <= 0.5 else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] >= 0.5 else False

if(1):
    nx = 4; ny = 6
    mesh = UnitSquare(nx, ny)
    #mesh = Mesh('nohole.xml')
    mesh = Mesh('hole.xml')

    # Define a MeshFunction over two subdomains
    subdomains = MeshFunction('uint', mesh, 2)

    # Mark subdomains with numbers 0 and 1
    subdomain0 = Omega0()
    subdomain0.mark(subdomains, 0)
    subdomain1 = Omega1()
    subdomain1.mark(subdomains, 1)

    # Discontinuous basis is used for boundary condition
    # LAter, continuous basis used for PDE solution
    V0 = FunctionSpace(mesh, 'DG', 0)
    k = Function(V0)

    print 'mesh:', mesh
    print 'subdomains:', subdomains
    print 'k:', k

    import numpy
    help = numpy.asarray(subdomains.array(), dtype=numpy.int32)
    #print help
    print "bad ", np.where(help>1)
    k.vector()[:] = numpy.choose(help, k_values)

Revision history for this message
Peter Kekenes-Huskey (huskeypm) wrote :
description: updated
To post a comment you must log in.
This report contains Public information  
Everyone can see this information.

Other bug subscribers

Remote bug watches

Bug watches keep track of this bug in other bug trackers.