Assembly

As explained in Section Complementary rotation the assembly procedure of the residual vector and of the Jacobian matrix must be tweaked in order to check whether an element has some but not all the unknown rotation vectors changed into their complementary vector.

To this end, one need to write a customized assembly class in C++, compile it, and make it visible to the python interface of Dolfin.

An excellent starting point is an existing assembly class. Thus, one can start copying the files dolfin/fem/Assembler.h and dolfin/fem/Assembler.cpp into an unique file, say RotationAssembly.cpp, and save it into the directory where this project lives.

Then, replace all the occurence of Assembler class name with RotationAssembler to avoid any name clash.

The standard assembler class defines an assemble method as

void assemble(GenericTensor& A,
              const Form& a);

where A is where Form a should be assembled into.

The method assemble performs some initialization and checks, and then calls

void assemble_cells(GenericTensor& A, const Form& a, UFC& ufc,
                    std::shared_ptr<const MeshFunction<std::size_t>> domains,
                    std::vector<double>* values);

void assemble_exterior_facets(GenericTensor& A, const Form& a,
                              UFC& ufc,
                              std::shared_ptr<const MeshFunction<std::size_t>> domains,
                              std::vector<double>* values);

void assemble_interior_facets(GenericTensor& A, const Form& a,
                              UFC& ufc,
                              std::shared_ptr<const MeshFunction<std::size_t>> domains,
                              std::shared_ptr<const MeshFunction<std::size_t>> cell_domains,
                              std::vector<double>* values);

and

void assemble_vertices(GenericTensor& A, const Form& a, UFC& ufc,
                       std::shared_ptr<const MeshFunction<std::size_t>> domains);

that loops over all the cells, exterior_facets, interior_facets and vertices of the mesh and ssemble into the GenricTensor A the relative contributions coming from Form a.

The aim is to leave as much as possible unchanged. Looking into e.g. assemble_cells one easily identifies a loop over all the cells of the mesh

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

Inside this loop the assembly is performed by collecting all the degree of freedom values used by a given cell into a double** stored inside the helper class variable UFC ufc with

ufc.update(*cell, coordinate_dofs, ufc_cell,
           integral->enabled_coefficients());

and then compute the contribution of the cell with this call

integral->tabulate_tensor(ufc.A.data(), ufc.w(),
                          coordinate_dofs.data(),
                          ufc_cell.orientation);

Here method ufc.w() passes the double** to tabulate_tensor. The result of the integral is stored into ufc.A.data(), regardless of the rank of the form that is being assembled. If the rank is 0, one would get only a scalar, for a rank one form one would get a (dense) vector, and for a rank two form ufc.A.data() the resulting vector would store a dense matrix with row-major ordering.

At last, the result of the integral is assembled into the global tensor A with

if (is_cell_functional)
  (*values)[cell->index()] = ufc.A[0];
else
  A.add_local(ufc.A.data(), dofs);

where dofs stores the dofmap between the current cell DOFs and the global numbering of the GenericTensor A.

The assembly procedure needs to be modified between these three calls:

• immediately after ufc.update it is possible to check whether not all of the cell rotation vectors has been transformed into its complement, and eventually modify in place the right rotation vectors of ufc.w() by changing them again into their complement so that the rotation vector field can be interpolated.

