### Piecewise values for a function

156

views

0

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:

Thank you so much in advance

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

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:

The MeshFunction

```
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.

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

You could try out this approach and see if it suits your needs:

Would you explain the MeshFunction's topological dimension argument? Here

```
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.

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

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.

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.

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.

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.