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 <vector>
11 
12 #include <Eigen/Dense>
13 #include <dolfinx/fem/utils.h>
14 #include <dolfinx/mesh/Mesh.h>
15 
16 namespace dolfinx::function
17 {
18 
19 template <typename T>
20 class Expression;
25 template <typename T>
26 void eval(
27  Eigen::Ref<Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
28  values,
29  const function::Expression<T>& e,
30  const std::vector<std::int32_t>& active_cells)
31 {
32  // Extract data from Expression
33  auto mesh = e.mesh();
34  assert(mesh);
35 
36  // Prepare coefficients
37  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
39 
40  // Prepare constants
41  const Eigen::Array<T, Eigen::Dynamic, 1> constant_values
43 
44  const auto& fn = e.get_tabulate_expression();
45 
46  // Prepare cell geometry
48  = mesh->geometry().dofmap();
49 
50  // FIXME: Add proper interface for num coordinate dofs
51  const int num_dofs_g = x_dofmap.num_links(0);
52  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
53  = mesh->geometry().x();
54 
55  // Create data structures used in evaluation
56  const int gdim = mesh->geometry().dim();
57  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
58  coordinate_dofs(num_dofs_g, gdim);
59 
60  Eigen::Array<T, Eigen::Dynamic, 1> values_e;
61  const Eigen::Index num_points = e.num_points();
62  const Eigen::Index value_size = e.value_size();
63  const Eigen::Index size = num_points * value_size;
64  values_e.setZero(size);
65 
66  // Iterate over cells and 'assemble' into values
67  Eigen::Index i = 0;
68  for (std::int32_t c : active_cells)
69  {
70  auto x_dofs = x_dofmap.links(c);
71  for (Eigen::Index j = 0; j < num_dofs_g; ++j)
72  {
73  const auto x_dof = x_dofs[j];
74  for (Eigen::Index k = 0; k < gdim; ++k)
75  coordinate_dofs(j, k) = x_g(x_dof, k);
76  }
77 
78  auto coeff_cell = coeffs.row(c);
79 
80  // Experimentally faster than .setZero().
81  for (Eigen::Index j = 0; j < size; j++)
82  values_e(j) = 0.0;
83 
84  fn(values_e.data(), coeff_cell.data(), constant_values.data(),
85  coordinate_dofs.data());
86 
87  values.row(i) = values_e;
88  ++i;
89  }
90 }
91 
92 } // namespace dolfinx::function
dolfinx::graph::AdjacencyList::num_links
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:143
dolfinx::function::Expression::value_size
const std::size_t value_size() const
Get value size.
Definition: Expression.h:139
dolfinx::fem::pack_coefficients
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:407
dolfinx::function::eval
void eval(Eigen::Ref< Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> values, const function::Expression< T > &e, const std::vector< std::int32_t > &active_cells)
Evaluate a UFC expression.
Definition: evaluate.h:26
dolfinx::function::Expression::get_tabulate_expression
const std::function< void(T *, const T *, const T *, const double *)> & get_tabulate_expression() const
Get function for tabulate_expression.
Definition: Expression.h:110
dolfinx::graph::AdjacencyList::links
Eigen::Array< T, Eigen::Dynamic, 1 >::SegmentReturnType links(int node)
Links (edges) for given node.
Definition: AdjacencyList.h:153
dolfinx::graph::AdjacencyList
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:28
dolfinx::function::Expression::mesh
std::shared_ptr< const mesh::Mesh > mesh() const
Get mesh.
Definition: Expression.h:127
dolfinx::function::Expression
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell....
Definition: Expression.h:39
dolfinx::function
Functions tools, including FEM functions and pointwise defined functions.
Definition: assembler.h:19
dolfinx::fem::pack_constants
Eigen::Array< typename U::scalar_type, Eigen::Dynamic, 1 > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:456
dolfinx::function::Expression::num_points
const Eigen::Index num_points() const
Get number of points.
Definition: Expression.h:143