• immediately after integral->tabulate_tensor is it instead possible to modify in place the values of ufc.A.data(), as explained in Section Complementary rotation:

  • if the form rank is equal to zero, do nothing;

  • if the form rank is equal to one then the equations that has as test function corresponding to the rotation vectors modified in the previous step should be multiplied on the left by \(\left(\partial\boldsymbol{\varphi}_C/\partial\boldsymbol{\varphi}\right)^T\);

  • if the form rank is equal to two then

    • the matrix rows that has a test function that corresponds to the rotation vectors modified in the previous step should be multiplied on the left by \(\left(\partial\boldsymbol{\varphi}_C/\partial\boldsymbol{\varphi}\right)^T\);

    • the matrix columns that has a trial function that corresponds to the rotation vectors modified in the previous step should be multiplied on the right by \(\left(\partial\boldsymbol{\varphi}_C/\partial\boldsymbol{\varphi}\right)\);

    • finally, one should add the term

      (1)\[\frac{2\pi}{\varphi^3}\left[ \left( \boldsymbol{f}_i\otimes\boldsymbol{\varphi} + \boldsymbol{\varphi}\otimes\boldsymbol{f}_i \right) + (\boldsymbol{\varphi}\boldsymbol{f}_i) \left( \boldsymbol{I} - 3\ \frac{\boldsymbol{\varphi}}{\varphi}\otimes\frac{\boldsymbol{\varphi}}{\varphi} \right) \right]\]

To do so, one needs to know where the components of the rotation coefficients are stored inside ufc.w(), and what are the corresponding positions, if any, inside the test and trial functions. The first step is to add a reference to the (sub)Function that approximates the rotation field as a third argument of assemble:

void assemble(GenericTensor& A,
              const Form& a,
              const Function& rotation_function);

Then, in the private section of the class, add the following declarations

bool walk_coefficients(std::size_t& n,
                       std::shared_ptr<const dolfin::FunctionSpace> func_space,
                       std::shared_ptr<const dolfin::FunctionSpace> rotation_space);

std::size_t value_size(std::shared_ptr<const dolfin::FunctionSpace> func_space);

struct rotation_coefficient_struct
{
  std::size_t coefficient_dimension = 0;
  std::size_t coef_idx = 0;
  std::size_t rotation_value_size = 0;
  std::size_t rotation_n_nodes = 0;
  std::size_t start_rotation_index;
};
rotation_coefficient_struct rot_coefs;
rotation_coefficient_struct rot_args[2];

void find_coef_rotation_indices(const std::vector<std::shared_ptr<const GenericFunction>>& coefficients,
                                const Function& rotation_function);
void find_rotation_indices(rotation_coefficient_struct & rot,
                           const std::size_t i,
                           const std::shared_ptr<const dolfin::FunctionSpace> func_space,
                           const Function& rotation_function);

The structure rotation_coefficient_struct is designed to sore the sought informations into rot_coefs for the coefficients of the form, rot_args[0] for the test functions and rot_args[1] for the trial functions:

• rotation_coefficient_struct::coefficient_dimension stores the space_dimension (the total number of DOfs) of the corresponding FiniteElement;

• rotation_coefficient_struct::coef_idx stores inside which sub element of the form’s function space the rotation_function is held;

• rotation_coefficient_struct::rotation_value_size is the value_size of the rotation function space (it should be equal to three)

• rotation_coefficient_struct::rotation_n_nodes is the number of nodes of the finite element; it is equal to rotation_coefficient_struct::coefficient_dimension/rotation_coefficient_struct::rotation_value_size;

• rotation_coefficient_struct::start_rotation_index is the starting index of the rotation coefficients inside the local vector.

The local vector will store, starting from index start_rotation_index, rotation_n_nodes values that corresponds to the first component of the rotation vector, then rotation_n_nodes values for the second components, and at last rotation_n_nodes values for the third component.

These helper functions below can populate variables rot_coefs and rot_args by walking the form coefficients and arguments

void RotationAssembler::find_coef_rotation_indices(
    const std::vector<std::shared_ptr<const GenericFunction>>& coefficients,
    const Function& rotation_function)
{
  for (std::size_t i = 0; i < coefficients.size(); i++)
  {
    const Function*const f = dynamic_cast<const Function*const>(coefficients[i].get());
    if (f != nullptr) {
      const std::shared_ptr<const dolfin::FunctionSpace> func_space = f->function_space();
      find_rotation_indices(rot_coefs, i, func_space, rotation_function);
    }
  }
}

