Conversion of minimal code snippet from Python to C++

Question

This is a question about understanding a piece of random code that does not necessarily require knowledge of it's theory. This is very specific and may not be of use to the community in general but some may find it interesting. I have shown a snippet of python code below. In my attempt to take this code from python to C++, I have marked the lines which can go directly to UFL and those which require C++ conversion and then go to main.cpp.

Du_elastic_energy = derivative(elastic_energy, u, w) # Goes as it is to UFL
Dv_kinetic_energy = derivative(kinetic_energy, v, w)  # Goes as it is to UFL
M, F_tot = lhs(replace(Dv_kinetic_energy, {v: du})), - Du_elastic_energy     # Goes as it is to UFL
M_lumped = action(M, Constant(1.0))  # Mass lumping  # Goes as it is to UFL
ML = PETScVector()          # Goes to main.cpp
FT = PETScVector()          # Goes to main.cpp
assemble(M_lumped, tensor=ML)        # Goes to main.cpp
def solve_a():                                       # Goes to main.cpp
    assemble(F_tot, tensor=FT)
    a.vector().set_local(FT.array() / ML.array())

Is the following a correct C++ conversion of the above FEniCS python code?

auto ML = std::make_shared<dolfin::PETScVector>();   # M_lumped.(new dolfin::PETScMatrix());
auto FT = std::make_shared<dolfin::PETScVector>();
dolfin::Assembler assembler;
assembler.assemble(ML, *M_lumped);   # assembler_->assemble(*A);

void solve_a()
{
assembler.assemble(FT, *F_tot);
a.vector() = FT.vector() / ML.vector()
}

Comments

What types are M_lumped, F_tot, and a?

---

Hello senseiwa,

M is mass matrix, F_tot is sum of internal and external forces and 'a' is acceleration vector.

The only relevant part of code before this snippet which may be relevant is:

elastic_energy = 0.5*inner(grad(u), grad(u))*dx
kinetic_energy = 0.5*inner(v, v)*dx

What the author is trying to do is to solve for the acceleration

Ma = Fext - Fint = Fext - Ku

Action() and Replace() are defined in UFL manual here:

https://fenicsproject.org/pub/documents/ufl/ufl-user-manual/ufl-user-manual.pdf

User tianyikillua has provided details of this same algorithm here:

https://fenicsproject.org/qa/8006/elementwise-multiplication-parallelisation-difference

Answer 1

Since I have no knowledge about types (not mathematical domains, I meant in the comment "C++ types"), this is my best guess:

// Somewhere you do this
{
    // Shared pointer    
    auto ML = std::make_shared<dolfin::PETScVector>();
    // Shared pointer
    auto FT = std::make_shared<dolfin::PETScVector>();

    dolfin::Assembler assembler;
    assembler.assemble(ML, *M_lumped);
}

// Then you define this function
void solve_a()
{
    assembler.assemble(FT, *F_tot);

    // I don't see anywhere what a is I suppose it is **NOT** a pointer

    // FT and ML are (shared) pointers, vector() returns a shared pointer
    *a.vector() = *FT->vector() / *ML->vector();
}

Obviously if a is indeed a pointer, you must use the syntax for FT.

Cheers!