Integrating FEniCS solution up to some point

variational formulation
1

As a part of a bigger problem there is a first order equation u'(x) = f(x) on [0,1] coupled with integration condition \int_{[0,1]}u(x)dx = 1 which has to be solved in every step. RHS f is a solution of some other equation. This equation can be solved by pure integration but question is how to do that in FEniCS.

As far as I know we can integrate on a whole domain or on subdomains but how can we say something like u(x)=\int_0^xf(s)ds?

Pick points from mesh, define step functions and then integrate product f\cdot\text{step} and later interpolate to starting function space?

Do you maybe know some tutorial about similar problem?

2

Did you find a solution for your problem? I’m facing a similar situation.

3

I would suggest implementing your own integration scheme with quadrature points. As you can evaluate a function at any point, i.e. u(x_p), where x_p is the point of interest, you have that \int_{0}^{x}u(x)\mathrm{ds}=\sum_q w_q u(p_q) where p_q are the quadrature points and w_q the weights on the interval [0, x]

4

Thank you for the reply. So, considering your comment, an attempt to include the term f(x)\cdot v \Bigg|^{x_1}_{x_0} in the variational form of the problem can be done by evaluating:

f(x)\cdot v \Bigg|^{x_1}_{x_0} = f(x_1)*v(x_1) - f(x_0)*v(x_0)

In that case, how could one access the values of the weights in x = x_1 and x = x_0?

5

You can find quadrature points and weights for the interval [0,1] at wikipedia: Gaussian quadrature - Wikipedia

You could also, as suggested in the original post, create a step function, (for instance using ufl.conditional), as shown in What is the syntax to import conditional operators from ufl package? - #4 by dokken

from dolfin import *
N = 10
mesh = UnitIntervalMesh(N)
V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)
f = project(x[0], V)
def u(x0, degree):    
    step = conditional(lt(x[0], x0), 1, 0)
    return  assemble(step*f*dx(domain=mesh, degree=degree))
for p in [0.05, 0.1, 0.23, 0.28]:
    for degree in [2**i for i in range(1,6)]:
        exact = 1/2*(p**2)
        approx = u(p, degree)
        print(p, degree, exact, approx, (exact-approx)/exact )

giving the following for N=10

0.05 2 0.0012500000000000002 0.0010566243270259353 0.15470053837925188
0.05 4 0.0012500000000000002 0.0003130601816090507 0.7495518547127595
0.05 8 0.0012500000000000002 0.0006078258433900407 0.5137393252879675
0.05 16 0.0012500000000000002 0.0008600886103124731 0.31192911175002164
0.05 32 0.0012500000000000002 0.0010324149186984053 0.17406806504127595
0.1 2 0.005000000000000001 0.0049999999999999975 6.938893903907227e-16
0.1 4 0.005000000000000001 0.004999999999999996 1.040834085586084e-15
0.1 8 0.005000000000000001 0.0049999999999999975 6.938893903907227e-16
0.1 16 0.005000000000000001 0.004999999999999997 8.6736173798840335e-16
0.1 32 0.005000000000000001 0.004999999999999999 3.4694469519536137e-16
0.23 2 0.02645 0.031056624327025925 -0.17416349062479863
0.23 4 0.02645 0.025868615737164594 0.021980501430450168
0.23 8 0.02645 0.027763381398945595 -0.04965525137790526
0.23 16 0.02645 0.025557752394857115 0.03373336881447584
0.23 32 0.02645 0.02605440254041938 0.014956425693029152
0.28 2 0.039200000000000006 0.044999999999999984 -0.14795918367346883
0.28 4 0.039200000000000006 0.036979726848275704 0.056639621217456665
0.28 8 0.039200000000000006 0.041501668016300265 -0.058716020823986206
0.28 16 0.039200000000000006 0.037494421258372934 0.043509661776200796
0.28 32 0.039200000000000006 0.03903467013137407 0.004217598689437088

For N=20

0.05 2 0.0012500000000000002 0.0012499999999999994 6.938893903907227e-16
0.05 4 0.0012500000000000002 0.001249999999999999 1.040834085586084e-15
0.05 8 0.0012500000000000002 0.0012499999999999994 6.938893903907227e-16
0.05 16 0.0012500000000000002 0.0012499999999999992 8.6736173798840335e-16
0.05 32 0.0012500000000000002 0.0012499999999999998 3.4694469519536137e-16
0.1 2 0.005000000000000001 0.004999999999999998 5.20417042793042e-16
0.1 4 0.005000000000000001 0.0049999999999999975 6.938893903907227e-16
0.1 8 0.005000000000000001 0.004999999999999999 3.4694469519536137e-16
0.1 16 0.005000000000000001 0.004999999999999998 5.20417042793042e-16
0.1 32 0.005000000000000001 0.005000000000000001 0.0
0.23 2 0.02645 0.025264156081756483 0.04483341845911222
0.23 4 0.02645 0.027856042823180036 -0.053158518834783904
0.23 8 0.02645 0.026929734238625287 -0.01813740032609774
0.23 16 0.02645 0.02624702013695705 0.007674096901434859
0.23 32 0.02645 0.026831735198445915 -0.01443233264445798
0.28 2 0.039200000000000006 0.03776415608175648 0.03662867138376349
0.28 4 0.039200000000000006 0.04091159837873558 -0.04366322394733602
0.28 8 0.039200000000000006 0.039785289794180834 -0.014930862096449692
0.28 16 0.039200000000000006 0.03895341973383283 0.0062903129124279285
0.28 32 0.039200000000000006 0.039664592624127276 -0.011851852656307919

Note that as long as the point x0 does not align with the interval, you do not obtain an exact value.

1 Like
6

Thank you, dokken. The model I’m developing relies on a second-order stress equation on a boundary of a 2D (rectangle) domain.

I integrate the equation by parts to reduce the order, and the process yields the term t \cdot \phi \Big |^{S_f}_{S_i}, in which t is the tangent vector, \phi is a test function from a vector function space and S_f and S_i are the end and begin of the boundary S, respectively.

Considering the value of \phi(S_f)= \phi(S_i) = 1, the terms t \cdot \phi (S_f) and - t \cdot \phi (S_i) should be included in the residual of the element located at S_f and S_i. I’ve been struggling to find a solution to this problem, as I don’t know how to mark an edge using a MeshFunction and I still don’t know how could I work with the residual of an element alongside the weighted residual from the problem formulation. I defined a dummy function t_edges that takes the tangent value on the edges and is zero everywhere else, and the I add the term inner(t_edges, phi)*ds(marker) into the variational formulation. However this approach is not working. Do you have any thoughts on how I could do this?