Feature request: modify exising matrix

[VikingScientist]

Sorry if this functionality already exist, but I could not find it. Say I have two functions, which both should contribute to my system matrix and/or force vectors. Is it possible to chain these, so the second would be added to the results from the first, instead of just overwriting. In the example below I was hoping to get 3 along the diagonal, but I only get 2.

element Vertex
end

extern verts : set{Vertex};

func assemble_one(v : Vertex) -> (K : matrix[verts,verts](float))
    K(v,v) = 1.0;
end

func assemble_two(v : Vertex) -> (K : matrix[verts,verts](float))
    K(v,v) = 2.0;
end

export func main()
    var K = map assemble_one to verts reduce +;
    K     = map assemble_two to verts reduce +;
    print K;
end

[fredrikbk]

We're trying to encourage a more functional style so creating two matrices and adding them together is preferred:

A = map assemble_one to verts reduce +;
B = map assemble_two to verts reduce +;
C = A + B;

or

C = (map assemble_one to verts reduce +) + (map assemble_two to verts reduce +);

Our intention (and belief) is that our compiler will eventually be able to fuse these together. We have some stuff in the pipeline that we think will lead to such advanced capabilities.

To directly answer your question, there is no such syntax at the moment. We could add it as an extension if it's really needed to get the performance you seek, but we think it's not desirable longer term.

[VikingScientist]

This works fine if the matrix has the same sparsity pattern. I get segmentation fault if I try with different patterns. A typical scenario here is to assemble all interior terms (body forces) as one function call and all boundary terms (external forces) as a second function call. The boundary terms will not create any new non-zero entries in the global matrix, but rather be a lot more sparse, i.e. edge functions already create internal contributions in addition to their now edge contributions. The following code will cause a crash:

element Vertex
end
element Line
end

extern verts : set{Vertex};
extern line : set{Line}(verts,verts);

func assemble(v : Vertex) -> (K : matrix[verts,verts](float))
    K(v,v) = 1.0;
end

func assemble_edge(l : Line, v : (Vertex*2)) -> (K : matrix[verts,verts](float))
    K(v(0),v(0)) = 1.0;
    K(v(1),v(1)) = 1.0;
end

export func main()
    K1 = map assemble to verts reduce +;
    K2 = map assemble_edge to line reduce +;
    var K = K1 + K2;
    print K;
end

Interestingly if you swap the sum, i.e. var K = K2 + K1, it will not crash, but it will create the wrong output.