Average a CellField over a dimension

[hgpeterson]

Hello,

Is there an efficient way to average a CellField over one or multiple dimensions? For instance, I have a SingleFieldFEFunction defined on a 3D mesh that I'd like to average in x and y to just get a 1D function of z.

I suppose one way to do this is by using evaluate on a uniform grid of points and then averaging over that, but I was wondering if there is a better way.

Thanks!

[JordiManyer]

There is no easy way of doing this. I think slicing up domains is something we could expand on. However, and as always, it is involved but not impossible.

The basic structure should go like this:

using Gridap
using Gridap.Geometry

model = CartesianDiscreteModel((0,1,0,1,0,1),(4,4,4))

order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe);

function create_face_mask(model::DiscreteModel{Dc,Dp},mask) where {Dc,Dp}
  topo = get_grid_topology(model)
  node_coords = Geometry.get_vertex_coordinates(topo)
  face_to_node = Geometry.get_faces(topo,Dc-1,0)
  face_coords = lazy_map(Broadcasting(Reindex(node_coords)),face_to_node)
  bgface_to_mask = map(mask,face_coords)
  return bgface_to_mask
end

function my_mean(z,f)
  mask(p::Point{3}) = abs(p[3]-z) < 10*eps(z)
  mask(nodes::AbstractVector) = all(mask,nodes)

  bgface_to_mask = create_face_mask(model,mask)
  Γ = BoundaryTriangulation(model,bgface_to_mask)
  dΓ = Measure(Γ,order)
  return sum(∫(f)dΓ)/sum(∫(1)dΓ)
end

uh = interpolate(x -> x[3],V)
my_mean(0.0,uh) # 0.0
my_mean(1.0,uh) # 1.0

[hgpeterson]

Hi Jordi,

Thanks for the quick reply! Being able to slice up domains quickly would definitely be super useful if implemented.

It looks like this example you sent works will for a structured grid where you know that abs(p[3]-z) < 10*eps(z) will return a valid slice of points. Is there an hope if one has an UnstructuredDiscreteModel()?

Henry

[JordiManyer]

This relies on selecting a subset of faces. If you don't have that, you cannot do things efficiently using regular FEM. You would have to go to embedded methods, which you can find in GridapEmbedded.jl, use that to embed your cut into your background mesh and then do the integrals.

[hgpeterson]

Ah got it. I wasn't aware of GridapEmbedded.jl, so I'll have to check that out. Thanks for the help!