void RotationAssembler::find_rotation_indices(
    rotation_coefficient_struct & rot,
    const std::size_t i,
    const std::shared_ptr<const FunctionSpace> func_space,
    const Function& rotation_function)
{
  rot.start_rotation_index = 0;
  if (walk_coefficients(rot.start_rotation_index, func_space, rotation_function.function_space()))
  {
    rot.coefficient_dimension = func_space->element()->space_dimension();
    rot.coef_idx = i;
    rot.rotation_value_size = value_size(rotation_function.function_space());
    rot.rotation_n_nodes = rotation_function.function_space()->element()->space_dimension() /
                               rot.rotation_value_size;
  }
}

std::size_t RotationAssembler::value_size(std::shared_ptr<const dolfin::FunctionSpace> func_space)
{
  std::size_t value_size = 1;
  for (std::size_t d = 0; d < func_space->element()->value_rank(); d++)
    value_size *= func_space->element()->value_dimension(d);
  return value_size;
};

bool RotationAssembler::walk_coefficients(std::size_t& n,
                       std::shared_ptr<const dolfin::FunctionSpace> func_space,
                       std::shared_ptr<const dolfin::FunctionSpace> rotation_space)
{
  bool found = false;
  for (std::size_t i = 0; i < func_space->element()->num_sub_elements(); i++) {
    if (func_space->sub(i) == rotation_space) {
      if (value_size(func_space->sub(i)) != 3)
         dolfin_error("RotationAssembler::walk_coefficients",
                               "correctly identify the rotation (sub)space",
                               "The size of rotation field is != 3: coef.sub(s_space).value_size() = %d != 3", value_size(func_space->sub(i)));
      found = true;
      break;
    }
    else if (func_space->sub(i)->element()->num_sub_elements() > 0)
    {
      found = walk_coefficients(n, func_space->sub(i), rotation_space);
    }
    else
    {
      n += func_space->sub(i)->dofmap()->max_element_dofs();
    }
  }
  return found;
};

The variables can be populated inside assemble by

// Find rotation (sub)space indices
find_coef_rotation_indices(coefficients, rotation_function);
for (std::size_t i = 0; i < ufc.form.rank(); ++i)
  find_rotation_indices(rot_args[i], i, a.function_space(i), rotation_function);

At this point assemble_cell can be modified. Let see first how can this be accomplished for forms up to rank 1.

First, immediately before tabulate_tensor save the values of the rotation vectors inside Eigen::Vector3d save_phi[rot_coefs.rotation_n_nodes] and check wheter some vector should be changed back to its complement. This is done by checking whether the scalar product between the first node rotation vector and the subsequent ones is negative or not; note that this test somewhat limits the relative rotation allowed within a single element. Whether a rotation is changed into its complement or not is stored in side the switched_node vector, while switched_cell is true only if at least one rotation was transformed

// If the cell has switched rotation dofs switch then ...
bool switched_cell = false;
double * wr = ufc.w()[rot_coefs.coef_idx] + rot_coefs.start_rotation_index;
Eigen::Vector3d save_phi[rot_coefs.rotation_n_nodes];
bool switched_node[rot_coefs.rotation_n_nodes];
switched_node[0] = false;
constexpr double half_pi = std::acos(-1.) / 2.;
constexpr double pi2 = 2. * std::acos(-1.);
if (rot_coefs.rotation_value_size)
{
  save_phi[0](0) = wr[0];
  save_phi[0](1) = wr[0 + rot_coefs.rotation_n_nodes];
  save_phi[0](2) = wr[0 + 2 * rot_coefs.rotation_n_nodes];
  if (std::sqrt(save_phi[0].dot(save_phi[0])) > half_pi)
  {
    for (std::size_t n = 1; n < rot_coefs.rotation_n_nodes; n++)
    {
      Eigen::Vector3d& phi = *(save_phi + n);
      phi(0) = wr[n];
      phi(1) = wr[n + rot_coefs.rotation_n_nodes];
      phi(2) = wr[n + 2 * rot_coefs.rotation_n_nodes];
      switched_node[n] = (save_phi[0].dot(phi) < 0);
      if (switched_node[n])
      {
        switched_cell = true;
        double normphi = std::sqrt(phi.dot(phi));
        double c = 1. - pi2 / normphi;
        wr[n] = phi(0) * c;
        wr[n + rot_coefs.rotation_n_nodes] = phi(1) * c;
        wr[n + 2 * rot_coefs.rotation_n_nodes] = phi(2) * c;
      }
    }
  }
}

