LocalSolver does not support facet integrals

[Tormod Landet]

It seems that LocalSolver does not support facet integrals. I will look at how hard it is to implement this tomorrow. Any hints are appreciated. I already have a pure Python implementation, but do not know my way around the dolfin C++ code

Ref: qa question

[Martin Sandve Alnæs]

Good luck. Look for these landmarks to navigate after in LocalSolver.cpp:

Iteration over cells in mesh:

 for (CellIterator cell(mesh); !cell.end(); ++cell)

Computation of the element matrix on a cell:

 integral->tabulate_tensor(A...)

Look at assemble_cells vs assemble_exterior_facets in Assembler.cpp for guidance on the differences.

[Prof Garth Wells]

It should be quite simple to implement. Do you want interior facet integrals or just exterior facet integrals?

[Marie Elisabeth Rognes]

You could also take a look at the LocalAssembler classes:

dolfin/adaptivity/LocalAssembler.cpp

[Tormod Landet]

Tanks for the pointers, I will look at them now.

I want to support at least all local facet integrals, ie uvds, u('+')v('+')dS and u('-')v('-')dS. I am currently not so interested in non-local integrals such as u('+')v('-')dS

[Martin Sandve Alnæs]

You can't see the difference between the dS integrals.

[Tormod Landet]

Some questions:

• I will try to use LocalAssembler::assemble instead of duplicating that code, does this seem reasonable?
• In solve_local_rhs the following code appears when putting the solution into the user provided function:

x.set_local(x_e.data(), dofs_a0.size(), dofs_a0.data());

Doesn't this assume that the function space of the trial and test functions are the same? I would think that it should be dofs_a1 that should be used?

[Martin Sandve Alnæs]

Reusing LocalAssembler seems very reasonable, without known the details of this code.

There's an assertion on line 147 in LocalSover.cpp that the test and trial spaces have the same dimension:

dolfin_assert(dofs_a0.size() == dofs_a1.size());

If they were not the same, you would have a non-square matrix so the solving part wouldn't work out so well, unless you want to do a local least squares solve or something?

[Martin Sandve Alnæs]

However, I agree that it should be dofs_a1, the trial space.

Either the spaces truly are exactly the same, and then only one of the dofs vectors is needed, or they just have the same size and can have e.g. different ordering, and dofs_a1 must be used.

[Tormod Landet]

I got it working (unit tests at least) by removing quite some code and using LocalAssembler. The only real change is that I needed to make the local A matrix column major, it was row major for some reason, probably efficiency of the code LocalSolver was originally used for? If it must remain row major then I would either have to reimplement LocalAssembler with row major matrices or change LocalAssembler.

With regards to the dofs_a0 vs dofs_a1 situation it is the latter case you mention that is important to me. The matrices are not square globally, but they are square locally with different function spaces for trial and test functions. As an example this is a unit test I have added to test_local_assembler.py (this unit test is the one I aim to make work):

def test_solve_local_rhs_facet_integrals():
    mesh = UnitSquareMesh(1, 1)

    Vu = VectorFunctionSpace(mesh, 'DG', 1)
    Vv = FunctionSpace(mesh, 'DGT', 1) 
    u = TrialFunction(Vu)
    v = TestFunction(Vv)

    n = FacetNormal(mesh)
    w = Constant([1, 1])

    a = dot(u, n)*v*ds
    L = dot(w, n)*v*ds
    for R in '+-':
        a += dot(u(R), n(R))*v(R)*dS
        L += dot(w(R), n(R))*v(R)*dS

    A = assemble(a)
    b = assemble(L)
    U = Function(Vu)

    ls = LocalSolver(a, L)
    ls.solve_local_rhs(U)

    print U.vector().array()
    err = abs(U.vector().array() - 1).sum()
    print err
    assert err < 1e-15

[Tormod Landet]

I moved this to a pull request: https://bitbucket.org/fenics-project/dolfin/pull-requests/244/support-facet-integrals-in-localassembler

[Prof Garth Wells]

You'll need to make the mesh bigger for the unit test, otherwise it will fail in parallel.

Also, don't use NumPy to manipulate arrays. Use the GenericVector interface, which works in parallel.

[Tormod Landet]

Will do. Is 4x4 enough? How many processors are you testing on?

Should I squash commits on the branch or just add more?

[Prof Garth Wells]

4x4 should be ok. The automated tests run on 3 processes, so I aim for at least 10 cells in total so that the partitioner puts at least one cell on each process (otherwise some partitioner crash).

[Tormod Landet]

Fix is in master after merging pull request #244

[Jack Hale]

Thanks for this @trlandet, I was working on fenics-shells today with a projection between functions in CG^2 to NED_1^1 and it will be nice to shift it over to LocalSolver instead of using KrylovSolver for the post-processing.

On a related note, there are some strange choices on row vs column major ordering in the LocalAssembler interface. Perhaps I will open a separate issue though for discussion.

[Tormod Landet]

Great that this is useful to you, @jackhale!

I changed the default ordering in LocalAssembler to match what was used in LocalSolver previously. By using the storage order that tabulate_tensor expects we can avoid copying data from one structure to the other when assembling pure dx integrals. I did not benchmark that change in isolation, but the goal was to not make a major impact on the speed of the LocalSolver for the use cases where it had been previously used (i.e only cell integral projections). For ds and dS integrals the copying must take place anyway so that previously assembled dx integrals for the same cell are not overwritten, so is only for dx integrals that the storage order can make a speed difference. One exception is that perhaps the LU/Cholesky factorization may be faster (or better tested at least, according to the documentation) when using the default storage order of Eigen (which we are not using now).

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)