DOLFIN-X
DOLFIN-X C++ interface
evaluate.h
1 // Copyright (C) 2020 Jack S. Hale
2 //
3 // This file is part of DOLFINX (https://www.fenicsproject.org)
4 //
5 // SPDX-License-Identifier: LGPL-3.0-or-later
6 //
7 
8 #pragma once
9 
10 #include <Eigen/Core>
11 #include <dolfinx/fem/utils.h>
12 #include <dolfinx/mesh/Mesh.h>
13 #include <vector>
14 
15 namespace dolfinx::fem
16 {
17 
18 template <typename T>
19 class Expression;
24 template <typename T>
25 void eval(
26  Eigen::Ref<Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
27  values,
28  const fem::Expression<T>& e, const std::vector<std::int32_t>& active_cells)
29 {
30  // Extract data from Expression
31  auto mesh = e.mesh();
32  assert(mesh);
33 
34  // Prepare coefficients
35  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
37 
38  // Prepare constants
39  const std::vector<T> constant_values = dolfinx::fem::pack_constants(e);
40 
41  const auto& fn = e.get_tabulate_expression();
42 
43  // Prepare cell geometry
45  = mesh->geometry().dofmap();
46  const fem::CoordinateElement& cmap = mesh->geometry().cmap();
47 
48  // Prepate cell permutation info
49  mesh->topology_mutable().create_entity_permutations();
50  const std::vector<std::uint32_t>& cell_info
51  = mesh->topology().get_cell_permutation_info();
52 
53  // FIXME: Add proper interface for num coordinate dofs
54  const int num_dofs_g = x_dofmap.num_links(0);
55  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
56  = mesh->geometry().x();
57 
58  // Create data structures used in evaluation
59  const int gdim = mesh->geometry().dim();
60  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
61  coordinate_dofs(num_dofs_g, gdim);
62 
63  // Iterate over cells and 'assemble' into values
64  std::vector<T> values_e(e.num_points() * e.value_size(), 0);
65  for (std::size_t i = 0; i < active_cells.size(); ++i)
66  {
67  const std::int32_t c = active_cells[i];
68  auto x_dofs = x_dofmap.links(c);
69  for (int j = 0; j < num_dofs_g; ++j)
70  {
71  const auto x_dof = x_dofs[j];
72  for (int k = 0; k < gdim; ++k)
73  coordinate_dofs(j, k) = x_g(x_dof, k);
74  }
75 
76  auto coeff_cell = coeffs.row(c);
77  std::fill(values_e.begin(), values_e.end(), 0.0);
78  fn(values_e.data(), coeff_cell.data(), constant_values.data(),
79  coordinate_dofs.data());
80 
81  for (std::size_t j = 0; j < values_e.size(); ++j)
82  values(i, j) = values_e[j];
83  }
84 }
85 
86 } // namespace dolfinx::fem
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:26
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell....
Definition: Expression.h:39
const std::size_t value_size() const
Get value size.
Definition: Expression.h:136
std::shared_ptr< const mesh::Mesh > mesh() const
Get mesh.
Definition: Expression.h:124
const std::function< void(T *, const T *, const T *, const double *)> & get_tabulate_expression() const
Get function for tabulate_expression.
Definition: Expression.h:108
const Eigen::Index num_points() const
Get number of points.
Definition: Expression.h:140
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:46
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:128
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:118
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:455
Eigen::Array< typename U::scalar_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:400
void eval(Eigen::Ref< Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> values, const fem::Expression< T > &e, const std::vector< std::int32_t > &active_cells)
Evaluate a UFC expression.
Definition: evaluate.h:25