It is now possible to call:

// Tabulate cell tensor
integral->tabulate_tensor(ufc.A.data(), ufc.w(),
                          coordinate_dofs.data(),
                          ufc_cell.orientation);

and transform the result

// Transform the tabulated values
if (switched_cell) {
  for (std::size_t n = 0; n < rot_coefs.rotation_n_nodes; n++)
  {
    if (switched_node[n])
    {
      Eigen::Vector3d& phi = *(save_phi + n);
      if (form_rank == 1 && rot_args[0].rotation_value_size)
      {
        const double normphi2 = phi.dot(phi);
        const double normphi = std::sqrt(normphi2);
        const double c1 = 1. - pi2 / normphi;

        double * r = ufc.A.data() + rot_args[0].start_rotation_index + n;

        const double c2 = (r[0] * phi(0) +
                    r[0 + rot_coefs.rotation_n_nodes] * phi(1) +
                    r[0 + 2 * rot_coefs.rotation_n_nodes] * phi(2)
                   ) * pi2 / normphi2 / normphi;
        r[0] = r[0] * c1 + phi(0) * c2;
        r[0 + rot_coefs.rotation_n_nodes] = r[0 + rot_coefs.rotation_n_nodes] * c1 + phi(1) * c2;
        r[0 + 2 * rot_coefs.rotation_n_nodes] = r[0 + 2 * rot_coefs.rotation_n_nodes] * c1 + phi(2) * c2;
      }
    }
  }
}

before assembling it

// Add entries to global tensor. Either store values cell-by-cell
// (currently only available for functionals)
if (is_cell_functional)
  (*values)[cell->index()] = ufc.A[0];
else
  A.add_local(ufc.A.data(), dofs);

In order to use this new assembly class from within Python copy the file assembling.py into the new file rotation_assembling.py.

In order to compile and link the new class add this code

from dolfin import compile_extension_module
from dolfin.fem.assembling import _create_dolfin_form, _create_tensor

pwd = os.path.dirname(os.path.abspath(__file__))
with open(pwd + "/RotationAssembler.cpp", "r") as f:
  rotation_assembler_code = f.read()
rotation_assembler = compile_extension_module(rotation_assembler_code)

__all__ = ["r_assemble"]

Then add the new required argument rotation_function by changing

def assemble(form,
             tensor=None,
             form_compiler_parameters=None,
             add_values=False,
             finalize_tensor=True,
             keep_diagonal=False,
             backend=None):

into

def r_assemble(form,
             rotation_function,
             tensor=None,
             form_compiler_parameters=None,
             add_values=False,
             finalize_tensor=True,
             keep_diagonal=False,
             backend=None):

Then, istead of calling the standard Assembler with

# Call C++ assemble function
assembler = cpp.Assembler()
assembler.add_values = add_values
assembler.finalize_tensor = finalize_tensor
assembler.keep_diagonal = keep_diagonal
assembler.assemble(tensor, dolfin_form)

one can call the new code with

# Call C++ assemble function
assembler = rotation_assembler.RotationAssembler()
assembler.add_values = add_values
assembler.finalize_tensor = finalize_tensor
assembler.keep_diagonal = keep_diagonal
assembler.assemble(tensor, dolfin_form, rotation_function)

# Convert to float for scalars
if dolfin_form.rank() == 0:
    tensor = tensor.get_scalar_value()

# Return value
return tensor

Note however that, in order to assemble a rank two form one should assemble also the corresponding rank one. This can be accomplished by defining a new C++ RotationAssembleSystem class that takes both forms as arguments

