Get average value of a function, cell per cell, and facet per facet, with non regular mesh (cell of different sizes)
3 months ago by
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.
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
3 months ago by
No one has an idea? I think that's a pretty common issue people might be interested with.
3 months ago by
Please login to add an answer/comment or follow this question.