Local assembling

Question

Hello,

for sensitivity analysis of eigenvalues respect to material parameters I need to assemble differentiated form on each cell. I use the following loop:

for i, cell in enumerate(cells(Omega)):
    cell_id = cell.index()
    # extract eigenvector for given cell
    eigvect_e = eigvects[-1, dofmap.cell_dofs(cell_id)].real
    # extract  differentiated stiffness matrix on cell
    dk_t = assemble_local(dk * dx, cell)
    # compute sensitivity
    edk[i] = dot(dot(eigvect_e, dk_t), eigvect_e.T)

However it is too slow with higher number of cells than a few thousands. Please, any idea how to get speed up?

Petr

Answer 1

I'm not aware of any straightforward ways speed this up, unless you write a function in C++ and let Instant to the wrapping.

You can find an example of wrapping a large function at https://bitbucket.org/unilucompmech/fenics-shells/src/76149165aaadfd2cfdc7b06225e35681143d492f/fenics_shells/fem/?at=master

Comments

Thank you, Mr. Wells,
As you suggested, I have tried to write down a simple c++ function, but I stopped on the following problem - I declared a function similarly to link you posted:

#ifndef __LOCAL_ASSEMBLER_H
#define __LOCAL_ASSEMBLER_H

namespace dolfin
{
    void local_assembler(const Form& a, std::vector<double>& tensor);
}

#endif

and in execution file:

header_file = open("local_assembler/local_assembler.h", "r")
code = header_file.read()
header_file.close()
cpp_module = compile_extension_module(
    code=code, source_directory="local_assembler", 
    sources=["local_assembler.cpp"],
    include_dirs=[".", os.path.abspath("local_assembler")])

The headers are in source file.

However, after calling the function like:

k = inner(grad(v), sigma(u)) * dx
cpp_module.local_assembler(k)

I get the error:
TypeError: in method 'local_assembler', argument 1 of type 'dolfin::Form const &'

Please, any idea why the data type is not correct?

I am not much familiar with c++, so please excuse if the error is trivial.

Petr

---

Hi,

a normal form is a ufl.Form. When you compile it you get a dolfin.Form

 a = ...
 a2 = dolfin.Form(a)

See also _create_dolfin_form here:
https://bitbucket.org/fenics-project/dolfin/src/07d3743c91637ff1c2f90474eac70a283c51101c/site-packages/dolfin/fem/assembling.py?at=master&fileviewer=file-view-default

and
https://bitbucket.org/fenics-project/dolfin/src/07d3743c91637ff1c2f90474eac70a283c51101c/site-packages/dolfin/fem/assembling_local.py?at=master&fileviewer=file-view-default

---

Hi,

Thanks, now it works correctly.