class RotationAssemblerSystem {
public:

void assemble(GenericTensor& A,
              const Form& a,
              GenericTensor& L,
              const Form& b,
              const Function& rotation_function);

/// Assemble tensor from given form over cells. This function is
/// provided for users who wish to build a customized assembler.
void assemble_cells(GenericTensor& A, const Form& a,
                    GenericTensor& L, const Form& b,
                    UFC& ufca,
                    UFC& ufcb,
                    std::shared_ptr<const MeshFunction<std::size_t>> domains,
                    std::vector<double>* values);

and so on.

The assemble_cell code now computes both forms

// Tabulate cell tensor
integrala->tabulate_tensor(ufca.A.data(), ufca.w(),
                          coordinate_dofs.data(),
                          ufc_cell.orientation);
integralb->tabulate_tensor(ufcb.A.data(), ufcb.w(),
                          coordinate_dofs.data(),
                          ufc_cell.orientation);

storing into ufca the local contribution from the rank two form a, into ufcb that from the rank 1 form b.

If the cell needs a transoformation

// Transform the tabulated values
if (switched_cell) {

loop over all the nodes and, if a node was switched

for (std::size_t n = 0; n < rot_coefs.rotation_n_nodes; n++)
{
  if (switched_node[n])
  {
    Eigen::Vector3d& phi = *(save_phi + n);
    const double normphi2 = phi.dot(phi);
    const double normphi = std::sqrt(normphi2);
    const double pi2_normphi3 = pi2 / normphi2 / normphi;
    const double c1 = 1. - pi2 / normphi;

If the trial function needs the transformation multiply the rank 2 tensor on the right by \(\left(\partial\boldsymbol{\varphi}_C/\partial\boldsymbol{\varphi}\right)\). Remember that the matrix is stored with a Row-Major ordering

// Rank 2 on the right (trial function)
if (rot_args[1].rotation_value_size)
{
  //Row Major
  for (std::size_t row = 0; row < dofs[0].size(); row++)
  {
    double * r = ufca.A.data() + row * dofs[1].size() + rot_args[1].start_rotation_index + n;
    const double c2 = (r[0] * phi(0) +
                 r[0 + rot_coefs.rotation_n_nodes] * phi(1) +
                 r[0 + 2 * rot_coefs.rotation_n_nodes] * phi(2)
                ) * pi2_normphi3;
    r[0] = r[0] * c1 + phi(0) * c2;
    r[0 + rot_coefs.rotation_n_nodes] = r[0 + rot_coefs.rotation_n_nodes] * c1 + phi(1) * c2;
    r[0 + 2 * rot_coefs.rotation_n_nodes] = r[0 + 2 * rot_coefs.rotation_n_nodes] * c1 + phi(2) * c2;
  }
}

If the test function needs the transformation multiply the rank 2 tensor on the left by \(\left(\partial\boldsymbol{\varphi}_C/\partial\boldsymbol{\varphi}\right)^T\)

if (rot_args[0].rotation_value_size)
{
  // Rank 2 on the left (test function)
  double * rb = ufcb.A.data() + rot_args[0].start_rotation_index + n;
  Eigen::Vector3d res(rb[0], rb[0 + rot_coefs.rotation_n_nodes], rb[0 + 2 * rot_coefs.rotation_n_nodes]);
  for (std::size_t col = 0; col < dofs[1].size(); col++)
  {
    double * r = ufca.A.data() + col + dofs[1].size() * (rot_args[0].start_rotation_index + n);

    double c2 = (r[0] * phi(0) +
                 r[0 + dofs[1].size() * rot_coefs.rotation_n_nodes] * phi(1) +
                 r[0 + 2 * dofs[1].size() * rot_coefs.rotation_n_nodes] * phi(2)
                ) * pi2_normphi3;
    r[0] = r[0] * c1 + phi(0) * c2;
    r[0 + dofs[1].size() * rot_coefs.rotation_n_nodes] =
        r[0 + dofs[1].size() * rot_coefs.rotation_n_nodes] * c1 + phi(1) * c2;
    r[0 + 2 * dofs[1].size() * rot_coefs.rotation_n_nodes] =
        r[0 + 2 * dofs[1].size() * rot_coefs.rotation_n_nodes] * c1 + phi(2) * c2;
  }

Still inside the test function block, having at hand the residual it is possible to compute the second order tensor (1) into the Eigen::Matrix3d K variable, and add it to the tensor at the right locations

Eigen::Vector3d n_phi = phi / normphi;
Eigen::Matrix3d K = (res * n_phi.transpose() + n_phi * res.transpose() +
                        n_phi.dot(res) * (Eigen::Matrix3d().setIdentity() - n_phi * n_phi.transpose() * 3.)
                    ) * (pi2 / normphi2);
double * r = ufca.A.data() + (rot_args[1].start_rotation_index + n) * dofs[1].size() + rot_args[1].start_rotation_index + n;
for (std::size_t ir = 0; ir < 3; ir++)
  for (std::size_t ic = 0; ic < 3; ic++)
    r[0 + dofs[1].size() * (ir * rot_coefs.rotation_n_nodes) + ic * rot_coefs.rotation_n_nodes] -= K(ir, ic);

Finally, tranform the rank one form

  // Rank 1 on the left (test function)
  const double c2 = res.dot(phi) * pi2 / normphi2 / normphi;
  rb[0] = res(0) * c1 + phi(0) * c2;
  rb[0 + rot_coefs.rotation_n_nodes] = res(1) * c1 + phi(1) * c2;
  rb[0 + 2 * rot_coefs.rotation_n_nodes] = res(2) * c1 + phi(2) * c2;
}

and assemble the result into the GenericTensor``s ``A and L

// Add entries to global tensor.
A.add_local(ufca.A.data(), dofs);
L.add_local(ufcb.A.data(), dofs);

The new class can be called from Python by extending again rotation_assembly.py

from dolfin import compile_extension_module
from dolfin.fem.assembling import _create_dolfin_form, _create_tensor

pwd = os.path.dirname(os.path.abspath(__file__))
with open(pwd + "/RotationAssembler.cpp", "r") as f:
  rotation_assembler_code = f.read()
rotation_assembler = compile_extension_module(rotation_assembler_code)

pwd = os.path.dirname(os.path.abspath(__file__))
with open(pwd + "/RotationAssemblerSystem.cpp", "r") as f:
    rotation_assembler_system_code = f.read()
rotation_assembler_system = compile_extension_module(rotation_assembler_system_code)

__all__ = ["r_assemble", "r_assemble_system"]

and adding a function

def r_assemble_system(forma, formb,
             rotation_function,
             tensora=None,
             tensorb=None,
             form_compiler_parameters=None,
             add_values=False,
             finalize_tensor=True,
             keep_diagonal=False,
             backend=None):

that calls rotation_assembler_system

# Create dolfin Form object referencing all data needed by assembler
dolfin_forma = _create_dolfin_form(forma, form_compiler_parameters)
dolfin_formb = _create_dolfin_form(formb, form_compiler_parameters)

# Create tensor
comma = dolfin_forma.mesh().mpi_comm()
commb = dolfin_formb.mesh().mpi_comm()
tensora = _create_tensor(comma, forma, dolfin_forma.rank(), backend, tensora)
tensorb = _create_tensor(commb, formb, dolfin_formb.rank(), backend, tensorb)

# Call C++ assemble function
assembler = rotation_assembler_system.RotationAssemblerSystem()
assembler.add_values = add_values
assembler.finalize_tensor = finalize_tensor
assembler.keep_diagonal = keep_diagonal
assembler.assemble(tensora, dolfin_forma, tensorb, dolfin_formb, rotation_function)

# Return value
return (tensora, tensorb)

Methods assemble_exterior_facets, assemble_interior_facets and assemble_vertices need to be modified too. Right now we skip them because either they are not used in this example or they do evaluate the rotation function only at a node.