Piecewise values for a function


156
views
0
10 weeks ago by
Nima  
Hi,

I am trying to change the values of a function corresponding to each sub-domain. So I created a toy problem. Suppose I have the following domain:


As you can see the marker on the left is '3' and the right is '4'.
I have the following code

from dolfin import *

mesh = Mesh('Square.xml')

Interior = MeshFunction('size_t' , mesh , 'Square_physical_region.xml' )  
boundary = MeshFunction('size_t', mesh , 'Square_facet_region.xml') 

V = FunctionSpace(mesh,'P',1)

f=TestFunction(V)

g=Function(V)

g_D = Constant(2)

a = dot(grad(g),grad(f))*dx+g*f*dx
L=Constant(0)

bc=DirichletBC(V,g_D,boundary)

solve(a == L, g, bc)

Now, I want the values of 'g' to be zero on the right domain (where Interior==4). Does anyone know a quick way to do this?

Thank you so much in advance
Community: FEniCS Project

It's very unclear to me what you're trying to do.

First, would you please provide an MWE that we can run by itself? For example, use fenics.UnitSquareMesh instead of loading an XML file. Also immediately after g=TrialFunction(V) you write g=Function(V), so there must be a mistake. It sort of looks to me like you're mixing up the notation for linear and nonlinear problems.

written 10 weeks ago by Alexander G. Zimmerman  
Thank you for the reply. About "g" it was a typo, so I edited my code. But the reason I did not use a fenics mesh was because I wanted to see if I can do the piecrwise values for the domain markers that come from GMSH.
written 10 weeks ago by Nima  
You could use this in place of the external mesh:
from dolfin import *

mesh = UnitSquareMesh(2,2)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= .5
        
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= .5
        
Interior = MeshFunction("size_t", mesh, mesh.topology().dim())

left = Left()
right = Right()
left.mark(meshFun, 3)
right.mark(meshFun, 4)​

The MeshFunction Interior you are going to use is essentially the same thing, no matter if your mesh comes from gmsh or whatever tool you want to use.
written 10 weeks ago by klunkean  
It seems I have done a good job confusing everyone. The domain marking is not my problem. I want to see if there is a quick way to assign for example 0 to "g" where "interior==3". I did not add the marking code because I have already done it via GMSH.

Thank you so much for the kind reply.
written 10 weeks ago by Nima  
We get that it's not the marking that is the problem. But we can't just copy your code and modify it because you did not supply us with your mesh xml file. That means it's not a minimal working example and that's why I added the snippet above.
written 10 weeks ago by klunkean  

1 Answer


5
10 weeks ago by
You could try out this approach and see if it suits your needs:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(9,9)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= .3
        
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= .3
        
meshFun = MeshFunction("size_t", mesh, mesh.topology().dim())

left = Left()
right = Right()
left.mark(meshFun, 3)
right.mark(meshFun, 4)

FE = FiniteElement('DG', mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, FE)

g = Function(V)

# first make g non-zero everywhere for demonstration purposes
g.interpolate(Constant(5.0)) 

subdmarker = np.where(meshFun.array()==4)
vals = g.vector().get_local()
vals[subdmarker] = 0
g.vector().set_local(vals)

plot(g)
plt.show()​
get_local() and set_local() might solve a related problem I've been digging into the past days. Thanks!

Would you explain the MeshFunction's topological dimension argument? Here len(meshFun) is 162, so I'm guessing there's a marker for each cell (because 162 = 2x9x9).
written 10 weeks ago by Alexander G. Zimmerman  
The third argument to the MeshFunction constructor refers to whether it is volumes (3), facets (2) or lines (1) you're marking.
In this case we mark facets (2D cells) so there is a value for each cell just as you pointed out.
written 10 weeks ago by klunkean  
Does this only work with DG0, since there is one degree of freedom per cell? I'm currently trying to do this with mixed higher order elements. Edit: Also OP has more than one DoF per cell.
written 10 weeks ago by Alexander G. Zimmerman  
I think this direct approach will only work with DG0. The function itself and the MeshFunction will have the same length and ordering so you can directly set the degrees of freedom corresponding to the entry in the MeshFunction.

I just thought of the following hack for more DoFs, nastily abusing the DirichletBC.apply() method:
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

mesh = UnitSquareMesh(9,9)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= .3
        
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= .3
        
meshFun = MeshFunction("size_t", mesh, mesh.topology().dim())

left = Left()
right = Right()
left.mark(meshFun, 3)
right.mark(meshFun, 4)

# More dofs for function u
FElin = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
Vlin = FunctionSpace(mesh, FElin)

u = Function(Vlin)
# non-zero values for demonstration
u.interpolate(Expression('x[0]', element=FElin))

hack = DirichletBC(Vlin, Constant(0.0), right)
hack.apply(u.vector())

plot(u)
plt.show()​

Notice, however, the possibly unwanted steep linear decrease over the boundary cells between the two subdomains which might result in errors when calculating integrals over these cells.
written 10 weeks ago by klunkean  
1
Brilliant. I won't blame you for the hack. It is very odd that FEniCS lets us apply "boundary conditions" on the interior of the domain, and this caused me some unexpected issues when first using FEniCS.

I asked a related question here that you might want to get credit for answering: https://www.allanswered.com/post/xzxlb/assign-mixed-element-function-values-on-subdomain/

I already tried this for the lid-driven cavity and it seems to work at first glance. I'll try for my real problem now.
written 10 weeks ago by Alexander G. Zimmerman  
Interesting. So no matter if you use cell by cell assignment or your method, we always see that weird error between the two subdomains. To be honest that was one of the reasons I asked this question, because I thought I was doing something wrong. But thank you so very much. It was enlightening in all sorts of ways.
written 10 weeks ago by Nima  
I actually don't understand why we would classify the behavior on the layer of cells between the domains as an error. For example with P1 elements, you're only setting degrees of freedom at vertices. Over a cell there will always be a linear interpolation of the vertex values. With DG0 you wouldn't have this behavior. This all seems as expected.
written 10 weeks ago by Alexander G. Zimmerman  
I totally agree with you and I wouldn't call it an error per se but it can lead to errors on a coarse mesh.
written 10 weeks ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »