[Gradient w.r.t MeshFunction]

General
1

Dear all:

I am new to fenics, and I wonder if it’s possible to find the gradient w.r.t the MeshFunction?

More specified:
For a given mesh, with stored simulation result from external numerical model, we denote the f_{i} \in R be the value of the i-th mesh center. So we have no explicit expression for f but only mesh-center value.

Consider for the nonlinear poisson problem (g is known function, boundary setting is same as demo nonlinear poisson case), my question is : how we formulate the || \nabla f||^{2}, or gradient of f in fenics, with only information of the mesh and the cell center value ?

\nabla \cdot \left( {||\nabla f|{|^2}\nabla u} \right) = g

Tried solution:
[1] In demo case ‘tensor weighted poisson’, an example of the meshFunction may solve this problem. However, it reports that '‘MeshFunctionDouble can not be converted to any UFL type’ .

Thanks for anyone watched this thread.
Best regards,
Anthony

2

You would have to convert your mesh function into a function in the (“DG”, 0) space, see for instance: How to define different materials / importet 3D geometry (gmsh) - #2 by dokken

4

Hello dokken, thanks for your kindly reply!

As your suggestion, did you mean the code should be like this?

f = [val_cell_1, val_cell_2, ..., val_cell_nxn]
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "DG", 0)
v = Function(V)
x = V.tabulate_dof_coordinates()

for i in range(x.shape[0]): # cell_number x dimension
    v.vector().vec().setValueLocal(i, f[i])
grad_of_v  = grad(v)
do something

Sorry for this kind of reply since something wrong with my docker and i currently not able to verified by myself.

Thanks for your help!

5

You have written sudo-code almost identical to the post I linked you to, so yes, it should do what you expect.

6

Yeah, the code above almost referenced from your post, just to verfify the code behave as i expect.

Thank you for your kindly help! It really helps me a lot.
Best regards!
Anthony

7

Hello dokken, I tried previous code but something wired happened,

  1. the mesh is uni-square mesh with 64 \times 64
  2. the function represent synthetic mesh value, here f(x,y):=x^{2}+y^{2}, so \nabla{f}=[2x, 2y]^{T}

The gradient function defined by following code produces 0 at each point, but it should not

import numpy as np
from dolfin import *
# Toy example 
n = 64
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "DG", 0) 
f = Function(V)
x = V.tabulate_dof_coordinates()

# Define f := x^2 + y^2 
for i in range(x.shape[0]):
    x_, y_ = x[i]
    f.vector().vec().setValueLocal(i, np.sum(np.power(x[i], 2)))

# Loop to check the gradient and coordinate
# Only first 5 res to check
Lx, Ly = grad(f)
for idx, crds in enumerate(x):
    if idx < 5:
        print("x, y, f(x,y)", crds, f(crds))
        print("residual at crds", f(crds) - np.sum(np.power(crds, 2)))
        print("grad at x is", Lx(crds))
        print("grad at y is", Ly(crds))
        print("\n")
    else:
        break

and the output

x, y, f(x,y) [ 0.01041667  0.00520833] 0.000135633680556
residual at crds 0.0

grad at x is 0.0
grad at y is 0.0


x, y, f(x,y) [ 0.00520833  0.01041667] 0.000135633680556
residual at crds 0.0
grad at x is 0.0
grad at y is 0.0


x, y, f(x,y) [ 0.02604167  0.00520833] 0.000705295138889
residual at crds 0.0
grad at x is 0.0
grad at y is 0.0


x, y, f(x,y) [ 0.02083333  0.01041667] 0.000542534722222
residual at crds 0.0
grad at x is 0.0
grad at y is 0.0


x, y, f(x,y) [ 0.04166667  0.00520833] 0.00176323784722
residual at crds 0.0
grad at x is 0.0
grad at y is 0.0

The function value is correctly assigned to mesh and the function is properly constructed, but obviouly the gradient should not be 0, so what may be the problem?

Thanks for your help

8

An oversight perhaps is using a DG space V^p_\text{DG} of degree p=0. Naturally the broken gradient operator defined on a mesh \mathcal{T}_h = \{\kappa\} is \nabla_h := \nabla|_{\kappa \in \mathcal{T}_h} when applied to a function \nabla_h u_h(x) = 0, \forall u_h \in V^0_\text{DG}. So the trick would be to decide how to interpolate your piecewise constant data to something which varies intra-cell. Or choose some other method for interpolation/projection of your cell-wise data which varies whilst conforming with continuity inter-cell.

1 Like
9

That might be the solution,I should check it more precicely and spend some time checking the document of Fenics.

Thans for your kindly help!
Have a nice day :slight_smile: