Why I can't assemble a matrix on a submesh?

Question

I encounter the problem in domain decomposition,
For Example:

from dolfin import *

mesh = UnitSquareMesh(4,4)
V = FunctionSpace(mesh,'CG',1)
u = TrialFunction(V)
v = TestFunction(V)

class My_domain(SubDomain):
    def inside(self,x,on_boundary):
        return x[0]<=0.5

sub_domain = CellFunction('size_t',mesh)
sub_domain.set_all(0)
my_domain = My_domain()
my_domain.mark( sub_domain,1 )

mesh1 = SubMesh( mesh,sub_domain,1 )

a = u*v*dx
A = assemble(a,mesh = mesh1)

Thanks !

Answer 1

Hi, consider

from dolfin import *

mesh = UnitSquareMesh(4,4)

class My_domain(SubDomain):
    def inside(self,x,on_boundary):
        return x[0]<=0.5

subdomain = CellFunction('size_t',mesh)
subdomain.set_all(0)
my_domain = My_domain()
my_domain.mark(subdomain,1 )

# dim V x dim V mass matrix with zero rows for basis functions
# with no support in subdomain
V = FunctionSpace(mesh,'CG',1)
u = TrialFunction(V)
v = TestFunction(V)

dx_sub = Measure('dx', domain=mesh, subdomain_id=1, subdomain_data=subdomain)
a = u*v*dx_sub
A = assemble(a)

A_ = A.array()
n_basis_in_sub = sum(1 for row in A_ if np.abs(row).sum() < DOLFIN_EPS)
print A_.shape, n_basis_in_sub

# mass matrix of space defined over subdomain
submesh = SubMesh(mesh, subdomain, 1)
V = FunctionSpace(submesh,'CG',1)
u = TrialFunction(V)
v = TestFunction(V)

a = u*v*dx
A = assemble(a)

A_ = A.array()
print A_.shape 

Comments

Would you know of a good way to determine which global dofs correspond to the dofs of the subproblem?

---

Are you asking for some map between dofs of the full problem that are "in the subdomain" and dofs of the problem defined only in the subdomain?

---

Yes, exactly. For instance for dd methods one would setup a global system from the assembled subdomain matrices and introduce additional interface couplings (e.g. Lagrange multipliers). In order to obtain the global solution, it is required to know which subdomain dofs correspond to which global dofs.
How can one obtain this mapping?

---

Okay, this is a nice topic so in order for it to get the attention it deserves I suggest you open it as a new/separate question.

---

Sure, can do. However, you already indicate by your suggestion that there is no obvious or simple way to get this information (which really is a pity). Otherwise you'd certainly just written it down ;)

---

You got me :). But seriously, separate question will attract more answers and once answered it will be easier to find for others with a similar problem.

---

Thanks for your answer, but with introduce a new function space on the sub mesh, I will can't know the relationship between local dofs and global ones !