How to compute the average value of u of a part of the domain and of a facet?

variational formulation
1

Hello everyone,

I try to solve a transient non-homogeneous poisson equation with Robin and Dirchlet boundary conditions.

I applied the source term f to a part of my 1x1 UnitSquareMesh domain by the following expression:

    class MyExpression(UserExpression):
        def eval(self, value, x):
            if x[0] > 0 and x[0] <= 0.1: 				# I
                alpha_t_plus_1_I = alpha_t_I - dt * sqrt(alpha_t_I)
                value[0] = alpha_t_plus_1_I * (constant - u_avg_I)         
            elif x[0] >= 0.9 and x[0] < 1.0:				# II
                alpha_t_plus_1_II = alpha_t_II - dt * sqrt(alpha_t_II)
                value[0] = alpha_t_plus_1_II * (constant - u_avg_II)
            elif x[1] > 0 and x[1] <= 0.1:				# III
                alpha_t_plus_1_III = alpha_t_III - dt * sqrt(alpha_t_III)
                value[0] = alpha_t_plus_1_III * (constant - u_avg_III)
            elif x[1] >= 0.9 and x[1] < 1.0:				# IV
                alpha_t_plus_1_IV = alpha_t_IV - dt * sqrt(alpha_t_IV)
                value[0] = alpha_t_plus_1_IV * (constant - u_avg_IV)
            #else:
                value[0] = 0
            
        def value_shape(self):
            return(1,)
                
    f = MyExpression()

, where u_{avg\,I/II/III/IV} is the average value of u for the previous time step in the region of the domain described by the if/elif condition. So for example, u_{avg\,I} is the average value of u for the part of the domain, that is described by x[0] > 0 and x[0] \leq 0.1.

So my first question is:
How do I compute u_{avg\,I}, u_{avg\,II}, u_{avg\,III} and u_{avg\,IV}?

I use the same set of formulars to describe the Dirchlet boundary condition for every facet:

    # define facets
    gamma00 = CompiledSubDomain("on_boundary")
    gamma01 = CompiledSubDomain("near(x[1], 0.0) && on_boundary")
    gamma02 = CompiledSubDomain("near(x[0], 1.0) && on_boundary")
    gamma03 = CompiledSubDomain("near(x[1], 1.0) && on_boundary")
    gamma04 = CompiledSubDomain("near(x[0], 0.0) && on_boundary")

    facet_marker = MeshFunction("size_t", mesh, space_dim - 1)
    facet_marker.set_all(0)
    gamma01.mark(facet_marker, 1)
    gamma02.mark(facet_marker, 2)
    gamma03.mark(facet_marker, 3)
    gamma04.mark(facet_marker, 4)

    alpha_t_plus_1_facet1 = alpha_t_facet1 - dt * sqrt(alpha_t­_facet1 )
    u_Dirichlet _facet1= alpha_t_plus_1_ facet1 * (constant – u_avg_facet1)
        
    alpha_t_plus_1_facet2 = alpha_t_ facet2 - dt * sqrt(alpha_t_facet2)
    u_Dirichlet _facet2= alpha_t_plus_1­_facet2 * (constant – u_avg_facet2)

    alpha_t_plus_1_facet3 = alpha_t_facet3 - dt * sqrt(alpha_t_facet3)
    u_Dirichlet _facet3= alpha_t_plus_1­_facet3 * (constant – u_avg_facet3)

    alpha_t_plus_1_facet4 = alpha_t_ facet4 - dt * sqrt(alpha_t_facet4)
    u_Dirichlet _facet4= alpha_t_plus_1_ facet4 * (constant - u_avg_facet4)

    dirichlet_bcs = []
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_facet1, facet_marker, 1))
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_facet2, facet_marker, 2))
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_facet3, facet_marker, 3))
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_facet4, facet_marker, 4))

Here, I‘ve got the same question:
How do i compute u_{avg\:facet1}, u_{avg\:facet2}, u_{avg\:facet3} and u_{avg\:facet4}?

Thanks,
Jan

P.S.: I have to apologize for the indexation, which is probably pretty confusing, if you don‘t work on that problem on a daily basis.

2

I am not entirely sure what you want to do, but if u is the function you want an average of on some mesh entities, you might want to calculate u_F = \frac{1}{|F|} \int_F u, which can be done e.g. for a boundary facet with a construction like:

u_F = assemble(u*ds(subdomain_data=mf_facets, subdomain_id=1)) / assemble(Constant(1.0)*ds(subdomain_data=mf_facets, subdomain_id=1))

where mf_facets is a MeshFunction object where facet F is marked with 1, similar to your second code example.

1 Like
3

Dear @volkerk,

I simplified my example, which I just made up for a code demonstration, rather than to be a reasonable model any physical system, as follows:

In order to solve a transient non-homogeneous poisson equation with Robin and Dirchlet boundary conditions, I applied the source term f to a part of my 1x1 UnitSquareMesh domain by the following expression:

    class MyExpression(UserExpression):
        def eval(self, value, x):
            if x[0] > 0 and x[0] <= 0.1: 				
                value[0] = Constant(2.0) * u_avg       
            else:
                value[0] = 0
            return value[0]
            
        def value_shape(self):
            return(1,)
                
    f = MyExpression()

, where u_{avg} is the average value of u for the previous time step in the region of the domain described by the if condition. So, u_{avg} is the average value of u for the part of the domain, that is described by x[0]>0 and x[0]≤0.1 .

With your help,

I was able to compute u_{avg} by the following code:

u_avg = assemble(u*dx(1)) / assemble(Constant(1.0)*dx(1))

, whereas dx(1) was implemented in the following manner

# Create classes for defining parts of the interior of the domain   
 class Obstacle_left(SubDomain):
        def inside(self, x, on_boundary):
            return (between(x[0], (0.0, 0.1)))

# Initialize sub-domain instances
obstacle_left = Obstacle_left()

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, space_dim-1)
domains.set_all(0)
obstacle_left.mark(domains, 1)

# Define new measures associated with the interior domains
dx = Measure('dx', domain=mesh, subdomain_data=domains)

Thanks a lot!