### Get average value of a function, cell per cell, and facet per facet, with non regular mesh (cell of different sizes)

359
views
2
9 months ago by
Hello everyone,
My question is quite generic, but i could not find a satisfactory solution yet.

Let's say we have a mesh whom cell and facet sizes are not unique: some cell are bigger than others (thus each dof doesnt have the same weight from a statistical point of view). We have a function solution u calculated on this mesh.
If i want to calculate the average value, i simply do an assemble of the solution u divided by the mesh volume. But i want to get average value cell per cell and facet per facet, for a subset of the facets, ds(id_)

- Question 1: how to extract an array, whom row i is equal to the average value of u within the cell i (i.e. assemble u*dx(this cell)/volume this cell), for all the cell of the mesh, with the best efficient way possible ?

Right now, i am doing the following (it works), but i would like to know if someone know a faster version:

# Calculate cell volume per cell
def cell_volume(mesh):
# We use the Discontinous Galerkin element
V_mesh_DG = FunctionSpace(mesh, 'DG', 0)
# Create associate function
Volume_cell=Function(V_mesh_DG)
test_function_DG = TestFunction(V_mesh_DG)
# Calculate volume cell per cell
c = Constant(1.0)
Volume_ = assemble(test_function_DG*c*dx)
Volume_cell_array = Volume_.array()
# Replace array
Volume_cell.vector()[:] = Volume_cell_array
return Volume_cell,Volume_cell_array

# Calculate function avergae value per cell
def fct_average_percell(u,u_DG,t_VG,V_DG,cell_volume_array):
# Calculate assemble cell per cell
u_=assemble(t_VG*u*dx)
# Divide by cell volume
u_array=u_.array()/cell_volume_array
# Replace array
u_DG.vector()[:] = u_array
return u_DG,u_array

### We have a mesh called mesh_
### We have a function solution called u_

### Create DG objects (so i dont have to define them each time)
V_mesh_DG = FunctionSpace(mesh_, 'DG', 0)
u_avg=Function(V_mesh_DG )
test_function_mesh_DG = TestFunction(V_mesh_DG )

### Calculate cell volume
Volume_cell,Volume_cell_array = cell_volume(mesh_)

### Calculate function average
u_avg,u_avg_array = fct_average_percell(u_,u_avg,test_function_mesh_DG ,V_mesh_DG ,Volume_cell_array)

​

- Question 2: how to extract an array, whom row i is equal to the average value of u for the facet i (i.e. assemble u*ds(this facet)/area this facet), ONLY for the facets marked with a certain value, with the best efficient way possible ?

I do not have a solution right now. If someone can help, that would be much appreciable.

I have seen this post that could help, but i didnot make it yet.
https://fenicsproject.org/qa/2822/how-to-calculate-facet-average

Best.

EDIT:

ok i found a method for the facets, but i feel a better solution exist out there.

def facet_surface(mesh):
# We use the Discontinous Galerkin element
V_mesh_DG = FunctionSpace(mesh, 'DG', 0)
# Create associate function
Area_face=Function(V_mesh_DG)
test_function_DG = TestFunction(V_mesh_DG)
# Calculate area face per facet
c = Constant(1.0)
Area_ = assemble(test_function_DG*c*ds)
Area_facet_array = Area_.array()
# Replace array
Area_face.vector()[:] = Area_facet_array
return Area_face,Area_facet_array

# Calculate function avergae value per facet
def fct_average_perfacet(u,u_DG,t_VG,V_DG,facet_area_array,ds_,id_interface):
# Calculate assemble cell per cell
u_=assemble(t_VG*u*ds_(id_interface))
# Divide by facet area
u_array=u_.array()/facet_area_array
# Replace NAN by 0, so you can plot it later
u_array = np.nan_to_num(u_array)
# Replace array
u_DG.vector()[:] = u_array
return u_DG,u_array​

# Example (you have to define your ds and mark the facetfunction with id_ first)
u_facetavg=Function(V_mesh_DG)
Area_facet,Area_facet_array = facet_surface(mesh)
u_facetavg,u_facetavg_array = fct_average_perfacet(u,u_facetavg,test_DG,V_mesh_DG,Area_facet_array,ds_,id_)

# You need to remove all the 0 for further analysis
index=np.nonzero(u_facetavg_array)
u_facetavg_array=u_facetavg_array[index]



Community: FEniCS Project