need a way to loop over quadrature points

[stevengj]

@hammy4815 and I were looking for a lower-level abstraction that gives us access to the quadrature points used when assembling a matrix. (This is for topology optimization where the degrees of freedom are specified on some external grid, not on the FE mesh.)

I think all we need is a API like the following, given a quadratic form a and two FEFunctions u,v on some mesh Ω:

a(u, v) = ...some quadratic form...
quadratureloop(a, u, v, Ω) do aₓ, x, wₓ
     ... do something with aₓ(u, v) at the quadrature point x with quadrature weight wₓ ...
end

so that, for example, ∫a(u,v)dΩ could be computed with:

integral = 0.0
quadratureloop(a, u, v, Ω) do aₓ, x, wₓ
     integral += aₓ * wₓ
end

(In our case we want to do something different with the quadrature points and weights, not just an integral)

[stevengj]

(You could also have the analogous API corresponding to assemble_vector instead of assemble_matrix, where you just pass 1 space/function instead of two, i.e. quadratureloop(f, a, u, Ω) corresponding to ∫a(u)dΩ. We don't need this ourselves, but it is good to think about.)

[fverdugo]

Hi @stevengj,

It is possible to get the values of the bilinear and linear forms at quadrature points, but it requires some effort.

Accessing the elemental matrices (once integrated) is much more easy. Is this perhaps enough for your needs?

using Gridap
using Gridap.FESpaces
U = ... # trial space
V = ... # test space
Ω = ... # triangulation
dv = get_fe_basis(V)
du = get_trial_fe_basis(U)

cell = 100
a(du,dv)[Ω][cell] # elemental matrix of the cell number 100 in the triangulation Ω

[stevengj]

We need the individual quadrature points, not just the cell-integrated elemental matrices.

(The reason is that backpropagating from assembly to an interpolation function from an external grid requires us to map each quadrature point to the external grid via the interpolation weights for that point, and these vary from point to point.)

[stevengj]

Note that I originally suggested a lower-level API, but then deleted it in favor of the simpler API above, but maybe it's clearer to give the lower-level API. In particular, we just need the same quadrature points and bilinear-form values that you use internally in assemble_matrix. In particular, it would suffice to have the lower-level API:

quadratureloop(a, U, V) do aₓ, x, wₓ, i, j
     # do something the value aₓ of a(uᵢ, vⱼ) at x, along with quadrature weight wₓ,
     # for overlapping basis functions i and j
end

so that assemble_matrix could conceptually be replaced by:

A = zeros(m, n)
quadratureloop(a, U, V) do aₓ, x, wₓ, i, j
    A[i,j] += aₓ * wₓ
end

(This might be useful for other reasons, e.g. to use custom matrix data structures.)

Then, for example, ∫a(u,v)dΩ could be computed by:

u_vec, v_vec = get_vector(u), get_vector(v)
U, V = get_space(U), get_space(V)
integral = 0.0
quadratureloop(a, U, V) do aₓ, x, wₓ, i, j
    integral += aₓ * wₓ * u_vec[i] * v_vec[j]
end

[fverdugo]

In this case, you can access the data you mention as follows.

module MWE

using Gridap

function get_cell_weights(dΩ::Gridap.CellData.Measure)
    get_cell_weights(dΩ.quad)
end

function get_cell_weights(quad::Gridap.CellData.CellQuadrature)
  if quad.data_domain_style == PhysicalDomain() && quad.integration_domain_style == PhysicalDomain()
    quad.cell_weight
  elseif quad.data_domain_style == ReferenceDomain() && quad.integration_domain_style == PhysicalDomain()
    cell_map = get_cell_map(quad.trian)
    cell_Jt = lazy_map(∇,cell_map)
    cell_Jtx = lazy_map(evaluate,cell_Jt,quad.cell_point)
    cell_m = lazy_map(Broadcasting(Gridap.TensorValues.meas),cell_Jtx)
    lazy_map(Broadcasting(*),quad.cell_weight,cell_m)
  elseif quad.data_domain_style == ReferenceDomain() && quad.integration_domain_style == ReferenceDomain()
    quad.cell_weight
  else
    Gridap.Helpers.@notimplemented
  end
end

domain = (0,1,0,1)
cells = (4,4)
model = CartesianDiscreteModel(domain,cells)
Ω = Triangulation(model)
k = 1
dΩ = Measure(Ω,4*k) # more points than needed to make number of points different than number of local shape functions
reffe = ReferenceFE(lagrangian,Float64,k)
V = FESpace(model,reffe)
U = V
a_integrand(u,v) = ∇(u)⋅∇(v)
dv = get_fe_basis(V)
du = get_trial_fe_basis(U)

p = get_cell_points(dΩ)
cell_ax = a_integrand(du,dv)(p)
cell_w = get_cell_weights(dΩ)
cell_x = Gridap.CellData.get_array(p)

c1 = array_cache(cell_ax)
c2 = array_cache(cell_w)
c3 = array_cache(cell_x)

ncells = length(cell_ax)
for cell in 1:ncells
    ax = getindex!(c1,cell_ax,cell)
    w = getindex!(c2,cell_w,cell)
    x = getindex!(c3,cell_x,cell)
    @show cell
    display(ax) # ax[q,:,:] local matrix at point q
    display(w) # w[q] integration weight at point q
    display(x) # x[q] physical coordinate of integration point q
end


end # module

[hammy4815]

This is great, Thank you Francesc! I think with this I can properly construct the gradient I need. Just to ensure I am parsing this correctly, if I wanted to compute the integral as in Steven's example above, it might look like:

integral = 0.0
for cell in 1:ncells
    ax = getindex!(c1,cell_ax,cell)
    w = getindex!(c2,cell_w,cell)
    x = getindex!(c3,cell_x,cell)
    a_mat = ax * w # tensor vector product
    integral += v_vec[...] * a_mat * u_vec[...] # quadratic form
end

[fverdugo]

Just to ensure I am parsing this correctly, if I wanted to compute the integral as in Steven's example above, it might look like:

You still need the loop over integration points.

for cell in 1:ncells
    ax = getindex!(c1,cell_ax,cell)
    w = getindex!(c2,cell_w,cell)
    x = getindex!(c3,cell_x,cell)
   for q in 1:length(w)
        ax[q,:,:] # This creates a matrix, you can use a view instead
        w[q]
        x[q]
   end
end

Feel free to add the new function get_cell_weights in this file. I it can be useful also for others.

https://github.com/gridap/Gridap.jl/blob/master/src/CellData/DomainContributions.jl

[hammy4815]

Feel free to add the new function get_cell_weights in this file. I it can be useful also for others.

https://github.com/gridap/Gridap.jl/blob/master/src/CellData/DomainContributions.jl

Gotcha, will do.

I'm trying to figure out how ax[q,:,:] maps to the global matrix acquired from assemble_matrix(). I need to find the indices in u_vec and v_vec that correspond to a given cell so I can perform something similar to the following:

integral = 0.0
for cell in 1:ncells
    ax = getindex!(c1,cell_ax,cell)
    w = getindex!(c2,cell_w,cell)
    x = getindex!(c3,cell_x,cell)
   for q in 1:length(w)
        a_mat = ax[q,:,:] * w[q]
        integral  += v_vec[indices] * a_mat * u_vec[indices]
   end
end

Basically, I'm unsure how to get indices given cell.

[JordiManyer]

I would start reading some of Gridap's asssembly-related code if I were you. The files containing most of what you want are CellQuadratures.jl (mostly for setting up integration), Assemblers.jl (mostly for assembly) and maybe DomainContributions.jl (which is the interface we have to manage multiple integrals in different triangulations).

The function that sets up the assembly for matrices can be found here. The rows and cols variables in that function are the ones you want.

[fverdugo]

Hi @hammy4815,

I am not sure what u_vec and v_vec are, but perhaps you are looking for something like this:

uh::FEFunction
vh::FEFunction
u_vec = get_free_dof_values(uh)
v_vec = get_free_dof_values(vh)
cell_dof_ids = get_cell_dof_ids(U,\Omega)
c4 = array_cache(cell_dof_ids)
integral = 0.0
for cell in 1:ncells
    ax = getindex!(c1,cell_ax,cell)
    w = getindex!(c2,cell_w,cell)
    x = getindex!(c3,cell_x,cell)
    indices = getindex!(c4,cell_dof_ids,cell)
   for q in 1:length(w)
        a_mat = ax[q,:,:] * w[q]
        # Following line assumes that you don't have Dirichlet BCs
        # Otherwise you need to handle negative values in indices
         integral  += transpose(v_vec[indices]) * a_mat * u_vec[indices]
   end
end

This would be equivalent to

A = assemble_matrix(a,U,V)
integral = transpose(v_vec)*A*u_vec

[hammy4815]

I am not sure what u_vec and v_vec are, but perhaps you are looking for something like this:

Great! This is what I needed. Thank you!

I will create a PR with the added functionality (get_cell_weights) as suggested in https://github.com/gridap/Gridap.jl/blob/master/src/CellData/DomainContributions.jl .