Finite Element Assembly using map over elements

Specific Domains Numerics
1

Hi,

I am wondering whether the loop over the elements in a finite element assembly can be replace by a flavor of map over the list of elements. I am happy to consider the linear elements on a one-dimensional mesh to start from.

Assume thus having on the interval (or computational domain) 0 <= x <= 1 three elements (or subintervals) 0 <= x <= 1/3, 1/3 <= x <= 2/3 and 2/3 <= x <= 1. We thus have as connectivity matrix

e = [1 2; 2 3; 3 4]

and a list of coordinates p = [0 1/3 2/3 1].

Is there a way to compute the length (area) of each subinterval by a map on the matrix e?

Thanks. Domenico.

2

Have you seen Gridap.jl? There’s no need to implement FEM from scratch.

3

It can be done. There may not be any advantage to re-doing something that works well.

4

Dear Petr and Steven,

Sincere thanks to both of you for getting in touch.

Yes, I have seen GridApp.jl (and similar such as JuliaFEM.jl and JuAFEM.jl) and the overview given here JuliaPDE · GitHub . And yes, I do agree that there might be no intrinsic need to replace the loop over elements by a map.

My interest, however, is in gaining a better understanding of Julia in general and of coding PDE solvers in Julia in particular. Replacing the for-loop by a map might constitute a fun coding exercise. Being able to go beyond a blind use of a package like GridApp.jl might put me in the position to contribute to the package in a next stage.

I therefore remains curious in the question as originally posted.

Greetings, Domenico.

5

Fair enough. Here is to get you started:

julia> map(i -> p[e[i, 2]]-p[e[i, 1]], 1:size(e, 1))
3-element Vector{Float64}:
 3.33333e-01
 3.33333e-01
 3.33333e-01
1 Like
6

Wonderful, very nice!

This, however, assumes that the variables p and e are global variables to the anonymous function in the input variables i. Is there a way to avoid this? I tried a map over both e and p, leading to the error that e and p are of different size.

Can mapreduce be exploited to assemble the stiffness matrix, that is, can reduce be instructed to take node-element connectivity into account?

Thanks again, Domenico.

7

@petr: your advice has been instrumental. Thx!

Below how to remove dependence on global variables.

Below can be extended to elementary matrix contribution. The set of matrix contributions can subsequently be assembled into the global matrix using the sparse function.

Can mapreduce be used for the global assembly?

How do map and explicit for loop compare in assembling matrices?

Does map allow easy extension to pmap for parellel assembly?

function compute_element_area(elem_id, e, p)
  area_id = p[e[elem_id,2]] - p[e[elem_id,1]]
  return area_id
end 

elem_ids = [1:3] 
e = [1 2; 2 3; 3 4]
p = [0 1/3 2/3 1]

map(elem_ids) do elem_id 
    compute_element_area(elem_id, e, p)
end
2 Likes