Integral of solution on (sub)domains

[fewagner]

Hi all,

Great project! I only recently started using skfem and am still getting familiar with the syntax. So here is a basic question:

I am working with 2D triangular meshes, and P1 basis functions. I would like to integrate my solution over

1. the whole mesh, or
2. a 2D subdomain, or
3. along a line.

For 3. I could just use a probe, and integrate with numpy.trapz. Is this the most elegant solution, or are there any built-in alternatives?

For 2. - is it possoble to define a 2D probe? Or what would be the best way to do this?

For 1., this should not be too tricky, I would just need to build the scalar product of the solution vector with the weight of each triangle? Is it possible to access these weights; or should I calculate them from the "dx" property of the mesh?

Thank you already!

[kinnala]

The basic idea is described here: https://scikit-fem.readthedocs.io/en/latest/howto.html#postprocessing-the-solution

To integrate over subdomain, you tag the respective elements using Mesh.with_subdomains or using external mesh generator and initialize intbasis=Basis(mesh, elem, elements='yourtag') to use together with Functional wrapper.

To integrate over facets you use FacetBasis. I am now on mobile phone so cannot answer in more detail.

[fewagner]

Hi, thanks for the quick response, I think that solves it! I'll come back if I don't manage from here.

[fewagner]

Hi, sorry, my original question was answered but I haven't fully managed to do what I wanted to, it's possible that I am just confused. I would appreciate any comments.

I am solving a standard diffusion equation with trivial Dirichlet boundary conditions on the left side of a quadratic mesh; with an implicit Euler algorithm in time:

mesh = fem.MeshTri().refined(6)  
basis = fem.Basis(mesh, fem.ElementTriP1())
fbasis = basis.boundary('left')
boundary = basis.get_dofs('left').all()

def LHS(u, v, kwargs): 
    out = u * v 
    out += kwargs['k'] * kwargs['D'] * dot(grad(u), grad(v)) 
    return out 

LHS_form = fem.BilinearForm(LHS)

def RHS(v, kwargs):
    out = kwargs['u_prev_t'] * v
    return out 

RHS_form = fem.LinearForm(RHS)

As starting values, I define a region in the center of the mesh thats filled with ones (to approximate a Dirac delta function). Then I iterate over the times, and find the FEM solution for each time step:

times = np.linspace(0,1,200)

def starting_values(x):
    rad = np.sqrt((x[0]-0.5)**2 + (x[1]-0.5)**2)
    return np.heaviside(0.05-rad, 0.)

u = basis.project(starting_values)
solutions = [u.copy()]

# define D
D = .1

# ITERATE TIMES 
for i in trange(1,len(times)):

    k = times[i] - times[i-1]

    RHS_vec = RHS_form.assemble(basis, u_prev_t=u)  
    LHS_matrix = LHS_form.assemble(basis, k=k, D=D)
    u = fem.solve(*fem.condense(LHS_matrix, RHS_vec, D=boundary)) 
    
    solutions.append(u.copy())

Now in the postprocessing of the solution, I want to track two things over time: the total density within the mesh (called norm in the code below) and the flux across the Dirichlet boundary (i.e. the density lost to the environment, called flux below). I am calculating these quantities as follows:

@fem.Functional
def integral(kwargs):
    return kwargs['uh'] 

norm = [integral.assemble(basis, uh=basis.interpolate(s)) for s in solutions] 

@fem.Functional
def bndintegral(kwargs):
    return kwargs['uh'].grad[0] 

flux = [bndintegral.assemble(fbasis, uh=fbasis.interpolate(s)) for s in solutions]

The last part is mostly copied from the tutorial. Now the part that I don't understand is, I would expect that the time-integrated flux over the boundary should make up for the decrease of the density within the mesh. But thats not the case here, in fact the flux over the boundary is much larger than the density in the mesh (see plot below). I think I have some misconception about how the facet basis and the integral over that works?

I would be very happy for suggestions what I am doing wrong.

[kinnala]

Perhaps you want to multiply the derivative by the diffusion coefficient to get the heat flux?

[fewagner]

Aaah, yes of course, because I was tracking the space-derivative, but was interested in the time-derivative. Thanks so much! :)

I actually have a second question, that is related, and for simplicity I would append it here:

If I understand correctly, I enforce trivial Neumann boundary conditions (du/dn=0) along all boundaries where I don't explicitly enforce Dirichlet boundary conditions.

If I wanted to enforce along one of the edges non-trivial Neumann boundary conditions, e.g. du/dn = 1, how would that work? Is there an as simple way to do this as for the Dirichlet conditions, or is it more complicated?

[kinnala]

For variational formulations, boundary conditions can be split into two categories: (i) essential and (ii) natural. In this case Dirichlet conditions are category (i) and those types of conditions are typically implemented at linear algebra level by constraining some DOFs to zero or other value. The Neumann conditions go to category (ii) meaning that they will be automatically satisfied by the variational formulation if correct terms have been added to the bilinear and linear forms. In this case another term must be added to the linear form which is (1, v) integrated over the respective boundary.

[fewagner]

Thank you again! You were a great help :)