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


272
views
2
7 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

2 Answers


0
7 months ago by
Hello there,
No one has an idea? I think that's a pretty common issue people might be interested with.
0
7 months ago by
I think ufl.cell_avg(), ufl.facet_avg() might be a good starting point.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »