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

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

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])
do something


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

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

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

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
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("\n")
else:
break


and the output

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

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

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

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

x, y, f(x,y) [ 0.04166667  0.00520833] 0.00176323784722
residual at crds 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?