Project marked mesh

mesh
1

Hi! I have a gmsh mesh and the separate array of markers for every point of this mesh. I want to create UserExpression and project this expression on the function space. The problem is next: how can I get the index of the current evaluated point in the markers array? I want to do something like this (simplified example):

class MyExpression(UserExpression):

    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
   
    # In reality we don't have idx variable
    def eval(self, values, x, idx):  
        point_marker = self.markers[idx]

        if point_marker == 0:
            values[0] = x[0] ** 2
        else:
            values[0] = 0
            
    def value_shape(self):
        return ()

mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'Lagrange', 1)
markers = [1, 0, 1, 0, 1, 0, 1, 0]
expr = MyExpression(markers)
result = project(expr, V)

Is it possible? Thank you in advance!

2

To me it seems easier to first define your custom UserExpression after which you apply the markers to select an appropriate part of your solution. You can use the FEniCS MeshFunction to define your markers, or use your own array. Try something like this:

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)

# Interpolate expression
class MyExpression(UserExpression):
    def eval(self, values, x):
        values[0] = x[0] ** 2
    def value_shape(self):
        return ()
expr = MyExpression(degree=1)
u = project(expr, V)

# Define markers
vertex_markers = MeshFunction("size_t", mesh, dim=0)
vertex_markers.set_all(0)
for vertex in vertices(mesh):
    if vertex.x(0) <= 0.5: # Example
        vertex_markers.set_value(vertex.index(), 1)

# Apply on expr
v2d = vertex_to_dof_map(V)
u.vector()[v2d] = u.vector()[v2d] * vertex_markers.array()
1 Like
3

Thank you very much!

4

Hi! Have you any ideas about how can I project function with multiple conditions (not only two: 1 or 0)? For example:

if (x[0] < 0.25):
    values[0] = x[0] ** 2
elif (0.25 <= x[0] < 0.5):
    values[0] = x[0] ** 3
else:
    values[0] = x[0] ** 4
5

I am not sure what exactly you are trying to do here. In case you would like a UserExpression with multiple conditions, define something as follows:

class MyExpression(UserExpression):
    def eval(self, values, x):
        if (x[0] < 0.25):
            values[0] = x[0] ** 2
        elif (0.25 <= x[0] < 0.5):
            values[0] = x[0] ** 3
        else:
            values[0] = x[0] ** 4
    def value_shape(self):
        return ()

In case you have more complex divisions of your domain, consider using the MeshFunction functionality in combination with the definition of SubDomains as in the following demo: https://fenicsproject.org/docs/dolfin/2018.1.0/python/demos/subdomains/documentation.html?highlight=subdomain. I hope this answers your question.

1 Like