how can I add local contributions to a system matrix in appropriate way?

Question

Hello,

I have been implemented a stabilized FEM whose stabilization term is the sum over the
interior facet of some facet integrals.
Until now I have the following code

for facet in facets(mesh):
      subdomains.set_all(0)
      facet_integration.set_all(0)
      if interior_facets[facet] == 1:
        facet_integration[facet] = 1
        if extremum_domain_parts[facet] == 1:
           t0 = 1.0
        else:
          t0 = alpha(u_h)

      AA1 = assemble(a_ld3, interior_facet_domains = facet_integration)
      AA2 = assemble(a_lps3, interior_facet_domains = facet_integration)
     AA = AA + t0*AA1 + (1-t0)*AA2

but it is very slow. Can I implement it in a better way?

thank you in advance
Fotini

Comments

Please provide a complete problem description and maybe even a code-snippet which can be used to reproduce your problem (see also). You did not even specify what kind of stabilization you want to use or what the code is supposed to do.

a stabilized FEM whose stabilization term is the sum over the 
interior facet of some facet integrals

This sounds like you want to implement some interior-penalty-type stabilization. This is covered in the demos and much easier to accomplish than looping manually over the facets.

---

thank you for your answer! it is a kind of interior penalty but the coefficients are defined on facets and not element-wise. If there was something like Coefficient(facet) it would be easy.

Answer 1

As Christian said, it is necessary to avoid your approach as it requires triggering assembler for every facet in the loop. Check Expression documentation.

User-defined expressions by subclassing allows you to define an expression depedent ufc_cell.local_facet. To speed it up you can additionally switch to Complex user-defined JIT-compiled Expressions.